Actual source code: ex9.c

  1: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
  2:   "Solves scalar and vector problems, choose the physical model with -physics\n"
  3:   "  advection   - Constant coefficient scalar advection\n"
  4:   "                u_t       + (a*u)_x               = 0\n"
  5:   "  burgers     - Burgers equation\n"
  6:   "                u_t       + (u^2/2)_x             = 0\n"
  7:   "  traffic     - Traffic equation\n"
  8:   "                u_t       + (u*(1-u))_x           = 0\n"
  9:   "  isogas      - Isothermal gas dynamics\n"
 10:   "                rho_t     + (rho*u)_x             = 0\n"
 11:   "                (rho*u)_t + (rho*u^2 + c^2*rho)_x = 0\n"
 12:   "  shallow     - Shallow water equations\n"
 13:   "                h_t       + (h*u)_x               = 0\n"
 14:   "                (h*u)_t   + (h*u^2 + g*h^2/2)_x   = 0\n"
 15:   "Some of these physical models have multiple Riemann solvers, select these with -physics_xxx_riemann\n"
 16:   "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
 17:   "                the states across shocks and rarefactions\n"
 18:   "  roe         - Linearized scheme, usually with an entropy fix inside sonic rarefactions\n"
 19:   "The systems provide a choice of reconstructions with -physics_xxx_reconstruct\n"
 20:   "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
 21:   "  conservative   - Limit the conservative variables directly, can cause undesired interaction of waves\n\n"
 22:   "A variety of limiters for high-resolution TVD limiters are available with -limit\n"
 23:   "  upwind,minmod,superbee,mc,vanleer,vanalbada,koren,cada-torillhon (last two are nominally third order)\n"
 24:   "  and non-TVD schemes lax-wendroff,beam-warming,fromm\n\n"
 25:   "To preserve the TVD property, one should time step with a strong stability preserving method.\n"
 26:   "The optimal high order explicit Runge-Kutta methods in TSSSP are recommended for non-stiff problems.\n\n"
 27:   "Several initial conditions can be chosen with -initial N\n\n"
 28:   "The problem size should be set with -da_grid_x M\n\n";

 30: /* To get isfinite in math.h */
 31: #define _XOPEN_SOURCE 600

 33: #include <petscts.h>
 34: #include <petscdmda.h>

 36: #include <../src/mat/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */

 38: static inline PetscReal Sgn(PetscReal a) { return (a<0) ? -1 : 1; }
 39: static inline PetscReal Abs(PetscReal a) { return (a<0) ? 0 : a; }
 40: static inline PetscReal Sqr(PetscReal a) { return a*a; }
 41: static inline PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; }
 42: static inline PetscReal MinAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) < PetscAbs(b)) ? a : b; }
 43: static inline PetscReal MinMod2(PetscReal a,PetscReal b)
 44: { return (a*b<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscAbs(b)); }
 45: static inline PetscReal MaxMod2(PetscReal a,PetscReal b)
 46: { return (a*b<0) ? 0 : Sgn(a)*PetscMax(PetscAbs(a),PetscAbs(b)); }
 47: static inline PetscReal MinMod3(PetscReal a,PetscReal b,PetscReal c)
 48: {return (a*b<0 || a*c<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscMin(PetscAbs(b),PetscAbs(c))); }

 50: static inline PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax)
 51: { PetscReal range = xmax-xmin; return xmin + fmod(range+fmod(a,range),range); }


 54: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
 55: typedef struct _LimitInfo {
 56:   PetscReal hx;
 57:   PetscInt m;
 58: } *LimitInfo;
 59: static void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
 60: {
 61:   PetscInt i;
 62:   for (i=0; i<info->m; i++) lmt[i] = 0;
 63: }
 64: static void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
 65: {
 66:   PetscInt i;
 67:   for (i=0; i<info->m; i++) lmt[i] = jR[i];
 68: }
 69: static void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
 70: {
 71:   PetscInt i;
 72:   for (i=0; i<info->m; i++) lmt[i] = jL[i];
 73: }
 74: static void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
 75: {
 76:   PetscInt i;
 77:   for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i] + jR[i]);
 78: }
 79: static void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
 80: {
 81:   PetscInt i;
 82:   for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]);
 83: }
 84: static void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
 85: {
 86:   PetscInt i;
 87:   for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]));
 88: }
 89: static void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
 90: {
 91:   PetscInt i;
 92:   for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]);
 93: }
 94: static void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
 95: { /* phi = (t + abs(t)) / (1 + abs(t)) */
 96:   PetscInt i;
 97:   for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Abs(jR[i]) + Abs(jL[i])*jR[i]) / (Abs(jL[i]) + Abs(jR[i]) + 1e-15);
 98: }
 99: static void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
100: { /* phi = (t + t^2) / (1 + t^2) */
101:   PetscInt i;
102:   for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
103: }
104: static void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
105: { /* phi = (t + t^2) / (1 + t^2) */
106:   PetscInt i;
107:   for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0
108:                         : (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
109: }
110: static void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
111: { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
112:   PetscInt i;
113:   for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i]) + 2*Sqr(jL[i])*jR[i])
114:                                 / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
115: }
116: static void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
117: { /* Symmetric version of above */
118:   PetscInt i;
119:   for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i])
120:                                 / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
121: }
122: static void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
123: { /* Eq 11 of Cada-Torrilhon 2009 */
124:   PetscInt i;
125:   for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]);
126: }

128: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)
129: { return PetscMax(0,PetscMin((L+2*R)/3,
130:                               PetscMax(-0.5*L,
131:                                        PetscMin(2*L,
132:                                                 PetscMin((L+2*R)/3,1.6*R)))));
133: }
134: static void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
135: { /* Cada-Torrilhon 2009, Eq 13 */
136:   PetscInt i;
137:   for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]);
138: }
139: static void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
140: { /* Cada-Torrilhon 2009, Eq 22 */
141:   /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
142:   const PetscReal eps = 1e-7,hx = info->hx;
143:   PetscInt i;
144:   for (i=0; i<info->m; i++) {
145:     const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r*hx);
146:     lmt[i] = ((eta < 1-eps)
147:               ? (jL[i] + 2*jR[i]) / 3
148:               : ((eta > 1+eps)
149:                  ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i])
150:                  : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3
151:                         + (1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]))));
152:   }
153: }
154: static void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
155: { Limit_CadaTorrilhon3R(0.1,info,jL,jR,lmt); }
156: static void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
157: { Limit_CadaTorrilhon3R(1,info,jL,jR,lmt); }
158: static void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
159: { Limit_CadaTorrilhon3R(10,info,jL,jR,lmt); }
160: static void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
161: { Limit_CadaTorrilhon3R(100,info,jL,jR,lmt); }


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

