Actual source code: ex7.c

  1: static const char help[] = "1D periodic Finite Volume solver by a particular slope limiter with semidiscrete time stepping.\n"
  2:   "  advection   - Constant coefficient scalar advection\n"
  3:   "                u_t       + (a*u)_x               = 0\n"
  4:   "  for this toy problem, we choose different meshsizes for different sub-domains, say\n"
  5:   "                hxs  = (xmax - xmin)/2.0*(hratio+1.0)/Mx, \n"
  6:   "                hxf  = (xmax - xmin)/2.0*(1.0+1.0/hratio)/Mx, \n"
  7:   "  with x belongs to (xmin,xmax), the number of total mesh points is Mx and the ratio between the meshsize of corse\n\n"
  8:   "  grids and fine grids is hratio.\n"
  9:   "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
 10:   "                the states across shocks and rarefactions\n"
 11:   "  simulation  - use reference solution which is generated by smaller time step size to be true solution,\n"
 12:   "                also the reference solution should be generated by user and stored in a binary file.\n"
 13:   "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
 14:   "Several initial conditions can be chosen with -initial N\n\n"
 15:   "The problem size should be set with -da_grid_x M\n\n"
 16:   "This script choose the slope limiter by biased second-order upwind procedure which is proposed by Van Leer in 1994\n"
 17:   "                             u(x_(k+1/2),t) = u(x_k,t) + phi(x_(k+1/2),t)*(u(x_k,t)-u(x_(k-1),t))                 \n"
 18:   "                     limiter phi(x_(k+1/2),t) = max(0,min(r(k+1/2),min(2,gamma(k+1/2)*r(k+1/2)+alpha(k+1/2))))    \n"
 19:   "                             r(k+1/2) = (u(x_(k+1))-u(x_k))/(u(x_k)-u(x_(k-1)))                                   \n"
 20:   "                             alpha(k+1/2) = (h_k*h_(k+1))/(h_(k-1)+h_k)/(h_(k-1)+h_k+h_(k+1))                     \n"
 21:   "                             gamma(k+1/2) = h_k*(h_(k-1)+h_k)/(h_k+h_(k+1))/(h_(k-1)+h_k+h_(k+1))                 \n";

 23: #include <petscts.h>
 24: #include <petscdm.h>
 25: #include <petscdmda.h>
 26: #include <petscdraw.h>
 27: #include <petscmath.h>

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

 31: /* --------------------------------- Finite Volume data structures ----------------------------------- */

 33: typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType;
 34: static const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0};

 36: typedef struct {
 37:   PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*);
 38:   PetscErrorCode (*flux)(void*,const PetscScalar*,PetscScalar*,PetscReal*);
 39:   PetscErrorCode (*destroy)(void*);
 40:   void           *user;
 41:   PetscInt       dof;
 42:   char           *fieldname[16];
 43: } PhysicsCtx;

 45: typedef struct {
 46:   PhysicsCtx  physics;
 47:   MPI_Comm    comm;
 48:   char        prefix[256];

 50:   /* Local work arrays */
 51:   PetscScalar *flux;            /* Flux across interface                                                      */
 52:   PetscReal   *speeds;          /* Speeds of each wave                                                        */
 53:   PetscScalar *u;               /* value at face                                                              */

 55:   PetscReal   cfl_idt;          /* Max allowable value of 1/Delta t                                           */
 56:   PetscReal   cfl;
 57:   PetscReal   xmin,xmax;
 58:   PetscInt    initial;
 59:   PetscBool   exact;
 60:   PetscBool   simulation;
 61:   FVBCType    bctype;
 62:   PetscInt    hratio;           /* hratio = hslow/hfast */
 63:   IS          isf,iss;
 64:   PetscInt    sf,fs;            /* slow-fast and fast-slow interfaces */
 65: } FVCtx;

 67: /* --------------------------------- Physics ----------------------------------- */
 68: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
 69: {

 73:   PetscFree(vctx);
 74:   return(0);
 75: }

 77: /* --------------------------------- Advection ----------------------------------- */
 78: typedef struct {
 79:   PetscReal a;                  /* advective velocity */
 80: } AdvectCtx;

 82: static PetscErrorCode PhysicsFlux_Advect(void *vctx,const PetscScalar *u,PetscScalar *flux,PetscReal *maxspeed)
 83: {
 84:   AdvectCtx *ctx = (AdvectCtx*)vctx;
 85:   PetscReal speed;

 88:   speed     = ctx->a;
 89:   flux[0]   = speed*u[0];
 90:   *maxspeed = speed;
 91:   return(0);
 92: }

 94: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
 95: {
 96:   AdvectCtx *ctx = (AdvectCtx*)vctx;
 97:   PetscReal a    = ctx->a,x0;

100:   switch (bctype) {
101:     case FVBC_OUTFLOW:   x0 = x-a*t; break;
102:     case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
103:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
104:   }
105:   switch (initial) {
106:     case 0: u[0] = (x0 < 0) ? 1 : -1; break;
107:     case 1: u[0] = (x0 < 0) ? -1 : 1; break;
108:     case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
109:     case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
110:     case 4: u[0] = PetscAbs(x0); break;
111:     case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
112:     case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
113:     case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
114:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
115:   }
116:   return(0);
117: }

