1 /* -----------------------------------------------------------------
2 * Programmer(s): Jimmy Almgren-Bell @ LLNL
3 * Based on prior version by: Scott D. Cohen, Alan C. Hindmarsh, and
4 * Radu Serban @ LLNL
5 * -----------------------------------------------------------------
6 * SUNDIALS Copyright Start
7 * Copyright (c) 2002-2020, 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 coding
19 * needed for its solution by CVODE. The problem is from
20 * chemical kinetics, and consists of the following three rate
21 * equations:
22 * dy1/dt = -.04*y1 + 1.e4*y2*y3
23 * dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2
24 * dy3/dt = 3.e7*(y2)^2
25 * on the interval from t = 0.0 to t = 4.e10, with initial
26 * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.
27 * While integrating the system, we also use the rootfinding
28 * feature to find the points at which y1 = 1e-4 or at which
29 * y3 = 0.01. This program solves the problem with the BDF method,
30 * Newton iteration with the SUNDENSE dense linear solver, and a
31 * user-supplied Jacobian routine. It uses a scalar relative tolerance
32 * and a vector absolute tolerance. The constraint y_i >= 0 is
33 * posed for all components. Output is printed in decades
34 * from t = .4 to t = 4.e10. Run statistics (optional outputs)
35 * are printed at the end.
36 * -----------------------------------------------------------------*/
37
38 #include <stdio.h>
39
40 #include <cvodes/cvodes.h> /* prototypes for CVODE fcts., consts. */
41 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
42 #include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
43 #include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
44 #include <cvodes/cvodes_direct.h> /* access to CVDls interface */
45 #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
46
47 #if defined(SUNDIALS_EXTENDED_PRECISION)
48 #define GSYM "Lg"
49 #define ESYM "Le"
50 #define FSYM "Lf"
51 #else
52 #define GSYM "g"
53 #define ESYM "e"
54 #define FSYM "f"
55 #endif
56
57 /* User-defined vector and matrix accessor macros: Ith, IJth */
58
59 /* These macros are defined in order to write code which exactly matches
60 the mathematical problem description given above.
61
62 Ith(v,i) references the ith component of the vector v, where i is in
63 the range [1..NEQ] and NEQ is defined below. The Ith macro is defined
64 using the N_VIth macro in nvector.h. N_VIth numbers the components of
65 a vector starting from 0.
66
67 IJth(A,i,j) references the (i,j)th element of the dense matrix A, where
68 i and j are in the range [1..NEQ]. The IJth macro is defined using the
69 SM_ELEMENT_D macro in dense.h. SM_ELEMENT_D numbers rows and columns of
70 a dense matrix starting from 0. */
71
72 #define Ith(v,i) NV_Ith_S(v,i-1) /* Ith numbers components 1..NEQ */
73 #define IJth(A,i,j) SM_ELEMENT_D(A,i-1,j-1) /* IJth numbers rows,cols 1..NEQ */
74
75
76 /* Problem Constants */
77
78 #define NEQ 3 /* number of equations */
79 #define Y1 RCONST(1.0) /* initial y components */
80 #define Y2 RCONST(0.0)
81 #define Y3 RCONST(0.0)
82 #define RTOL RCONST(1.0e-4) /* scalar relative tolerance */
83 #define ATOL1 RCONST(1.0e-6) /* vector absolute tolerance components */
84 #define ATOL2 RCONST(1.0e-11)
85 #define ATOL3 RCONST(1.0e-5)
86 #define T0 RCONST(0.0) /* initial time */
87 #define T1 RCONST(0.4) /* first output time */
88 #define TMULT RCONST(10.0) /* output time factor */
89 #define NOUT 12 /* number of output times */
90
91 #define ZERO RCONST(0.0)
92 #define ONE RCONST(1.0)
93
94 /* Functions Called by the Solver */
95
96 static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
97
98 static int g(realtype t, N_Vector y, realtype *gout, void *user_data);
99
100 static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
101 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
102
103 /* Private functions to output results */
104
105 static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3);
106 static void PrintRootInfo(int root_f1, int root_f2);
107
108 /* Private function to print final statistics */
109
110 static void PrintFinalStats(void *cvode_mem);
111
112 /* Private function to check function return values */
113
114 static int check_retval(void *returnvalue, const char *funcname, int opt);
115
116 /* Private function to check computed solution */
117
118 static int check_ans(N_Vector y, realtype t, realtype rtol, N_Vector atol);
119
120
121 /*
122 *-------------------------------
123 * Main Program
124 *-------------------------------
125 */
126
main()127 int main()
128 {
129 realtype reltol, t, tout;
130 N_Vector y, abstol, constraints;
131 SUNMatrix A;
132 SUNLinearSolver LS;
133 void *cvode_mem;
134 int retval, retvalr, iout;
135 int rootsfound[2];
136
137 y = abstol = NULL;
138 A = NULL;
139 constraints = NULL;
140 LS = NULL;
141 cvode_mem = NULL;
142
143 /* Create serial vector of length NEQ for I.C., abstol, and constraints */
144 y = N_VNew_Serial(NEQ);
145 if (check_retval((void *)y, "N_VNew_Serial", 0)) return(1);
146 abstol = N_VNew_Serial(NEQ);
147 if (check_retval((void *)abstol, "N_VNew_Serial", 0)) return(1);
148
149 constraints = N_VNew_Serial(NEQ);
150 if (check_retval((void *)constraints, "N_VNew_Serial", 0)) return(1);
151
152 /* Initialize y */
153 Ith(y,1) = Y1;
154 Ith(y,2) = Y2;
155 Ith(y,3) = Y3;
156
157 /* Set the scalar relative tolerance */
158 reltol = RTOL;
159 /* Set the vector absolute tolerance */
160 Ith(abstol,1) = ATOL1;
161 Ith(abstol,2) = ATOL2;
162 Ith(abstol,3) = ATOL3;
163
164 /* Set constraints to all 1's for nonnegative solution values. */
165 N_VConst(ONE, constraints);
166
167 /* Call CVodeCreate to create the solver memory and specify the
168 * Backward Differentiation Formula */
169 cvode_mem = CVodeCreate(CV_BDF);
170 if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
171
172 /* Call CVodeInit to initialize the integrator memory and specify the
173 * user's right hand side function in y'=f(t,y), the inital time T0, and
174 * the initial dependent variable vector y. */
175 retval = CVodeInit(cvode_mem, f, T0, y);
176 if (check_retval(&retval, "CVodeInit", 1)) return(1);
177
178 /* Call CVodeSVtolerances to specify the scalar relative tolerance
179 * and vector absolute tolerances */
180 retval = CVodeSVtolerances(cvode_mem, reltol, abstol);
181 if (check_retval(&retval, "CVodeSVtolerances", 1)) return(1);
182
183 /* Call CVodeRootInit to specify the root function g with 2 components */
184 retval = CVodeRootInit(cvode_mem, 2, g);
185 if (check_retval(&retval, "CVodeRootInit", 1)) return(1);
186
187 /* Create dense SUNMatrix for use in linear solves */
188 A = SUNDenseMatrix(NEQ, NEQ);
189 if(check_retval((void *)A, "SUNDenseMatrix", 0)) return(1);
190
191 /* Create dense SUNLinearSolver object for use by CVode */
192 LS = SUNLinSol_Dense(y, A);
193 if(check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1);
194
195 /* Call CVDlsSetLinearSolver to attach the matrix and linear solver to CVode */
196 retval = CVDlsSetLinearSolver(cvode_mem, LS, A);
197 if(check_retval(&retval, "CVDlsSetLinearSolver", 1)) return(1);
198
199 /* Set the user-supplied Jacobian routine Jac */
200 retval = CVDlsSetJacFn(cvode_mem, Jac);
201 if(check_retval(&retval, "CVDlsSetJacFn", 1)) return(1);
202
203 /* Call CVodeSetConstraints to initialize constraints */
204 retval = CVodeSetConstraints(cvode_mem, constraints);
205 if(check_retval(&retval, "CVodeSetConstraints", 1)) return(1);
206 N_VDestroy(constraints);
207
208 /* In loop, call CVode, print results, and test for error.
209 Break out of loop when NOUT preset output times have been reached. */
210 printf(" \n3-species kinetics problem\n\n");
211
212 iout = 0; tout = T1;
213 while(1) {
214 retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
215 PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));
216
217 if (retval == CV_ROOT_RETURN) {
218 retvalr = CVodeGetRootInfo(cvode_mem, rootsfound);
219 if (check_retval(&retvalr, "CVodeGetRootInfo", 1)) return(1);
220 PrintRootInfo(rootsfound[0],rootsfound[1]);
221 }
222
223 if (check_retval(&retval, "CVode", 1)) break;
224 if (retval == CV_SUCCESS) {
225 iout++;
226 tout *= TMULT;
227 }
228
229 if (iout == NOUT) break;
230 }
231
232 /* Print some final statistics */
233 PrintFinalStats(cvode_mem);
234
235 /* check the solution error */
236 retval = check_ans(y, t, reltol, abstol);
237
238 /* Free y and abstol vectors */
239 N_VDestroy(y);
240 N_VDestroy(abstol);
241
242 /* Free integrator memory */
243 CVodeFree(&cvode_mem);
244
245 /* Free the linear solver memory */
246 SUNLinSolFree(LS);
247
248 /* Free the matrix memory */
249 SUNMatDestroy(A);
250
251 return(retval);
252 }
253
254
255 /*
256 *-------------------------------
257 * Functions called by the solver
258 *-------------------------------
259 */
260
261 /*
262 * f routine. Compute function f(t,y).
263 */
264
f(realtype t,N_Vector y,N_Vector ydot,void * user_data)265 static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
266 {
267 realtype y1, y2, y3, yd1, yd3;
268
269 y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
270
271 yd1 = Ith(ydot,1) = RCONST(-0.04)*y1 + RCONST(1.0e4)*y2*y3;
272 yd3 = Ith(ydot,3) = RCONST(3.0e7)*y2*y2;
273 Ith(ydot,2) = -yd1 - yd3;
274
275 return(0);
276 }
277
278 /*
279 * g routine. Compute functions g_i(t,y) for i = 0,1.
280 */
281
g(realtype t,N_Vector y,realtype * gout,void * user_data)282 static int g(realtype t, N_Vector y, realtype *gout, void *user_data)
283 {
284 realtype y1, y3;
285
286 y1 = Ith(y,1); y3 = Ith(y,3);
287 gout[0] = y1 - RCONST(0.0001);
288 gout[1] = y3 - RCONST(0.01);
289
290 return(0);
291 }
292
293 /*
294 * Jacobian routine. Compute J(t,y) = df/dy. *
295 */
296
Jac(realtype t,N_Vector y,N_Vector fy,SUNMatrix J,void * user_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)297 static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
298 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
299 {
300 realtype y2, y3;
301
302 y2 = Ith(y,2); y3 = Ith(y,3);
303
304 IJth(J,1,1) = RCONST(-0.04);
305 IJth(J,1,2) = RCONST(1.0e4)*y3;
306 IJth(J,1,3) = RCONST(1.0e4)*y2;
307
308 IJth(J,2,1) = RCONST(0.04);
309 IJth(J,2,2) = RCONST(-1.0e4)*y3-RCONST(6.0e7)*y2;
310 IJth(J,2,3) = RCONST(-1.0e4)*y2;
311
312 IJth(J,3,1) = ZERO;
313 IJth(J,3,2) = RCONST(6.0e7)*y2;
314 IJth(J,3,3) = ZERO;
315
316 return(0);
317 }
318
319 /*
320 *-------------------------------
321 * Private helper functions
322 *-------------------------------
323 */
324
PrintOutput(realtype t,realtype y1,realtype y2,realtype y3)325 static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3)
326 {
327 #if defined(SUNDIALS_EXTENDED_PRECISION)
328 printf("At t = %0.4Le y =%14.6Le %14.6Le %14.6Le\n", t, y1, y2, y3);
329 #elif defined(SUNDIALS_DOUBLE_PRECISION)
330 printf("At t = %0.4e y =%14.6e %14.6e %14.6e\n", t, y1, y2, y3);
331 #else
332 printf("At t = %0.4e y =%14.6e %14.6e %14.6e\n", t, y1, y2, y3);
333 #endif
334
335 return;
336 }
337
PrintRootInfo(int root_f1,int root_f2)338 static void PrintRootInfo(int root_f1, int root_f2)
339 {
340 printf(" rootsfound[] = %3d %3d\n", root_f1, root_f2);
341
342 return;
343 }
344
345 /*
346 * Get and print some final statistics
347 */
348
PrintFinalStats(void * cvode_mem)349 static void PrintFinalStats(void *cvode_mem)
350 {
351 long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
352 int retval;
353
354 retval = CVodeGetNumSteps(cvode_mem, &nst);
355 check_retval(&retval, "CVodeGetNumSteps", 1);
356 retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);
357 check_retval(&retval, "CVodeGetNumRhsEvals", 1);
358 retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
359 check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);
360 retval = CVodeGetNumErrTestFails(cvode_mem, &netf);
361 check_retval(&retval, "CVodeGetNumErrTestFails", 1);
362 retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
363 check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);
364 retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
365 check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);
366
367 retval = CVDlsGetNumJacEvals(cvode_mem, &nje);
368 check_retval(&retval, "CVDlsGetNumJacEvals", 1);
369 retval = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);
370 check_retval(&retval, "CVDlsGetNumRhsEvals", 1);
371
372 retval = CVodeGetNumGEvals(cvode_mem, &nge);
373 check_retval(&retval, "CVodeGetNumGEvals", 1);
374
375 printf("\nFinal Statistics:\n");
376 printf("nst = %-6ld nfe = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld\n",
377 nst, nfe, nsetups, nfeLS, nje);
378 printf("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n \n",
379 nni, ncfn, netf, nge);
380 }
381
382 /*
383 * Check function return value...
384 * opt == 0 means SUNDIALS function allocates memory so check if
385 * returned NULL pointer
386 * opt == 1 means SUNDIALS function returns an integer value so check if
387 * retval < 0
388 * opt == 2 means function allocates memory so check if returned
389 * NULL pointer
390 */
391
check_retval(void * returnvalue,const char * funcname,int opt)392 static int check_retval(void *returnvalue, const char *funcname, int opt)
393 {
394 int *retval;
395
396 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
397 if (opt == 0 && returnvalue == NULL) {
398 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
399 funcname);
400 return(1); }
401
402 /* Check if retval < 0 */
403 else if (opt == 1) {
404 retval = (int *) returnvalue;
405 if (*retval < 0) {
406 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
407 funcname, *retval);
408 return(1); }}
409
410 /* Check if function returned NULL pointer - no memory allocated */
411 else if (opt == 2 && returnvalue == NULL) {
412 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
413 funcname);
414 return(1); }
415
416 return(0);
417 }
418
419 /* compare the solution at the final time 4e10s to a reference solution computed
420 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)421 static int check_ans(N_Vector y, realtype t, realtype rtol, N_Vector atol)
422 {
423 int passfail=0; /* answer pass (0) or fail (1) retval */
424 N_Vector ref; /* reference solution vector */
425 N_Vector ewt; /* error weight vector */
426 realtype err; /* wrms error */
427
428 /* create reference solution and error weight vectors */
429 ref = N_VClone(y);
430 ewt = N_VClone(y);
431
432 /* set the reference solution data */
433 NV_Ith_S(ref,0) = RCONST(5.2083495894337328e-08);
434 NV_Ith_S(ref,1) = RCONST(2.0833399429795671e-13);
435 NV_Ith_S(ref,2) = RCONST(9.9999994791629776e-01);
436
437 /* compute the error weight vector, loosen atol */
438 N_VAbs(ref, ewt);
439 N_VLinearSum(rtol, ewt, RCONST(10.0), atol, ewt);
440 if (N_VMin(ewt) <= ZERO) {
441 fprintf(stderr, "\nSUNDIALS_ERROR: check_ans failed - ewt <= 0\n\n");
442 return(-1);
443 }
444 N_VInv(ewt, ewt);
445
446 /* compute the solution error */
447 N_VLinearSum(ONE, y, -ONE, ref, ref);
448 err = N_VWrmsNorm(ref, ewt);
449
450 /* is the solution within the tolerances? */
451 passfail = (err < ONE) ? 0 : 1;
452
453 if (passfail) {
454 fprintf(stdout, "\nSUNDIALS_WARNING: check_ans error=%"GSYM"\n\n", err);
455 }
456
457 /* Free vectors */
458 N_VDestroy(ref);
459 N_VDestroy(ewt);
460
461 return(passfail);
462 }
463