Actual source code: ex6.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  1: /*
  2:   Note:
  3:     -hratio is the ratio between mesh size of corse grids and fine grids
  4:     -ts_rk_dtratio is the ratio between the large stepsize and the small stepsize
  5: */

  7: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
  8:   "  advection   - Constant coefficient scalar advection\n"
  9:   "                u_t       + (a*u)_x               = 0\n"
 10:   "  for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
 11:   "                hxs  = hratio*hxf \n"
 12:   "  where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n"
 13:   "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
 14:   "                the states across shocks and rarefactions\n"
 15:   "  simulation  - use reference solution which is generated by smaller time step size to be true solution,\n"
 16:   "                also the reference solution should be generated by user and stored in a binary file.\n"
 17:   "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
 18:   "Several initial conditions can be chosen with -initial N\n\n"
 19:   "The problem size should be set with -da_grid_x M\n\n";

 21:  #include <petscts.h>
 22:  #include <petscdm.h>
 23:  #include <petscdmda.h>
 24:  #include <petscdraw.h>
 25: #include "finitevolume1d.h"

 27: PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }

 29: /* --------------------------------- Advection ----------------------------------- */
 30: typedef struct {
 31:   PetscReal a;                  /* advective velocity */
 32: } AdvectCtx;

 34: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
 35: {
 36:   AdvectCtx *ctx = (AdvectCtx*)vctx;
 37:   PetscReal speed;

 40:   speed     = ctx->a;
 41:   flux[0]   = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
 42:   *maxspeed = speed;
 43:   return(0);
 44: }

 46: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
 47: {
 48:   AdvectCtx *ctx = (AdvectCtx*)vctx;

 51:   X[0]      = 1.;
 52:   Xi[0]     = 1.;
 53:   speeds[0] = ctx->a;
 54:   return(0);
 55: }

 57: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
 58: {
 59:   AdvectCtx *ctx = (AdvectCtx*)vctx;
 60:   PetscReal a    = ctx->a,x0;

 63:   switch (bctype) {
 64:     case FVBC_OUTFLOW:   x0 = x-a*t; break;
 65:     case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
 66:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
 67:   }
 68:   switch (initial) {
 69:     case 0: u[0] = (x0 < 0) ? 1 : -1; break;
 70:     case 1: u[0] = (x0 < 0) ? -1 : 1; break;
 71:     case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
 72:     case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
 73:     case 4: u[0] = PetscAbs(x0); break;
 74:     case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
 75:     case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
 76:     case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
 77:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
 78:   }
 79:   return(0);
 80: }

 82: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
 83: {
 85:   AdvectCtx      *user;

 88:   PetscNew(&user);
 89:   ctx->physics2.sample2         = PhysicsSample_Advect;
 90:   ctx->physics2.riemann2        = PhysicsRiemann_Advect;
 91:   ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
 92:   ctx->physics2.destroy         = PhysicsDestroy_SimpleFree;
 93:   ctx->physics2.user            = user;
 94:   ctx->physics2.dof             = 1;
 95:   PetscStrallocpy("u",&ctx->physics2.fieldname[0]);
 96:   user->a = 1;
 97:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
 98:   {
 99:     PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
100:   }
101:   PetscOptionsEnd();
102:   return(0);
103: }

