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