Actual source code: ex5.c

  1: /*
  2:   Note:
  3:   To performance a convergence study, a reference solution can be generated with a small stepsize and stored into a binary file, e.g. -ts_type ssp -ts_dt 1e-5 -ts_max_steps 30000 -ts_view_solution binary:reference.bin
  4:   Errors can be computed in the following runs with -simulation -f reference.bin

  6:   Multirate RK options:
  7:   -ts_rk_dtratio is the ratio between larger time step size and small time step size
  8:   -ts_rk_multirate_type has three choices: none (for single step RK)
  9:                                            combined (for multirate method and user just need to provide the same RHS with the single step RK)
 10:                                            partitioned (for multiraet method and user need to provide two RHS, one is for fast components and the orther is for slow components
 11: */

 13: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
 14:   " advection   - Variable coefficient scalar advection\n"
 15:   "                u_t       + (a*u)_x               = 0\n"
 16:   " a is a piecewise function with two values say a0 and a1 (both a0 and a1 are constant).\n"
 17:   " in this toy problem we assume a=a0 when x<0 and a=a1 when x>0 with a0<a1 which has the same discontinuous point with initial condition.\n"
 18:   " we don't provide an exact solution, so you need to generate reference solution to do convergnece staudy,\n"
 19:   " more precisely, you need firstly generate a reference to a binary file say file.bin, then on commend line,\n"
 20:   " you should type -simulation -f file.bin.\n"
 21:   " you can choose the number of grids by -da_grid_x.\n"
 22:   " you can choose the value of a by -physics_advect_a1 and -physics_advect_a2.\n";

 24: #include <petscts.h>
 25: #include <petscdm.h>
 26: #include <petscdmda.h>
 27: #include <petscdraw.h>
 28: #include <petsc/private/tsimpl.h>

 30: #include <petsc/private/kernels/blockinvert.h>

 32: #include "finitevolume1d.h"

 34: static inline PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }

 36: /* --------------------------------- Advection ----------------------------------- */

 38: typedef struct {
 39:   PetscReal a[2];                  /* advective velocity */
 40: } AdvectCtx;

 42: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed,PetscReal x,PetscReal xmin, PetscReal xmax)
 43: {
 44:   AdvectCtx *ctx = (AdvectCtx*)vctx;
 45:   PetscReal *speed;

 48:   speed = ctx->a;
 49:   if (x==0 || x == xmin || x == xmax) flux[0] = PetscMax(0,(speed[0]+speed[1])/2.0)*uL[0] + PetscMin(0,(speed[0]+speed[1])/2.0)*uR[0]; /* if condition need to be changed base on different problem, '0' is the discontinuous point of a */
 50:   else if (x<0) flux[0] = PetscMax(0,speed[0])*uL[0] + PetscMin(0,speed[0])*uR[0];  /* else if condition need to be changed based on diferent problem, 'x = 0' is discontinuous point of a */
 51:   else flux[0] = PetscMax(0,speed[1])*uL[0] + PetscMin(0,speed[1])*uR[0];
 52:   *maxspeed = *speed;
 53:   return 0;
 54: }

 56: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds,PetscReal x)
 57: {
 58:   AdvectCtx *ctx = (AdvectCtx*)vctx;

 61:   X[0]      = 1.;
 62:   Xi[0]     = 1.;
 63:   if (x<0) speeds[0] = ctx->a[0];  /* x is discontinuous point of a */
 64:   else    speeds[0] = ctx->a[1];
 65:   return 0;
 66: }

 68: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
 69: {
 70:   AdvectCtx *ctx = (AdvectCtx*)vctx;
 71:   PetscReal *a    = ctx->a,x0;

 74:   if (x<0){   /* x is cell center */
 75:     switch (bctype) {
 76:       case FVBC_OUTFLOW:  x0 = x-a[0]*t; break;
 77:       case FVBC_PERIODIC: x0 = RangeMod(x-a[0]*t,xmin,xmax); break;
 78:       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
 79:     }
 80:     switch (initial) {
 81:       case 0: u[0] = (x0 < 0) ? 1 : -1; break;
 82:       case 1: u[0] = (x0 < 0) ? -1 : 1; break;
 83:       case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
 84:       case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
 85:       case 4: u[0] = PetscAbs(x0); break;
 86:       case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
 87:       case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
 88:       case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
 89:       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
 90:     }
 91:   }
 92:   else{
 93:     switch (bctype) {
 94:       case FVBC_OUTFLOW:  x0 = x-a[1]*t; break;
 95:       case FVBC_PERIODIC: x0 = RangeMod(x-a[1]*t,xmin,xmax); break;
 96:       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
 97:     }
 98:     switch (initial) {
 99:       case 0: u[0] = (x0 < 0) ? 1 : -1; break;
100:       case 1: u[0] = (x0 < 0) ? -1 : 1; break;
101:       case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
102:       case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
103:       case 4: u[0] = PetscAbs(x0); break;
104:       case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
105:       case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
106:       case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
107:       default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
108:     }
109:   }
110:   return 0;
111: }