166: typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType;
167: static const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0};
168: typedef PetscErrorCode (*RiemannFunction)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*);
169: typedef PetscErrorCode (*ReconstructFunction)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*);

171: typedef struct {
172:   PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*);
173:   RiemannFunction riemann;
174:   ReconstructFunction characteristic;
175:   PetscErrorCode (*destroy)(void*);
176:   void *user;
177:   PetscInt dof;
178:   char *fieldname[16];
179: } PhysicsCtx;

181: typedef struct {
182:   void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*);
183:   PhysicsCtx physics;

185:   MPI_Comm comm;
186:   char prefix[256];
187:   DM da;
188:   /* Local work arrays */
189:   PetscScalar *R,*Rinv;         /* Characteristic basis, and it's inverse.  COLUMN-MAJOR */
190:   PetscScalar *cjmpLR;          /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
191:   PetscScalar *cslope;          /* Limited slope, written in characteristic basis */
192:   PetscScalar *uLR;             /* Solution at left and right of interface, conservative variables, len=2*dof */
193:   PetscScalar *flux;            /* Flux across interface */

195:   PetscReal cfl_idt;            /* Max allowable value of 1/Delta t */
196:   PetscReal cfl;
197:   PetscReal xmin,xmax;
198:   PetscInt initial;
199:   PetscBool  exact;
200:   FVBCType bctype;
201: } FVCtx;


204: /* Utility */

208: PetscErrorCode RiemannListAdd(PetscFList *flist,const char *name,RiemannFunction rsolve)
209: {

213:   PetscFListAdd(flist,name,"",(void(*)(void))rsolve);
214:   return(0);
215: }

219: PetscErrorCode RiemannListFind(PetscFList flist,const char *name,RiemannFunction *rsolve)
220: {

224:   PetscFListFind(flist,PETSC_COMM_WORLD,name,PETSC_FALSE,(void(**)(void))rsolve);
225:   if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,1,"Riemann solver \"%s\" could not be found",name);
226:   return(0);
227: }

231: PetscErrorCode ReconstructListAdd(PetscFList *flist,const char *name,ReconstructFunction r)
232: {

236:   PetscFListAdd(flist,name,"",(void(*)(void))r);
237:   return(0);
238: }

242: PetscErrorCode ReconstructListFind(PetscFList flist,const char *name,ReconstructFunction *r)
243: {

247:   PetscFListFind(flist,PETSC_COMM_WORLD,name,PETSC_FALSE,(void(**)(void))r);
248:   if (!*r) SETERRQ1(PETSC_COMM_SELF,1,"Reconstruction \"%s\" could not be found",name);
249:   return(0);
250: }


253: /* --------------------------------- Physics ----------------------------------- */
254: /**
255: * Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction.  These
256: * are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the
257: * number of fields and their names, and a function to deallocate private storage.
258: **/

260: /* First a few functions useful to several different physics */
263: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi)
264: {
265:   PetscInt i,j;

268:   for (i=0; i<m; i++) {
269:     for (j=0; j<m; j++) {
270:       Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j);
271:     }
272:   }
273:   return(0);
274: }

278: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
279: {

283:   PetscFree(vctx);
284:   return(0);
285: }



289: /* --------------------------------- Advection ----------------------------------- */

291: typedef struct {
292:   PetscReal a;                  /* advective velocity */
293: } AdvectCtx;

297: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
298: {
299:   AdvectCtx *ctx = (AdvectCtx*)vctx;
300:   PetscReal speed;

303:   speed = ctx->a;
304:   flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
305:   *maxspeed = speed;
306:   return(0);
307: }

311: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
312: {
313:   AdvectCtx *ctx = (AdvectCtx*)vctx;
314:   PetscReal a = ctx->a,x0;

317:   switch (bctype) {
318:     case FVBC_OUTFLOW: x0 = x-a*t; break;
319:     case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
320:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
321:   }
322:   switch (initial) {
323:     case 0: u[0] = (x0 < 0) ? 1 : -1; break;
324:     case 1: u[0] = (x0 < 0) ? -1 : 1; break;
325:     case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
326:     case 3: u[0] = sin(2*M_PI*x0); break;
327:     case 4: u[0] = PetscAbs(x0); break;
328:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
329:   }
330:   return(0);
331: }

335: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
336: {
338:   AdvectCtx *user;

341:   PetscNew(*user,&user);
342:   ctx->physics.sample         = PhysicsSample_Advect;
343:   ctx->physics.riemann        = PhysicsRiemann_Advect;
344:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
345:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
346:   ctx->physics.user           = user;
347:   ctx->physics.dof            = 1;
348:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
349:   user->a = 1;
350:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
351:   {
352:     PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,PETSC_NULL);
353:   }
354:   PetscOptionsEnd();
355:   return(0);
356: }



360: /* --------------------------------- Burgers ----------------------------------- */

362: typedef struct {
363:   PetscReal lxf_speed;
364: } BurgersCtx;

368: static PetscErrorCode PhysicsSample_Burgers(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
369: {

372:   if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
373:   switch (initial) {
374:     case 0: u[0] = (x < 0) ? 1 : -1; break;
375:     case 1:
376:       if       (x < -t) u[0] = -1;
377:       else if  (x < t)  u[0] = x/t;
378:       else              u[0] = 1;
379:       break;
380:     case 2:
381:       if      (x < 0)       u[0] = 0;
382:       else if (x <= t)      u[0] = x/t;
383:       else if (x < 1+0.5*t) u[0] = 1;
384:       else                  u[0] = 0;
385:       break;
386:     case 3:
387:       if       (x < 0.2*t) u[0] = 0.2;
388:       else if  (x < t) u[0] = x/t;
389:       else             u[0] = 1;
390:       break;
391:     case 4:
392:       if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
393:       u[0] = 0.7 + 0.3*sin(2*M_PI*((x-xmin)/(xmax-xmin)));
394:       break;
395:     case 5:                     /* Pure shock solution */
396:       if (x < 0.5*t) u[0] = 1;
397:       else u[0] = 0;
398:       break;
399:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
400:   }
401:   return(0);
402: }

406: static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
407: {

410:   if (uL[0] < uR[0]) {          /* rarefaction */
411:     flux[0] = (uL[0]*uR[0] < 0)
412:       ? 0                       /* sonic rarefaction */
413:       : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0]));
414:   } else {                      /* shock */
415:     flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0]));
416:   }
417:   *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
418:   return(0);
419: }