105: PetscErrorCode FVSample_2WaySplit(FVCtx *ctx,DM da,PetscReal time,Vec U)
106: {
107:   PetscErrorCode  ierr;
108:   PetscScalar     *u,*uj,xj,xi;
109:   PetscInt        i,j,k,dof,xs,xm,Mx;
110:   const PetscInt  N = 200;
111:   PetscReal       hs,hf;

114:   if (!ctx->physics2.sample2) SETERRQ(PETSC_COMM_SELF,1,"Physics has not provided a sampling function");
115:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
116:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
117:   DMDAVecGetArray(da,U,&u);
118:   PetscMalloc1(dof,&uj);
119:   hs   = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
120:   hf   = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
121:   for (i=xs; i<xs+xm; i++) {
122:     if (i < ctx->sf) {
123:       xi = ctx->xmin+0.5*hs+i*hs;
124:       /* Integrate over cell i using trapezoid rule with N points. */
125:       for (k=0; k<dof; k++) u[i*dof+k] = 0;
126:       for (j=0; j<N+1; j++) {
127:         xj = xi+hs*(j-N/2)/(PetscReal)N;
128:         (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
129:         for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
130:       }
131:     } else if (i < ctx->fs) {
132:       xi = ctx->xmin+ctx->sf*hs+0.5*hf+(i-ctx->sf)*hf;
133:       /* Integrate over cell i using trapezoid rule with N points. */
134:       for (k=0; k<dof; k++) u[i*dof+k] = 0;
135:       for (j=0; j<N+1; j++) {
136:         xj = xi+hf*(j-N/2)/(PetscReal)N;
137:         (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
138:         for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
139:       }
140:     } else {
141:       xi = ctx->xmin+ctx->sf*hs+(ctx->fs-ctx->sf)*hf+0.5*hs+(i-ctx->fs)*hs;
142:       /* Integrate over cell i using trapezoid rule with N points. */
143:       for (k=0; k<dof; k++) u[i*dof+k] = 0;
144:       for (j=0; j<N+1; j++) {
145:         xj = xi+hs*(j-N/2)/(PetscReal)N;
146:         (*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
147:         for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
148:       }
149:     }
150:   }
151:   DMDAVecRestoreArray(da,U,&u);
152:   PetscFree(uj);
153:   return(0);
154: }

156: static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1)
157: {
158:   PetscErrorCode    ierr;
159:   Vec               Y;
160:   PetscInt          i,Mx;
161:   const PetscScalar *ptr_X,*ptr_Y;
162:   PetscReal         hs,hf;

165:   VecGetSize(X,&Mx);
166:   VecDuplicate(X,&Y);
167:   FVSample_2WaySplit(ctx,da,t,Y);
168:   hs   = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
169:   hf   = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
170:   VecGetArrayRead(X,&ptr_X);
171:   VecGetArrayRead(Y,&ptr_Y);
172:   for (i=0; i<Mx; i++) {
173:     if (i < ctx->sf || i > ctx->fs -1) *nrm1 +=  hs*PetscAbs(ptr_X[i]-ptr_Y[i]);
174:     else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]);
175:   }
176:   VecRestoreArrayRead(X,&ptr_X);
177:   VecRestoreArrayRead(Y,&ptr_Y);
178:   VecDestroy(&Y);
179:   return(0);
180: }

182: PetscErrorCode FVRHSFunction_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
183: {
184:   FVCtx          *ctx = (FVCtx*)vctx;
186:   PetscInt       i,j,k,Mx,dof,xs,xm,sf = ctx->sf,fs = ctx->fs;
187:   PetscReal      hxf,hxs,cfl_idt = 0;
188:   PetscScalar    *x,*f,*slope;
189:   Vec            Xloc;
190:   DM             da;

193:   TSGetDM(ts,&da);
194:   DMGetLocalVector(da,&Xloc);                          /* Xloc contains ghost points                                     */
195:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);   /* Mx is the number of center points                              */
196:   hxs  = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
197:   hxf  = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
198:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);       /* X is solution vector which does not contain ghost points       */
199:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);

201:   VecZeroEntries(F);                                   /* F is the right hand side function corresponds to center points */

203:   DMDAVecGetArray(da,Xloc,&x);
204:   DMDAVecGetArray(da,F,&f);
205:   DMDAGetArray(da,PETSC_TRUE,&slope);                  /* contains ghost points                                           */

207:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

