source: branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/TestCvodes.cs @ 16227

using System;
using HeuristicLab.Problems.DynamicalSystemsModelling;
using Microsoft.VisualStudio.TestTools.UnitTesting;

namespace AutoDiffForDynamicalModelsTest {
  [TestClass]
  public class TestCvodes {
    [TestMethod]
    public unsafe void TestCvodesMethod() {
      // test vectors
      var arr = new double[] { 3.14, 2.71 };
      var vec = Problem.N_VMake_Serial(2, arr);
      Problem.N_VDestroy_Serial(vec);

      vec = Problem.N_VNew_Serial(10);

      unsafe {
        int* content = (int*)*(int*)vec.ToPointer();
        long length = *content;
        Console.WriteLine(*content);
        int own_data = *(content + 2);
        Console.WriteLine(own_data);
        double* data = (double*) *(content + 3);
        double data0 = *data;
      }
      Problem.N_VConst_Serial(2.0, vec);
      Problem.N_VPrint_Serial(vec);
      Assert.AreEqual(20, Problem.N_VL1Norm_Serial(vec));
      Problem.N_VDestroy_Serial(vec);


      // linear oscillator
      int numberOfEquations = 2;
      var y = Problem.N_VNew_Serial(numberOfEquations);
      // y must be initialized before calling CVodeInit
      // Problem.N_VConst_Serial(100.0, y);
      Problem.NV_Set_Ith_S(y, 0, 0.5); // x
      Problem.NV_Set_Ith_S(y, 1, 1);  // v

      var cvode_mem = Problem.CVodeCreate(Problem.MultistepMethod.CV_ADAMS, Problem.NonlinearSolverIteration.CV_FUNCTIONAL);

     var flag = Problem.CVodeInit(cvode_mem, F, 0.0, y);
      Assert.AreEqual(Problem.CV_SUCCESS, flag);

      flag = Problem.CVodeSStolerances(cvode_mem, 1E-4, 1.0);
      Assert.AreEqual(Problem.CV_SUCCESS, flag);

      var A = Problem.SUNDenseMatrix(numberOfEquations, numberOfEquations);
      Assert.AreNotSame(A, IntPtr.Zero);


      var linearSolver = Problem.SUNDenseLinearSolver(y, A); 
      Assert.AreNotSame(linearSolver, IntPtr.Zero);

      flag = Problem.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
      Assert.AreEqual(Problem.CV_SUCCESS, flag);

      flag = Problem.CVDlsSetJacFn(cvode_mem, JacF);
      Assert.AreEqual(Problem.CV_SUCCESS, flag);


      double t = 0.0;
      double tout = 0.1; // first output time
      int nout = 100; // number of output times
      for (int iout = 0; iout < nout; iout++) {
        flag = Problem.CVode(cvode_mem, tout, y, ref t, Problem.CV_NORMAL);
        Assert.AreEqual(Problem.CV_SUCCESS, flag);
        Console.WriteLine("{0} {1} {2}", t, Problem.NV_Get_Ith_S(y, 0), Problem.NV_Get_Ith_S(y, 1));
        // Problem.N_VPrint_Serial(y);
        tout += 0.1; 
      }

      Problem.N_VDestroy_Serial(y);
      Problem.CVodeFree(cvode_mem);
      Problem.SUNLinSolFree(linearSolver);
      Problem.SUNMatDestroy(A);
    }

    private int JacF(double t, IntPtr y, IntPtr fy, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
      // TODO
      return 0;
    }

    public static int F(
      double t, // realtype
      IntPtr y, // N_Vector
      IntPtr ydot, // N_Vector
      IntPtr user_data) {

      Problem.NV_Set_Ith_S(ydot, 0, Problem.NV_Get_Ith_S(y, 1));
      Problem.NV_Set_Ith_S(ydot, 1, -0.3 * Problem.NV_Get_Ith_S(y, 0));
      return 0;
      ;
    }
  }
}