423: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
424: {
425:   PetscReal speed;

428:   speed = 0.5*(uL[0] + uR[0]);
429:   flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
430:   if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
431:   *maxspeed = speed;
432:   return(0);
433: }

437: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
438: {
439:   PetscReal c;
440:   PetscScalar fL,fR;

443:   c = ((BurgersCtx*)vctx)->lxf_speed;
444:   fL = 0.5*PetscSqr(uL[0]);
445:   fR = 0.5*PetscSqr(uR[0]);
446:   flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
447:   *maxspeed = c;
448:   return(0);
449: }

453: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
454: {
455:   PetscReal c;
456:   PetscScalar fL,fR;

459:   c = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0]));
460:   fL = 0.5*PetscSqr(uL[0]);
461:   fR = 0.5*PetscSqr(uR[0]);
462:   flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
463:   *maxspeed = c;
464:   return(0);
465: }

469: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
470: {
471:   BurgersCtx *user;
473:   RiemannFunction r;
474:   PetscFList rlist = 0;
475:   char rname[256] = "exact";

478:   PetscNew(*user,&user);
479:   ctx->physics.sample         = PhysicsSample_Burgers;
480:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
481:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
482:   ctx->physics.user           = user;
483:   ctx->physics.dof            = 1;
484:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
485:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Burgers_Exact);
486:   RiemannListAdd(&rlist,"roe",    PhysicsRiemann_Burgers_Roe);
487:   RiemannListAdd(&rlist,"lxf",    PhysicsRiemann_Burgers_LxF);
488:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Burgers_Rusanov);
489:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
490:   {
491:     PetscOptionsList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
492:   }
493:   PetscOptionsEnd();
494:   RiemannListFind(rlist,rname,&r);
495:   PetscFListDestroy(&rlist);
496:   ctx->physics.riemann = r;

498:   /* *
499:   * Hack to deal with LxF in semi-discrete form
500:   * max speed is 1 for the basic initial conditions (where |u| <= 1)
501:   * */
502:   if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
503:   return(0);
504: }



508: /* --------------------------------- Traffic ----------------------------------- */

510: typedef struct {
511:   PetscReal lxf_speed;
512:   PetscReal a;
513: } TrafficCtx;

515: static inline PetscScalar TrafficFlux(PetscScalar a,PetscScalar u) { return a*u*(1-u); }

519: static PetscErrorCode PhysicsSample_Traffic(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
520: {
521:   PetscReal a = ((TrafficCtx*)vctx)->a;

524:   if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
525:   switch (initial) {
526:     case 0:
527:       u[0] = (-a*t < x) ? 2 : 0; break;
528:     case 1:
529:       if      (x < PetscMin(2*a*t,0.5+a*t)) u[0] = -1;
530:       else if (x < 1)                       u[0] = 0;
531:       else                                  u[0] = 1;
532:       break;
533:     case 2:
534:       if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
535:       u[0] = 0.7 + 0.3*sin(2*M_PI*((x-xmin)/(xmax-xmin)));
536:       break;
537:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
538:   }
539:   return(0);
540: }

544: static PetscErrorCode PhysicsRiemann_Traffic_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
545: {
546:   PetscReal a = ((TrafficCtx*)vctx)->a;

549:   if (uL[0] < uR[0]) {
550:     flux[0] = PetscMin(TrafficFlux(a,uL[0]),
551:                        TrafficFlux(a,uR[0]));
552:   } else {
553:     flux[0] = (uR[0] < 0.5 && 0.5 < uL[0])
554:       ? TrafficFlux(a,0.5)
555:       : PetscMax(TrafficFlux(a,uL[0]),
556:                  TrafficFlux(a,uR[0]));
557:   }
558:   *maxspeed = a*MaxAbs(1-2*uL[0],1-2*uR[0]);
559:   return(0);
560: }

564: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
565: {
566:   PetscReal a = ((TrafficCtx*)vctx)->a;
567:   PetscReal speed;

570:   speed = a*(1 - (uL[0] + uR[0]));
571:   flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
572:   *maxspeed = speed;
573:   return(0);
574: }

578: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
579: {
580:   TrafficCtx *phys = (TrafficCtx*)vctx;
581:   PetscReal a = phys->a;
582:   PetscReal speed;

585:   speed = a*(1 - (uL[0] + uR[0]));
586:   flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]);
587:   *maxspeed = speed;
588:   return(0);
589: }

593: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
594: {
595:   PetscReal a = ((TrafficCtx*)vctx)->a;
596:   PetscReal speed;

599:   speed = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0]));
600:   flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]);
601:   *maxspeed = speed;
602:   return(0);
603: }

607: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
608: {
610:   TrafficCtx *user;
611:   RiemannFunction r;
612:   PetscFList rlist = 0;
613:   char rname[256] = "exact";

616:   PetscNew(*user,&user);
617:   ctx->physics.sample         = PhysicsSample_Traffic;
618:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
619:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
620:   ctx->physics.user           = user;
621:   ctx->physics.dof            = 1;
622:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
623:   user->a = 0.5;
624:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Traffic_Exact);
625:   RiemannListAdd(&rlist,"roe",    PhysicsRiemann_Traffic_Roe);
626:   RiemannListAdd(&rlist,"lxf",    PhysicsRiemann_Traffic_LxF);
627:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Traffic_Rusanov);
628:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");
629:   {
630:     PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,PETSC_NULL);
631:     PetscOptionsList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
632:   }
633:   PetscOptionsEnd();

635:   RiemannListFind(rlist,rname,&r);
636:   PetscFListDestroy(&rlist);
637:   ctx->physics.riemann = r;