119: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
120: {
122:   AdvectCtx      *user;

125:   PetscNew(&user);
126:   ctx->physics.sample         = PhysicsSample_Advect;
127:   ctx->physics.flux           = PhysicsFlux_Advect;
128:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
129:   ctx->physics.user           = user;
130:   ctx->physics.dof            = 1;
131:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
132:   user->a = 1;
133:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
134:   {
135:     PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
136:   }
137:   PetscOptionsEnd();
138:   return(0);
139: }

141: /* --------------------------------- Finite Volume Solver ----------------------------------- */

143: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
144: {
145:   FVCtx          *ctx = (FVCtx*)vctx;
147:   PetscInt       i,j,Mx,dof,xs,xm,sf = ctx->sf,fs = ctx->fs;
148:   PetscReal      hf,hs,cfl_idt = 0;
149:   PetscScalar    *x,*f,*r,*min,*alpha,*gamma;
150:   Vec            Xloc;
151:   DM             da;

154:   TSGetDM(ts,&da);
155:   DMGetLocalVector(da,&Xloc);                          /* Xloc contains ghost points                                     */
156:   DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);   /* Mx is the number of center points                              */
157:   hs   = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
158:   hf   = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
159:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);       /* X is solution vector which does not contain ghost points       */
160:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
161:   VecZeroEntries(F);                                   /* F is the right hand side function corresponds to center points */
162:   DMDAVecGetArray(da,Xloc,&x);
163:   DMDAVecGetArray(da,F,&f);
164:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
165:   PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);

167:   if (ctx->bctype == FVBC_OUTFLOW) {
168:     for (i=xs-2; i<0; i++) {
169:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
170:     }
171:     for (i=Mx; i<xs+xm+2; i++) {
172:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
173:     }
174:   }

