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: 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: {
 71:   PetscFree(vctx);
 72:   return 0;
 73: }

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

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

 86:   speed     = ctx->a;
 87:   flux[0]   = speed*u[0];
 88:   *maxspeed = speed;
 89:   return 0;
 90: }

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

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

117: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
118: {
120:   AdvectCtx      *user;

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

139: /* --------------------------------- Finite Volume Solver ----------------------------------- */

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

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

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

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

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

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

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

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

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

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

437:   if (ctx->bctype == FVBC_OUTFLOW) {
438:     for (i=xs-2; i<0; i++) {
439:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
440:     }
441:     for (i=Mx; i<xs+xm+2; i++) {
442:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
443:     }
444:   }

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

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

521: PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
522: {
523:   PetscScalar    *u,*uj,xj,xi;
524:   PetscInt       i,j,k,dof,xs,xm,Mx,count_slow,count_fast;
525:   const PetscInt N=200;

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

572: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
573: {
574:   PetscReal         xmin,xmax;
575:   PetscScalar       sum,tvsum,tvgsum;
576:   const PetscScalar *x;
577:   PetscInt          imin,imax,Mx,i,j,xs,xm,dof;
578:   Vec               Xloc;
579:   PetscBool         iascii;

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

599:     VecMin(X,&imin,&xmin);
600:     VecMax(X,&imax,&xmax);
601:     VecSum(X,&sum);
602:     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));
603:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported");
604:   return 0;
605: }

607: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1)
608: {
609:   Vec               Y;
610:   PetscInt          i,Mx,count_slow=0,count_fast=0;
611:   const PetscScalar *ptr_X,*ptr_Y;

614:   VecGetSize(X,&Mx);
615:   VecDuplicate(X,&Y);
616:   FVSample(ctx,da,t,Y);
617:   const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx;
618:   const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx;
619:   count_slow = (PetscReal)Mx/(1.0+ctx->hratio);
620:   count_fast = Mx-count_slow;
621:   VecGetArrayRead(X,&ptr_X);
622:   VecGetArrayRead(Y,&ptr_Y);
623:   for (i=0; i<Mx; i++) {
624:     if (i < count_slow/2 || i > count_slow/2+count_fast-1) *nrm1 +=  hs*PetscAbs(ptr_X[i]-ptr_Y[i]);
625:     else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]);
626:   }
627:   VecRestoreArrayRead(X,&ptr_X);
628:   VecRestoreArrayRead(Y,&ptr_Y);
629:   VecDestroy(&Y);
630:   return 0;
631: }

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

647:   PetscInitialize(&argc,&argv,0,help);
648:   comm = PETSC_COMM_WORLD;
649:   PetscMemzero(&ctx,sizeof(ctx));

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

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

672:   /* Choose the physics from the list of registered models */
673:   {
674:     PetscErrorCode (*r)(FVCtx*);
675:     PetscFunctionListFind(physics,physname,&r);
677:     /* Create the physics, will set the number of fields and their names */
678:     (*r)(&ctx);
679:   }

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

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

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

699:   /* Create a vector to store the solution and to save the initial state */
700:   DMCreateGlobalVector(da,&X);
701:   VecDuplicate(X,&X0);
702:   VecDuplicate(X,&R);

704:   /* create index for slow parts and fast parts*/
705:   count_slow = Mx/(1+ctx.hratio);
707:   count_fast = Mx-count_slow;
708:   ctx.sf = count_slow/2;
709:   ctx.fs = ctx.sf + count_fast;
710:   PetscMalloc1(xm*dof,&index_slow);
711:   PetscMalloc1(xm*dof,&index_fast);
712:   for (i=xs; i<xs+xm; i++) {
713:     if (i < count_slow/2 || i > count_slow/2+count_fast-1)
714:       for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k;
715:     else
716:       for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k;
717:   }
718:   ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);
719:   ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);

721:   /* Create a time-stepping object */
722:   TSCreate(comm,&ts);
723:   TSSetDM(ts,da);
724:   TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
725:   TSRHSSplitSetIS(ts,"slow",ctx.iss);
726:   TSRHSSplitSetIS(ts,"fast",ctx.isf);
727:   TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);
728:   TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);

730:   TSSetType(ts,TSMPRK);
731:   TSSetMaxTime(ts,10);
732:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

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

805:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
806:   if (draw & 0x1) VecView(X0,PETSC_VIEWER_DRAW_WORLD);
807:   if (draw & 0x2) VecView(X,PETSC_VIEWER_DRAW_WORLD);
808:   if (draw & 0x4) {
809:     Vec Y;
810:     VecDuplicate(X,&Y);
811:     FVSample(&ctx,da,ptime,Y);
812:     VecAYPX(Y,-1,X);
813:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
814:     VecDestroy(&Y);
815:   }

817:   if (view_final) {
818:     PetscViewer viewer;
819:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
820:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
821:     VecView(X,viewer);
822:     PetscViewerPopFormat(viewer);
823:     PetscViewerDestroy(&viewer);
824:   }

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

844: /*TEST

846:     build:
847:       requires: !complex

849:     test:
850:       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

852:     test:
853:       suffix: 2
854:       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
855:       output_file: output/ex7_1.out

857:     test:
858:       suffix: 3
859:       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

861:     test:
862:       suffix: 4
863:       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
864:       output_file: output/ex7_3.out

866:     test:
867:       suffix: 5
868:       nsize: 2
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 1
870:       output_file: output/ex7_3.out
871: TEST*/