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