176:   for (i=xs; i<xs+xm+1; i++) {
177:     PetscReal   maxspeed;
178:     PetscScalar *u;
179:     if (i < sf || i > fs+1) {
180:       u = &ctx->u[0];
181:       alpha[0] = 1.0/6.0;
182:       gamma[0] = 1.0/3.0;
183:       for (j=0; j<dof; j++) {
184:         r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
185:         min[j] = PetscMin(r[j],2.0);
186:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
187:       }
188:        (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
189:       cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hs));
190:       if (i > xs) {
191:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
192:       }
193:       if (i < xs+xm) {
194:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
195:       }
196:     } else if (i == sf) {
197:       u = &ctx->u[0];
198:       alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
199:       gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
200:       for (j=0; j<dof; j++) {
201:         r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
202:         min[j] = PetscMin(r[j],2.0);
203:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
204:       }
205:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
206:       if (i > xs) {
207:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
208:       }
209:       if (i < xs+xm) {
210:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hf;
211:       }
212:     } else if (i == sf+1) {
213:       u = &ctx->u[0];
214:       alpha[0] = hf*hf/(hs+hf)/(hs+hf+hf);
215:       gamma[0] = hf*(hs+hf)/(hf+hf)/(hs+hf+hf);
216:       for (j=0; j<dof; j++) {
217:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
218:         min[j] = PetscMin(r[j],2.0);
219:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
220:       }
221:        (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
222:       if (i > xs) {
223:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hf;
224:       }
225:       if (i < xs+xm) {
226:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hf;
227:       }
228:     } else if (i > sf+1 && i < fs) {
229:       u = &ctx->u[0];
230:       alpha[0] = 1.0/6.0;
231:       gamma[0] = 1.0/3.0;
232:       for (j=0; j<dof; j++) {
233:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
234:         min[j] = PetscMin(r[j],2.0);
235:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
236:       }
237:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
238:       if (i > xs) {
239:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hf;
240:       }
241:       if (i < xs+xm) {
242:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hf;
243:       }
244:     } else if (i == fs) {
245:       u = &ctx->u[0];
246:       alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
247:       gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
248:       for (j=0; j<dof; j++) {
249:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
250:         min[j] = PetscMin(r[j],2.0);
251:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
252:       }
253:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
254:       if (i > xs) {
255:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hf;
256:       }
257:       if (i < xs+xm) {
258:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
259:       }
260:     } else if (i == fs+1) {
261:       u = &ctx->u[0];
262:       alpha[0] = hs*hs/(hf+hs)/(hf+hs+hs);
263:       gamma[0] = hs*(hf+hs)/(hs+hs)/(hf+hs+hs);
264:       for (j=0; j<dof; j++) {
265:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
266:         min[j] = PetscMin(r[j],2.0);
267:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
268:       }
269:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
270:       if (i > xs) {
271:         for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hs;
272:       }
273:       if (i < xs+xm) {
274:         for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hs;
275:       }
276:     }
277:   }
278:   DMDAVecRestoreArray(da,Xloc,&x);
279:   DMDAVecRestoreArray(da,F,&f);
280:   DMRestoreLocalVector(da,&Xloc);
281:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da));
282:   if (0) {
283:     /* We need a way to inform the TS of a CFL constraint, this is a debugging fragment */
284:     PetscReal dt,tnow;
285:     TSGetTimeStep(ts,&dt);
286:     TSGetTime(ts,&tnow);
287:     if (dt > 0.5/ctx->cfl_idt) {
288:       PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
289:     }
290:   }
291:   PetscFree4(r,min,alpha,gamma);
292:   return(0);
293:  }

295: static PetscErrorCode FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
296: {
297:   FVCtx             *ctx = (FVCtx*)vctx;
298:   PetscErrorCode    ierr;
299:   PetscInt          i,j,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs;
300:   PetscReal         hf,hs;
301:   PetscScalar       *x,*f,*r,*min,*alpha,*gamma;
302:   Vec               Xloc;
303:   DM                da;

306:   TSGetDM(ts,&da);
307:   DMGetLocalVector(da,&Xloc);                          /* Xloc contains ghost points                                     */
308:   DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);   /* Mx is the number of center points                              */
309:   hs   = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
310:   hf   = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
311:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);       /* X is solution vector which does not contain ghost points       */
312:   DMGlobalToLocalEnd  (da,X,INSERT_VALUES,Xloc);
313:   VecZeroEntries(F);                                   /* F is the right hand side function corresponds to center points */
314:   DMDAVecGetArray(da,Xloc,&x);
315:   VecGetArray(F,&f);
316:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
317:   PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);

319:   if (ctx->bctype == FVBC_OUTFLOW) {
320:     for (i=xs-2; i<0; i++) {
321:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
322:     }
323:     for (i=Mx; i<xs+xm+2; i++) {
324:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
325:     }
326:   }

328:   for (i=xs; i<xs+xm+1; i++) {
329:     PetscReal   maxspeed;
330:     PetscScalar *u;
331:     if (i < sf) {
332:       u = &ctx->u[0];
333:       alpha[0] = 1.0/6.0;
334:       gamma[0] = 1.0/3.0;
335:       for (j=0; j<dof; j++) {
336:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
337:         min[j] = PetscMin(r[j],2.0);
338:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
339:       }
340:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
341:       if (i > xs) {
342:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
343:       }
344:       if (i < xs+xm) {
345:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
346:         islow++;
347:       }
348:     } else if (i == sf) {
349:       u = &ctx->u[0];
350:       alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
351:       gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
352:       for (j=0; j<dof; j++) {
353:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
354:         min[j] = PetscMin(r[j],2.0);
355:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
356:       }
357:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
358:       if (i > xs) {
359:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
360:       }
361:     } else if (i == fs) {
362:       u = &ctx->u[0];
363:       alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
364:       gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
365:       for (j=0; j<dof; j++) {
366:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
367:         min[j] = PetscMin(r[j],2.0);
368:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
369:       }
370:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
371:       if (i < xs+xm) {
372:         for (j=0; j<dof; j++)  f[islow*dof+j] += ctx->flux[j]/hs;
373:         islow++;
374:       }
375:     } else if (i == fs+1) {
376:       u = &ctx->u[0];
377:       alpha[0] = hs*hs/(hf+hs)/(hf+hs+hs);
378:       gamma[0] = hs*(hf+hs)/(hs+hs)/(hf+hs+hs);
379:       for (j=0; j<dof; j++) {
380:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
381:         min[j] = PetscMin(r[j],2.0);
382:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
383:       }
384:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
385:       if (i > xs) {
386:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
387:       }
388:       if (i < xs+xm) {
389:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
390:         islow++;
391:       }
392:     } else if (i > fs+1) {
393:       u = &ctx->u[0];
394:       alpha[0] = 1.0/6.0;
395:       gamma[0] = 1.0/3.0;
396:       for (j=0; j<dof; j++) {
397:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
398:         min[j] = PetscMin(r[j],2.0);
399:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
400:       }
401:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
402:       if (i > xs) {
403:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
404:       }
405:       if (i < xs+xm) {
406:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
407:         islow++;
408:       }
409:     }
410:   }
411:   DMDAVecRestoreArray(da,Xloc,&x);
412:   VecRestoreArray(F,&f);
413:   DMRestoreLocalVector(da,&Xloc);
414:   PetscFree4(r,min,alpha,gamma);
415:   return(0);
416:  }

