SUNDIALS的C++使用例子

SUNDIALS的C++使用例子

本例子基于cvAdvDiff_bnd.c改造。可以作为sundials使用的模板

本例子采用CMake进行编译，在ubuntun18下测试通过

特点

1. 回调的函数是类的静态成员函数。不能够是成员函数，因为接口需要的是函数指针，而成员函数需要对象.
2. 可以支持不需要jac
3. 起始时间可以是任意时间,更改T0的值，以及main.cpp的迭代起始值即可
4. 测试的结果和原始代码一致

---

##头文件

#ifndef CV_ADV_DIFF_BND_H_
#define CV_ADV_DIFF_BND_H_

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <cvode/cvode.h>              /* prototypes for CVODE fcts., consts.  */
#include <nvector/nvector_serial.h>   /* access to serial N_Vector            */
#include <sunmatrix/sunmatrix_band.h> /* access to band SUNMatrix             */
#include <sunlinsol/sunlinsol_band.h> /* access to band SUNLinearSolver       */
#include <sundials/sundials_types.h>  /* definition of type realtype          */
#include <sundials/sundials_math.h>   /* definition of ABS and EXP            */

class CvAdvDiffBnd
{

public:
    CvAdvDiffBnd();
    virtual ~CvAdvDiffBnd();

 

public:
    int Init();
    int Run(double tout);
    void Finish();
    void Release();

public:
    /**
     * 1. Initialize parallel or multi-threaded environment, if appropriate.
     */
    int init_environment_1();

    /**
     *2. Defining the length of the problem.
     */
    int define_problem_length_2();

    /**
     *  3. Set vector of initial values.
     */
    int set_vector_initial_val_3();

    /**
     * 4. Create CVODE Object.
     */
    int create_vode_obj_4();

    /**
     * 5. Initialize CVODE solver.
     */
    int init_cvode_solver_5();

    /**
     *  6. Specify integration tolerances.
     */
    int specify_integration_tolerances_6();

    /**
     * 7. Set Optional inputs.
     */
    int set_optional_inputs_7();

    /**
     * 8. Create Matrix Object.
     */
    int create_matrix_ob_8();

    /**
     * 9. Create Linear Solver Object.
     */
    int create_linear_solver_obj_9();

    /**
     * 10. Set linear solver optional inputs.
     */
    int set_linear_solver_optional_inputs_10();

    /**
     *11. Attach linear solver module. 
     */
    int attach_linear_solver_module_11();

    /**
     *  12. Set linear solver interface optional inputs.
     */
    int set_linear_solver_interface_optional_inputs_12();

    /**
     * 13. Specify rootfinding problem.
     */
    int specify_rootfinding_problem_13();

    /**
     * 14. Advance solution in time.
     */
    int advance_solution_in_time_14();

    //----------------------------------------------------------------
    /**
     * 15. Get optional outputs.
     */
    int get_optional_output_15();

    /**
     * 16. Deallocate memory for solution vector.
     */
    int deallocate_mem_16();

    /**
     * 17. Free solver memory.
     */
    int free_solver_mem_17();

    /**
     * 18. Free linear solver and matrix memory.
     */
    int free_linear_solver_and_matrix_mem_18();


public:
    static int f_sub(realtype t, N_Vector u, N_Vector udot, void *user_data);

    static int Jac_sub(realtype t, N_Vector u, N_Vector fu,
        SUNMatrix J, void *user_data,
        N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);


public:
    /* Functions Called by the Solver */

    int f_func(realtype t, N_Vector u, N_Vector udot);
    int Jac_func(realtype t, N_Vector u, N_Vector fu,
            SUNMatrix J,
            N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);

public:
    /* Private Helper Functions */
    void SetIC_func(N_Vector u);
    void PrintHeader_func(realtype reltol, realtype abstol, realtype umax);
    void PrintOutput_func(realtype t, realtype umax, long int nst);
    void PrintFinalStats_func(void *cvode_mem);

