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