418: static PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
419: {
420:   FVCtx          *ctx = (FVCtx*)vctx;
422:   PetscInt       i,j,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs;
423:   PetscReal      hf,hs;
424:   PetscScalar    *x,*f,*r,*min,*alpha,*gamma;
425:   Vec            Xloc;
426:   DM             da;

429:   TSGetDM(ts,&da);
430:   DMGetLocalVector(da,&Xloc);                          /* Xloc contains ghost points                                     */
431:   DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);   /* Mx is the number of center points                              */
432:   hs   = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
433:   hf   = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
434:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);       /* X is solution vector which does not contain ghost points       */
435:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
436:   VecZeroEntries(F);                                   /* F is the right hand side function corresponds to center points */
437:   DMDAVecGetArray(da,Xloc,&x);
438:   VecGetArray(F,&f);
439:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
440:   PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);

442:   if (ctx->bctype == FVBC_OUTFLOW) {
443:     for (i=xs-2; i<0; i++) {
444:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
445:     }
446:     for (i=Mx; i<xs+xm+2; i++) {
447:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
448:     }
449:   }

451:   for (i=xs; i<xs+xm+1; i++) {
452:     PetscReal   maxspeed;
453:     PetscScalar *u;
454:     if (i == sf) {
455:       u = &ctx->u[0];
456:       alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf);
457:       gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf);
458:       for (j=0; j<dof; j++) {
459:         r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
460:         min[j] = PetscMin(r[j],2.0);
461:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
462:       }
463:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
464:       if (i < xs+xm) {
465:         for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf;
466:         ifast++;
467:       }
468:     } else if (i == sf+1) {
469:       u = &ctx->u[0];
470:       alpha[0] = hf*hf/(hs+hf)/(hs+hf+hf);
471:       gamma[0] = hf*(hs+hf)/(hf+hf)/(hs+hf+hf);
472:       for (j=0; j<dof; j++) {
473:         r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
474:         min[j] = PetscMin(r[j],2.0);
475:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
476:       }
477:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
478:       if (i > xs) {
479:         for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf;
480:       }
481:       if (i < xs+xm) {
482:         for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf;
483:         ifast++;
484:       }
485:     } else if (i > sf+1 && i < fs) {
486:       u = &ctx->u[0];
487:       alpha[0] = 1.0/6.0;
488:       gamma[0] = 1.0/3.0;
489:       for (j=0; j<dof; j++) {
490:         r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
491:         min[j] = PetscMin(r[j],2.0);
492:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
493:       }
494:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
495:       if (i > xs) {
496:         for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf;
497:       }
498:       if (i < xs+xm) {
499:         for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf;
500:         ifast++;
501:       }
502:     } else if (i == fs) {
503:       u = &ctx->u[0];
504:       alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
505:       gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
506:       for (j=0; j<dof; j++) {
507:         r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
508:         min[j] = PetscMin(r[j],2.0);
509:         u[j] = x[(i-1)*dof+j]+PetscMax(0,PetscMin(min[j],alpha[0]+gamma[0]*r[j]))*(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
510:       }
511:        (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
512:       if (i > xs) {
513:         for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf;
514:       }
515:     }
516:   }
517:   DMDAVecRestoreArray(da,Xloc,&x);
518:   VecRestoreArray(F,&f);
519:   DMRestoreLocalVector(da,&Xloc);
520:   PetscFree4(r,min,alpha,gamma);
521:   return(0);
522:  }