209:   if (ctx->bctype == FVBC_OUTFLOW) {
210:     for (i=xs-2; i<0; i++) {
211:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
212:     }
213:     for (i=Mx; i<xs+xm+2; i++) {
214:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
215:     }
216:   }
217:   for (i=xs-1; i<xs+xm+1; i++) {
218:     struct _LimitInfo info;
219:     PetscScalar       *cjmpL,*cjmpR;
220:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
221:     (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
222:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
223:     PetscArrayzero(ctx->cjmpLR,2*dof);
224:     cjmpL = &ctx->cjmpLR[0];
225:     cjmpR = &ctx->cjmpLR[dof];
226:     for (j=0; j<dof; j++) {
227:       PetscScalar jmpL,jmpR;
228:       jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
229:       jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
230:       for (k=0; k<dof; k++) {
231:         cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
232:         cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
233:       }
234:     }
235:     /* Apply limiter to the left and right characteristic jumps */
236:     info.m  = dof;
237:     info.hxs = hxs;
238:     info.hxf = hxf;
239:     (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
240:     for (j=0; j<dof; j++) {
241:       PetscScalar tmp = 0;
242:       for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
243:       slope[i*dof+j] = tmp;
244:     }
245:   }

247:   for (i=xs; i<xs+xm+1; i++) {
248:     PetscReal   maxspeed;
249:     PetscScalar *uL,*uR;
250:     uL = &ctx->uLR[0];
251:     uR = &ctx->uLR[dof];
252:     if (i < sf) { /* slow region */
253:       for (j=0; j<dof; j++) {
254:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
255:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
256:       }
257:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
258:       if (i > xs) {
259:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs;
260:       }
261:       if (i < xs+xm) {
262:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs;
263:       }
264:     } else if (i == sf) { /* interface between the slow region and the fast region */
265:       for (j=0; j<dof; j++) {
266:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
267:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
268:       }
269:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
270:       if (i > xs) {
271:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs;
272:       }
273:       if (i < xs+xm) {
274:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf;
275:       }
276:     } else if (i > sf && i < fs) { /* fast region */
277:       for (j=0; j<dof; j++) {
278:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
279:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
280:       }
281:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
282:       if (i > xs) {
283:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf;
284:       }
285:       if (i < xs+xm) {
286:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf;
287:       }
288:     } else if (i == fs) { /* interface between the fast region and the slow region */
289:       for (j=0; j<dof; j++) {
290:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
291:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
292:       }
293:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
294:       if (i > xs) {
295:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf;
296:       }
297:       if (i < xs+xm) {
298:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs;
299:       }
300:     } else { /* slow region */
301:       for (j=0; j<dof; j++) {
302:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
303:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
304:       }
305:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
306:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */
307:       if (i > xs) {
308:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs;
309:       }
310:       if (i < xs+xm) {
311:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs;
312:       }
313:     }
314:   }
315:   DMDAVecRestoreArray(da,Xloc,&x);
316:   DMDAVecRestoreArray(da,F,&f);
317:   DMDARestoreArray(da,PETSC_TRUE,&slope);
318:   DMRestoreLocalVector(da,&Xloc);
319:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));
320:   if (0) {
321:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
322:     PetscReal dt,tnow;
323:     TSGetTimeStep(ts,&dt);
324:     TSGetTime(ts,&tnow);
325:     if (dt > 0.5/ctx->cfl_idt) {
326:       if (1) {
327:         PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
328:       } else SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt));
329:     }
330:   }
331:   return(0);
332: }

