1 /*
2 * -----------------------------------------------------------------
3 * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, George D. Byrne,
4 * and 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 program for
19 * its solution by CVODES. The problem is the semi-discrete form of
20 * the advection-diffusion equation in 1-D:
21 * du/dt = q1 * d^2 u / dx^2 + q2 * du/dx
22 * on the interval 0 <= x <= 2, and the time interval 0 <= t <= 5.
23 * Homogeneous Dirichlet boundary conditions are posed, and the
24 * initial condition is:
25 * u(x,y,t=0) = x(2-x)exp(2x).
26 * The PDE is discretized on a uniform grid of size MX+2 with
27 * central differencing, and with boundary values eliminated,
28 * leaving an ODE system of size NEQ = MX.
29 * This program solves the problem with the option for nonstiff
30 * systems: ADAMS method and functional iteration.
31 * It uses scalar relative and absolute tolerances.
32 * Output is printed at t = .5, 1.0, ..., 5.
33 * Run statistics (optional outputs) are printed at the end.
34 *
35 * Optionally, CVODES can compute sensitivities with respect to the
36 * problem parameters q1 and q2.
37 * Any of three sensitivity methods (SIMULTANEOUS, STAGGERED, and
38 * STAGGERED1) can be used and sensitivities may be included in the
39 * error test or not (error control set on FULL or PARTIAL,
40 * respectively).
41 *
42 * Execution:
43 *
44 * Note: This version uses MPI for user routines, and the CVODES
45 * solver. In what follows, N is the number of processors,
46 * N = NPEX*NPEY (see constants below) and it is assumed that
47 * the MPI script mpirun is used to run a parallel
48 * application.
49 * If no sensitivities are desired:
50 * % mpirun -np N cvsAdvDiff_FSA_non_p -nosensi
51 * If sensitivities are to be computed:
52 * % mpirun -np N cvsAdvDiff_FSA_non_p -sensi sensi_meth err_con
53 * where sensi_meth is one of {sim, stg, stg1} and err_con is one of
54 * {t, f}.
55 * -----------------------------------------------------------------
56 */
57
58 #include <stdio.h>
59 #include <stdlib.h>
60 #include <math.h>
61 #include <string.h>
62
63 #include <cvodes/cvodes.h>
64 #include <nvector/nvector_parallel.h>
65 #include <sundials/sundials_types.h>
66 #include "sunnonlinsol/sunnonlinsol_fixedpoint.h" /* access to the fixed point SUNNonlinearSolver */
67
68 #include <mpi.h>
69
70 /* Problem Constants */
71 #define XMAX RCONST(2.0) /* domain boundary */
72 #define MX 10 /* mesh dimension */
73 #define NEQ MX /* number of equations */
74 #define ATOL RCONST(1.e-5) /* scalar absolute tolerance */
75 #define T0 RCONST(0.0) /* initial time */
76 #define T1 RCONST(0.5) /* first output time */
77 #define DTOUT RCONST(0.5) /* output time increment */
78 #define NOUT 10 /* number of output times */
79
80 #define NP 2
81 #define NS 2
82
83 #define ZERO RCONST(0.0)
84
85 /* Type : UserData
86 contains problem parameters, grid constants, work array. */
87
88 typedef struct {
89 realtype *p;
90 realtype dx;
91 int npes, my_pe;
92 MPI_Comm comm;
93 realtype z[100];
94 } *UserData;
95
96
97 /* Prototypes of user-supplied functins */
98
99 static int f(realtype t, N_Vector u, N_Vector udot, void *user_data);
100
101 /* Prototypes of private functions */
102
103 static void ProcessArgs(int argc, char *argv[], int my_pe,
104 booleantype *sensi, int *sensi_meth, booleantype *err_con);
105 static void WrongArgs(int my_pe, char *name);
106 static void SetIC(N_Vector u, realtype dx, sunindextype my_length, sunindextype my_base);
107 static void PrintOutput(void *cvode_mem, int my_pe, realtype t, N_Vector u);
108 static void PrintOutputS(int my_pe, N_Vector *uS);
109 static void PrintFinalStats(void *cvode_mem, booleantype sensi,
110 booleantype err_con, int sensi_meth);
111 static int check_retval(void *returnvalue, const char *funcname, int opt, int id);
112
113 /*
114 *--------------------------------------------------------------------
115 * MAIN PROGRAM
116 *--------------------------------------------------------------------
117 */
118
main(int argc,char * argv[])119 int main(int argc, char *argv[])
120 {
121 realtype dx, reltol, abstol, t, tout;
122 N_Vector u;
123 UserData data;
124 void *cvode_mem;
125 int iout, retval, my_pe, npes;
126 sunindextype local_N, nperpe, nrem, my_base;
127
128 realtype *pbar;
129 int is, *plist;
130 N_Vector *uS;
131 booleantype sensi, err_con;
132 int sensi_meth;
133
134 SUNNonlinearSolver NLS, NLSsens;
135
136 MPI_Comm comm;
137
138 u = NULL;
139 data = NULL;
140 cvode_mem = NULL;
141 pbar = NULL;
142 plist = NULL;
143 uS = NULL;
144 NLS = NULL;
145 NLSsens = NULL;
146
147 /* Get processor number, total number of pe's, and my_pe. */
148 MPI_Init(&argc, &argv);
149 comm = MPI_COMM_WORLD;
150 MPI_Comm_size(comm, &npes);
151 MPI_Comm_rank(comm, &my_pe);
152
153 /* Process arguments */
154 ProcessArgs(argc, argv, my_pe, &sensi, &sensi_meth, &err_con);
155
156 /* Set local vector length. */
157 nperpe = NEQ/npes;
158 nrem = NEQ - npes*nperpe;
159 local_N = (my_pe < nrem) ? nperpe+1 : nperpe;
160 my_base = (my_pe < nrem) ? my_pe*local_N : my_pe*nperpe + nrem;
161
162 /* USER DATA STRUCTURE */
163 data = (UserData) malloc(sizeof *data); /* Allocate data memory */
164 data->p = NULL;
165 if(check_retval((void *)data, "malloc", 2, my_pe)) MPI_Abort(comm, 1);
166 data->comm = comm;
167 data->npes = npes;
168 data->my_pe = my_pe;
169 data->p = (realtype *) malloc(NP * sizeof(realtype));
170 if(check_retval((void *)data->p, "malloc", 2, my_pe)) MPI_Abort(comm, 1);
171 dx = data->dx = XMAX/((realtype)(MX+1));
172 data->p[0] = RCONST(1.0);
173 data->p[1] = RCONST(0.5);
174
175 /* INITIAL STATES */
176 u = N_VNew_Parallel(comm, local_N, NEQ); /* Allocate u vector */
177 if(check_retval((void *)u, "N_VNew_Parallel", 0, my_pe)) MPI_Abort(comm, 1);
178 SetIC(u, dx, local_N, my_base); /* Initialize u vector */
179
180 /* TOLERANCES */
181 reltol = ZERO; /* Set the tolerances */
182 abstol = ATOL;
183
184 /* CVODE_CREATE & CVODE_MALLOC */
185 cvode_mem = CVodeCreate(CV_ADAMS);
186 if(check_retval((void *)cvode_mem, "CVodeCreate", 0, my_pe)) MPI_Abort(comm, 1);
187
188 retval = CVodeSetUserData(cvode_mem, data);
189 if(check_retval(&retval, "CVodeSetUserData", 1, my_pe)) MPI_Abort(comm, 1);
190
191 retval = CVodeInit(cvode_mem, f, T0, u);
192 if(check_retval(&retval, "CVodeInit", 1, my_pe)) MPI_Abort(comm, 1);
193
194 retval = CVodeSStolerances(cvode_mem, reltol, abstol);
195 if(check_retval(&retval, "CVodeSStolerances", 1, my_pe)) MPI_Abort(comm, 1);
196
197 /* create fixed point nonlinear solver object */
198 NLS = SUNNonlinSol_FixedPoint(u, 0);
199 if(check_retval((void *)NLS, "SUNNonlinSol_FixedPoint", 0, my_pe)) MPI_Abort(comm, 1);
200
201 /* attach nonlinear solver object to CVode */
202 retval = CVodeSetNonlinearSolver(cvode_mem, NLS);
203 if(check_retval(&retval, "CVodeSetNonlinearSolver", 1, my_pe)) MPI_Abort(comm, 1);
204
205 if (my_pe == 0) {
206 printf("\n1-D advection-diffusion equation, mesh size =%3d \n", MX);
207 printf("\nNumber of PEs = %3d \n",npes);
208 }
209
210 if(sensi) {
211
212 plist = (int *) malloc(NS * sizeof(int));
213 if(check_retval((void *)plist, "malloc", 2, my_pe)) MPI_Abort(comm, 1);
214 for(is=0; is<NS; is++)
215 plist[is] = is; /* sensitivity w.r.t. i-th parameter */
216
217 pbar = (realtype *) malloc(NS * sizeof(realtype));
218 if(check_retval((void *)pbar, "malloc", 2, my_pe)) MPI_Abort(comm, 1);
219 for(is=0; is<NS; is++) pbar[is] = data->p[plist[is]];
220
221 uS = N_VCloneVectorArray_Parallel(NS, u);
222 if(check_retval((void *)uS, "N_VCloneVectorArray_Parallel", 0, my_pe))
223 MPI_Abort(comm, 1);
224 for(is=0;is<NS;is++)
225 N_VConst(ZERO,uS[is]);
226
227 retval = CVodeSensInit1(cvode_mem, NS, sensi_meth, NULL, uS);
228 if(check_retval(&retval, "CVodeSensInit1", 1, my_pe)) MPI_Abort(comm, 1);
229
230 retval = CVodeSensEEtolerances(cvode_mem);
231 if(check_retval(&retval, "CVodeSensEEtolerances", 1, my_pe)) MPI_Abort(comm, 1);
232
233 retval = CVodeSetSensErrCon(cvode_mem, err_con);
234 if(check_retval(&retval, "CVodeSetSensErrCon", 1, my_pe)) MPI_Abort(comm, 1);
235
236 retval = CVodeSetSensDQMethod(cvode_mem, CV_CENTERED, ZERO);
237 if(check_retval(&retval, "CVodeSetSensDQMethod", 1, my_pe)) MPI_Abort(comm, 1);
238
239 retval = CVodeSetSensParams(cvode_mem, data->p, pbar, plist);
240 if(check_retval(&retval, "CVodeSetSensParams", 1, my_pe)) MPI_Abort(comm, 1);
241
242 /* create sensitivity fixed point nonlinear solver object */
243 if (sensi_meth == CV_SIMULTANEOUS)
244 NLSsens = SUNNonlinSol_FixedPointSens(NS+1, u, 0);
245 else if(sensi_meth == CV_STAGGERED)
246 NLSsens = SUNNonlinSol_FixedPointSens(NS, u, 0);
247 else
248 NLSsens = SUNNonlinSol_FixedPoint(u, 0);
249 if(check_retval((void *)NLS, "SUNNonlinSol_FixedPoint", 0, my_pe)) MPI_Abort(comm, 1);
250
251 /* attach nonlinear solver object to CVode */
252 if (sensi_meth == CV_SIMULTANEOUS)
253 retval = CVodeSetNonlinearSolverSensSim(cvode_mem, NLSsens);
254 else if(sensi_meth == CV_STAGGERED)
255 retval = CVodeSetNonlinearSolverSensStg(cvode_mem, NLSsens);
256 else
257 retval = CVodeSetNonlinearSolverSensStg1(cvode_mem, NLSsens);
258 if(check_retval(&retval, "CVodeSetNonlinearSolver", 1, my_pe)) MPI_Abort(comm, 1);
259
260 if(my_pe == 0) {
261 printf("Sensitivity: YES ");
262 if(sensi_meth == CV_SIMULTANEOUS)
263 printf("( SIMULTANEOUS +");
264 else
265 if(sensi_meth == CV_STAGGERED) printf("( STAGGERED +");
266 else printf("( STAGGERED1 +");
267 if(err_con) printf(" FULL ERROR CONTROL )");
268 else printf(" PARTIAL ERROR CONTROL )");
269 }
270
271 } else {
272
273 if(my_pe == 0) printf("Sensitivity: NO ");
274
275 }
276
277 /* In loop over output points, call CVode, print results, test for error */
278
279 if(my_pe == 0) {
280 printf("\n\n");
281 printf("============================================================\n");
282 printf(" T Q H NST Max norm \n");
283 printf("============================================================\n");
284 }
285
286 for (iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {
287
288 retval = CVode(cvode_mem, tout, u, &t, CV_NORMAL);
289 if(check_retval(&retval, "CVode", 1, my_pe)) break;
290 PrintOutput(cvode_mem, my_pe, t, u);
291 if (sensi) {
292 retval = CVodeGetSens(cvode_mem, &t, uS);
293 if(check_retval(&retval, "CVodeGetSens", 1, my_pe)) break;
294 PrintOutputS(my_pe, uS);
295 }
296 if (my_pe == 0)
297 printf("------------------------------------------------------------\n");
298
299 }
300
301 /* Print final statistics */
302 if (my_pe == 0)
303 PrintFinalStats(cvode_mem, sensi, err_con, sensi_meth);
304
305 /* Free memory */
306 N_VDestroy(u); /* Free the u vector */
307 if (sensi)
308 N_VDestroyVectorArray(uS, NS); /* Free the uS vectors */
309 free(data->p); /* Free the p vector */
310 free(data); /* Free block of UserData */
311 CVodeFree(&cvode_mem); /* Free the CVODES problem memory */
312 SUNNonlinSolFree(NLS);
313 free(pbar);
314 if(sensi) free(plist);
315 if(sensi) SUNNonlinSolFree(NLSsens);
316
317 MPI_Finalize();
318
319 return(0);
320 }
321
322 /*
323 *--------------------------------------------------------------------
324 * FUNCTIONS CALLED BY CVODES
325 *--------------------------------------------------------------------
326 */
327
328 /*
329 * f routine. Compute f(t,u).
330 */
331
f(realtype t,N_Vector u,N_Vector udot,void * user_data)332 static int f(realtype t, N_Vector u, N_Vector udot, void *user_data)
333 {
334 realtype ui, ult, urt, hordc, horac, hdiff, hadv;
335 realtype *udata, *dudata, *z;
336 realtype dx;
337 int npes, my_pe, my_pe_m1, my_pe_p1, last_pe;
338 sunindextype i, my_length;
339 UserData data;
340 MPI_Status status;
341 MPI_Comm comm;
342
343 udata = N_VGetArrayPointer_Parallel(u);
344 dudata = N_VGetArrayPointer_Parallel(udot);
345
346 /* Extract needed problem constants from data */
347 data = (UserData) user_data;
348 dx = data->dx;
349 hordc = data->p[0]/(dx*dx);
350 horac = data->p[1]/(RCONST(2.0)*dx);
351
352 /* Extract parameters for parallel computation. */
353 comm = data->comm;
354 npes = data->npes; /* Number of processes. */
355 my_pe = data->my_pe; /* Current process number. */
356 my_length = N_VGetLocalLength_Parallel(u); /* Number of local elements of u. */
357 z = data->z;
358
359 /* Compute related parameters. */
360 my_pe_m1 = my_pe - 1;
361 my_pe_p1 = my_pe + 1;
362 last_pe = npes - 1;
363
364 /* Store local segment of u in the working array z. */
365 for (i = 1; i <= my_length; i++)
366 z[i] = udata[i - 1];
367
368 /* Pass needed data to processes before and after current process. */
369 if (my_pe != 0)
370 MPI_Send(&z[1], 1, MPI_SUNREALTYPE, my_pe_m1, 0, comm);
371 if (my_pe != last_pe)
372 MPI_Send(&z[my_length], 1, MPI_SUNREALTYPE, my_pe_p1, 0, comm);
373
374 /* Receive needed data from processes before and after current process. */
375 if (my_pe != 0)
376 MPI_Recv(&z[0], 1, MPI_SUNREALTYPE, my_pe_m1, 0, comm, &status);
377 else z[0] = ZERO;
378 if (my_pe != last_pe)
379 MPI_Recv(&z[my_length+1], 1, MPI_SUNREALTYPE, my_pe_p1, 0, comm,
380 &status);
381 else z[my_length + 1] = ZERO;
382
383 /* Loop over all grid points in current process. */
384 for (i=1; i<=my_length; i++) {
385
386 /* Extract u at x_i and two neighboring points */
387 ui = z[i];
388 ult = z[i-1];
389 urt = z[i+1];
390
391 /* Set diffusion and advection terms and load into udot */
392 hdiff = hordc*(ult - RCONST(2.0)*ui + urt);
393 hadv = horac*(urt - ult);
394 dudata[i-1] = hdiff + hadv;
395 }
396
397 return(0);
398 }
399
400 /*
401 *--------------------------------------------------------------------
402 * PRIVATE FUNCTIONS
403 *--------------------------------------------------------------------
404 */
405
406 /*
407 * Process and verify arguments to cvsfwdnonx_p.
408 */
409
ProcessArgs(int argc,char * argv[],int my_pe,booleantype * sensi,int * sensi_meth,booleantype * err_con)410 static void ProcessArgs(int argc, char *argv[], int my_pe,
411 booleantype *sensi, int *sensi_meth, booleantype *err_con)
412 {
413 *sensi = SUNFALSE;
414 *sensi_meth = -1;
415 *err_con = SUNFALSE;
416
417 if (argc < 2) WrongArgs(my_pe, argv[0]);
418
419 if (strcmp(argv[1],"-nosensi") == 0)
420 *sensi = SUNFALSE;
421 else if (strcmp(argv[1],"-sensi") == 0)
422 *sensi = SUNTRUE;
423 else
424 WrongArgs(my_pe, argv[0]);
425
426 if (*sensi) {
427
428 if (argc != 4)
429 WrongArgs(my_pe, argv[0]);
430
431 if (strcmp(argv[2],"sim") == 0)
432 *sensi_meth = CV_SIMULTANEOUS;
433 else if (strcmp(argv[2],"stg") == 0)
434 *sensi_meth = CV_STAGGERED;
435 else if (strcmp(argv[2],"stg1") == 0)
436 *sensi_meth = CV_STAGGERED1;
437 else
438 WrongArgs(my_pe, argv[0]);
439
440 if (strcmp(argv[3],"t") == 0)
441 *err_con = SUNTRUE;
442 else if (strcmp(argv[3],"f") == 0)
443 *err_con = SUNFALSE;
444 else
445 WrongArgs(my_pe, argv[0]);
446 }
447
448 }
449
WrongArgs(int my_pe,char * name)450 static void WrongArgs(int my_pe, char *name)
451 {
452 if (my_pe == 0) {
453 printf("\nUsage: %s [-nosensi] [-sensi sensi_meth err_con]\n",name);
454 printf(" sensi_meth = sim, stg, or stg1\n");
455 printf(" err_con = t or f\n");
456 }
457 MPI_Finalize();
458 exit(0);
459 }
460
461 /*
462 * Set initial conditions in u vector
463 */
464
SetIC(N_Vector u,realtype dx,sunindextype my_length,sunindextype my_base)465 static void SetIC(N_Vector u, realtype dx, sunindextype my_length,
466 sunindextype my_base)
467 {
468 int i;
469 sunindextype iglobal;
470 realtype x;
471 realtype *udata;
472
473 /* Set pointer to data array and get local length of u. */
474 udata = N_VGetArrayPointer_Parallel(u);
475 my_length = N_VGetLocalLength_Parallel(u);
476
477 /* Load initial profile into u vector */
478 for (i=1; i<=my_length; i++) {
479 iglobal = my_base + i;
480 x = iglobal*dx;
481 udata[i-1] = x*(XMAX - x)*exp(2.0*x);
482 }
483 }
484
485 /*
486 * Print current t, step count, order, stepsize, and max norm of solution
487 */
488
PrintOutput(void * cvode_mem,int my_pe,realtype t,N_Vector u)489 static void PrintOutput(void *cvode_mem, int my_pe, realtype t, N_Vector u)
490 {
491 long int nst;
492 int qu, retval;
493 realtype hu, umax;
494
495 retval = CVodeGetNumSteps(cvode_mem, &nst);
496 check_retval(&retval, "CVodeGetNumSteps", 1, my_pe);
497 retval = CVodeGetLastOrder(cvode_mem, &qu);
498 check_retval(&retval, "CVodeGetLastOrder", 1, my_pe);
499 retval = CVodeGetLastStep(cvode_mem, &hu);
500 check_retval(&retval, "CVodeGetLastStep", 1, my_pe);
501
502 umax = N_VMaxNorm(u);
503
504 if (my_pe == 0) {
505
506 #if defined(SUNDIALS_EXTENDED_PRECISION)
507 printf("%8.3Le %2d %8.3Le %5ld\n", t,qu,hu,nst);
508 #elif defined(SUNDIALS_DOUBLE_PRECISION)
509 printf("%8.3e %2d %8.3e %5ld\n", t,qu,hu,nst);
510 #else
511 printf("%8.3e %2d %8.3e %5ld\n", t,qu,hu,nst);
512 #endif
513
514 printf(" Solution ");
515
516 #if defined(SUNDIALS_EXTENDED_PRECISION)
517 printf("%12.4Le \n", umax);
518 #elif defined(SUNDIALS_DOUBLE_PRECISION)
519 printf("%12.4e \n", umax);
520 #else
521 printf("%12.4e \n", umax);
522 #endif
523
524 }
525
526 }
527
528 /*
529 * Print max norm of sensitivities
530 */
531
PrintOutputS(int my_pe,N_Vector * uS)532 static void PrintOutputS(int my_pe, N_Vector *uS)
533 {
534 realtype smax;
535
536 smax = N_VMaxNorm(uS[0]);
537 if (my_pe == 0) {
538 printf(" Sensitivity 1 ");
539 #if defined(SUNDIALS_EXTENDED_PRECISION)
540 printf("%12.4Le \n", smax);
541 #elif defined(SUNDIALS_DOUBLE_PRECISION)
542 printf("%12.4e \n", smax);
543 #else
544 printf("%12.4e \n", smax);
545 #endif
546 }
547
548 smax = N_VMaxNorm(uS[1]);
549 if (my_pe == 0) {
550 printf(" Sensitivity 2 ");
551 #if defined(SUNDIALS_EXTENDED_PRECISION)
552 printf("%12.4Le \n", smax);
553 #elif defined(SUNDIALS_DOUBLE_PRECISION)
554 printf("%12.4e \n", smax);
555 #else
556 printf("%12.4e \n", smax);
557 #endif
558 }
559
560 }
561
562 /*
563 * Print some final statistics located in the iopt array
564 */
565
PrintFinalStats(void * cvode_mem,booleantype sensi,booleantype err_con,int sensi_meth)566 static void PrintFinalStats(void *cvode_mem, booleantype sensi,
567 booleantype err_con, int sensi_meth)
568 {
569 long int nst;
570 long int nfe, nsetups, nni, ncfn, netf;
571 long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;
572 int retval;
573
574 retval = CVodeGetNumSteps(cvode_mem, &nst);
575 check_retval(&retval, "CVodeGetNumSteps", 1, 0);
576 retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);
577 check_retval(&retval, "CVodeGetNumRhsEvals", 1, 0);
578 retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
579 check_retval(&retval, "CVodeGetNumLinSolvSetups", 1, 0);
580 retval = CVodeGetNumErrTestFails(cvode_mem, &netf);
581 check_retval(&retval, "CVodeGetNumErrTestFails", 1, 0);
582 retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
583 check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1, 0);
584 retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
585 check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1, 0);
586
587 if (sensi) {
588 retval = CVodeGetSensNumRhsEvals(cvode_mem, &nfSe);
589 check_retval(&retval, "CVodeGetSensNumRhsEvals", 1, 0);
590 retval = CVodeGetNumRhsEvalsSens(cvode_mem, &nfeS);
591 check_retval(&retval, "CVodeGetNumRhsEvalsSens", 1, 0);
592 retval = CVodeGetSensNumLinSolvSetups(cvode_mem, &nsetupsS);
593 check_retval(&retval, "CVodeGetSensNumLinSolvSetups", 1, 0);
594 retval = CVodeGetSensNumErrTestFails(cvode_mem, &netfS);
595 if (err_con) {
596 retval = CVodeGetSensNumErrTestFails(cvode_mem, &netfS);
597 check_retval(&retval, "CVodeGetSensNumErrTestFails", 1, 0);
598 } else {
599 netfS = 0;
600 }
601 if ((sensi_meth == CV_STAGGERED) || (sensi_meth == CV_STAGGERED1)) {
602 retval = CVodeGetSensNumNonlinSolvIters(cvode_mem, &nniS);
603 check_retval(&retval, "CVodeGetSensNumNonlinSolvIters", 1, 0);
604 retval = CVodeGetSensNumNonlinSolvConvFails(cvode_mem, &ncfnS);
605 check_retval(&retval, "CVodeGetSensNumNonlinSolvConvFails", 1, 0);
606 } else {
607 nniS = 0;
608 ncfnS = 0;
609 }
610 }
611
612 printf("\nFinal Statistics\n\n");
613 printf("nst = %5ld\n\n", nst);
614 printf("nfe = %5ld\n", nfe);
615 printf("netf = %5ld nsetups = %5ld\n", netf, nsetups);
616 printf("nni = %5ld ncfn = %5ld\n", nni, ncfn);
617
618 if(sensi) {
619 printf("\n");
620 printf("nfSe = %5ld nfeS = %5ld\n", nfSe, nfeS);
621 printf("netfs = %5ld nsetupsS = %5ld\n", netfS, nsetupsS);
622 printf("nniS = %5ld ncfnS = %5ld\n", nniS, ncfnS);
623 }
624
625 }
626
627 /*
628 * Check function return value...
629 * opt == 0 means SUNDIALS function allocates memory so check if
630 * returned NULL pointer
631 * opt == 1 means SUNDIALS function returns an integer value so check if
632 * retval < 0
633 * opt == 2 means function allocates memory so check if returned
634 * NULL pointer
635 */
636
check_retval(void * returnvalue,const char * funcname,int opt,int id)637 static int check_retval(void *returnvalue, const char *funcname, int opt, int id)
638 {
639 int *retval;
640
641 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
642 if (opt == 0 && returnvalue == NULL) {
643 fprintf(stderr, "\nSUNDIALS_ERROR(%d): %s() failed - returned NULL pointer\n\n",
644 id, funcname);
645 return(1); }
646
647 /* Check if retval < 0 */
648 else if (opt == 1) {
649 retval = (int *) returnvalue;
650 if (*retval < 0) {
651 fprintf(stderr, "\nSUNDIALS_ERROR(%d): %s() failed with retval = %d\n\n",
652 id, funcname, *retval);
653 return(1); }}
654
655 /* Check if function returned NULL pointer - no memory allocated */
656 else if (opt == 2 && returnvalue == NULL) {
657 fprintf(stderr, "\nMEMORY_ERROR(%d): %s() failed - returned NULL pointer\n\n",
658 id, funcname);
659 return(1); }
660
661 return(0);
662 }
663