113: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
114: {
116:   AdvectCtx      *user;

119:   PetscNew(&user);
120:   ctx->physics.sample         = PhysicsSample_Advect;
121:   ctx->physics.riemann        = PhysicsRiemann_Advect;
122:   ctx->physics.characteristic = PhysicsCharacteristic_Advect;
123:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
124:   ctx->physics.user           = user;
125:   ctx->physics.dof            = 1;
126:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
127:   user->a[0] = 1;
128:   user->a[1] = 1;
129:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
130:   {
131:     PetscOptionsReal("-physics_advect_a1","Speed1","",user->a[0],&user->a[0],NULL);
132:     PetscOptionsReal("-physics_advect_a2","Speed2","",user->a[1],&user->a[1],NULL);
133:   }
134:   PetscOptionsEnd();
135:   return 0;
136: }

138: /* --------------------------------- Finite Volume Solver for slow parts ----------------------------------- */

140: PetscErrorCode FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
141: {
142:   FVCtx          *ctx = (FVCtx*)vctx;
143:   PetscInt       i,j,k,Mx,dof,xs,xm,len_slow;
144:   PetscReal      hx,cfl_idt = 0;
145:   PetscScalar    *x,*f,*slope;
146:   Vec            Xloc;
147:   DM             da;

150:   TSGetDM(ts,&da);
151:   DMGetLocalVector(da,&Xloc);
152:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
153:   hx   = (ctx->xmax-ctx->xmin)/Mx;
154:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
155:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
156:   VecZeroEntries(F);
157:   DMDAVecGetArray(da,Xloc,&x);
158:   VecGetArray(F,&f);
159:   DMDAGetArray(da,PETSC_TRUE,&slope);
160:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
161:   ISGetSize(ctx->iss,&len_slow);

163:   if (ctx->bctype == FVBC_OUTFLOW) {
164:     for (i=xs-2; i<0; i++) {
165:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
166:     }
167:     for (i=Mx; i<xs+xm+2; i++) {
168:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
169:     }
170:   }
171:   for (i=xs-1; i<xs+xm+1; i++) {
172:     struct _LimitInfo info;
173:     PetscScalar       *cjmpL,*cjmpR;
174:     if (i < len_slow+1) {
175:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
176:       (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);
177:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
178:       PetscArrayzero(ctx->cjmpLR,2*dof);
179:       cjmpL = &ctx->cjmpLR[0];
180:       cjmpR = &ctx->cjmpLR[dof];
181:       for (j=0; j<dof; j++) {
182:         PetscScalar jmpL,jmpR;
183:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
184:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
185:         for (k=0; k<dof; k++) {
186:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
187:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
188:         }
189:       }
190:       /* Apply limiter to the left and right characteristic jumps */
191:       info.m  = dof;
192:       info.hx = hx;
193:       (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
194:       for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
195:       for (j=0; j<dof; j++) {
196:         PetscScalar tmp = 0;
197:         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
198:         slope[i*dof+j] = tmp;
199:       }
200:     }
201:   }

203:   for (i=xs; i<xs+xm+1; i++) {
204:     PetscReal   maxspeed;
205:     PetscScalar *uL,*uR;
206:     if (i < len_slow) { /* slow parts can be changed based on a */
207:       uL = &ctx->uLR[0];
208:       uR = &ctx->uLR[dof];
209:       for (j=0; j<dof; j++) {
210:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
211:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
212:       }
213:       (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
214:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
215:       if (i > xs) {
216:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
217:       }
218:       if (i < xs+xm) {
219:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
220:       }
221:     }
222:     if (i == len_slow) { /* slow parts can be changed based on a */
223:       uL = &ctx->uLR[0];
224:       uR = &ctx->uLR[dof];
225:       for (j=0; j<dof; j++) {
226:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
227:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
228:       }
229:       (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
230:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
231:       if (i > xs) {
232:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
233:       }
234:     }
235:   }
236:   DMDAVecRestoreArray(da,Xloc,&x);
237:   VecRestoreArray(F,&f);
238:   DMDARestoreArray(da,PETSC_TRUE,&slope);
239:   DMRestoreLocalVector(da,&Xloc);
240:   return 0;
241: }

243: /* --------------------------------- Finite Volume Solver for fast  parts ----------------------------------- */

245: PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
246: {
247:   FVCtx          *ctx = (FVCtx*)vctx;
248:   PetscInt       i,j,k,Mx,dof,xs,xm,len_slow;
249:   PetscReal      hx,cfl_idt = 0;
250:   PetscScalar    *x,*f,*slope;
251:   Vec            Xloc;
252:   DM             da;

255:   TSGetDM(ts,&da);
256:   DMGetLocalVector(da,&Xloc);
257:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
258:   hx   = (ctx->xmax-ctx->xmin)/Mx;
259:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
260:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
261:   VecZeroEntries(F);
262:   DMDAVecGetArray(da,Xloc,&x);
263:   VecGetArray(F,&f);
264:   DMDAGetArray(da,PETSC_TRUE,&slope);
265:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
266:   ISGetSize(ctx->iss,&len_slow);

268:   if (ctx->bctype == FVBC_OUTFLOW) {
269:     for (i=xs-2; i<0; i++) {
270:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
271:     }
272:     for (i=Mx; i<xs+xm+2; i++) {
273:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
274:     }
275:   }
276:   for (i=xs-1; i<xs+xm+1; i++) {
277:     struct _LimitInfo info;
278:     PetscScalar       *cjmpL,*cjmpR;
279:     if (i > len_slow-2) {
280:       (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);
281:       PetscArrayzero(ctx->cjmpLR,2*dof);
282:       cjmpL = &ctx->cjmpLR[0];
283:       cjmpR = &ctx->cjmpLR[dof];
284:       for (j=0; j<dof; j++) {
285:         PetscScalar jmpL,jmpR;
286:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
287:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
288:         for (k=0; k<dof; k++) {
289:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
290:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
291:         }
292:       }
293:       /* Apply limiter to the left and right characteristic jumps */
294:       info.m  = dof;
295:       info.hx = hx;
296:       (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
297:       for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
298:       for (j=0; j<dof; j++) {
299:         PetscScalar tmp = 0;
300:         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
301:         slope[i*dof+j] = tmp;
302:       }
303:     }
304:   }

306:   for (i=xs; i<xs+xm+1; i++) {
307:     PetscReal   maxspeed;
308:     PetscScalar *uL,*uR;
309:     if (i > len_slow) {
310:       uL = &ctx->uLR[0];
311:       uR = &ctx->uLR[dof];
312:       for (j=0; j<dof; j++) {
313:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
314:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
315:       }
316:       (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
317:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
318:       if (i > xs) {
319:         for (j=0; j<dof; j++) f[(i-len_slow-1)*dof+j] -= ctx->flux[j]/hx;
320:       }
321:       if (i < xs+xm) {
322:         for (j=0; j<dof; j++) f[(i-len_slow)*dof+j] += ctx->flux[j]/hx;
323:       }
324:     }
325:     if (i == len_slow) {
326:       uL = &ctx->uLR[0];
327:       uR = &ctx->uLR[dof];
328:       for (j=0; j<dof; j++) {
329:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
330:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
331:       }
332:       (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
333:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
334:       if (i < xs+xm) {
335:         for (j=0; j<dof; j++) f[(i-len_slow)*dof+j] += ctx->flux[j]/hx;
336:       }
337:     }
338:   }
339:   DMDAVecRestoreArray(da,Xloc,&x);
340:   VecRestoreArray(F,&f);
341:   DMDARestoreArray(da,PETSC_TRUE,&slope);
342:   DMRestoreLocalVector(da,&Xloc);
343:   return 0;
344: }

346: PetscErrorCode FVRHSFunctionslow2(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
347: {
348:   FVCtx          *ctx = (FVCtx*)vctx;
349:   PetscInt       i,j,k,Mx,dof,xs,xm,len_slow1,len_slow2;
350:   PetscReal      hx,cfl_idt = 0;
351:   PetscScalar    *x,*f,*slope;
352:   Vec            Xloc;
353:   DM             da;

356:   TSGetDM(ts,&da);
357:   DMGetLocalVector(da,&Xloc);
358:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
359:   hx   = (ctx->xmax-ctx->xmin)/Mx;
360:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
361:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
362:   VecZeroEntries(F);
363:   DMDAVecGetArray(da,Xloc,&x);
364:   VecGetArray(F,&f);
365:   DMDAGetArray(da,PETSC_TRUE,&slope);
366:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
367:   ISGetSize(ctx->iss,&len_slow1);
368:   ISGetSize(ctx->iss2,&len_slow2);
369:   if (ctx->bctype == FVBC_OUTFLOW) {
370:     for (i=xs-2; i<0; i++) {
371:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
372:     }
373:     for (i=Mx; i<xs+xm+2; i++) {
374:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
375:     }
376:   }
377:   for (i=xs-1; i<xs+xm+1; i++) {
378:     struct _LimitInfo info;
379:     PetscScalar       *cjmpL,*cjmpR;
380:     if (i < len_slow1+len_slow2+1 && i > len_slow1-2) {
381:       (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);
382:       PetscArrayzero(ctx->cjmpLR,2*dof);
383:       cjmpL = &ctx->cjmpLR[0];
384:       cjmpR = &ctx->cjmpLR[dof];
385:       for (j=0; j<dof; j++) {
386:         PetscScalar jmpL,jmpR;
387:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
388:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
389:         for (k=0; k<dof; k++) {
390:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
391:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
392:         }
393:       }
394:       /* Apply limiter to the left and right characteristic jumps */
395:       info.m  = dof;
396:       info.hx = hx;
397:       (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
398:       for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
399:       for (j=0; j<dof; j++) {
400:         PetscScalar tmp = 0;
401:         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
402:         slope[i*dof+j] = tmp;
403:       }
404:     }
405:   }

407:   for (i=xs; i<xs+xm+1; i++) {
408:     PetscReal   maxspeed;
409:     PetscScalar *uL,*uR;
410:     if (i < len_slow1+len_slow2 && i > len_slow1) { /* slow parts can be changed based on a */
411:       uL = &ctx->uLR[0];
412:       uR = &ctx->uLR[dof];
413:       for (j=0; j<dof; j++) {
414:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
415:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
416:       }
417:       (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
418:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
419:       if (i > xs) {
420:         for (j=0; j<dof; j++) f[(i-len_slow1-1)*dof+j] -= ctx->flux[j]/hx;
421:       }
422:       if (i < xs+xm) {
423:         for (j=0; j<dof; j++) f[(i-len_slow1)*dof+j] += ctx->flux[j]/hx;
424:       }
425:     }
426:     if (i == len_slow1+len_slow2) { /* slow parts can be changed based on a */
427:       uL = &ctx->uLR[0];
428:       uR = &ctx->uLR[dof];
429:       for (j=0; j<dof; j++) {
430:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
431:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
432:       }
433:       (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
434:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
435:       if (i > xs) {
436:         for (j=0; j<dof; j++) f[(i-len_slow1-1)*dof+j] -= ctx->flux[j]/hx;
437:       }
438:     }
439:     if (i == len_slow1) { /* slow parts can be changed based on a */
440:       uL = &ctx->uLR[0];
441:       uR = &ctx->uLR[dof];
442:       for (j=0; j<dof; j++) {
443:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
444:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
445:       }
446:       (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
447:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
448:       if (i < xs+xm) {
449:         for (j=0; j<dof; j++) f[(i-len_slow1)*dof+j] += ctx->flux[j]/hx;
450:       }
451:     }
452:   }

454:   DMDAVecRestoreArray(da,Xloc,&x);
455:   VecRestoreArray(F,&f);
456:   DMDARestoreArray(da,PETSC_TRUE,&slope);
457:   DMRestoreLocalVector(da,&Xloc);
458:   return 0;
459: }

461: PetscErrorCode FVRHSFunctionfast2(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
462: {
463:   FVCtx          *ctx = (FVCtx*)vctx;
464:   PetscInt       i,j,k,Mx,dof,xs,xm,len_slow1,len_slow2;
465:   PetscReal      hx,cfl_idt = 0;
466:   PetscScalar    *x,*f,*slope;
467:   Vec            Xloc;
468:   DM             da;

471:   TSGetDM(ts,&da);
472:   DMGetLocalVector(da,&Xloc);
473:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
474:   hx   = (ctx->xmax-ctx->xmin)/Mx;
475:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
476:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
477:   VecZeroEntries(F);
478:   DMDAVecGetArray(da,Xloc,&x);
479:   VecGetArray(F,&f);
480:   DMDAGetArray(da,PETSC_TRUE,&slope);
481:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
482:   ISGetSize(ctx->iss,&len_slow1);
483:   ISGetSize(ctx->iss2,&len_slow2);

485:   if (ctx->bctype == FVBC_OUTFLOW) {
486:     for (i=xs-2; i<0; i++) {
487:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
488:     }
489:     for (i=Mx; i<xs+xm+2; i++) {
490:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
491:     }
492:   }
493:   for (i=xs-1; i<xs+xm+1; i++) {
494:     struct _LimitInfo info;
495:     PetscScalar       *cjmpL,*cjmpR;
496:     if (i > len_slow1+len_slow2-2) {
497:       (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);
498:       PetscArrayzero(ctx->cjmpLR,2*dof);
499:       cjmpL = &ctx->cjmpLR[0];
500:       cjmpR = &ctx->cjmpLR[dof];
501:       for (j=0; j<dof; j++) {
502:         PetscScalar jmpL,jmpR;
503:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
504:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
505:         for (k=0; k<dof; k++) {
506:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
507:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
508:         }
509:       }
510:       /* Apply limiter to the left and right characteristic jumps */
511:       info.m  = dof;
512:       info.hx = hx;
513:       (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
514:       for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
515:       for (j=0; j<dof; j++) {
516:         PetscScalar tmp = 0;
517:         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
518:         slope[i*dof+j] = tmp;
519:       }
520:     }
521:   }

523:   for (i=xs; i<xs+xm+1; i++) {
524:     PetscReal   maxspeed;
525:     PetscScalar *uL,*uR;
526:     if (i > len_slow1+len_slow2) {
527:       uL = &ctx->uLR[0];
528:       uR = &ctx->uLR[dof];
529:       for (j=0; j<dof; j++) {
530:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
531:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
532:       }
533:       (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
534:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
535:       if (i > xs) {
536:         for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2-1)*dof+j] -= ctx->flux[j]/hx;
537:       }
538:       if (i < xs+xm) {
539:         for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2)*dof+j] += ctx->flux[j]/hx;
540:       }
541:     }
542:     if (i == len_slow1+len_slow2) {
543:       uL = &ctx->uLR[0];
544:       uR = &ctx->uLR[dof];
545:       for (j=0; j<dof; j++) {
546:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
547:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
548:       }
549:       (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);
550:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
551:       if (i < xs+xm) {
552:         for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2)*dof+j] += ctx->flux[j]/hx;
553:       }
554:     }
555:   }
556:   DMDAVecRestoreArray(da,Xloc,&x);
557:   VecRestoreArray(F,&f);
558:   DMDARestoreArray(da,PETSC_TRUE,&slope);
559:   DMRestoreLocalVector(da,&Xloc);
560:   return 0;
561: }

563: int main(int argc,char *argv[])
564: {
565:   char              lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
566:   PetscFunctionList limiters   = 0,physics = 0;
567:   MPI_Comm          comm;
568:   TS                ts;
569:   DM                da;
570:   Vec               X,X0,R;
571:   FVCtx             ctx;
572:   PetscInt          i,k,dof,xs,xm,Mx,draw = 0,*index_slow,*index_fast,islow = 0,ifast = 0;
573:   PetscBool         view_final = PETSC_FALSE;
574:   PetscReal         ptime;
575:   PetscErrorCode    ierr;

577:   PetscInitialize(&argc,&argv,0,help);
578:   comm = PETSC_COMM_WORLD;
579:   PetscMemzero(&ctx,sizeof(ctx));

581:   /* Register limiters to be available on the command line */
582:   PetscFunctionListAdd(&limiters,"upwind"              ,Limit_Upwind);
583:   PetscFunctionListAdd(&limiters,"lax-wendroff"        ,Limit_LaxWendroff);
584:   PetscFunctionListAdd(&limiters,"beam-warming"        ,Limit_BeamWarming);
585:   PetscFunctionListAdd(&limiters,"fromm"               ,Limit_Fromm);
586:   PetscFunctionListAdd(&limiters,"minmod"              ,Limit_Minmod);
587:   PetscFunctionListAdd(&limiters,"superbee"            ,Limit_Superbee);
588:   PetscFunctionListAdd(&limiters,"mc"                  ,Limit_MC);
589:   PetscFunctionListAdd(&limiters,"vanleer"             ,Limit_VanLeer);
590:   PetscFunctionListAdd(&limiters,"vanalbada"           ,Limit_VanAlbada);
591:   PetscFunctionListAdd(&limiters,"vanalbadatvd"        ,Limit_VanAlbadaTVD);
592:   PetscFunctionListAdd(&limiters,"koren"               ,Limit_Koren);
593:   PetscFunctionListAdd(&limiters,"korensym"            ,Limit_KorenSym);
594:   PetscFunctionListAdd(&limiters,"koren3"              ,Limit_Koren3);
595:   PetscFunctionListAdd(&limiters,"cada-torrilhon2"     ,Limit_CadaTorrilhon2);
596:   PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1);
597:   PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1"  ,Limit_CadaTorrilhon3R1);
598:   PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10);
599:   PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100);

601:   /* Register physical models to be available on the command line */
602:   PetscFunctionListAdd(&physics,"advect"          ,PhysicsCreate_Advect);

604:   ctx.comm = comm;
605:   ctx.cfl  = 0.9;
606:   ctx.bctype = FVBC_PERIODIC;
607:   ctx.xmin = -1.0;
608:   ctx.xmax = 1.0;
609:   PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
610:   PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
611:   PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
612:   PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);
613:   PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
614:   PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
615:   PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
616:   PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);
617:   PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
618:   PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
619:   PetscOptionsBool("-recursive_split","Split the domain recursively","",ctx.recursive,&ctx.recursive,NULL);
620:   PetscOptionsEnd();

622:   /* Choose the limiter from the list of registered limiters */
623:   PetscFunctionListFind(limiters,lname,&ctx.limit);

626:   /* Choose the physics from the list of registered models */
627:   {
628:     PetscErrorCode (*r)(FVCtx*);
629:     PetscFunctionListFind(physics,physname,&r);
631:     /* Create the physics, will set the number of fields and their names */
632:     (*r)(&ctx);
633:   }

635:   /* Create a DMDA to manage the parallel grid */
636:   DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);
637:   DMSetFromOptions(da);
638:   DMSetUp(da);
639:   /* Inform the DMDA of the field names provided by the physics. */
640:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
641:   for (i=0; i<ctx.physics.dof; i++) {
642:     DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
643:   }
644:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
645:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

647:   /* Set coordinates of cell centers */
648:   DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);

650:   /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
651:   PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope);
652:   PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds);

