1 /* -----------------------------------------------------------------
2 * Programmer(s): Cosmin Petra and Radu Serban @ LLNL
3 * -----------------------------------------------------------------
4 * SUNDIALS Copyright Start
5 * Copyright (c) 2002-2021, Lawrence Livermore National Security
6 * and Southern Methodist University.
7 * All rights reserved.
8 *
9 * See the top-level LICENSE and NOTICE files for details.
10 *
11 * SPDX-License-Identifier: BSD-3-Clause
12 * SUNDIALS Copyright End
13 * -----------------------------------------------------------------
14 * Example problem:
15 *
16 * This simple example problem for IDA, due to Robertson,
17 * is from chemical kinetics, and consists of the following three
18 * equations:
19 *
20 * dy1/dt = -p1*y1 + p2*y2*y3
21 * dy2/dt = p1*y1 - p2*y2*y3 - p3*y2**2
22 * 0 = y1 + y2 + y3 - 1
23 *
24 * on the interval from t = 0.0 to t = 4.e10, with initial
25 * conditions: y1 = 1, y2 = y3 = 0.The reaction rates are: p1=0.04,
26 * p2=1e4, and p3=3e7
27 *
28 * Optionally, IDAS can compute sensitivities with respect to the
29 * problem parameters p1, p2, and p3.
30 * The sensitivity right hand side is given analytically through the
31 * user routine fS (of type SensRhs1Fn).
32 * Any of two sensitivity methods (SIMULTANEOUS and STAGGERED can be
33 * used and sensitivities may be included in the error test or not
34 *(error control set on SUNTRUE or SUNFALSE, respectively).
35 *
36 * Execution:
37 *
38 * If no sensitivities are desired:
39 * % idasRoberts_FSA_dns -nosensi
40 * If sensitivities are to be computed:
41 * % idasRoberts_FSA_dns -sensi sensi_meth err_con
42 * where sensi_meth is one of {sim, stg} and err_con is one of
43 * {t, f}.
44 * -----------------------------------------------------------------*/
45
46 #include <stdio.h>
47 #include <stdlib.h>
48 #include <string.h>
49
50 #include <idas/idas.h> /* prototypes for IDA fcts., consts. */
51 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
52 #include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
53 #include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
54 #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
55 #include <sundials/sundials_math.h> /* defs. of SUNRabs, SUNRexp, etc. */
56
57 /* Accessor macros */
58
59 #define Ith(v,i) NV_Ith_S(v,i-1) /* i-th vector component i=1..NEQ */
60
61 /* Problem Constants */
62
63 #define NEQ 3 /* number of equations */
64 #define T0 RCONST(0.0) /* initial time */
65 #define T1 RCONST(0.4) /* first output time */
66 #define TMULT RCONST(10.0) /* output time factor */
67 #define NOUT 12 /* number of output times */
68
69 #define NP 3 /* number of problem parameters */
70 #define NS 3 /* number of sensitivities computed */
71
72 #define ZERO RCONST(0.0)
73 #define ONE RCONST(1.0)
74
75 /* Type : UserData */
76
77 typedef struct {
78 realtype p[3]; /* problem parameters */
79 realtype coef;
80 } *UserData;
81
82 /* Prototypes of functions by IDAS */
83
84 static int res(realtype t, N_Vector y, N_Vector yp, N_Vector resval, void *user_data);
85
86 static int resS(int Ns, realtype t,
87 N_Vector y, N_Vector yp, N_Vector resval,
88 N_Vector *yyS, N_Vector *ypS, N_Vector *resvalS,
89 void *user_data,
90 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
91
92 static int rhsQ(realtype tres, N_Vector yy, N_Vector yp,
93 N_Vector rrQ, void *user_data);
94
95 /* Prototypes of private functions */
96
97 static void ProcessArgs(int argc, char *argv[],
98 booleantype *sensi, int *sensi_meth,
99 booleantype *err_con);
100 static void WrongArgs(char *name);
101
102 static void PrintIC(N_Vector y, N_Vector yp);
103 static void PrintSensIC(N_Vector y, N_Vector yp, N_Vector* yS, N_Vector* ypS);
104
105 static void PrintOutput(void *ida_mem, realtype t, N_Vector u);
106 static void PrintSensOutput(N_Vector *uS);
107
108 static void PrintFinalStats(void *ida_mem, booleantype sensi);
109
110 static int check_retval(void *returnvalue, const char *funcname, int opt);
111 /*
112 *--------------------------------------------------------------------
113 * MAIN PROGRAM
114 *--------------------------------------------------------------------
115 */
116
main(int argc,char * argv[])117 int main(int argc, char *argv[])
118 {
119 void *ida_mem;
120 SUNMatrix A;
121 SUNLinearSolver LS;
122 UserData data;
123 realtype reltol, t, tout;
124 N_Vector y, yp, abstol, id;
125 int iout, retval;
126
127 realtype pbar[NS];
128 int is;
129 N_Vector *yS, *ypS;
130 booleantype sensi, err_con;
131 int sensi_meth;
132
133 N_Vector yQ, *yQS;
134
135 ida_mem = NULL;
136 A = NULL;
137 LS = NULL;
138 data = NULL;
139 y = NULL;
140 yS = NULL;
141 ypS = NULL;
142 yQS = NULL;
143
144 /* Process arguments */
145 ProcessArgs(argc, argv, &sensi, &sensi_meth, &err_con);
146
147 /* User data structure */
148 data = (UserData) malloc(sizeof *data);
149 if (check_retval((void *)data, "malloc", 2)) return(1);
150 data->p[0] = RCONST(0.040);
151 data->p[1] = RCONST(1.0e4);
152 data->p[2] = RCONST(3.0e7);
153 data->coef = 0.5;
154
155 /* Initial conditions */
156 y = N_VNew_Serial(NEQ);
157 if (check_retval((void *)y, "N_VNew_Serial", 0)) return(1);
158
159 Ith(y,1) = ONE;
160 Ith(y,2) = ZERO;
161 Ith(y,3) = ZERO;
162
163 yp = N_VNew_Serial(NEQ);
164 if(check_retval((void *)yp, "N_VNew_Serial", 0)) return(1);
165
166 /* These initial conditions are NOT consistent. See IDACalcIC below. */
167 Ith(yp,1) = RCONST(0.1);
168 Ith(yp,2) = ZERO;
169 Ith(yp,3) = ZERO;
170
171 /* Create IDAS object */
172 ida_mem = IDACreate();
173 if (check_retval((void *)ida_mem, "IDACreate", 0)) return(1);
174
175 /* Allocate space for IDAS */
176 retval = IDAInit(ida_mem, res, T0, y, yp);
177 if (check_retval(&retval, "IDAInit", 1)) return(1);
178
179 /* Specify scalar relative tol. and vector absolute tol. */
180 reltol = RCONST(1.0e-6);
181 abstol = N_VNew_Serial(NEQ);
182 Ith(abstol,1) = RCONST(1.0e-8);
183 Ith(abstol,2) = RCONST(1.0e-14);
184 Ith(abstol,3) = RCONST(1.0e-6);
185 retval = IDASVtolerances(ida_mem, reltol, abstol);
186 if (check_retval(&retval, "IDASVtolerances", 1)) return(1);
187
188 /* Set ID vector */
189 id = N_VNew_Serial(NEQ);
190 Ith(id,1) = 1.0;
191 Ith(id,2) = 1.0;
192 Ith(id,3) = 0.0;
193 retval = IDASetId(ida_mem, id);
194 if (check_retval(&retval, "IDASetId", 1)) return(1);
195
196 /* Attach user data */
197 retval = IDASetUserData(ida_mem, data);
198 if (check_retval(&retval, "IDASetUserData", 1)) return(1);
199
200 /* Create dense SUNMatrix for use in linear solves */
201 A = SUNDenseMatrix(NEQ, NEQ);
202 if(check_retval((void *)A, "SUNDenseMatrix", 0)) return(1);
203
204 /* Create dense SUNLinearSolver object */
205 LS = SUNLinSol_Dense(y, A);
206 if(check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1);
207
208 /* Attach the matrix and linear solver */
209 retval = IDASetLinearSolver(ida_mem, LS, A);
210 if(check_retval(&retval, "IDASetLinearSolver", 1)) return(1);
211
212 printf("\n3-species chemical kinetics problem\n");
213
214 /* Sensitivity-related settings */
215 if (sensi) {
216
217 pbar[0] = data->p[0];
218 pbar[1] = data->p[1];
219 pbar[2] = data->p[2];
220
221 yS = N_VCloneVectorArray(NS, y);
222 if (check_retval((void *)yS, "N_VCloneVectorArray", 0)) return(1);
223 for (is=0;is<NS;is++) N_VConst(ZERO, yS[is]);
224
225 ypS = N_VCloneVectorArray(NS, y);
226 if (check_retval((void *)ypS, "N_VCloneVectorArray", 0)) return(1);
227 for (is=0;is<NS;is++) N_VConst(ZERO, ypS[is]);
228
229 /*
230 * Only non-zero sensitivity I.C. are ypS[0]:
231 * - Ith(ypS[0],1) = -ONE;
232 * - Ith(ypS[0],2) = ONE;
233 *
234 * They are not set. IDACalcIC also computes consistent IC for sensitivities.
235 */
236
237 retval = IDASensInit(ida_mem, NS, sensi_meth, resS, yS, ypS);
238 if(check_retval(&retval, "IDASensInit", 1)) return(1);
239
240 retval = IDASensEEtolerances(ida_mem);
241 if(check_retval(&retval, "IDASensEEtolerances", 1)) return(1);
242
243 retval = IDASetSensErrCon(ida_mem, err_con);
244 if (check_retval(&retval, "IDASetSensErrCon", 1)) return(1);
245
246 retval = IDASetSensParams(ida_mem, data->p, pbar, NULL);
247 if (check_retval(&retval, "IDASetSensParams", 1)) return(1);
248
249 printf("Sensitivity: YES ");
250 if(sensi_meth == IDA_SIMULTANEOUS)
251 printf("( SIMULTANEOUS +");
252 else
253 printf("( STAGGERED +");
254 if(err_con) printf(" FULL ERROR CONTROL )");
255 else printf(" PARTIAL ERROR CONTROL )");
256
257 } else {
258
259 printf("Sensitivity: NO ");
260
261 }
262
263 /*----------------------------------------------------------
264 * Q U A D R A T U R E S
265 * ---------------------------------------------------------*/
266 yQ = N_VNew_Serial(2);
267
268 Ith(yQ,1) = 0;
269 Ith(yQ,2) = 0;
270
271 IDAQuadInit(ida_mem, rhsQ, yQ);
272
273 if (sensi) {
274 yQS = N_VCloneVectorArray(NS, yQ);
275 for (is=0;is<NS;is++) N_VConst(ZERO, yQS[is]);
276
277 IDAQuadSensInit(ida_mem, NULL, yQS);
278 }
279
280 /* Call IDACalcIC to compute consistent initial conditions. If sensitivity is
281 enabled, this function also try to find consistent IC for the sensitivities. */
282
283 retval = IDACalcIC(ida_mem, IDA_YA_YDP_INIT, T1);;
284 if (check_retval(&retval, "IDACalcIC", 1)) return(1);
285
286 retval = IDAGetConsistentIC(ida_mem, y, yp);
287 if (check_retval(&retval, "IDAGetConsistentIC", 1)) return(1);
288
289 PrintIC(y, yp);
290
291 if(sensi) {
292 IDAGetSensConsistentIC(ida_mem, yS, ypS);
293 PrintSensIC(y, yp, yS, ypS);
294 }
295
296 /* In loop over output points, call IDA, print results, test for error */
297
298 printf("\n\n");
299 printf("===========================================");
300 printf("============================\n");
301 printf(" T Q H NST y1");
302 printf(" y2 y3 \n");
303 printf("===========================================");
304 printf("============================\n");
305
306 for (iout=1, tout=T1; iout <= NOUT; iout++, tout *= TMULT) {
307
308 retval = IDASolve(ida_mem, tout, &t, y, yp, IDA_NORMAL);
309 if (check_retval(&retval, "IDASolve", 1)) break;
310
311 PrintOutput(ida_mem, t, y);
312
313 if (sensi) {
314 retval = IDAGetSens(ida_mem, &t, yS);
315 if (check_retval(&retval, "IDAGetSens", 1)) break;
316 PrintSensOutput(yS);
317 }
318 printf("-----------------------------------------");
319 printf("------------------------------\n");
320
321 }
322
323 printf("\nQuadrature:\n");
324 IDAGetQuad(ida_mem, &t, yQ);
325 #if defined(SUNDIALS_EXTENDED_PRECISION)
326 printf("G: %10.4Le\n", Ith(yQ,1));
327 #else
328 printf("G: %10.4e\n", Ith(yQ,1));
329 #endif
330
331 if(sensi) {
332 IDAGetQuadSens(ida_mem, &t, yQS);
333 #if defined(SUNDIALS_EXTENDED_PRECISION)
334 printf("\nSensitivities at t=%Lg:\n",t);
335 printf("dG/dp1: %11.4Le\n", Ith(yQS[0], 1));
336 printf("dG/dp1: %11.4Le\n", Ith(yQS[1], 1));
337 printf("dG/dp1: %11.4Le\n", Ith(yQS[2], 1));
338 #else
339 printf("\nSensitivities at t=%g:\n",t);
340 printf("dG/dp1: %11.4e\n", Ith(yQS[0], 1));
341 printf("dG/dp1: %11.4e\n", Ith(yQS[1], 1));
342 printf("dG/dp1: %11.4e\n", Ith(yQS[2], 1));
343 #endif
344 }
345
346 /* Print final statistics */
347 PrintFinalStats(ida_mem, sensi);
348
349 /* Free memory */
350 N_VDestroy(y);
351 N_VDestroy(yp);
352 N_VDestroy(abstol);
353 N_VDestroy(id);
354 N_VDestroy(yQ);
355 if (sensi) {
356 N_VDestroyVectorArray(yS, NS);
357 N_VDestroyVectorArray(ypS, NS);
358 N_VDestroyVectorArray(yQS, NS);
359 }
360 free(data);
361 IDAFree(&ida_mem);
362 SUNLinSolFree(LS);
363 SUNMatDestroy(A);
364
365 return(0);
366 }
367
368 /*
369 *--------------------------------------------------------------------
370 * FUNCTIONS CALLED BY IDAS
371 *--------------------------------------------------------------------
372 */
373
374 /*
375 * Residual routine. Compute F(t,y,y',p).
376 */
res(realtype t,N_Vector yy,N_Vector yp,N_Vector resval,void * user_data)377 static int res(realtype t, N_Vector yy, N_Vector yp, N_Vector resval, void *user_data)
378 {
379 UserData data;
380 realtype p1, p2, p3;
381 realtype y1, y2, y3;
382 realtype yp1, yp2;
383
384 data = (UserData) user_data;
385 p1 = data->p[0];
386 p2 = data->p[1];
387 p3 = data->p[2];
388
389 y1 = Ith(yy,1);
390 y2 = Ith(yy,2);
391 y3 = Ith(yy,3);
392
393 yp1 = Ith(yp,1);
394 yp2 = Ith(yp,2);
395
396 Ith(resval,1) = yp1 + p1*y1 - p2*y2*y3;
397 Ith(resval,2) = yp2 - p1*y1 + p2*y2*y3 + p3*y2*y2;
398 Ith(resval,3) = y1 + y2 + y3 - ONE;
399
400 return(0);
401 }
402
403
404 /*
405 * resS routine. Compute sensitivity r.h.s.
406 */
407
resS(int Ns,realtype t,N_Vector yy,N_Vector yp,N_Vector resval,N_Vector * yyS,N_Vector * ypS,N_Vector * resvalS,void * user_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)408 static int resS(int Ns, realtype t,
409 N_Vector yy, N_Vector yp, N_Vector resval,
410 N_Vector *yyS, N_Vector *ypS, N_Vector *resvalS,
411 void *user_data,
412 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
413 {
414 UserData data;
415 realtype p1, p2, p3;
416 realtype y1, y2, y3;
417 realtype s1, s2, s3;
418 realtype sd1, sd2;
419 realtype rs1, rs2, rs3;
420 int is;
421
422 data = (UserData) user_data;
423 p1 = data->p[0];
424 p2 = data->p[1];
425 p3 = data->p[2];
426
427 y1 = Ith(yy,1);
428 y2 = Ith(yy,2);
429 y3 = Ith(yy,3);
430
431 for (is=0; is<NS; is++) {
432
433 s1 = Ith(yyS[is],1);
434 s2 = Ith(yyS[is],2);
435 s3 = Ith(yyS[is],3);
436
437 sd1 = Ith(ypS[is],1);
438 sd2 = Ith(ypS[is],2);
439
440 rs1 = sd1 + p1*s1 - p2*y3*s2 - p2*y2*s3;
441 rs2 = sd2 - p1*s1 + p2*y3*s2 + p2*y2*s3 + 2*p3*y2*s2;
442 rs3 = s1 + s2 + s3;
443
444 switch (is) {
445 case 0:
446 rs1 += y1;
447 rs2 -= y1;
448 break;
449 case 1:
450 rs1 -= y2*y3;
451 rs2 += y2*y3;
452 break;
453 case 2:
454 rs2 += y2*y2;
455 break;
456 }
457
458 Ith(resvalS[is],1) = rs1;
459 Ith(resvalS[is],2) = rs2;
460 Ith(resvalS[is],3) = rs3;
461
462 }
463
464 return(0);
465 }
466
rhsQ(realtype t,N_Vector y,N_Vector yp,N_Vector ypQ,void * user_data)467 static int rhsQ(realtype t, N_Vector y, N_Vector yp,
468 N_Vector ypQ, void* user_data)
469 {
470 UserData data;
471
472 data = (UserData) user_data;
473
474 Ith(ypQ,1) = Ith(y,3);
475
476 Ith(ypQ,2) = data->coef*( Ith(y,1)*Ith(y,1)+
477 Ith(y,2)*Ith(y,2)+
478 Ith(y,3)*Ith(y,3) );
479
480 return(0);
481 }
482
483
484 /*
485 *--------------------------------------------------------------------
486 * PRIVATE FUNCTIONS
487 *--------------------------------------------------------------------
488 */
489
490 /*
491 * Process and verify arguments to idasfwddenx.
492 */
493
ProcessArgs(int argc,char * argv[],booleantype * sensi,int * sensi_meth,booleantype * err_con)494 static void ProcessArgs(int argc, char *argv[],
495 booleantype *sensi, int *sensi_meth, booleantype *err_con)
496 {
497 *sensi = SUNFALSE;
498 *sensi_meth = -1;
499 *err_con = SUNFALSE;
500
501 if (argc < 2) WrongArgs(argv[0]);
502
503 if (strcmp(argv[1],"-nosensi") == 0)
504 *sensi = SUNFALSE;
505 else if (strcmp(argv[1],"-sensi") == 0)
506 *sensi = SUNTRUE;
507 else
508 WrongArgs(argv[0]);
509
510 if (*sensi) {
511
512 if (argc != 4)
513 WrongArgs(argv[0]);
514
515 if (strcmp(argv[2],"sim") == 0)
516 *sensi_meth = IDA_SIMULTANEOUS;
517 else if (strcmp(argv[2],"stg") == 0)
518 *sensi_meth = IDA_STAGGERED;
519 else
520 WrongArgs(argv[0]);
521
522 if (strcmp(argv[3],"t") == 0)
523 *err_con = SUNTRUE;
524 else if (strcmp(argv[3],"f") == 0)
525 *err_con = SUNFALSE;
526 else
527 WrongArgs(argv[0]);
528 }
529
530 }
531
WrongArgs(char * name)532 static void WrongArgs(char *name)
533 {
534 printf("\nUsage: %s [-nosensi] [-sensi sensi_meth err_con]\n",name);
535 printf(" sensi_meth = sim or stg\n");
536 printf(" err_con = t or f\n");
537
538 exit(0);
539 }
540
541
PrintIC(N_Vector y,N_Vector yp)542 static void PrintIC(N_Vector y, N_Vector yp)
543 {
544 realtype* data;
545
546 data = N_VGetArrayPointer(y);
547 printf("\n\nConsistent IC:\n");
548 printf("\ty = ");
549 #if defined(SUNDIALS_EXTENDED_PRECISION)
550 printf("%12.4Le %12.4Le %12.4Le \n", data[0], data[1], data[2]);
551 #elif defined(SUNDIALS_DOUBLE_PRECISION)
552 printf("%12.4e %12.4e %12.4e \n", data[0], data[1], data[2]);
553 #else
554 printf("%12.4e %12.4e %12.4e \n", data[0], data[1], data[2]);
555 #endif
556
557 data = N_VGetArrayPointer(yp);
558 printf("\typ= ");
559 #if defined(SUNDIALS_EXTENDED_PRECISION)
560 printf("%12.4Le %12.4Le %12.4Le \n", data[0], data[1], data[2]);
561 #elif defined(SUNDIALS_DOUBLE_PRECISION)
562 printf("%12.4e %12.4e %12.4e \n", data[0], data[1], data[2]);
563 #else
564 printf("%12.4e %12.4e %12.4e \n", data[0], data[1], data[2]);
565 #endif
566
567 }
PrintSensIC(N_Vector y,N_Vector yp,N_Vector * yS,N_Vector * ypS)568 static void PrintSensIC(N_Vector y, N_Vector yp, N_Vector* yS, N_Vector* ypS)
569 {
570 realtype *sdata;
571
572 sdata = N_VGetArrayPointer(yS[0]);
573 printf(" Sensitivity 1 ");
574
575 printf("\n\ts1 = ");
576 #if defined(SUNDIALS_EXTENDED_PRECISION)
577 printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
578 #elif defined(SUNDIALS_DOUBLE_PRECISION)
579 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
580 #else
581 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
582 #endif
583 sdata = N_VGetArrayPointer(ypS[0]);
584 printf("\ts1'= ");
585 #if defined(SUNDIALS_EXTENDED_PRECISION)
586 printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
587 #elif defined(SUNDIALS_DOUBLE_PRECISION)
588 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
589 #else
590 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
591 #endif
592
593
594 printf(" Sensitivity 2 ");
595 sdata = N_VGetArrayPointer(yS[1]);
596 printf("\n\ts2 = ");
597 #if defined(SUNDIALS_EXTENDED_PRECISION)
598 printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
599 #elif defined(SUNDIALS_DOUBLE_PRECISION)
600 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
601 #else
602 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
603 #endif
604 sdata = N_VGetArrayPointer(ypS[1]);
605 printf("\ts2'= ");
606 #if defined(SUNDIALS_EXTENDED_PRECISION)
607 printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
608 #elif defined(SUNDIALS_DOUBLE_PRECISION)
609 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
610 #else
611 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
612 #endif
613
614
615 printf(" Sensitivity 3 ");
616 sdata = N_VGetArrayPointer(yS[2]);
617 printf("\n\ts3 = ");
618 #if defined(SUNDIALS_EXTENDED_PRECISION)
619 printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
620 #elif defined(SUNDIALS_DOUBLE_PRECISION)
621 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
622 #else
623 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
624 #endif
625 sdata = N_VGetArrayPointer(ypS[2]);
626 printf("\ts3'= ");
627 #if defined(SUNDIALS_EXTENDED_PRECISION)
628 printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
629 #elif defined(SUNDIALS_DOUBLE_PRECISION)
630 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
631 #else
632 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
633 #endif
634
635
636 }
637
638 /*
639 * Print current t, step count, order, stepsize, and solution.
640 */
641
PrintOutput(void * ida_mem,realtype t,N_Vector u)642 static void PrintOutput(void *ida_mem, realtype t, N_Vector u)
643 {
644 long int nst;
645 int qu, retval;
646 realtype hu, *udata;
647
648 udata = N_VGetArrayPointer(u);
649
650 retval = IDAGetNumSteps(ida_mem, &nst);
651 check_retval(&retval, "IDAGetNumSteps", 1);
652 retval = IDAGetLastOrder(ida_mem, &qu);
653 check_retval(&retval, "IDAGetLastOrder", 1);
654 retval = IDAGetLastStep(ida_mem, &hu);
655 check_retval(&retval, "IDAGetLastStep", 1);
656
657 #if defined(SUNDIALS_EXTENDED_PRECISION)
658 printf("%8.3Le %2d %8.3Le %5ld\n", t, qu, hu, nst);
659 #elif defined(SUNDIALS_DOUBLE_PRECISION)
660 printf("%8.3e %2d %8.3e %5ld\n", t, qu, hu, nst);
661 #else
662 printf("%8.3e %2d %8.3e %5ld\n", t, qu, hu, nst);
663 #endif
664
665 printf(" Solution ");
666
667 #if defined(SUNDIALS_EXTENDED_PRECISION)
668 printf("%12.4Le %12.4Le %12.4Le \n", udata[0], udata[1], udata[2]);
669 #elif defined(SUNDIALS_DOUBLE_PRECISION)
670 printf("%12.4e %12.4e %12.4e \n", udata[0], udata[1], udata[2]);
671 #else
672 printf("%12.4e %12.4e %12.4e \n", udata[0], udata[1], udata[2]);
673 #endif
674
675 }
676
677 /*
678 * Print sensitivities.
679 */
680
PrintSensOutput(N_Vector * uS)681 static void PrintSensOutput(N_Vector *uS)
682 {
683 realtype *sdata;
684
685 sdata = N_VGetArrayPointer(uS[0]);
686 printf(" Sensitivity 1 ");
687
688 #if defined(SUNDIALS_EXTENDED_PRECISION)
689 printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
690 #elif defined(SUNDIALS_DOUBLE_PRECISION)
691 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
692 #else
693 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
694 #endif
695
696 sdata = N_VGetArrayPointer(uS[1]);
697 printf(" Sensitivity 2 ");
698
699 #if defined(SUNDIALS_EXTENDED_PRECISION)
700 printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
701 #elif defined(SUNDIALS_DOUBLE_PRECISION)
702 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
703 #else
704 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
705 #endif
706
707 sdata = N_VGetArrayPointer(uS[2]);
708 printf(" Sensitivity 3 ");
709
710 #if defined(SUNDIALS_EXTENDED_PRECISION)
711 printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
712 #elif defined(SUNDIALS_DOUBLE_PRECISION)
713 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
714 #else
715 printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
716 #endif
717 }
718
719 /*
720 * Print some final statistics from the IDAS memory.
721 */
722
PrintFinalStats(void * ida_mem,booleantype sensi)723 static void PrintFinalStats(void *ida_mem, booleantype sensi)
724 {
725 long int nst;
726 long int nfe, nsetups, nni, ncfn, netf;
727 long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;
728 long int nje, nfeLS;
729 int retval;
730
731 retval = IDAGetNumSteps(ida_mem, &nst);
732 check_retval(&retval, "IDAGetNumSteps", 1);
733 retval = IDAGetNumResEvals(ida_mem, &nfe);
734 check_retval(&retval, "IDAGetNumRhsEvals", 1);
735 retval = IDAGetNumLinSolvSetups(ida_mem, &nsetups);
736 check_retval(&retval, "IDAGetNumLinSolvSetups", 1);
737 retval = IDAGetNumErrTestFails(ida_mem, &netf);
738 check_retval(&retval, "IDAGetNumErrTestFails", 1);
739 retval = IDAGetNumNonlinSolvIters(ida_mem, &nni);
740 check_retval(&retval, "IDAGetNumNonlinSolvIters", 1);
741 retval = IDAGetNumNonlinSolvConvFails(ida_mem, &ncfn);
742 check_retval(&retval, "IDAGetNumNonlinSolvConvFails", 1);
743
744 if (sensi) {
745 retval = IDAGetSensNumResEvals(ida_mem, &nfSe);
746 check_retval(&retval, "IDAGetSensNumRhsEvals", 1);
747 retval = IDAGetNumResEvalsSens(ida_mem, &nfeS);
748 check_retval(&retval, "IDAGetNumResEvalsSens", 1);
749 retval = IDAGetSensNumLinSolvSetups(ida_mem, &nsetupsS);
750 check_retval(&retval, "IDAGetSensNumLinSolvSetups", 1);
751 retval = IDAGetSensNumErrTestFails(ida_mem, &netfS);
752 check_retval(&retval, "IDAGetSensNumErrTestFails", 1);
753 retval = IDAGetSensNumNonlinSolvIters(ida_mem, &nniS);
754 check_retval(&retval, "IDAGetSensNumNonlinSolvIters", 1);
755 retval = IDAGetSensNumNonlinSolvConvFails(ida_mem, &ncfnS);
756 check_retval(&retval, "IDAGetSensNumNonlinSolvConvFails", 1);
757 }
758
759 retval = IDAGetNumJacEvals(ida_mem, &nje);
760 check_retval(&retval, "IDAGetNumJacEvals", 1);
761 retval = IDAGetNumLinResEvals(ida_mem, &nfeLS);
762 check_retval(&retval, "IDAGetNumLinResEvals", 1);
763
764 printf("\nFinal Statistics\n\n");
765 printf("nst = %5ld\n\n", nst);
766 printf("nfe = %5ld\n", nfe);
767 printf("netf = %5ld nsetups = %5ld\n", netf, nsetups);
768 printf("nni = %5ld ncfn = %5ld\n", nni, ncfn);
769
770 if(sensi) {
771 printf("\n");
772 printf("nfSe = %5ld nfeS = %5ld\n", nfSe, nfeS);
773 printf("netfs = %5ld nsetupsS = %5ld\n", netfS, nsetupsS);
774 printf("nniS = %5ld ncfnS = %5ld\n", nniS, ncfnS);
775 }
776
777 printf("\n");
778 printf("nje = %5ld nfeLS = %5ld\n", nje, nfeLS);
779
780 }
781
782 /*
783 * Check function return value.
784 * opt == 0 means SUNDIALS function allocates memory so check if
785 * returned NULL pointer
786 * opt == 1 means SUNDIALS function returns an integer value so check if
787 * retval < 0
788 * opt == 2 means function allocates memory so check if returned
789 * NULL pointer
790 */
791
check_retval(void * returnvalue,const char * funcname,int opt)792 static int check_retval(void *returnvalue, const char *funcname, int opt)
793 {
794 int *retval;
795
796 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
797 if (opt == 0 && returnvalue == NULL) {
798 fprintf(stderr,
799 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
800 funcname);
801 return(1); }
802
803 /* Check if retval < 0 */
804 else if (opt == 1) {
805 retval = (int *) returnvalue;
806 if (*retval < 0) {
807 fprintf(stderr,
808 "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
809 funcname, *retval);
810 return(1); }}
811
812 /* Check if function returned NULL pointer - no memory allocated */
813 else if (opt == 2 && returnvalue == NULL) {
814 fprintf(stderr,
815 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
816 funcname);
817 return(1); }
818
819 return(0);
820 }
821