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