654:   /* Create a vector to store the solution and to save the initial state */
655:   DMCreateGlobalVector(da,&X);
656:   VecDuplicate(X,&X0);
657:   VecDuplicate(X,&R);

659:   /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/
660:   PetscMalloc1(xm*dof,&index_slow);
661:   PetscMalloc1(xm*dof,&index_fast);
662:   for (i=xs; i<xs+xm; i++) {
663:     if (ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx < 0)
664:       for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
665:     else
666:       for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
667:   }  /* this step need to be changed based on discontinuous point of a */
668:   ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
669:   ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);

671:   /* Create a time-stepping object */
672:   TSCreate(comm,&ts);
673:   TSSetDM(ts,da);
674:   TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
675:   TSRHSSplitSetIS(ts,"slow",ctx.iss);
676:   TSRHSSplitSetIS(ts,"fast",ctx.isf);
677:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);
678:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);

680:   if (ctx.recursive) {
681:     TS subts;
682:     islow = 0;
683:     ifast = 0;
684:     for (i=xs; i<xs+xm; i++) {
685:       PetscReal coord = ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx;
686:       if (coord >= 0 && coord < ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.)
687:         for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
688:       if (coord >= ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.)
689:         for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
690:     }  /* this step need to be changed based on discontinuous point of a */
691:     ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss2);
692:     ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf2);