334: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
335: PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
336: {
337:   FVCtx          *ctx = (FVCtx*)vctx;
339:   PetscInt       i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth;
340:   PetscReal      hxs,hxf,cfl_idt = 0;
341:   PetscScalar    *x,*f,*slope;
342:   Vec            Xloc;
343:   DM             da;

346:   TSGetDM(ts,&da);
347:   DMGetLocalVector(da,&Xloc);
348:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
349:   hxs  = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
350:   hxf  = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
351:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
352:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
353:   VecZeroEntries(F);
354:   DMDAVecGetArray(da,Xloc,&x);
355:   VecGetArray(F,&f);
356:   DMDAGetArray(da,PETSC_TRUE,&slope);
357:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

359:   if (ctx->bctype == FVBC_OUTFLOW) {
360:     for (i=xs-2; i<0; i++) {
361:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
362:     }
363:     for (i=Mx; i<xs+xm+2; i++) {
364:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
365:     }
366:   }
367:   for (i=xs-1; i<xs+xm+1; i++) {
368:     struct _LimitInfo info;
369:     PetscScalar       *cjmpL,*cjmpR;
370:     if (i < sf-lsbwidth+1 || i > fs+rsbwidth-2) { /* slow components and the first and last fast components */
371:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
372:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
373:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
374:       PetscArrayzero(ctx->cjmpLR,2*dof);
375:       cjmpL = &ctx->cjmpLR[0];
376:       cjmpR = &ctx->cjmpLR[dof];
377:       for (j=0; j<dof; j++) {
378:         PetscScalar jmpL,jmpR;
379:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
380:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
381:         for (k=0; k<dof; k++) {
382:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
383:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
384:         }
385:       }
386:       /* Apply limiter to the left and right characteristic jumps */
387:       info.m  = dof;
388:       info.hxs = hxs;
389:       info.hxf = hxf;
390:       (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
391:       for (j=0; j<dof; j++) {
392:         PetscScalar tmp = 0;
393:         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
394:           slope[i*dof+j] = tmp;
395:       }
396:     }
397:   }

399:   for (i=xs; i<xs+xm+1; i++) {
400:     PetscReal   maxspeed;
401:     PetscScalar *uL,*uR;
402:     uL = &ctx->uLR[0];
403:     uR = &ctx->uLR[dof];
404:     if (i < sf-lsbwidth) { /* slow region */
405:       for (j=0; j<dof; j++) {
406:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
407:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
408:       }
409:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
410:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */
411:       if (i > xs) {
412:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
413:       }
414:       if (i < xs+xm) {
415:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
416:         islow++;
417:       }
418:     }
419:     if (i == sf-lsbwidth) { /* interface between the slow region and the fast region */
420:       for (j=0; j<dof; j++) {
421:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
422:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
423:       }
424:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
425:       if (i > xs) {
426:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
427:       }
428:     }
429:     if (i == fs+rsbwidth) { /* slow region */
430:       for (j=0; j<dof; j++) {
431:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
432:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
433:       }
434:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
435:       if (i < xs+xm) {
436:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
437:         islow++;
438:       }
439:     }
440:     if (i > fs+rsbwidth) { /* slow region */
441:       for (j=0; j<dof; j++) {
442:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
443:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
444:       }
445:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
446:       if (i > xs) {
447:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
448:       }
449:       if (i < xs+xm) {
450:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
451:         islow++;
452:       }
453:     }
454:   }
455:   DMDAVecRestoreArray(da,Xloc,&x);
456:   VecRestoreArray(F,&f);
457:   DMDARestoreArray(da,PETSC_TRUE,&slope);
458:   DMRestoreLocalVector(da,&Xloc);
459:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));
460:   return(0);
461: }

463: PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
464: {
465:   FVCtx          *ctx = (FVCtx*)vctx;
467:   PetscInt       i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth;
468:   PetscReal      hxs,hxf;
469:   PetscScalar    *x,*f,*slope;
470:   Vec            Xloc;
471:   DM             da;

474:   TSGetDM(ts,&da);
475:   DMGetLocalVector(da,&Xloc);
476:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
477:   hxs  = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
478:   hxf  = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
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);

