API tools faq
paste
Advertisement
Guest User

Untitled

a guest
May 25th, 2019
133
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 20.50 KB | None | 0 0
raw download clone embed print report
  1. /* -----------------------------------------------------------------
  2. * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, and
  3. * Radu Serban @ LLNL
  4. * -----------------------------------------------------------------
  5. * LLNS Copyright Start
  6. * Copyright (c) 2017, Lawrence Livermore National Security
  7. * This work was performed under the auspices of the U.S. Department
  8. * of Energy by Lawrence Livermore National Laboratory in part under
  9. * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344.
  10. * Produced at the Lawrence Livermore National Laboratory.
  11. * All rights reserved.
  12. * For details, see the LICENSE file.
  13. * LLNS Copyright End
  14. * -----------------------------------------------------------------
  15. * Example problem:
  16. *
  17. * The following is a simple example problem, with the coding
  18. * needed for its solution by CVODES for Forward Sensitivity
  19. * Analysis. The problem is from chemical kinetics, and consists
  20. * of the following three rate equations:
  21. * dy1/dt = -p1*y1 + p2*y2*y3
  22. * dy2/dt = p1*y1 - p2*y2*y3 - p3*(y2)^2
  23. * dy3/dt = p3*(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 reaction rates are: p1=0.04,
  26. * p2=1e4, and p3=3e7. The problem is stiff.
  27. * This program solves the problem with the BDF method, Newton
  28. * iteration with the DENSE linear solver, and a
  29. * user-supplied Jacobian routine.
  30. * It uses a scalar relative tolerance and a vector absolute
  31. * tolerance.
  32. * Output is printed in decades from t = .4 to t = 4.e10.
  33. * Run statistics (optional outputs) are printed at the end.
  34. *
  35. * Optionally, CVODES can compute sensitivities with respect to the
  36. * problem parameters p1, p2, and p3.
  37. * The sensitivity right hand side is given analytically through the
  38. * user routine fS (of type SensRhs1Fn).
  39. * Any of three sensitivity methods (SIMULTANEOUS, STAGGERED, and
  40. * STAGGERED1) can be used and sensitivities may be included in the
  41. * error test or not (error control set on SUNTRUE or SUNFALSE,
  42. * respectively).
  43. *
  44. * Execution:
  45. *
  46. * If no sensitivities are desired:
  47. * % cvsRoberts_FSA_dns -nosensi
  48. * If sensitivities are to be computed:
  49. * % cvsRoberts_FSA_dns -sensi sensi_meth err_con
  50. * where sensi_meth is one of {sim, stg, stg1} and err_con is one of
  51. * {t, f}.
  52. * -----------------------------------------------------------------*/
  53.  
  54. #include <stdio.h>
  55. #include <stdlib.h>
  56. #include <string.h>
  57.  
  58. #include <cvodes/cvodes.h> /* prototypes for CVODE fcts., consts. */
  59. #include <nvector/nvector_serial.h> /* access to serial N_Vector */
  60. #include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
  61. #include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
  62. #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
  63. #include <sundials/sundials_math.h> /* definition of ABS */
  64.  
  65.  
  66. #include <iostream>
  67. #include <vector>
  68.  
  69. /* Accessor macros */
  70.  
  71. #define Ith(v,i) NV_Ith_S(v,i-1) /* i-th vector component i=1..NEQ */
  72. #define IJth(A,i,j) SM_ELEMENT_D(A,i-1,j-1) /* (i,j)-th matrix component i,j=1..NEQ */
  73.  
  74. /* Problem Constants */
  75.  
  76. #define NEQ 3 /* number of equations */
  77. #define Y1 RCONST(1.0) /* initial y components */
  78. #define Y2 RCONST(0.0)
  79. #define Y3 RCONST(0.0)
  80. #define RTOL RCONST(1e-4) /* scalar relative tolerance */
  81. #define ATOL1 RCONST(1e-8) /* vector absolute tolerance components */
  82. #define ATOL2 RCONST(1e-14)
  83. #define ATOL3 RCONST(1e-6)
  84. #define T0 RCONST(0.0) /* initial time */
  85. #define T1 RCONST(0.4) /* first output time */
  86. #define TMULT RCONST(10.0) /* output time factor */
  87. #define NOUT 12 /* number of output times */
  88.  
  89. #define NP 3 /* number of problem parameters */
  90. #define NS 3 /* number of sensitivities computed */
  91.  
  92. #define ZERO RCONST(0.0)
  93.  
  94. /* Type : UserData */
  95.  
  96. typedef struct {
  97. // realtype p[3]; /* problem parameters */
  98. std::vector<double> p;
  99. } *UserData;
  100.  
  101. /* Prototypes of functions by CVODES */
  102.  
  103. static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
  104.  
  105. static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
  106. void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
  107.  
  108. static int fS(int Ns, realtype t, N_Vector y, N_Vector ydot,
  109. int iS, N_Vector yS, N_Vector ySdot,
  110. void *user_data, N_Vector tmp1, N_Vector tmp2);
  111.  
  112. static int ewt(N_Vector y, N_Vector w, void *user_data);
  113.  
  114. /* Prototypes of private functions */
  115.  
  116. static void ProcessArgs(int argc, char *argv[],
  117. booleantype *sensi, int *sensi_meth,
  118. booleantype *err_con);
  119. static void WrongArgs(char *name);
  120. static void PrintOutput(void *cvode_mem, realtype t, N_Vector u);
  121. static void PrintOutputS(N_Vector *uS);
  122. static void PrintFinalStats(void *cvode_mem, booleantype sensi);
  123. static int check_retval(void *returnvalue, const char *funcname, int opt);
  124.  
  125. /*
  126. *--------------------------------------------------------------------
  127. * MAIN PROGRAM
  128. *--------------------------------------------------------------------
  129. */
  130.  
  131. int main(int argc, char *argv[])
  132. {
  133. SUNMatrix A;
  134. SUNLinearSolver LS;
  135. void *cvode_mem;
  136. UserData data;
  137. realtype t, tout;
  138. N_Vector y;
  139. int iout, retval;
  140.  
  141. realtype pbar[NS];
  142. int is;
  143. N_Vector *yS;
  144. booleantype sensi, err_con;
  145. int sensi_meth;
  146.  
  147. cvode_mem = NULL;
  148. data = NULL;
  149. y = NULL;
  150. yS = NULL;
  151. A = NULL;
  152. LS = NULL;
  153.  
  154. /* Process arguments */
  155. ProcessArgs(argc, argv, &sensi, &sensi_meth, &err_con);
  156.  
  157. /* User data structure */
  158. data = (UserData) malloc(sizeof *data);
  159. if (check_retval((void *)data, "malloc", 2)) return(1);
  160. // data->p[0] = RCONST(0.04);
  161. // data->p[1] = RCONST(1.0e4);
  162. // data->p[2] = RCONST(3.0e7);
  163.  
  164. data->p.push_back(RCONST(0.04));
  165. data->p.push_back(RCONST(1.0e4));
  166. data->p.push_back(RCONST(3.0e7));
  167.  
  168. /* Initial conditions */
  169. y = N_VNew_Serial(NEQ);
  170. if (check_retval((void *)y, "N_VNew_Serial", 0)) return(1);
  171.  
  172. Ith(y,1) = Y1;
  173. Ith(y,2) = Y2;
  174. Ith(y,3) = Y3;
  175.  
  176. /* Create CVODES object */
  177. cvode_mem = CVodeCreate(CV_BDF);
  178. if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
  179.  
  180. /* Allocate space for CVODES */
  181. retval = CVodeInit(cvode_mem, f, T0, y);
  182. if (check_retval(&retval, "CVodeInit", 1)) return(1);
  183.  
  184. /* Use private function to compute error weights */
  185. retval = CVodeWFtolerances(cvode_mem, ewt);
  186. if (check_retval(&retval, "CVodeSetEwtFn", 1)) return(1);
  187.  
  188. /* Attach user data */
  189. retval = CVodeSetUserData(cvode_mem, data);
  190. if (check_retval(&retval, "CVodeSetUserData", 1)) return(1);
  191.  
  192. /* Create dense SUNMatrix */
  193. A = SUNDenseMatrix(NEQ, NEQ);
  194. if (check_retval((void *)A, "SUNDenseMatrix", 0)) return(1);
  195.  
  196. /* Create dense SUNLinearSolver */
  197. LS = SUNLinSol_Dense(y, A);
  198. if (check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1);
  199.  
  200. /* Attach the matrix and linear solver */
  201. retval = CVodeSetLinearSolver(cvode_mem, LS, A);
  202. if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);
  203.  
  204. /* Set the user-supplied Jacobian routine Jac */
  205. /*retval = CVodeSetJacFn(cvode_mem, Jac);
  206. if (check_retval(&retval, "CVodeSetJacFn", 1)) return(1);
  207. */
  208. printf("\n3-species chemical kinetics problem\n");
  209.  
  210. /* Sensitivity-related settings */
  211. if (sensi) {
  212.  
  213. /* Set parameter scaling factor */
  214. pbar[0] = data->p.at(0);
  215. pbar[1] = data->p.at(1);
  216. pbar[2] = data->p.at(2);
  217.  
  218. /* Set sensitivity initial conditions */
  219. yS = N_VCloneVectorArray(NS, y);
  220. if (check_retval((void *)yS, "N_VCloneVectorArray", 0)) return(1);
  221. for (is=0;is<NS;is++) N_VConst(ZERO, yS[is]);
  222.  
  223. /* Call CVodeSensInit1 to activate forward sensitivity computations
  224. and allocate internal memory for COVEDS related to sensitivity
  225. calculations. Computes the right-hand sides of the sensitivity
  226. ODE, one at a time */
  227. retval = CVodeSensInit1(cvode_mem, NS, sensi_meth, NULL, yS); // fS is third argument
  228. if(check_retval(&retval, "CVodeSensInit", 1)) return(1);
  229.  
  230.  
  231. /* Call CVodeSensEEtolerances to estimate tolerances for sensitivity
  232. variables based on the rolerances supplied for states variables and
  233. the scaling factor pbar */
  234. retval = CVodeSensEEtolerances(cvode_mem);
  235. if(check_retval(&retval, "CVodeSensEEtolerances", 1)) return(1);
  236.  
  237. /* Set sensitivity analysis optional inputs */
  238. /* Call CVodeSetSensErrCon to specify the error control strategy for
  239. sensitivity variables */
  240. retval = CVodeSetSensErrCon(cvode_mem, err_con);
  241. if (check_retval(&retval, "CVodeSetSensErrCon", 1)) return(1);
  242.  
  243. /* Call CVodeSetSensParams to specify problem parameter information for
  244. sensitivity calculations */
  245. /*
  246. int parameter_len = 3;
  247. realtype *param_array;
  248. param_array = (realtype*)malloc(parameter_len * sizeof(realtype));
  249. param_array[0] = RCONST(0.04);
  250. param_array[1] = RCONST(1.0e4);
  251. param_array[2] = RCONST(3.0e7);
  252. */
  253. retval = CVodeSetSensParams(cvode_mem, data->p[0], pbar, NULL); // data->p
  254. if (check_retval(&retval, "CVodeSetSensParams", 1)) return(1);
  255.  
  256. printf("Sensitivity: YES ");
  257. if(sensi_meth == CV_SIMULTANEOUS)
  258. printf("( SIMULTANEOUS +");
  259. else
  260. if(sensi_meth == CV_STAGGERED) printf("( STAGGERED +");
  261. else printf("( STAGGERED1 +");
  262. if(err_con) printf(" FULL ERROR CONTROL )");
  263. else printf(" PARTIAL ERROR CONTROL )");
  264.  
  265. } else {
  266.  
  267. printf("Sensitivity: NO ");
  268.  
  269. }
  270.  
  271. /* In loop over output points, call CVode, print results, test for error */
  272.  
  273. printf("\n\n");
  274. printf("===========================================");
  275. printf("============================\n");
  276. printf(" T Q H NST y1");
  277. printf(" y2 y3 \n");
  278. printf("===========================================");
  279. printf("============================\n");
  280.  
  281. for (iout=1, tout=T1; iout <= NOUT; iout++, tout *= TMULT) {
  282.  
  283. retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
  284. if (check_retval(&retval, "CVode", 1)) break;
  285.  
  286. PrintOutput(cvode_mem, t, y);
  287.  
  288. /* Call CVodeGetSens to get the sensitivity solution vector after a
  289. successful return from CVode */
  290. if (sensi) {
  291. retval = CVodeGetSens(cvode_mem, &t, yS);
  292. if (check_retval(&retval, "CVodeGetSens", 1)) break;
  293. PrintOutputS(yS);
  294. }
  295. printf("-----------------------------------------");
  296. printf("------------------------------\n");
  297.  
  298. }
  299.  
  300. /* Print final statistics */
  301. PrintFinalStats(cvode_mem, sensi);
  302.  
  303. /* Free memory */
  304.  
  305. N_VDestroy(y); /* Free y vector */
  306. if (sensi) {
  307. N_VDestroyVectorArray(yS, NS); /* Free yS vector */
  308. }
  309. free(data); /* Free user data */
  310. CVodeFree(&cvode_mem); /* Free CVODES memory */
  311. SUNLinSolFree(LS); /* Free the linear solver memory */
  312. SUNMatDestroy(A); /* Free the matrix memory */
  313.  
  314. return(0);
  315. }
  316.  
  317. /*
  318. *--------------------------------------------------------------------
  319. * FUNCTIONS CALLED BY CVODES
  320. *--------------------------------------------------------------------
  321. */
  322.  
  323. /*
  324. * f routine. Compute f(t,y).
  325. */
  326.  
  327. static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
  328. {
  329. realtype y1, y2, y3, yd1, yd3;
  330. UserData data;
  331. realtype p1, p2, p3;
  332.  
  333. y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
  334. data = (UserData) user_data;
  335. p1 = data->p[0]; p2 = data->p[1]; p3 = data->p[2];
  336.  
  337. yd1 = Ith(ydot,1) = -p1*y1 + p2*y2*y3;
  338. yd3 = Ith(ydot,3) = p3*y2*y2;
  339. Ith(ydot,2) = -yd1 - yd3;
  340.  
  341. return(0);
  342. }
  343.  
  344.  
  345. /*
  346. * Jacobian routine. Compute J(t,y).
  347. */
  348.  
  349. static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
  350. void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
  351. {
  352. realtype y2, y3;
  353. UserData data;
  354. realtype p1, p2, p3;
  355.  
  356. y2 = Ith(y,2); y3 = Ith(y,3);
  357. data = (UserData) user_data;
  358. p1 = data->p[0]; p2 = data->p[1]; p3 = data->p[2];
  359.  
  360. IJth(J,1,1) = -p1; IJth(J,1,2) = p2*y3; IJth(J,1,3) = p2*y2;
  361. IJth(J,2,1) = p1; IJth(J,2,2) = -p2*y3-2*p3*y2; IJth(J,2,3) = -p2*y2;
  362. IJth(J,3,2) = 2*p3*y2;
  363.  
  364. return(0);
  365. }
  366.  
  367. /*
  368. * fS routine. Compute sensitivity r.h.s.
  369. */
  370.  
  371. static int fS(int Ns, realtype t, N_Vector y, N_Vector ydot,
  372. int iS, N_Vector yS, N_Vector ySdot,
  373. void *user_data, N_Vector tmp1, N_Vector tmp2)
  374. {
  375. UserData data;
  376. realtype p1, p2, p3;
  377. realtype y1, y2, y3;
  378. realtype s1, s2, s3;
  379. realtype sd1, sd2, sd3;
  380.  
  381. data = (UserData) user_data;
  382. p1 = data->p[0]; p2 = data->p[1]; p3 = data->p[2];
  383.  
  384. y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
  385. s1 = Ith(yS,1); s2 = Ith(yS,2); s3 = Ith(yS,3);
  386.  
  387. sd1 = -p1*s1 + p2*y3*s2 + p2*y2*s3;
  388. sd3 = 2*p3*y2*s2;
  389. sd2 = -sd1-sd3;
  390.  
  391. switch (iS) {
  392. case 0:
  393. sd1 += -y1;
  394. sd2 += y1;
  395. break;
  396. case 1:
  397. sd1 += y2*y3;
  398. sd2 += -y2*y3;
  399. break;
  400. case 2:
  401. sd2 += -y2*y2;
  402. sd3 += y2*y2;
  403. break;
  404. }
  405.  
  406. Ith(ySdot,1) = sd1;
  407. Ith(ySdot,2) = sd2;
  408. Ith(ySdot,3) = sd3;
  409.  
  410. return(0);
  411. }
  412.  
  413. /*
  414. * EwtSet function. Computes the error weights at the current solution.
  415. */
  416.  
  417. static int ewt(N_Vector y, N_Vector w, void *user_data)
  418. {
  419. int i;
  420. realtype yy, ww, rtol, atol[3];
  421.  
  422. rtol = RTOL;
  423. atol[0] = ATOL1;
  424. atol[1] = ATOL2;
  425. atol[2] = ATOL3;
  426.  
  427. for (i=1; i<=3; i++) {
  428. yy = Ith(y,i);
  429. ww = rtol * SUNRabs(yy) + atol[i-1];
  430. if (ww <= 0.0) return (-1);
  431. Ith(w,i) = 1.0/ww;
  432. }
  433.  
  434. return(0);
  435. }
  436.  
  437. /*
  438. *--------------------------------------------------------------------
  439. * PRIVATE FUNCTIONS
  440. *--------------------------------------------------------------------
  441. */
  442.  
  443. /*
  444. * Process and verify arguments to cvsfwddenx.
  445. */
  446.  
  447. static void ProcessArgs(int argc, char *argv[],
  448. booleantype *sensi, int *sensi_meth, booleantype *err_con)
  449. {
  450. *sensi = SUNFALSE;
  451. *sensi_meth = -1;
  452. *err_con = SUNFALSE;
  453.  
  454. if (argc < 2) WrongArgs(argv[0]);
  455.  
  456. if (strcmp(argv[1],"-nosensi") == 0)
  457. *sensi = SUNFALSE;
  458. else if (strcmp(argv[1],"-sensi") == 0)
  459. *sensi = SUNTRUE;
  460. else
  461. WrongArgs(argv[0]);
  462.  
  463. if (*sensi) {
  464.  
  465. if (argc != 4)
  466. WrongArgs(argv[0]);
  467.  
  468. if (strcmp(argv[2],"sim") == 0)
  469. *sensi_meth = CV_SIMULTANEOUS;
  470. else if (strcmp(argv[2],"stg") == 0)
  471. *sensi_meth = CV_STAGGERED;
  472. else if (strcmp(argv[2],"stg1") == 0)
  473. *sensi_meth = CV_STAGGERED1;
  474. else
  475. WrongArgs(argv[0]);
  476.  
  477. if (strcmp(argv[3],"t") == 0)
  478. *err_con = SUNTRUE;
  479. else if (strcmp(argv[3],"f") == 0)
  480. *err_con = SUNFALSE;
  481. else
  482. WrongArgs(argv[0]);
  483. }
  484.  
  485. }
  486.  
  487. static void WrongArgs(char *name)
  488. {
  489. printf("\nUsage: %s [-nosensi] [-sensi sensi_meth err_con]\n",name);
  490. printf(" sensi_meth = sim, stg, or stg1\n");
  491. printf(" err_con = t or f\n");
  492.  
  493. exit(0);
  494. }
  495.  
  496. /*
  497. * Print current t, step count, order, stepsize, and solution.
  498. */
  499.  
  500. static void PrintOutput(void *cvode_mem, realtype t, N_Vector u)
  501. {
  502. long int nst;
  503. int qu, retval;
  504. realtype hu, *udata;
  505.  
  506. udata = N_VGetArrayPointer(u);
  507.  
  508. retval = CVodeGetNumSteps(cvode_mem, &nst);
  509. check_retval(&retval, "CVodeGetNumSteps", 1);
  510. retval = CVodeGetLastOrder(cvode_mem, &qu);
  511. check_retval(&retval, "CVodeGetLastOrder", 1);
  512. retval = CVodeGetLastStep(cvode_mem, &hu);
  513. check_retval(&retval, "CVodeGetLastStep", 1);
  514.  
  515. #if defined(SUNDIALS_EXTENDED_PRECISION)
  516. printf("%8.3Le %2d %8.3Le %5ld\n", t, qu, hu, nst);
  517. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  518. printf("%8.3e %2d %8.3e %5ld\n", t, qu, hu, nst);
  519. #else
  520. printf("%8.3e %2d %8.3e %5ld\n", t, qu, hu, nst);
  521. #endif
  522.  
  523. printf(" Solution ");
  524.  
  525. #if defined(SUNDIALS_EXTENDED_PRECISION)
  526. printf("%12.4Le %12.4Le %12.4Le \n", udata[0], udata[1], udata[2]);
  527. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  528. printf("%12.4e %12.4e %12.4e \n", udata[0], udata[1], udata[2]);
  529. #else
  530. printf("%12.4e %12.4e %12.4e \n", udata[0], udata[1], udata[2]);
  531. #endif
  532.  
  533. }
  534.  
  535. /*
  536. * Print sensitivities.
  537. */
  538.  
  539. static void PrintOutputS(N_Vector *uS)
  540. {
  541. realtype *sdata;
  542.  
  543. sdata = N_VGetArrayPointer(uS[0]);
  544. printf(" Sensitivity 1 ");
  545.  
  546. #if defined(SUNDIALS_EXTENDED_PRECISION)
  547. printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
  548. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  549. printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  550. #else
  551. printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  552. #endif
  553.  
  554. sdata = N_VGetArrayPointer(uS[1]);
  555. printf(" Sensitivity 2 ");
  556.  
  557. #if defined(SUNDIALS_EXTENDED_PRECISION)
  558. printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
  559. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  560. printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  561. #else
  562. printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  563. #endif
  564.  
  565. sdata = N_VGetArrayPointer(uS[2]);
  566. printf(" Sensitivity 3 ");
  567.  
  568. #if defined(SUNDIALS_EXTENDED_PRECISION)
  569. printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
  570. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  571. printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  572. #else
  573. printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  574. #endif
  575. }
  576.  
  577. /*
  578. * Print some final statistics from the CVODES memory.
  579. */
  580.  
  581. static void PrintFinalStats(void *cvode_mem, booleantype sensi)
  582. {
  583. long int nst;
  584. long int nfe, nsetups, nni, ncfn, netf;
  585. long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;
  586. long int nje, nfeLS;
  587. int retval;
  588.  
  589. retval = CVodeGetNumSteps(cvode_mem, &nst);
  590. check_retval(&retval, "CVodeGetNumSteps", 1);
  591. retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);
  592. check_retval(&retval, "CVodeGetNumRhsEvals", 1);
  593. retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
  594. check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);
  595. retval = CVodeGetNumErrTestFails(cvode_mem, &netf);
  596. check_retval(&retval, "CVodeGetNumErrTestFails", 1);
  597. retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
  598. check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);
  599. retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
  600. check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);
  601.  
  602. if (sensi) {
  603. retval = CVodeGetSensNumRhsEvals(cvode_mem, &nfSe);
  604. check_retval(&retval, "CVodeGetSensNumRhsEvals", 1);
  605. retval = CVodeGetNumRhsEvalsSens(cvode_mem, &nfeS);
  606. check_retval(&retval, "CVodeGetNumRhsEvalsSens", 1);
  607. retval = CVodeGetSensNumLinSolvSetups(cvode_mem, &nsetupsS);
  608. check_retval(&retval, "CVodeGetSensNumLinSolvSetups", 1);
  609. retval = CVodeGetSensNumErrTestFails(cvode_mem, &netfS);
  610. check_retval(&retval, "CVodeGetSensNumErrTestFails", 1);
  611. retval = CVodeGetSensNumNonlinSolvIters(cvode_mem, &nniS);
  612. check_retval(&retval, "CVodeGetSensNumNonlinSolvIters", 1);
  613. retval = CVodeGetSensNumNonlinSolvConvFails(cvode_mem, &ncfnS);
  614. check_retval(&retval, "CVodeGetSensNumNonlinSolvConvFails", 1);
  615. }
  616.  
  617. retval = CVodeGetNumJacEvals(cvode_mem, &nje);
  618. check_retval(&retval, "CVodeGetNumJacEvals", 1);
  619. retval = CVodeGetNumLinRhsEvals(cvode_mem, &nfeLS);
  620. check_retval(&retval, "CVodeGetNumLinRhsEvals", 1);
  621.  
  622. printf("\nFinal Statistics\n\n");
  623. printf("nst = %5ld\n\n", nst);
  624. printf("nfe = %5ld\n", nfe);
  625. printf("netf = %5ld nsetups = %5ld\n", netf, nsetups);
  626. printf("nni = %5ld ncfn = %5ld\n", nni, ncfn);
  627.  
  628. if(sensi) {
  629. printf("\n");
  630. printf("nfSe = %5ld nfeS = %5ld\n", nfSe, nfeS);
  631. printf("netfs = %5ld nsetupsS = %5ld\n", netfS, nsetupsS);
  632. printf("nniS = %5ld ncfnS = %5ld\n", nniS, ncfnS);
  633. }
  634.  
  635. printf("\n");
  636. printf("nje = %5ld nfeLS = %5ld\n", nje, nfeLS);
  637.  
  638. }
  639.  
  640. /*
  641. * Check function return value.
  642. * opt == 0 means SUNDIALS function allocates memory so check if
  643. * returned NULL pointer
  644. * opt == 1 means SUNDIALS function returns an integer value so check if
  645. * retval < 0
  646. * opt == 2 means function allocates memory so check if returned
  647. * NULL pointer
  648. */
  649.  
  650. static int check_retval(void *returnvalue, const char *funcname, int opt)
  651. {
  652. int *retval;
  653.  
  654. /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
  655. if (opt == 0 && returnvalue == NULL) {
  656. fprintf(stderr,
  657. "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
  658. funcname);
  659. return(1); }
  660.  
  661. /* Check if retval < 0 */
  662. else if (opt == 1) {
  663. retval = (int *) returnvalue;
  664. if (*retval < 0) {
  665. fprintf(stderr,
  666. "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
  667. funcname, *retval);
  668. return(1); }}
  669.  
  670. /* Check if function returned NULL pointer - no memory allocated */
  671. else if (opt == 2 && returnvalue == NULL) {
  672. fprintf(stderr,
  673. "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
  674. funcname);
  675. return(1); }
  676.  
  677. return(0);
  678. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!