694:     TSRHSSplitGetSubTS(ts,"fast",&subts);
695:     TSRHSSplitSetIS(subts,"slow",ctx.iss2);
696:     TSRHSSplitSetIS(subts,"fast",ctx.isf2);
697:     TSRHSSplitSetRHSFunction(subts,"slow",NULL,FVRHSFunctionslow2,&ctx);
698:     TSRHSSplitSetRHSFunction(subts,"fast",NULL,FVRHSFunctionfast2,&ctx);
699:   }

701:   /*TSSetType(ts,TSSSP);*/
702:   TSSetType(ts,TSMPRK);
703:   TSSetMaxTime(ts,10);
704:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

706:   /* Compute initial conditions and starting time step */
707:   FVSample(&ctx,da,0,X0);
708:   FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
709:   VecCopy(X0,X);                        /* The function value was not used so we set X=X0 again */
710:   TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
711:   TSSetFromOptions(ts); /* Take runtime options */
712:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
713:   {
714:     PetscInt    steps;
715:     PetscScalar mass_initial,mass_final,mass_difference;

717:     TSSolve(ts,X);
718:     TSGetSolveTime(ts,&ptime);
719:     TSGetStepNumber(ts,&steps);
720:     PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
721:     /* calculate the total mass at initial time and final time */
722:     mass_initial = 0.0;
723:     mass_final   = 0.0;
724:     VecSum(X0,&mass_initial);
725:     VecSum(X,&mass_final);
726:     mass_difference = (ctx.xmax-ctx.xmin)/(PetscScalar)Mx*(mass_final - mass_initial);
727:     PetscPrintf(comm,"Mass difference %g\n",(double)mass_difference);
728:     if (ctx.simulation) {
729:       PetscViewer  fd;
730:       char         filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
731:       Vec          XR;
732:       PetscReal    nrm1,nrmsup;
733:       PetscBool    flg;

735:       PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
737:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
738:       VecDuplicate(X0,&XR);
739:       VecLoad(XR,fd);
740:       PetscViewerDestroy(&fd);
741:       VecAYPX(XR,-1,X);
742:       VecNorm(XR,NORM_1,&nrm1);
743:       VecNorm(XR,NORM_INFINITY,&nrmsup);
744:       nrm1 /= Mx;
745:       PetscPrintf(comm,"Error ||x-x_e||_1 %g  ||x-x_e||_sup %g\n",(double)nrm1,(double)nrmsup);
746:       VecDestroy(&XR);
747:     }
748:   }

750:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
751:   if (draw & 0x1) VecView(X0,PETSC_VIEWER_DRAW_WORLD);
752:   if (draw & 0x2) VecView(X,PETSC_VIEWER_DRAW_WORLD);
753:   if (draw & 0x4) {
754:     Vec Y;
755:     VecDuplicate(X,&Y);
756:     FVSample(&ctx,da,ptime,Y);
757:     VecAYPX(Y,-1,X);
758:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
759:     VecDestroy(&Y);
760:   }

762:   if (view_final) {
763:     PetscViewer viewer;
764:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
765:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
766:     VecView(X,viewer);
767:     PetscViewerPopFormat(viewer);
768:     PetscViewerDestroy(&viewer);
769:   }

771:   /* Clean up */
772:   (*ctx.physics.destroy)(ctx.physics.user);
773:   for (i=0; i<ctx.physics.dof; i++) PetscFree(ctx.physics.fieldname[i]);
774:   PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
775:   PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
776:   ISDestroy(&ctx.iss);
777:   ISDestroy(&ctx.isf);
778:   ISDestroy(&ctx.iss2);
779:   ISDestroy(&ctx.isf2);
780:   VecDestroy(&X);
781:   VecDestroy(&X0);
782:   VecDestroy(&R);
783:   DMDestroy(&da);
784:   TSDestroy(&ts);
785:   PetscFree(index_slow);
786:   PetscFree(index_fast);
787:   PetscFunctionListDestroy(&limiters);
788:   PetscFunctionListDestroy(&physics);
789:   PetscFinalize();
790:   return 0;
791: }

793: /*TEST
794:     build:
795:       requires: !complex
796:       depends: finitevolume1d.c

798:     test:
799:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0

801:     test:
802:       suffix: 2
803:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1
804:       output_file: output/ex5_1.out

806:     test:
807:       suffix: 3
808:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0

810:     test:
811:       suffix: 4
812:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
813:       output_file: output/ex5_3.out

815:     test:
816:       suffix: 5
817:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0 -recursive_split

819:     test:
820:       suffix: 6
821:       args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1 -recursive_split
822: TEST*/