639:   /* *
640:   * Hack to deal with LxF in semi-discrete form
641:   * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
642:   * */
643:   if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3*user->a;

645:   return(0);
646: }




651: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */

653: typedef struct {
654:   PetscReal acoustic_speed;
655: } IsoGasCtx;

657: static inline void IsoGasFlux(PetscReal c,const PetscScalar *u,PetscScalar *f)
658: {
659:   f[0] = u[1];
660:   f[1] = PetscSqr(u[1])/u[0] + c*c*u[0];
661: }

665: static PetscErrorCode PhysicsSample_IsoGas(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
666: {

669:   if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solutions not implemented for t > 0");
670:   switch (initial) {
671:     case 0:
672:       u[0] = (x < 0) ? 1 : 0.5;
673:       u[1] = (x < 0) ? 1 : 0.7;
674:       break;
675:     case 1:
676:       u[0] = 1+0.5*sin(2*M_PI*x);
677:       u[1] = 1*u[0];
678:       break;
679:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
680:   }
681:   return(0);
682: }

686: static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
687: {
688:   IsoGasCtx *phys = (IsoGasCtx*)vctx;
689:   PetscReal c = phys->acoustic_speed;
690:   PetscScalar ubar,du[2],a[2],fL[2],fR[2],lam[2],ustar[2],R[2][2];
691:   PetscInt i;

694:   ubar = (uL[1]/PetscSqrtScalar(uL[0]) + uR[1]/PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0]));
695:   /* write fluxuations in characteristic basis */
696:   du[0] = uR[0] - uL[0];
697:   du[1] = uR[1] - uL[1];
698:   a[0] = (1/(2*c)) * ((ubar + c)*du[0] - du[1]);
699:   a[1] = (1/(2*c)) * ((-ubar + c)*du[0] + du[1]);
700:   /* wave speeds */
701:   lam[0] = ubar - c;
702:   lam[1] = ubar + c;
703:   /* Right eigenvectors */
704:   R[0][0] = 1; R[0][1] = ubar-c;
705:   R[1][0] = 1; R[1][1] = ubar+c;
706:   /* Compute state in star region (between the 1-wave and 2-wave) */
707:   for (i=0; i<2; i++) ustar[i] = uL[i] + a[0]*R[0][i];
708:   if (uL[1]/uL[0] < c && c < ustar[1]/ustar[0]) { /* 1-wave is sonic rarefaction */
709:     PetscScalar ufan[2];
710:     ufan[0] = uL[0]*PetscExpScalar(uL[1]/(uL[0]*c) - 1);
711:     ufan[1] = c*ufan[0];
712:     IsoGasFlux(c,ufan,flux);
713:   } else if (ustar[1]/ustar[0] < -c && -c < uR[1]/uR[0]) { /* 2-wave is sonic rarefaction */
714:     PetscScalar ufan[2];
715:     ufan[0] = uR[0]*PetscExpScalar(-uR[1]/(uR[0]*c) - 1);
716:     ufan[1] = -c*ufan[0];
717:     IsoGasFlux(c,ufan,flux);
718:   } else {                      /* Centered form */
719:     IsoGasFlux(c,uL,fL);
720:     IsoGasFlux(c,uR,fR);
721:     for (i=0; i<2; i++) {
722:       PetscScalar absdu = PetscAbsScalar(lam[0])*a[0]*R[0][i] + PetscAbsScalar(lam[1])*a[1]*R[1][i];
723:       flux[i] = 0.5*(fL[i]+fR[i]) - 0.5*absdu;
724:     }
725:   }
726:   *maxspeed = MaxAbs(lam[0],lam[1]);
727:   return(0);
728: }

732: static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
733: {
734:   IsoGasCtx *phys = (IsoGasCtx*)vctx;
735:   PetscReal c = phys->acoustic_speed;
736:   PetscScalar ustar[2];
737:   struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
738:   PetscInt i;

741:   if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
742:   {
743:     /* Solve for star state */
744:     PetscScalar res,tmp,rho = 0.5*(L.rho + R.rho); /* initial guess */
745:     for (i=0; i<20; i++) {
746:       PetscScalar fr,fl,dfr,dfl;
747:       fl = (L.rho < rho)
748:         ? (rho-L.rho)/PetscSqrtScalar(L.rho*rho)       /* shock */
749:         : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
750:       fr = (R.rho < rho)
751:         ? (rho-R.rho)/PetscSqrtScalar(R.rho*rho)       /* shock */
752:         : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
753:       res = R.u-L.u + c*(fr+fl);
754:       if (!isfinite(res)) SETERRQ1(PETSC_COMM_SELF,1,"non-finite residual=%g",res);
755:       if (PetscAbsScalar(res) < 1e-10) {
756:         star.rho = rho;
757:         star.u   = L.u - c*fl;
758:         goto converged;
759:       }
760:       dfl = (L.rho < rho)
761:         ? 1/PetscSqrtScalar(L.rho*rho)*(1 - 0.5*(rho-L.rho)/rho)
762:         : 1/rho;
763:       dfr = (R.rho < rho)
764:         ? 1/PetscSqrtScalar(R.rho*rho)*(1 - 0.5*(rho-R.rho)/rho)
765:         : 1/rho;
766:       tmp = rho - res/(c*(dfr+dfl));
767:       if (tmp <= 0) rho /= 2;   /* Guard against Newton shooting off to a negative density */
768:       else rho = tmp;
769:       if (!((rho > 0) && isnormal(rho))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate rho=%g",rho);
770:     }
771:     SETERRQ1(PETSC_COMM_SELF,1,"Newton iteration for star.rho diverged after %d iterations",i);
772:   }
773:   converged:
774:   if (L.u-c < 0 && 0 < star.u-c) { /* 1-wave is sonic rarefaction */
775:     PetscScalar ufan[2];
776:     ufan[0] = L.rho*PetscExpScalar(L.u/c - 1);
777:     ufan[1] = c*ufan[0];
778:     IsoGasFlux(c,ufan,flux);
779:   } else if (star.u+c < 0 && 0 < R.u+c) { /* 2-wave is sonic rarefaction */
780:     PetscScalar ufan[2];
781:     ufan[0] = R.rho*PetscExpScalar(-R.u/c - 1);
782:     ufan[1] = -c*ufan[0];
783:     IsoGasFlux(c,ufan,flux);
784:   } else if ((L.rho >= star.rho && L.u-c >= 0)
785:              || (L.rho < star.rho && (star.rho*star.u-L.rho*L.u)/(star.rho-L.rho) > 0)) {
786:     /* 1-wave is supersonic rarefaction, or supersonic shock */
787:     IsoGasFlux(c,uL,flux);
788:   } else if ((star.rho <= R.rho && R.u+c <= 0)
789:              || (star.rho > R.rho && (R.rho*R.u-star.rho*star.u)/(R.rho-star.rho) < 0)) {
790:     /* 2-wave is supersonic rarefaction or supersonic shock */
791:     IsoGasFlux(c,uR,flux);
792:   } else {
793:     ustar[0] = star.rho;
794:     ustar[1] = star.rho*star.u;
795:     IsoGasFlux(c,ustar,flux);
796:   }
797:   *maxspeed = MaxAbs(MaxAbs(star.u-c,star.u+c),MaxAbs(L.u-c,R.u+c));
798:   return(0);
799: }

803: static PetscErrorCode PhysicsRiemann_IsoGas_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
804: {
805:   IsoGasCtx *phys = (IsoGasCtx*)vctx;
806:   PetscScalar c = phys->acoustic_speed,fL[2],fR[2],s;
807:   struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};

810:   if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
811:   IsoGasFlux(c,uL,fL);
812:   IsoGasFlux(c,uR,fR);
813:   s = PetscMax(PetscAbs(L.u),PetscAbs(R.u))+c;
814:   flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
815:   flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
816:   *maxspeed = s;
817:   return(0);
818: }

