Actual source code: ex6.c

  1: /*
  2:   Note:
  3:     -hratio is the ratio between mesh size of coarse 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: 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,PETSC_ERR_ARG_UNKNOWN_TYPE,"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,PETSC_ERR_ARG_UNKNOWN_TYPE,"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:   PetscScalar     *u,*uj,xj,xi;
108:   PetscInt        i,j,k,dof,xs,xm,Mx;
109:   const PetscInt  N = 200;
110:   PetscReal       hs,hf;

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

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

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

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

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

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

200:   DMDAVecGetArray(da,Xloc,&x);
201:   DMDAVecGetArray(da,F,&f);
202:   DMDAGetArray(da,PETSC_TRUE,&slope);                  /* contains ghost points                                           */

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

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

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

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

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

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

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

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

469:   TSGetDM(ts,&da);
470:   DMGetLocalVector(da,&Xloc);
471:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
472:   hxs  = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf;
473:   hxf  = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf);
474:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
475:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
476:   VecZeroEntries(F);
477:   DMDAVecGetArray(da,Xloc,&x);
478:   VecGetArray(F,&f);
479:   DMDAGetArray(da,PETSC_TRUE,&slope);
480:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

482:   if (ctx->bctype == FVBC_OUTFLOW) {
483:     for (i=xs-2; i<0; i++) {
484:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
485:     }
486:     for (i=Mx; i<xs+xm+2; i++) {
487:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
488:     }
489:   }
490:   for (i=xs-1; i<xs+xm+1; i++) {
491:     struct _LimitInfo info;
492:     PetscScalar       *cjmpL,*cjmpR;
493:     if ((i > sf-lsbwidth-2 && i < sf+1) || (i > fs-2 && i < fs+rsbwidth+1)) {
494:       /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
495:       (*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
496:       /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
497:       PetscArrayzero(ctx->cjmpLR,2*dof);
498:       cjmpL = &ctx->cjmpLR[0];
499:       cjmpR = &ctx->cjmpLR[dof];
500:       for (j=0; j<dof; j++) {
501:         PetscScalar jmpL,jmpR;
502:         jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
503:         jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
504:         for (k=0; k<dof; k++) {
505:           cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
506:           cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
507:         }
508:       }
509:       /* Apply limiter to the left and right characteristic jumps */
510:       info.m  = dof;
511:       info.hxs = hxs;
512:       info.hxf = hxf;
513:       (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope);
514:       for (j=0; j<dof; j++) {
515:         PetscScalar tmp = 0;
516:         for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
517:           slope[i*dof+j] = tmp;
518:       }
519:     }
520:   }

522:   for (i=xs; i<xs+xm+1; i++) {
523:     PetscReal   maxspeed;
524:     PetscScalar *uL,*uR;
525:     uL = &ctx->uLR[0];
526:     uR = &ctx->uLR[dof];
527:     if (i == sf-lsbwidth) {
528:       for (j=0; j<dof; j++) {
529:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
530:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
531:       }
532:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
533:       if (i < xs+xm) {
534:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
535:         islow++;
536:       }
537:     }
538:     if (i > sf-lsbwidth && i < sf) {
539:       for (j=0; j<dof; j++) {
540:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
541:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
542:       }
543:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
544:       if (i > xs) {
545:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
546:       }
547:       if (i < xs+xm) {
548:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
549:         islow++;
550:       }
551:     }
552:     if (i == sf) { /* interface between the slow region and the fast region */
553:       for (j=0; j<dof; j++) {
554:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
555:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2;
556:       }
557:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
558:       if (i > xs) {
559:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
560:       }
561:     }
562:     if (i == fs) { /* interface between the fast region and the slow region */
563:       for (j=0; j<dof; j++) {
564:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2;
565:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
566:       }
567:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
568:       if (i < xs+xm) {
569:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
570:         islow++;
571:       }
572:     }
573:     if (i > fs && i < fs+rsbwidth) {
574:       for (j=0; j<dof; j++) {
575:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
576:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
577:       }
578:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
579:       if (i > xs) {
580:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
581:       }
582:       if (i < xs+xm) {
583:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs;
584:         islow++;
585:       }
586:     }
587:     if (i == fs+rsbwidth) {
588:       for (j=0; j<dof; j++) {
589:         uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2;
590:         uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2;
591:       }
592:       (*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed);
593:       if (i > xs) {
594:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs;
595:       }
596:     }
597:   }
598:   DMDAVecRestoreArray(da,Xloc,&x);
599:   VecRestoreArray(F,&f);
600:   DMDARestoreArray(da,PETSC_TRUE,&slope);
601:   DMRestoreLocalVector(da,&Xloc);
602:   return 0;
603: }

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

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

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

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

