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