822: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi)
823: {
824:   IsoGasCtx *phys = (IsoGasCtx*)vctx;
825:   PetscReal c = phys->acoustic_speed;

829:   X[0*2+0] = 1;
830:   X[0*2+1] = u[1]/u[0] - c;
831:   X[1*2+0] = 1;
832:   X[1*2+1] = u[1]/u[0] + c;
833:   PetscMemcpy(Xi,X,4*sizeof(X[0]));
834:   Kernel_A_gets_inverse_A_2(Xi,0);
835:   return(0);
836: }

840: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
841: {
843:   IsoGasCtx *user;
844:   PetscFList rlist = 0,rclist = 0;
845:   char rname[256] = "exact",rcname[256] = "characteristic";

848:   PetscNew(*user,&user);
849:   ctx->physics.sample         = PhysicsSample_IsoGas;
850:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
851:   ctx->physics.user           = user;
852:   ctx->physics.dof            = 2;
853:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
854:   PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
855:   user->acoustic_speed = 1;
856:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_IsoGas_Exact);
857:   RiemannListAdd(&rlist,"roe",    PhysicsRiemann_IsoGas_Roe);
858:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_IsoGas_Rusanov);
859:   ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_IsoGas);
860:   ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
861:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");
862:   {
863:     PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,PETSC_NULL);
864:     PetscOptionsList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
865:     PetscOptionsList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
866:   }
867:   PetscOptionsEnd();
868:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
869:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
870:   PetscFListDestroy(&rlist);
871:   PetscFListDestroy(&rclist);
872:   return(0);
873: }



877: /* --------------------------------- Shallow Water ----------------------------------- */

879: typedef struct {
880:   PetscReal gravity;
881: } ShallowCtx;

883: static inline void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
884: {
885:   f[0] = u[1];
886:   f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
887: }

891: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
892: {
893:   ShallowCtx *phys = (ShallowCtx*)vctx;
894:   PetscScalar g = phys->gravity,ustar[2],cL,cR,c,cstar;
895:   struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
896:   PetscInt i;

899:   if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed thickness is negative");
900:   cL = PetscSqrtScalar(g*L.h);
901:   cR = PetscSqrtScalar(g*R.h);
902:   c = PetscMax(cL,cR);
903:   {
904:     /* Solve for star state */
905:     const PetscInt maxits = 50;
906:     PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */
907:     h0 = h;
908:     for (i=0; i<maxits; i++) {
909:       PetscScalar fr,fl,dfr,dfl;
910:       fl = (L.h < h)
911:         ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */
912:         : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h);   /* rarefaction */
913:       fr = (R.h < h)
914:         ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */
915:         : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h);   /* rarefaction */
916:       res = R.u - L.u + fr + fl;
917:       if (!isfinite(res)) SETERRQ1(PETSC_COMM_SELF,1,"non-finite residual=%g",res);
918:       //PetscPrintf(PETSC_COMM_WORLD,"h=%g, res[%d] = %g\n",h,i,res);
919:       if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h-h0) < 1e-8)) {
920:         star.h = h;
921:         star.u = L.u - fl;
922:         goto converged;
923:       } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) {        /* Line search */
924:         h = 0.8*h0 + 0.2*h;
925:         continue;
926:       }
927:       /* Accept the last step and take another */
928:       res0 = res;
929:       h0 = h;
930:       dfl = (L.h < h)
931:         ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h)
932:         : PetscSqrtScalar(g/h);
933:       dfr = (R.h < h)
934:         ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h)
935:         : PetscSqrtScalar(g/h);
936:       tmp = h - res/(dfr+dfl);
937:       if (tmp <= 0) h /= 2;   /* Guard against Newton shooting off to a negative thickness */
938:       else h = tmp;
939:       if (!((h > 0) && isnormal(h))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate h=%g",h);
940:     }
941:     SETERRQ1(PETSC_COMM_SELF,1,"Newton iteration for star.h diverged after %d iterations",i);
942:   }
943:   converged:
944:   cstar = PetscSqrtScalar(g*star.h);
945:   if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */
946:     PetscScalar ufan[2];
947:     ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL);
948:     ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0];
949:     ShallowFlux(phys,ufan,flux);
950:   } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */
951:     PetscScalar ufan[2];
952:     ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR);
953:     ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0];
954:     ShallowFlux(phys,ufan,flux);
955:   } else if ((L.h >= star.h && L.u-c >= 0)
956:              || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) {
957:     /* 1-wave is right-travelling shock (supersonic) */
958:     ShallowFlux(phys,uL,flux);
959:   } else if ((star.h <= R.h && R.u+c <= 0)
960:              || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) {
961:     /* 2-wave is left-travelling shock (supersonic) */
962:     ShallowFlux(phys,uR,flux);
963:   } else {
964:     ustar[0] = star.h;
965:     ustar[1] = star.h*star.u;
966:     ShallowFlux(phys,ustar,flux);
967:   }
968:   *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR));
969:   return(0);
970: }

