1 /* -----------------------------------------------------------------
2 * Programmer(s): Scott D. Cohen, Alan C. 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 problem:
16 *
17 * An ODE system is generated from the following 2-species diurnal
18 * kinetics advection-diffusion PDE system in 2 space dimensions:
19 *
20 * dc(i)/dt = Kh*(d/dx)^2 c(i) + V*dc(i)/dx + (d/dy)(Kv(y)*dc(i)/dy)
21 * + Ri(c1,c2,t) for i = 1,2, where
22 * R1(c1,c2,t) = -q1*c1*c3 - q2*c1*c2 + 2*q3(t)*c3 + q4(t)*c2 ,
23 * R2(c1,c2,t) = q1*c1*c3 - q2*c1*c2 - q4(t)*c2 ,
24 * Kv(y) = Kv0*exp(y/5) ,
25 * Kh, V, Kv0, q1, q2, and c3 are constants, and q3(t) and q4(t)
26 * vary diurnally. The problem is posed on the square
27 * 0 <= x <= 20, 30 <= y <= 50 (all in km),
28 * with homogeneous Neumann boundary conditions, and for time t in
29 * 0 <= t <= 86400 sec (1 day).
30 * The PDE system is treated by central differences on a uniform
31 * 10 x 10 mesh, with simple polynomial initial profiles.
32 * The problem is solved with CVODE, with the BDF/GMRES
33 * method (i.e. using the SUNLinSol_SPGMR linear solver) and a banded
34 * preconditioner, generated by difference quotients, using the
35 * module CVBANDPRE. The problem is solved with left and right
36 * preconditioning.
37 * -----------------------------------------------------------------*/
38
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <math.h>
42
43 #include <cvode/cvode.h> /* prototypes for CVODE fcts., consts. */
44 #include <nvector/nvector_serial.h> /* access to serial N_Vector */
45 #include <sunlinsol/sunlinsol_spgmr.h> /* access to SPGMR SUNLinearSolver */
46 #include <cvode/cvode_bandpre.h> /* access to CVBANDPRE module */
47 #include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
48
49 /* helpful macros */
50
51 #ifndef SQR
52 #define SQR(A) ((A)*(A))
53 #endif
54
55 /* Problem Constants */
56
57 #define ZERO RCONST(0.0)
58 #define ONE RCONST(1.0)
59 #define TWO RCONST(2.0)
60
61 #define NUM_SPECIES 2 /* number of species */
62 #define KH RCONST(4.0e-6) /* horizontal diffusivity Kh */
63 #define VEL RCONST(0.001) /* advection velocity V */
64 #define KV0 RCONST(1.0e-8) /* coefficient in Kv(y) */
65 #define Q1 RCONST(1.63e-16) /* coefficients q1, q2, c3 */
66 #define Q2 RCONST(4.66e-16)
67 #define C3 RCONST(3.7e16)
68 #define A3 RCONST(22.62) /* coefficient in expression for q3(t) */
69 #define A4 RCONST(7.601) /* coefficient in expression for q4(t) */
70 #define C1_SCALE RCONST(1.0e6) /* coefficients in initial profiles */
71 #define C2_SCALE RCONST(1.0e12)
72
73 #define T0 ZERO /* initial time */
74 #define NOUT 12 /* number of output times */
75 #define TWOHR RCONST(7200.0) /* number of seconds in two hours */
76 #define HALFDAY RCONST(4.32e4) /* number of seconds in a half day */
77 #define PI RCONST(3.1415926535898) /* pi */
78
79 #define XMIN ZERO /* grid boundaries in x */
80 #define XMAX RCONST(20.0)
81 #define YMIN RCONST(30.0) /* grid boundaries in y */
82 #define YMAX RCONST(50.0)
83 #define XMID RCONST(10.0) /* grid midpoints in x,y */
84 #define YMID RCONST(40.0)
85
86 #define MX 10 /* MX = number of x mesh points */
87 #define MY 10 /* MY = number of y mesh points */
88 #define NSMX 20 /* NSMX = NUM_SPECIES*MX */
89 #define MM (MX*MY) /* MM = MX*MY */
90
91 /* CVodeInit Constants */
92
93 #define RTOL RCONST(1.0e-5) /* scalar relative tolerance */
94 #define FLOOR RCONST(100.0) /* value of C1 or C2 at which tolerances */
95 /* change from relative to absolute */
96 #define ATOL (RTOL*FLOOR) /* scalar absolute tolerance */
97 #define NEQ (NUM_SPECIES*MM) /* NEQ = number of equations */
98
99 /* User-defined vector and matrix accessor macros: IJKth, IJth */
100
101 /* IJKth is defined in order to isolate the translation from the
102 mathematical 3-dimensional structure of the dependent variable vector
103 to the underlying 1-dimensional storage. IJth is defined in order to
104 write code which indexes into small dense matrices with a (row,column)
105 pair, where 1 <= row, column <= NUM_SPECIES.
106
107 IJKth(vdata,i,j,k) references the element in the vdata array for
108 species i at mesh point (j,k), where 1 <= i <= NUM_SPECIES,
109 0 <= j <= MX-1, 0 <= k <= MY-1. The vdata array is obtained via
110 the call vdata = N_VGetArrayPointer(v), where v is an N_Vector.
111 For each mesh point (j,k), the elements for species i and i+1 are
112 contiguous within vdata.
113
114 IJth(a,i,j) references the (i,j)th entry of the small matrix realtype **a,
115 where 1 <= i,j <= NUM_SPECIES. The small matrix routines in cvode_bandpre.h
116 work with matrices stored by column in a 2-dimensional array. In C,
117 arrays are indexed starting at 0, not 1. */
118
119 #define IJKth(vdata,i,j,k) (vdata[i-1 + (j)*NUM_SPECIES + (k)*NSMX])
120 #define IJth(a,i,j) (a[j-1][i-1])
121
122 /* Type : UserData
123 contains preconditioner blocks, pivot arrays, and problem constants */
124
125 typedef struct {
126 realtype q4, om, dx, dy, hdco, haco, vdco;
127 } *UserData;
128
129 /* Private Helper Functions */
130
131 static void InitUserData(UserData data);
132 static void SetInitialProfiles(N_Vector u, realtype dx, realtype dy);
133 static void PrintIntro(sunindextype mu, sunindextype ml);
134 static void PrintOutput(void *cvode_mem, N_Vector u, realtype t);
135 static void PrintFinalStats(void *cvode_mem);
136
137 /* Private function to check function return values */
138 static int check_retval(void *returnvalue, const char *funcname, int opt);
139
140 /* Function Called by the Solver */
141
142 static int f(realtype t, N_Vector u, N_Vector udot, void *user_data);
143
144 /*
145 *-------------------------------
146 * Main Program
147 *-------------------------------
148 */
149
main()150 int main()
151 {
152 realtype abstol, reltol, t, tout;
153 N_Vector u;
154 UserData data;
155 SUNLinearSolver LS;
156 void *cvode_mem;
157 int retval, iout, jpre;
158 sunindextype ml, mu;
159
160 u = NULL;
161 data = NULL;
162 LS = NULL;
163 cvode_mem = NULL;
164
165 /* Allocate and initialize u, and set problem data and tolerances */
166 u = N_VNew_Serial(NEQ);
167 if(check_retval((void *)u, "N_VNew_Serial", 0)) return(1);
168 data = (UserData) malloc(sizeof *data);
169 if(check_retval((void *)data, "malloc", 2)) return(1);
170 InitUserData(data);
171 SetInitialProfiles(u, data->dx, data->dy);
172 abstol = ATOL;
173 reltol = RTOL;
174
175 /* Call CVodeCreate to create the solver memory and specify the
176 * Backward Differentiation Formula */
177 cvode_mem = CVodeCreate(CV_BDF);
178 if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
179
180 /* Set the pointer to user-defined data */
181 retval = CVodeSetUserData(cvode_mem, data);
182 if(check_retval(&retval, "CVodeSetUserData", 1)) return(1);
183
184 /* Call CVodeInit to initialize the integrator memory and specify the
185 * user's right hand side function in u'=f(t,u), the inital time T0, and
186 * the initial dependent variable vector u. */
187 retval = CVodeInit(cvode_mem, f, T0, u);
188 if(check_retval(&retval, "CVodeInit", 1)) return(1);
189
190 /* Call CVodeSStolerances to specify the scalar relative tolerance
191 * and scalar absolute tolerances */
192 retval = CVodeSStolerances(cvode_mem, reltol, abstol);
193 if (check_retval(&retval, "CVodeSStolerances", 1)) return(1);
194
195 /* Call SUNLinSol_SPGMR to specify the linear solver SPGMR
196 * with left preconditioning and the default Krylov dimension */
197 LS = SUNLinSol_SPGMR(u, PREC_LEFT, 0);
198 if(check_retval((void *)LS, "SUNLinSol_SPGMR", 0)) return(1);
199
200 /* Call CVodeSetLinearSolver to attach the linear sovler to CVode */
201 retval = CVodeSetLinearSolver(cvode_mem, LS, NULL);
202 if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1;
203
204 /* Call CVBandPreInit to initialize band preconditioner */
205 ml = mu = 2;
206 retval = CVBandPrecInit(cvode_mem, NEQ, mu, ml);
207 if(check_retval(&retval, "CVBandPrecInit", 0)) return(1);
208
209 PrintIntro(mu, ml);
210
211 /* Loop over jpre (= PREC_LEFT, PREC_RIGHT), and solve the problem */
212
213 for (jpre = PREC_LEFT; jpre <= PREC_RIGHT; jpre++) {
214
215 /* On second run, re-initialize u, the solver, and SPGMR */
216
217 if (jpre == PREC_RIGHT) {
218
219 SetInitialProfiles(u, data->dx, data->dy);
220
221 retval = CVodeReInit(cvode_mem, T0, u);
222 if(check_retval(&retval, "CVodeReInit", 1)) return(1);
223
224 retval = SUNLinSol_SPGMRSetPrecType(LS, PREC_RIGHT);
225 if(check_retval(&retval, "SUNLinSol_SPGMRSetPrecType", 1)) return(1);
226
227 retval = CVBandPrecInit(cvode_mem, NEQ, mu, ml);
228 if(check_retval(&retval, "CVBandPrecInit", 0)) return(1);
229
230 printf("\n\n-------------------------------------------------------");
231 printf("------------\n");
232 }
233
234 printf("\n\nPreconditioner type is: jpre = %s\n\n",
235 (jpre == PREC_LEFT) ? "PREC_LEFT" : "PREC_RIGHT");
236
237 /* In loop over output points, call CVode, print results, test for error */
238
239 for (iout = 1, tout = TWOHR; iout <= NOUT; iout++, tout += TWOHR) {
240 retval = CVode(cvode_mem, tout, u, &t, CV_NORMAL);
241 check_retval(&retval, "CVode", 1);
242 PrintOutput(cvode_mem, u, t);
243 if (retval != CV_SUCCESS) {
244 break;
245 }
246 }
247
248 /* Print final statistics */
249
250 PrintFinalStats(cvode_mem);
251
252 } /* End of jpre loop */
253
254 /* Free memory */
255 N_VDestroy(u);
256 free(data);
257 CVodeFree(&cvode_mem);
258 SUNLinSolFree(LS);
259
260 return(0);
261 }
262
263 /*
264 *-------------------------------
265 * Private helper functions
266 *-------------------------------
267 */
268
269 /* Load problem constants in data */
270
InitUserData(UserData data)271 static void InitUserData(UserData data)
272 {
273 data->om = PI/HALFDAY;
274 data->dx = (XMAX-XMIN)/(MX-1);
275 data->dy = (YMAX-YMIN)/(MY-1);
276 data->hdco = KH/SQR(data->dx);
277 data->haco = VEL/(TWO*data->dx);
278 data->vdco = (ONE/SQR(data->dy))*KV0;
279 }
280
281 /* Set initial conditions in u */
282
SetInitialProfiles(N_Vector u,realtype dx,realtype dy)283 static void SetInitialProfiles(N_Vector u, realtype dx, realtype dy)
284 {
285 int jx, jy;
286 realtype x, y, cx, cy;
287 realtype *udata;
288
289 /* Set pointer to data array in vector u. */
290
291 udata = N_VGetArrayPointer(u);
292
293 /* Load initial profiles of c1 and c2 into u vector */
294
295 for (jy=0; jy < MY; jy++) {
296 y = YMIN + jy*dy;
297 cy = SQR(RCONST(0.1)*(y - YMID));
298 cy = ONE - cy + RCONST(0.5)*SQR(cy);
299 for (jx=0; jx < MX; jx++) {
300 x = XMIN + jx*dx;
301 cx = SQR(RCONST(0.1)*(x - XMID));
302 cx = ONE - cx + RCONST(0.5)*SQR(cx);
303 IJKth(udata,1,jx,jy) = C1_SCALE*cx*cy;
304 IJKth(udata,2,jx,jy) = C2_SCALE*cx*cy;
305 }
306 }
307 }
308
PrintIntro(sunindextype mu,sunindextype ml)309 static void PrintIntro(sunindextype mu, sunindextype ml)
310 {
311 printf("2-species diurnal advection-diffusion problem, %d by %d mesh\n",
312 MX, MY);
313 printf("SPGMR solver; band preconditioner; mu = %ld, ml = %ld\n\n",
314 (long int) mu, (long int) ml);
315
316 return;
317 }
318
319 /* Print current t, step count, order, stepsize, and sampled c1,c2 values */
320
PrintOutput(void * cvode_mem,N_Vector u,realtype t)321 static void PrintOutput(void *cvode_mem, N_Vector u, realtype t)
322 {
323 long int nst;
324 int qu, retval;
325 realtype hu, *udata;
326 int mxh = MX/2 - 1, myh = MY/2 - 1, mx1 = MX - 1, my1 = MY - 1;
327
328 udata = N_VGetArrayPointer(u);
329
330 retval = CVodeGetNumSteps(cvode_mem, &nst);
331 check_retval(&retval, "CVodeGetNumSteps", 1);
332 retval = CVodeGetLastOrder(cvode_mem, &qu);
333 check_retval(&retval, "CVodeGetLastOrder", 1);
334 retval = CVodeGetLastStep(cvode_mem, &hu);
335 check_retval(&retval, "CVodeGetLastStep", 1);
336
337 #if defined(SUNDIALS_EXTENDED_PRECISION)
338 printf("t = %.2Le no. steps = %ld order = %d stepsize = %.2Le\n",
339 t, nst, qu, hu);
340 printf("c1 (bot.left/middle/top rt.) = %12.3Le %12.3Le %12.3Le\n",
341 IJKth(udata,1,0,0), IJKth(udata,1,mxh,myh), IJKth(udata,1,mx1,my1));
342 printf("c2 (bot.left/middle/top rt.) = %12.3Le %12.3Le %12.3Le\n\n",
343 IJKth(udata,2,0,0), IJKth(udata,2,mxh,myh), IJKth(udata,2,mx1,my1));
344 #elif defined(SUNDIALS_DOUBLE_PRECISION)
345 printf("t = %.2e no. steps = %ld order = %d stepsize = %.2e\n",
346 t, nst, qu, hu);
347 printf("c1 (bot.left/middle/top rt.) = %12.3e %12.3e %12.3e\n",
348 IJKth(udata,1,0,0), IJKth(udata,1,mxh,myh), IJKth(udata,1,mx1,my1));
349 printf("c2 (bot.left/middle/top rt.) = %12.3e %12.3e %12.3e\n\n",
350 IJKth(udata,2,0,0), IJKth(udata,2,mxh,myh), IJKth(udata,2,mx1,my1));
351 #else
352 printf("t = %.2e no. steps = %ld order = %d stepsize = %.2e\n",
353 t, nst, qu, hu);
354 printf("c1 (bot.left/middle/top rt.) = %12.3e %12.3e %12.3e\n",
355 IJKth(udata,1,0,0), IJKth(udata,1,mxh,myh), IJKth(udata,1,mx1,my1));
356 printf("c2 (bot.left/middle/top rt.) = %12.3e %12.3e %12.3e\n\n",
357 IJKth(udata,2,0,0), IJKth(udata,2,mxh,myh), IJKth(udata,2,mx1,my1));
358 #endif
359 }
360
361 /* Get and print final statistics */
362
PrintFinalStats(void * cvode_mem)363 static void PrintFinalStats(void *cvode_mem)
364 {
365 long int lenrw, leniw ;
366 long int lenrwLS, leniwLS;
367 long int lenrwBP, leniwBP;
368 long int nst, nfe, nsetups, nni, ncfn, netf;
369 long int nli, npe, nps, ncfl, nfeLS;
370 long int nfeBP;
371 int retval;
372
373 retval = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw);
374 check_retval(&retval, "CVodeGetWorkSpace", 1);
375 retval = CVodeGetNumSteps(cvode_mem, &nst);
376 check_retval(&retval, "CVodeGetNumSteps", 1);
377 retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);
378 check_retval(&retval, "CVodeGetNumRhsEvals", 1);
379 retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
380 check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);
381 retval = CVodeGetNumErrTestFails(cvode_mem, &netf);
382 check_retval(&retval, "CVodeGetNumErrTestFails", 1);
383 retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
384 check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);
385 retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
386 check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);
387
388 retval = CVodeGetLinWorkSpace(cvode_mem, &lenrwLS, &leniwLS);
389 check_retval(&retval, "CVodeGetLinWorkSpace", 1);
390 retval = CVodeGetNumLinIters(cvode_mem, &nli);
391 check_retval(&retval, "CVodeGetNumLinIters", 1);
392 retval = CVodeGetNumPrecEvals(cvode_mem, &npe);
393 check_retval(&retval, "CVodeGetNumPrecEvals", 1);
394 retval = CVodeGetNumPrecSolves(cvode_mem, &nps);
395 check_retval(&retval, "CVodeGetNumPrecSolves", 1);
396 retval = CVodeGetNumLinConvFails(cvode_mem, &ncfl);
397 check_retval(&retval, "CVodeGetNumLinConvFails", 1);
398 retval = CVodeGetNumLinRhsEvals(cvode_mem, &nfeLS);
399 check_retval(&retval, "CVodeGetNumLinRhsEvals", 1);
400
401 retval = CVBandPrecGetWorkSpace(cvode_mem, &lenrwBP, &leniwBP);
402 check_retval(&retval, "CVBandPrecGetWorkSpace", 1);
403 retval = CVBandPrecGetNumRhsEvals(cvode_mem, &nfeBP);
404 check_retval(&retval, "CVBandPrecGetNumRhsEvals", 1);
405
406 printf("\nFinal Statistics.. \n\n");
407 printf("lenrw = %5ld leniw = %5ld\n" , lenrw, leniw);
408 printf("lenrwls = %5ld leniwls = %5ld\n" , lenrwLS, leniwLS);
409 printf("lenrwbp = %5ld leniwbp = %5ld\n" , lenrwBP, leniwBP);
410 printf("nst = %5ld\n" , nst);
411 printf("nfe = %5ld nfetot = %5ld\n" , nfe, nfe+nfeLS+nfeBP);
412 printf("nfeLS = %5ld nfeBP = %5ld\n" , nfeLS, nfeBP);
413 printf("nni = %5ld nli = %5ld\n" , nni, nli);
414 printf("nsetups = %5ld netf = %5ld\n" , nsetups, netf);
415 printf("npe = %5ld nps = %5ld\n" , npe, nps);
416 printf("ncfn = %5ld ncfl = %5ld\n\n", ncfn, ncfl);
417 }
418
419 /* Check function return value...
420 opt == 0 means SUNDIALS function allocates memory so check if
421 returned NULL pointer
422 opt == 1 means SUNDIALS function returns an integer value so check if
423 retval < 0
424 opt == 2 means function allocates memory so check if returned
425 NULL pointer */
426
check_retval(void * returnvalue,const char * funcname,int opt)427 static int check_retval(void *returnvalue, const char *funcname, int opt)
428 {
429 int *retval;
430
431 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
432 if (opt == 0 && returnvalue == NULL) {
433 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
434 funcname);
435 return(1); }
436
437 /* Check if retval < 0 */
438 else if (opt == 1) {
439 retval = (int *) returnvalue;
440 if (*retval < 0) {
441 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
442 funcname, *retval);
443 return(1); }}
444
445 /* Check if function returned NULL pointer - no memory allocated */
446 else if (opt == 2 && returnvalue == NULL) {
447 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
448 funcname);
449 return(1); }
450
451 return(0);
452 }
453
454 /*
455 *-------------------------------
456 * Function called by the solver
457 *-------------------------------
458 */
459
460 /* f routine. Compute RHS function f(t,u). */
461
f(realtype t,N_Vector u,N_Vector udot,void * user_data)462 static int f(realtype t, N_Vector u, N_Vector udot, void *user_data)
463 {
464 realtype q3, c1, c2, c1dn, c2dn, c1up, c2up, c1lt, c2lt;
465 realtype c1rt, c2rt, cydn, cyup, hord1, hord2, horad1, horad2;
466 realtype qq1, qq2, qq3, qq4, rkin1, rkin2, s, vertd1, vertd2, ydn, yup;
467 realtype q4coef, dely, verdco, hordco, horaco;
468 realtype *udata, *dudata;
469 int jx, jy, idn, iup, ileft, iright;
470 UserData data;
471
472 data = (UserData) user_data;
473 udata = N_VGetArrayPointer(u);
474 dudata = N_VGetArrayPointer(udot);
475
476 /* Set diurnal rate coefficients. */
477
478 s = sin(data->om*t);
479 if (s > ZERO) {
480 q3 = exp(-A3/s);
481 data->q4 = exp(-A4/s);
482 } else {
483 q3 = ZERO;
484 data->q4 = ZERO;
485 }
486
487 /* Make local copies of problem variables, for efficiency. */
488
489 q4coef = data->q4;
490 dely = data->dy;
491 verdco = data->vdco;
492 hordco = data->hdco;
493 horaco = data->haco;
494
495 /* Loop over all grid points. */
496
497 for (jy=0; jy < MY; jy++) {
498
499 /* Set vertical diffusion coefficients at jy +- 1/2 */
500
501 ydn = YMIN + (jy - RCONST(0.5))*dely;
502 yup = ydn + dely;
503 cydn = verdco*exp(RCONST(0.2)*ydn);
504 cyup = verdco*exp(RCONST(0.2)*yup);
505 idn = (jy == 0) ? 1 : -1;
506 iup = (jy == MY-1) ? -1 : 1;
507 for (jx=0; jx < MX; jx++) {
508
509 /* Extract c1 and c2, and set kinetic rate terms. */
510
511 c1 = IJKth(udata,1,jx,jy);
512 c2 = IJKth(udata,2,jx,jy);
513 qq1 = Q1*c1*C3;
514 qq2 = Q2*c1*c2;
515 qq3 = q3*C3;
516 qq4 = q4coef*c2;
517 rkin1 = -qq1 - qq2 + TWO*qq3 + qq4;
518 rkin2 = qq1 - qq2 - qq4;
519
520 /* Set vertical diffusion terms. */
521
522 c1dn = IJKth(udata,1,jx,jy+idn);
523 c2dn = IJKth(udata,2,jx,jy+idn);
524 c1up = IJKth(udata,1,jx,jy+iup);
525 c2up = IJKth(udata,2,jx,jy+iup);
526 vertd1 = cyup*(c1up - c1) - cydn*(c1 - c1dn);
527 vertd2 = cyup*(c2up - c2) - cydn*(c2 - c2dn);
528
529 /* Set horizontal diffusion and advection terms. */
530
531 ileft = (jx == 0) ? 1 : -1;
532 iright =(jx == MX-1) ? -1 : 1;
533 c1lt = IJKth(udata,1,jx+ileft,jy);
534 c2lt = IJKth(udata,2,jx+ileft,jy);
535 c1rt = IJKth(udata,1,jx+iright,jy);
536 c2rt = IJKth(udata,2,jx+iright,jy);
537 hord1 = hordco*(c1rt - TWO*c1 + c1lt);
538 hord2 = hordco*(c2rt - TWO*c2 + c2lt);
539 horad1 = horaco*(c1rt - c1lt);
540 horad2 = horaco*(c2rt - c2lt);
541
542 /* Load all terms into udot. */
543
544 IJKth(dudata, 1, jx, jy) = vertd1 + hord1 + horad1 + rkin1;
545 IJKth(dudata, 2, jx, jy) = vertd2 + hord2 + horad2 + rkin2;
546 }
547 }
548
549 return(0);
550 }
551