Sundial的c++例子
宋华灿
sundials自带的例子都是c接口，不太适合我的项目中使用，所以
尝试一下怎么用c++的形式来包装一下，方便使用
ODE原始来源
由OdeInt的例子修改而来：https://www.boost.org/doc/libs/1_67_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/tutorial/stiff_systems.html
后续目标
1 测试一下没有jac的效果
2 和OdeInt做比较
3 性能的评估，之前有做过简单的评估，SUNDIALS大概是OdeInt的3-4倍。
4 SUNDIALS的封装
头文件
#ifndef SimpleCvodeEample_h
#define SimpleCvodeEample_h
#include <iostream>
#include <cvode/cvode.h> // prototypes for CVODE fcts., consts.
#include <nvector/nvector_serial.h> // access to serial N_Vector
#include <sunlinsol/sunlinsol_spgmr.h> //access to SPGMR SUNLinearSolver
#include <cvode/cvode_spils.h> // access to CVSpils interface
#include <sundials/sundials_types.h> // defs. of realtype, sunindextype
#include <sundials/sundials_math.h> // contains the macros ABS, SUNSQR, EXP
class SimpleCvodeEample
{
public:
SimpleCvodeEample();
virtual ~SimpleCvodeEample();
public:
static SimpleCvodeEample *GetInstance();
public:
int Init();
int Run(double tout);
void Release();
public:
int f_func(realtype t, N_Vector u, N_Vector u_dot, void *user_data);
int jtv_func(N_Vector v, N_Vector Jv, realtype t, N_Vector u, N_Vector fu,
void *user_data, N_Vector tmp);
int check_flag_func(void *flagvalue, const char *funcname, int opt);
public:
int flag; // For checking if functions have run properly
realtype abstol; // real tolerance of system
realtype reltol; // absolute tolerance of system
sunindextype N;
N_Vector y; // Problem vector.
void *cvode_mem;
realtype t0;
SUNLinearSolver LS;
int print_steps;
realtype tout;
realtype end_time;
realtype step_length;
realtype t;
};
#endif
cpp文件
#include "SimpleCvodeEample.h"
#define NV_Ith_S(v,i) ( NV_DATA_S(v)[i] )
// Simple function that calculates the differential equation.
static int f(realtype t, N_Vector u, N_Vector u_dot, void *user_data)
{
return SimpleCvodeEample::GetInstance()->f_func(t,u,u_dot,user_data);
}
// Jacobian function vector routine.
static int jtv(N_Vector v, N_Vector Jv, realtype t, N_Vector u, N_Vector fu,
void *user_data, N_Vector tmp)
{
return SimpleCvodeEample::GetInstance()->jtv_func(v,Jv,t,u,fu,user_data,tmp);
}
static int check_flag(void *flagvalue, const char *funcname, int opt)
{
return SimpleCvodeEample::GetInstance()->check_flag_func(flagvalue,funcname,opt);
}
SimpleCvodeEample::SimpleCvodeEample()
{
}
SimpleCvodeEample::~SimpleCvodeEample()
{
Release();
}
SimpleCvodeEample *SimpleCvodeEample::GetInstance()
{
static SimpleCvodeEample *ins = new SimpleCvodeEample();
return ins;
}
int SimpleCvodeEample::Init()
{
abstol = 1e-5; // real tolerance of system
reltol = 1e-5; // absolute tolerance of system
// 1. Initialize parallel or multi-threaded environment, if appropriate.
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// 2. Defining the length of the problem.
N = 2;
// 3. Set vector of initial values.
// realtype y_0[N] = {2.0, 1.0};
// y = N_VMake_Serial(N, y_0);
y = N_VNew_Serial(N);
NV_Ith_S(y, 0) = 2.0;
NV_Ith_S(y, 1) = 1.0;
// 4. Create CVODE Object.
cvode_mem = CVodeCreate(CV_BDF);
// 5. Initialize CVODE solver.
t0=0;
flag = CVodeInit(cvode_mem, f, t0, y);
if(check_flag(&flag, "CVodeSetUserData", 1)) return(1);
// 6. Specify integration tolerances.
flag = CVodeSStolerances(cvode_mem, reltol, abstol);
if (check_flag(&flag, "CVodeSStolerances", 1)) return(1);
// 7. Set Optional inputs.
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// 8. Create Matrix Object.
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// 9. Create Linear Solver Object.
LS = SUNSPGMR(y, 0, 0);
if(check_flag((void *)LS, "SUNSPGMR", 0)) return(1);
// 10. Set linear solver optional inputs.
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// 11. Attach linear solver module.
flag = CVSpilsSetLinearSolver(cvode_mem, LS);
if (check_flag(&flag, "CVSpilsSetLinearSolver", 1)) return 1;
// 12. Set linear solver interface optional inputs.
flag = CVSpilsSetJacTimes(cvode_mem, NULL, jtv);
if(check_flag(&flag, "CVSpilsSetJacTimes", 1)) return(1);
// 13. Specify rootfinding problem.
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// 14. Advance solution in time.
// ---------------------------------------------------------------------------
// Have the solution advance over time, but stop to log 100 of the steps.
print_steps = 100;
end_time = 50;
step_length = 0.5;
t = 0;
}
int SimpleCvodeEample::Run(double tout)
{
realtype t = 0;
flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
std::cout << "t: " << t;
std::cout << "\ny:";
N_VPrint_Serial(y);
if(check_flag(&flag, "CVode", 1))
{
return 1;
}
else
{
return 0;
}
}
void SimpleCvodeEample::Release()
{
// 15. Get optional outputs.
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// 16. Deallocate memory for solution vector.
// ---------------------------------------------------------------------------
N_VDestroy(y);
// ---------------------------------------------------------------------------
// 17. Free solver memory.
// ---------------------------------------------------------------------------
CVodeFree(&cvode_mem);
// ---------------------------------------------------------------------------
// 18. Free linear solver and matrix memory.
// ---------------------------------------------------------------------------
SUNLinSolFree(LS);
// ---------------------------------------------------------------------------
}
// Simple function that calculates the differential equation.
int SimpleCvodeEample::f_func(realtype t, N_Vector u, N_Vector u_dot, void *user_data) {
// N_VGetArrayPointer returns a pointer to the data in the N_Vector class.
realtype *udata = N_VGetArrayPointer(u); // pointer u vector data
realtype *dudata = N_VGetArrayPointer(u_dot); // pointer to udot vector data
dudata[0] = -101.0 * udata[0] - 100.0 * udata[1];
dudata[1] = udata[0];
return(0);
}
// Jacobian function vector routine.
int SimpleCvodeEample::jtv_func(N_Vector v, N_Vector Jv, realtype t, N_Vector u, N_Vector fu,
void *user_data, N_Vector tmp) {
realtype *udata = N_VGetArrayPointer(u);
realtype *vdata = N_VGetArrayPointer(v);
realtype *Jvdata = N_VGetArrayPointer(Jv);
realtype *fudata = N_VGetArrayPointer(fu);
Jvdata[0] = -101.0 * vdata[0] + -100.0 * vdata[1];
Jvdata[1] = vdata[0] + 0 * vdata[1];
fudata[0] = 0;
fudata[1] = 0;
return(0);
}
// check_flag function is from the cvDiurnals_ky.c example from the CVODE
// package.
/* Check function return value...
opt == 0 means SUNDIALS function allocates memory so check if
returned NULL pointer
opt == 1 means SUNDIALS function returns a flag so check if
flag >= 0
opt == 2 means function allocates memory so check if returned
NULL pointer */
int SimpleCvodeEample::check_flag_func(void *flagvalue, const char *funcname, int opt) {
int *errflag;
/* Check if SUNDIALS function returned NULL pointer - no memory allocated */
if (opt == 0 && flagvalue == NULL) {
fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
funcname);
return(1); }
/* Check if flag < 0 */
else if (opt == 1) {
errflag = (int *) flagvalue;
if (*errflag < 0) {
fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
funcname, *errflag);
return(1); }}
/* Check if function returned NULL pointer - no memory allocated */
else if (opt == 2 && flagvalue == NULL) {
fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
funcname);
return(1); }
return(0);
}
测试入口代码
#include <stdio.h>
#include <iostream>
#include "SimpleCvodeEample.h"
using namespace std;
int main()
{
SimpleCvodeEample test;
if(1 == test.Init())
{
printf("Init error\n");
return 0;
}
int print_steps = 100;
realtype tout;
realtype end_time = 50;
realtype step_length = 0.5;
realtype t = 0;
// loop over output points, call CVode, print results, test for error
for (tout = step_length; tout <= end_time; tout += step_length)
{
if(1 == test.Run(tout))
{
break;
}
}
return 0;
}
输出
t: 0.5
y:-0.6249119745010792
0.6249119551816967
t: 1
y: -0.379070277786671
0.3790702791895351
t: 1.5
y:-0.2299103165625554
0.2299103126197803
t: 2
y:-0.1394434400073003
0.1394434406695035
t: 2.5
y:-0.08454829657003644
0.08454830412509731
t: 3
y:-0.05127930281061206
0.05127931029361821
t: 3.5
y:-0.0311024737657117
0.03110245162704426
t: 4
y:-0.01886629458285826
0.01886628385578891
t: 4.5
y:-0.01144272056713856
0.01144273456266413
t: 5
y:-0.006945105331771737
0.006945035771704891
t: 5.5
y:-0.004224473244729604
0.004224327316890079
t: 6
y:-0.00257470647381695
0.002574668593885272
t: 6.5
y:-0.001566067716433651
0.001566053134294854
t: 7
y:-0.0009501147598849761
0.0009501161651030121
t: 7.5
y:-0.0005788732835267033
0.0005788745302760797
t: 8
y:-0.000359009562134631
0.0003590106839948549
t: 8.5
y:-0.0002332064937038593
0.0002332110070028788
t: 9
y:-0.0001551129067039888
0.0001551148745547571
t: 9.5
y:-0.0001000625783188643
0.0001000624987545105
t: 10
y:-6.12786020735907e-05
6.127161930100425e-05
t: 10.5
y:-3.560524006424445e-05
3.553294750139934e-05
t: 11
y:-1.912367516965566e-05
1.90409034985556e-05
t: 11.5
y:-1.034810017343243e-05
1.033000446221329e-05
t: 12
y:-7.662977402627386e-06
7.785951319782404e-06
t: 12.5
y:-7.831120566794994e-06
7.99844861699368e-06
t: 13
y:-7.984808393223872e-06
7.988672221302979e-06
t: 13.5
y:-9.419213635552393e-06
9.413215116337815e-06
t: 14
y:-7.934132298185239e-06
7.934582029360176e-06
t: 14.5
y:-6.333997220997139e-06
6.343042253911147e-06
t: 15
y:-4.130913358365324e-06
4.136779805322092e-06
t: 15.5
y:-1.614847840609949e-06
1.618436063094846e-06
t: 16
y:3.135490166074421e-07
-3.09500983956612e-07
t: 16.5
y:8.434891282761371e-07
-8.443113781444554e-07
t: 17
y:8.001513512591652e-07
-8.026519065980999e-07
t: 17.5
y:3.345496151746862e-07
-3.346430265247854e-07
t: 18
y:-7.666923683566892e-07
7.669878252988782e-07
t: 18.5
y:-1.309930352692327e-06
1.31014411485376e-06
t: 19
y:-1.450518602771783e-06
1.450627246789752e-06
t: 19.5
y:-1.284636103687557e-06
1.284742638624243e-06
t: 20
y:-7.075641747285998e-07
7.075229907795396e-07
t: 20.5
y:4.850466910265154e-08
-4.861626913565485e-08
t: 21
y:4.511555942860418e-07
-4.511763079351868e-07
t: 21.5
y:1.659404448308286e-06
-1.659588122931253e-06
t: 22
y:2.893645364426666e-06
-2.894000402230771e-06
t: 22.5
y:3.442550364388173e-06
-3.442852240580473e-06
t: 23
y:3.463250744704413e-06
-3.46334900436318e-06
t: 23.5
y:4.54118590846849e-06
-4.541509623084763e-06
t: 24
y:5.214854199393966e-06
-5.21530899991576e-06
t: 24.5
y:4.725129202436564e-06
-4.725337336614224e-06
t: 25
y:4.31237501592474e-06
-4.312442925367661e-06
t: 25.5
y:4.243449496053212e-06
-4.2435884707076e-06
t: 26
y:3.823008111093036e-06
-3.823158931093899e-06
t: 26.5
y:2.815965134658726e-06
-2.815970585340204e-06
t: 27
y:2.136042066505748e-06
-2.136050040487766e-06
t: 27.5
y:1.456611463632213e-06
-1.456621988608491e-06
t: 28
y:8.786691385639032e-07
-8.786775331306455e-07
t: 28.5
y:4.656137682549884e-07
-4.658737829217555e-07
t: 29
y:1.303755732246936e-07
-1.321699331737769e-07
t: 29.5
y:-1.404659476286379e-08
9.194509644209812e-09
t: 30
y:5.906672705554527e-08
-6.923315876038694e-08
t: 30.5
y:2.052038861745693e-07
-2.195608081560045e-07
t: 31
y:3.330027490467403e-07
-3.483117619200677e-07
t: 31.5
y:4.732636112468342e-07
-4.878119016440941e-07
t: 32
y:5.787515782694418e-07
-5.905740338896326e-07
t: 32.5
y:5.946973676299159e-07
-6.016510299291738e-07
t: 33
y:4.711642902749203e-07
-4.711405009677789e-07
t: 33.5
y:5.452992207773848e-07
-5.40691745696408e-07
t: 34
y:5.94631109831263e-07
-5.862637006687925e-07
t: 34.5
y:6.057233973632193e-07
-5.953355194394757e-07
t: 35
y:5.643866168910427e-07
-5.547834067724029e-07
t: 35.5
y:4.556783955236465e-07
-4.508806186409208e-07
t: 36
y:3.290980931783813e-07
-3.294852942750648e-07
t: 36.5
y:2.440505478124803e-07
-2.456924677087564e-07
t: 37
y:1.570242588679222e-07
-1.596767225193197e-07
t: 37.5
y:7.835240245731255e-08
-8.141430037570089e-08
t: 38
y:2.035569275301578e-08
-2.280847842172226e-08
t: 38.5
y:-2.657618012845719e-09
2.310430723919075e-09
t: 39
y:-5.718491990723776e-08
5.723974974957211e-08
t: 39.5
y:-1.062677963708564e-07
1.064632767693695e-07
t: 40
y:-1.391367786493883e-07
1.394419566773716e-07
t: 40.5
y:-1.50661988710465e-07
1.509924261927111e-07
t: 41
y:-1.353183232947853e-07
1.35524698823497e-07
t: 41.5
y:-1.119209274682172e-07
1.120307715351344e-07
t: 42
y:-1.229066245901732e-07
1.23358701155367e-07
t: 42.5
y:-1.324338000230351e-07
1.332696823365114e-07
t: 43
y:-1.399817184904001e-07
1.411998956872731e-07
t: 43.5
y:-1.44672746220611e-07
1.462176175577227e-07
t: 44
y:-1.452723509467553e-07
1.470232200392959e-07
t: 44.5
y:-1.40189101906666e-07
1.419491709647932e-07
t: 45
y:-1.27474669842921e-07
1.289600339083802e-07
t: 45.5
y:-1.048238270028436e-07
1.056524681855875e-07
t: 46
y:-8.16327463180977e-08
8.184339979968385e-08
t: 46.5
y:-8.890388157078401e-08
8.986945764476941e-08
t: 47
y:-9.653315694094676e-08
9.827263446907575e-08
t: 47.5
y:-1.041000991891444e-07
1.065999869389798e-07
t: 48
y:-1.111842350759358e-07
1.143985717208587e-07
t: 48.5
y:-1.173650913618795e-07
1.212154454810894e-07
t: 49
y:-1.222221948075341e-07
1.26597664886049e-07
t: 49.5
y:-1.253350721734585e-07
1.300922866021146e-07
t: 50
y:-1.262832502202113e-07
1.312463672956633e-07
类似资料：