1 /* -----------------------------------------------------------------
2 * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and
3 * Radu Serban, Cody J. Balos @ LLNL
4 * -----------------------------------------------------------------
5 * SUNDIALS Copyright Start
6 * Copyright (c) 2002-2021, 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 * This example loops through the available iterative linear solvers:
16 * SPGMR, SPFGMR, SPBCG and SPTFQMR.
17 *
18 * Example problem:
19 *
20 * An ODE system is generated from the following 2-species diurnal
21 * kinetics advection-diffusion PDE system in 2 space dimensions:
22 *
23 * dc(i)/dt = Kh*(d/dx)^2 c(i) + V*dc(i)/dx + (d/dy)(Kv(y)*dc(i)/dy)
24 * + Ri(c1,c2,t) for i = 1,2, where
25 * R1(c1,c2,t) = -q1*c1*c3 - q2*c1*c2 + 2*q3(t)*c3 + q4(t)*c2 ,
26 * R2(c1,c2,t) = q1*c1*c3 - q2*c1*c2 - q4(t)*c2 ,
27 * Kv(y) = Kv0*exp(y/5) ,
28 * Kh, V, Kv0, q1, q2, and c3 are constants, and q3(t) and q4(t)
29 * vary diurnally. The problem is posed on the square
30 * 0 <= x <= 20, 30 <= y <= 50 (all in km),
31 * with homogeneous Neumann boundary conditions, and for time t in
32 * 0 <= t <= 86400 sec (1 day).
33 * The PDE system is treated by central differences on a uniform
34 * 10 x 10 mesh, with simple polynomial initial profiles.
35 * The problem is solved with CVODE, with the BDF/GMRES, BDF/FGMRES
36 * BDF/Bi-CGStab, and BDF/TFQMR methods (i.e. using the SUNLinSol_SPGMR,
37 * SUNLinSol_SPFGMR, SUNLinSol_SPBCGS, and SUNLinSol_SPTFQMR linear solvers)
38 * and the block-diagonal part of the Newton matrix as a left preconditioner.
39 * A copy of the block-diagonal part of the Jacobian is saved and
40 * conditionally reused within the Precond routine.
41 * -----------------------------------------------------------------*/
42
43 #include <stdio.h>
44 #include <stdlib.h>
45 #include <math.h>
46
47 #include <cvode/cvode.h> /* main integrator header file */
48 #include <sunlinsol/sunlinsol_spgmr.h> /* access to SPGMR SUNLinearSolver */
49 #include <sunlinsol/sunlinsol_spfgmr.h> /* access to SPFGMR SUNLinearSolver */
50 #include <sunlinsol/sunlinsol_spbcgs.h> /* access to SPBCGS SUNLinearSolver */
51 #include <sunlinsol/sunlinsol_sptfqmr.h> /* access to SPTFQMR SUNLinearSolver */
52 #include <sunnonlinsol/sunnonlinsol_newton.h> /* access to Newton SUNNonlinearSolver */
53 #include <nvector/nvector_serial.h> /* serial N_Vector types, fct. and macros */
54 #include <sundials/sundials_dense.h> /* use generic DENSE solver in preconditioning */
55 #include <sundials/sundials_types.h> /* definition of realtype */
56
57 /* helpful macros */
58
59 #ifndef SQR
60 #define SQR(A) ((A)*(A))
61 #endif
62
63 #ifndef SQRT
64 #if defined(SUNDIALS_DOUBLE_PRECISION)
65 #define SQRT(x) (sqrt((x)))
66 #elif defined(SUNDIALS_SINGLE_PRECISION)
67 #define SQRT(x) (sqrtf((x)))
68 #elif defined(SUNDIALS_EXTENDED_PRECISION)
69 #define SQRT(x) (sqrtl((x)))
70 #endif
71 #endif
72
73 /* Problem Constants */
74
75 #define ZERO RCONST(0.0)
76 #define ONE RCONST(1.0)
77 #define TWO RCONST(2.0)
78
79 #define NUM_SPECIES 2 /* number of species */
80 #define KH RCONST(4.0e-6) /* horizontal diffusivity Kh */
81 #define VEL RCONST(0.001) /* advection velocity V */
82 #define KV0 RCONST(1.0e-8) /* coefficient in Kv(y) */
83 #define Q1 RCONST(1.63e-16) /* coefficients q1, q2, c3 */
84 #define Q2 RCONST(4.66e-16)
85 #define C3 RCONST(3.7e16)
86 #define A3 RCONST(22.62) /* coefficient in expression for q3(t) */
87 #define A4 RCONST(7.601) /* coefficient in expression for q4(t) */
88 #define C1_SCALE RCONST(1.0e6) /* coefficients in initial profiles */
89 #define C2_SCALE RCONST(1.0e12)
90
91 #define T0 ZERO /* initial time */
92 #define NOUT 12 /* number of output times */
93 #define TWOHR RCONST(7200.0) /* number of seconds in two hours */
94 #define HALFDAY RCONST(4.32e4) /* number of seconds in a half day */
95 #define PI RCONST(3.1415926535898) /* pi */
96
97 #define XMIN ZERO /* grid boundaries in x */
98 #define XMAX RCONST(20.0)
99 #define YMIN RCONST(30.0) /* grid boundaries in y */
100 #define YMAX RCONST(50.0)
101 #define XMID RCONST(10.0) /* grid midpoints in x,y */
102 #define YMID RCONST(40.0)
103
104 #define MX 10 /* MX = number of x mesh points */
105 #define MY 10 /* MY = number of y mesh points */
106 #define NSMX 20 /* NSMX = NUM_SPECIES*MX */
107 #define MM (MX*MY) /* MM = MX*MY */
108
109 /* CVodeInit Constants */
110
111 #define RTOL RCONST(1.0e-5) /* scalar relative tolerance */
112 #define FLOOR RCONST(100.0) /* value of C1 or C2 at which tolerances */
113 /* change from relative to absolute */
114 #define ATOL (RTOL*FLOOR) /* scalar absolute tolerance */
115 #define NEQ (NUM_SPECIES*MM) /* NEQ = number of equations */
116
117 /* Linear Solver Loop Constants */
118
119 #define USE_SPGMR 0
120 #define USE_SPFGMR 1
121 #define USE_SPBCG 2
122 #define USE_SPTFQMR 3
123
124 /* User-defined vector and matrix accessor macros: IJKth, IJth */
125
126 /* IJKth is defined in order to isolate the translation from the
127 mathematical 3-dimensional structure of the dependent variable vector
128 to the underlying 1-dimensional storage. IJth is defined in order to
129 write code which indexes into dense matrices with a (row,column)
130 pair, where 1 <= row, column <= NUM_SPECIES.
131
132 IJKth(vdata,i,j,k) references the element in the vdata array for
133 species i at mesh point (j,k), where 1 <= i <= NUM_SPECIES,
134 0 <= j <= MX-1, 0 <= k <= MY-1. The vdata array is obtained via
135 the call vdata = N_VGetArrayPointer(v), where v is an N_Vector.
136 For each mesh point (j,k), the elements for species i and i+1 are
137 contiguous within vdata.
138
139 IJth(a,i,j) references the (i,j)th entry of the matrix realtype **a,
140 where 1 <= i,j <= NUM_SPECIES. The small matrix routines in
141 sundials_dense.h work with matrices stored by column in a 2-dimensional
142 array. In C, arrays are indexed starting at 0, not 1. */
143
144 #define IJKth(vdata,i,j,k) (vdata[i-1 + (j)*NUM_SPECIES + (k)*NSMX])
145 #define IJth(a,i,j) (a[j-1][i-1])
146
147 /* Type : UserData
148 contains preconditioner blocks, pivot arrays, and problem constants,
149 solution vector, and linsolver type */
150
151 typedef struct {
152 realtype **P[MX][MY], **Jbd[MX][MY];
153 sunindextype *pivot[MX][MY];
154 realtype q4, om, dx, dy, hdco, haco, vdco;
155 N_Vector u;
156 int linsolver;
157 } *UserData;
158
159 /* Private Helper Functions */
160
161 static UserData AllocUserData(void);
162 static void InitUserData(UserData data, N_Vector u);
163 static void FreeUserData(UserData data);
164 static void SetInitialProfiles(N_Vector u, realtype dx, realtype dy);
165 static void PrintOutput(void *cvode_mem, N_Vector u, realtype t);
166 static void PrintStats(void *cvode_mem, int linsolver, int stats);
167 static int check_retval(void *returnvalue, const char *funcname, int opt);
168
169 /* Functions Called by the Solver */
170
171 static int f(realtype t, N_Vector u, N_Vector udot, void *user_data);
172
173 static int Precond(realtype tn, N_Vector u, N_Vector fu, booleantype jok,
174 booleantype *jcurPtr, realtype gamma, void *user_data);
175
176 static int PSolve(realtype tn, N_Vector u, N_Vector fu,
177 N_Vector r, N_Vector z,
178 realtype gamma, realtype delta,
179 int lr, void *user_data);
180
181 static int myMonitorFunction(void *cvode_mem, void *user_data);
182
183 /*
184 *-------------------------------
185 * Main Program
186 *-------------------------------
187 */
188
main(int argc,char * argv[])189 int main(int argc, char* argv[])
190 {
191 realtype abstol, reltol, t, tout;
192 N_Vector u;
193 UserData data;
194 SUNLinearSolver LS;
195 SUNNonlinearSolver NLS;
196 void *cvode_mem;
197 int linsolver, iout, retval;
198 FILE* infofp;
199 int nrmfactor; /* LS norm conversion factor flag */
200 realtype nrmfac; /* LS norm conversion factor */
201 int monitor; /* LS resiudal monitoring flag */
202
203 u = NULL;
204 data = NULL;
205 LS = NULL;
206 cvode_mem = NULL;
207 infofp = NULL;
208 nrmfactor = 0;
209 monitor = 0;
210
211 /* Retrieve the command-line options */
212 if (argc > 1) nrmfactor = atoi(argv[1]);
213 if (argc > 2) monitor = atoi(argv[2]);
214
215 /* Open info file if monitoring is turned on */
216 if (monitor) {
217 infofp = fopen("cvKrylovDemo_ls-info.txt", "w+");
218 if (check_retval((void *)infofp, "fopen", 0)) return(1);
219 }
220
221 /* Allocate memory, and set problem data, initial values, tolerances */
222 u = N_VNew_Serial(NEQ);
223 if (check_retval((void *)u, "N_VNew_Serial", 0)) return(1);
224 data = AllocUserData();
225 if (check_retval((void *)data, "AllocUserData", 2)) return(1);
226 InitUserData(data, u);
227 SetInitialProfiles(u, data->dx, data->dy);
228 abstol=ATOL;
229 reltol=RTOL;
230
231 /* Call CVodeCreate to create the solver memory and specify the
232 * Backward Differentiation Formula */
233 cvode_mem = CVodeCreate(CV_BDF);
234 if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
235
236 /* Set the pointer to user-defined data */
237 retval = CVodeSetUserData(cvode_mem, data);
238 if (check_retval(&retval, "CVodeSetUserData", 1)) return(1);
239
240 /* Call CVodeInit to initialize the integrator memory and specify the
241 * user's right hand side function in u'=f(t,u), the inital time T0, and
242 * the initial dependent variable vector u. */
243 retval = CVodeInit(cvode_mem, f, T0, u);
244 if (check_retval(&retval, "CVodeInit", 1)) return(1);
245
246 /* Call CVodeSStolerances to specify the scalar relative tolerance
247 * and scalar absolute tolerances */
248 retval = CVodeSStolerances(cvode_mem, reltol, abstol);
249 if (check_retval(&retval, "CVodeSStolerances", 1)) return(1);
250
251 /* Set a function that CVode will call every 50 successful time steps.
252 * This will be used to monitor the solution and integrator statistics. */
253 if (monitor) {
254 retval = CVodeSetMonitorFn(cvode_mem, myMonitorFunction);
255 if (check_retval(&retval, "CVodeSetMonitorFn", 1)) return(1);
256 retval = CVodeSetMonitorFrequency(cvode_mem, 50);
257 if (check_retval(&retval, "CVodeSetMonitorFrequency", 1)) return(1);
258 }
259
260 /* Create the SUNNonlinearSolver */
261 NLS = SUNNonlinSol_Newton(u);
262 if (check_retval(&retval, "SUNNonlinSol_Newton", 0)) return(1);
263 if (monitor) {
264 /* Set the print level set to 1, so that the nonlinear residual
265 is printed every newton iteration. */
266 retval = SUNNonlinSolSetPrintLevel_Newton(NLS, 1);
267 if (check_retval(&retval, "SUNNonlinSolSetPrintLevel_Newton", 1)) return(1);
268 retval = SUNNonlinSolSetInfoFile_Newton(NLS, infofp);
269 if (check_retval(&retval, "SUNNonlinSolSetInfoFile_Newton", 1)) return(1);
270 }
271
272 /* Call CVodeSetNonlinearSolver to attach the nonlinear solver to CVode */
273 retval = CVodeSetNonlinearSolver(cvode_mem, NLS);
274 if (check_retval(&retval, "CVodeSetNonlinearSolver", 1)) return(1);
275
276 /* START: Loop through SPGMR, SPFGMR, SPBCG and SPTFQMR linear solver modules */
277 for (linsolver = 0; linsolver < 4; ++linsolver) {
278
279 if (linsolver != 0) {
280
281 /* Re-initialize user data */
282 InitUserData(data, u);
283 SetInitialProfiles(u, data->dx, data->dy);
284
285 /* Re-initialize CVode for the solution of the same problem, but
286 using a different linear solver module */
287 retval = CVodeReInit(cvode_mem, T0, u);
288 if (check_retval(&retval, "CVodeReInit", 1)) return(1);
289
290 }
291
292 /* Free previous linear solver and attach a new linear solver module */
293 SUNLinSolFree(LS);
294
295 /* Set the linear sovler type in user data */
296 data->linsolver = linsolver;
297
298 switch(linsolver) {
299
300 /* (a) SPGMR */
301 case(USE_SPGMR):
302
303 /* Print header */
304 printf(" -------");
305 printf(" \n| SPGMR |\n");
306 printf(" -------\n");
307 if (monitor) {
308 fprintf(infofp, " ---------");
309 fprintf(infofp, " \n| SPGMR |\n");
310 fprintf(infofp, " ---------\n");
311 }
312
313 /* Call SUNLinSol_SPGMR to specify the linear solver SPGMR with
314 left preconditioning and the default maximum Krylov dimension */
315 LS = SUNLinSol_SPGMR(u, PREC_LEFT, 0);
316 if (check_retval((void *)LS, "SUNLinSol_SPGMR", 0)) return(1);
317 if (monitor) {
318 retval = SUNLinSolSetPrintLevel_SPGMR(LS, 1);
319 if (check_retval(&retval, "SUNLinSolSetPrintLevel_SPGMR", 1)) return(1);
320 retval = SUNLinSolSetInfoFile_SPGMR(LS, infofp);
321 if (check_retval(&retval, "SUNLinSolSetInfoFile_SPGMR", 1)) return(1);
322 }
323 retval = CVodeSetLinearSolver(cvode_mem, LS, NULL);
324 if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1;
325
326 break;
327
328 /* (b) SPFGMR */
329 case(USE_SPFGMR):
330
331 /* Print header */
332 printf(" ---------");
333 printf(" \n| SPFGMR |\n");
334 printf(" ---------\n");
335 if (monitor) {
336 fprintf(infofp, " ---------");
337 fprintf(infofp, " \n| SPFGMR |\n");
338 fprintf(infofp, " ---------\n");
339 }
340
341 /* Call SUNLinSol_SPFGMR to specify the linear solver SPFGMR with
342 left preconditioning and the default maximum Krylov dimension */
343 LS = SUNLinSol_SPFGMR(u, PREC_LEFT, 0);
344 if (check_retval((void *)LS, "SUNLinSol_SPFGMR", 0)) return(1);
345 if (monitor) {
346 retval = SUNLinSolSetPrintLevel_SPFGMR(LS, 1);
347 if (check_retval(&retval, "SUNLinSolSetPrintLevel_SPFGMR", 1)) return(1);
348 retval = SUNLinSolSetInfoFile_SPFGMR(LS, infofp);
349 if (check_retval(&retval, "SUNLinSolSetInfoFile_SPFGMR", 1)) return(1);
350 }
351 retval = CVodeSetLinearSolver(cvode_mem, LS, NULL);
352 if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1;
353
354 break;
355
356 /* (c) SPBCG */
357 case(USE_SPBCG):
358
359 /* Print header */
360 printf(" -------");
361 printf(" \n| SPBCGS |\n");
362 printf(" -------\n");
363 if (monitor) {
364 fprintf(infofp, " ---------");
365 fprintf(infofp, " \n| SPBCGS |\n");
366 fprintf(infofp, " ---------\n");
367 }
368
369 /* Call SUNLinSol_SPBCGS to specify the linear solver SPBCGS with
370 left preconditioning and the default maximum Krylov dimension */
371 LS = SUNLinSol_SPBCGS(u, PREC_LEFT, 0);
372 if (check_retval((void *)LS, "SUNLinSol_SPBCGS", 0)) return(1);
373 if (monitor) {
374 retval = SUNLinSolSetPrintLevel_SPBCGS(LS, 1);
375 if (check_retval(&retval, "SUNLinSolSetPrintLevel_SPBCGS", 1)) return(1);
376 retval = SUNLinSolSetInfoFile_SPBCGS(LS, infofp);
377 if (check_retval(&retval, "SUNLinSolSetInfoFile_SPBCGS", 1)) return(1);
378 }
379 retval = CVodeSetLinearSolver(cvode_mem, LS, NULL);
380 if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1;
381
382 break;
383
384 /* (d) SPTFQMR */
385 case(USE_SPTFQMR):
386
387 /* Print header */
388 printf(" ---------");
389 printf(" \n| SPTFQMR |\n");
390 printf(" ---------\n");
391 if (monitor) {
392 fprintf(infofp, " ---------");
393 fprintf(infofp, " \n| SPTFQMR |\n");
394 fprintf(infofp, " ---------\n");
395 }
396
397 /* Call SUNLinSol_SPTFQMR to specify the linear solver SPTFQMR with
398 left preconditioning and the default maximum Krylov dimension */
399 LS = SUNLinSol_SPTFQMR(u, PREC_LEFT, 0);
400 if (check_retval((void *)LS, "SUNLinSol_SPTFQMR", 0)) return(1);
401 if (monitor) {
402 retval = SUNLinSolSetPrintLevel_SPTFQMR(LS, 1);
403 if (check_retval(&retval, "SUNLinSolSetPrintLevel_SPTFQMR", 1)) return(1);
404 retval = SUNLinSolSetInfoFile_SPTFQMR(LS, infofp);
405 if (check_retval(&retval, "SUNLinSolSetInfoFile_SPTFQMR", 1)) return(1);
406 }
407 retval = CVodeSetLinearSolver(cvode_mem, LS, NULL);
408 if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1;
409
410 break;
411 }
412
413
414 /* Set preconditioner setup and solve routines Precond and PSolve,
415 and the pointer to the user-defined block data */
416 retval = CVodeSetPreconditioner(cvode_mem, Precond, PSolve);
417 if (check_retval(&retval, "CVodeSetPreconditioner", 1)) return(1);
418
419 /* Set the linear solver tolerance conversion factor */
420 switch(nrmfactor) {
421
422 case(1):
423 /* use the square root of the vector length */
424 nrmfac = SQRT((realtype)NEQ);
425 break;
426 case(2):
427 /* compute with dot product */
428 nrmfac = -ONE;
429 break;
430 default:
431 /* use the default */
432 nrmfac = ZERO;
433 break;
434 }
435
436 retval = CVodeSetLSNormFactor(cvode_mem, nrmfac);
437 if (check_retval(&retval, "CVodeSetLSNormFactor", 1)) return(1);
438
439 /* In loop over output points, call CVode, print results, and test for error */
440 printf(" \n2-species diurnal advection-diffusion problem\n\n");
441 for (iout=1, tout = TWOHR; iout <= NOUT; iout++, tout += TWOHR) {
442 retval = CVode(cvode_mem, tout, u, &t, CV_NORMAL);
443 if (!monitor) PrintOutput(cvode_mem, u, t);
444 if (check_retval(&retval, "CVode", 1)) break;
445 }
446 if (monitor) PrintOutput(cvode_mem, u, t);
447 PrintStats(cvode_mem, linsolver, 1);
448
449 } /* END: Loop through SPGMR, SPBCG and SPTFQMR linear solver modules */
450
451 /* Free memory */
452 if (monitor) fclose(infofp);
453 N_VDestroy(u);
454 FreeUserData(data);
455 CVodeFree(&cvode_mem);
456 SUNLinSolFree(LS);
457
458 return(0);
459 }
460
461 /*
462 *-------------------------------
463 * Private helper functions
464 *-------------------------------
465 */
466
467 /* Allocate memory for data structure of type UserData */
468
AllocUserData(void)469 static UserData AllocUserData(void)
470 {
471 int jx, jy;
472 UserData data;
473
474 data = (UserData) malloc(sizeof *data);
475
476 for (jx=0; jx < MX; jx++) {
477 for (jy=0; jy < MY; jy++) {
478 (data->P)[jx][jy] = newDenseMat(NUM_SPECIES, NUM_SPECIES);
479 (data->Jbd)[jx][jy] = newDenseMat(NUM_SPECIES, NUM_SPECIES);
480 (data->pivot)[jx][jy] = newIndexArray(NUM_SPECIES);
481 }
482 }
483
484 return(data);
485 }
486
487 /* Load problem constants in data */
488
InitUserData(UserData data,N_Vector u)489 static void InitUserData(UserData data, N_Vector u)
490 {
491 data->u = u;
492 data->om = PI/HALFDAY;
493 data->dx = (XMAX-XMIN)/(MX-1);
494 data->dy = (YMAX-YMIN)/(MY-1);
495 data->hdco = KH/SQR(data->dx);
496 data->haco = VEL/(TWO*data->dx);
497 data->vdco = (ONE/SQR(data->dy))*KV0;
498 }
499
500 /* Free data memory */
501
FreeUserData(UserData data)502 static void FreeUserData(UserData data)
503 {
504 int jx, jy;
505
506 for (jx=0; jx < MX; jx++) {
507 for (jy=0; jy < MY; jy++) {
508 destroyMat((data->P)[jx][jy]);
509 destroyMat((data->Jbd)[jx][jy]);
510 destroyArray((data->pivot)[jx][jy]);
511 }
512 }
513
514 free(data);
515 }
516
517 /* Set initial conditions in u */
518
SetInitialProfiles(N_Vector u,realtype dx,realtype dy)519 static void SetInitialProfiles(N_Vector u, realtype dx, realtype dy)
520 {
521 int jx, jy;
522 realtype x, y, cx, cy;
523 realtype *udata;
524
525 /* Set pointer to data array in vector u. */
526
527 udata = N_VGetArrayPointer(u);
528
529 /* Load initial profiles of c1 and c2 into u vector */
530
531 for (jy=0; jy < MY; jy++) {
532 y = YMIN + jy*dy;
533 cy = SQR(RCONST(0.1)*(y - YMID));
534 cy = ONE - cy + RCONST(0.5)*SQR(cy);
535 for (jx=0; jx < MX; jx++) {
536 x = XMIN + jx*dx;
537 cx = SQR(RCONST(0.1)*(x - XMID));
538 cx = ONE - cx + RCONST(0.5)*SQR(cx);
539 IJKth(udata,1,jx,jy) = C1_SCALE*cx*cy;
540 IJKth(udata,2,jx,jy) = C2_SCALE*cx*cy;
541 }
542 }
543 }
544
545 /* Print current t, step count, order, stepsize, and sampled c1,c2 values */
546
PrintOutput(void * cvode_mem,N_Vector u,realtype t)547 static void PrintOutput(void *cvode_mem, N_Vector u, realtype t)
548 {
549 long int nst;
550 int qu, retval;
551 realtype hu, *udata;
552 int mxh = MX/2 - 1, myh = MY/2 - 1, mx1 = MX - 1, my1 = MY - 1;
553
554 udata = N_VGetArrayPointer(u);
555
556 retval = CVodeGetNumSteps(cvode_mem, &nst);
557 check_retval(&retval, "CVodeGetNumSteps", 1);
558 retval = CVodeGetLastOrder(cvode_mem, &qu);
559 check_retval(&retval, "CVodeGetLastOrder", 1);
560 retval = CVodeGetLastStep(cvode_mem, &hu);
561 check_retval(&retval, "CVodeGetLastStep", 1);
562
563 #if defined(SUNDIALS_EXTENDED_PRECISION)
564 printf("t = %.2Le no. steps = %ld order = %d stepsize = %.2Le\n",
565 t, nst, qu, hu);
566 printf("c1 (bot.left/middle/top rt.) = %12.3Le %12.3Le %12.3Le\n",
567 IJKth(udata,1,0,0), IJKth(udata,1,mxh,myh), IJKth(udata,1,mx1,my1));
568 printf("c2 (bot.left/middle/top rt.) = %12.3Le %12.3Le %12.3Le\n\n",
569 IJKth(udata,2,0,0), IJKth(udata,2,mxh,myh), IJKth(udata,2,mx1,my1));
570 #elif defined(SUNDIALS_DOUBLE_PRECISION)
571 printf("t = %.2e no. steps = %ld order = %d stepsize = %.2e\n",
572 t, nst, qu, hu);
573 printf("c1 (bot.left/middle/top rt.) = %12.3e %12.3e %12.3e\n",
574 IJKth(udata,1,0,0), IJKth(udata,1,mxh,myh), IJKth(udata,1,mx1,my1));
575 printf("c2 (bot.left/middle/top rt.) = %12.3e %12.3e %12.3e\n\n",
576 IJKth(udata,2,0,0), IJKth(udata,2,mxh,myh), IJKth(udata,2,mx1,my1));
577 #else
578 printf("t = %.2e no. steps = %ld order = %d stepsize = %.2e\n",
579 t, nst, qu, hu);
580 printf("c1 (bot.left/middle/top rt.) = %12.3e %12.3e %12.3e\n",
581 IJKth(udata,1,0,0), IJKth(udata,1,mxh,myh), IJKth(udata,1,mx1,my1));
582 printf("c2 (bot.left/middle/top rt.) = %12.3e %12.3e %12.3e\n\n",
583 IJKth(udata,2,0,0), IJKth(udata,2,mxh,myh), IJKth(udata,2,mx1,my1));
584 #endif
585 }
586
587 /* Get and print final statistics */
588
PrintStats(void * cvode_mem,int linsolver,int final)589 static void PrintStats(void *cvode_mem, int linsolver, int final)
590 {
591 long int lenrw, leniw ;
592 long int lenrwLS, leniwLS;
593 long int nst, nfe, nsetups, nni, ncfn, netf;
594 long int nje, nli, npe, nps, ncfl, nfeLS, njts, njte;
595 int retval;
596
597 retval = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);
598 check_retval(&retval, "CVodeGetWorkSpace", 1);
599 retval = CVodeGetNumSteps(cvode_mem, &nst);
600 check_retval(&retval, "CVodeGetNumSteps", 1);
601 retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);
602 check_retval(&retval, "CVodeGetNumRhsEvals", 1);
603 retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
604 check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);
605 retval = CVodeGetNumErrTestFails(cvode_mem, &netf);
606 check_retval(&retval, "CVodeGetNumErrTestFails", 1);
607 retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
608 check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);
609 retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
610 check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);
611
612 retval = CVodeGetLinWorkSpace(cvode_mem, &lenrwLS, &leniwLS);
613 check_retval(&retval, "CVodeGetLinWorkSpace", 1);
614 CVodeGetLinSolveStats(cvode_mem, &nje, &nfeLS, &nli, &ncfl, &npe,
615 &nps, &njts, &njte);
616 check_retval(&retval, "CVodeGetLinWorkSpace", 1);
617
618 if (final) printf("\nFinal Statistics.. \n\n");
619 else printf("\nIntermediate Statistics.. \n\n");
620 printf("lenrw = %5ld leniw = %5ld\n" , lenrw, leniw);
621 printf("lenrwLS = %5ld leniwLS = %5ld\n" , lenrwLS, leniwLS);
622 printf("nst = %5ld\n" , nst);
623 printf("nfe = %5ld nfeLS = %5ld\n" , nfe, nfeLS);
624 printf("nni = %5ld nli = %5ld\n" , nni, nli);
625 printf("nsetups = %5ld netf = %5ld\n" , nsetups, netf);
626 printf("npe = %5ld nps = %5ld\n" , npe, nps);
627 printf("ncfn = %5ld ncfl = %5ld\n\n", ncfn, ncfl);
628
629 if (linsolver < 2)
630 printf("======================================================================\n\n");
631 }
632
633 /* Check function return value...
634 opt == 0 means SUNDIALS function allocates memory so check if
635 returned NULL pointer
636 opt == 1 means SUNDIALS function returns an integer value so check if
637 retval < 0
638 opt == 2 means function allocates memory so check if returned
639 NULL pointer */
640
check_retval(void * returnvalue,const char * funcname,int opt)641 static int check_retval(void *returnvalue, const char *funcname, int opt)
642 {
643 int *retval;
644
645 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
646 if (opt == 0 && returnvalue == NULL) {
647 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
648 funcname);
649 return(1); }
650
651 /* Check if retval < 0 */
652 else if (opt == 1) {
653 retval = (int *) returnvalue;
654 if (*retval < 0) {
655 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
656 funcname, *retval);
657 return(1); }}
658
659 /* Check if function returned NULL pointer - no memory allocated */
660 else if (opt == 2 && returnvalue == NULL) {
661 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
662 funcname);
663 return(1); }
664
665 return(0);
666 }
667
668 /*
669 *-------------------------------
670 * Functions called by the solver
671 *-------------------------------
672 */
673
674 /* f routine. Compute RHS function f(t,u). */
675
f(realtype t,N_Vector u,N_Vector udot,void * user_data)676 static int f(realtype t, N_Vector u, N_Vector udot, void *user_data)
677 {
678 realtype q3, c1, c2, c1dn, c2dn, c1up, c2up, c1lt, c2lt;
679 realtype c1rt, c2rt, cydn, cyup, hord1, hord2, horad1, horad2;
680 realtype qq1, qq2, qq3, qq4, rkin1, rkin2, s, vertd1, vertd2, ydn, yup;
681 realtype q4coef, dely, verdco, hordco, horaco;
682 realtype *udata, *dudata;
683 int jx, jy, idn, iup, ileft, iright;
684 UserData data;
685
686 data = (UserData) user_data;
687 udata = N_VGetArrayPointer(u);
688 dudata = N_VGetArrayPointer(udot);
689
690 /* Set diurnal rate coefficients. */
691
692 s = sin(data->om*t);
693 if (s > ZERO) {
694 q3 = exp(-A3/s);
695 data->q4 = exp(-A4/s);
696 } else {
697 q3 = ZERO;
698 data->q4 = ZERO;
699 }
700
701 /* Make local copies of problem variables, for efficiency. */
702
703 q4coef = data->q4;
704 dely = data->dy;
705 verdco = data->vdco;
706 hordco = data->hdco;
707 horaco = data->haco;
708
709 /* Loop over all grid points. */
710
711 for (jy=0; jy < MY; jy++) {
712
713 /* Set vertical diffusion coefficients at jy +- 1/2 */
714
715 ydn = YMIN + (jy - RCONST(0.5))*dely;
716 yup = ydn + dely;
717 cydn = verdco*exp(RCONST(0.2)*ydn);
718 cyup = verdco*exp(RCONST(0.2)*yup);
719 idn = (jy == 0) ? 1 : -1;
720 iup = (jy == MY-1) ? -1 : 1;
721 for (jx=0; jx < MX; jx++) {
722
723 /* Extract c1 and c2, and set kinetic rate terms. */
724
725 c1 = IJKth(udata,1,jx,jy);
726 c2 = IJKth(udata,2,jx,jy);
727 qq1 = Q1*c1*C3;
728 qq2 = Q2*c1*c2;
729 qq3 = q3*C3;
730 qq4 = q4coef*c2;
731 rkin1 = -qq1 - qq2 + TWO*qq3 + qq4;
732 rkin2 = qq1 - qq2 - qq4;
733
734 /* Set vertical diffusion terms. */
735
736 c1dn = IJKth(udata,1,jx,jy+idn);
737 c2dn = IJKth(udata,2,jx,jy+idn);
738 c1up = IJKth(udata,1,jx,jy+iup);
739 c2up = IJKth(udata,2,jx,jy+iup);
740 vertd1 = cyup*(c1up - c1) - cydn*(c1 - c1dn);
741 vertd2 = cyup*(c2up - c2) - cydn*(c2 - c2dn);
742
743 /* Set horizontal diffusion and advection terms. */
744
745 ileft = (jx == 0) ? 1 : -1;
746 iright =(jx == MX-1) ? -1 : 1;
747 c1lt = IJKth(udata,1,jx+ileft,jy);
748 c2lt = IJKth(udata,2,jx+ileft,jy);
749 c1rt = IJKth(udata,1,jx+iright,jy);
750 c2rt = IJKth(udata,2,jx+iright,jy);
751 hord1 = hordco*(c1rt - TWO*c1 + c1lt);
752 hord2 = hordco*(c2rt - TWO*c2 + c2lt);
753 horad1 = horaco*(c1rt - c1lt);
754 horad2 = horaco*(c2rt - c2lt);
755
756 /* Load all terms into udot. */
757
758 IJKth(dudata, 1, jx, jy) = vertd1 + hord1 + horad1 + rkin1;
759 IJKth(dudata, 2, jx, jy) = vertd2 + hord2 + horad2 + rkin2;
760 }
761 }
762
763 return(0);
764 }
765
766 /* Preconditioner setup routine. Generate and preprocess P. */
767
Precond(realtype tn,N_Vector u,N_Vector fu,booleantype jok,booleantype * jcurPtr,realtype gamma,void * user_data)768 static int Precond(realtype tn, N_Vector u, N_Vector fu, booleantype jok,
769 booleantype *jcurPtr, realtype gamma, void *user_data)
770 {
771 realtype c1, c2, cydn, cyup, diag, ydn, yup, q4coef, dely, verdco, hordco;
772 realtype **(*P)[MY], **(*Jbd)[MY];
773 sunindextype *(*pivot)[MY], retval;
774 int jx, jy;
775 realtype *udata, **a, **j;
776 UserData data;
777
778 /* Make local copies of pointers in user_data, and of pointer to u's data */
779
780 data = (UserData) user_data;
781 P = data->P;
782 Jbd = data->Jbd;
783 pivot = data->pivot;
784 udata = N_VGetArrayPointer(u);
785
786 if (jok) {
787
788 /* jok = SUNTRUE: Copy Jbd to P */
789
790 for (jy=0; jy < MY; jy++)
791 for (jx=0; jx < MX; jx++)
792 denseCopy(Jbd[jx][jy], P[jx][jy], NUM_SPECIES, NUM_SPECIES);
793
794 *jcurPtr = SUNFALSE;
795
796 }
797
798 else {
799 /* jok = SUNFALSE: Generate Jbd from scratch and copy to P */
800
801 /* Make local copies of problem variables, for efficiency. */
802
803 q4coef = data->q4;
804 dely = data->dy;
805 verdco = data->vdco;
806 hordco = data->hdco;
807
808 /* Compute 2x2 diagonal Jacobian blocks (using q4 values
809 computed on the last f call). Load into P. */
810
811 for (jy=0; jy < MY; jy++) {
812 ydn = YMIN + (jy - RCONST(0.5))*dely;
813 yup = ydn + dely;
814 cydn = verdco*exp(RCONST(0.2)*ydn);
815 cyup = verdco*exp(RCONST(0.2)*yup);
816 diag = -(cydn + cyup + TWO*hordco);
817 for (jx=0; jx < MX; jx++) {
818 c1 = IJKth(udata,1,jx,jy);
819 c2 = IJKth(udata,2,jx,jy);
820 j = Jbd[jx][jy];
821 a = P[jx][jy];
822 IJth(j,1,1) = (-Q1*C3 - Q2*c2) + diag;
823 IJth(j,1,2) = -Q2*c1 + q4coef;
824 IJth(j,2,1) = Q1*C3 - Q2*c2;
825 IJth(j,2,2) = (-Q2*c1 - q4coef) + diag;
826 denseCopy(j, a, NUM_SPECIES, NUM_SPECIES);
827 }
828 }
829
830 *jcurPtr = SUNTRUE;
831
832 }
833
834 /* Scale by -gamma */
835
836 for (jy=0; jy < MY; jy++)
837 for (jx=0; jx < MX; jx++)
838 denseScale(-gamma, P[jx][jy], NUM_SPECIES, NUM_SPECIES);
839
840 /* Add identity matrix and do LU decompositions on blocks in place. */
841
842 for (jx=0; jx < MX; jx++) {
843 for (jy=0; jy < MY; jy++) {
844 denseAddIdentity(P[jx][jy], NUM_SPECIES);
845 retval =denseGETRF(P[jx][jy], NUM_SPECIES, NUM_SPECIES, pivot[jx][jy]);
846 if (retval != 0) return(1);
847 }
848 }
849
850 return(0);
851 }
852
853 /* Preconditioner solve routine */
854
PSolve(realtype tn,N_Vector u,N_Vector fu,N_Vector r,N_Vector z,realtype gamma,realtype delta,int lr,void * user_data)855 static int PSolve(realtype tn, N_Vector u, N_Vector fu,
856 N_Vector r, N_Vector z,
857 realtype gamma, realtype delta,
858 int lr, void *user_data)
859 {
860 realtype **(*P)[MY];
861 sunindextype *(*pivot)[MY];
862 int jx, jy;
863 realtype *zdata, *v;
864 UserData data;
865
866 /* Extract the P and pivot arrays from user_data. */
867
868 data = (UserData) user_data;
869 P = data->P;
870 pivot = data->pivot;
871 zdata = N_VGetArrayPointer(z);
872
873 N_VScale(ONE, r, z);
874
875 /* Solve the block-diagonal system Px = r using LU factors stored
876 in P and pivot data in pivot, and return the solution in z. */
877
878 for (jx=0; jx < MX; jx++) {
879 for (jy=0; jy < MY; jy++) {
880 v = &(IJKth(zdata, 1, jx, jy));
881 denseGETRS(P[jx][jy], NUM_SPECIES, pivot[jx][jy], v);
882 }
883 }
884
885 return(0);
886 }
887
888
889 /* Function that is called at some step interval by CVODE */
890
myMonitorFunction(void * cvode_mem,void * user_data)891 static int myMonitorFunction(void* cvode_mem, void* user_data)
892 {
893 UserData data = (UserData) user_data;
894 realtype t = 0;
895
896 CVodeGetCurrentTime(cvode_mem, &t);
897 PrintOutput(cvode_mem, data->u, t);
898 PrintStats(cvode_mem, data->linsolver, 0);
899
900 return(0);
901 }
902