487:   if (ctx->bctype == FVBC_OUTFLOW) {
488:     for (i=xs-2; i<0; i++) {
489:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
490:     }
491:     for (i=Mx; i<xs+xm+2; i++) {
492:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
493:     }
494:   }
495:   for (i=xs-1; i<xs+xm+1; i++) {
496:     struct _LimitInfo info;
497:     PetscScalar       *cjmpL,*cjmpR;
498:     if ((i > sf-lsbwidth-2 && i < sf+1) || (i > fs-2 && i < fs+rsbwidth+1)) {
499:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
500:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
501:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
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.hxs = hxs;
517:       info.hxf = hxf;
518:       (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
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:     uL = &ctx->uLR[0];
531:     uR = &ctx->uLR[dof];
532:     if (i == sf-lsbwidth) {
533:       for (j=0; j<dof; j++) {
534:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
535:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
536:       }
537:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
538:       if (i < xs+xm) {
539:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
540:         islow++;
541:       }
542:     }
543:     if (i > sf-lsbwidth && i < sf) {
544:       for (j=0; j<dof; j++) {
545:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
546:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
547:       }
548:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
549:       if (i > xs) {
550:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
551:       }
552:       if (i < xs+xm) {
553:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
554:         islow++;
555:       }
556:     }
557:     if (i == sf) { /* interface between the slow region and the fast region */
558:       for (j=0; j<dof; j++) {
559:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
560:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
561:       }
562:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
563:       if (i > xs) {
564:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
565:       }
566:     }
567:     if (i == fs) { /* interface between the fast region and the slow region */
568:       for (j=0; j<dof; j++) {
569:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
570:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
571:       }
572:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
573:       if (i < xs+xm) {
574:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
575:         islow++;
576:       }
577:     }
578:     if (i > fs && i < fs+rsbwidth) {
579:       for (j=0; j<dof; j++) {
580:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
581:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
582:       }
583:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
584:       if (i > xs) {
585:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
586:       }
587:       if (i < xs+xm) {
588:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
589:         islow++;
590:       }
591:     }
592:     if (i == fs+rsbwidth) {
593:       for (j=0; j<dof; j++) {
594:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
595:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
596:       }
597:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
598:       if (i > xs) {
599:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
600:       }
601:     }
602:   }
603:   DMDAVecRestoreArray(da,Xloc,&x);
604:   VecRestoreArray(F,&f);
605:   DMDARestoreArray(da,PETSC_TRUE,&slope);
606:   DMRestoreLocalVector(da,&Xloc);
607:   return(0);
608: }

610: /* --------------------------------- Finite Volume Solver for fast  parts ----------------------------------- */
611: PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
612: {
613:   FVCtx          *ctx = (FVCtx*)vctx;
615:   PetscInt       i,j,k,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs;
616:   PetscReal      hxs,hxf;
617:   PetscScalar    *x,*f,*slope;
618:   Vec            Xloc;
619:   DM             da;

622:   TSGetDM(ts,&da);
623:   DMGetLocalVector(da,&Xloc);
624:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
625:   hxs  = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
626:   hxf  = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
627:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
628:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
629:   VecZeroEntries(F);
630:   DMDAVecGetArray(da,Xloc,&x);
631:   VecGetArray(F,&f);
632:   DMDAGetArray(da,PETSC_TRUE,&slope);
633:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

635:   if (ctx->bctype == FVBC_OUTFLOW) {
636:     for (i=xs-2; i<0; i++) {
637:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
638:     }
639:     for (i=Mx; i<xs+xm+2; i++) {
640:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
641:     }
642:   }
643:   for (i=xs-1; i<xs+xm+1; i++) {
644:     struct _LimitInfo info;
645:     PetscScalar       *cjmpL,*cjmpR;
646:     if (i > sf-2 && i < fs+1) {
647:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
648:       PetscArrayzero(ctx->cjmpLR,2*dof);
649:       cjmpL = &ctx->cjmpLR[0];
650:       cjmpR = &ctx->cjmpLR[dof];
651:       for (j=0; j<dof; j++) {
652:         PetscScalar jmpL,jmpR;
653:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
654:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
655:         for (k=0; k<dof; k++) {
656:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
657:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
658:         }
659:       }
660:       /* Apply limiter to the left and right characteristic jumps */
661:       info.m  = dof;
662:       info.hxs = hxs;
663:       info.hxf = hxf;
664:       (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
665:       for (j=0; j<dof; j++) {
666:       PetscScalar tmp = 0;
667:       for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
668:         slope[i*dof+j] = tmp;
669:       }
670:     }
671:   }

673:   for (i=xs; i<xs+xm+1; i++) {
674:     PetscReal   maxspeed;
675:     PetscScalar *uL,*uR;
676:     uL = &ctx->uLR[0];
677:     uR = &ctx->uLR[dof];
678:     if (i == sf) { /* interface between the slow region and the fast region */
679:       for (j=0; j<dof; j++) {
680:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
681:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
682:       }
683:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
684:       if (i < xs+xm) {
685:         for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf;
686:         ifast++;
687:       }
688:     }
689:     if (i > sf && i < fs) { /* fast region */
690:       for (j=0; j<dof; j++) {
691:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
692:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
693:       }
694:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
695:       if (i > xs) {
696:         for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf;
697:       }
698:       if (i < xs+xm) {
699:         for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf;
700:         ifast++;
701:       }
702:     }
703:     if (i == fs) { /* interface between the fast region and the slow region */
704:       for (j=0; j<dof; j++) {
705:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
706:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
707:       }
708:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
709:       if (i > xs) {
710:         for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf;
711:       }
712:     }
713:   }
714:   DMDAVecRestoreArray(da,Xloc,&x);
715:   VecRestoreArray(F,&f);
716:   DMDARestoreArray(da,PETSC_TRUE,&slope);
717:   DMRestoreLocalVector(da,&Xloc);
718:   return(0);
719: }

721: int main(int argc,char *argv[])
722: {
723:   char              lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
724:   PetscFunctionList limiters   = 0,physics = 0;
725:   MPI_Comm          comm;
726:   TS                ts;
727:   DM                da;
728:   Vec               X,X0,R;
729:   FVCtx             ctx;
730:   PetscInt          i,k,dof,xs,xm,Mx,draw = 0,count_slow,count_fast,islow = 0,ifast =0,islowbuffer = 0,*index_slow,*index_fast,*index_slowbuffer;
731:   PetscBool         view_final = PETSC_FALSE;
732:   PetscReal         ptime;
733:   PetscErrorCode    ierr;

735:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
736:   comm = PETSC_COMM_WORLD;
737:   PetscMemzero(&ctx,sizeof(ctx));

739:   /* Register limiters to be available on the command line */
740:   PetscFunctionListAdd(&limiters,"upwind"              ,Limit2_Upwind);
741:   PetscFunctionListAdd(&limiters,"lax-wendroff"        ,Limit2_LaxWendroff);
742:   PetscFunctionListAdd(&limiters,"beam-warming"        ,Limit2_BeamWarming);
743:   PetscFunctionListAdd(&limiters,"fromm"               ,Limit2_Fromm);
744:   PetscFunctionListAdd(&limiters,"minmod"              ,Limit2_Minmod);
745:   PetscFunctionListAdd(&limiters,"superbee"            ,Limit2_Superbee);
746:   PetscFunctionListAdd(&limiters,"mc"                  ,Limit2_MC);
747:   PetscFunctionListAdd(&limiters,"koren3"              ,Limit2_Koren3);

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

752:   ctx.comm = comm;
753:   ctx.cfl  = 0.9;
754:   ctx.bctype = FVBC_PERIODIC;
755:   ctx.xmin = -1.0;
756:   ctx.xmax = 1.0;
757:   PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
758:   PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
759:   PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
760:   PetscOptionsFList("-limit","Name of flux imiter to use","",limiters,lname,lname,sizeof(lname),NULL);
761:   PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
762:   PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
763:   PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
764:   PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
765:   PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);
766:   PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
767:   PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
768:   PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL);
769:   PetscOptionsEnd();

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

775:   /* Choose the physics from the list of registered models */
776:   {
777:     PetscErrorCode (*r)(FVCtx*);
778:     PetscFunctionListFind(physics,physname,&r);
779:     if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
780:     /* Create the physics, will set the number of fields and their names */
781:     (*r)(&ctx);
782:   }

784:   /* Create a DMDA to manage the parallel grid */
785:   DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics2.dof,2,NULL,&da);
786:   DMSetFromOptions(da);
787:   DMSetUp(da);
788:   /* Inform the DMDA of the field names provided by the physics. */
789:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
790:   for (i=0; i<ctx.physics2.dof; i++) {
791:     DMDASetFieldName(da,i,ctx.physics2.fieldname[i]);
792:   }
793:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
794:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

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

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

