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