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