524: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */

526: PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
527: {
529:   PetscScalar    *u,*uj,xj,xi;
530:   PetscInt       i,j,k,dof,xs,xm,Mx,count_slow,count_fast;
531:   const PetscInt N=200;

534:   if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function");
535:   DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
536:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
537:   DMDAVecGetArray(da,U,&u);
538:   PetscMalloc1(dof,&uj);
539:   const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
540:   const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
541:   count_slow = Mx/(1+ctx->hratio);
542:   count_fast = Mx-count_slow;
543:   for (i=xs; i<xs+xm; i++) {
544:     if (i*hs+0.5*hs<(ctx->xmax-ctx->xmin)*0.25) {
545:       xi = ctx->xmin+0.5*hs+i*hs;
546:       /* Integrate over cell i using trapezoid rule with N points. */
547:       for (k=0; k<dof; k++) u[i*dof+k] = 0;
548:       for (j=0; j<N+1; j++) {
549:         xj = xi+hs*(j-N/2)/(PetscReal)N;
550:         (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
551:         for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
552:       }
553:     } else if ((ctx->xmax-ctx->xmin)*0.25+(i-count_slow/2)*hf+0.5*hf<(ctx->xmax-ctx->xmin)*0.75) {
554:       xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.25+0.5*hf+(i-count_slow/2)*hf;
555:       /* Integrate over cell i using trapezoid rule with N points. */
556:       for (k=0; k<dof; k++) u[i*dof+k] = 0;
557:       for (j=0; j<N+1; j++) {
558:         xj = xi+hf*(j-N/2)/(PetscReal)N;
559:         (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
560:         for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
561:       }
562:     } else {
563:       xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.75+0.5*hs+(i-count_slow/2-count_fast)*hs;
564:       /* Integrate over cell i using trapezoid rule with N points. */
565:       for (k=0; k<dof; k++) u[i*dof+k] = 0;
566:       for (j=0; j<N+1; j++) {
567:         xj = xi+hs*(j-N/2)/(PetscReal)N;
568:         (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
569:         for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
570:       }
571:     }
572:   }
573:   DMDAVecRestoreArray(da,U,&u);
574:   PetscFree(uj);
575:   return(0);
576: }

578: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
579: {
580:   PetscErrorCode    ierr;
581:   PetscReal         xmin,xmax;
582:   PetscScalar       sum,tvsum,tvgsum;
583:   const PetscScalar *x;
584:   PetscInt          imin,imax,Mx,i,j,xs,xm,dof;
585:   Vec               Xloc;
586:   PetscBool         iascii;

589:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
590:   if (iascii) {
591:     /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
592:     DMGetLocalVector(da,&Xloc);
593:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
594:     DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
595:     DMDAVecGetArrayRead(da,Xloc,(void*)&x);
596:     DMDAGetCorners(da,&xs,0,0,&xm,0,0);
597:     DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
598:     tvsum = 0;
599:     for (i=xs; i<xs+xm; i++) {
600:       for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j]-x[(i-1)*dof+j]);
601:     }
602:     MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)da));
603:     DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);
604:     DMRestoreLocalVector(da,&Xloc);

606:     VecMin(X,&imin,&xmin);
607:     VecMax(X,&imax,&xmax);
608:     VecSum(X,&sum);
609:     PetscViewerASCIIPrintf(viewer,"Solution range [%g,%g] with minimum at %D, mean %g, ||x||_TV %g\n",(double)xmin,(double)xmax,imin,(double)(sum/Mx),(double)(tvgsum/Mx));
610:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported");
611:   return(0);
612: }

