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 * Demonstration program for CVODE - direct linear solvers.
16 * Two separate problems are solved using both the CV_ADAMS and CV_BDF
17 * linear multistep methods in combination with the
18 * SUNNONLINSOL_FIXEDPOINT and SUNNONLINSOL_NEWTON nonlinear solver
19 * modules:
20 *
21 * Problem 1: Van der Pol oscillator
22 * xdotdot - 3*(1 - x^2)*xdot + x = 0, x(0) = 2, xdot(0) = 0.
23 * This second-order ODE is converted to a first-order system by
24 * defining y0 = x and y1 = xdot.
25 * The NEWTON iteration cases use the following types of Jacobian
26 * approximation: (1) dense, user-supplied, (2) dense, difference
27 * quotient approximation, (3) diagonal approximation.
28 *
29 * Problem 2: ydot = A * y, where A is a banded lower triangular
30 * matrix derived from 2-D advection PDE.
31 * The NEWTON iteration cases use the following types of Jacobian
32 * approximation: (1) band, user-supplied, (2) band, difference
33 * quotient approximation, (3) diagonal approximation.
34 *
35 * For each problem, in the series of eight runs, CVodeInit is
36 * called only once, for the first run, whereas CVodeReInit is
37 * called for each of the remaining seven runs.
38 *
39 * Notes: This program demonstrates the usage of the sequential
40 * macros NV_Ith_S, SM_ELEMENT_D, SM_COLUMN_B, and
41 * SM_COLUMN_ELEMENT_B. The NV_Ith_S macro is used to reference the
42 * components of an N_Vector. It works for any size N=NEQ, but
43 * due to efficiency concerns it should only by used when the
44 * problem size is small. The Problem 1 right hand side and
45 * Jacobian functions f1 and Jac1 both use NV_Ith_S. The
46 * N_VGetArrayPointer function gives the user access to the
47 * memory used for the component storage of an N_Vector. In the
48 * sequential case, the user may assume that this is one contiguous
49 * array of reals. The N_VGetArrayPointer function
50 * gives a more efficient means (than the NV_Ith_S macro) to
51 * access the components of an N_Vector and should be used when the
52 * problem size is large. The Problem 2 right hand side function f2
53 * uses the N_VGetArrayPointer function. The SM_ELEMENT_D macro
54 * used in Jac1 gives access to an element of a dense SUNMatrix. It
55 * should be used only when the problem size is small (the
56 * size of a Dense SUNMatrix is NEQ x NEQ) due to efficiency concerns. For
57 * larger problem sizes, the macro SM_COLUMN_D can be used in order
58 * to work directly with a column of a Dense SUNMatrix. The SM_COLUMN_B and
59 * SM_COLUMN_ELEMENT_B allow efficient columnwise access to the elements
60 * of a Banded SUNMatix. These macros are used in the Jac2 function.
61 * -----------------------------------------------------------------*/
62
63 #include <stdio.h>
64 #include <stdlib.h>
65 #include <math.h>
66
67 #include <cvodes/cvodes.h> /* prototypes for CVODE fcts., consts. */
68 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
69 #include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
70 #include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
71 #include <sunmatrix/sunmatrix_band.h> /* access to band SUNMatrix */
72 #include <sunlinsol/sunlinsol_band.h> /* access to band SUNLinearSolver */
73 #include <cvodes/cvodes_diag.h> /* access to CVDIAG linear solver */
74 #include "sunnonlinsol/sunnonlinsol_newton.h" /* access to the newton SUNNonlinearSolver */
75 #include "sunnonlinsol/sunnonlinsol_fixedpoint.h" /* access to the fixed point SUNNonlinearSolver */
76 #include <sundials/sundials_types.h> /* definition of realtype */
77
78 /* helpful macros */
79
80 #ifndef SQR
81 #define SQR(A) ((A)*(A))
82 #endif
83
84 /* Shared Problem Constants */
85
86 #define ATOL RCONST(1.0e-6)
87 #define RTOL RCONST(0.0)
88
89 #define ZERO RCONST(0.0)
90 #define ONE RCONST(1.0)
91 #define TWO RCONST(2.0)
92 #define THIRTY RCONST(30.0)
93
94 /* Problem #1 Constants */
95
96 #define P1_NEQ 2
97 #define P1_ETA RCONST(3.0)
98 #define P1_NOUT 4
99 #define P1_T0 RCONST(0.0)
100 #define P1_T1 RCONST(1.39283880203)
101 #define P1_DTOUT RCONST(2.214773875)
102 #define P1_TOL_FACTOR RCONST(1.0e4)
103
104 /* Problem #2 Constants */
105
106 #define P2_MESHX 5
107 #define P2_MESHY 5
108 #define P2_NEQ P2_MESHX*P2_MESHY
109 #define P2_ALPH1 RCONST(1.0)
110 #define P2_ALPH2 RCONST(1.0)
111 #define P2_NOUT 5
112 #define P2_ML 5
113 #define P2_MU 0
114 #define P2_T0 RCONST(0.0)
115 #define P2_T1 RCONST(0.01)
116 #define P2_TOUT_MULT RCONST(10.0)
117 #define P2_TOL_FACTOR RCONST(1.0e3)
118
119 /* Linear Solver Options */
120
121 enum {FUNC, DENSE_USER, DENSE_DQ, DIAG, BAND_USER, BAND_DQ};
122
123 /* Private Helper Functions */
124
125 static int Problem1(void);
126 static void PrintIntro1(void);
127 static void PrintHeader1(void);
128 static void PrintOutput1(realtype t, realtype y0, realtype y1, int qu, realtype hu);
129 static int Problem2(void);
130 static void PrintIntro2(void);
131 static void PrintHeader2(void);
132 static void PrintOutput2(realtype t, realtype erm, int qu, realtype hu);
133 static realtype MaxError(N_Vector y, realtype t);
134 static int PrepareNextRun(void *cvode_mem, int lmm, int miter, N_Vector y,
135 SUNMatrix* A, sunindextype mu, sunindextype ml,
136 SUNLinearSolver* LS, SUNNonlinearSolver* NLS);
137 static void PrintErrOutput(realtype tol_factor);
138 static void PrintFinalStats(void *cvode_mem, int miter, realtype ero);
139 static void PrintErrInfo(int nerr);
140
141 /* Functions Called by the Solver */
142
143 static int f1(realtype t, N_Vector y, N_Vector ydot, void *user_data);
144 static int Jac1(realtype tn, N_Vector y, N_Vector fy, SUNMatrix J,
145 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
146
147 static int f2(realtype t, N_Vector y, N_Vector ydot, void *user_data);
148 static int Jac2(realtype tn, N_Vector y, N_Vector fy, SUNMatrix J, void *user_data,
149 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
150
151 /* Private function to check function return values */
152
153 static int check_retval(void *returnvalue, const char *funcname, int opt);
154
155 /* Implementation */
156
main(void)157 int main(void)
158 {
159 int nerr;
160
161 nerr = Problem1();
162 nerr += Problem2();
163 PrintErrInfo(nerr);
164
165 return(0);
166 }
167
Problem1(void)168 static int Problem1(void)
169 {
170 realtype reltol=RTOL, abstol=ATOL, t, tout, ero, er;
171 int miter, retval, temp_retval, iout, nerr=0;
172 N_Vector y;
173 SUNMatrix A;
174 SUNLinearSolver LS;
175 SUNNonlinearSolver NLS;
176 void *cvode_mem;
177 booleantype firstrun;
178 int qu;
179 realtype hu;
180
181 y = NULL;
182 A = NULL;
183 LS = NULL;
184 NLS = NULL;
185 cvode_mem = NULL;
186
187 y = N_VNew_Serial(P1_NEQ);
188 if(check_retval((void *)y, "N_VNew_Serial", 0)) return(1);
189 PrintIntro1();
190
191 cvode_mem = CVodeCreate(CV_ADAMS);
192 if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
193
194 for (miter=FUNC; miter <= DIAG; miter++) {
195 ero = ZERO;
196 NV_Ith_S(y,0) = TWO;
197 NV_Ith_S(y,1) = ZERO;
198
199 firstrun = (miter==FUNC);
200 if (firstrun) {
201
202 /* initialize CVode */
203 retval = CVodeInit(cvode_mem, f1, P1_T0, y);
204 if(check_retval(&retval, "CVodeInit", 1)) return(1);
205
206 /* set scalar tolerances */
207 retval = CVodeSStolerances(cvode_mem, reltol, abstol);
208 if(check_retval(&retval, "CVodeSStolerances", 1)) return(1);
209
210 } else {
211
212 /* reinitialize CVode */
213 retval = CVodeReInit(cvode_mem, P1_T0, y);
214 if(check_retval(&retval, "CVodeReInit", 1)) return(1);
215
216 }
217
218 retval = PrepareNextRun(cvode_mem, CV_ADAMS, miter, y, &A, 0, 0, &LS, &NLS);
219 if(check_retval(&retval, "PrepareNextRun", 1)) return(1);
220
221 PrintHeader1();
222
223 for(iout=1, tout=P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT) {
224 retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
225 check_retval(&retval, "CVode", 1);
226 temp_retval = CVodeGetLastOrder(cvode_mem, &qu);
227 if(check_retval(&temp_retval, "CVodeGetLastOrder", 1)) ++nerr;
228 temp_retval = CVodeGetLastStep(cvode_mem, &hu);
229 if(check_retval(&temp_retval, "CVodeGetLastStep", 1)) ++nerr;
230 PrintOutput1(t, NV_Ith_S(y,0), NV_Ith_S(y,1), qu, hu);
231 if (retval != CV_SUCCESS) {
232 nerr++;
233 break;
234 }
235 if (iout%2 == 0) {
236 er = fabs(NV_Ith_S(y,0)) / abstol;
237 if (er > ero) ero = er;
238 if (er > P1_TOL_FACTOR) {
239 nerr++;
240 PrintErrOutput(P1_TOL_FACTOR);
241 }
242 }
243 }
244
245 PrintFinalStats(cvode_mem, miter, ero);
246 }
247
248 CVodeFree(&cvode_mem);
249 SUNNonlinSolFree(NLS);
250 NLS = NULL;
251 LS = NULL;
252 A = NULL;
253
254 cvode_mem = CVodeCreate(CV_BDF);
255 if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
256
257 for (miter=FUNC; miter <= DIAG; miter++) {
258 ero = ZERO;
259 NV_Ith_S(y,0) = TWO;
260 NV_Ith_S(y,1) = ZERO;
261
262 firstrun = (miter==FUNC);
263 if (firstrun) {
264
265 /* initialize CVode */
266 retval = CVodeInit(cvode_mem, f1, P1_T0, y);
267 if(check_retval(&retval, "CVodeInit", 1)) return(1);
268
269 /* set scalar tolerances */
270 retval = CVodeSStolerances(cvode_mem, reltol, abstol);
271 if(check_retval(&retval, "CVodeSStolerances", 1)) return(1);
272
273 } else {
274
275 /* reinitialize CVode */
276 retval = CVodeReInit(cvode_mem, P1_T0, y);
277 if(check_retval(&retval, "CVodeReInit", 1)) return(1);
278
279 }
280
281 retval = PrepareNextRun(cvode_mem, CV_BDF, miter, y, &A, 0, 0, &LS, &NLS);
282 if(check_retval(&retval, "PrepareNextRun", 1)) return(1);
283
284 PrintHeader1();
285
286 for(iout=1, tout=P1_T1; iout <= P1_NOUT; iout++, tout += P1_DTOUT) {
287 retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
288 check_retval(&retval, "CVode", 1);
289 temp_retval = CVodeGetLastOrder(cvode_mem, &qu);
290 if(check_retval(&temp_retval, "CVodeGetLastOrder", 1)) ++nerr;
291 temp_retval = CVodeGetLastStep(cvode_mem, &hu);
292 if(check_retval(&temp_retval, "CVodeGetLastStep", 1)) ++nerr;
293 PrintOutput1(t, NV_Ith_S(y,0), NV_Ith_S(y,1), qu, hu);
294 if (retval != CV_SUCCESS) {
295 nerr++;
296 break;
297 }
298 if (iout%2 == 0) {
299 er = fabs(NV_Ith_S(y,0)) / abstol;
300 if (er > ero) ero = er;
301 if (er > P1_TOL_FACTOR) {
302 nerr++;
303 PrintErrOutput(P1_TOL_FACTOR);
304 }
305 }
306 }
307
308 PrintFinalStats(cvode_mem, miter, ero);
309 }
310
311 CVodeFree(&cvode_mem);
312 SUNNonlinSolFree(NLS);
313 N_VDestroy(y);
314
315 return(nerr);
316 }
317
PrintIntro1(void)318 static void PrintIntro1(void)
319 {
320 printf("Demonstration program for CVODE package - direct linear solvers\n");
321 printf("\n\n");
322 printf("Problem 1: Van der Pol oscillator\n");
323 printf(" xdotdot - 3*(1 - x^2)*xdot + x = 0, x(0) = 2, xdot(0) = 0\n");
324 #if defined(SUNDIALS_EXTENDED_PRECISION)
325 printf(" neq = %d, reltol = %.2Lg, abstol = %.2Lg",
326 P1_NEQ, RTOL, ATOL);
327 #elif defined(SUNDIALS_DOUBLE_PRECISION)
328 printf(" neq = %d, reltol = %.2g, abstol = %.2g",
329 P1_NEQ, RTOL, ATOL);
330 #else
331 printf(" neq = %d, reltol = %.2g, abstol = %.2g",
332 P1_NEQ, RTOL, ATOL);
333 #endif
334 }
335
PrintHeader1(void)336 static void PrintHeader1(void)
337 {
338 printf("\n t x xdot qu hu \n");
339
340 return;
341 }
342
PrintOutput1(realtype t,realtype y0,realtype y1,int qu,realtype hu)343 static void PrintOutput1(realtype t, realtype y0, realtype y1, int qu, realtype hu)
344 {
345 #if defined(SUNDIALS_EXTENDED_PRECISION)
346 printf("%10.5Lf %12.5Le %12.5Le %2d %6.4Le\n", t, y0, y1, qu, hu);
347 #elif defined(SUNDIALS_DOUBLE_PRECISION)
348 printf("%10.5f %12.5e %12.5e %2d %6.4e\n", t, y0, y1, qu, hu);
349 #else
350 printf("%10.5f %12.5e %12.5e %2d %6.4e\n", t, y0, y1, qu, hu);
351 #endif
352
353 return;
354 }
355
f1(realtype t,N_Vector y,N_Vector ydot,void * user_data)356 static int f1(realtype t, N_Vector y, N_Vector ydot, void *user_data)
357 {
358 realtype y0, y1;
359
360 y0 = NV_Ith_S(y,0);
361 y1 = NV_Ith_S(y,1);
362
363 NV_Ith_S(ydot,0) = y1;
364 NV_Ith_S(ydot,1) = (ONE - SQR(y0))* P1_ETA * y1 - y0;
365
366 return(0);
367 }
368
Jac1(realtype tn,N_Vector y,N_Vector fy,SUNMatrix J,void * user_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)369 static int Jac1(realtype tn, N_Vector y, N_Vector fy, SUNMatrix J,
370 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
371 {
372 realtype y0, y1;
373
374 y0 = NV_Ith_S(y,0);
375 y1 = NV_Ith_S(y,1);
376
377 SM_ELEMENT_D(J,0,1) = ONE;
378 SM_ELEMENT_D(J,1,0) = -TWO * P1_ETA * y0 * y1 - ONE;
379 SM_ELEMENT_D(J,1,1) = P1_ETA * (ONE - SQR(y0));
380
381 return(0);
382 }
383
Problem2(void)384 static int Problem2(void)
385 {
386 realtype reltol=RTOL, abstol=ATOL, t, tout, er, erm, ero;
387 int miter, retval, temp_retval, nerr=0;
388 N_Vector y;
389 SUNMatrix A;
390 SUNLinearSolver LS;
391 SUNNonlinearSolver NLS;
392 void *cvode_mem;
393 booleantype firstrun;
394 int qu, iout;
395 realtype hu;
396
397 y = NULL;
398 A = NULL;
399 LS = NULL;
400 NLS = NULL;
401 cvode_mem = NULL;
402
403 y = N_VNew_Serial(P2_NEQ);
404 if(check_retval((void *)y, "N_VNew_Serial", 0)) return(1);
405
406 PrintIntro2();
407
408 cvode_mem = CVodeCreate(CV_ADAMS);
409 if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
410
411 for (miter=FUNC; miter <= BAND_DQ; miter++) {
412 if ((miter==DENSE_USER) || (miter==DENSE_DQ)) continue;
413 ero = ZERO;
414 N_VConst(ZERO, y);
415 NV_Ith_S(y,0) = ONE;
416
417 firstrun = (miter==FUNC);
418 if (firstrun) {
419
420 /* initialize CVode */
421 retval = CVodeInit(cvode_mem, f2, P2_T0, y);
422 if(check_retval(&retval, "CVodeInit", 1)) return(1);
423
424 /* set scalar tolerances */
425 retval = CVodeSStolerances(cvode_mem, reltol, abstol);
426 if(check_retval(&retval, "CVodeSStolerances", 1)) return(1);
427
428 } else {
429
430 /* reinitialize CVode */
431 retval = CVodeReInit(cvode_mem, P2_T0, y);
432 if(check_retval(&retval, "CVodeReInit", 1)) return(1);
433
434 }
435
436 retval = PrepareNextRun(cvode_mem, CV_ADAMS, miter, y, &A, P2_MU, P2_ML, &LS, &NLS);
437 if(check_retval(&retval, "PrepareNextRun", 1)) return(1);
438
439 PrintHeader2();
440
441 for(iout=1, tout=P2_T1; iout <= P2_NOUT; iout++, tout*=P2_TOUT_MULT) {
442 retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
443 check_retval(&retval, "CVode", 1);
444 erm = MaxError(y, t);
445 temp_retval = CVodeGetLastOrder(cvode_mem, &qu);
446 if(check_retval(&temp_retval, "CVodeGetLastOrder", 1)) ++nerr;
447 temp_retval = CVodeGetLastStep(cvode_mem, &hu);
448 if(check_retval(&temp_retval, "CVodeGetLastStep", 1)) ++nerr;
449 PrintOutput2(t, erm, qu, hu);
450 if (retval != CV_SUCCESS) {
451 nerr++;
452 break;
453 }
454 er = erm / abstol;
455 if (er > ero) ero = er;
456 if (er > P2_TOL_FACTOR) {
457 nerr++;
458 PrintErrOutput(P2_TOL_FACTOR);
459 }
460 }
461
462 PrintFinalStats(cvode_mem, miter, ero);
463 }
464
465 CVodeFree(&cvode_mem);
466 SUNNonlinSolFree(NLS);
467 SUNLinSolFree(LS);
468 SUNMatDestroy(A);
469 NLS = NULL;
470 LS = NULL;
471 A = NULL;
472
473 cvode_mem = CVodeCreate(CV_BDF);
474 if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
475
476 for (miter=FUNC; miter <= BAND_DQ; miter++) {
477 if ((miter==DENSE_USER) || (miter==DENSE_DQ)) continue;
478 ero = ZERO;
479 N_VConst(ZERO, y);
480 NV_Ith_S(y,0) = ONE;
481
482 firstrun = (miter==FUNC);
483 if (firstrun) {
484
485 /* initialize CVode */
486 retval = CVodeInit(cvode_mem, f2, P2_T0, y);
487 if(check_retval(&retval, "CVodeInit", 1)) return(1);
488
489 /* set scalar tolerances */
490 retval = CVodeSStolerances(cvode_mem, reltol, abstol);
491 if(check_retval(&retval, "CVodeSStolerances", 1)) return(1);
492
493 } else {
494
495 /* reinitialize CVode */
496 retval = CVodeReInit(cvode_mem, P2_T0, y);
497 if(check_retval(&retval, "CVodeReInit", 1)) return(1);
498
499 }
500
501 retval = PrepareNextRun(cvode_mem, CV_BDF, miter, y, &A, P2_MU, P2_ML, &LS, &NLS);
502 if(check_retval(&retval, "PrepareNextRun", 1)) return(1);
503
504 PrintHeader2();
505
506 for(iout=1, tout=P2_T1; iout <= P2_NOUT; iout++, tout*=P2_TOUT_MULT) {
507 retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
508 check_retval(&retval, "CVode", 1);
509 erm = MaxError(y, t);
510 temp_retval = CVodeGetLastOrder(cvode_mem, &qu);
511 if(check_retval(&temp_retval, "CVodeGetLastOrder", 1)) ++nerr;
512 temp_retval = CVodeGetLastStep(cvode_mem, &hu);
513 if(check_retval(&temp_retval, "CVodeGetLastStep", 1)) ++nerr;
514 PrintOutput2(t, erm, qu, hu);
515 if (retval != CV_SUCCESS) {
516 nerr++;
517 break;
518 }
519 er = erm / abstol;
520 if (er > ero) ero = er;
521 if (er > P2_TOL_FACTOR) {
522 nerr++;
523 PrintErrOutput(P2_TOL_FACTOR);
524 }
525 }
526
527 PrintFinalStats(cvode_mem, miter, ero);
528 }
529
530 CVodeFree(&cvode_mem);
531 SUNNonlinSolFree(NLS);
532 SUNLinSolFree(LS);
533 SUNMatDestroy(A);
534 N_VDestroy(y);
535
536 return(nerr);
537 }
538
PrintIntro2(void)539 static void PrintIntro2(void)
540 {
541 printf("\n\n-------------------------------------------------------------");
542 printf("\n-------------------------------------------------------------");
543 printf("\n\nProblem 2: ydot = A * y, where A is a banded lower\n");
544 printf("triangular matrix derived from 2-D advection PDE\n\n");
545 printf(" neq = %d, ml = %d, mu = %d\n", P2_NEQ, P2_ML, P2_MU);
546 #if defined(SUNDIALS_EXTENDED_PRECISION)
547 printf(" itol = %s, reltol = %.2Lg, abstol = %.2Lg", "CV_SS", RTOL, ATOL);
548 #elif defined(SUNDIALS_DOUBLE_PRECISION)
549 printf(" itol = %s, reltol = %.2g, abstol = %.2g", "CV_SS", RTOL, ATOL);
550 #else
551 printf(" itol = %s, reltol = %.2g, abstol = %.2g", "CV_SS", RTOL, ATOL);
552 #endif
553 printf("\n t max.err qu hu \n");
554 }
555
PrintHeader2(void)556 static void PrintHeader2(void)
557 {
558 printf("\n t max.err qu hu \n");
559
560 return;
561 }
562
PrintOutput2(realtype t,realtype erm,int qu,realtype hu)563 static void PrintOutput2(realtype t, realtype erm, int qu, realtype hu)
564 {
565 #if defined(SUNDIALS_EXTENDED_PRECISION)
566 printf("%10.3Lf %12.4Le %2d %12.4Le\n", t, erm, qu, hu);
567 #elif defined(SUNDIALS_DOUBLE_PRECISION)
568 printf("%10.3f %12.4e %2d %12.4e\n", t, erm, qu, hu);
569 #else
570 printf("%10.3f %12.4e %2d %12.4e\n", t, erm, qu, hu);
571 #endif
572
573 return;
574 }
575
f2(realtype t,N_Vector y,N_Vector ydot,void * user_data)576 static int f2(realtype t, N_Vector y, N_Vector ydot, void *user_data)
577 {
578 sunindextype i, j, k;
579 realtype d, *ydata, *dydata;
580
581 ydata = N_VGetArrayPointer(y);
582 dydata = N_VGetArrayPointer(ydot);
583
584 /*
585 Excluding boundaries,
586
587 ydot = f = -2 y + alpha1 * y + alpha2 * y
588 i,j i,j i,j i-1,j i,j-1
589 */
590
591 for (j=0; j < P2_MESHY; j++) {
592 for (i=0; i < P2_MESHX; i++) {
593 k = i + j * P2_MESHX;
594 d = -TWO*ydata[k];
595 if (i != 0) d += P2_ALPH1 * ydata[k-1];
596 if (j != 0) d += P2_ALPH2 * ydata[k-P2_MESHX];
597 dydata[k] = d;
598 }
599 }
600
601 return(0);
602 }
603
Jac2(realtype tn,N_Vector y,N_Vector fy,SUNMatrix J,void * user_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)604 static int Jac2(realtype tn, N_Vector y, N_Vector fy, SUNMatrix J,
605 void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
606 {
607 int i, j, k;
608 realtype *kthCol;
609
610 /*
611 The components of f(t,y) which depend on y are
612 i,j
613 f , f , and f :
614 i,j i+1,j i,j+1
615
616 f = -2 y + alpha1 * y + alpha2 * y
617 i,j i,j i-1,j i,j-1
618
619 f = -2 y + alpha1 * y + alpha2 * y
620 i+1,j i+1,j i,j i+1,j-1
621
622 f = -2 y + alpha1 * y + alpha2 * y
623 i,j+1 i,j+1 i-1,j+1 i,j
624 */
625
626 for (j=0; j < P2_MESHY; j++) {
627 for (i=0; i < P2_MESHX; i++) {
628 k = i + j * P2_MESHX;
629 kthCol = SM_COLUMN_B(J,k);
630 SM_COLUMN_ELEMENT_B(kthCol,k,k) = -TWO;
631 if (i != P2_MESHX-1) SM_COLUMN_ELEMENT_B(kthCol,k+1,k) = P2_ALPH1;
632 if (j != P2_MESHY-1) SM_COLUMN_ELEMENT_B(kthCol,k+P2_MESHX,k) = P2_ALPH2;
633 }
634 }
635
636 return(0);
637 }
638
MaxError(N_Vector y,realtype t)639 static realtype MaxError(N_Vector y, realtype t)
640 {
641 sunindextype i, j, k;
642 realtype *ydata, er, ex=ZERO, yt, maxError=ZERO, ifact_inv, jfact_inv=ONE;
643
644 if (t == ZERO) return(ZERO);
645
646 ydata = N_VGetArrayPointer(y);
647 if (t <= THIRTY) ex = exp(-TWO*t);
648
649 for (j = 0; j < P2_MESHY; j++) {
650 ifact_inv = ONE;
651 for (i = 0; i < P2_MESHX; i++) {
652 k = i + j * P2_MESHX;
653 yt = pow(t, i+j) * ex * ifact_inv * jfact_inv;
654 er = fabs(ydata[k] - yt);
655 if (er > maxError) maxError = er;
656 ifact_inv /= (i+1);
657 }
658 jfact_inv /= (j+1);
659 }
660 return(maxError);
661 }
662
PrepareNextRun(void * cvode_mem,int lmm,int miter,N_Vector y,SUNMatrix * A,sunindextype mu,sunindextype ml,SUNLinearSolver * LS,SUNNonlinearSolver * NLS)663 static int PrepareNextRun(void *cvode_mem, int lmm, int miter, N_Vector y,
664 SUNMatrix* A, sunindextype mu, sunindextype ml,
665 SUNLinearSolver* LS, SUNNonlinearSolver* NLS)
666 {
667 int retval = CV_SUCCESS;
668
669 if (*NLS)
670 SUNNonlinSolFree(*NLS);
671 if (*LS)
672 SUNLinSolFree(*LS);
673 if (*A)
674 SUNMatDestroy(*A);
675
676 printf("\n\n-------------------------------------------------------------");
677
678 printf("\n\nLinear Multistep Method : ");
679 if (lmm == CV_ADAMS) {
680 printf("ADAMS\n");
681 } else {
682 printf("BDF\n");
683 }
684
685 printf("Iteration : ");
686 if (miter == FUNC) {
687 printf("FIXEDPOINT\n");
688
689 /* create fixed point nonlinear solver object */
690 *NLS = SUNNonlinSol_FixedPoint(y, 0);
691 if(check_retval((void *)*NLS, "SUNNonlinSol_FixedPoint", 0)) return(1);
692
693 /* attach nonlinear solver object to CVode */
694 retval = CVodeSetNonlinearSolver(cvode_mem, *NLS);
695 if(check_retval(&retval, "CVodeSetNonlinearSolver", 1)) return(1);
696
697 } else {
698 printf("NEWTON\n");
699
700 /* create Newton nonlinear solver object */
701 *NLS = SUNNonlinSol_Newton(y);
702 if(check_retval((void *)NLS, "SUNNonlinSol_Newton", 0)) return(1);
703
704 /* attach nonlinear solver object to CVode */
705 retval = CVodeSetNonlinearSolver(cvode_mem, *NLS);
706 if(check_retval(&retval, "CVodeSetNonlinearSolver", 1)) return(1);
707
708 printf("Linear Solver : ");
709
710 switch(miter) {
711
712 case DENSE_USER :
713 printf("Dense, User-Supplied Jacobian\n");
714
715 /* Create dense SUNMatrix for use in linear solves */
716 *A = SUNDenseMatrix(P1_NEQ, P1_NEQ);
717 if(check_retval((void *)*A, "SUNDenseMatrix", 0)) return(1);
718
719 /* Create dense SUNLinearSolver object for use by CVode */
720 *LS = SUNLinSol_Dense(y, *A);
721 if(check_retval((void *)*LS, "SUNLinSol_Dense", 0)) return(1);
722
723 /* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
724 retval = CVodeSetLinearSolver(cvode_mem, *LS, *A);
725 if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);
726
727 /* Set the user-supplied Jacobian routine Jac */
728 retval = CVodeSetJacFn(cvode_mem, Jac1);
729 if(check_retval(&retval, "CVodeSetJacFn", 1)) return(1);
730 break;
731
732 case DENSE_DQ :
733 printf("Dense, Difference Quotient Jacobian\n");
734
735 /* Create dense SUNMatrix for use in linear solves */
736 *A = SUNDenseMatrix(P1_NEQ, P1_NEQ);
737 if(check_retval((void *)*A, "SUNDenseMatrix", 0)) return(1);
738
739 /* Create dense SUNLinearSolver object for use by CVode */
740 *LS = SUNLinSol_Dense(y, *A);
741 if(check_retval((void *)*LS, "SUNLinSol_Dense", 0)) return(1);
742
743 /* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
744 retval = CVodeSetLinearSolver(cvode_mem, *LS, *A);
745 if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);
746
747 /* Use a difference quotient Jacobian */
748 retval = CVodeSetJacFn(cvode_mem, NULL);
749 if(check_retval(&retval, "CVodeSetJacFn", 1)) return(1);
750 break;
751
752 case DIAG :
753 printf("Diagonal Jacobian\n");
754
755 /* Call CVDiag to create/attach the CVODE-specific diagonal solver */
756 retval = CVDiag(cvode_mem);
757 if(check_retval(&retval, "CVDiag", 1)) return(1);
758 break;
759
760 case BAND_USER :
761 printf("Band, User-Supplied Jacobian\n");
762
763 /* Create band SUNMatrix for use in linear solves */
764 *A = SUNBandMatrix(P2_NEQ, mu, ml);
765 if(check_retval((void *)*A, "SUNBandMatrix", 0)) return(1);
766
767 /* Create banded SUNLinearSolver object for use by CVode */
768 *LS = SUNLinSol_Band(y, *A);
769 if(check_retval((void *)*LS, "SUNLinSol_Band", 0)) return(1);
770
771 /* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
772 retval = CVodeSetLinearSolver(cvode_mem, *LS, *A);
773 if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);
774
775 /* Set the user-supplied Jacobian routine Jac */
776 retval = CVodeSetJacFn(cvode_mem, Jac2);
777 if(check_retval(&retval, "CVodeSetJacFn", 1)) return(1);
778 break;
779
780 case BAND_DQ :
781 printf("Band, Difference Quotient Jacobian\n");
782
783 /* Create band SUNMatrix for use in linear solves */
784 *A = SUNBandMatrix(P2_NEQ, mu, ml);
785 if(check_retval((void *)*A, "SUNBandMatrix", 0)) return(1);
786
787 /* Create banded SUNLinearSolver object for use by CVode */
788 *LS = SUNLinSol_Band(y, *A);
789 if(check_retval((void *)*LS, "SUNLinSol_Band", 0)) return(1);
790
791 /* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
792 retval = CVodeSetLinearSolver(cvode_mem, *LS, *A);
793 if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);
794
795 /* Use a difference quotient Jacobian */
796 retval = CVodeSetJacFn(cvode_mem, NULL);
797 if(check_retval(&retval, "CVodeSetJacFn", 1)) return(1);
798 break;
799 }
800 }
801
802 return(retval);
803 }
804
PrintErrOutput(realtype tol_factor)805 static void PrintErrOutput(realtype tol_factor)
806 {
807 #if defined(SUNDIALS_EXTENDED_PRECISION)
808 printf("\n\n Error exceeds %Lg * tolerance \n\n", tol_factor);
809 #elif defined(SUNDIALS_DOUBLE_PRECISION)
810 printf("\n\n Error exceeds %g * tolerance \n\n", tol_factor);
811 #else
812 printf("\n\n Error exceeds %g * tolerance \n\n", tol_factor);
813 #endif
814
815 return;
816 }
817
PrintFinalStats(void * cvode_mem,int miter,realtype ero)818 static void PrintFinalStats(void *cvode_mem, int miter, realtype ero)
819 {
820 long int lenrw, leniw, lenrwLS, leniwLS;
821 long int nst, nfe, nsetups, nni, ncfn, netf, nje, nfeLS;
822 int retval;
823
824 retval = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);
825 check_retval(&retval, "CVodeGetWorkSpace", 1);
826 retval = CVodeGetNumSteps(cvode_mem, &nst);
827 check_retval(&retval, "CVodeGetNumSteps", 1);
828 retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);
829 check_retval(&retval, "CVodeGetNumRhsEvals", 1);
830 retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
831 check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);
832 retval = CVodeGetNumErrTestFails(cvode_mem, &netf);
833 check_retval(&retval, "CVodeGetNumErrTestFails", 1);
834 retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
835 check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);
836 retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
837 check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);
838
839 printf("\n Final statistics for this run:\n\n");
840 printf(" CVode real workspace length = %4ld \n", lenrw);
841 printf(" CVode integer workspace length = %4ld \n", leniw);
842 printf(" Number of steps = %4ld \n", nst);
843 printf(" Number of f-s = %4ld \n", nfe);
844 printf(" Number of setups = %4ld \n", nsetups);
845 printf(" Number of nonlinear iterations = %4ld \n", nni);
846 printf(" Number of nonlinear convergence failures = %4ld \n", ncfn);
847 printf(" Number of error test failures = %4ld \n\n",netf);
848
849 if (miter != FUNC) {
850 if (miter != DIAG) {
851 retval = CVodeGetNumJacEvals(cvode_mem, &nje);
852 check_retval(&retval, "CVodeGetNumJacEvals", 1);
853 retval = CVodeGetNumLinRhsEvals(cvode_mem, &nfeLS);
854 check_retval(&retval, "CVodeGetNumLinRhsEvals", 1);
855 retval = CVodeGetLinWorkSpace(cvode_mem, &lenrwLS, &leniwLS);
856 check_retval(&retval, "CVodeGetLinWorkSpace", 1);
857 } else {
858 nje = nsetups;
859 retval = CVDiagGetNumRhsEvals(cvode_mem, &nfeLS);
860 check_retval(&retval, "CVDiagGetNumRhsEvals", 1);
861 retval = CVDiagGetWorkSpace(cvode_mem, &lenrwLS, &leniwLS);
862 check_retval(&retval, "CVDiagGetWorkSpace", 1);
863 }
864 printf(" Linear solver real workspace length = %4ld \n", lenrwLS);
865 printf(" Linear solver integer workspace length = %4ld \n", leniwLS);
866 printf(" Number of Jacobian evaluations = %4ld \n", nje);
867 printf(" Number of f evals. in linear solver = %4ld \n\n", nfeLS);
868 }
869
870 #if defined(SUNDIALS_EXTENDED_PRECISION)
871 printf(" Error overrun = %.3Lf \n", ero);
872 #else
873 printf(" Error overrun = %.3f \n", ero);
874 #endif
875 }
876
PrintErrInfo(int nerr)877 static void PrintErrInfo(int nerr)
878 {
879 printf("\n\n-------------------------------------------------------------");
880 printf("\n-------------------------------------------------------------");
881 printf("\n\n Number of errors encountered = %d \n", nerr);
882
883 return;
884 }
885
886 /* Check function return value...
887 opt == 0 means SUNDIALS function allocates memory so check if
888 returned NULL pointer
889 opt == 1 means SUNDIALS function returns an integer value so check if
890 retval < 0
891 opt == 2 means function allocates memory so check if returned
892 NULL pointer */
893
check_retval(void * returnvalue,const char * funcname,int opt)894 static int check_retval(void *returnvalue, const char *funcname, int opt)
895 {
896 int *retval;
897
898 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
899 if (opt == 0 && returnvalue == NULL) {
900 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
901 funcname);
902 return(1); }
903
904 /* Check if retval < 0 */
905 else if (opt == 1) {
906 retval = (int *) returnvalue;
907 if (*retval < 0) {
908 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
909 funcname, *retval);
910 return(1); }}
911
912 /* Check if function returned NULL pointer - no memory allocated */
913 else if (opt == 2 && returnvalue == NULL) {
914 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
915 funcname);
916 return(1); }
917
918 return(0);
919 }
920