1 /* -----------------------------------------------------------------
2 * Programmer(s): Allan Taylor, Alan Hindmarsh and
3 * Radu Serban @ 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 * Example program for IDA: Food web problem.
16 *
17 * This example program (serial version) uses the banded linear
18 * solver, and IDACalcIC for initial condition calculation.
19 *
20 * The mathematical problem solved in this example is a DAE system
21 * that arises from a system of partial differential equations after
22 * spatial discretization. The PDE system is a food web population
23 * model, with predator-prey interaction and diffusion on the unit
24 * square in two dimensions. The dependent variable vector is:
25 *
26 * 1 2 ns
27 * c = (c , c , ..., c ) , ns = 2 * np
28 *
29 * and the PDE's are as follows:
30 *
31 * i i i
32 * dc /dt = d(i)*(c + c ) + R (x,y,c) (i = 1,...,np)
33 * xx yy i
34 *
35 * i i
36 * 0 = d(i)*(c + c ) + R (x,y,c) (i = np+1,...,ns)
37 * xx yy i
38 *
39 * where the reaction terms R are:
40 *
41 * i ns j
42 * R (x,y,c) = c * (b(i) + sum a(i,j)*c )
43 * i j=1
44 *
45 * The number of species is ns = 2 * np, with the first np being
46 * prey and the last np being predators. The coefficients a(i,j),
47 * b(i), d(i) are:
48 *
49 * a(i,i) = -AA (all i)
50 * a(i,j) = -GG (i <= np , j > np)
51 * a(i,j) = EE (i > np, j <= np)
52 * all other a(i,j) = 0
53 * b(i) = BB*(1+ alpha * x*y + beta*sin(4 pi x)*sin(4 pi y)) (i <= np)
54 * b(i) =-BB*(1+ alpha * x*y + beta*sin(4 pi x)*sin(4 pi y)) (i > np)
55 * d(i) = DPREY (i <= np)
56 * d(i) = DPRED (i > np)
57 *
58 * The various scalar parameters required are set using '#define'
59 * statements or directly in routine InitUserData. In this program,
60 * np = 1, ns = 2. The boundary conditions are homogeneous Neumann:
61 * normal derivative = 0.
62 *
63 * A polynomial in x and y is used to set the initial values of the
64 * first np variables (the prey variables) at each x,y location,
65 * while initial values for the remaining (predator) variables are
66 * set to a flat value, which is corrected by IDACalcIC.
67 *
68 * The PDEs are discretized by central differencing on a MX by MY
69 * mesh.
70 *
71 * The DAE system is solved by IDA using the banded linear solver.
72 * Output is printed at t = 0, .001, .01, .1, .4, .7, 1.
73 * -----------------------------------------------------------------
74 * References:
75 * [1] Peter N. Brown and Alan C. Hindmarsh,
76 * Reduced Storage Matrix Methods in Stiff ODE systems, Journal
77 * of Applied Mathematics and Computation, Vol. 31 (May 1989),
78 * pp. 40-91.
79 *
80 * [2] Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold,
81 * Using Krylov Methods in the Solution of Large-Scale
82 * Differential-Algebraic Systems, SIAM J. Sci. Comput., 15
83 * (1994), pp. 1467-1488.
84 *
85 * [3] Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold,
86 * Consistent Initial Condition Calculation for Differential-
87 * Algebraic Systems, SIAM J. Sci. Comput., 19 (1998),
88 * pp. 1495-1512.
89 * -----------------------------------------------------------------*/
90
91 #include <stdio.h>
92 #include <stdlib.h>
93 #include <math.h>
94
95 #include <ida/ida.h> /* prototypes for IDA fcts., consts. */
96 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
97 #include <sunmatrix/sunmatrix_band.h> /* access to band SUNMatrix */
98 #include <sunlinsol/sunlinsol_band.h> /* access to band SUNLinearSolver */
99 #include <sundials/sundials_types.h> /* definition of type realtype */
100
101 /* Problem Constants. */
102
103 #define NPREY 1 /* No. of prey (= no. of predators). */
104 #define NUM_SPECIES 2*NPREY
105
106 #define PI RCONST(3.1415926535898)
107 #define FOURPI (RCONST(4.0)*PI)
108
109 #define MX 20 /* MX = number of x mesh points */
110 #define MY 20 /* MY = number of y mesh points */
111 #define NSMX (NUM_SPECIES * MX)
112 #define NEQ (NUM_SPECIES*MX*MY)
113 #define AA RCONST(1.0) /* Coefficient in above eqns. for a */
114 #define EE RCONST(10000.) /* Coefficient in above eqns. for a */
115 #define GG RCONST(0.5e-6) /* Coefficient in above eqns. for a */
116 #define BB RCONST(1.0) /* Coefficient in above eqns. for b */
117 #define DPREY RCONST(1.0) /* Coefficient in above eqns. for d */
118 #define DPRED RCONST(0.05) /* Coefficient in above eqns. for d */
119 #define ALPHA RCONST(50.) /* Coefficient alpha in above eqns. */
120 #define BETA RCONST(1000.) /* Coefficient beta in above eqns. */
121 #define AX RCONST(1.0) /* Total range of x variable */
122 #define AY RCONST(1.0) /* Total range of y variable */
123 #define RTOL RCONST(1.e-5) /* Relative tolerance */
124 #define ATOL RCONST(1.e-5) /* Absolute tolerance */
125 #define NOUT 6 /* Number of output times */
126 #define TMULT RCONST(10.0) /* Multiplier for tout values */
127 #define TADD RCONST(0.3) /* Increment for tout values */
128 #define ZERO RCONST(0.)
129 #define ONE RCONST(1.0)
130
131 /*
132 * User-defined vector and accessor macro: IJ_Vptr.
133 * IJ_Vptr is defined in order to express the underlying 3-D structure of
134 * the dependent variable vector from its underlying 1-D storage (an N_Vector).
135 * IJ_Vptr(vv,i,j) returns a pointer to the location in vv corresponding to
136 * species index is = 0, x-index ix = i, and y-index jy = j.
137 */
138
139 #define IJ_Vptr(vv,i,j) (&NV_Ith_S(vv, (i)*NUM_SPECIES + (j)*NSMX))
140
141 /* Type: UserData. Contains problem constants, etc. */
142
143 typedef struct {
144 sunindextype Neq, ns, np, mx, my;
145 realtype dx, dy, **acoef;
146 realtype cox[NUM_SPECIES], coy[NUM_SPECIES], bcoef[NUM_SPECIES];
147 N_Vector rates;
148 } *UserData;
149
150 /* Prototypes for functions called by the IDA Solver. */
151
152 static int resweb(realtype time, N_Vector cc, N_Vector cp, N_Vector resval,
153 void *user_data);
154
155 /* Prototypes for private Helper Functions. */
156
157 static void InitUserData(UserData webdata);
158 static void SetInitialProfiles(N_Vector cc, N_Vector cp, N_Vector id,
159 UserData webdata);
160 static void PrintHeader(sunindextype mu, sunindextype ml, realtype rtol, realtype atol);
161 static void PrintOutput(void *mem, N_Vector c, realtype t);
162 static void PrintFinalStats(void *mem);
163 static void Fweb(realtype tcalc, N_Vector cc, N_Vector crate, UserData webdata);
164 static void WebRates(realtype xx, realtype yy, realtype *cxy, realtype *ratesxy,
165 UserData webdata);
166 static realtype dotprod(sunindextype size, realtype *x1, realtype *x2);
167 static int check_retval(void *returnvalue, const char *funcname, int opt);
168
169 /*
170 *--------------------------------------------------------------------
171 * MAIN PROGRAM
172 *--------------------------------------------------------------------
173 */
174
main()175 int main()
176 {
177 void *mem;
178 UserData webdata;
179 N_Vector cc, cp, id;
180 int iout, retval;
181 sunindextype mu, ml;
182 realtype rtol, atol, t0, tout, tret;
183 SUNMatrix A;
184 SUNLinearSolver LS;
185
186 mem = NULL;
187 webdata = NULL;
188 cc = cp = id = NULL;
189 A = NULL;
190 LS = NULL;
191
192 /* Allocate and initialize user data block webdata. */
193
194 webdata = (UserData) malloc(sizeof *webdata);
195 webdata->rates = N_VNew_Serial(NEQ);
196 webdata->acoef = newDenseMat(NUM_SPECIES, NUM_SPECIES);
197
198 InitUserData(webdata);
199
200 /* Allocate N-vectors and initialize cc, cp, and id. */
201
202 cc = N_VNew_Serial(NEQ);
203 if(check_retval((void *)cc, "N_VNew_Serial", 0)) return(1);
204
205 cp = N_VNew_Serial(NEQ);
206 if(check_retval((void *)cp, "N_VNew_Serial", 0)) return(1);
207
208 id = N_VNew_Serial(NEQ);
209 if(check_retval((void *)id, "N_VNew_Serial", 0)) return(1);
210
211 SetInitialProfiles(cc, cp, id, webdata);
212
213 /* Set remaining inputs to IDAMalloc. */
214
215 t0 = ZERO;
216 rtol = RTOL;
217 atol = ATOL;
218
219 /* Call IDACreate and IDAMalloc to initialize IDA. */
220
221 mem = IDACreate();
222 if(check_retval((void *)mem, "IDACreate", 0)) return(1);
223
224 retval = IDASetUserData(mem, webdata);
225 if(check_retval(&retval, "IDASetUserData", 1)) return(1);
226
227 retval = IDASetId(mem, id);
228 if(check_retval(&retval, "IDASetId", 1)) return(1);
229
230 retval = IDAInit(mem, resweb, t0, cc, cp);
231 if(check_retval(&retval, "IDAInit", 1)) return(1);
232
233 retval = IDASStolerances(mem, rtol, atol);
234 if(check_retval(&retval, "IDASStolerances", 1)) return(1);
235
236 /* Create banded SUNMatrix for use in linear solves */
237 mu = ml = NSMX;
238 A = SUNBandMatrix(NEQ, mu, ml);
239 if(check_retval((void *)A, "SUNBandMatrix", 0)) return(1);
240
241 /* Create banded SUNLinearSolver object */
242 LS = SUNLinSol_Band(cc, A);
243 if(check_retval((void *)LS, "SUNLinSol_Band", 0)) return(1);
244
245 /* Attach the matrix and linear solver */
246 retval = IDASetLinearSolver(mem, LS, A);
247 if(check_retval(&retval, "IDASetLinearSolver", 1)) return(1);
248
249 /* Call IDACalcIC (with default options) to correct the initial values. */
250
251 tout = RCONST(0.001);
252 retval = IDACalcIC(mem, IDA_YA_YDP_INIT, tout);
253 if(check_retval(&retval, "IDACalcIC", 1)) return(1);
254
255 /* Print heading, basic parameters, and initial values. */
256
257 PrintHeader(mu, ml, rtol, atol);
258 PrintOutput(mem, cc, ZERO);
259
260 /* Loop over iout, call IDASolve (normal mode), print selected output. */
261
262 for (iout = 1; iout <= NOUT; iout++) {
263
264 retval = IDASolve(mem, tout, &tret, cc, cp, IDA_NORMAL);
265 if(check_retval(&retval, "IDASolve", 1)) return(retval);
266
267 PrintOutput(mem, cc, tret);
268
269 if (iout < 3) tout *= TMULT; else tout += TADD;
270
271 }
272
273 /* Print final statistics and free memory. */
274
275 PrintFinalStats(mem);
276
277 /* Free memory */
278
279 IDAFree(&mem);
280 SUNLinSolFree(LS);
281 SUNMatDestroy(A);
282
283 N_VDestroy(cc);
284 N_VDestroy(cp);
285 N_VDestroy(id);
286
287
288 destroyMat(webdata->acoef);
289 N_VDestroy(webdata->rates);
290 free(webdata);
291
292 return(0);
293 }
294
295 /* Define lines for readability in later routines */
296
297 #define acoef (webdata->acoef)
298 #define bcoef (webdata->bcoef)
299 #define cox (webdata->cox)
300 #define coy (webdata->coy)
301
302 /*
303 *--------------------------------------------------------------------
304 * FUNCTIONS CALLED BY IDA
305 *--------------------------------------------------------------------
306 */
307
308 /*
309 * resweb: System residual function for predator-prey system.
310 * This routine calls Fweb to get all the right-hand sides of the
311 * equations, then loads the residual vector accordingly,
312 * using cp in the case of prey species.
313 */
314
resweb(realtype tt,N_Vector cc,N_Vector cp,N_Vector res,void * user_data)315 static int resweb(realtype tt, N_Vector cc, N_Vector cp,
316 N_Vector res, void *user_data)
317 {
318 sunindextype jx, jy, is, yloc, loc, np;
319 realtype *resv, *cpv;
320 UserData webdata;
321
322 webdata = (UserData)user_data;
323
324 cpv = N_VGetArrayPointer(cp);
325 resv = N_VGetArrayPointer(res);
326 np = webdata->np;
327
328 /* Call Fweb to set res to vector of right-hand sides. */
329 Fweb(tt, cc, res, webdata);
330
331 /* Loop over all grid points, setting residual values appropriately
332 for differential or algebraic components. */
333
334 for (jy = 0; jy < MY; jy++) {
335 yloc = NSMX * jy;
336 for (jx = 0; jx < MX; jx++) {
337 loc = yloc + NUM_SPECIES * jx;
338 for (is = 0; is < NUM_SPECIES; is++) {
339 if (is < np)
340 resv[loc+is] = cpv[loc+is] - resv[loc+is];
341 else
342 resv[loc+is] = -resv[loc+is];
343 }
344 }
345 }
346
347 return(0);
348
349 }
350
351 /*
352 *--------------------------------------------------------------------
353 * PRIVATE FUNCTIONS
354 *--------------------------------------------------------------------
355 */
356
357 /*
358 * InitUserData: Load problem constants in webdata (of type UserData).
359 */
360
InitUserData(UserData webdata)361 static void InitUserData(UserData webdata)
362 {
363 sunindextype i, j, np;
364 realtype *a1,*a2, *a3, *a4, dx2, dy2;
365
366 webdata->mx = MX;
367 webdata->my = MY;
368 webdata->ns = NUM_SPECIES;
369 webdata->np = NPREY;
370 webdata->dx = AX/(MX-1);
371 webdata->dy = AY/(MY-1);
372 webdata->Neq= NEQ;
373
374 /* Set up the coefficients a and b, and others found in the equations. */
375 np = webdata->np;
376 dx2 = (webdata->dx)*(webdata->dx); dy2 = (webdata->dy)*(webdata->dy);
377
378 for (i = 0; i < np; i++) {
379 a1 = &(acoef[i][np]);
380 a2 = &(acoef[i+np][0]);
381 a3 = &(acoef[i][0]);
382 a4 = &(acoef[i+np][np]);
383 /* Fill in the portion of acoef in the four quadrants, row by row. */
384 for (j = 0; j < np; j++) {
385 *a1++ = -GG;
386 *a2++ = EE;
387 *a3++ = ZERO;
388 *a4++ = ZERO;
389 }
390
391 /* Reset the diagonal elements of acoef to -AA. */
392 acoef[i][i] = -AA; acoef[i+np][i+np] = -AA;
393
394 /* Set coefficients for b and diffusion terms. */
395 bcoef[i] = BB; bcoef[i+np] = -BB;
396 cox[i] = DPREY/dx2; cox[i+np] = DPRED/dx2;
397 coy[i] = DPREY/dy2; coy[i+np] = DPRED/dy2;
398 }
399
400 }
401
402 /*
403 * SetInitialProfiles: Set initial conditions in cc, cp, and id.
404 * A polynomial profile is used for the prey cc values, and a constant
405 * (1.0e5) is loaded as the initial guess for the predator cc values.
406 * The id values are set to 1 for the prey and 0 for the predators.
407 * The prey cp values are set according to the given system, and
408 * the predator cp values are set to zero.
409 */
410
SetInitialProfiles(N_Vector cc,N_Vector cp,N_Vector id,UserData webdata)411 static void SetInitialProfiles(N_Vector cc, N_Vector cp, N_Vector id,
412 UserData webdata)
413 {
414 sunindextype loc, yloc, is, jx, jy, np;
415 realtype xx, yy, xyfactor;
416 realtype *ccv, *cpv, *idv;
417
418 ccv = N_VGetArrayPointer(cc);
419 cpv = N_VGetArrayPointer(cp);
420 idv = N_VGetArrayPointer(id);
421 np = webdata->np;
422
423 /* Loop over grid, load cc values and id values. */
424 for (jy = 0; jy < MY; jy++) {
425 yy = jy * webdata->dy;
426 yloc = NSMX * jy;
427 for (jx = 0; jx < MX; jx++) {
428 xx = jx * webdata->dx;
429 xyfactor = RCONST(16.0)*xx*(ONE-xx)*yy*(ONE-yy);
430 xyfactor *= xyfactor;
431 loc = yloc + NUM_SPECIES*jx;
432
433 for (is = 0; is < NUM_SPECIES; is++) {
434 if (is < np) {
435 ccv[loc+is] = RCONST(10.0) + (realtype)(is+1) * xyfactor;
436 idv[loc+is] = ONE;
437 }
438 else {
439 ccv[loc+is] = RCONST(1.0e5);
440 idv[loc+is] = ZERO;
441 }
442 }
443 }
444 }
445
446 /* Set c' for the prey by calling the function Fweb. */
447 Fweb(ZERO, cc, cp, webdata);
448
449 /* Set c' for predators to 0. */
450 for (jy = 0; jy < MY; jy++) {
451 yloc = NSMX * jy;
452 for (jx = 0; jx < MX; jx++) {
453 loc = yloc + NUM_SPECIES * jx;
454 for (is = np; is < NUM_SPECIES; is++) {
455 cpv[loc+is] = ZERO;
456 }
457 }
458 }
459 }
460
461 /*
462 * Print first lines of output (problem description)
463 */
464
PrintHeader(sunindextype mu,sunindextype ml,realtype rtol,realtype atol)465 static void PrintHeader(sunindextype mu, sunindextype ml, realtype rtol, realtype atol)
466 {
467 printf("\nidaFoodWeb_bnd: Predator-prey DAE serial example problem for IDA \n\n");
468 printf("Number of species ns: %d", NUM_SPECIES);
469 printf(" Mesh dimensions: %d x %d", MX, MY);
470 printf(" System size: %d\n", NEQ);
471 #if defined(SUNDIALS_EXTENDED_PRECISION)
472 printf("Tolerance parameters: rtol = %Lg atol = %Lg\n", rtol, atol);
473 #elif defined(SUNDIALS_DOUBLE_PRECISION)
474 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
475 #else
476 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
477 #endif
478 printf("Linear solver: BAND, Band parameters mu = %ld, ml = %ld\n", (long int) mu, (long int) ml);
479 printf("CalcIC called to correct initial predator concentrations.\n\n");
480 printf("-----------------------------------------------------------\n");
481 printf(" t bottom-left top-right");
482 printf(" | nst k h\n");
483 printf("-----------------------------------------------------------\n\n");
484
485 }
486
487 /*
488 * PrintOutput: Print output values at output time t = tt.
489 * Selected run statistics are printed. Then values of the concentrations
490 * are printed for the bottom left and top right grid points only.
491 */
492
PrintOutput(void * mem,N_Vector c,realtype t)493 static void PrintOutput(void *mem, N_Vector c, realtype t)
494 {
495 int i, kused, retval;
496 long int nst;
497 realtype *c_bl, *c_tr, hused;
498
499 retval = IDAGetLastOrder(mem, &kused);
500 check_retval(&retval, "IDAGetLastOrder", 1);
501 retval = IDAGetNumSteps(mem, &nst);
502 check_retval(&retval, "IDAGetNumSteps", 1);
503 retval = IDAGetLastStep(mem, &hused);
504 check_retval(&retval, "IDAGetLastStep", 1);
505
506 c_bl = IJ_Vptr(c,0,0);
507 c_tr = IJ_Vptr(c,MX-1,MY-1);
508
509 #if defined(SUNDIALS_EXTENDED_PRECISION)
510 printf("%8.2Le %12.4Le %12.4Le | %3ld %1d %12.4Le\n",
511 t, c_bl[0], c_tr[0], nst, kused, hused);
512 for (i=1;i<NUM_SPECIES;i++)
513 printf(" %12.4Le %12.4Le |\n",c_bl[i],c_tr[i]);
514 #elif defined(SUNDIALS_DOUBLE_PRECISION)
515 printf("%8.2e %12.4e %12.4e | %3ld %1d %12.4e\n",
516 t, c_bl[0], c_tr[0], nst, kused, hused);
517 for (i=1;i<NUM_SPECIES;i++)
518 printf(" %12.4e %12.4e |\n",c_bl[i],c_tr[i]);
519 #else
520 printf("%8.2e %12.4e %12.4e | %3ld %1d %12.4e\n",
521 t, c_bl[0], c_tr[0], nst, kused, hused);
522 for (i=1;i<NUM_SPECIES;i++)
523 printf(" %12.4e %12.4e |\n",c_bl[i],c_tr[i]);
524 #endif
525
526 printf("\n");
527 }
528
529 /*
530 * PrintFinalStats: Print final run data contained in iopt.
531 */
532
PrintFinalStats(void * mem)533 static void PrintFinalStats(void *mem)
534 {
535 long int nst, nre, nreLS, nni, nje, netf, ncfn;
536 int retval;
537
538 retval = IDAGetNumSteps(mem, &nst);
539 check_retval(&retval, "IDAGetNumSteps", 1);
540 retval = IDAGetNumNonlinSolvIters(mem, &nni);
541 check_retval(&retval, "IDAGetNumNonlinSolvIters", 1);
542 retval = IDAGetNumResEvals(mem, &nre);
543 check_retval(&retval, "IDAGetNumResEvals", 1);
544 retval = IDAGetNumErrTestFails(mem, &netf);
545 check_retval(&retval, "IDAGetNumErrTestFails", 1);
546 retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
547 check_retval(&retval, "IDAGetNumNonlinSolvConvFails", 1);
548 retval = IDAGetNumJacEvals(mem, &nje);
549 check_retval(&retval, "IDAGetNumJacEvals", 1);
550 retval = IDAGetNumLinResEvals(mem, &nreLS);
551 check_retval(&retval, "IDAGetNumLinResEvals", 1);
552
553 printf("-----------------------------------------------------------\n");
554 printf("Final run statistics: \n\n");
555 printf("Number of steps = %ld\n", nst);
556 printf("Number of residual evaluations = %ld\n", nre+nreLS);
557 printf("Number of Jacobian evaluations = %ld\n", nje);
558 printf("Number of nonlinear iterations = %ld\n", nni);
559 printf("Number of error test failures = %ld\n", netf);
560 printf("Number of nonlinear conv. failures = %ld\n", ncfn);
561
562 }
563
564 /*
565 * Fweb: Rate function for the food-web problem.
566 * This routine computes the right-hand sides of the system equations,
567 * consisting of the diffusion term and interaction term.
568 * The interaction term is computed by the function WebRates.
569 */
570
Fweb(realtype tcalc,N_Vector cc,N_Vector crate,UserData webdata)571 static void Fweb(realtype tcalc, N_Vector cc, N_Vector crate,
572 UserData webdata)
573 {
574 sunindextype jx, jy, is, idyu, idyl, idxu, idxl;
575 realtype xx, yy, *cxy, *ratesxy, *cratexy, dcyli, dcyui, dcxli, dcxui;
576
577 /* Loop over grid points, evaluate interaction vector (length ns),
578 form diffusion difference terms, and load crate. */
579
580 for (jy = 0; jy < MY; jy++) {
581 yy = (webdata->dy) * jy ;
582 idyu = (jy!=MY-1) ? NSMX : -NSMX;
583 idyl = (jy!= 0 ) ? NSMX : -NSMX;
584
585 for (jx = 0; jx < MX; jx++) {
586 xx = (webdata->dx) * jx;
587 idxu = (jx!= MX-1) ? NUM_SPECIES : -NUM_SPECIES;
588 idxl = (jx!= 0 ) ? NUM_SPECIES : -NUM_SPECIES;
589 cxy = IJ_Vptr(cc,jx,jy);
590 ratesxy = IJ_Vptr(webdata->rates,jx,jy);
591 cratexy = IJ_Vptr(crate,jx,jy);
592
593 /* Get interaction vector at this grid point. */
594 WebRates(xx, yy, cxy, ratesxy, webdata);
595
596 /* Loop over species, do differencing, load crate segment. */
597 for (is = 0; is < NUM_SPECIES; is++) {
598
599 /* Differencing in y. */
600 dcyli = *(cxy+is) - *(cxy - idyl + is) ;
601 dcyui = *(cxy + idyu + is) - *(cxy+is);
602
603 /* Differencing in x. */
604 dcxli = *(cxy+is) - *(cxy - idxl + is);
605 dcxui = *(cxy + idxu +is) - *(cxy+is);
606
607 /* Compute the crate values at (xx,yy). */
608 cratexy[is] = coy[is] * (dcyui - dcyli) +
609 cox[is] * (dcxui - dcxli) + ratesxy[is];
610
611 } /* End is loop */
612 } /* End of jx loop */
613 } /* End of jy loop */
614
615 }
616
617 /*
618 * WebRates: Evaluate reaction rates at a given spatial point.
619 * At a given (x,y), evaluate the array of ns reaction terms R.
620 */
621
WebRates(realtype xx,realtype yy,realtype * cxy,realtype * ratesxy,UserData webdata)622 static void WebRates(realtype xx, realtype yy, realtype *cxy, realtype *ratesxy,
623 UserData webdata)
624 {
625 int is;
626 realtype fac;
627
628 for (is = 0; is < NUM_SPECIES; is++)
629 ratesxy[is] = dotprod(NUM_SPECIES, cxy, acoef[is]);
630
631 fac = ONE + ALPHA*xx*yy + BETA*sin(FOURPI*xx)*sin(FOURPI*yy);
632
633 for (is = 0; is < NUM_SPECIES; is++)
634 ratesxy[is] = cxy[is]*( bcoef[is]*fac + ratesxy[is] );
635
636 }
637
638 /*
639 * dotprod: dot product routine for realtype arrays, for use by WebRates.
640 */
641
dotprod(sunindextype size,realtype * x1,realtype * x2)642 static realtype dotprod(sunindextype size, realtype *x1, realtype *x2)
643 {
644 sunindextype i;
645 realtype *xx1, *xx2, temp = ZERO;
646
647 xx1 = x1; xx2 = x2;
648 for (i = 0; i < size; i++) temp += (*xx1++) * (*xx2++);
649 return(temp);
650
651 }
652
653 /*
654 * Check function return value...
655 * opt == 0 means SUNDIALS function allocates memory so check if
656 * returned NULL pointer
657 * opt == 1 means SUNDIALS function returns an integer value so check if
658 * retval < 0
659 * opt == 2 means function allocates memory so check if returned
660 * NULL pointer
661 */
662
check_retval(void * returnvalue,const char * funcname,int opt)663 static int check_retval(void *returnvalue, const char *funcname, int opt)
664 {
665 int *retval;
666
667 if (opt == 0 && returnvalue == NULL) {
668 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
669 fprintf(stderr,
670 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
671 funcname);
672 return(1);
673 } else if (opt == 1) {
674 /* Check if retval < 0 */
675 retval = (int *) returnvalue;
676 if (*retval < 0) {
677 fprintf(stderr,
678 "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
679 funcname, *retval);
680 return(1);
681 }
682 } else if (opt == 2 && returnvalue == NULL) {
683 /* Check if function returned NULL pointer - no memory allocated */
684 fprintf(stderr,
685 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
686 funcname);
687 return(1);
688 }
689
690 return(0);
691 }
692