1 /* -----------------------------------------------------------------
2 * Programmer(s): Carol Woodward @ LLNL
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 * Example problem:
15 *
16 * The following is a simple example problem, with the coding
17 * needed for its solution by the accelerated fixed point solver in
18 * KINSOL.
19 * The problem is from chemical kinetics, and consists of solving
20 * the first time step in a Backward Euler solution for the
21 * following three rate equations:
22 * dy1/dt = -.04*y1 + 1.e4*y2*y3
23 * dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e2*(y2)^2
24 * dy3/dt = 3.e2*(y2)^2
25 * on the interval from t = 0.0 to t = 0.1, with initial
26 * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.
27 * Run statistics (optional outputs) are printed at the end.
28 * -----------------------------------------------------------------
29 */
30
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <math.h>
34
35 #include <kinsol/kinsol.h> /* access to KINSOL func., consts. */
36 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
37 #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
38
39 #if defined(SUNDIALS_EXTENDED_PRECISION)
40 #define GSYM "Lg"
41 #define ESYM "Le"
42 #define FSYM "Lf"
43 #else
44 #define GSYM "g"
45 #define ESYM "e"
46 #define FSYM "f"
47 #endif
48
49 /* Problem Constants */
50
51 #define NEQ 3 /* number of equations */
52 #define Y10 RCONST(1.0) /* initial y components */
53 #define Y20 RCONST(0.0)
54 #define Y30 RCONST(0.0)
55 #define TOL RCONST(1.e-10) /* function tolerance */
56 #define DSTEP RCONST(0.1) /* Size of the single time step used */
57
58 #define PRIORS 2
59
60 #define ZERO RCONST(0.0)
61 #define ONE RCONST(1.0)
62
63 /* User-defined vector accessor macro: Ith */
64
65 /* This macro is defined in order to write code which exactly matches
66 the mathematical problem description given above.
67
68 Ith(v,i) references the ith component of the vector v, where i is in
69 the range [1..NEQ] and NEQ is defined above. The Ith macro is defined
70 using the N_VIth macro in nvector.h. N_VIth numbers the components of
71 a vector starting from 0.
72 */
73
74 #define Ith(v,i) NV_Ith_S(v,i-1) /* Ith numbers components 1..NEQ */
75
76
77 /* Private functions */
78
79 static int funcRoberts(N_Vector u, N_Vector f, void *user_data);
80 static void PrintOutput(N_Vector u);
81 static void PrintFinalStats(void *kmem);
82 static int check_flag(void *flagvalue, const char *funcname, int opt);
83 static int check_ans(N_Vector u, realtype rtol, realtype atol);
84
85 /*
86 *--------------------------------------------------------------------
87 * MAIN PROGRAM
88 *--------------------------------------------------------------------
89 */
90
main()91 int main()
92 {
93 realtype fnormtol, fnorm;
94 N_Vector y, scale;
95 int flag;
96 void *kmem;
97
98 fnorm = 0.0;
99 y = scale = NULL;
100 kmem = NULL;
101
102 /* -------------------------
103 * Print problem description
104 * ------------------------- */
105
106 printf("Example problem from chemical kinetics solving\n");
107 printf("the first time step in a Backward Euler solution for the\n");
108 printf("following three rate equations:\n");
109 printf(" dy1/dt = -.04*y1 + 1.e4*y2*y3\n");
110 printf(" dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e2*(y2)^2\n");
111 printf(" dy3/dt = 3.e2*(y2)^2\n");
112 printf("on the interval from t = 0.0 to t = 0.1, with initial\n");
113 printf("conditions: y1 = 1.0, y2 = y3 = 0.\n");
114 printf("Solution method: Anderson accelerated fixed point iteration.\n");
115
116 /* --------------------------------------
117 * Create vectors for solution and scales
118 * -------------------------------------- */
119
120 y = N_VNew_Serial(NEQ);
121 if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);
122
123 scale = N_VNew_Serial(NEQ);
124 if (check_flag((void *)scale, "N_VNew_Serial", 0)) return(1);
125
126 /* -----------------------------------------
127 * Initialize and allocate memory for KINSOL
128 * ----------------------------------------- */
129
130 kmem = KINCreate();
131 if (check_flag((void *)kmem, "KINCreate", 0)) return(1);
132
133 /* y is used as a template */
134
135 /* Set number of prior residuals used in Anderson acceleration */
136 flag = KINSetMAA(kmem, PRIORS);
137
138 flag = KINInit(kmem, funcRoberts, y);
139 if (check_flag(&flag, "KINInit", 1)) return(1);
140
141 /* -------------------
142 * Set optional inputs
143 * ------------------- */
144
145 /* Specify stopping tolerance based on residual */
146
147 fnormtol = TOL;
148 flag = KINSetFuncNormTol(kmem, fnormtol);
149 if (check_flag(&flag, "KINSetFuncNormTol", 1)) return(1);
150
151 /* -------------
152 * Initial guess
153 * ------------- */
154
155 N_VConst(ZERO, y);
156 Ith(y,1) = ONE;
157
158 /* ----------------------------
159 * Call KINSol to solve problem
160 * ---------------------------- */
161
162 /* No scaling used */
163 N_VConst(ONE,scale);
164
165 /* Call main solver */
166 flag = KINSol(kmem, /* KINSol memory block */
167 y, /* initial guess on input; solution vector */
168 KIN_FP, /* global strategy choice */
169 scale, /* scaling vector, for the variable cc */
170 scale); /* scaling vector for function values fval */
171 if (check_flag(&flag, "KINSol", 1)) return(1);
172
173
174 /* ------------------------------------
175 * Print solution and solver statistics
176 * ------------------------------------ */
177
178 /* Get scaled norm of the system function */
179
180 flag = KINGetFuncNorm(kmem, &fnorm);
181 if (check_flag(&flag, "KINGetfuncNorm", 1)) return(1);
182
183 #if defined(SUNDIALS_EXTENDED_PRECISION)
184 printf("\nComputed solution (||F|| = %Lg):\n\n",fnorm);
185 #else
186 printf("\nComputed solution (||F|| = %g):\n\n",fnorm);
187 #endif
188 PrintOutput(y);
189
190 PrintFinalStats(kmem);
191
192 /* check the solution error */
193 flag = check_ans(y, RCONST(1e-4), RCONST(1e-6));
194
195 /* -----------
196 * Free memory
197 * ----------- */
198
199 N_VDestroy(y);
200 N_VDestroy(scale);
201 KINFree(&kmem);
202
203 return(flag);
204 }
205
206 /*
207 *--------------------------------------------------------------------
208 * PRIVATE FUNCTIONS
209 *--------------------------------------------------------------------
210 */
211
212 /*
213 * System function
214 */
215
funcRoberts(N_Vector y,N_Vector g,void * user_data)216 static int funcRoberts(N_Vector y, N_Vector g, void *user_data)
217 {
218 realtype y1, y2, y3;
219 realtype yd1, yd3;
220
221 y1 = Ith(y,1);
222 y2 = Ith(y,2);
223 y3 = Ith(y,3);
224
225 yd1 = DSTEP * ( RCONST(-0.04)*y1 + RCONST(1.0e4)*y2*y3 );
226 yd3 = DSTEP * RCONST(3.0e2)*y2*y2;
227
228 Ith(g,1) = yd1 + Y10;
229 Ith(g,2) = -yd1 - yd3 + Y20;
230 Ith(g,3) = yd3 + Y30;
231
232 return(0);
233 }
234
235 /*
236 * Print solution at selected points
237 */
238
PrintOutput(N_Vector y)239 static void PrintOutput(N_Vector y)
240 {
241 realtype y1, y2, y3;
242
243 y1 = Ith(y,1);
244 y2 = Ith(y,2);
245 y3 = Ith(y,3);
246
247 #if defined(SUNDIALS_EXTENDED_PRECISION)
248 printf("y =%14.6Le %14.6Le %14.6Le\n", y1, y2, y3);
249 #elif defined(SUNDIALS_DOUBLE_PRECISION)
250 printf("y =%14.6e %14.6e %14.6e\n", y1, y2, y3);
251 #else
252 printf("y =%14.6e %14.6e %14.6e\n", y1, y2, y3);
253 #endif
254
255 return;
256 }
257
258 /*
259 * Print final statistics
260 */
261
PrintFinalStats(void * kmem)262 static void PrintFinalStats(void *kmem)
263 {
264 long int nni, nfe;
265 int flag;
266
267 /* Main solver statistics */
268
269 flag = KINGetNumNonlinSolvIters(kmem, &nni);
270 check_flag(&flag, "KINGetNumNonlinSolvIters", 1);
271 flag = KINGetNumFuncEvals(kmem, &nfe);
272 check_flag(&flag, "KINGetNumFuncEvals", 1);
273
274 printf("\nFinal Statistics.. \n\n");
275 printf("nni = %6ld nfe = %6ld \n", nni, nfe);
276 }
277
278 /*
279 * Check function return value...
280 * opt == 0 means SUNDIALS function allocates memory so check if
281 * returned NULL pointer
282 * opt == 1 means SUNDIALS function returns a flag so check if
283 * flag >= 0
284 * opt == 2 means function allocates memory so check if returned
285 * NULL pointer
286 */
287
check_flag(void * flagvalue,const char * funcname,int opt)288 static int check_flag(void *flagvalue, const char *funcname, int opt)
289 {
290 int *errflag;
291
292 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
293 if (opt == 0 && flagvalue == NULL) {
294 fprintf(stderr,
295 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
296 funcname);
297 return(1);
298 }
299
300 /* Check if flag < 0 */
301 else if (opt == 1) {
302 errflag = (int *) flagvalue;
303 if (*errflag < 0) {
304 fprintf(stderr,
305 "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
306 funcname, *errflag);
307 return(1);
308 }
309 }
310
311 /* Check if function returned NULL pointer - no memory allocated */
312 else if (opt == 2 && flagvalue == NULL) {
313 fprintf(stderr,
314 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
315 funcname);
316 return(1);
317 }
318
319 return(0);
320 }
321
322 /* compare the solution to a reference solution computed with a
323 tolerance of 1e-14 */
check_ans(N_Vector u,realtype rtol,realtype atol)324 static int check_ans(N_Vector u, realtype rtol, realtype atol)
325 {
326 int passfail=0; /* answer pass (0) or fail (1) flag */
327 N_Vector ref; /* reference solution vector */
328 N_Vector ewt; /* error weight vector */
329 realtype err; /* wrms error */
330
331 /* create reference solution and error weight vectors */
332 ref = N_VClone(u);
333 ewt = N_VClone(u);
334
335 /* set the reference solution data */
336 NV_Ith_S(ref,0) = RCONST(9.9678538655358029e-01);
337 NV_Ith_S(ref,1) = RCONST(2.9530060962800345e-03);
338 NV_Ith_S(ref,2) = RCONST(2.6160735013975683e-04);
339
340 /* compute the error weight vector */
341 N_VAbs(ref, ewt);
342 N_VScale(rtol, ewt, ewt);
343 N_VAddConst(ewt, atol, ewt);
344 if (N_VMin(ewt) <= ZERO) {
345 fprintf(stderr, "\nSUNDIALS_ERROR: check_ans failed - ewt <= 0\n\n");
346 return(-1);
347 }
348 N_VInv(ewt, ewt);
349
350 /* compute the solution error */
351 N_VLinearSum(ONE, u, -ONE, ref, ref);
352 err = N_VWrmsNorm(ref, ewt);
353
354 /* is the solution within the tolerances? */
355 passfail = (err < ONE) ? 0 : 1;
356
357 if (passfail) {
358 fprintf(stdout, "\nSUNDIALS_WARNING: check_ans error=%"GSYM"\n\n", err);
359 }
360
361 /* Free vectors */
362 N_VDestroy(ref);
363 N_VDestroy(ewt);
364
365 return(passfail);
366 }
367