1 /* -----------------------------------------------------------------
2 * Programmer(s): Ting Yan @ SMU
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 * Example program for IDA: Food web problem, OpenMP, GMRES,
15 * user-supplied preconditioner
16 *
17 * This example program uses the SPGMR as the 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 SPGMR 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 <sunlinsol/sunlinsol_spgmr.h> /* access to spgmr SUNLinearSolver */
98 #include <sundials/sundials_dense.h> /* use generic dense solver in precond. */
99 #include <sundials/sundials_types.h> /* definition of type realtype */
100
101 /* helpful macros */
102
103 #ifndef MAX
104 #define MAX(A, B) ((A) > (B) ? (A) : (B))
105 #endif
106
107 /* Problem Constants. */
108
109 #define NPREY 1 /* No. of prey (= no. of predators). */
110 #define NUM_SPECIES 2*NPREY
111
112 #define PI RCONST(3.1415926535898)
113 #define FOURPI (RCONST(4.0)*PI)
114
115 #define MX 20 /* MX = number of x mesh points */
116 #define MY 20 /* MY = number of y mesh points */
117 #define NSMX (NUM_SPECIES * MX)
118 #define NEQ (NUM_SPECIES*MX*MY)
119 #define AA RCONST(1.0) /* Coefficient in above eqns. for a */
120 #define EE RCONST(10000.) /* Coefficient in above eqns. for a */
121 #define GG RCONST(0.5e-6) /* Coefficient in above eqns. for a */
122 #define BB RCONST(1.0) /* Coefficient in above eqns. for b */
123 #define DPREY RCONST(1.0) /* Coefficient in above eqns. for d */
124 #define DPRED RCONST(0.05) /* Coefficient in above eqns. for d */
125 #define ALPHA RCONST(50.) /* Coefficient alpha in above eqns. */
126 #define BETA RCONST(1000.) /* Coefficient beta in above eqns. */
127 #define AX RCONST(1.0) /* Total range of x variable */
128 #define AY RCONST(1.0) /* Total range of y variable */
129 #define RTOL RCONST(1.e-5) /* Relative tolerance */
130 #define ATOL RCONST(1.e-5) /* Absolute tolerance */
131 #define NOUT 6 /* Number of output times */
132 #define TMULT RCONST(10.0) /* Multiplier for tout values */
133 #define TADD RCONST(0.3) /* Increment for tout values */
134 #define ZERO RCONST(0.)
135 #define ONE RCONST(1.0)
136
137 /*
138 * User-defined vector and accessor macro: IJ_Vptr.
139 * IJ_Vptr is defined in order to express the underlying 3-D structure of
140 * the dependent variable vector from its underlying 1-D storage (an N_Vector).
141 * IJ_Vptr(vv,i,j) returns a pointer to the location in vv corresponding to
142 * species index is = 0, x-index ix = i, and y-index jy = j.
143 */
144
145 #define IJ_Vptr(vv,i,j) (&NV_Ith_S(vv, (i)*NUM_SPECIES + (j)*NSMX))
146
147 /* Type: UserData. Contains problem constants, etc. */
148
149 typedef struct {
150 sunindextype Neq, ns, np, mx, my;
151 realtype dx, dy, **acoef;
152 realtype cox[NUM_SPECIES], coy[NUM_SPECIES], bcoef[NUM_SPECIES];
153 realtype **PP[MX][MY];
154 sunindextype *pivot[MX][MY];
155 N_Vector rates;
156 N_Vector ewt;
157 void *ida_mem;
158 } *UserData;
159
160 /* Prototypes for functions called by the IDA Solver. */
161
162 static int resweb(realtype time, N_Vector cc, N_Vector cp, N_Vector resval,
163 void *user_data);
164
165 static int Precond(realtype tt,
166 N_Vector cc, N_Vector cp, N_Vector rr,
167 realtype cj, void *user_data);
168
169 static int PSolve(realtype tt,
170 N_Vector cc, N_Vector cp, N_Vector rr,
171 N_Vector rvec, N_Vector zvec,
172 realtype cj, realtype delta, void *user_data);
173
174 /* Prototypes for private Helper Functions. */
175
176 static void InitUserData(UserData webdata);
177 static void SetInitialProfiles(N_Vector cc, N_Vector cp, N_Vector id,
178 UserData webdata);
179 static void PrintHeader(int maxl, realtype rtol, realtype atol);
180 static void PrintOutput(void *mem, N_Vector c, realtype t);
181 static void PrintFinalStats(void *mem);
182 static void Fweb(realtype tcalc, N_Vector cc, N_Vector crate, UserData webdata);
183 static void WebRates(realtype xx, realtype yy, realtype *cxy, realtype *ratesxy,
184 UserData webdata);
185 static realtype dotprod(sunindextype size, realtype *x1, realtype *x2);
186 static int check_retval(void *returnvalue, char *funcname, int opt);
187
188 /*
189 *--------------------------------------------------------------------
190 * MAIN PROGRAM
191 *--------------------------------------------------------------------
192 */
193
main()194 int main()
195 {
196 void *mem;
197 UserData webdata;
198 N_Vector cc, cp, id;
199 int iout, jx, jy, retval;
200 int maxl;
201 realtype rtol, atol, t0, tout, tret;
202 SUNLinearSolver LS;
203
204 mem = NULL;
205 webdata = NULL;
206 cc = cp = id = NULL;
207 LS = NULL;
208
209 /* Allocate and initialize user data block webdata. */
210
211 webdata = (UserData) malloc(sizeof *webdata);
212 webdata->rates = N_VNew_Serial(NEQ);
213 webdata->acoef = newDenseMat(NUM_SPECIES, NUM_SPECIES);
214 webdata->ewt = N_VNew_Serial(NEQ);
215 for (jx = 0; jx < MX; jx++) {
216 for (jy = 0; jy < MY; jy++) {
217 (webdata->pivot)[jx][jy] = newIndexArray(NUM_SPECIES);
218 (webdata->PP)[jx][jy] = newDenseMat(NUM_SPECIES, NUM_SPECIES);
219 }
220 }
221
222 InitUserData(webdata);
223
224 /* Allocate N-vectors and initialize cc, cp, and id. */
225
226 cc = N_VNew_Serial(NEQ);
227 if(check_retval((void *)cc, "N_VNew_Serial", 0)) return(1);
228
229 cp = N_VNew_Serial(NEQ);
230 if(check_retval((void *)cp, "N_VNew_Serial", 0)) return(1);
231
232 id = N_VNew_Serial(NEQ);
233 if(check_retval((void *)id, "N_VNew_Serial", 0)) return(1);
234
235 SetInitialProfiles(cc, cp, id, webdata);
236
237 /* Set remaining inputs to IDAMalloc. */
238
239 t0 = ZERO;
240 rtol = RTOL;
241 atol = ATOL;
242
243 /* Call IDACreate and IDAMalloc to initialize IDA. */
244
245 mem = IDACreate();
246 if(check_retval((void *)mem, "IDACreate", 0)) return(1);
247
248 retval = IDASetUserData(mem, webdata);
249 if(check_retval(&retval, "IDASetUserData", 1)) return(1);
250
251 retval = IDASetId(mem, id);
252 if(check_retval(&retval, "IDASetId", 1)) return(1);
253
254 retval = IDAInit(mem, resweb, t0, cc, cp);
255 if(check_retval(&retval, "IDAInit", 1)) return(1);
256
257 retval = IDASStolerances(mem, rtol, atol);
258 if(check_retval(&retval, "IDASStolerances", 1)) return(1);
259
260 webdata->ida_mem = mem;
261
262 /* Create the linear solver SUNLinSol_SPGMR with left preconditioning
263 and maximum Krylov dimension maxl */
264 maxl = 16;
265 LS = SUNLinSol_SPGMR(cc, PREC_LEFT, maxl);
266 if(check_retval((void *)LS, "SUNLinSol_SPGMR", 0)) return(1);
267
268 /* IDA recommends allowing up to 5 restarts (default is 0) */
269 retval = SUNLinSol_SPGMRSetMaxRestarts(LS, 5);
270 if(check_retval(&retval, "SUNLinSol_SPGMRSetMaxRestarts", 1)) return(1);
271
272 /* Attach the linear sovler */
273 retval = IDASetLinearSolver(mem, LS, NULL);
274 if(check_retval(&retval, "IDASetLinearSolver", 1)) return(1);
275
276 /* Set the preconditioner solve and setup functions */
277 retval = IDASetPreconditioner(mem, Precond, PSolve);
278 if(check_retval(&retval, "IDASetPreconditioner", 1)) return(1);
279
280 /* Call IDACalcIC (with default options) to correct the initial values. */
281
282 tout = RCONST(0.001);
283 retval = IDACalcIC(mem, IDA_YA_YDP_INIT, tout);
284 if(check_retval(&retval, "IDACalcIC", 1)) return(1);
285
286 /* Print heading, basic parameters, and initial values. */
287
288 PrintHeader(maxl, rtol, atol);
289 PrintOutput(mem, cc, ZERO);
290
291 /* Loop over iout, call IDASolve (normal mode), print selected output. */
292
293 for (iout = 1; iout <= NOUT; iout++) {
294
295 retval = IDASolve(mem, tout, &tret, cc, cp, IDA_NORMAL);
296 if(check_retval(&retval, "IDASolve", 1)) return(retval);
297
298 PrintOutput(mem, cc, tret);
299
300 if (iout < 3) tout *= TMULT; else tout += TADD;
301
302 }
303
304 /* Print final statistics and free memory. */
305
306 PrintFinalStats(mem);
307
308 /* Free memory */
309
310 IDAFree(&mem);
311 SUNLinSolFree(LS);
312
313 N_VDestroy(cc);
314 N_VDestroy(cp);
315 N_VDestroy(id);
316
317
318 destroyMat(webdata->acoef);
319 N_VDestroy(webdata->rates);
320 N_VDestroy(webdata->ewt);
321 for (jx = 0; jx < MX; jx++) {
322 for (jy = 0; jy < MY; jy ++) {
323 destroyArray((webdata->pivot)[jx][jy]);
324 destroyMat((webdata->PP)[jx][jy]);
325 }
326 }
327 free(webdata);
328
329 return(0);
330 }
331
332 /* Define lines for readability in later routines */
333
334 #define acoef (webdata->acoef)
335 #define bcoef (webdata->bcoef)
336 #define cox (webdata->cox)
337 #define coy (webdata->coy)
338
339 /*
340 *--------------------------------------------------------------------
341 * FUNCTIONS CALLED BY IDA
342 *--------------------------------------------------------------------
343 */
344
345 /*
346 * resweb: System residual function for predator-prey system.
347 * This routine calls Fweb to get all the right-hand sides of the
348 * equations, then loads the residual vector accordingly,
349 * using cp in the case of prey species.
350 */
351
resweb(realtype tt,N_Vector cc,N_Vector cp,N_Vector res,void * user_data)352 static int resweb(realtype tt, N_Vector cc, N_Vector cp,
353 N_Vector res, void *user_data)
354 {
355 sunindextype jx, jy, is, yloc, loc, np;
356 realtype *resv, *cpv;
357 UserData webdata;
358
359 webdata = (UserData)user_data;
360
361 cpv = N_VGetArrayPointer(cp);
362 resv = N_VGetArrayPointer(res);;
363 np = webdata->np;
364
365 /* Call Fweb to set res to vector of right-hand sides. */
366 Fweb(tt, cc, res, webdata);
367
368 /* Loop over all grid points, setting residual values appropriately
369 for differential or algebraic components. */
370
371 for (jy = 0; jy < MY; jy++) {
372 yloc = NSMX * jy;
373 for (jx = 0; jx < MX; jx++) {
374 loc = yloc + NUM_SPECIES * jx;
375 for (is = 0; is < NUM_SPECIES; is++) {
376 if (is < np)
377 resv[loc+is] = cpv[loc+is] - resv[loc+is];
378 else
379 resv[loc+is] = -resv[loc+is];
380 }
381 }
382 }
383
384 return(0);
385
386 }
387
388
Precond(realtype tt,N_Vector cc,N_Vector cp,N_Vector rr,realtype cj,void * user_data)389 static int Precond(realtype tt,
390 N_Vector cc, N_Vector cp, N_Vector rr,
391 realtype cj, void *user_data)
392 {
393 int retval;
394 sunindextype ret;
395 realtype uround, xx, yy, del_x, del_y;
396 realtype **Pxy, *ratesxy, *Pxycol, *cxy, *cpxy, *ewtxy, cctmp;
397 realtype inc, fac, sqru, perturb_rates[NUM_SPECIES];
398 int is, js, jx, jy;
399 void *mem;
400 N_Vector ewt;
401 realtype hh;
402 UserData webdata;
403
404 webdata = (UserData) user_data;
405 del_x = webdata->dx;
406 del_y = webdata->dy;
407
408 uround = UNIT_ROUNDOFF;
409 sqru = sqrt(uround);
410
411 mem = webdata->ida_mem;
412 ewt = webdata->ewt;
413 retval = IDAGetErrWeights(mem, ewt);
414 if(check_retval(&retval, "IDAGetErrWeights", 1)) return(1);
415 retval = IDAGetCurrentStep(mem, &hh);
416 if(check_retval(&retval, "IDAGetCurrentStep", 1)) return(1);
417
418 for (jy = 0; jy < MY; jy++) {
419 yy = jy * del_y;
420
421 for (jx = 0; jx < MX; jx++) {
422 xx = jx * del_x;
423 Pxy = (webdata->PP)[jx][jy];
424 cxy = IJ_Vptr(cc, jx, jy);
425 cpxy = IJ_Vptr(cp, jx, jy);
426 ewtxy = IJ_Vptr(ewt, jx, jy);
427 ratesxy = IJ_Vptr((webdata->rates), jx, jy);
428
429 for (js = 0; js < NUM_SPECIES; js++) {
430 inc = sqru*(MAX(fabs(cxy[js]), MAX(hh*fabs(cpxy[js]), ONE/ewtxy[js])));
431 cctmp = cxy[js];
432 cxy[js] += inc;
433 fac = -ONE/inc;
434
435 WebRates(xx, yy, cxy, perturb_rates, webdata);
436
437 Pxycol = Pxy[js];
438
439 for (is = 0; is < NUM_SPECIES; is++)
440 Pxycol[is] = (perturb_rates[is] - ratesxy[is])*fac;
441
442 if (js < 1) Pxycol[js] += cj;
443
444 cxy[js] = cctmp;
445 }
446
447 ret = denseGETRF(Pxy, NUM_SPECIES, NUM_SPECIES, (webdata->pivot)[jx][jy]);
448
449 if (ret != 0) return(1);
450 }
451 }
452
453 return(0);
454
455 }
456
457
PSolve(realtype tt,N_Vector cc,N_Vector cp,N_Vector rr,N_Vector rvec,N_Vector zvec,realtype cj,realtype dalta,void * user_data)458 static int PSolve(realtype tt,
459 N_Vector cc, N_Vector cp, N_Vector rr,
460 N_Vector rvec, N_Vector zvec,
461 realtype cj, realtype dalta,
462 void *user_data)
463 {
464 realtype **Pxy, *zxy;
465 sunindextype *pivot;
466 int jx, jy;
467 UserData webdata;
468
469 webdata = (UserData) user_data;
470
471 N_VScale(ONE, rvec, zvec);
472
473 for (jx = 0; jx < MX; jx++) {
474 for (jy = 0; jy <MY; jy++) {
475
476 zxy = IJ_Vptr(zvec, jx, jy);
477 Pxy = (webdata->PP)[jx][jy];
478 pivot = (webdata->pivot)[jx][jy];
479 denseGETRS(Pxy, NUM_SPECIES, pivot, zxy);
480 }
481 }
482
483 return(0);
484
485 }
486
487
488 /*
489 *--------------------------------------------------------------------
490 * PRIVATE FUNCTIONS
491 *--------------------------------------------------------------------
492 */
493
494 /*
495 * InitUserData: Load problem constants in webdata (of type UserData).
496 */
497
InitUserData(UserData webdata)498 static void InitUserData(UserData webdata)
499 {
500 sunindextype i, j, np;
501 realtype *a1,*a2, *a3, *a4, dx2, dy2;
502
503 webdata->mx = MX;
504 webdata->my = MY;
505 webdata->ns = NUM_SPECIES;
506 webdata->np = NPREY;
507 webdata->dx = AX/(MX-1);
508 webdata->dy = AY/(MY-1);
509 webdata->Neq= NEQ;
510
511 /* Set up the coefficients a and b, and others found in the equations. */
512 np = webdata->np;
513 dx2 = (webdata->dx)*(webdata->dx); dy2 = (webdata->dy)*(webdata->dy);
514
515 for (i = 0; i < np; i++) {
516 a1 = &(acoef[i][np]);
517 a2 = &(acoef[i+np][0]);
518 a3 = &(acoef[i][0]);
519 a4 = &(acoef[i+np][np]);
520 /* Fill in the portion of acoef in the four quadrants, row by row. */
521 for (j = 0; j < np; j++) {
522 *a1++ = -GG;
523 *a2++ = EE;
524 *a3++ = ZERO;
525 *a4++ = ZERO;
526 }
527
528 /* Reset the diagonal elements of acoef to -AA. */
529 acoef[i][i] = -AA; acoef[i+np][i+np] = -AA;
530
531 /* Set coefficients for b and diffusion terms. */
532 bcoef[i] = BB; bcoef[i+np] = -BB;
533 cox[i] = DPREY/dx2; cox[i+np] = DPRED/dx2;
534 coy[i] = DPREY/dy2; coy[i+np] = DPRED/dy2;
535 }
536
537 }
538
539 /*
540 * SetInitialProfiles: Set initial conditions in cc, cp, and id.
541 * A polynomial profile is used for the prey cc values, and a constant
542 * (1.0e5) is loaded as the initial guess for the predator cc values.
543 * The id values are set to 1 for the prey and 0 for the predators.
544 * The prey cp values are set according to the given system, and
545 * the predator cp values are set to zero.
546 */
547
SetInitialProfiles(N_Vector cc,N_Vector cp,N_Vector id,UserData webdata)548 static void SetInitialProfiles(N_Vector cc, N_Vector cp, N_Vector id,
549 UserData webdata)
550 {
551 sunindextype loc, yloc, is, jx, jy, np;
552 realtype xx, yy, xyfactor;
553 realtype *ccv, *cpv, *idv;
554
555 ccv = N_VGetArrayPointer(cc);
556 cpv = N_VGetArrayPointer(cp) ;
557 idv = N_VGetArrayPointer(id);
558 np = webdata->np;
559
560 /* Loop over grid, load cc values and id values. */
561 for (jy = 0; jy < MY; jy++) {
562 yy = jy * webdata->dy;
563 yloc = NSMX * jy;
564 for (jx = 0; jx < MX; jx++) {
565 xx = jx * webdata->dx;
566 xyfactor = RCONST(16.0)*xx*(ONE-xx)*yy*(ONE-yy);
567 xyfactor *= xyfactor;
568 loc = yloc + NUM_SPECIES*jx;
569
570 for (is = 0; is < NUM_SPECIES; is++) {
571 if (is < np) {
572 ccv[loc+is] = RCONST(10.0) + (realtype)(is+1) * xyfactor;
573 idv[loc+is] = ONE;
574 }
575 else {
576 ccv[loc+is] = RCONST(1.0e5);
577 idv[loc+is] = ZERO;
578 }
579 }
580 }
581 }
582
583 /* Set c' for the prey by calling the function Fweb. */
584 Fweb(ZERO, cc, cp, webdata);
585
586 /* Set c' for predators to 0. */
587 for (jy = 0; jy < MY; jy++) {
588 yloc = NSMX * jy;
589 for (jx = 0; jx < MX; jx++) {
590 loc = yloc + NUM_SPECIES * jx;
591 for (is = np; is < NUM_SPECIES; is++) {
592 cpv[loc+is] = ZERO;
593 }
594 }
595 }
596 }
597
598 /*
599 * Print first lines of output (problem description)
600 */
601
PrintHeader(int maxl,realtype rtol,realtype atol)602 static void PrintHeader(int maxl, realtype rtol, realtype atol)
603 {
604 printf("\nidaFoodWeb_kry: Predator-prey DAE serial example problem for IDA \n\n");
605 printf("Number of species ns: %d", NUM_SPECIES);
606 printf(" Mesh dimensions: %d x %d", MX, MY);
607 printf(" System size: %d\n", NEQ);
608 #if defined(SUNDIALS_EXTENDED_PRECISION)
609 printf("Tolerance parameters: rtol = %Lg atol = %Lg\n", rtol, atol);
610 #elif defined(SUNDIALS_DOUBLE_PRECISION)
611 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
612 #else
613 printf("Tolerance parameters: rtol = %g atol = %g\n", rtol, atol);
614 #endif
615 printf("Linear solver: SPGMR, SPGMR parameters maxl = %d\n", maxl);
616 printf("CalcIC called to correct initial predator concentrations.\n\n");
617 printf("-----------------------------------------------------------\n");
618 printf(" t bottom-left top-right");
619 printf(" | nst k h\n");
620 printf("-----------------------------------------------------------\n\n");
621
622 }
623
624 /*
625 * PrintOutput: Print output values at output time t = tt.
626 * Selected run statistics are printed. Then values of the concentrations
627 * are printed for the bottom left and top right grid points only.
628 */
629
PrintOutput(void * mem,N_Vector c,realtype t)630 static void PrintOutput(void *mem, N_Vector c, realtype t)
631 {
632 int i, kused, retval;
633 long int nst;
634 realtype *c_bl, *c_tr, hused;
635
636 retval = IDAGetLastOrder(mem, &kused);
637 check_retval(&retval, "IDAGetLastOrder", 1);
638 retval = IDAGetNumSteps(mem, &nst);
639 check_retval(&retval, "IDAGetNumSteps", 1);
640 retval = IDAGetLastStep(mem, &hused);
641 check_retval(&retval, "IDAGetLastStep", 1);
642
643 c_bl = IJ_Vptr(c,0,0);
644 c_tr = IJ_Vptr(c,MX-1,MY-1);
645
646 #if defined(SUNDIALS_EXTENDED_PRECISION)
647 printf("%8.2Le %12.4Le %12.4Le | %3ld %1d %12.4Le\n",
648 t, c_bl[0], c_tr[0], nst, kused, hused);
649 for (i=1;i<NUM_SPECIES;i++)
650 printf(" %12.4Le %12.4Le |\n",c_bl[i],c_tr[i]);
651 #elif defined(SUNDIALS_DOUBLE_PRECISION)
652 printf("%8.2e %12.4e %12.4e | %3ld %1d %12.4e\n",
653 t, c_bl[0], c_tr[0], nst, kused, hused);
654 for (i=1;i<NUM_SPECIES;i++)
655 printf(" %12.4e %12.4e |\n",c_bl[i],c_tr[i]);
656 #else
657 printf("%8.2e %12.4e %12.4e | %3ld %1d %12.4e\n",
658 t, c_bl[0], c_tr[0], nst, kused, hused);
659 for (i=1;i<NUM_SPECIES;i++)
660 printf(" %12.4e %12.4e |\n",c_bl[i],c_tr[i]);
661 #endif
662
663 printf("\n");
664 }
665
666 /*
667 * PrintFinalStats: Print final run data contained in iopt.
668 */
669
PrintFinalStats(void * mem)670 static void PrintFinalStats(void *mem)
671 {
672 long int nst, nre, sli, netf, nps, npevals, nrevalsLS;
673 int retval;
674
675 retval = IDAGetNumSteps(mem, &nst);
676 check_retval(&retval, "IDAGetNumSteps", 1);
677 retval = IDAGetNumLinIters(mem, &sli);
678 check_retval(&retval, "IDAGetNumLinIters", 1);
679 retval = IDAGetNumResEvals(mem, &nre);
680 check_retval(&retval, "IDAGetNumResEvals", 1);
681 retval = IDAGetNumErrTestFails(mem, &netf);
682 check_retval(&retval, "IDAGetNumErrTestFails", 1);
683 retval = IDAGetNumPrecSolves(mem, &nps);
684 check_retval(&retval, "IDAGetNumPrecSolves", 1);
685 retval = IDAGetNumPrecEvals(mem, &npevals);
686 check_retval(&retval, "IDAGetNumPrecEvals", 1);
687 retval = IDAGetNumLinResEvals(mem, &nrevalsLS);
688 check_retval(&retval, "IDAGetNumLinResEvals", 1);
689
690 printf("-----------------------------------------------------------\n");
691 printf("Final run statistics: \n\n");
692 printf("Number of steps = %ld\n", nst);
693 printf("Number of residual evaluations = %ld\n", nre);
694 printf("Number of Preconditioner evaluations = %ld\n", npevals);
695 printf("Number of linear iterations = %ld\n", sli);
696 printf("Number of error test failures = %ld\n", netf);
697 printf("Number of precond solve fun called = %ld\n", nps);
698
699 }
700
701 /*
702 * Fweb: Rate function for the food-web problem.
703 * This routine computes the right-hand sides of the system equations,
704 * consisting of the diffusion term and interaction term.
705 * The interaction term is computed by the function WebRates.
706 */
707
Fweb(realtype tcalc,N_Vector cc,N_Vector crate,UserData webdata)708 static void Fweb(realtype tcalc, N_Vector cc, N_Vector crate,
709 UserData webdata)
710 {
711 sunindextype jx, jy, is, idyu, idyl, idxu, idxl;
712 realtype xx, yy, *cxy, *ratesxy, *cratexy, dcyli, dcyui, dcxli, dcxui;
713
714 /* Loop over grid points, evaluate interaction vector (length ns),
715 form diffusion difference terms, and load crate. */
716
717 for (jy = 0; jy < MY; jy++) {
718 yy = (webdata->dy) * jy ;
719 idyu = (jy!=MY-1) ? NSMX : -NSMX;
720 idyl = (jy!= 0 ) ? NSMX : -NSMX;
721
722 for (jx = 0; jx < MX; jx++) {
723 xx = (webdata->dx) * jx;
724 idxu = (jx!= MX-1) ? NUM_SPECIES : -NUM_SPECIES;
725 idxl = (jx!= 0 ) ? NUM_SPECIES : -NUM_SPECIES;
726 cxy = IJ_Vptr(cc,jx,jy);
727 ratesxy = IJ_Vptr(webdata->rates,jx,jy);
728 cratexy = IJ_Vptr(crate,jx,jy);
729
730 /* Get interaction vector at this grid point. */
731 WebRates(xx, yy, cxy, ratesxy, webdata);
732
733 /* Loop over species, do differencing, load crate segment. */
734 for (is = 0; is < NUM_SPECIES; is++) {
735
736 /* Differencing in y. */
737 dcyli = *(cxy+is) - *(cxy - idyl + is) ;
738 dcyui = *(cxy + idyu + is) - *(cxy+is);
739
740 /* Differencing in x. */
741 dcxli = *(cxy+is) - *(cxy - idxl + is);
742 dcxui = *(cxy + idxu +is) - *(cxy+is);
743
744 /* Compute the crate values at (xx,yy). */
745 cratexy[is] = coy[is] * (dcyui - dcyli) +
746 cox[is] * (dcxui - dcxli) + ratesxy[is];
747
748 } /* End is loop */
749 } /* End of jx loop */
750 } /* End of jy loop */
751
752 }
753
754 /*
755 * WebRates: Evaluate reaction rates at a given spatial point.
756 * At a given (x,y), evaluate the array of ns reaction terms R.
757 */
758
WebRates(realtype xx,realtype yy,realtype * cxy,realtype * ratesxy,UserData webdata)759 static void WebRates(realtype xx, realtype yy, realtype *cxy, realtype *ratesxy,
760 UserData webdata)
761 {
762 int is;
763 realtype fac;
764
765 for (is = 0; is < NUM_SPECIES; is++)
766 ratesxy[is] = dotprod(NUM_SPECIES, cxy, acoef[is]);
767
768 fac = ONE + ALPHA*xx*yy + BETA*sin(FOURPI*xx)*sin(FOURPI*yy);
769
770 for (is = 0; is < NUM_SPECIES; is++)
771 ratesxy[is] = cxy[is]*( bcoef[is]*fac + ratesxy[is] );
772
773 }
774
775 /*
776 * dotprod: dot product routine for realtype arrays, for use by WebRates.
777 */
778
dotprod(sunindextype size,realtype * x1,realtype * x2)779 static realtype dotprod(sunindextype size, realtype *x1, realtype *x2)
780 {
781 sunindextype i;
782 realtype *xx1, *xx2, temp = ZERO;
783
784 xx1 = x1; xx2 = x2;
785 for (i = 0; i < size; i++) temp += (*xx1++) * (*xx2++);
786 return(temp);
787
788 }
789
790 /*
791 * Check function return value...
792 * opt == 0 means SUNDIALS function allocates memory so check if
793 * returned NULL pointer
794 * opt == 1 means SUNDIALS function returns an integer value so check if
795 * retval < 0
796 * opt == 2 means function allocates memory so check if returned
797 * NULL pointer
798 */
799
check_retval(void * returnvalue,char * funcname,int opt)800 static int check_retval(void *returnvalue, char *funcname, int opt)
801 {
802 int *retval;
803
804 if (opt == 0 && returnvalue == NULL) {
805 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
806 fprintf(stderr,
807 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
808 funcname);
809 return(1);
810 } else if (opt == 1) {
811 /* Check if retval < 0 */
812 retval = (int *) returnvalue;
813 if (*retval < 0) {
814 fprintf(stderr,
815 "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
816 funcname, *retval);
817 return(1);
818 }
819 } else if (opt == 2 && returnvalue == NULL) {
820 /* Check if function returned NULL pointer - no memory allocated */
821 fprintf(stderr,
822 "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
823 funcname);
824 return(1);
825 }
826
827 return(0);
828 }
829