974: static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
975: {
976:   ShallowCtx *phys = (ShallowCtx*)vctx;
977:   PetscScalar g = phys->gravity,fL[2],fR[2],s;
978:   struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};

981:   if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed thickness is negative");
982:   ShallowFlux(phys,uL,fL);
983:   ShallowFlux(phys,uR,fR);
984:   s = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h));
985:   flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
986:   flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
987:   *maxspeed = s;
988:   return(0);
989: }

993: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi)
994: {
995:   ShallowCtx *phys = (ShallowCtx*)vctx;
996:   PetscReal c;

1000:   c = PetscSqrtScalar(u[0]*phys->gravity);
1001:   X[0*2+0] = 1;
1002:   X[0*2+1] = u[1]/u[0] - c;
1003:   X[1*2+0] = 1;
1004:   X[1*2+1] = u[1]/u[0] + c;
1005:   PetscMemcpy(Xi,X,4*sizeof(X[0]));
1006:   Kernel_A_gets_inverse_A_2(Xi,0);
1007:   return(0);
1008: }

1012: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
1013: {
1015:   ShallowCtx *user;
1016:   PetscFList rlist = 0,rclist = 0;
1017:   char rname[256] = "exact",rcname[256] = "characteristic";

1020:   PetscNew(*user,&user);
1021:   /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
1022:   ctx->physics.sample         = PhysicsSample_IsoGas;
1023:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
1024:   ctx->physics.user           = user;
1025:   ctx->physics.dof            = 2;
1026:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1027:   PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1028:   user->gravity = 1;
1029:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Shallow_Exact);
1030:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov);
1031:   ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Shallow);
1032:   ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1033:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");
1034:   {
1035:     PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,PETSC_NULL);
1036:     PetscOptionsList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
1037:     PetscOptionsList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
1038:   }
1039:   PetscOptionsEnd();
1040:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
1041:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1042:   PetscFListDestroy(&rlist);
1043:   PetscFListDestroy(&rclist);
1044:   return(0);
1045: }



1049: /* --------------------------------- Finite Volume Solver ----------------------------------- */

1053: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
1054: {
1055:   FVCtx          *ctx = (FVCtx*)vctx;
1056:   PetscErrorCode  ierr;
1057:   PetscInt        i,j,k,Mx,dof,xs,xm;
1058:   PetscReal       hx,cfl_idt = 0;
1059:   PetscScalar    *x,*f,*slope;
1060:   Vec             Xloc;

1063:   DMGetLocalVector(ctx->da,&Xloc);
1064:   DMDAGetInfo(ctx->da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1065:   hx = (ctx->xmax - ctx->xmin)/Mx;
1066:   DMGlobalToLocalBegin(ctx->da,X,INSERT_VALUES,Xloc);
1067:   DMGlobalToLocalEnd  (ctx->da,X,INSERT_VALUES,Xloc);

1069:   VecZeroEntries(F);

1071:   DMDAVecGetArray(ctx->da,Xloc,&x);
1072:   DMDAVecGetArray(ctx->da,F,&f);
1073:   DMDAGetArray(ctx->da,PETSC_TRUE,&slope);

1075:   DMDAGetCorners(ctx->da,&xs,0,0,&xm,0,0);

1077:   if (ctx->bctype == FVBC_OUTFLOW) {
1078:     for (i=xs-2; i<0; i++) {
1079:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
1080:     }
1081:     for (i=Mx; i<xs+xm+2; i++) {
1082:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
1083:     }
1084:   }
1085:   for (i=xs-1; i<xs+xm+1; i++) {
1086:     struct _LimitInfo info;
1087:     PetscScalar *cjmpL,*cjmpR;
1088:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
1089:     (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv);
1090:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
1091:     PetscMemzero(ctx->cjmpLR,2*dof*sizeof(ctx->cjmpLR[0]));
1092:     cjmpL = &ctx->cjmpLR[0];
1093:     cjmpR = &ctx->cjmpLR[dof];
1094:     for (j=0; j<dof; j++) {
1095:       PetscScalar jmpL,jmpR;
1096:       jmpL = x[(i+0)*dof+j] - x[(i-1)*dof+j];
1097:       jmpR = x[(i+1)*dof+j] - x[(i+0)*dof+j];
1098:       for (k=0; k<dof; k++) {
1099:         cjmpL[k] += ctx->Rinv[k+j*dof] * jmpL;
1100:         cjmpR[k] += ctx->Rinv[k+j*dof] * jmpR;
1101:       }
1102:     }
1103:     /* Apply limiter to the left and right characteristic jumps */
1104:     info.m = dof;
1105:     info.hx = hx;
1106:     (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
1107:     for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
1108:     for (j=0; j<dof; j++) {
1109:       PetscScalar tmp = 0;
1110:       for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof] * ctx->cslope[k];
1111:       slope[i*dof+j] = tmp;
1112:     }
1113:   }

1115:   for (i=xs; i<xs+xm+1; i++) {
1116:     PetscReal maxspeed;
1117:     PetscScalar *uL,*uR;
1118:     uL = &ctx->uLR[0];
1119:     uR = &ctx->uLR[dof];
1120:     for (j=0; j<dof; j++) {
1121:       uL[j] = x[(i-1)*dof+j] + slope[(i-1)*dof+j]*hx/2;
1122:       uR[j] = x[(i-0)*dof+j] - slope[(i-0)*dof+j]*hx/2;
1123:     }
1124:     (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed);
1125:     cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */

1127:     if (i > xs) {
1128:       for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
1129:     }
1130:     if (i < xs+xm) {
1131:       for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
1132:     }
1133:   }

1135:   DMDAVecRestoreArray(ctx->da,Xloc,&x);
1136:   DMDAVecRestoreArray(ctx->da,F,&f);
1137:   DMDARestoreArray(ctx->da,PETSC_TRUE,&slope);
1138:   DMRestoreLocalVector(ctx->da,&Xloc);

1140:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,((PetscObject)ctx->da)->comm);
1141:   if (0) {
1142:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
1143:     PetscReal dt,tnow;
1144:     TSGetTimeStep(ts,&dt);
1145:     TSGetTime(ts,&tnow);
1146:     if (dt > 0.5/ctx->cfl_idt) {
1147:       if (1) {
1148:         PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",tnow,dt,0.5/ctx->cfl_idt);
1149:       } else {
1150:         SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",dt,ctx->cfl/ctx->cfl_idt);
1151:       }
1152:     }
1153:   }
1154:   return(0);
1155: }