803:   /* Create a vector to store the solution and to save the initial state */
804:   DMCreateGlobalVector(da,&X);
805:   VecDuplicate(X,&X0);
806:   VecDuplicate(X,&R);

808:   /* create index for slow parts and fast parts,
809:      count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
810:   count_slow = Mx/(1.0+ctx.hratio/3.0);
811:   if (count_slow%2) SETERRQ(PETSC_COMM_WORLD,1,"Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio/3) is even");
812:   count_fast = Mx-count_slow;
813:   ctx.sf = count_slow/2;
814:   ctx.fs = ctx.sf+count_fast;
815:   PetscMalloc1(xm*dof,&index_slow);
816:   PetscMalloc1(xm*dof,&index_fast);
817:   PetscMalloc1(6*dof,&index_slowbuffer);
818:   if (((AdvectCtx*)ctx.physics2.user)->a > 0) {
819:     ctx.lsbwidth = 2;
820:     ctx.rsbwidth = 4;
821:   } else {
822:     ctx.lsbwidth = 4;
823:     ctx.rsbwidth = 2;
824:   }
825:   for (i=xs; i<xs+xm; i++) {
826:     if (i < ctx.sf-ctx.lsbwidth || i > ctx.fs+ctx.rsbwidth-1)
827:       for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
828:     else if ((i >= ctx.sf-ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs-1 && i <= ctx.fs+ctx.rsbwidth-1))
829:       for (k=0; k<dof; k++) index_slowbuffer[islowbuffer++] = i*dof+k;
830:     else
831:       for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
832:   }
833:   ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
834:   ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);
835:   ISCreateGeneral(PETSC_COMM_WORLD,islowbuffer,index_slowbuffer,PETSC_COPY_VALUES,&ctx.issb);

837:   /* Create a time-stepping object */
838:   TSCreate(comm,&ts);
839:   TSSetDM(ts,da);
840:   TSSetRHSFunction(ts,R,FVRHSFunction_2WaySplit,&ctx);
841:   TSRHSSplitSetIS(ts,"slow",ctx.iss);
842:   TSRHSSplitSetIS(ts,"slowbuffer",ctx.issb);
843:   TSRHSSplitSetIS(ts,"fast",ctx.isf);
844:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow_2WaySplit,&ctx);
845:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast_2WaySplit,&ctx);
846:   TSRHSSplitSetRHSFunction(ts,"slowbuffer",NULL,FVRHSFunctionslowbuffer_2WaySplit,&ctx);

