1 /* -----------------------------------------------------------------
2 * Programmer(s): Allan Taylor, Alan Hindmarsh and
3 * Radu Serban @ LLNL
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 * This simple example problem for IDA, due to Robertson,
16 * is from chemical kinetics, and consists of the following three
17 * equations:
18 *
19 * dy1/dt = -.04*y1 + 1.e4*y2*y3
20 * dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
21 * 0 = y1 + y2 + y3 - 1
22 *
23 * on the interval from t = 0.0 to t = 4.e10, with initial
24 * conditions: y1 = 1, y2 = y3 = 0.
25 *
26 * While integrating the system, we also use the rootfinding
27 * feature to find the points at which y1 = 1e-4 or at which
28 * y3 = 0.01.
29 *
30 * The problem is solved with IDA using the DENSE linear
31 * solver, with a user-supplied Jacobian. Output is printed at
32 * t = .4, 4, 40, ..., 4e10.
33 * -----------------------------------------------------------------*/
34
35 #include <stdio.h>
36 #include <math.h>
37
38 #include <ida/ida.h> /* prototypes for IDA fcts., consts. */
39 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
40 #include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
41 #include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
42 #include <sunnonlinsol/sunnonlinsol_newton.h> /* access to Newton SUNNonlinearSolver */
43 #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
44 #include <sundials/sundials_math.h> /* defs. of SUNRabs, SUNRexp, etc. */
45
46 #if defined(SUNDIALS_EXTENDED_PRECISION)
47 #define GSYM "Lg"
48 #define ESYM "Le"
49 #define FSYM "Lf"
50 #else
51 #define GSYM "g"
52 #define ESYM "e"
53 #define FSYM "f"
54 #endif
55
56 /* Problem Constants */
57
58 #define NEQ 3
59 #define NOUT 12
60
61 #define ZERO RCONST(0.0)
62 #define ONE RCONST(1.0)
63
64 /* Macro to define dense matrix elements, indexed from 1. */
65
66 #define IJth(A,i,j) SM_ELEMENT_D(A,i-1,j-1)
67
68 /* Prototypes of functions called by IDA */
69
70 int resrob(realtype tres, N_Vector yy, N_Vector yp,
71 N_Vector resval, void *user_data);
72
73 static int grob(realtype t, N_Vector yy, N_Vector yp,
74 realtype *gout, void *user_data);
75
76 int jacrob(realtype tt, realtype cj,
77 N_Vector yy, N_Vector yp, N_Vector resvec,
78 SUNMatrix JJ, void *user_data,
79 N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
80
81 /* Prototypes of private functions */
82 static void PrintHeader(realtype rtol, N_Vector avtol, N_Vector y);
83 static void PrintOutput(void *mem, realtype t, N_Vector y);
84 static void PrintRootInfo(int root_f1, int root_f2);
85 static void PrintFinalStats(void *mem);
86 static int check_retval(void *returnvalue, const char *funcname, int opt);
87 static int check_ans(N_Vector y, realtype t, realtype rtol, N_Vector atol);
88
89 /*
90 *--------------------------------------------------------------------
91 * Main Program
92 *--------------------------------------------------------------------
93 */
94
main(void)95 int main(void)
96 {
97 void *mem;
98 N_Vector yy, yp, avtol;
99 realtype rtol, *yval, *ypval, *atval;
100 realtype t0, tout1, tout, tret;
101 int iout, retval, retvalr;
102 int rootsfound[2];
103 SUNMatrix A;
104 SUNLinearSolver LS;
105 SUNNonlinearSolver NLS;
106
107 mem = NULL;
108 yy = yp = avtol = NULL;
109 yval = ypval = atval = NULL;
110 A = NULL;
111 LS = NULL;
112 NLS = NULL;
113
114 /* Allocate N-vectors. */
115 yy = N_VNew_Serial(NEQ);
116 if(check_retval((void *)yy, "N_VNew_Serial", 0)) return(1);
117 yp = N_VNew_Serial(NEQ);
118 if(check_retval((void *)yp, "N_VNew_Serial", 0)) return(1);
119 avtol = N_VNew_Serial(NEQ);
120 if(check_retval((void *)avtol, "N_VNew_Serial", 0)) return(1);
121
122 /* Create and initialize y, y', and absolute tolerance vectors. */
123 yval = N_VGetArrayPointer(yy);
124 yval[0] = ONE;
125 yval[1] = ZERO;
126 yval[2] = ZERO;
127
128 ypval = N_VGetArrayPointer(yp);
129 ypval[0] = RCONST(-0.04);
130 ypval[1] = RCONST(0.04);
131 ypval[2] = ZERO;
132
133 rtol = RCONST(1.0e-4);
134
135 atval = N_VGetArrayPointer(avtol);
136 atval[0] = RCONST(1.0e-8);
137 atval[1] = RCONST(1.0e-6);
138 atval[2] = RCONST(1.0e-6);
139
140 /* Integration limits */
141 t0 = ZERO;
142 tout1 = RCONST(0.4);
143
144 PrintHeader(rtol, avtol, yy);
145
146 /* Call IDACreate and IDAInit to initialize IDA memory */
147 mem = IDACreate();
148 if(check_retval((void *)mem, "IDACreate", 0)) return(1);
149 retval = IDAInit(mem, resrob, t0, yy, yp);
150 if(check_retval(&retval, "IDAInit", 1)) return(1);
151 /* Call IDASVtolerances to set tolerances */
152 retval = IDASVtolerances(mem, rtol, avtol);
153 if(check_retval(&retval, "IDASVtolerances", 1)) return(1);
154
155 /* Call IDARootInit to specify the root function grob with 2 components */
156 retval = IDARootInit(mem, 2, grob);
157 if (check_retval(&retval, "IDARootInit", 1)) return(1);
158
159 /* Create dense SUNMatrix for use in linear solves */
160 A = SUNDenseMatrix(NEQ, NEQ);
161 if(check_retval((void *)A, "SUNDenseMatrix", 0)) return(1);
162
163 /* Create dense SUNLinearSolver object */
164 LS = SUNLinSol_Dense(yy, A);
165 if(check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1);
166
167 /* Attach the matrix and linear solver */
168 retval = IDASetLinearSolver(mem, LS, A);
169 if(check_retval(&retval, "IDASetLinearSolver", 1)) return(1);
170
171 /* Set the user-supplied Jacobian routine */
172 retval = IDASetJacFn(mem, jacrob);
173 if(check_retval(&retval, "IDASetJacFn", 1)) return(1);
174
175 /* Create Newton SUNNonlinearSolver object. IDA uses a
176 * Newton SUNNonlinearSolver by default, so it is unecessary
177 * to create it and attach it. It is done in this example code
178 * solely for demonstration purposes. */
179 NLS = SUNNonlinSol_Newton(yy);
180 if(check_retval((void *)NLS, "SUNNonlinSol_Newton", 0)) return(1);
181
182 /* Attach the nonlinear solver */
183 retval = IDASetNonlinearSolver(mem, NLS);
184 if(check_retval(&retval, "IDASetNonlinearSolver", 1)) return(1);
185
186 /* In loop, call IDASolve, print results, and test for error.
187 Break out of loop when NOUT preset output times have been reached. */
188
189 iout = 0; tout = tout1;
190 while(1) {
191
192 retval = IDASolve(mem, tout, &tret, yy, yp, IDA_NORMAL);
193
194 PrintOutput(mem,tret,yy);
195
196 if(check_retval(&retval, "IDASolve", 1)) return(1);
197
198 if (retval == IDA_ROOT_RETURN) {
199 retvalr = IDAGetRootInfo(mem, rootsfound);
200 check_retval(&retvalr, "IDAGetRootInfo", 1);
201 PrintRootInfo(rootsfound[0],rootsfound[1]);
202 }
203
204 if (retval == IDA_SUCCESS) {
205 iout++;
206 tout *= RCONST(10.0);
207 }
208
209 if (iout == NOUT) break;
210 }
211
212 PrintFinalStats(mem);
213
214 /* check the solution error */
215 retval = check_ans(yy, tret, rtol, avtol);
216
217 /* Free memory */
218 IDAFree(&mem);
219 SUNNonlinSolFree(NLS);
220 SUNLinSolFree(LS);
221 SUNMatDestroy(A);
222 N_VDestroy(avtol);
223 N_VDestroy(yy);
224 N_VDestroy(yp);
225
226 return(retval);
227
228 }
229
230 /*
231 *--------------------------------------------------------------------
232 * Functions called by IDA
233 *--------------------------------------------------------------------
234 */
235
236 /*
237 * Define the system residual function.
238 */
239
resrob(realtype tres,N_Vector yy,N_Vector yp,N_Vector rr,void * user_data)240 int resrob(realtype tres, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
241 {
242 realtype *yval, *ypval, *rval;
243
244 yval = N_VGetArrayPointer(yy);
245 ypval = N_VGetArrayPointer(yp);
246 rval = N_VGetArrayPointer(rr);
247
248 rval[0] = RCONST(-0.04)*yval[0] + RCONST(1.0e4)*yval[1]*yval[2];
249 rval[1] = -rval[0] - RCONST(3.0e7)*yval[1]*yval[1] - ypval[1];
250 rval[0] -= ypval[0];
251 rval[2] = yval[0] + yval[1] + yval[2] - ONE;
252
253 return(0);
254 }
255
256 /*
257 * Root function routine. Compute functions g_i(t,y) for i = 0,1.
258 */
259
grob(realtype t,N_Vector yy,N_Vector yp,realtype * gout,void * user_data)260 static int grob(realtype t, N_Vector yy, N_Vector yp, realtype *gout,
261 void *user_data)
262 {
263 realtype *yval, y1, y3;
264
265 yval = N_VGetArrayPointer(yy);
266 y1 = yval[0]; y3 = yval[2];
267 gout[0] = y1 - RCONST(0.0001);
268 gout[1] = y3 - RCONST(0.01);
269
270 return(0);
271 }
272
273 /*
274 * Define the Jacobian function.
275 */
276
jacrob(realtype tt,realtype cj,N_Vector yy,N_Vector yp,N_Vector resvec,SUNMatrix JJ,void * user_data,N_Vector tempv1,N_Vector tempv2,N_Vector tempv3)277 int jacrob(realtype tt, realtype cj,
278 N_Vector yy, N_Vector yp, N_Vector resvec,
279 SUNMatrix JJ, void *user_data,
280 N_Vector tempv1, N_Vector tempv2, N_Vector tempv3)
281 {
282 realtype *yval;
283
284 yval = N_VGetArrayPointer(yy);
285
286 IJth(JJ,1,1) = RCONST(-0.04) - cj;
287 IJth(JJ,2,1) = RCONST(0.04);
288 IJth(JJ,3,1) = ONE;
289 IJth(JJ,1,2) = RCONST(1.0e4)*yval[2];
290 IJth(JJ,2,2) = RCONST(-1.0e4)*yval[2] - RCONST(6.0e7)*yval[1] - cj;
291 IJth(JJ,3,2) = ONE;
292 IJth(JJ,1,3) = RCONST(1.0e4)*yval[1];
293 IJth(JJ,2,3) = RCONST(-1.0e4)*yval[1];
294 IJth(JJ,3,3) = ONE;
295
296 return(0);
297 }
298
299 /*
300 *--------------------------------------------------------------------
301 * Private functions
302 *--------------------------------------------------------------------
303 */
304
305 /*
306 * Print first lines of output (problem description)
307 */
308
PrintHeader(realtype rtol,N_Vector avtol,N_Vector y)309 static void PrintHeader(realtype rtol, N_Vector avtol, N_Vector y)
310 {
311 realtype *atval, *yval;
312
313 atval = N_VGetArrayPointer(avtol);
314 yval = N_VGetArrayPointer(y);
315
316 printf("\nidaRoberts_dns: Robertson kinetics DAE serial example problem for IDA\n");
317 printf(" Three equation chemical kinetics problem.\n\n");
318 printf("Linear solver: DENSE, with user-supplied Jacobian.\n");
319 #if defined(SUNDIALS_EXTENDED_PRECISION)
320 printf("Tolerance parameters: rtol = %Lg atol = %Lg %Lg %Lg \n",
321 rtol, atval[0],atval[1],atval[2]);
322 printf("Initial conditions y0 = (%Lg %Lg %Lg)\n",
323 yval[0], yval[1], yval[2]);
324 #elif defined(SUNDIALS_DOUBLE_PRECISION)
325 printf("Tolerance parameters: rtol = %g atol = %g %g %g \n",
326 rtol, atval[0],atval[1],atval[2]);
327 printf("Initial conditions y0 = (%g %g %g)\n",
328 yval[0], yval[1], yval[2]);
329 #else
330 printf("Tolerance parameters: rtol = %g atol = %g %g %g \n",
331 rtol, atval[0],atval[1],atval[2]);
332 printf("Initial conditions y0 = (%g %g %g)\n",
333 yval[0], yval[1], yval[2]);
334 #endif
335 printf("Constraints and id not used.\n\n");
336 printf("-----------------------------------------------------------------------\n");
337 printf(" t y1 y2 y3");
338 printf(" | nst k h\n");
339 printf("-----------------------------------------------------------------------\n");
340 }
341
342 /*
343 * Print Output
344 */
345
PrintOutput(void * mem,realtype t,N_Vector y)346 static void PrintOutput(void *mem, realtype t, N_Vector y)
347 {
348 realtype *yval;
349 int retval, kused;
350 long int nst;
351 realtype hused;
352
353 yval = N_VGetArrayPointer(y);
354
355 retval = IDAGetLastOrder(mem, &kused);
356 check_retval(&retval, "IDAGetLastOrder", 1);
357 retval = IDAGetNumSteps(mem, &nst);
358 check_retval(&retval, "IDAGetNumSteps", 1);
359 retval = IDAGetLastStep(mem, &hused);
360 check_retval(&retval, "IDAGetLastStep", 1);
361 #if defined(SUNDIALS_EXTENDED_PRECISION)
362 printf("%10.4Le %12.4Le %12.4Le %12.4Le | %3ld %1d %12.4Le\n",
363 t, yval[0], yval[1], yval[2], nst, kused, hused);
364 #elif defined(SUNDIALS_DOUBLE_PRECISION)
365 printf("%10.4e %12.4e %12.4e %12.4e | %3ld %1d %12.4e\n",
366 t, yval[0], yval[1], yval[2], nst, kused, hused);
367 #else
368 printf("%10.4e %12.4e %12.4e %12.4e | %3ld %1d %12.4e\n",
369 t, yval[0], yval[1], yval[2], nst, kused, hused);
370 #endif
371 }
372
PrintRootInfo(int root_f1,int root_f2)373 static void PrintRootInfo(int root_f1, int root_f2)
374 {
375 printf(" rootsfound[] = %3d %3d\n", root_f1, root_f2);
376 return;
377 }
378
379 /*
380 * Print final integrator statistics
381 */
382
PrintFinalStats(void * mem)383 static void PrintFinalStats(void *mem)
384 {
385 int retval;
386 long int nst, nni, nje, nre, nreLS, netf, ncfn, nge;
387
388 retval = IDAGetNumSteps(mem, &nst);
389 check_retval(&retval, "IDAGetNumSteps", 1);
390 retval = IDAGetNumResEvals(mem, &nre);
391 check_retval(&retval, "IDAGetNumResEvals", 1);
392 retval = IDAGetNumJacEvals(mem, &nje);
393 check_retval(&retval, "IDAGetNumJacEvals", 1);
394 retval = IDAGetNumNonlinSolvIters(mem, &nni);
395 check_retval(&retval, "IDAGetNumNonlinSolvIters", 1);
396 retval = IDAGetNumErrTestFails(mem, &netf);
397 check_retval(&retval, "IDAGetNumErrTestFails", 1);
398 retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
399 check_retval(&retval, "IDAGetNumNonlinSolvConvFails", 1);
400 retval = IDAGetNumLinResEvals(mem, &nreLS);
401 check_retval(&retval, "IDAGetNumLinResEvals", 1);
402 retval = IDAGetNumGEvals(mem, &nge);
403 check_retval(&retval, "IDAGetNumGEvals", 1);
404
405 printf("\nFinal Run Statistics: \n\n");
406 printf("Number of steps = %ld\n", nst);
407 printf("Number of residual evaluations = %ld\n", nre+nreLS);
408 printf("Number of Jacobian evaluations = %ld\n", nje);
409 printf("Number of nonlinear iterations = %ld\n", nni);
410 printf("Number of error test failures = %ld\n", netf);
411 printf("Number of nonlinear conv. failures = %ld\n", ncfn);
412 printf("Number of root fn. evaluations = %ld\n", nge);
413 }
414
415 /*
416 * Check function return value...
417 * opt == 0 means SUNDIALS function allocates memory so check if
418 * returned NULL pointer
419 * opt == 1 means SUNDIALS function returns an integer value so check if
420 * retval < 0
421 * opt == 2 means function allocates memory so check if returned
422 * NULL pointer
423 */
424
check_retval(void * returnvalue,const char * funcname,int opt)425 static int check_retval(void *returnvalue, const char *funcname, int opt)
426 {
427 int *retval;
428 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
429 if (opt == 0 && returnvalue == NULL) {
430 fprintf(stderr,
431 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
432 funcname);
433 return(1);
434 } else if (opt == 1) {
435 /* Check if retval < 0 */
436 retval = (int *) returnvalue;
437 if (*retval < 0) {
438 fprintf(stderr,
439 "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
440 funcname, *retval);
441 return(1);
442 }
443 } else if (opt == 2 && returnvalue == NULL) {
444 /* Check if function returned NULL pointer - no memory allocated */
445 fprintf(stderr,
446 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
447 funcname);
448 return(1);
449 }
450
451 return(0);
452 }
453
454 /* compare the solution at the final time 4e10s to a reference solution computed
455 using a relative tolerance of 1e-8 and absoltue tolerance of 1e-14 */
check_ans(N_Vector y,realtype t,realtype rtol,N_Vector atol)456 static int check_ans(N_Vector y, realtype t, realtype rtol, N_Vector atol)
457 {
458 int passfail=0; /* answer pass (0) or fail (1) retval */
459 N_Vector ref; /* reference solution vector */
460 N_Vector ewt; /* error weight vector */
461 realtype err; /* wrms error */
462
463 /* create reference solution and error weight vectors */
464 ref = N_VClone(y);
465 ewt = N_VClone(y);
466
467 /* set the reference solution data */
468 NV_Ith_S(ref,0) = RCONST(5.2083474251394888e-08);
469 NV_Ith_S(ref,1) = RCONST(2.0833390772616859e-13);
470 NV_Ith_S(ref,2) = RCONST(9.9999994791631752e-01);
471
472 /* compute the error weight vector, loosen atol */
473 N_VAbs(ref, ewt);
474 N_VLinearSum(rtol, ewt, RCONST(10.0), atol, ewt);
475 if (N_VMin(ewt) <= ZERO) {
476 fprintf(stderr, "\nSUNDIALS_ERROR: check_ans failed - ewt <= 0\n\n");
477 return(-1);
478 }
479 N_VInv(ewt, ewt);
480
481 /* compute the solution error */
482 N_VLinearSum(ONE, y, -ONE, ref, ref);
483 err = N_VWrmsNorm(ref, ewt);
484
485 /* is the solution within the tolerances? */
486 passfail = (err < ONE) ? 0 : 1;
487
488 if (passfail) {
489 fprintf(stdout, "\nSUNDIALS_WARNING: check_ans error=%"GSYM"\n\n", err);
490 }
491
492 /* Free vectors */
493 N_VDestroy(ref);
494 N_VDestroy(ewt);
495
496 return(passfail);
497 }
498