Actual source code: ex7.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  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,1,"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,1,"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 to 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:       if (1) {
289:         PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
290:       } else SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt));
291:     }
292:   }
293:   PetscFree4(r,min,alpha,gamma);
294:   return(0);
295:  }

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

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

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

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

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

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

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

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

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

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

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

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

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

607:     VecMin(X,&imin,&xmin);
608:     VecMax(X,&imax,&xmax);
609:     VecSum(X,&sum);
610:     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));
611:   } else SETERRQ(PETSC_COMM_SELF,1,"Viewer type not supported");
612:   return(0);
613: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

849: /*TEST

851:     build:
852:       requires: !complex c99

854:     test:
855:       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

857:     test:
858:       suffix: 2
859:       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
860:       output_file: output/ex7_1.out

862:     test:
863:       suffix: 3
864:       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

866:     test:
867:       suffix: 4
868:       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
869:       output_file: output/ex7_3.out
870: TEST*/