• CVODE入门

    1、安装

    cvode是sundials软件包的一部分。在windows+visualstudio下可以用vcpkg安装使用。

    git clone https://github.com/microsoft/vcpkg

    然后运行目录中的Bootstrap_vcpkg.bat。

    然后进入vcpkg目录中，vcpkg install sundials:x64-windows，这样安装64位dll。

    然后管理员身份运行cmd，vcpkg integrate install。

    2 使用

    新建项目后，在项目属性中可以看到vcpkg选项，证明安装好了。

    3 示例程序

    程序的主要步骤

    1、初始化多线程环境（可选）
    2、创建sundials上下文对象
    SUNContext_Create()
    3、设置问题的大小？？？
    4、设置初值向量
    如果已有数据存储在ydata中，可以调用y0 = N_VMake_***(..., ydata)
    否则，调用N_VNew_***(...)创建新的向量
    然后用ydata = N_VGetArrayPointer(y0)设置单独数据
    5、创建CVODE对象
    CVodeCreate()返回一个内存结构的指针
    6、初始化CVODE求解器
    CVodeInit()，定义方程组信息，分配内部的内存，初始化。
    7、设置积分误差
    CVodeSStolerances() or CVodeSVtolerances() 设置相对误差或者绝对误差
    8、创建矩阵对象

    matrix_var = SUNDenseMatrix(x, x, context_var)；

    9、创建求解器

    ls = SUNxxx_xxx(y, matrix_var, context_var)

    10、把矩阵加入求解器

    CVodeSetLinearSolver(mem, ls, matrix_var);

    11、设置用户定义的雅克比矩阵

    CVodeSetJacFn(mem, Jac), jac是用户定义函数
    13、时间步长增加
    ier = CVode(cvode_mem, tout, yout, tret itask)
    14、释放解向量内存
    15、释放求解器内存
    CVodeFree()
    16、释放线性求解器和矩阵的内存
    SUNLinSolFree()
    SUNMatDestroy()
    17、释放SUNContext对象
     

    其他

    Ith和IJth是自定义函数，方便函数形式符合数学公式的形式，即下标序号从1开始，而不是c语言的从0开始。
    函数f是ode表达式的直接转写,怎么写详细见https://sundials.readthedocs.io/en/latest/cvode/Usage/index.html#user-supplied-functions。函数g定义了两个数学函数g0和g1，用于求特定值时的参数（比如y1=1,此时求其他的y，可选)。Jac为Jacobian矩阵，该矩阵初始值为0，所以可以只设置非0值。

    CVodePrintAllStats输出所有统计信息

    1. /* -----------------------------------------------------------------
    2.  * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and
    3.  *                Radu Serban @ LLNL
    4.  * -----------------------------------------------------------------
    5.  * SUNDIALS Copyright Start
    6.  * Copyright (c) 2002-2022, Lawrence Livermore National Security
    7.  * and Southern Methodist University.
    8.  * All rights reserved.
    9.  *
    10.  * See the top-level LICENSE and NOTICE files for details.
    11.  *
    12.  * SPDX-License-Identifier: BSD-3-Clause
    13.  * SUNDIALS Copyright End
    14.  * -----------------------------------------------------------------
    15.  * Example problem:
    16.  *
    17.  * The following is a simple example problem, with the coding
    18.  * needed for its solution by CVODE. The problem is from
    19.  * chemical kinetics, and consists of the following three rate
    20.  * equations:
    21.  *    dy1/dt = -.04*y1 + 1.e4*y2*y3
    22.  *    dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2
    23.  *    dy3/dt = 3.e7*(y2)^2
    24.  * on the interval from t = 0.0 to t = 4.e10, with initial
    25.  * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.
    26.  * While integrating the system, we also use the rootfinding
    27.  * feature to find the points at which y1 = 1e-4 or at which
    28.  * y3 = 0.01. This program solves the problem with the BDF method,
    29.  * Newton iteration with the dense linear solver, and a
    30.  * user-supplied Jacobian routine.
    31.  * It uses a scalar relative tolerance and a vector absolute
    32.  * tolerance. Output is printed in decades from t = .4 to t = 4.e10.
    33.  * Run statistics (optional outputs) are printed at the end.
    34.  * -----------------------------------------------------------------*/
    35.  
    36. #include 
    37. //#include 
    38.  
    39. #include                /* prototypes for CVODE fcts., consts.  */
    40. #include     /* access to serial N_Vector            */
    41. #include  /* access to dense SUNMatrix            */
    42. #include  /* access to dense SUNLinearSolver      */
    43.  
    44. #if defined(SUNDIALS_EXTENDED_PRECISION)
    45. #define GSYM "Lg"
    46. #define ESYM "Le"
    47. #define FSYM "Lf"
    48. #else
    49. #define GSYM "g"
    50. #define ESYM "e"
    51. #define FSYM "f"
    52. #endif
    53.  
    54. /* User-defined vector and matrix accessor macros: Ith, IJth */
    55.  
    56. /* These macros are defined in order to write code which exactly matches
    57.    the mathematical problem description given above.
    59.    Ith(v,i) references the ith component of the vector v, where i is in
    60.    the range [1..NEQ] and NEQ is defined below. The Ith macro is defined
    61.    using the N_VIth macro in nvector.h. N_VIth numbers the components of
    62.    a vector starting from 0.
    64.    IJth(A,i,j) references the (i,j)th element of the dense matrix A, where
    65.    i and j are in the range [1..NEQ]. The IJth macro is defined using the
    66.    SM_ELEMENT_D macro. SM_ELEMENT_D numbers rows and columns of
    67.    a dense matrix starting from 0. */
    68.  
    69. #define Ith(v,i)    NV_Ith_S(v,i-1)         /* i-th vector component i=1..NEQ */
    70. #define IJth(A,i,j) SM_ELEMENT_D(A,i-1,j-1) /* (i,j)-th matrix component i,j=1..NEQ */
    71.  
    72. /* Problem Constants */
    73.  
    74. #define NEQ   3                /* number of equations  */
    75. #define Y1    RCONST(1.0)      /* initial y components */
    76. #define Y2    RCONST(0.0)
    77. #define Y3    RCONST(0.0)
    78. #define RTOL  RCONST(1.0e-4)   /* scalar relative tolerance            */
    79. #define ATOL1 RCONST(1.0e-8)   /* vector absolute tolerance components */
    80. #define ATOL2 RCONST(1.0e-14)
    81. #define ATOL3 RCONST(1.0e-6)
    82. #define T0    RCONST(0.0)      /* initial time           */
    83. #define T1    RCONST(0.4)      /* first output time      */
    84. #define TMULT RCONST(10.0)     /* output time factor     */
    85. #define NOUT  12               /* number of output times */
    86.  
    87. #define ZERO  RCONST(0.0)
    88.  
    89. /* Functions Called by the Solver */
    90.  
    91. static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
    92.  
    93. static int g(realtype t, N_Vector y, realtype *gout, void *user_data);
    94.  
    95. static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
    96.                void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
    97.  
    98. /* Private functions to output results */
    99.  
    100. static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3);
    101. static void PrintRootInfo(int root_f1, int root_f2);
    102.  
    103. /* Private function to check function return values */
    104.  
    105. static int check_retval(void *returnvalue, const char *funcname, int opt);
    106.  
    107. /* Private function to check computed solution */
    108.  
    109. static int check_ans(N_Vector y, realtype t, realtype rtol, N_Vector atol);
    110.  
    111.  
    112. /*
    113.  *-------------------------------
    114.  * Main Program
    115.  *-------------------------------
    116.  */
    117.  
    118. int main()
    119. {
    120.   SUNContext sunctx;
    121.   realtype t, tout;
    122.   N_Vector y;
    123.   N_Vector abstol;
    124.   SUNMatrix A;
    125.   SUNLinearSolver LS;
    126.   void *cvode_mem;
    127.   int retval, iout;
    128.   int retvalr;
    129.   int rootsfound[2];
    130.   FILE* FID;
    131.  
    132.   y = NULL;
    133.   abstol = NULL;
    134.   A = NULL;
    135.   LS = NULL;
    136.   cvode_mem = NULL;
    137.  
    138.   /* Create the SUNDIALS context */
    139.   retval = SUNContext_Create(NULL, &sunctx);
    140.   if (check_retval(&retval, "SUNContext_Create", 1)) return(1);
    141.  
    142.   /* Initial conditions */
    143.   y = N_VNew_Serial(NEQ, sunctx);
    144.   if (check_retval((void *)y, "N_VNew_Serial", 0)) return(1);
    145.  
    146.   /* Initialize y */
    147.   Ith(y,1) = Y1;
    148.   Ith(y,2) = Y2;
    149.   Ith(y,3) = Y3;
    150.  
    151.   /* Set the vector absolute tolerance */
    152.   abstol = N_VNew_Serial(NEQ, sunctx);
    153.   if (check_retval((void *)abstol, "N_VNew_Serial", 0)) return(1);
    154.  
    155.   Ith(abstol,1) = ATOL1;
    156.   Ith(abstol,2) = ATOL2;
    157.   Ith(abstol,3) = ATOL3;
    158.  
    159.   /* Call CVodeCreate to create the solver memory and specify the
    160.    * Backward Differentiation Formula */
    161.   cvode_mem = CVodeCreate(CV_BDF, sunctx);
    162.   if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
    163.  
    164.   /* Call CVodeInit to initialize the integrator memory and specify the
    165.    * user's right hand side function in y'=f(t,y), the initial time T0, and
    166.    * the initial dependent variable vector y. */
    167.   retval = CVodeInit(cvode_mem, f, T0, y);
    168.   if (check_retval(&retval, "CVodeInit", 1)) return(1);
    169.  
    170.   /* Call CVodeSVtolerances to specify the scalar relative tolerance
    171.    * and vector absolute tolerances */
    172.   retval = CVodeSVtolerances(cvode_mem, RTOL, abstol);
    173.   if (check_retval(&retval, "CVodeSVtolerances", 1)) return(1);
    174.  
    175.   /* Call CVodeRootInit to specify the root function g with 2 components */
    176.   //retval = CVodeRootInit(cvode_mem, 2, g);
    177.   //if (check_retval(&retval, "CVodeRootInit", 1)) return(1);
    178.  
    179.   /* Create dense SUNMatrix for use in linear solves */
    180.   A = SUNDenseMatrix(3, 3, sunctx);
    181.   if (check_retval((void *)A, "SUNDenseMatrix", 0)) return(1);
    182.  
    183.   /* Create dense SUNLinearSolver object for use by CVode */
    184.   LS = SUNLinSol_Dense(y, A, sunctx);
    185.   if (check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1);
    186.  
    187.   /* Attach the matrix and linear solver */
    188.   retval = CVodeSetLinearSolver(cvode_mem, LS, A);
    189.   if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);
    190.  
    191.   /* Set the user-supplied Jacobian routine Jac */
    192.   retval = CVodeSetJacFn(cvode_mem, Jac);
    193.   if (check_retval(&retval, "CVodeSetJacFn", 1)) return(1);
    194.  
    195.   /* In loop, call CVode, print results, and test for error.
    196.      Break out of loop when NOUT preset output times have been reached.  */
    197.   printf(" \n3-species kinetics problem\n\n");
    198.  
    199.   FILE* vstream;
    200.  
    201.   /* Open file for printing statistics */
    202.   //FID = fopen_s(vstream, "cvRoberts_dns_stats.csv", "w");
    203.   errno_t ee = fopen_s(&vstream, "cvRoberts_dns_stats.csv", "w");
    204.  
    205.   iout = 0;  tout = T1;
    206.   while(1) {
    207.     retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
    208.     PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));
    209.  
    210.     if (retval == CV_ROOT_RETURN) {
    211.       retvalr = CVodeGetRootInfo(cvode_mem, rootsfound);
    212.       if (check_retval(&retvalr, "CVodeGetRootInfo", 1)) return(1);
    213.       PrintRootInfo(rootsfound[0],rootsfound[1]);
    214.     }
    215.  
    216.     if (check_retval(&retval, "CVode", 1)) break;
    217.     if (retval == CV_SUCCESS) {
    218.       iout++;
    219.       tout *= TMULT;
    220.     }
    221.  
    222.     retval = CVodePrintAllStats(cvode_mem, vstream, SUN_OUTPUTFORMAT_CSV);
    223.  
    224.     if (iout == NOUT) break;
    225.   }
    226.   fclose(vstream);
    227.  
    228.   /* Print final statistics to the screen */
    229.   printf("\nFinal Statistics:\n");
    230.   retval = CVodePrintAllStats(cvode_mem, stdout, SUN_OUTPUTFORMAT_TABLE);
    231.  
    232.   /* check the solution error */
    233.   retval = check_ans(y, t, RTOL, abstol);
    234.  
    235.   /* Free memory */
    236.   N_VDestroy(y);                            /* Free y vector */
    237.   N_VDestroy(abstol);                       /* Free abstol vector */
    238.   CVodeFree(&cvode_mem);                    /* Free CVODE memory */
    239.   SUNLinSolFree(LS);                        /* Free the linear solver memory */
    240.   SUNMatDestroy(A);                         /* Free the matrix memory */
    241.   SUNContext_Free(&sunctx);                 /* Free the SUNDIALS context */
    242.  
    243.   return(retval);
    244. }
    245.  
    246.  
    247. /*
    248.  *-------------------------------
    249.  * Functions called by the solver
    250.  *-------------------------------
    251.  */
    252.  
    253. /*
    254.  * f routine. Compute function f(t,y).
    255.  */
    256.  
    257. static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
    258. {
    259.   realtype y1, y2, y3, yd1, yd3;
    260.  
    261.   y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
    262.  
    263.   yd1 = Ith(ydot,1) = RCONST(-0.04)*y1 + RCONST(1.0e4)*y2*y3;
    264.   yd3 = Ith(ydot,3) = RCONST(3.0e7)*y2*y2;
    265.         Ith(ydot,2) = -yd1 - yd3;
    266.  
    267.   return(0);
    268. }
    269.  
    270. /*
    271.  * g routine. Compute functions g_i(t,y) for i = 0,1.
    272.  */
    273.  
    274. static int g(realtype t, N_Vector y, realtype *gout, void *user_data)
    275. {
    276.   realtype y1, y3;
    277.  
    278.   y1 = Ith(y,1); y3 = Ith(y,3);
    279.   gout[0] = y1 - RCONST(0.0001);
    280.   gout[1] = y3 - RCONST(0.01);
    281.  
    282.   return(0);
    283. }
    284.  
    285. /*
    286.  * Jacobian routine. Compute J(t,y) = df/dy. *
    287.  */
    288.  
    289. static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
    290.                void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
    291. {
    292.   realtype y2, y3;
    293.  
    294.   y2 = Ith(y,2); y3 = Ith(y,3);
    295.  
    296.   IJth(J,1,1) = RCONST(-0.04);
    297.   IJth(J,1,2) = RCONST(1.0e4)*y3;
    298.   IJth(J,1,3) = RCONST(1.0e4)*y2;
    299.  
    300.   IJth(J,2,1) = RCONST(0.04);
    301.   IJth(J,2,2) = RCONST(-1.0e4)*y3-RCONST(6.0e7)*y2;
    302.   IJth(J,2,3) = RCONST(-1.0e4)*y2;
    303.  
    304.   IJth(J,3,1) = ZERO;
    305.   IJth(J,3,2) = RCONST(6.0e7)*y2;
    306.   IJth(J,3,3) = ZERO;
    307.  
    308.   return(0);
    309. }
    310.  
    311. /*
    312.  *-------------------------------
    313.  * Private helper functions
    314.  *-------------------------------
    315.  */
    316.  
    317. static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3)
    318. {
    319. #if defined(SUNDIALS_EXTENDED_PRECISION)
    320.   printf("At t = %0.4Le      y =%14.6Le  %14.6Le  %14.6Le\n", t, y1, y2, y3);
    321. #elif defined(SUNDIALS_DOUBLE_PRECISION)
    322.   printf("At t = %0.4e      y =%14.6e  %14.6e  %14.6e\n", t, y1, y2, y3);
    323. #else
    324.   printf("At t = %0.4e      y =%14.6e  %14.6e  %14.6e\n", t, y1, y2, y3);
    325. #endif
    326.  
    327.   return;
    328. }
    329.  
    330. static void PrintRootInfo(int root_f1, int root_f2)
    331. {
    332.   printf("    rootsfound[] = %3d %3d\n", root_f1, root_f2);
    333.  
    334.   return;
    335. }
    336.  
    337. /*
    338.  * Check function return value...
    339.  *   opt == 0 means SUNDIALS function allocates memory so check if
    340.  *            returned NULL pointer
    341.  *   opt == 1 means SUNDIALS function returns an integer value so check if
    342.  *            retval < 0
    343.  *   opt == 2 means function allocates memory so check if returned
    344.  *            NULL pointer
    345.  */
    346.  
    347. static int check_retval(void *returnvalue, const char *funcname, int opt)
    348. {
    349.   int *retval;
    350.  
    351.   /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
    352.   if (opt == 0 && returnvalue == NULL) {
    353.     fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
    354.             funcname);
    355.     return(1); }
    356.  
    357.   /* Check if retval < 0 */
    358.   else if (opt == 1) {
    359.     retval = (int *) returnvalue;
    360.     if (*retval < 0) {
    361.       fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
    362.               funcname, *retval);
    363.       return(1); }}
    364.  
    365.   /* Check if function returned NULL pointer - no memory allocated */
    366.   else if (opt == 2 && returnvalue == NULL) {
    367.     fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
    368.             funcname);
    369.     return(1); }
    370.  
    371.   return(0);
    372. }
    373.  
    374. /* compare the solution at the final time 4e10s to a reference solution computed
    375.    using a relative tolerance of 1e-8 and absoltue tolerance of 1e-14 */
    376. static int check_ans(N_Vector y, realtype t, realtype rtol, N_Vector atol)
    377. {
    378.   int      passfail=0;        /* answer pass (0) or fail (1) flag */
    379.   N_Vector ref;               /* reference solution vector        */
    380.   N_Vector ewt;               /* error weight vector              */
    381.   realtype err;               /* wrms error                       */
    382.   realtype ONE=RCONST(1.0);
    383.  
    384.   /* create reference solution and error weight vectors */
    385.   ref = N_VClone(y);
    386.   ewt = N_VClone(y);
    387.  
    388.   /* set the reference solution data */
    389.   NV_Ith_S(ref,0) = RCONST(5.2083495894337328e-08);
    390.   NV_Ith_S(ref,1) = RCONST(2.0833399429795671e-13);
    391.   NV_Ith_S(ref,2) = RCONST(9.9999994791629776e-01);
    392.  
    393.   /* compute the error weight vector, loosen atol */
    394.   N_VAbs(ref, ewt);
    395.   N_VLinearSum(rtol, ewt, RCONST(10.0), atol, ewt);
    396.   if (N_VMin(ewt) <= ZERO) {
    397.     fprintf(stderr, "\nSUNDIALS_ERROR: check_ans failed - ewt <= 0\n\n");
    398.     return(-1);
    399.   }
    400.   N_VInv(ewt, ewt);
    401.  
    402.   /* compute the solution error */
    403.   N_VLinearSum(ONE, y, -ONE, ref, ref);
    404.   err = N_VWrmsNorm(ref, ewt);
    405.  
    406.   /* is the solution within the tolerances? */
    407.   passfail = (err < ONE) ? 0 : 1;
    408.  
    409.   if (passfail) {
    410.     fprintf(stdout, "\nSUNDIALS_WARNING: check_ans error=%"GSYM"\n\n", err);
    411.   }
    412.  
    413.   /* Free vectors */
    414.   N_VDestroy(ref);
    415.   N_VDestroy(ewt);
    416.  
    417.   return(passfail);
    418. }

  • 相关阅读:
    一键自动化博客发布工具,用过的人都说好(腾讯云篇)
    2023-09-14 mysql-Item_subselect-分析
    【C++进阶】：哈希
    数据治理：误区梳理篇
    SpringBoot基础入门
    接口自动化测试小结
    Qt扫盲-QPen 理论使用总结
    基于flowable的upp(统一流程平台)运行性能优化（2）
    Spread.NET 16.0[Winform|WPF|ASP.NET]
    某电商网站的数据库设计（3）
  • 原文地址：https://blog.csdn.net/novanova2009/article/details/126798087