614: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1)
615: {
616:   PetscErrorCode    ierr;
617:   Vec               Y;
618:   PetscInt          i,Mx,count_slow=0,count_fast=0;
619:   const PetscScalar *ptr_X,*ptr_Y;

622:   VecGetSize(X,&Mx);
623:   VecDuplicate(X,&Y);
624:   FVSample(ctx,da,t,Y);
625:   const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
626:   const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
627:   count_slow = (PetscReal)Mx/(1.0+ctx->hratio);
628:   count_fast = Mx-count_slow;
629:   VecGetArrayRead(X,&ptr_X);
630:   VecGetArrayRead(Y,&ptr_Y);
631:   for (i=0; i<Mx; i++) {
632:     if (i < count_slow/2 || i > count_slow/2+count_fast-1) *nrm1 +=  hs*PetscAbs(ptr_X[i]-ptr_Y[i]);
633:     else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]);
634:   }
635:   VecRestoreArrayRead(X,&ptr_X);
636:   VecRestoreArrayRead(Y,&ptr_Y);
637:   VecDestroy(&Y);
638:   return(0);
639: }

641: int main(int argc,char *argv[])
642: {
643:   char              physname[256] = "advect",final_fname[256] = "solution.m";
644:   PetscFunctionList physics = 0;
645:   MPI_Comm          comm;
646:   TS                ts;
647:   DM                da;
648:   Vec               X,X0,R;
649:   FVCtx             ctx;
650:   PetscInt          i,k,dof,xs,xm,Mx,draw = 0,count_slow,count_fast,islow = 0,ifast = 0,*index_slow,*index_fast;
651:   PetscBool         view_final = PETSC_FALSE;
652:   PetscReal         ptime;
653:   PetscErrorCode    ierr;

655:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
656:   comm = PETSC_COMM_WORLD;
657:   PetscMemzero(&ctx,sizeof(ctx));

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

662:   ctx.comm = comm;
663:   ctx.cfl  = 0.9;
664:   ctx.bctype = FVBC_PERIODIC;
665:   ctx.xmin = -1.0;
666:   ctx.xmax = 1.0;
667:   PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
668:   PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
669:   PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
670:   PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
671:   PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
672:   PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
673:   PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
674:   PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);
675:   PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
676:   PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
677:   PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL);
678:   PetscOptionsEnd();

680:   /* Choose the physics from the list of registered models */
681:   {
682:     PetscErrorCode (*r)(FVCtx*);
683:     PetscFunctionListFind(physics,physname,&r);
684:     if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Physics '%s' not found",physname);
685:     /* Create the physics, will set the number of fields and their names */
686:     (*r)(&ctx);
687:   }

689:   /* Create a DMDA to manage the parallel grid */
690:   DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);
691:   DMSetFromOptions(da);
692:   DMSetUp(da);
693:   /* Inform the DMDA of the field names provided by the physics. */
694:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
695:   for (i=0; i<ctx.physics.dof; i++) {
696:     DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
697:   }
698:   DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);
699:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

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

704:   /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
705:   PetscMalloc3(dof,&ctx.u,dof,&ctx.flux,dof,&ctx.speeds);

707:   /* Create a vector to store the solution and to save the initial state */
708:   DMCreateGlobalVector(da,&X);
709:   VecDuplicate(X,&X0);
710:   VecDuplicate(X,&R);

712:   /* create index for slow parts and fast parts*/
713:   count_slow = Mx/(1+ctx.hratio);
714:   if (count_slow%2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio) is even");
715:   count_fast = Mx-count_slow;
716:   ctx.sf = count_slow/2;
717:   ctx.fs = ctx.sf + count_fast;
718:   PetscMalloc1(xm*dof,&index_slow);
719:   PetscMalloc1(xm*dof,&index_fast);
720:   for (i=xs; i<xs+xm; i++) {
721:     if (i < count_slow/2 || i > count_slow/2+count_fast-1)
722:       for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
723:     else
724:       for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
725:   }
726:   ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
727:   ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);

729:   /* Create a time-stepping object */
730:   TSCreate(comm,&ts);
731:   TSSetDM(ts,da);
732:   TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
733:   TSRHSSplitSetIS(ts,"slow",ctx.iss);
734:   TSRHSSplitSetIS(ts,"fast",ctx.isf);
735:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);
736:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);

738:   TSSetType(ts,TSMPRK);
739:   TSSetMaxTime(ts,10);
740:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

