Actual source code: ex7.c

petsc-3.13.4 2020-08-01
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,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:       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[(islow-1)*dof+j] -= ctx->flux[j]/hs;
345:       }
346:       if (i < xs+xm) {
347:         for (j=0; j<dof; j++) f[islow*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) {
361:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
362:       }
363:     } else if (i == fs) {
364:       u = &ctx->u[0];
365:       alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs);
366:       gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs);
367:       for (j=0; j<dof; j++) {
368:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
369:         min[j] = PetscMin(r[j],2.0);
370:         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]);
371:       }
372:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
373:       if (i < xs+xm) {
374:         for (j=0; j<dof; j++)  f[islow*dof+j] += ctx->flux[j]/hs;
375:         islow++;
376:       }
377:     } else if (i == fs+1) {
378:       u = &ctx->u[0];
379:       alpha[0] = hs*hs/(hf+hs)/(hf+hs+hs);
380:       gamma[0] = hs*(hf+hs)/(hs+hs)/(hf+hs+hs);
381:       for (j=0; j<dof; j++) {
382:         r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]);
383:         min[j] = PetscMin(r[j],2.0);
384:         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]);
385:       }
386:       (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);
387:       if (i > xs) {
388:         for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs;
389:       }
390:       if (i < xs+xm) {
391:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
392:         islow++;
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[(islow-1)*dof+j] -= ctx->flux[j]/hs;
406:       }
407:       if (i < xs+xm) {
408:         for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs;
409:         islow++;
410:       }
411:     }
412:   }
413:   DMDAVecRestoreArray(da,Xloc,&x);
414:   VecRestoreArray(F,&f);
415:   DMRestoreLocalVector(da,&Xloc);
416:   PetscFree4(r,min,alpha,gamma);
417:   return(0);
418:  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

854: /*TEST

856:     build:
857:       requires: !complex c99

859:     test:
860:       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

862:     test:
863:       suffix: 2
864:       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
865:       output_file: output/ex7_1.out

867:     test:
868:       suffix: 3
869:       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

871:     test:
872:       suffix: 4
873:       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
874:       output_file: output/ex7_3.out

876:     test:
877:       suffix: 5
878:       nsize: 2
879:       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
880:       output_file: output/ex7_3.out
881: TEST*/