1 /*---------------------------------------------------------------
2 * Programmer(s): Daniel R. Reynolds @ SMU
3 *---------------------------------------------------------------
4 * SUNDIALS Copyright Start
5 * Copyright (c) 2002-2020, 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 * Example problem:
15 *
16 * The following test simulates the Robertson problem,
17 * corresponding to the kinetics of an autocatalytic reaction.
18 * This is an ODE system with 3 components, Y = [u,v,w], satisfying
19 * the equations,
20 * du/dt = -0.04*u + 1e4*v*w
21 * dv/dt = 0.04*u - 1e4*v*w - 3e7*v^2
22 * dw/dt = 3e7*v^2
23 * for t in the interval [0.0, 1e11], with initial conditions
24 * Y0 = [1,0,0].
25 *
26 * This program solves the problem with one of the solvers, ERK,
27 * DIRK or ARK. For DIRK and ARK, implicit subsystems are solved
28 * using a Newton iteration with the dense SUNLinearSolver, and a
29 * user-supplied Jacobian routine. The constraint y_i >= 0 is
30 * posed for all components.
31 *
32 * 100 outputs are printed at equal intervals, and run statistics
33 * are printed at the end.
34 *---------------------------------------------------------------*/
35
36 /* Header files */
37 #include <stdio.h>
38 #include <math.h>
39 #include <arkode/arkode_arkstep.h> /* prototypes for ARKStep fcts., consts */
40 #include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
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
45 #if defined(SUNDIALS_EXTENDED_PRECISION)
46 #define GSYM "Lg"
47 #define ESYM "Le"
48 #define FSYM "Lf"
49 #else
50 #define GSYM "g"
51 #define ESYM "e"
52 #define FSYM "f"
53 #endif
54
55 /* User-supplied Functions Called by the Solver */
56 static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
57 static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
58 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
59
60 /* Private function to check function return values */
61 static int check_flag(void *flagvalue, const char *funcname, int opt);
62
63 /* Private function to check computed solution */
64 static int check_ans(N_Vector y, realtype t, realtype rtol, realtype atol);
65
66 /* Main Program */
main()67 int main()
68 {
69 realtype ONE = RCONST(1.0);
70
71 /* general problem parameters */
72 realtype T0 = RCONST(0.0); /* initial time */
73 realtype Tf = RCONST(1.e11); /* final time */
74 realtype dTout = (Tf-T0)/100; /* time between outputs */
75 int Nt = (int) ceil(Tf/dTout); /* number of output times */
76 sunindextype NEQ = 3; /* number of dependent vars. */
77
78 /* general problem variables */
79 int flag; /* reusable error-checking flag */
80 N_Vector y = NULL; /* empty vector for storing solution */
81 N_Vector constraints = NULL; /* empty vector for storing constraints */
82 SUNMatrix A = NULL; /* empty matrix for linear solver */
83 SUNLinearSolver LS = NULL; /* empty linear solver object */
84 void *arkode_mem = NULL; /* empty ARKode memory structure */
85 FILE *UFID;
86 realtype t, tout;
87 int iout;
88 long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf, nctf;
89
90 /* set up the initial conditions, tolerances, initial time step size */
91 realtype u0 = RCONST(1.0);
92 realtype v0 = RCONST(0.0);
93 realtype w0 = RCONST(0.0);
94 realtype reltol = 1.e-3;
95 realtype abstol = 1.e-7;
96 realtype h0 = 1.e-4 * reltol;
97
98 /* Initial problem output */
99 printf("\nRobertson ODE test problem:\n");
100 printf(" initial conditions: u0 = %"GSYM", v0 = %"GSYM", w0 = %"GSYM"\n",u0,v0,w0);
101
102 /* Initialize data structures */
103 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
104 if (check_flag((void *) y, "N_VNew_Serial", 0)) return 1;
105 NV_Ith_S(y,0) = u0; /* Set initial conditions into y */
106 NV_Ith_S(y,1) = v0;
107 NV_Ith_S(y,2) = w0;
108
109 constraints = N_VNew_Serial(NEQ);
110 if (check_flag((void *) constraints, "N_VNew_Serial", 0)) return 1;
111 /* Set constraints to all 1's for nonnegative solution values. */
112 N_VConst(ONE, constraints);
113
114 /* Call ARKStepCreate to initialize the ARK timestepper module and
115 specify the right-hand side function in y'=f(t,y), the inital time
116 T0, and the initial dependent variable vector y. Note: since this
117 problem is fully implicit, we set f_E to NULL and f_I to f. */
118 arkode_mem = ARKStepCreate(NULL, f, T0, y);
119 if (check_flag((void *)arkode_mem, "ARKStepCreate", 0)) return 1;
120
121 /* Set routines */
122 flag = ARKStepSetInitStep(arkode_mem, h0); /* Set custom initial step */
123 if (check_flag(&flag, "ARKStepSetInitStep", 1)) return 1;
124 flag = ARKStepSetMaxErrTestFails(arkode_mem, 20); /* Increase max error test fails */
125 if (check_flag(&flag, "ARKStepSetMaxErrTestFails", 1)) return 1;
126 flag = ARKStepSetMaxNonlinIters(arkode_mem, 8); /* Increase max nonlin iters */
127 if (check_flag(&flag, "ARKStepSetMaxNonlinIters", 1)) return 1;
128 flag = ARKStepSetNonlinConvCoef(arkode_mem, 1.e-7); /* Set nonlinear convergence coeff. */
129 if (check_flag(&flag, "ARKStepSetNonlinConvCoef", 1)) return 1;
130 flag = ARKStepSetMaxNumSteps(arkode_mem, 100000); /* Increase max num steps */
131 if (check_flag(&flag, "ARKStepSetMaxNumSteps", 1)) return 1;
132 flag = ARKStepSetPredictorMethod(arkode_mem, 1); /* Specify maximum-order predictor */
133 if (check_flag(&flag, "ARKStepSetPredictorMethod", 1)) return 1;
134 flag = ARKStepSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */
135 if (check_flag(&flag, "ARKStepSStolerances", 1)) return 1;
136 flag = ARKStepSetConstraints(arkode_mem, constraints); /* Set constraints */
137 if (check_flag(&flag, "ARKStepSetConstraints", 1)) return 1;
138
139 /* Initialize dense matrix data structure and solver */
140 A = SUNDenseMatrix(NEQ, NEQ);
141 if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
142 LS = SUNLinSol_Dense(y, A);
143 if (check_flag((void *)LS, "SUNLinSol_Dense", 0)) return 1;
144
145 /* Linear solver interface */
146 flag = ARKStepSetLinearSolver(arkode_mem, LS, A); /* Attach matrix and linear solver */
147 if (check_flag(&flag, "ARKStepSetLinearSolver", 1)) return 1;
148 flag = ARKStepSetJacFn(arkode_mem, Jac); /* Set the Jacobian routine */
149 if (check_flag(&flag, "ARKStepSetJacFn", 1)) return 1;
150
151 /* Open output stream for results, output comment line */
152 UFID = fopen("solution.txt","w");
153 fprintf(UFID,"# t u v w\n");
154
155 /* output initial condition to disk */
156 fprintf(UFID," %.16"ESYM" %.16"ESYM" %.16"ESYM" %.16"ESYM"\n",
157 T0, NV_Ith_S(y,0), NV_Ith_S(y,1), NV_Ith_S(y,2));
158
159 /* Main time-stepping loop: calls ARKStepEvolve to perform the integration, then
160 prints results. Stops when the final time has been reached */
161 t = T0;
162 tout = T0+dTout;
163 printf(" t u v w\n");
164 printf(" --------------------------------------------------\n");
165 printf(" %10.3"ESYM" %12.5"ESYM" %12.5"ESYM" %12.5"ESYM"\n",
166 t, NV_Ith_S(y,0), NV_Ith_S(y,1), NV_Ith_S(y,2));
167 for (iout=0; iout<Nt; iout++) {
168
169 flag = ARKStepEvolve(arkode_mem, tout, y, &t, ARK_NORMAL); /* call integrator */
170 if (check_flag(&flag, "ARKStepEvolve", 1)) break;
171 printf(" %10.3"ESYM" %12.5"ESYM" %12.5"ESYM" %12.5"ESYM"\n", /* access/print solution */
172 t, NV_Ith_S(y,0), NV_Ith_S(y,1), NV_Ith_S(y,2));
173 fprintf(UFID," %.16"ESYM" %.16"ESYM" %.16"ESYM" %.16"ESYM"\n",
174 t, NV_Ith_S(y,0), NV_Ith_S(y,1), NV_Ith_S(y,2));
175 if (flag >= 0) { /* successful solve: update time */
176 tout += dTout;
177 tout = (tout > Tf) ? Tf : tout;
178 } else { /* unsuccessful solve: break */
179 fprintf(stderr,"Solver failure, stopping integration\n");
180 break;
181 }
182 }
183 printf(" --------------------------------------------------\n");
184 fclose(UFID);
185
186 /* Print some final statistics */
187 flag = ARKStepGetNumSteps(arkode_mem, &nst);
188 check_flag(&flag, "ARKStepGetNumSteps", 1);
189 flag = ARKStepGetNumStepAttempts(arkode_mem, &nst_a);
190 check_flag(&flag, "ARKStepGetNumStepAttempts", 1);
191 flag = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi);
192 check_flag(&flag, "ARKStepGetNumRhsEvals", 1);
193 flag = ARKStepGetNumLinSolvSetups(arkode_mem, &nsetups);
194 check_flag(&flag, "ARKStepGetNumLinSolvSetups", 1);
195 flag = ARKStepGetNumErrTestFails(arkode_mem, &netf);
196 check_flag(&flag, "ARKStepGetNumErrTestFails", 1);
197 flag = ARKStepGetNumNonlinSolvIters(arkode_mem, &nni);
198 check_flag(&flag, "ARKStepGetNumNonlinSolvIters", 1);
199 flag = ARKStepGetNumNonlinSolvConvFails(arkode_mem, &ncfn);
200 check_flag(&flag, "ARKStepGetNumNonlinSolvConvFails", 1);
201 flag = ARKStepGetNumJacEvals(arkode_mem, &nje);
202 check_flag(&flag, "ARKStepGetNumJacEvals", 1);
203 flag = ARKStepGetNumLinRhsEvals(arkode_mem, &nfeLS);
204 check_flag(&flag, "ARKStepGetNumLinRhsEvals", 1);
205 flag = ARKStepGetNumConstrFails(arkode_mem, &nctf);
206 check_flag(&flag, "ARKStepGetNumConstrFails", 1);
207
208 printf("\nFinal Solver Statistics:\n");
209 printf(" Internal solver steps = %li (attempted = %li)\n",
210 nst, nst_a);
211 printf(" Total RHS evals: Fe = %li, Fi = %li\n", nfe, nfi);
212 printf(" Total linear solver setups = %li\n", nsetups);
213 printf(" Total RHS evals for setting up the linear system = %li\n", nfeLS);
214 printf(" Total number of Jacobian evaluations = %li\n", nje);
215 printf(" Total number of Newton iterations = %li\n", nni);
216 printf(" Total number of nonlinear solver convergence failures = %li\n", ncfn);
217 printf(" Total number of error test failures = %li\n", netf);
218 printf(" Total number of constraint test failures = %li\n", nctf);
219
220 /* check the solution error */
221 flag = check_ans(y, t, reltol, abstol);
222
223 /* Clean up and return with successful completion */
224 N_VDestroy(y); /* Free y vector */
225 N_VDestroy(constraints); /* Free constraints vector */
226 ARKStepFree(&arkode_mem); /* Free integrator memory */
227 SUNLinSolFree(LS); /* Free linear solver */
228 SUNMatDestroy(A); /* Free A matrix */
229
230 return flag;
231 }
232
233 /*-------------------------------
234 * Functions called by the solver
235 *-------------------------------*/
236
237 /* f routine to compute the ODE RHS function f(t,y). */
f(realtype t,N_Vector y,N_Vector ydot,void * user_data)238 static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
239 {
240 realtype u = NV_Ith_S(y,0); /* access current solution */
241 realtype v = NV_Ith_S(y,1);
242 realtype w = NV_Ith_S(y,2);
243
244 /* Fill in ODE RHS function */
245 NV_Ith_S(ydot,0) = -0.04*u + 1.e4*v*w;
246 NV_Ith_S(ydot,1) = 0.04*u - 1.e4*v*w - 3.e7*v*v;
247 NV_Ith_S(ydot,2) = 3.e7*v*v;
248
249 return 0; /* Return with success */
250 }
251
252 /* Jacobian routine to compute J(t,y) = df/dy. */
Jac(realtype t,N_Vector y,N_Vector fy,SUNMatrix J,void * user_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)253 static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
254 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
255 {
256 realtype v = NV_Ith_S(y,1); /* access current solution */
257 realtype w = NV_Ith_S(y,2);
258 SUNMatZero(J); /* initialize Jacobian to zero */
259
260 /* Fill in the Jacobian of the ODE RHS function */
261 SM_ELEMENT_D(J,0,0) = -0.04;
262 SM_ELEMENT_D(J,0,1) = 1.e4*w;
263 SM_ELEMENT_D(J,0,2) = 1.e4*v;
264
265 SM_ELEMENT_D(J,1,0) = 0.04;
266 SM_ELEMENT_D(J,1,1) = -1.e4*w - 6.e7*v;
267 SM_ELEMENT_D(J,1,2) = -1.e4*v;
268
269 SM_ELEMENT_D(J,2,1) = 6.e7*v;
270
271 return 0; /* Return with success */
272 }
273
274 /*-------------------------------
275 * Private helper functions
276 *-------------------------------*/
277
278 /* Check function return value...
279 opt == 0 means SUNDIALS function allocates memory so check if
280 returned NULL pointer
281 opt == 1 means SUNDIALS function returns a flag so check if
282 flag >= 0
283 opt == 2 means function allocates memory so check if returned
284 NULL pointer
285 */
check_flag(void * flagvalue,const char * funcname,int opt)286 static int check_flag(void *flagvalue, const char *funcname, int opt)
287 {
288 int *errflag;
289
290 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
291 if (opt == 0 && flagvalue == NULL) {
292 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
293 funcname);
294 return 1; }
295
296 /* Check if flag < 0 */
297 else if (opt == 1) {
298 errflag = (int *) flagvalue;
299 if (*errflag < 0) {
300 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
301 funcname, *errflag);
302 return 1; }}
303
304 /* Check if function returned NULL pointer - no memory allocated */
305 else if (opt == 2 && flagvalue == NULL) {
306 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
307 funcname);
308 return 1; }
309
310 return 0;
311 }
312
313 /* compare the solution at the final time 1e11s to a reference solution computed
314 using a relative tolerance of 1e-8 and absoltue tolerance of 1e-14 */
check_ans(N_Vector y,realtype t,realtype rtol,realtype atol)315 static int check_ans(N_Vector y, realtype t, realtype rtol, realtype atol)
316 {
317 int passfail=0; /* answer pass (0) or fail (1) flag */
318 N_Vector ref; /* reference solution vector */
319 N_Vector ewt; /* error weight vector */
320 realtype err; /* wrms error */
321 realtype ZERO=RCONST(0.0);
322 realtype ONE=RCONST(1.0);
323
324 /* create reference solution and error weight vectors */
325 ref = N_VClone(y);
326 ewt = N_VClone(y);
327
328 /* set the reference solution data */
329 NV_Ith_S(ref,0) = RCONST(2.0833403356917897e-08);
330 NV_Ith_S(ref,1) = RCONST(8.1470714598028223e-14);
331 NV_Ith_S(ref,2) = RCONST(9.9999997916651040e-01);
332
333 /* compute the error weight vector */
334 N_VAbs(ref, ewt);
335 N_VScale(rtol, ewt, ewt);
336 N_VAddConst(ewt, atol, ewt);
337 if (N_VMin(ewt) <= ZERO) {
338 fprintf(stderr, "\nSUNDIALS_ERROR: check_ans failed - ewt <= 0\n\n");
339 return(-1);
340 }
341 N_VInv(ewt, ewt);
342
343 /* compute the solution error */
344 N_VLinearSum(ONE, y, -ONE, ref, ref);
345 err = N_VWrmsNorm(ref, ewt);
346
347 /* is the solution within the tolerances? */
348 passfail = (err < ONE) ? 0 : 1;
349
350 if (passfail) {
351 fprintf(stdout, "\nSUNDIALS_WARNING: check_ans error=%"GSYM"\n\n", err);
352 }
353
354 /* Free vectors */
355 N_VDestroy(ref);
356 N_VDestroy(ewt);
357
358 return(passfail);
359 }
360
361 /*---- end of file ----*/
362