742:   /* Compute initial conditions and starting time step */
743:   FVSample(&ctx,da,0,X0);
744:   FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
745:   VecCopy(X0,X);                        /* The function value was not used so we set X=X0 again */
746:   TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
747:   TSSetFromOptions(ts); /* Take runtime options */
748:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
749:   {
750:     PetscInt          steps;
751:     PetscScalar       mass_initial,mass_final,mass_difference,mass_differenceg;
752:     const PetscScalar *ptr_X,*ptr_X0;
753:     const PetscReal   hs  = (ctx.xmax-ctx.xmin)/2.0/count_slow;
754:     const PetscReal   hf  = (ctx.xmax-ctx.xmin)/2.0/count_fast;
755:     TSSolve(ts,X);
756:     TSGetSolveTime(ts,&ptime);
757:     TSGetStepNumber(ts,&steps);
758:     /* calculate the total mass at initial time and final time */
759:     mass_initial = 0.0;
760:     mass_final   = 0.0;
761:     DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0);
762:     DMDAVecGetArrayRead(da,X,(void*)&ptr_X);
763:     for (i=xs; i<xs+xm; i++) {
764:       if (i < ctx.sf || i > ctx.fs-1) {
765:         for (k=0; k<dof; k++) {
766:           mass_initial = mass_initial+hs*ptr_X0[i*dof+k];
767:           mass_final = mass_final+hs*ptr_X[i*dof+k];
768:         }
769:       } else {
770:         for (k=0; k<dof; k++) {
771:           mass_initial = mass_initial+hf*ptr_X0[i*dof+k];
772:           mass_final = mass_final+hf*ptr_X[i*dof+k];
773:         }
774:       }
775:     }
776:     DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0);
777:     DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X);
778:     mass_difference = mass_final-mass_initial;
779:     MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm);
780:     PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg);
781:     PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);
782:     if (ctx.exact) {
783:       PetscReal nrm1 = 0;
784:       SolutionErrorNorms(&ctx,da,ptime,X,&nrm1);
785:       PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
786:     }
787:     if (ctx.simulation) {
788:       PetscReal         nrm1 = 0;
789:       PetscViewer       fd;
790:       char              filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
791:       Vec               XR;
792:       PetscBool         flg;
793:       const PetscScalar *ptr_XR;
794:       PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);
795:       if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
796:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
797:       VecDuplicate(X0,&XR);
798:       VecLoad(XR,fd);
799:       PetscViewerDestroy(&fd);
800:       VecGetArrayRead(X,&ptr_X);
801:       VecGetArrayRead(XR,&ptr_XR);
802:       for (i=0; i<Mx; i++) {
803:         if (i < count_slow/2 || i > count_slow/2+count_fast-1) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i]-ptr_XR[i]);
804:         else nrm1 = nrm1 + hf*PetscAbs(ptr_X[i]-ptr_XR[i]);
805:       }
806:       VecRestoreArrayRead(X,&ptr_X);
807:       VecRestoreArrayRead(XR,&ptr_XR);
808:       PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);
809:       VecDestroy(&XR);
810:     }
811:   }

813:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
814:   if (draw & 0x1) { VecView(X0,PETSC_VIEWER_DRAW_WORLD); }
815:   if (draw & 0x2) { VecView(X,PETSC_VIEWER_DRAW_WORLD); }
816:   if (draw & 0x4) {
817:     Vec Y;
818:     VecDuplicate(X,&Y);
819:     FVSample(&ctx,da,ptime,Y);
820:     VecAYPX(Y,-1,X);
821:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
822:     VecDestroy(&Y);
823:   }

825:   if (view_final) {
826:     PetscViewer viewer;
827:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
828:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
829:     VecView(X,viewer);
830:     PetscViewerPopFormat(viewer);
831:     PetscViewerDestroy(&viewer);
832:   }

834:   /* Clean up */
835:   (*ctx.physics.destroy)(ctx.physics.user);
836:   for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
837:   PetscFree3(ctx.u,ctx.flux,ctx.speeds);
838:   ISDestroy(&ctx.iss);
839:   ISDestroy(&ctx.isf);
840:   VecDestroy(&X);
841:   VecDestroy(&X0);
842:   VecDestroy(&R);
843:   DMDestroy(&da);
844:   TSDestroy(&ts);
845:   PetscFree(index_slow);
846:   PetscFree(index_fast);
847:   PetscFunctionListDestroy(&physics);
848:   PetscFinalize();
849:   return ierr;
850: }

852: /*TEST

854:     build:
855:       requires: !complex

857:     test:
858:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 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

860:     test:
861:       suffix: 2
862:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 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
863:       output_file: output/ex7_1.out

865:     test:
866:       suffix: 3
867:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0

869:     test:
870:       suffix: 4
871:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
872:       output_file: output/ex7_3.out

874:     test:
875:       suffix: 5
876:       nsize: 2
877:       args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
878:       output_file: output/ex7_3.out
879: TEST*/