1159: static PetscErrorCode FVSample(FVCtx *ctx,PetscReal time,Vec U)
1160: {
1162:   PetscScalar *u,*uj;
1163:   PetscInt i,j,k,dof,xs,xm,Mx;

1166:   if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,1,"Physics has not provided a sampling function");
1167:   DMDAGetInfo(ctx->da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1168:   DMDAGetCorners(ctx->da,&xs,0,0,&xm,0,0);
1169:   DMDAVecGetArray(ctx->da,U,&u);
1170:   PetscMalloc(dof*sizeof uj[0],&uj);
1171:   for (i=xs; i<xs+xm; i++) {
1172:     const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h;
1173:     const PetscInt N = 200;
1174:     /* Integrate over cell i using trapezoid rule with N points. */
1175:     for (k=0; k<dof; k++) u[i*dof+k] = 0;
1176:     for (j=0; j<N+1; j++) {
1177:       PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N;
1178:       (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
1179:       for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N)?0.5:1.0)*uj[k]/N;
1180:     }
1181:   }
1182:   DMDAVecRestoreArray(ctx->da,U,&u);
1183:   PetscFree(uj);
1184:   return(0);
1185: }

1189: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
1190: {
1192:   PetscReal xmin,xmax;
1193:   PetscScalar sum,*x,tvsum,tvgsum;
1194:   PetscInt imin,imax,Mx,i,j,xs,xm,dof;
1195:   Vec Xloc;
1196:   PetscBool  iascii;

1199:   PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1200:   if (iascii) {
1201:     /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
1202:     DMGetLocalVector(da,&Xloc);
1203:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1204:     DMGlobalToLocalEnd  (da,X,INSERT_VALUES,Xloc);
1205:     DMDAVecGetArray(da,Xloc,&x);
1206:     DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1207:     DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1208:     tvsum = 0;
1209:     for (i=xs; i<xs+xm; i++) {
1210:       for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j] - x[(i-1)*dof+j]);
1211:     }
1212:     MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_MAX,((PetscObject)da)->comm);
1213:     DMDAVecRestoreArray(da,Xloc,&x);
1214:     DMRestoreLocalVector(da,&Xloc);

1216:     VecMin(X,&imin,&xmin);
1217:     VecMax(X,&imax,&xmax);
1218:     VecSum(X,&sum);
1219:     PetscViewerASCIIPrintf(viewer,"Solution range [%8.5f,%8.5f] with extrema at %d and %d, mean %8.5f, ||x||_TV %8.5f\n",xmin,xmax,imin,imax,sum/Mx,tvgsum/Mx);
1220:   } else {
1221:     SETERRQ(PETSC_COMM_SELF,1,"Viewer type not supported");
1222:   }
1223:   return(0);
1224: }

1228: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,PetscReal t,Vec X,PetscReal *nrm1,PetscReal *nrmsup)
1229: {
1231:   Vec Y;
1232:   PetscInt Mx;

1235:   VecGetSize(X,&Mx);
1236:   VecDuplicate(X,&Y);
1237:   FVSample(ctx,t,Y);
1238:   VecAYPX(Y,-1,X);
1239:   VecNorm(Y,NORM_1,nrm1);
1240:   VecNorm(Y,NORM_INFINITY,nrmsup);
1241:   *nrm1 /= Mx;
1242:   VecDestroy(&Y);
1243:   return(0);
1244: }


