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