848:   TSSetType(ts,TSSSP);
849:   /*TSSetType(ts,TSMPRK);*/
850:   TSSetMaxTime(ts,10);
851:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

853:   /* Compute initial conditions and starting time step */
854:   FVSample_2WaySplit(&ctx,da,0,X0);
855:   FVRHSFunction_2WaySplit(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
856:   VecCopy(X0,X);                        /* The function value was not used so we set X=X0 again */
857:   TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
858:   TSSetFromOptions(ts); /* Take runtime options */
859:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
860:   {
861:     PetscInt          steps;
862:     PetscScalar       mass_initial,mass_final,mass_difference,mass_differenceg;
863:     const PetscScalar *ptr_X,*ptr_X0;
864:     const PetscReal   hs = (ctx.xmax-ctx.xmin)*3.0/4.0/count_slow;
865:     const PetscReal   hf = (ctx.xmax-ctx.xmin)/4.0/count_fast;

867:     TSSolve(ts,X);
868:     TSGetSolveTime(ts,&ptime);
869:     TSGetStepNumber(ts,&steps);
870:     /* calculate the total mass at initial time and final time */
871:     mass_initial = 0.0;
872:     mass_final   = 0.0;
873:     DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0);
874:     DMDAVecGetArrayRead(da,X,(void*)&ptr_X);
875:     for (i=xs;i<xs+xm;i++) {
876:       if (i < ctx.sf || i > ctx.fs-1) {
877:         for (k=0; k<dof; k++) {
878:           mass_initial = mass_initial + hs*ptr_X0[i*dof+k];
879:           mass_final = mass_final + hs*ptr_X[i*dof+k];
880:         }
881:       } else {
882:         for (k=0; k<dof; k++) {
883:           mass_initial = mass_initial + hf*ptr_X0[i*dof+k];
884:           mass_final = mass_final + hf*ptr_X[i*dof+k];
885:         }
886:       }
887:     }
888:     DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0);
889:     DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X);
890:     mass_difference = mass_final - mass_initial;
891:     MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm);
892:     PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg);
893:     PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
894:     PetscPrintf(comm,"Maximum allowable stepsize according to CFL %g\n",(double)(1.0/ctx.cfl_idt));
895:     if (ctx.exact) {
896:       PetscReal nrm1=0;
897:       SolutionErrorNorms_2WaySplit(&ctx,da,ptime,X,&nrm1);
898:       PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
899:     }
900:     if (ctx.simulation) {
901:       PetscReal    nrm1=0;
902:       PetscViewer  fd;
903:       char         filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
904:       Vec          XR;
905:       PetscBool    flg;
906:       const PetscScalar  *ptr_XR;
907:       PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
908:       if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
909:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
910:       VecDuplicate(X0,&XR);
911:       VecLoad(XR,fd);
912:       PetscViewerDestroy(&fd);
913:       VecGetArrayRead(X,&ptr_X);
914:       VecGetArrayRead(XR,&ptr_XR);
915:       for(i=xs;i<xs+xm;i++) {
916:         if(i < ctx.sf || i > ctx.fs-1)
917:           for (k=0; k<dof; k++) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]);
918:         else
919:           for (k=0; k<dof; k++) nrm1 = nrm1 + hf*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]);
920:       }
921:       VecRestoreArrayRead(X,&ptr_X);
922:       VecRestoreArrayRead(XR,&ptr_XR);
923:       PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
924:       VecDestroy(&XR);
925:     }
926:   }

928:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
929:   if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
930:   if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
931:   if (draw & 0x4) {
932:     Vec Y;
933:     VecDuplicate(X,&Y);
934:     FVSample_2WaySplit(&ctx,da,ptime,Y);
935:     VecAYPX(Y,-1,X);
936:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
937:     VecDestroy(&Y);
938:   }

940:   if (view_final) {
941:     PetscViewer viewer;
942:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
943:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
944:     VecView(X,viewer);
945:     PetscViewerPopFormat(viewer);
946:     PetscViewerDestroy(&viewer);
947:   }

949:   /* Clean up */
950:   (*ctx.physics2.destroy)(ctx.physics2.user);
951:   for (i=0; i<ctx.physics2.dof; i++) {PetscFree(ctx.physics2.fieldname[i]);}
952:   PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
953:   PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
954:   VecDestroy(&X);
955:   VecDestroy(&X0);
956:   VecDestroy(&R);
957:   DMDestroy(&da);
958:   TSDestroy(&ts);
959:   ISDestroy(&ctx.iss);
960:   ISDestroy(&ctx.isf);
961:   ISDestroy(&ctx.issb);
962:   PetscFree(index_slow);
963:   PetscFree(index_fast);
964:   PetscFree(index_slowbuffer);
965:   PetscFunctionListDestroy(&limiters);
966:   PetscFunctionListDestroy(&physics);
967:   PetscFinalize();
968:   return ierr;
969: }

971: /*TEST

973:     build:
974:       requires: !complex c99
975:       depends: finitevolume1d.c

977:     test:
978:       suffix: 1
979:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22

981:     test:
982:       suffix: 2
983:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0
984:       output_file: output/ex6_1.out

986: TEST*/