1 /* -----------------------------------------------------------------
2 * Programmer(s): Radu Serban @ 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 * This program solves a stiff ODE system that arises from a system
15 * of partial differential equations. The PDE system is a food web
16 * population model, with predator-prey interaction and diffusion on
17 * the unit square in two dimensions. The dependent variable vector
18 * is the following:
19 *
20 * 1 2 ns
21 * c = (c , c , ..., c )
22 *
23 * and the PDEs are as follows:
24 *
25 * i i i
26 * dc /dt = d(i)*(c + c ) + f (x,y,c) (i=1,...,ns)
27 * xx yy i
28 *
29 * where
30 *
31 * i ns j
32 * f (x,y,c) = c *(b(i) + sum a(i,j)*c )
33 * i j=1
34 *
35 * The number of species is ns = 2*np, with the first np being prey
36 * and the last np being predators. The coefficients a(i,j), b(i),
37 * d(i) are:
38 *
39 * a(i,i) = -a (all i)
40 * a(i,j) = -g (i <= np, j > np)
41 * a(i,j) = e (i > np, j <= np)
42 * b(i) = b*(1 + alpha*x*y) (i <= np)
43 * b(i) = -b*(1 + alpha*x*y) (i > np)
44 * d(i) = Dprey (i <= np)
45 * d(i) = Dpred (i > np)
46 *
47 * The spatial domain is the unit square. The final time is 10.
48 * The boundary conditions are: normal derivative = 0.
49 * A polynomial in x and y is used to set the initial conditions.
50 *
51 * The PDEs are discretized by central differencing on an MX by
52 * MY mesh. The resulting ODE system is stiff.
53 *
54 * The ODE system is solved by CVODES using Newton iteration and
55 * the SUNLinSol_SPGMR linear solver (scaled preconditioned GMRES).
56 *
57 * The preconditioner matrix used is the product of two matrices:
58 * (1) A matrix, only defined implicitly, based on a fixed number
59 * of Gauss-Seidel iterations using the diffusion terms only.
60 * (2) A block-diagonal matrix based on the partial derivatives
61 * of the interaction terms f only, using block-grouping (computing
62 * only a subset of the ns by ns blocks).
63 *
64 * Additionally, CVODES can integrate backwards in time the
65 * the semi-discrete form of the adjoint PDE:
66 * d(lambda)/dt = - D^T ( lambda_xx + lambda_yy )
67 * - F_c^T lambda
68 * with homogeneous Neumann boundary conditions and final conditions
69 * lambda(x,y,t=t_final) = - g_c^T(t_final)
70 * whose solution at t = 0 represents the sensitivity of
71 * int_x int _y g(t_final,c) dx dy dt
72 * with respect to the initial conditions of the original problem.
73 *
74 * In this example,
75 * g(t,c) = c(ISPEC), with ISPEC defined below.
76 * -----------------------------------------------------------------
77 * Reference: Peter N. Brown and Alan C. Hindmarsh, Reduced Storage
78 * Matrix Methods in Stiff ODE Systems, J. Appl. Math. & Comp., 31
79 * (1989), pp. 40-91. Also available as Lawrence Livermore National
80 * Laboratory Report UCRL-95088, Rev. 1, June 1987.
81 * -----------------------------------------------------------------
82 */
83
84 #include <stdio.h>
85 #include <stdlib.h>
86 #include <math.h>
87
88 #include <cvodes/cvodes.h> /* main integrator header file */
89 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
90 #include <sunlinsol/sunlinsol_spgmr.h> /* access to SPGMR SUNLinearSolver */
91 #include <sundials/sundials_dense.h> /* use generic dense solver in precond. */
92 #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
93
94 /* helpful macros */
95
96 #ifndef MAX
97 #define MAX(A, B) ((A) > (B) ? (A) : (B))
98 #endif
99
100 #ifndef SQR
101 #define SQR(A) ((A)*(A))
102 #endif
103
104 /* Constants */
105
106 #define ZERO RCONST(0.0)
107 #define ONE RCONST(1.0)
108 #define TWO RCONST(2.0)
109
110 /* Problem Specification Constants */
111
112 #define AA ONE /* AA = a */
113 #define EE RCONST(1.0e4) /* EE = e */
114 #define GG RCONST(0.5e-6) /* GG = g */
115 #define BB ONE /* BB = b */
116 #define DPREY ONE
117 #define DPRED RCONST(0.5)
118 #define ALPH ONE
119 #define NP 3
120 #define NS (2*NP)
121
122 /* Method Constants */
123
124 #define MX 20
125 #define MY 20
126 #define MXNS (MX*NS)
127 #define AX ONE
128 #define AY ONE
129 #define DX (AX/(realtype)(MX-1))
130 #define DY (AY/(realtype)(MY-1))
131 #define MP NS
132 #define MQ (MX*MY)
133 #define MXMP (MX*MP)
134 #define NGX 2
135 #define NGY 2
136 #define NGRP (NGX*NGY)
137 #define ITMAX 5
138
139 /* CVodeInit Constants */
140
141 #define NEQ (NS*MX*MY)
142 #define T0 ZERO
143 #define RTOL RCONST(1.0e-5)
144 #define ATOL RCONST(1.0e-5)
145
146 /* Output Constants */
147
148 #define TOUT RCONST(10.0)
149
150 /* Note: The value for species i at mesh point (j,k) is stored in */
151 /* component number (i-1) + j*NS + k*NS*MX of an N_Vector, */
152 /* where 1 <= i <= NS, 0 <= j < MX, 0 <= k < MY. */
153
154 /* Structure for user data */
155
156 typedef struct {
157 realtype **P[NGRP];
158 sunindextype *pivot[NGRP];
159 int ns, mxns, mp, mq, mx, my, ngrp, ngx, ngy, mxmp;
160 int jgx[NGX+1], jgy[NGY+1], jigx[MX], jigy[MY];
161 int jxr[NGX], jyr[NGY];
162 realtype acoef[NS][NS], bcoef[NS], diff[NS];
163 realtype cox[NS], coy[NS], dx, dy, srur;
164 realtype fsave[NEQ];
165 realtype fBsave[NEQ];
166 N_Vector rewt;
167 N_Vector vtemp;
168 void *cvode_mem;
169 int indexB;
170 } *WebData;
171
172 /* Adjoint calculation constants */
173 /* g = int_x int_y c(ISPEC) dy dx at t = Tfinal */
174 #define NSTEPS 80 /* check points every NSTEPS steps */
175 #define ISPEC 6 /* species # in objective */
176
177 /* Prototypes for user-supplied functions */
178
179 static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
180
181 static int Precond(realtype t, N_Vector c, N_Vector fc,
182 booleantype jok, booleantype *jcurPtr,
183 realtype gamma, void *user_data);
184
185 static int PSolve(realtype t, N_Vector c, N_Vector fc,
186 N_Vector r, N_Vector z,
187 realtype gamma, realtype delta,
188 int lr, void *user_data);
189
190 static int fB(realtype t, N_Vector c, N_Vector cB,
191 N_Vector cBdot, void *user_data);
192
193 static int PrecondB(realtype t, N_Vector c,
194 N_Vector cB, N_Vector fcB, booleantype jok,
195 booleantype *jcurPtr, realtype gamma,
196 void *user_data);
197
198 static int PSolveB(realtype t, N_Vector c,
199 N_Vector cB, N_Vector fcB,
200 N_Vector r, N_Vector z,
201 realtype gamma, realtype delta,
202 int lr, void *user_data);
203
204 /* Prototypes for private functions */
205
206 static WebData AllocUserData(void);
207 static void InitUserData(WebData wdata);
208 static void SetGroups(int m, int ng, int jg[], int jig[], int jr[]);
209 static void CInit(N_Vector c, WebData wdata);
210 static void CbInit(N_Vector c, int is, WebData wdata);
211 static void PrintOutput(N_Vector c, int ns, int mxns, WebData wdata);
212 static void FreeUserData(WebData wdata);
213 static void WebRates(realtype x, realtype y, realtype t, realtype c[], realtype rate[],
214 WebData wdata);
215 static void WebRatesB(realtype x, realtype y, realtype t, realtype c[], realtype cB[],
216 realtype rate[], realtype rateB[], WebData wdata);
217 static void fblock (realtype t, realtype cdata[], int jx, int jy, realtype cdotdata[],
218 WebData wdata);
219 static void GSIter(realtype gamma, N_Vector z, N_Vector x, WebData wdata);
220 static realtype doubleIntgr(N_Vector c, int i, WebData wdata);
221 static int check_retval(void *returnvalue, const char *funcname, int opt);
222
223 /* Small Vector Kernels */
224
225 static void v_inc_by_prod(realtype u[], realtype v[], realtype w[], int n);
226 static void v_sum_prods(realtype u[], realtype p[], realtype q[], realtype v[],
227 realtype w[], int n);
228 static void v_prod(realtype u[], realtype v[], realtype w[], int n);
229 static void v_zero(realtype u[], int n);
230
231 /*
232 *--------------------------------------------------------------------
233 * MAIN PROGRAM
234 *--------------------------------------------------------------------
235 */
236
main(int argc,char * argv[])237 int main(int argc, char *argv[])
238 {
239 realtype abstol=ATOL, reltol=RTOL, t;
240 N_Vector c;
241 WebData wdata;
242 void *cvode_mem;
243 SUNLinearSolver LS, LSB;
244
245 int retval, ncheck;
246
247 int indexB;
248
249 realtype reltolB=RTOL, abstolB=ATOL;
250 N_Vector cB;
251
252 c = cB = NULL;
253 wdata = NULL;
254 cvode_mem = NULL;
255 LS = LSB = NULL;
256
257 /* Allocate and initialize user data */
258
259 wdata = AllocUserData();
260 if(check_retval((void *)wdata, "AllocUserData", 2)) return(1);
261 InitUserData(wdata);
262
263 /* Set-up forward problem */
264
265 /* Initializations */
266 c = N_VNew_Serial(NEQ);
267 if(check_retval((void *)c, "N_VNew_Serial", 0)) return(1);
268 CInit(c, wdata);
269
270 /* Call CVodeCreate/CVodeInit for forward run */
271 printf("\nCreate and allocate CVODES memory for forward run\n");
272 cvode_mem = CVodeCreate(CV_BDF);
273 if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
274 wdata->cvode_mem = cvode_mem; /* Used in Precond */
275 retval = CVodeSetUserData(cvode_mem, wdata);
276 if(check_retval(&retval, "CVodeSetUserData", 1)) return(1);
277 retval = CVodeInit(cvode_mem, f, T0, c);
278 if(check_retval(&retval, "CVodeInit", 1)) return(1);
279 retval = CVodeSStolerances(cvode_mem, reltol, abstol);
280 if(check_retval(&retval, "CVodeSStolerances", 1)) return(1);
281
282 /* Create SUNLinSol_SPGMR linear solver for forward run */
283 LS = SUNLinSol_SPGMR(c, PREC_LEFT, 0);
284 if(check_retval((void *)LS, "SUNLinSol_SPGMR", 0)) return(1);
285
286 /* Attach the linear sovler */
287 retval = CVodeSetLinearSolver(cvode_mem, LS, NULL);
288 if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1;
289
290 /* Set the preconditioner solve and setup functions */
291 retval = CVodeSetPreconditioner(cvode_mem, Precond, PSolve);
292 if(check_retval(&retval, "CVodeSetPreconditioner", 1)) return(1);
293
294 /* Call CVodeSetMaxNumSteps to set the maximum number of steps the
295 * solver will take in an attempt to reach the next output time
296 * during forward integration. */
297 retval = CVodeSetMaxNumSteps(cvode_mem, 2500);
298 if(check_retval(&retval, "CVodeSetMaxNumSteps", 1)) return(1);
299
300 /* Set-up adjoint calculations */
301
302 printf("\nAllocate global memory\n");
303 retval = CVodeAdjInit(cvode_mem, NSTEPS, CV_HERMITE);
304 if(check_retval(&retval, "CVodeAdjInit", 1)) return(1);
305
306 /* Perform forward run */
307
308 printf("\nForward integration\n");
309 retval = CVodeF(cvode_mem, TOUT, c, &t, CV_NORMAL, &ncheck);
310 if(check_retval(&retval, "CVodeF", 1)) return(1);
311
312 printf("\nncheck = %d\n", ncheck);
313
314
315 #if defined(SUNDIALS_EXTENDED_PRECISION)
316 printf("\n g = int_x int_y c%d(Tfinal,x,y) dx dy = %Lf \n\n",
317 ISPEC, doubleIntgr(c,ISPEC,wdata));
318 #else
319 printf("\n g = int_x int_y c%d(Tfinal,x,y) dx dy = %f \n\n",
320 ISPEC, doubleIntgr(c,ISPEC,wdata));
321 #endif
322
323 /* Set-up backward problem */
324
325 /* Allocate cB */
326 cB = N_VNew_Serial(NEQ);
327 if(check_retval((void *)cB, "N_VNew_Serial", 0)) return(1);
328 /* Initialize cB = 0 */
329 CbInit(cB, ISPEC, wdata);
330
331 /* Create and allocate CVODES memory for backward run */
332 printf("\nCreate and allocate CVODES memory for backward run\n");
333 retval = CVodeCreateB(cvode_mem, CV_BDF, &indexB);
334 if(check_retval(&retval, "CVodeCreateB", 1)) return(1);
335 retval = CVodeSetUserDataB(cvode_mem, indexB, wdata);
336 if(check_retval(&retval, "CVodeSetUserDataB", 1)) return(1);
337 retval = CVodeInitB(cvode_mem, indexB, fB, TOUT, cB);
338 if(check_retval(&retval, "CVodeInitB", 1)) return(1);
339 retval = CVodeSStolerancesB(cvode_mem, indexB, reltolB, abstolB);
340 if(check_retval(&retval, "CVodeSStolerancesB", 1)) return(1);
341
342 wdata->indexB = indexB;
343
344 /* Create SUNLinSol_SPGMR linear solver for backward run */
345 LSB = SUNLinSol_SPGMR(cB, PREC_LEFT, 0);
346 if(check_retval((void *)LSB, "SUNLinSol_SPGMR", 0)) return(1);
347
348 /* Attach the linear sovler */
349 retval = CVodeSetLinearSolverB(cvode_mem, indexB, LSB, NULL);
350 if (check_retval(&retval, "CVodeSetLinearSolverB", 1)) return 1;
351
352 /* Set the preconditioner solve and setup functions */
353 retval = CVodeSetPreconditionerB(cvode_mem, indexB, PrecondB, PSolveB);
354 if(check_retval(&retval, "CVodeSetPreconditionerB", 1)) return(1);
355
356 /* Perform backward integration */
357
358 printf("\nBackward integration\n");
359 retval = CVodeB(cvode_mem, T0, CV_NORMAL);
360 if(check_retval(&retval, "CVodeB", 1)) return(1);
361
362 retval = CVodeGetB(cvode_mem, indexB, &t, cB);
363 if(check_retval(&retval, "CVodeGetB", 1)) return(1);
364
365 PrintOutput(cB, NS, MXNS, wdata);
366
367 /* Free all memory */
368 CVodeFree(&cvode_mem);
369
370 N_VDestroy(c);
371 N_VDestroy(cB);
372 SUNLinSolFree(LS);
373 SUNLinSolFree(LSB);
374
375 FreeUserData(wdata);
376
377 return(0);
378 }
379
380 /*
381 *--------------------------------------------------------------------
382 * FUNCTIONS CALLED BY CVODES
383 *--------------------------------------------------------------------
384 */
385
386 /*
387 * This routine computes the right-hand side of the ODE system and
388 * returns it in cdot. The interaction rates are computed by calls to WebRates,
389 * and these are saved in fsave for use in preconditioning.
390 */
391
f(realtype t,N_Vector c,N_Vector cdot,void * user_data)392 static int f(realtype t, N_Vector c, N_Vector cdot, void *user_data)
393 {
394 int i, ic, ici, idxl, idxu, idyl, idyu, iyoff, jx, jy, ns, mxns;
395 realtype dcxli, dcxui, dcyli, dcyui, x, y, *cox, *coy, *fsave, dx, dy;
396 realtype *cdata, *cdotdata;
397 WebData wdata;
398
399 wdata = (WebData) user_data;
400 cdata = N_VGetArrayPointer(c);
401 cdotdata = N_VGetArrayPointer(cdot);
402
403 mxns = wdata->mxns;
404 ns = wdata->ns;
405 fsave = wdata->fsave;
406 cox = wdata->cox;
407 coy = wdata->coy;
408 mxns = wdata->mxns;
409 dx = wdata->dx;
410 dy = wdata->dy;
411
412 for (jy = 0; jy < MY; jy++) {
413 y = jy*dy;
414 iyoff = mxns*jy;
415 idyu = (jy == MY-1) ? -mxns : mxns;
416 idyl = (jy == 0) ? -mxns : mxns;
417 for (jx = 0; jx < MX; jx++) {
418 x = jx*dx;
419 ic = iyoff + ns*jx;
420 /* Get interaction rates at one point (x,y). */
421 WebRates(x, y, t, cdata+ic, fsave+ic, wdata);
422 idxu = (jx == MX-1) ? -ns : ns;
423 idxl = (jx == 0) ? -ns : ns;
424 for (i = 1; i <= ns; i++) {
425 ici = ic + i-1;
426 /* Do differencing in y. */
427 dcyli = cdata[ici] - cdata[ici-idyl];
428 dcyui = cdata[ici+idyu] - cdata[ici];
429 /* Do differencing in x. */
430 dcxli = cdata[ici] - cdata[ici-idxl];
431 dcxui = cdata[ici+idxu] - cdata[ici];
432 /* Collect terms and load cdot elements. */
433 cdotdata[ici] = coy[i-1]*(dcyui - dcyli) +
434 cox[i-1]*(dcxui - dcxli) +
435 fsave[ici];
436 }
437 }
438 }
439
440 return(0);
441 }
442
443 /*
444 * This routine generates the block-diagonal part of the Jacobian
445 * corresponding to the interaction rates, multiplies by -gamma, adds
446 * the identity matrix, and calls denseGETRF to do the LU decomposition of
447 * each diagonal block. The computation of the diagonal blocks uses
448 * the preset block and grouping information. One block per group is
449 * computed. The Jacobian elements are generated by difference
450 * quotients using calls to the routine fblock.
451 *
452 * This routine can be regarded as a prototype for the general case
453 * of a block-diagonal preconditioner. The blocks are of size mp, and
454 * there are ngrp=ngx*ngy blocks computed in the block-grouping scheme.
455 */
456
Precond(realtype t,N_Vector c,N_Vector fc,booleantype jok,booleantype * jcurPtr,realtype gamma,void * user_data)457 static int Precond(realtype t, N_Vector c, N_Vector fc,
458 booleantype jok, booleantype *jcurPtr,
459 realtype gamma, void *user_data)
460 {
461 realtype ***P;
462 sunindextype **pivot;
463 int i, if0, if00, ig, igx, igy, j, jj, jx, jy;
464 int *jxr, *jyr, ngrp, ngx, ngy, mxmp, mp, retval;
465 sunindextype denseretval;
466 realtype uround, fac, r, r0, save, srur;
467 realtype *f1, *fsave, *cdata, *rewtdata;
468 WebData wdata;
469 N_Vector rewt;
470
471 wdata = (WebData) user_data;
472 rewt = wdata->rewt;
473 retval = CVodeGetErrWeights(wdata->cvode_mem, rewt);
474 if(check_retval(&retval, "CVodeGetErrWeights", 1)) return(1);
475
476 cdata = N_VGetArrayPointer(c);
477 rewtdata = N_VGetArrayPointer(rewt);
478
479 uround = UNIT_ROUNDOFF;
480
481 P = wdata->P;
482 pivot = wdata->pivot;
483 jxr = wdata->jxr;
484 jyr = wdata->jyr;
485 mp = wdata->mp;
486 srur = wdata->srur;
487 ngrp = wdata->ngrp;
488 ngx = wdata->ngx;
489 ngy = wdata->ngy;
490 mxmp = wdata->mxmp;
491 fsave = wdata->fsave;
492
493 /* Make mp calls to fblock to approximate each diagonal block of Jacobian.
494 Here, fsave contains the base value of the rate vector and
495 r0 is a minimum increment factor for the difference quotient. */
496
497 f1 = N_VGetArrayPointer(wdata->vtemp);
498
499 fac = N_VWrmsNorm (fc, rewt);
500 r0 = RCONST(1000.0)*fabs(gamma)*uround*NEQ*fac;
501 if (r0 == ZERO) r0 = ONE;
502
503 for (igy = 0; igy < ngy; igy++) {
504 jy = jyr[igy];
505 if00 = jy*mxmp;
506 for (igx = 0; igx < ngx; igx++) {
507 jx = jxr[igx];
508 if0 = if00 + jx*mp;
509 ig = igx + igy*ngx;
510 /* Generate ig-th diagonal block */
511 for (j = 0; j < mp; j++) {
512 /* Generate the jth column as a difference quotient */
513 jj = if0 + j;
514 save = cdata[jj];
515 r = MAX(srur*fabs(save),r0/rewtdata[jj]);
516 cdata[jj] += r;
517 fac = -gamma/r;
518 fblock (t, cdata, jx, jy, f1, wdata);
519 for (i = 0; i < mp; i++) {
520 P[ig][j][i] = (f1[i] - fsave[if0+i])*fac;
521 }
522 cdata[jj] = save;
523 }
524 }
525 }
526
527 /* Add identity matrix and do LU decompositions on blocks. */
528
529 for (ig = 0; ig < ngrp; ig++) {
530 denseAddIdentity(P[ig], mp);
531 denseretval = denseGETRF(P[ig], mp, mp, pivot[ig]);
532 if (denseretval != 0) return(1);
533 }
534
535 *jcurPtr = SUNTRUE;
536 return(0);
537 }
538
539 /*
540 * This routine applies two inverse preconditioner matrices
541 * to the vector r, using the interaction-only block-diagonal Jacobian
542 * with block-grouping, denoted Jr, and Gauss-Seidel applied to the
543 * diffusion contribution to the Jacobian, denoted Jd.
544 * It first calls GSIter for a Gauss-Seidel approximation to
545 * ((I - gamma*Jd)-inverse)*r, and stores the result in z.
546 * Then it computes ((I - gamma*Jr)-inverse)*z, using LU factors of the
547 * blocks in P, and pivot information in pivot, and returns the result in z.
548 */
549
PSolve(realtype t,N_Vector c,N_Vector fc,N_Vector r,N_Vector z,realtype gamma,realtype delta,int lr,void * user_data)550 static int PSolve(realtype t, N_Vector c, N_Vector fc,
551 N_Vector r, N_Vector z,
552 realtype gamma, realtype delta,
553 int lr, void *user_data)
554 {
555 realtype ***P;
556 sunindextype **pivot;
557 int jx, jy, igx, igy, iv, ig, *jigx, *jigy, mx, my, ngx, mp;
558 WebData wdata;
559
560 wdata = (WebData) user_data;
561
562 N_VScale(ONE, r, z);
563
564 /* call GSIter for Gauss-Seidel iterations */
565
566 GSIter(gamma, z, wdata->vtemp, wdata);
567
568 /* Do backsolves for inverse of block-diagonal preconditioner factor */
569
570 P = wdata->P;
571 pivot = wdata->pivot;
572 mx = wdata->mx;
573 my = wdata->my;
574 ngx = wdata->ngx;
575 mp = wdata->mp;
576 jigx = wdata->jigx;
577 jigy = wdata->jigy;
578
579 iv = 0;
580 for (jy = 0; jy < my; jy++) {
581 igy = jigy[jy];
582 for (jx = 0; jx < mx; jx++) {
583 igx = jigx[jx];
584 ig = igx + igy*ngx;
585 denseGETRS(P[ig], mp, pivot[ig], &(N_VGetArrayPointer(z)[iv]));
586 iv += mp;
587 }
588 }
589
590 return(0);
591 }
592
593 /*
594 * This routine computes the right-hand side of the adjoint ODE system and
595 * returns it in cBdot. The interaction rates are computed by calls to WebRates,
596 * and these are saved in fsave for use in preconditioning. The adjoint
597 * interaction rates are computed by calls to WebRatesB.
598 */
599
fB(realtype t,N_Vector c,N_Vector cB,N_Vector cBdot,void * user_data)600 static int fB(realtype t, N_Vector c, N_Vector cB,
601 N_Vector cBdot, void *user_data)
602 {
603 int i, ic, ici, idxl, idxu, idyl, idyu, iyoff, jx, jy, ns, mxns;
604 realtype dcxli, dcxui, dcyli, dcyui, x, y, *cox, *coy, *fsave, *fBsave, dx, dy;
605 realtype *cdata, *cBdata, *cBdotdata;
606 WebData wdata;
607
608 wdata = (WebData) user_data;
609 cdata = N_VGetArrayPointer(c);
610 cBdata = N_VGetArrayPointer(cB);
611 cBdotdata = N_VGetArrayPointer(cBdot);
612
613 mxns = wdata->mxns;
614 ns = wdata->ns;
615 fsave = wdata->fsave;
616 fBsave = wdata->fBsave;
617 cox = wdata->cox;
618 coy = wdata->coy;
619 mxns = wdata->mxns;
620 dx = wdata->dx;
621 dy = wdata->dy;
622
623 for (jy = 0; jy < MY; jy++) {
624 y = jy*dy;
625 iyoff = mxns*jy;
626 idyu = (jy == MY-1) ? -mxns : mxns;
627 idyl = (jy == 0) ? -mxns : mxns;
628 for (jx = 0; jx < MX; jx++) {
629 x = jx*dx;
630 ic = iyoff + ns*jx;
631 /* Get interaction rates at one point (x,y). */
632 WebRatesB(x, y, t, cdata+ic, cBdata+ic, fsave+ic, fBsave+ic, wdata);
633 idxu = (jx == MX-1) ? -ns : ns;
634 idxl = (jx == 0) ? -ns : ns;
635 for (i = 1; i <= ns; i++) {
636 ici = ic + i-1;
637 /* Do differencing in y. */
638 dcyli = cBdata[ici] - cBdata[ici-idyl];
639 dcyui = cBdata[ici+idyu] - cBdata[ici];
640 /* Do differencing in x. */
641 dcxli = cBdata[ici] - cBdata[ici-idxl];
642 dcxui = cBdata[ici+idxu] - cBdata[ici];
643 /* Collect terms and load cdot elements. */
644 cBdotdata[ici] = - coy[i-1]*(dcyui - dcyli)
645 - cox[i-1]*(dcxui - dcxli)
646 - fBsave[ici];
647 }
648 }
649 }
650
651 return(0);
652 }
653
654 /*
655 * Preconditioner setup function for the backward problem
656 */
657
PrecondB(realtype t,N_Vector c,N_Vector cB,N_Vector fcB,booleantype jok,booleantype * jcurPtr,realtype gamma,void * user_data)658 static int PrecondB(realtype t, N_Vector c,
659 N_Vector cB, N_Vector fcB, booleantype jok,
660 booleantype *jcurPtr, realtype gamma,
661 void *user_data)
662 {
663 realtype ***P;
664 sunindextype **pivot;
665 sunindextype denseretval;
666 int i, if0, if00, ig, igx, igy, j, jj, jx, jy;
667 int *jxr, *jyr, mp, ngrp, ngx, ngy, mxmp, retval;
668 realtype uround, fac, r, r0, save, srur;
669 realtype *f1, *fsave, *cdata, *rewtdata;
670 void *cvode_mem;
671 WebData wdata;
672 N_Vector rewt;
673
674 wdata = (WebData) user_data;
675 cvode_mem = CVodeGetAdjCVodeBmem(wdata->cvode_mem, wdata->indexB);
676 if(check_retval((void *)cvode_mem, "CVadjGetCVodeBmem", 0)) return(1);
677 rewt = wdata->rewt;
678 retval = CVodeGetErrWeights(cvode_mem, rewt);
679 if(check_retval(&retval, "CVodeGetErrWeights", 1)) return(1);
680
681 cdata = N_VGetArrayPointer(c);
682 rewtdata = N_VGetArrayPointer(rewt);
683
684 uround = UNIT_ROUNDOFF;
685
686 P = wdata->P;
687 pivot = wdata->pivot;
688 jxr = wdata->jxr;
689 jyr = wdata->jyr;
690 mp = wdata->mp;
691 srur = wdata->srur;
692 ngrp = wdata->ngrp;
693 ngx = wdata->ngx;
694 ngy = wdata->ngy;
695 mxmp = wdata->mxmp;
696 fsave = wdata->fsave;
697
698 /* Make mp calls to fblock to approximate each diagonal block of Jacobian.
699 Here, fsave contains the base value of the rate vector and
700 r0 is a minimum increment factor for the difference quotient. */
701
702 f1 = N_VGetArrayPointer(wdata->vtemp);
703
704 fac = N_VWrmsNorm (fcB, rewt);
705 r0 = RCONST(1000.0)*fabs(gamma)*uround*NEQ*fac;
706 if (r0 == ZERO) r0 = ONE;
707
708 for (igy = 0; igy < ngy; igy++) {
709 jy = jyr[igy];
710 if00 = jy*mxmp;
711 for (igx = 0; igx < ngx; igx++) {
712 jx = jxr[igx];
713 if0 = if00 + jx*mp;
714 ig = igx + igy*ngx;
715 /* Generate ig-th diagonal block */
716 for (j = 0; j < mp; j++) {
717 /* Generate the jth column as a difference quotient */
718 jj = if0 + j;
719 save = cdata[jj];
720 r = MAX(srur*fabs(save),r0/rewtdata[jj]);
721 cdata[jj] += r;
722 fac = gamma/r;
723 fblock (t, cdata, jx, jy, f1, wdata);
724 for (i = 0; i < mp; i++) {
725 P[ig][i][j] = (f1[i] - fsave[if0+i])*fac;
726 }
727 cdata[jj] = save;
728 }
729 }
730 }
731
732 /* Add identity matrix and do LU decompositions on blocks. */
733
734 for (ig = 0; ig < ngrp; ig++) {
735 denseAddIdentity(P[ig], mp);
736 denseretval = denseGETRF(P[ig], mp, mp, pivot[ig]);
737 if (denseretval != 0) return(1);
738 }
739
740 *jcurPtr = SUNTRUE;
741 return(0);
742 }
743
744 /*
745 * Preconditioner solve function for the backward problem
746 */
747
PSolveB(realtype t,N_Vector c,N_Vector cB,N_Vector fcB,N_Vector r,N_Vector z,realtype gamma,realtype delta,int lr,void * user_data)748 static int PSolveB(realtype t, N_Vector c,
749 N_Vector cB, N_Vector fcB,
750 N_Vector r, N_Vector z,
751 realtype gamma, realtype delta,
752 int lr, void *user_data)
753 {
754 realtype ***P;
755 sunindextype **pivot;
756 int jx, jy, igx, igy, iv, ig, *jigx, *jigy, mx, my, ngx, mp;
757 WebData wdata;
758
759 wdata = (WebData) user_data;
760
761 N_VScale(ONE, r, z);
762
763 /* call GSIter for Gauss-Seidel iterations (same routine but with gamma=-gamma) */
764
765 GSIter(-gamma, z, wdata->vtemp, wdata);
766
767 /* Do backsolves for inverse of block-diagonal preconditioner factor */
768
769 P = wdata->P;
770 pivot = wdata->pivot;
771 mx = wdata->mx;
772 my = wdata->my;
773 ngx = wdata->ngx;
774 mp = wdata->mp;
775 jigx = wdata->jigx;
776 jigy = wdata->jigy;
777
778 iv = 0;
779 for (jy = 0; jy < my; jy++) {
780 igy = jigy[jy];
781 for (jx = 0; jx < mx; jx++) {
782 igx = jigx[jx];
783 ig = igx + igy*ngx;
784 denseGETRS(P[ig], mp, pivot[ig], &(N_VGetArrayPointer(z)[iv]));
785 iv += mp;
786 }
787 }
788
789 return(0);
790 }
791
792 /*
793 *--------------------------------------------------------------------
794 * PRIVATE FUNCTIONS
795 *--------------------------------------------------------------------
796 */
797
798 /*
799 * Allocate space for user data structure
800 */
801
AllocUserData(void)802 static WebData AllocUserData(void)
803 {
804 int i, ngrp = NGRP;
805 sunindextype ns = NS;
806 WebData wdata;
807
808 wdata = (WebData) malloc(sizeof *wdata);
809 for(i=0; i < ngrp; i++) {
810 (wdata->P)[i] = newDenseMat(ns, ns);
811 (wdata->pivot)[i] = newIndexArray(ns);
812 }
813 wdata->rewt = N_VNew_Serial(NEQ);
814 wdata->vtemp = N_VNew_Serial(NEQ);
815
816 return(wdata);
817 }
818
819 /*
820 * Initialize user data structure
821 */
822
InitUserData(WebData wdata)823 static void InitUserData(WebData wdata)
824 {
825 int i, j, ns;
826 realtype *bcoef, *diff, *cox, *coy, dx, dy;
827 realtype (*acoef)[NS];
828
829 acoef = wdata->acoef;
830 bcoef = wdata->bcoef;
831 diff = wdata->diff;
832 cox = wdata->cox;
833 coy = wdata->coy;
834 ns = wdata->ns = NS;
835
836 for (j = 0; j < NS; j++) { for (i = 0; i < NS; i++) acoef[i][j] = ZERO; }
837 for (j = 0; j < NP; j++) {
838 for (i = 0; i < NP; i++) {
839 acoef[NP+i][j] = EE;
840 acoef[i][NP+j] = -GG;
841 }
842 acoef[j][j] = -AA;
843 acoef[NP+j][NP+j] = -AA;
844 bcoef[j] = BB;
845 bcoef[NP+j] = -BB;
846 diff[j] = DPREY;
847 diff[NP+j] = DPRED;
848 }
849
850 /* Set remaining problem parameters */
851
852 wdata->mxns = MXNS;
853 dx = wdata->dx = DX;
854 dy = wdata->dy = DY;
855 for (i = 0; i < ns; i++) {
856 cox[i] = diff[i]/SQR(dx);
857 coy[i] = diff[i]/SQR(dy);
858 }
859
860 /* Set remaining method parameters */
861
862 wdata->mp = MP;
863 wdata->mq = MQ;
864 wdata->mx = MX;
865 wdata->my = MY;
866 wdata->srur = sqrt(UNIT_ROUNDOFF);
867 wdata->mxmp = MXMP;
868 wdata->ngrp = NGRP;
869 wdata->ngx = NGX;
870 wdata->ngy = NGY;
871 SetGroups(MX, NGX, wdata->jgx, wdata->jigx, wdata->jxr);
872 SetGroups(MY, NGY, wdata->jgy, wdata->jigy, wdata->jyr);
873 }
874
875 /*
876 * This routine sets arrays jg, jig, and jr describing
877 * a uniform partition of (0,1,2,...,m-1) into ng groups.
878 * The arrays set are:
879 * jg = length ng+1 array of group boundaries.
880 * Group ig has indices j = jg[ig],...,jg[ig+1]-1.
881 * jig = length m array of group indices vs node index.
882 * Node index j is in group jig[j].
883 * jr = length ng array of indices representing the groups.
884 * The index for group ig is j = jr[ig].
885 */
886
SetGroups(int m,int ng,int jg[],int jig[],int jr[])887 static void SetGroups(int m, int ng, int jg[], int jig[], int jr[])
888 {
889 int ig, j, len1, mper, ngm1;
890
891 mper = m/ng; /* does integer division */
892 for (ig=0; ig < ng; ig++) jg[ig] = ig*mper;
893 jg[ng] = m;
894
895 ngm1 = ng - 1;
896 len1 = ngm1*mper;
897 for (j = 0; j < len1; j++) jig[j] = j/mper;
898 for (j = len1; j < m; j++) jig[j] = ngm1;
899
900 for (ig = 0; ig < ngm1; ig++) jr[ig] = ((2*ig+1)*mper-1)/2;
901 jr[ngm1] = (ngm1*mper+m-1)/2;
902 }
903
904 /*
905 * This routine computes and loads the vector of initial values.
906 */
907
CInit(N_Vector c,WebData wdata)908 static void CInit(N_Vector c, WebData wdata)
909 {
910 int i, ici, ioff, iyoff, jx, jy, ns, mxns;
911 realtype argx, argy, x, y, dx, dy, x_factor, y_factor, *cdata;
912
913 cdata = N_VGetArrayPointer(c);
914 ns = wdata->ns;
915 mxns = wdata->mxns;
916 dx = wdata->dx;
917 dy = wdata->dy;
918
919 x_factor = RCONST(4.0)/SQR(AX);
920 y_factor = RCONST(4.0)/SQR(AY);
921 for (jy = 0; jy < MY; jy++) {
922 y = jy*dy;
923 argy = SQR(y_factor*y*(AY-y));
924 iyoff = mxns*jy;
925 for (jx = 0; jx < MX; jx++) {
926 x = jx*dx;
927 argx = SQR(x_factor*x*(AX-x));
928 ioff = iyoff + ns*jx;
929 for (i = 1; i <= ns; i++) {
930 ici = ioff + i-1;
931 cdata[ici] = RCONST(10.0) + i*argx*argy;
932 }
933 }
934 }
935 }
936
937 /*
938 * This function computes and loads the final values for the adjoint variables
939 */
940
CbInit(N_Vector c,int is,WebData wdata)941 static void CbInit(N_Vector c, int is, WebData wdata)
942 {
943 int i, ici, ioff, iyoff, jx, jy, ns, mxns;
944 realtype *cdata;
945
946 realtype gu[NS];
947
948 cdata = N_VGetArrayPointer(c);
949 ns = wdata->ns;
950 mxns = wdata->mxns;
951
952 for ( i = 1; i <= ns; i++ ) gu[i-1] = ZERO;
953 gu[ISPEC-1] = ONE;
954
955 for (jy = 0; jy < MY; jy++) {
956 iyoff = mxns*jy;
957 for (jx = 0; jx < MX; jx++) {
958 ioff = iyoff + ns*jx;
959 for (i = 1; i <= ns; i++) {
960 ici = ioff + i-1;
961 cdata[ici] = gu[i-1];
962 }
963 }
964 }
965 }
966
967 /*
968 * This routine computes the interaction rates for the species
969 * c_1, ... ,c_ns (stored in c[0],...,c[ns-1]), at one spatial point
970 * and at time t.
971 */
972
WebRates(realtype x,realtype y,realtype t,realtype c[],realtype rate[],WebData wdata)973 static void WebRates(realtype x, realtype y, realtype t, realtype c[],
974 realtype rate[], WebData wdata)
975 {
976 int i, j, ns;
977 realtype fac, *bcoef;
978 realtype (*acoef)[NS];
979
980 ns = wdata->ns;
981 acoef = wdata->acoef;
982 bcoef = wdata->bcoef;
983
984 for (i = 0; i < ns; i++)
985 rate[i] = ZERO;
986
987 for (j = 0; j < ns; j++)
988 for (i = 0; i < ns; i++)
989 rate[i] += c[j] * acoef[i][j];
990
991 fac = ONE + ALPH*x*y;
992 for (i = 0; i < ns; i++)
993 rate[i] = c[i]*(bcoef[i]*fac + rate[i]);
994 }
995
996 /*
997 * This routine computes the interaction rates for the backward problem
998 */
999
WebRatesB(realtype x,realtype y,realtype t,realtype c[],realtype cB[],realtype rate[],realtype rateB[],WebData wdata)1000 static void WebRatesB(realtype x, realtype y, realtype t, realtype c[], realtype cB[],
1001 realtype rate[], realtype rateB[], WebData wdata)
1002 {
1003 int i, j, ns;
1004 realtype fac, *bcoef;
1005 realtype (*acoef)[NS];
1006
1007 ns = wdata->ns;
1008 acoef = wdata->acoef;
1009 bcoef = wdata->bcoef;
1010
1011 fac = ONE + ALPH*x*y;
1012
1013 for (i = 0; i < ns; i++)
1014 rate[i] = bcoef[i]*fac;
1015
1016 for (j = 0; j < ns; j++)
1017 for (i = 0; i < ns; i++)
1018 rate[i] += acoef[i][j]*c[j];
1019
1020 for (i = 0; i < ns; i++) {
1021 rateB[i] = cB[i]*rate[i];
1022 rate[i] = c[i]*rate[i];
1023 }
1024
1025 for (j = 0; j < ns; j++)
1026 for (i = 0; i < ns; i++)
1027 rateB[i] += acoef[j][i]*c[j]*cB[j];
1028
1029 }
1030
1031 /*
1032 * This routine computes one block of the interaction terms of the
1033 * system, namely block (jx,jy), for use in preconditioning.
1034 * Here jx and jy count from 0.
1035 */
1036
fblock(realtype t,realtype cdata[],int jx,int jy,realtype cdotdata[],WebData wdata)1037 static void fblock(realtype t, realtype cdata[], int jx, int jy,
1038 realtype cdotdata[], WebData wdata)
1039 {
1040 int iblok, ic;
1041 realtype x, y;
1042
1043 iblok = jx + jy*(wdata->mx);
1044 y = jy*(wdata->dy);
1045 x = jx*(wdata->dx);
1046 ic = (wdata->ns)*(iblok);
1047 WebRates(x, y, t, cdata+ic, cdotdata, wdata);
1048 }
1049
1050 /*
1051 * This routine performs ITMAX=5 Gauss-Seidel iterations to compute an
1052 * approximation to (P-inverse)*z, where P = I - gamma*Jd, and
1053 * Jd represents the diffusion contributions to the Jacobian.
1054 * The answer is stored in z on return, and x is a temporary vector.
1055 * The dimensions below assume a global constant NS >= ns.
1056 * Some inner loops of length ns are implemented with the small
1057 * vector kernels v_sum_prods, v_prod, v_inc_by_prod.
1058 */
1059
GSIter(realtype gamma,N_Vector z,N_Vector x,WebData wdata)1060 static void GSIter(realtype gamma, N_Vector z, N_Vector x, WebData wdata)
1061 {
1062 int i, ic, iter, iyoff, jx, jy, ns, mxns, mx, my, x_loc, y_loc;
1063 realtype beta[NS], beta2[NS], cof1[NS], gam[NS], gam2[NS];
1064 realtype temp, *cox, *coy, *xd, *zd;
1065
1066 xd = N_VGetArrayPointer(x);
1067 zd = N_VGetArrayPointer(z);
1068 ns = wdata->ns;
1069 mx = wdata->mx;
1070 my = wdata->my;
1071 mxns = wdata->mxns;
1072 cox = wdata->cox;
1073 coy = wdata->coy;
1074
1075 /* Write matrix as P = D - L - U.
1076 Load local arrays beta, beta2, gam, gam2, and cof1. */
1077
1078 for (i = 0; i < ns; i++) {
1079 temp = ONE/(ONE + TWO*gamma*(cox[i] + coy[i]));
1080 beta[i] = gamma*cox[i]*temp;
1081 beta2[i] = TWO*beta[i];
1082 gam[i] = gamma*coy[i]*temp;
1083 gam2[i] = TWO*gam[i];
1084 cof1[i] = temp;
1085 }
1086
1087 /* Begin iteration loop.
1088 Load vector x with (D-inverse)*z for first iteration. */
1089
1090 for (jy = 0; jy < my; jy++) {
1091 iyoff = mxns*jy;
1092 for (jx = 0; jx < mx; jx++) {
1093 ic = iyoff + ns*jx;
1094 v_prod(xd+ic, cof1, zd+ic, ns); /* x[ic+i] = cof1[i]z[ic+i] */
1095 }
1096 }
1097 N_VConst(ZERO, z);
1098
1099 /* Looping point for iterations. */
1100
1101 for (iter=1; iter <= ITMAX; iter++) {
1102
1103 /* Calculate (D-inverse)*U*x if not the first iteration. */
1104
1105 if (iter > 1) {
1106 for (jy=0; jy < my; jy++) {
1107 iyoff = mxns*jy;
1108 for (jx=0; jx < mx; jx++) { /* order of loops matters */
1109 ic = iyoff + ns*jx;
1110 x_loc = (jx == 0) ? 0 : ((jx == mx-1) ? 2 : 1);
1111 y_loc = (jy == 0) ? 0 : ((jy == my-1) ? 2 : 1);
1112 switch (3*y_loc+x_loc) {
1113 case 0 : /* jx == 0, jy == 0 */
1114 /* x[ic+i] = beta2[i]x[ic+ns+i] + gam2[i]x[ic+mxns+i] */
1115 v_sum_prods(xd+ic, beta2, xd+ic+ns, gam2, xd+ic+mxns, ns);
1116 break;
1117 case 1 : /* 1 <= jx <= mx-2, jy == 0 */
1118 /* x[ic+i] = beta[i]x[ic+ns+i] + gam2[i]x[ic+mxns+i] */
1119 v_sum_prods(xd+ic, beta, xd+ic+ns, gam2, xd+ic+mxns, ns);
1120 break;
1121 case 2 : /* jx == mx-1, jy == 0 */
1122 /* x[ic+i] = gam2[i]x[ic+mxns+i] */
1123 v_prod(xd+ic, gam2, xd+ic+mxns, ns);
1124 break;
1125 case 3 : /* jx == 0, 1 <= jy <= my-2 */
1126 /* x[ic+i] = beta2[i]x[ic+ns+i] + gam[i]x[ic+mxns+i] */
1127 v_sum_prods(xd+ic, beta2, xd+ic+ns, gam, xd+ic+mxns, ns);
1128 break;
1129 case 4 : /* 1 <= jx <= mx-2, 1 <= jy <= my-2 */
1130 /* x[ic+i] = beta[i]x[ic+ns+i] + gam[i]x[ic+mxns+i] */
1131 v_sum_prods(xd+ic, beta, xd+ic+ns, gam, xd+ic+mxns, ns);
1132 break;
1133 case 5 : /* jx == mx-1, 1 <= jy <= my-2 */
1134 /* x[ic+i] = gam[i]x[ic+mxns+i] */
1135 v_prod(xd+ic, gam, xd+ic+mxns, ns);
1136 break;
1137 case 6 : /* jx == 0, jy == my-1 */
1138 /* x[ic+i] = beta2[i]x[ic+ns+i] */
1139 v_prod(xd+ic, beta2, xd+ic+ns, ns);
1140 break;
1141 case 7 : /* 1 <= jx <= mx-2, jy == my-1 */
1142 /* x[ic+i] = beta[i]x[ic+ns+i] */
1143 v_prod(xd+ic, beta, xd+ic+ns, ns);
1144 break;
1145 case 8 : /* jx == mx-1, jy == my-1 */
1146 /* x[ic+i] = ZERO */
1147 v_zero(xd+ic, ns);
1148 break;
1149 }
1150 }
1151 }
1152 } /* end if (iter > 1) */
1153
1154 /* Overwrite x with [(I - (D-inverse)*L)-inverse]*x. */
1155
1156 for (jy=0; jy < my; jy++) {
1157 iyoff = mxns*jy;
1158 for (jx=0; jx < mx; jx++) { /* order of loops matters */
1159 ic = iyoff + ns*jx;
1160 x_loc = (jx == 0) ? 0 : ((jx == mx-1) ? 2 : 1);
1161 y_loc = (jy == 0) ? 0 : ((jy == my-1) ? 2 : 1);
1162 switch (3*y_loc+x_loc) {
1163 case 0 : /* jx == 0, jy == 0 */
1164 break;
1165 case 1 : /* 1 <= jx <= mx-2, jy == 0 */
1166 /* x[ic+i] += beta[i]x[ic-ns+i] */
1167 v_inc_by_prod(xd+ic, beta, xd+ic-ns, ns);
1168 break;
1169 case 2 : /* jx == mx-1, jy == 0 */
1170 /* x[ic+i] += beta2[i]x[ic-ns+i] */
1171 v_inc_by_prod(xd+ic, beta2, xd+ic-ns, ns);
1172 break;
1173 case 3 : /* jx == 0, 1 <= jy <= my-2 */
1174 /* x[ic+i] += gam[i]x[ic-mxns+i] */
1175 v_inc_by_prod(xd+ic, gam, xd+ic-mxns, ns);
1176 break;
1177 case 4 : /* 1 <= jx <= mx-2, 1 <= jy <= my-2 */
1178 /* x[ic+i] += beta[i]x[ic-ns+i] + gam[i]x[ic-mxns+i] */
1179 v_inc_by_prod(xd+ic, beta, xd+ic-ns, ns);
1180 v_inc_by_prod(xd+ic, gam, xd+ic-mxns, ns);
1181 break;
1182 case 5 : /* jx == mx-1, 1 <= jy <= my-2 */
1183 /* x[ic+i] += beta2[i]x[ic-ns+i] + gam[i]x[ic-mxns+i] */
1184 v_inc_by_prod(xd+ic, beta2, xd+ic-ns, ns);
1185 v_inc_by_prod(xd+ic, gam, xd+ic-mxns, ns);
1186 break;
1187 case 6 : /* jx == 0, jy == my-1 */
1188 /* x[ic+i] += gam2[i]x[ic-mxns+i] */
1189 v_inc_by_prod(xd+ic, gam2, xd+ic-mxns, ns);
1190 break;
1191 case 7 : /* 1 <= jx <= mx-2, jy == my-1 */
1192 /* x[ic+i] += beta[i]x[ic-ns+i] + gam2[i]x[ic-mxns+i] */
1193 v_inc_by_prod(xd+ic, beta, xd+ic-ns, ns);
1194 v_inc_by_prod(xd+ic, gam2, xd+ic-mxns, ns);
1195 break;
1196 case 8 : /* jx == mx-1, jy == my-1 */
1197 /* x[ic+i] += beta2[i]x[ic-ns+i] + gam2[i]x[ic-mxns+i] */
1198 v_inc_by_prod(xd+ic, beta2, xd+ic-ns, ns);
1199 v_inc_by_prod(xd+ic, gam2, xd+ic-mxns, ns);
1200 break;
1201 }
1202 }
1203 }
1204
1205 /* Add increment x to z : z <- z+x */
1206
1207 N_VLinearSum(ONE, z, ONE, x, z);
1208
1209 }
1210 }
1211
v_inc_by_prod(realtype u[],realtype v[],realtype w[],int n)1212 static void v_inc_by_prod(realtype u[], realtype v[], realtype w[], int n)
1213 {
1214 int i;
1215 for (i=0; i < n; i++) u[i] += v[i]*w[i];
1216 }
1217
v_sum_prods(realtype u[],realtype p[],realtype q[],realtype v[],realtype w[],int n)1218 static void v_sum_prods(realtype u[], realtype p[], realtype q[],
1219 realtype v[], realtype w[], int n)
1220 {
1221 int i;
1222 for (i=0; i < n; i++) u[i] = p[i]*q[i] + v[i]*w[i];
1223 }
1224
v_prod(realtype u[],realtype v[],realtype w[],int n)1225 static void v_prod(realtype u[], realtype v[], realtype w[], int n)
1226 {
1227 int i;
1228 for (i=0; i < n; i++) u[i] = v[i]*w[i];
1229 }
1230
v_zero(realtype u[],int n)1231 static void v_zero(realtype u[], int n)
1232 {
1233 int i;
1234 for (i=0; i < n; i++) u[i] = ZERO;
1235 }
1236
1237 /*
1238 * Print maximum sensitivity of G for each species
1239 */
1240
PrintOutput(N_Vector cB,int ns,int mxns,WebData wdata)1241 static void PrintOutput(N_Vector cB, int ns, int mxns, WebData wdata)
1242 {
1243 int i, jx, jy;
1244 realtype *cdata, cij, cmax, x, y;
1245
1246 x = y = ZERO;
1247
1248 cdata = N_VGetArrayPointer(cB);
1249
1250 for (i=1; i <= ns; i++) {
1251
1252 cmax = ZERO;
1253
1254 for (jy=MY-1; jy >= 0; jy--) {
1255 for (jx=0; jx < MX; jx++) {
1256 cij = cdata[(i-1) + jx*ns + jy*mxns];
1257 if (fabs(cij) > cmax) {
1258 cmax = cij;
1259 x = jx * wdata->dx;
1260 y = jy * wdata->dy;
1261 }
1262 }
1263 }
1264
1265 printf("\nMaximum sensitivity with respect to I.C. of species %d\n", i);
1266 #if defined(SUNDIALS_EXTENDED_PRECISION)
1267 printf(" mu max = %Le\n",cmax);
1268 #elif defined(SUNDIALS_DOUBLE_PRECISION)
1269 printf(" mu max = %e\n",cmax);
1270 #else
1271 printf(" mu max = %e\n",cmax);
1272 #endif
1273 printf("at\n");
1274 #if defined(SUNDIALS_EXTENDED_PRECISION)
1275 printf(" x = %Le\n y = %Le\n", x, y);
1276 #elif defined(SUNDIALS_DOUBLE_PRECISION)
1277 printf(" x = %e\n y = %e\n", x, y);
1278 #else
1279 printf(" x = %e\n y = %e\n", x, y);
1280 #endif
1281
1282 }
1283
1284 }
1285
1286 /*
1287 * Compute double space integral
1288 */
1289
doubleIntgr(N_Vector c,int i,WebData wdata)1290 static realtype doubleIntgr(N_Vector c, int i, WebData wdata)
1291 {
1292 realtype *cdata;
1293 int ns, mx, my, mxns;
1294 realtype dx, dy;
1295 realtype intgr_xy, intgr_x;
1296 int jx, jy;
1297
1298 cdata = N_VGetArrayPointer(c);
1299
1300 ns = wdata->ns;
1301 mx = wdata->mx;
1302 my = wdata->my;
1303 mxns = wdata->mxns;
1304 dx = wdata->dx;
1305 dy = wdata->dy;
1306
1307 jy = 0;
1308 intgr_x = cdata[(i-1)+jy*mxns];
1309 for (jx = 1; jx < mx-1; jx++) {
1310 intgr_x += TWO*cdata[(i-1) + jx*ns + jy*mxns];
1311 }
1312 intgr_x += cdata[(i-1)+(mx-1)*ns+jy*mxns];
1313 intgr_x *= RCONST(0.5)*dx;
1314
1315 intgr_xy = intgr_x;
1316
1317 for (jy = 1; jy < my-1; jy++) {
1318
1319 intgr_x = cdata[(i-1)+jy*mxns];
1320 for (jx = 1; jx < mx-1; jx++) {
1321 intgr_x += TWO*cdata[(i-1) + jx*ns + jy*mxns];
1322 }
1323 intgr_x += cdata[(i-1)+(mx-1)*ns+jy*mxns];
1324 intgr_x *= RCONST(0.5)*dx;
1325
1326 intgr_xy += TWO*intgr_x;
1327
1328 }
1329
1330 jy = my-1;
1331 intgr_x = cdata[(i-1)+jy*mxns];
1332 for (jx = 1; jx < mx-1; jx++) {
1333 intgr_x += TWO*cdata[(i-1) + jx*ns + jy*mxns];
1334 }
1335 intgr_x += cdata[(i-1)+(mx-1)*ns+jy*mxns];
1336 intgr_x *= RCONST(0.5)*dx;
1337
1338 intgr_xy += intgr_x;
1339
1340 intgr_xy *= RCONST(0.5)*dy;
1341
1342 return(intgr_xy);
1343 }
1344
1345 /*
1346 * Free space allocated for the user data structure
1347 */
1348
FreeUserData(WebData wdata)1349 static void FreeUserData(WebData wdata)
1350 {
1351 int i, ngrp;
1352
1353 ngrp = wdata->ngrp;
1354 for(i=0; i < ngrp; i++) {
1355 destroyMat((wdata->P)[i]);
1356 destroyArray((wdata->pivot)[i]);
1357 }
1358 N_VDestroy(wdata->rewt);
1359 N_VDestroy(wdata->vtemp);
1360 free(wdata);
1361 }
1362
1363 /*
1364 * Check function return value.
1365 * opt == 0 means SUNDIALS function allocates memory so check if
1366 * returned NULL pointer
1367 * opt == 1 means SUNDIALS function returns an integer value so check if
1368 * retval < 0
1369 * opt == 2 means function allocates memory so check if returned
1370 * NULL pointer
1371 */
1372
check_retval(void * returnvalue,const char * funcname,int opt)1373 static int check_retval(void *returnvalue, const char *funcname, int opt)
1374 {
1375 int *retval;
1376
1377 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
1378 if (opt == 0 && returnvalue == NULL) {
1379 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
1380 funcname);
1381 return(1); }
1382
1383 /* Check if retval < 0 */
1384 else if (opt == 1) {
1385 retval = (int *) returnvalue;
1386 if (*retval < 0) {
1387 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
1388 funcname, *retval);
1389 return(1); }}
1390
1391 /* Check if function returned NULL pointer - no memory allocated */
1392 else if (opt == 2 && returnvalue == NULL) {
1393 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
1394 funcname);
1395 return(1); }
1396
1397 return(0);
1398 }
1399