    /* Private function to check function return values */
    int check_retval(void *returnvalue, const char *funcname, int opt);

 
private:
    double XMAX;  // RCONST(2.0)    /* domain boundaries         */
    double YMAX;  //RCONST(1.0)
    int MX;       //    10             /* mesh dimensions           */
    int MY;       //   5
    int NEQ;      //  MX*MY          /* number of equations       */
    double ATOL;  //  RCONST(1.0e-5) /* scalar absolute tolerance */
    double T0;    //  RCONST(0.0)    /* initial time              */
    double T1;    //   RCONST(0.1)    /* first output time         */
    double DTOUT; // RCONST(0.1)    /* output time increment     */
    int NOUT;     //  10             /* number of output times    */

    double ZERO; // RCONST(0.0)
    double HALF; // RCONST(0.5)
    double ONE;  //  RCONST(1.0)
    double TWO;  // RCONST(2.0)
    double FIVE; // RCONST(5.0)

private:
    realtype dx, dy, reltol, abstol, t, tout, umax;
    N_Vector u;
    SUNMatrix A;
    SUNLinearSolver LS;
    void *cvode_mem;
    int iout, retval;
    long int nst;

private:
    //user data
    realtype m_dx, m_dy, m_hdcoef, m_hacoef, m_vdcoef;
};

#endif

---

cpp文件

#include "CvAdvDiffBnd.h"

#define IJth(vdata, i, j) (vdata[(j - 1) + (i - 1) * MY])


int CvAdvDiffBnd::f_sub(realtype t, N_Vector u, N_Vector udot, void *user_data)
{
    CvAdvDiffBnd *pt = (CvAdvDiffBnd *)user_data;
    return pt->f_func(t, u, udot);
}