1249: int main(int argc,char *argv[])
1250: {
1251:   char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
1252:   PetscFList limiters = 0,physics = 0;
1253:   MPI_Comm comm;
1254:   TS ts;
1255:   Vec X,X0,R;
1256:   FVCtx ctx;
1257:   PetscInt i,dof,xs,xm,Mx,draw = 0;
1258:   PetscBool  view_final = PETSC_FALSE;
1259:   PetscReal ptime;

1262:   PetscInitialize(&argc,&argv,0,help);
1263:   comm = PETSC_COMM_WORLD;
1264:   PetscMemzero(&ctx,sizeof(ctx));

1266:   /* Register limiters to be available on the command line */
1267:   PetscFListAdd(&limiters,"upwind"          ,"",(void(*)(void))Limit_Upwind);
1268:   PetscFListAdd(&limiters,"lax-wendroff"    ,"",(void(*)(void))Limit_LaxWendroff);
1269:   PetscFListAdd(&limiters,"beam-warming"    ,"",(void(*)(void))Limit_BeamWarming);
1270:   PetscFListAdd(&limiters,"fromm"           ,"",(void(*)(void))Limit_Fromm);
1271:   PetscFListAdd(&limiters,"minmod"          ,"",(void(*)(void))Limit_Minmod);
1272:   PetscFListAdd(&limiters,"superbee"        ,"",(void(*)(void))Limit_Superbee);
1273:   PetscFListAdd(&limiters,"mc"              ,"",(void(*)(void))Limit_MC);
1274:   PetscFListAdd(&limiters,"vanleer"         ,"",(void(*)(void))Limit_VanLeer);
1275:   PetscFListAdd(&limiters,"vanalbada"       ,"",(void(*)(void))Limit_VanAlbada);
1276:   PetscFListAdd(&limiters,"vanalbadatvd"    ,"",(void(*)(void))Limit_VanAlbadaTVD);
1277:   PetscFListAdd(&limiters,"koren"           ,"",(void(*)(void))Limit_Koren);
1278:   PetscFListAdd(&limiters,"korensym"        ,"",(void(*)(void))Limit_KorenSym);
1279:   PetscFListAdd(&limiters,"koren3"          ,"",(void(*)(void))Limit_Koren3);
1280:   PetscFListAdd(&limiters,"cada-torrilhon2" ,"",(void(*)(void))Limit_CadaTorrilhon2);
1281:   PetscFListAdd(&limiters,"cada-torrilhon3-r0p1","",(void(*)(void))Limit_CadaTorrilhon3R0p1);
1282:   PetscFListAdd(&limiters,"cada-torrilhon3-r1"  ,"",(void(*)(void))Limit_CadaTorrilhon3R1);
1283:   PetscFListAdd(&limiters,"cada-torrilhon3-r10" ,"",(void(*)(void))Limit_CadaTorrilhon3R10);
1284:   PetscFListAdd(&limiters,"cada-torrilhon3-r100","",(void(*)(void))Limit_CadaTorrilhon3R100);

1286:   /* Register physical models to be available on the command line */
1287:   PetscFListAdd(&physics,"advect"          ,"",(void(*)(void))PhysicsCreate_Advect);
1288:   PetscFListAdd(&physics,"burgers"         ,"",(void(*)(void))PhysicsCreate_Burgers);
1289:   PetscFListAdd(&physics,"traffic"         ,"",(void(*)(void))PhysicsCreate_Traffic);
1290:   PetscFListAdd(&physics,"isogas"          ,"",(void(*)(void))PhysicsCreate_IsoGas);
1291:   PetscFListAdd(&physics,"shallow"         ,"",(void(*)(void))PhysicsCreate_Shallow);

1293:   ctx.cfl = 0.9; ctx.bctype = FVBC_PERIODIC;
1294:   ctx.xmin = 0; ctx.xmax = 1;
1295:   PetscOptionsBegin(comm,PETSC_NULL,"Finite Volume solver options","");
1296:   {
1297:     PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,PETSC_NULL);
1298:     PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,PETSC_NULL);
1299:     PetscOptionsList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),PETSC_NULL);
1300:     PetscOptionsList("-physics","Name of physics (Riemann solver and characteristics) to use","",physics,physname,physname,sizeof(physname),PETSC_NULL);
1301:     PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,PETSC_NULL);
1302:     PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof final_fname,&view_final);
1303:     PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,PETSC_NULL);
1304:     PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,PETSC_NULL);
1305:     PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,PETSC_NULL);
1306:     PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,PETSC_NULL);
1307:   }
1308:   PetscOptionsEnd();

1310:   /* Choose the limiter from the list of registered limiters */
1311:   PetscFListFind(limiters,comm,lname,PETSC_FALSE,(void(**)(void))&ctx.limit);
1312:   if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);

1314:   /* Choose the physics from the list of registered models */
1315:   {
1316:     PetscErrorCode (*r)(FVCtx*);
1317:     PetscFListFind(physics,comm,physname,PETSC_FALSE,(void(**)(void))&r);
1318:     if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
1319:     /* Create the physics, will set the number of fields and their names */
1320:     (*r)(&ctx);
1321:   }

1323:   /* Create a DMDA to manage the parallel grid */
1324:   DMDACreate1d(comm,DMDA_BOUNDARY_PERIODIC,-50,ctx.physics.dof,2,PETSC_NULL,&ctx.da);
1325:   /* Inform the DMDA of the field names provided by the physics. */
1326:   /* The names will be shown in the title bars when run with -ts_monitor_solution */
1327:   for (i=0; i<ctx.physics.dof; i++) {
1328:     DMDASetFieldName(ctx.da,i,ctx.physics.fieldname[i]);
1329:   }
1330:   DMDAGetInfo(ctx.da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1331:   DMDAGetCorners(ctx.da,&xs,0,0,&xm,0,0);

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

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

1340:   /* Create a vector to store the solution and to save the initial state */
1341:   DMCreateGlobalVector(ctx.da,&X);
1342:   VecDuplicate(X,&X0);
1343:   VecDuplicate(X,&R);

1345:   /* Create a time-stepping object */
1346:   TSCreate(comm,&ts);
1347:   TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
1348:   TSSetType(ts,TSSSP);
1349:   TSSetDuration(ts,1000,10);

1351:   /* Compute initial conditions and starting time step */
1352:   FVSample(&ctx,0,X0);
1353:   FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
1354:   VecCopy(X0,X);                        /* The function value was not used so we set X=X0 again */
1355:   TSSetInitialTimeStep(ts,0,ctx.cfl/ctx.cfl_idt);

1357:   TSSetFromOptions(ts); /* Take runtime options */

1359:   SolutionStatsView(ctx.da,X,PETSC_VIEWER_STDOUT_WORLD);

1361:   {
1362:     PetscReal nrm1,nrmsup;
1363:     PetscInt steps;

1365:     TSSolve(ts,X,&ptime);
1366:     TSGetTimeStepNumber(ts,&steps);

1368:     PetscPrintf(comm,"Final time %8.5f, steps %d\n",ptime,steps);
1369:     if (ctx.exact) {
1370:       SolutionErrorNorms(&ctx,ptime,X,&nrm1,&nrmsup);
1371:       PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e  ||x-x_e||_sup %8.4e\n",nrm1,nrmsup);
1372:     }
1373:   }

1375:   SolutionStatsView(ctx.da,X,PETSC_VIEWER_STDOUT_WORLD);

1377:   if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
1378:   if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
1379:   if (draw & 0x4) {
1380:     Vec Y;
1381:     VecDuplicate(X,&Y);
1382:     FVSample(&ctx,ptime,Y);
1383:     VecAYPX(Y,-1,X);
1384:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1385:     VecDestroy(&Y);
1386:   }

1388:   if (view_final) {
1389:     PetscViewer viewer;
1390:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1391:     PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1392:     VecView(X,viewer);
1393:     PetscViewerDestroy(&viewer);
1394:   }

1396:   /* Clean up */
1397:   (*ctx.physics.destroy)(ctx.physics.user);
1398:   for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
1399:   PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
1400:   PetscFree2(ctx.uLR,ctx.flux);
1401:   VecDestroy(&X);
1402:   VecDestroy(&X0);
1403:   VecDestroy(&R);
1404:   DMDestroy(&ctx.da);
1405:   TSDestroy(&ts);
1406:   PetscFinalize();
1407:   return 0;
1408: }