Actual source code: ex5.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  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: PETSC_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;
144:   PetscInt       i,j,k,Mx,dof,xs,xm,len_slow;
145:   PetscReal      hx,cfl_idt = 0;
146:   PetscScalar    *x,*f,*slope;
147:   Vec            Xloc;
148:   DM             da;

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

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

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

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

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

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

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

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

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

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

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

457:   DMDAVecRestoreArray(da,Xloc,&x);
458:   VecRestoreArray(F,&f);
459:   DMDARestoreArray(da,PETSC_TRUE,&slope);
460:   DMRestoreLocalVector(da,&Xloc);
461:   return(0);
462: }

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

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

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

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

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

581:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
582:   comm = PETSC_COMM_WORLD;
583:   PetscMemzero(&ctx,sizeof(ctx));

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

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

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

626:   /* Choose the limiter from the list of registered limiters */
627:   PetscFunctionListFind(limiters,lname,&ctx.limit);
628:   if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);

630:   /* Choose the physics from the list of registered models */
631:   {
632:     PetscErrorCode (*r)(FVCtx*);
633:     PetscFunctionListFind(physics,physname,&r);
634:     if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
635:     /* Create the physics, will set the number of fields and their names */
636:     (*r)(&ctx);
637:   }

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

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

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

658:   /* Create a vector to store the solution and to save the initial state */
659:   DMCreateGlobalVector(da,&X);
660:   VecDuplicate(X,&X0);
661:   VecDuplicate(X,&R);

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

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

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

698:     TSRHSSplitGetSubTS(ts,"fast",&subts);
699:     TSRHSSplitSetIS(subts,"slow",ctx.iss2);
700:     TSRHSSplitSetIS(subts,"fast",ctx.isf2);
701:     TSRHSSplitSetRHSFunction(subts,"slow",NULL,FVRHSFunctionslow2,&ctx);
702:     TSRHSSplitSetRHSFunction(subts,"fast",NULL,FVRHSFunctionfast2,&ctx);
703:   }

705:   /*TSSetType(ts,TSSSP);*/
706:   TSSetType(ts,TSMPRK);
707:   TSSetMaxTime(ts,10);
708:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

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

721:     TSSolve(ts,X);
722:     TSGetSolveTime(ts,&ptime);
723:     TSGetStepNumber(ts,&steps);
724:     PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
725:     /* calculate the total mass at initial time and final time */
726:     mass_initial = 0.0;
727:     mass_final   = 0.0;
728:     VecSum(X0,&mass_initial);
729:     VecSum(X,&mass_final);
730:     mass_difference = (ctx.xmax-ctx.xmin)/(PetscScalar)Mx*(mass_final - mass_initial);
731:     PetscPrintf(comm,"Mass difference %g\n",(double)mass_difference);
732:     if (ctx.simulation) {
733:       PetscViewer  fd;
734:       char         filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
735:       Vec          XR;
736:       PetscReal    nrm1,nrmsup;
737:       PetscBool    flg;
738:       PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
739:       if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
740:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
741:       VecDuplicate(X0,&XR);
742:       VecLoad(XR,fd);
743:       PetscViewerDestroy(&fd);
744:       VecAYPX(XR,-1,X);
745:       VecNorm(XR,NORM_1,&nrm1);
746:       VecNorm(XR,NORM_INFINITY,&nrmsup);
747:       nrm1 /= Mx;
748:       PetscPrintf(comm,"Error ||x-x_e||_1 %g  ||x-x_e||_sup %g\n",(double)nrm1,(double)nrmsup);
749:       VecDestroy(&XR);
750:     }
751:   }

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

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

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

796: /*TEST
797:     build:
798:       requires: !complex c99
799:       depends: finitevolume1d.c

801:     test:
802:       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

804:     test:
805:       suffix: 2
806:       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
807:       output_file: output/ex5_1.out

809:     test:
810:       suffix: 3
811:       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

813:     test:
814:       suffix: 4
815:       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
816:       output_file: output/ex5_3.out

818:     test:
819:       suffix: 5
820:       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

822:     test:
823:       suffix: 6
824:       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
825: TEST*/