715: int main(int argc,char *argv[])
716: {
717:   char              lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
718:   PetscFunctionList limiters   = 0,physics = 0;
719:   MPI_Comm          comm;
720:   TS                ts;
721:   DM                da;
722:   Vec               X,X0,R;
723:   FVCtx             ctx;
724:   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;
725:   PetscBool         view_final = PETSC_FALSE;
726:   PetscReal         ptime;
727:   PetscErrorCode    ierr;

729:   PetscInitialize(&argc,&argv,0,help);
730:   comm = PETSC_COMM_WORLD;
731:   PetscMemzero(&ctx,sizeof(ctx));

733:   /* Register limiters to be available on the command line */
734:   PetscFunctionListAdd(&limiters,"upwind"              ,Limit2_Upwind);
735:   PetscFunctionListAdd(&limiters,"lax-wendroff"        ,Limit2_LaxWendroff);
736:   PetscFunctionListAdd(&limiters,"beam-warming"        ,Limit2_BeamWarming);
737:   PetscFunctionListAdd(&limiters,"fromm"               ,Limit2_Fromm);
738:   PetscFunctionListAdd(&limiters,"minmod"              ,Limit2_Minmod);
739:   PetscFunctionListAdd(&limiters,"superbee"            ,Limit2_Superbee);
740:   PetscFunctionListAdd(&limiters,"mc"                  ,Limit2_MC);
741:   PetscFunctionListAdd(&limiters,"koren3"              ,Limit2_Koren3);

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

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

765:   /* Choose the limiter from the list of registered limiters */
766:   PetscFunctionListFind(limiters,lname,&ctx.limit2);

769:   /* Choose the physics from the list of registered models */
770:   {
771:     PetscErrorCode (*r)(FVCtx*);
772:     PetscFunctionListFind(physics,physname,&r);
774:     /* Create the physics, will set the number of fields and their names */
775:     (*r)(&ctx);
776:   }

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

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

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

797:   /* Create a vector to store the solution and to save the initial state */
798:   DMCreateGlobalVector(da,&X);
799:   VecDuplicate(X,&X0);
800:   VecDuplicate(X,&R);

802:   /* create index for slow parts and fast parts,
803:      count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
804:   count_slow = Mx/(1.0+ctx.hratio/3.0);
806:   count_fast = Mx-count_slow;
807:   ctx.sf = count_slow/2;
808:   ctx.fs = ctx.sf+count_fast;
809:   PetscMalloc1(xm*dof,&index_slow);
810:   PetscMalloc1(xm*dof,&index_fast);
811:   PetscMalloc1(6*dof,&index_slowbuffer);
812:   if (((AdvectCtx*)ctx.physics2.user)->a > 0) {
813:     ctx.lsbwidth = 2;
814:     ctx.rsbwidth = 4;
815:   } else {
816:     ctx.lsbwidth = 4;
817:     ctx.rsbwidth = 2;
818:   }
819:   for (i=xs; i<xs+xm; i++) {
820:     if (i < ctx.sf-ctx.lsbwidth || i > ctx.fs+ctx.rsbwidth-1)
821:       for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
822:     else if ((i >= ctx.sf-ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs-1 && i <= ctx.fs+ctx.rsbwidth-1))
823:       for (k=0; k<dof; k++) index_slowbuffer[islowbuffer++] = i*dof+k;
824:     else
825:       for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
826:   }
827:   ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
828:   ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);
829:   ISCreateGeneral(PETSC_COMM_WORLD,islowbuffer,index_slowbuffer,PETSC_COPY_VALUES,&ctx.issb);

831:   /* Create a time-stepping object */
832:   TSCreate(comm,&ts);
833:   TSSetDM(ts,da);
834:   TSSetRHSFunction(ts,R,FVRHSFunction_2WaySplit,&ctx);
835:   TSRHSSplitSetIS(ts,"slow",ctx.iss);
836:   TSRHSSplitSetIS(ts,"slowbuffer",ctx.issb);
837:   TSRHSSplitSetIS(ts,"fast",ctx.isf);
838:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow_2WaySplit,&ctx);
839:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast_2WaySplit,&ctx);
840:   TSRHSSplitSetRHSFunction(ts,"slowbuffer",NULL,FVRHSFunctionslowbuffer_2WaySplit,&ctx);

842:   TSSetType(ts,TSSSP);
843:   /*TSSetType(ts,TSMPRK);*/
844:   TSSetMaxTime(ts,10);
845:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

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

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

922:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
923:   if (draw & 0x1) VecView(X0,PETSC_VIEWER_DRAW_WORLD);
924:   if (draw & 0x2) VecView(X,PETSC_VIEWER_DRAW_WORLD);
925:   if (draw & 0x4) {
926:     Vec Y;
927:     VecDuplicate(X,&Y);
928:     FVSample_2WaySplit(&ctx,da,ptime,Y);
929:     VecAYPX(Y,-1,X);
930:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
931:     VecDestroy(&Y);
932:   }

934:   if (view_final) {
935:     PetscViewer viewer;
936:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
937:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
938:     VecView(X,viewer);
939:     PetscViewerPopFormat(viewer);
940:     PetscViewerDestroy(&viewer);
941:   }

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

965: /*TEST

967:     build:
968:       requires: !complex
969:       depends: finitevolume1d.c

971:     test:
972:       suffix: 1
973:       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

975:     test:
976:       suffix: 2
977:       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
978:       output_file: output/ex6_1.out

980: TEST*/