• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest May 25th, 2019 69 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. /* -----------------------------------------------------------------
2.  * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, and
3.  *                Radu Serban @ LLNL
4.  * -----------------------------------------------------------------
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.
12.  * For details, see the LICENSE file.
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.  */
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. }
RAW Paste Data
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.

Top