int CvAdvDiffBnd::Jac_sub(realtype t, N_Vector u, N_Vector fu,
        SUNMatrix J, void *user_data,
        N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
{
    printf("i am jac func\n");
    CvAdvDiffBnd *pt = (CvAdvDiffBnd *)user_data;
    return pt->Jac_func(t, u, fu, J, tmp1, tmp2, tmp3);
}




CvAdvDiffBnd::CvAdvDiffBnd()
{
}

CvAdvDiffBnd::~CvAdvDiffBnd()
{
    Release();
}

int CvAdvDiffBnd::Init()
{

    init_environment_1();
    /**
     * 1. Initialize parallel or multi-threaded environment, if appropriate.
     */
      init_environment_1();

    /**
     *2. Defining the length of the problem.
     */
      define_problem_length_2();

    /**
     *  3. Set vector of initial values.
     */
      set_vector_initial_val_3();

    /**
     * 4. Create CVODE Object.
     */
      create_vode_obj_4();

    /**
     * 5. Initialize CVODE solver.
     */
      init_cvode_solver_5();

    /**
     *  6. Specify integration tolerances.
     */
      specify_integration_tolerances_6();

    /**
     * 7. Set Optional inputs.
     */
      set_optional_inputs_7();

    /**
     * 8. Create Matrix Object.
     */
      create_matrix_ob_8();

    /**
     * 9. Create Linear Solver Object.
     */
      create_linear_solver_obj_9();

    /**
     * 10. Set linear solver optional inputs.
     */
      set_linear_solver_optional_inputs_10();

    /**
     *11. Attach linear solver module. 
     */
      attach_linear_solver_module_11();

    /**
     *  12. Set linear solver interface optional inputs.
     */
      set_linear_solver_interface_optional_inputs_12();

    /**
     * 13. Specify rootfinding problem.
     */
      specify_rootfinding_problem_13();

    /**
     * 14. Advance solution in time.
     */
      advance_solution_in_time_14();
}

int CvAdvDiffBnd::Run(double tout)
{
    retval = CVode(cvode_mem, tout, u, &t, CV_NORMAL);
    if (check_retval(&retval, "CVode", 1))
    {
        return -1;
    }

    umax = N_VMaxNorm(u);
    retval = CVodeGetNumSteps(cvode_mem, &nst);
    check_retval(&retval, "CVodeGetNumSteps", 1);
    PrintOutput_func(t, umax, nst);

    return 0;
}

void CvAdvDiffBnd::Finish()
{
    //----------------------------------------------------------------
    /**
     * 15. Get optional outputs.
     */
    get_optional_output_15();
    PrintFinalStats_func(cvode_mem); /* Print some final statistics   */
}

void CvAdvDiffBnd::Release()
{

    /**
     * 16. Deallocate memory for solution vector.
     */
    deallocate_mem_16();

    /**
     * 17. Free solver memory.
     */
    free_solver_mem_17();

    /**
     * 18. Free linear solver and matrix memory.
     */
    free_linear_solver_and_matrix_mem_18();
}

int CvAdvDiffBnd::init_environment_1()
{
    XMAX = RCONST(2.0); /* domain boundaries         */
    YMAX = RCONST(1.0);
    MX = 10; /* mesh dimensions           */
    MY = 5;

    ATOL = RCONST(1.0e-5); /* scalar absolute tolerance */
    T0 = RCONST(0.5);      /* initial time              */
    T1 = RCONST(0.1);      /* first output time         */
    DTOUT = RCONST(0.1);   /* output time increment     */
    NOUT = 10;             /* number of output times    */

    ZERO = RCONST(0.0);
    HALF = RCONST(0.5);
    ONE = RCONST(1.0);
    TWO = RCONST(2.0);
    FIVE = RCONST(5.0);

    u = NULL;
    A = NULL;
    LS = NULL;
    cvode_mem = NULL;
    return 0;
}

/**
     *2. Defining the length of the problem.
     */
int CvAdvDiffBnd::define_problem_length_2()
{
    NEQ = MX * MY; /* number of equations       */
    return 0;
}

/**
     *  3. Set vector of initial values.
     */
int CvAdvDiffBnd::set_vector_initial_val_3()
{
    u = N_VNew_Serial(NEQ); /* Allocate u vector */
    if (check_retval((void *)u, "N_VNew_Serial", 0))
        return (1);

    dx = m_dx = XMAX / (MX + 1); /* Set grid coefficients in data */
    dy = m_dy = YMAX / (MY + 1);
    m_hdcoef = ONE / (dx * dx);
    m_hacoef = HALF / (TWO * dx);
    m_vdcoef = ONE / (dy * dy);

    SetIC_func(u);
    return 0;
}

/**
     * 4. Create CVODE Object.
     */
int CvAdvDiffBnd::create_vode_obj_4()
{
    /* Call CVodeCreate to create the solver memory and specify the 
   * Backward Differentiation Formula */
    cvode_mem = CVodeCreate(CV_BDF);
    if (check_retval((void *)cvode_mem, "CVodeCreate", 0))
        return (1);

    return 0;
}

/**
     * 5. Initialize CVODE solver.
     */
int CvAdvDiffBnd::init_cvode_solver_5()
{
    /* Call CVodeInit to initialize the integrator memory and specify the
   * user's right hand side function in u'=f(t,u), the inital time T0, and
   * the initial dependent variable vector u. */
    //retval = CVodeInit(cvode_mem, f, T0, u);
    retval = CVodeInit(cvode_mem, CvAdvDiffBnd::f_sub, T0, u);
    if (check_retval(&retval, "CVodeInit", 1))
        return (1);

    return 0;
}

/**
     *  6. Specify integration tolerances.
     */
int CvAdvDiffBnd::specify_integration_tolerances_6()
{
    reltol = ZERO; /* Set the tolerances */
    abstol = ATOL;
    /* Call CVodeSStolerances to specify the scalar relative tolerance
   * and scalar absolute tolerance */
    retval = CVodeSStolerances(cvode_mem, reltol, abstol);
    if (check_retval(&retval, "CVodeSStolerances", 1))
        return (1);

    return 0;
}

/**
     * 7. Set Optional inputs.
     */
int CvAdvDiffBnd::set_optional_inputs_7()
{

    return 0;
}

/**
     * 8. Create Matrix Object.
     */
int CvAdvDiffBnd::create_matrix_ob_8()
{
    /* Set the pointer to user-defined data */
    retval = CVodeSetUserData(cvode_mem, this);
    if (check_retval(&retval, "CVodeSetUserData", 1))
        return (1);

    /* Create banded SUNMatrix for use in linear solves -- since this will be factored, 
     set the storage bandwidth to be the sum of upper and lower bandwidths */
    A = SUNBandMatrix(NEQ, MY, MY);
    if (check_retval((void *)A, "SUNBandMatrix", 0))
        return (1);

    return 0;
}

/**
     * 9. Create Linear Solver Object.
     */
int CvAdvDiffBnd::create_linear_solver_obj_9()
{
    /* Create banded SUNLinearSolver object for use by CVode */
    LS = SUNLinSol_Band(u, A);
    if (check_retval((void *)LS, "SUNLinSol_Band", 0))
        return (1);

    return 0;
}

/**
     * 10. Set linear solver optional inputs.
     */
int CvAdvDiffBnd::set_linear_solver_optional_inputs_10()
{
    return 0;
}

/**
     *11. Attach linear solver module. 
     */
int CvAdvDiffBnd::attach_linear_solver_module_11()
{
    /* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
    retval = CVodeSetLinearSolver(cvode_mem, LS, A);
    if (check_retval(&retval, "CVodeSetLinearSolver", 1))
        return (1);

    umax = N_VMaxNorm(u);
    PrintHeader_func(reltol, abstol, umax);

    return 0;
}

/**
     *  12. Set linear solver interface optional inputs.
     */
int CvAdvDiffBnd::set_linear_solver_interface_optional_inputs_12()
{
    /* Set the user-supplied Jacobian routine Jac */
    //retval = CVodeSetJacFn(cvode_mem, Jac);
    
    retval = CVodeSetJacFn(cvode_mem, Jac_sub);
    if (check_retval(&retval, "CVodeSetJacFn", 1))
        return (1);

    return 0;
}

/**
     * 13. Specify rootfinding problem.
     */
int CvAdvDiffBnd::specify_rootfinding_problem_13()
{
    return 0;
}

/**
     * 14. Advance solution in time.
     */
int CvAdvDiffBnd::advance_solution_in_time_14()
{
    return 0;
}

//----------------------------------------------------------------
/**
     * 15. Get optional outputs.
     */
int CvAdvDiffBnd::get_optional_output_15()
{
    return 0;
}

/**
     * 16. Deallocate memory for solution vector.
     */
int CvAdvDiffBnd::deallocate_mem_16()
{
    N_VDestroy(u); /* Free the u vector    */
    return 0;
}

/**
     * 17. Free solver memory.
     */
int CvAdvDiffBnd::free_solver_mem_17()
{
    CVodeFree(&cvode_mem); /* Free the integrator memory */
    SUNLinSolFree(LS);     /* Free linear solver memory  */
    SUNMatDestroy(A);      /* Free the matrix memory     */
    return 0;
}

/**
     * 18. Free linear solver and matrix memory.
     */
int CvAdvDiffBnd::free_linear_solver_and_matrix_mem_18()
{
    return 0;
}

/*
 *-------------------------------
 * Functions called by the solver
 *-------------------------------
 */

/* f routine. Compute f(t,u). */

int CvAdvDiffBnd::f_func(realtype t, N_Vector u, N_Vector udot)
{
    realtype uij, udn, uup, ult, urt, hordc, horac, verdc, hdiff, hadv, vdiff;
    realtype *udata, *dudata;
    int i, j;
    //UserData data;

    udata = N_VGetArrayPointer(u);
    dudata = N_VGetArrayPointer(udot);

    /* Extract needed constants from data */
    //data = (UserData)user_data;
    hordc = m_hdcoef;
    horac = m_hacoef;
    verdc = m_vdcoef;

    /* Loop over all grid points. */

    for (j = 1; j <= MY; j++)
    {

        for (i = 1; i <= MX; i++)
        {

            /* Extract u at x_i, y_j and four neighboring points */

            uij = IJth(udata, i, j);
            udn = (j == 1) ? ZERO : IJth(udata, i, j - 1);
            uup = (j == MY) ? ZERO : IJth(udata, i, j + 1);
            ult = (i == 1) ? ZERO : IJth(udata, i - 1, j);
            urt = (i == MX) ? ZERO : IJth(udata, i + 1, j);

            /* Set diffusion and advection terms and load into udot */

            hdiff = hordc * (ult - TWO * uij + urt);
            hadv = horac * (urt - ult);
            vdiff = verdc * (uup - TWO * uij + udn);
            IJth(dudata, i, j) = hdiff + hadv + vdiff;
        }
    }

    return (0);
}

/* Jacobian routine. Compute J(t,u). */

int CvAdvDiffBnd::Jac_func(realtype t, N_Vector u, N_Vector fu,
                           SUNMatrix J,
                           N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
{
    sunindextype i, j, k;
    realtype *kthCol, hordc, horac, verdc;
    //UserData data;

    /*
   * The components of f = udot that depend on u(i,j) are
   * f(i,j), f(i-1,j), f(i+1,j), f(i,j-1), f(i,j+1), with
   *   df(i,j)/du(i,j) = -2 (1/dx^2 + 1/dy^2)
   *   df(i-1,j)/du(i,j) = 1/dx^2 + .25/dx  (if i > 1)
   *   df(i+1,j)/du(i,j) = 1/dx^2 - .25/dx  (if i < MX)
   *   df(i,j-1)/du(i,j) = 1/dy^2           (if j > 1)
   *   df(i,j+1)/du(i,j) = 1/dy^2           (if j < MY)
   */

    //data = (UserData)user_data;
    hordc = m_hdcoef;
    horac = m_hacoef;
    verdc = m_vdcoef;

    /* set non-zero Jacobian entries */
    for (j = 1; j <= MY; j++)
    {
        for (i = 1; i <= MX; i++)
        {
            k = j - 1 + (i - 1) * MY;
            kthCol = SUNBandMatrix_Column(J, k);

            /* set the kth column of J */

            SM_COLUMN_ELEMENT_B(kthCol, k, k) = -TWO * (verdc + hordc);
            if (i != 1)
                SM_COLUMN_ELEMENT_B(kthCol, k - MY, k) = hordc + horac;
            if (i != MX)
                SM_COLUMN_ELEMENT_B(kthCol, k + MY, k) = hordc - horac;
            if (j != 1)
                SM_COLUMN_ELEMENT_B(kthCol, k - 1, k) = verdc;
            if (j != MY)
                SM_COLUMN_ELEMENT_B(kthCol, k + 1, k) = verdc;
        }
    }

    return (0);
}

/*
 *-------------------------------
 * Private helper functions
 *-------------------------------
 */

/* Set initial conditions in u vector */

void CvAdvDiffBnd::SetIC_func(N_Vector u)
{
    int i, j;
    realtype x, y, dx, dy;
    realtype *udata;

    /* Extract needed constants from data */

    dx = m_dx;
    dy = m_dy;

    /* Set pointer to data array in vector u. */

    udata = N_VGetArrayPointer(u);

    /* Load initial profile into u vector */

    for (j = 1; j <= MY; j++)
    {
        y = j * dy;
        for (i = 1; i <= MX; i++)
        {
            x = i * dx;
            IJth(udata, i, j) = x * (XMAX - x) * y * (YMAX - y) * SUNRexp(FIVE * x * y);
        }
    }
}

/* Print first lines of output (problem description) */

void CvAdvDiffBnd::PrintHeader_func(realtype reltol, realtype abstol, realtype umax)
{
    printf("\n2-D Advection-Diffusion Equation\n");
    printf("Mesh dimensions = %d X %d\n", MX, MY);
    printf("Total system size = %d\n", NEQ);
#if defined(SUNDIALS_EXTENDED_PRECISION)
    printf("Tolerance parameters: reltol = %Lg   abstol = %Lg\n\n", reltol, abstol);
    printf("At t = %Lg      max.norm(u) =%14.6Le \n", T0, umax);
#elif defined(SUNDIALS_DOUBLE_PRECISION)
    printf("Tolerance parameters: reltol = %g   abstol = %g\n\n", reltol, abstol);
    printf("At t = %g      max.norm(u) =%14.6e \n", T0, umax);
#else
    printf("Tolerance parameters: reltol = %g   abstol = %g\n\n", reltol, abstol);
    printf("At t = %g      max.norm(u) =%14.6e \n", T0, umax);
#endif

    return;
}

/* Print current value */

void CvAdvDiffBnd::PrintOutput_func(realtype t, realtype umax, long int nst)
{
#if defined(SUNDIALS_EXTENDED_PRECISION)
    printf("At t = %4.2Lf   max.norm(u) =%14.6Le   nst = %4ld\n", t, umax, nst);
#elif defined(SUNDIALS_DOUBLE_PRECISION)
    printf("At t = %4.2f   max.norm(u) =%14.6e   nst = %4ld\n", t, umax, nst);
#else
    printf("At t = %4.2f   max.norm(u) =%14.6e   nst = %4ld\n", t, umax, nst);
#endif

    return;
}

/* Get and print some final statistics */

void CvAdvDiffBnd::PrintFinalStats_func(void *cvode_mem)
{
    int retval;
    long int nst, nfe, nsetups, netf, nni, ncfn, nje, nfeLS;

    retval = CVodeGetNumSteps(cvode_mem, &nst);
    check_retval(&retval, "CVodeGetNumSteps", 1);
    retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);
    check_retval(&retval, "CVodeGetNumRhsEvals", 1);
    retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
    check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);
    retval = CVodeGetNumErrTestFails(cvode_mem, &netf);
    check_retval(&retval, "CVodeGetNumErrTestFails", 1);
    retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
    check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);
    retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
    check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);

    retval = CVodeGetNumJacEvals(cvode_mem, &nje);
    check_retval(&retval, "CVodeGetNumJacEvals", 1);
    retval = CVodeGetNumLinRhsEvals(cvode_mem, &nfeLS);
    check_retval(&retval, "CVodeGetNumLinRhsEvals", 1);

    printf("\nFinal Statistics:\n");
    printf("nst = %-6ld nfe  = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld\n",
           nst, nfe, nsetups, nfeLS, nje);
    printf("nni = %-6ld ncfn = %-6ld netf = %ld\n",
           nni, ncfn, netf);

    return;
}

