SUNDIALS的C++使用例子
彭星津
SUNDIALS的C++使用例子
本例子基于cvAdvDiff_bnd.c改造。可以作为sundials使用的模板
本例子采用CMake进行编译，在ubuntun18下测试通过
特点
- 回调的函数是类的静态成员函数。不能够是成员函数，因为接口需要的是函数指针，而成员函数需要对象.
- 可以支持不需要jac
- 起始时间可以是任意时间,更改T0的值，以及main.cpp的迭代起始值即可
- 测试的结果和原始代码一致
##头文件
#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;
}
##测试的结果
跟原始一致
类似资料：