/* Check function return value...
     opt == 0 means SUNDIALS function allocates memory so check if
              returned NULL pointer
     opt == 1 means SUNDIALS function returns an integer value so check if
              retval < 0
     opt == 2 means function allocates memory so check if returned
              NULL pointer */

int CvAdvDiffBnd::check_retval(void *returnvalue, const char *funcname, int opt)
{
    int *retval;

    /* Check if SUNDIALS function returned NULL pointer - no memory allocated */

    if (opt == 0 && returnvalue == NULL)
    {
        fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
                funcname);
        return (1);
    }

    /* Check if retval < 0 */

    else if (opt == 1)
    {
        retval = (int *)returnvalue;
        if (*retval < 0)
        {
            fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
                    funcname, *retval);
            return (1);
        }
    }

    /* Check if function returned NULL pointer - no memory allocated */

    else if (opt == 2 && returnvalue == NULL)
    {
        fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
                funcname);
        return (1);
    }

    return (0);
}

---

main函数

#include "CvAdvDiffBnd.h"



int main()
{
	CvAdvDiffBnd test;
    test.Init();
 
 int NOUT = 10;
 double T1 = 0.6;
 double DTOUT = 0.1;
 int iout = 0;
 double tout = 0.0;
 for(iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {
     int ret = test.Run(tout);
     if(0 != ret)
     {
         break;
     }
 }

 test.Finish();

  return 0;
}

---

##测试的结果

跟原始一致