Actual source code: ex9.c

petsc-dev 2014-07-21
Report Typos and Errors
  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:   "  acoustics   - Acoustic wave propagation\n"
 10:   "                u_t       + (c*z*v)_x             = 0\n"
 11:   "                v_t       + (c/z*u)_x             = 0\n"
 12:   "  isogas      - Isothermal gas dynamics\n"
 13:   "                rho_t     + (rho*u)_x             = 0\n"
 14:   "                (rho*u)_t + (rho*u^2 + c^2*rho)_x = 0\n"
 15:   "  shallow     - Shallow water equations\n"
 16:   "                h_t       + (h*u)_x               = 0\n"
 17:   "                (h*u)_t   + (h*u^2 + g*h^2/2)_x   = 0\n"
 18:   "Some of these physical models have multiple Riemann solvers, select these with -physics_xxx_riemann\n"
 19:   "  exact       - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
 20:   "                the states across shocks and rarefactions\n"
 21:   "  roe         - Linearized scheme, usually with an entropy fix inside sonic rarefactions\n"
 22:   "The systems provide a choice of reconstructions with -physics_xxx_reconstruct\n"
 23:   "  characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
 24:   "  conservative   - Limit the conservative variables directly, can cause undesired interaction of waves\n\n"
 25:   "A variety of limiters for high-resolution TVD limiters are available with -limit\n"
 26:   "  upwind,minmod,superbee,mc,vanleer,vanalbada,koren,cada-torillhon (last two are nominally third order)\n"
 27:   "  and non-TVD schemes lax-wendroff,beam-warming,fromm\n\n"
 28:   "To preserve the TVD property, one should time step with a strong stability preserving method.\n"
 29:   "The optimal high order explicit Runge-Kutta methods in TSSSP are recommended for non-stiff problems.\n\n"
 30:   "Several initial conditions can be chosen with -initial N\n\n"
 31:   "The problem size should be set with -da_grid_x M\n\n";

 33: #include <petscts.h>
 34: #include <petscdm.h>
 35: #include <petscdmda.h>
 36: #include <petscdraw.h>

 38: #include <petsc-private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */

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

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


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

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


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

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

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

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

184:   MPI_Comm comm;
185:   char     prefix[256];

187:   /* Local work arrays */
188:   PetscScalar *R,*Rinv;         /* Characteristic basis, and it's inverse.  COLUMN-MAJOR */
189:   PetscScalar *cjmpLR;          /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
190:   PetscScalar *cslope;          /* Limited slope, written in characteristic basis */
191:   PetscScalar *uLR;             /* Solution at left and right of interface, conservative variables, len=2*dof */
192:   PetscScalar *flux;            /* Flux across interface */
193:   PetscReal   *speeds;          /* Speeds of each wave */

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(PetscFunctionList *flist,const char *name,RiemannFunction rsolve)
209: {

213:   PetscFunctionListAdd(flist,name,rsolve);
214:   return(0);
215: }

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

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

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

236:   PetscFunctionListAdd(flist,name,r);
237:   return(0);
238: }

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

247:   PetscFunctionListFind(flist,name,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,PetscReal *speeds)
264: {
265:   PetscInt i,j;

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

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

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



288: /* --------------------------------- Advection ----------------------------------- */

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

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

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

310: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
311: {
312:   AdvectCtx *ctx = (AdvectCtx*)vctx;

315:   X[0]      = 1.;
316:   Xi[0]     = 1.;
317:   speeds[0] = ctx->a;
318:   return(0);
319: }

323: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
324: {
325:   AdvectCtx *ctx = (AdvectCtx*)vctx;
326:   PetscReal a    = ctx->a,x0;

329:   switch (bctype) {
330:     case FVBC_OUTFLOW:   x0 = x-a*t; break;
331:     case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
332:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
333:   }
334:   switch (initial) {
335:     case 0: u[0] = (x0 < 0) ? 1 : -1; break;
336:     case 1: u[0] = (x0 < 0) ? -1 : 1; break;
337:     case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
338:     case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
339:     case 4: u[0] = PetscAbs(x0); break;
340:     case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
341:     case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
342:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
343:   }
344:   return(0);
345: }

349: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
350: {
352:   AdvectCtx      *user;

355:   PetscNew(&user);
356:   ctx->physics.sample         = PhysicsSample_Advect;
357:   ctx->physics.riemann        = PhysicsRiemann_Advect;
358:   ctx->physics.characteristic = PhysicsCharacteristic_Advect;
359:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
360:   ctx->physics.user           = user;
361:   ctx->physics.dof            = 1;
362:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
363:   user->a = 1;
364:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
365:   {
366:     PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
367:   }
368:   PetscOptionsEnd();
369:   return(0);
370: }



374: /* --------------------------------- Burgers ----------------------------------- */

376: typedef struct {
377:   PetscReal lxf_speed;
378: } BurgersCtx;

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

386:   if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
387:   switch (initial) {
388:     case 0: u[0] = (x < 0) ? 1 : -1; break;
389:     case 1:
390:       if       (x < -t) u[0] = -1;
391:       else if  (x < t)  u[0] = x/t;
392:       else              u[0] = 1;
393:       break;
394:     case 2:
395:       if      (x < 0)       u[0] = 0;
396:       else if (x <= t)      u[0] = x/t;
397:       else if (x < 1+0.5*t) u[0] = 1;
398:       else                  u[0] = 0;
399:       break;
400:     case 3:
401:       if       (x < 0.2*t) u[0] = 0.2;
402:       else if  (x < t) u[0] = x/t;
403:       else             u[0] = 1;
404:       break;
405:     case 4:
406:       if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
407:       u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin)));
408:       break;
409:     case 5:                     /* Pure shock solution */
410:       if (x < 0.5*t) u[0] = 1;
411:       else u[0] = 0;
412:       break;
413:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
414:   }
415:   return(0);
416: }

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

424:   if (uL[0] < uR[0]) {          /* rarefaction */
425:     flux[0] = (uL[0]*uR[0] < 0)
426:       ? 0                       /* sonic rarefaction */
427:       : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0]));
428:   } else {                      /* shock */
429:     flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0]));
430:   }
431:   *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
432:   return(0);
433: }

437: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
438: {
439:   PetscReal speed;

442:   speed   = 0.5*(uL[0] + uR[0]);
443:   flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
444:   if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
445:   *maxspeed = speed;
446:   return(0);
447: }

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

457:   c         = ((BurgersCtx*)vctx)->lxf_speed;
458:   fL        = 0.5*PetscSqr(uL[0]);
459:   fR        = 0.5*PetscSqr(uR[0]);
460:   flux[0]   = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
461:   *maxspeed = c;
462:   return(0);
463: }

467: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
468: {
469:   PetscReal   c;
470:   PetscScalar fL,fR;

473:   c         = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0]));
474:   fL        = 0.5*PetscSqr(uL[0]);
475:   fR        = 0.5*PetscSqr(uR[0]);
476:   flux[0]   = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
477:   *maxspeed = c;
478:   return(0);
479: }

483: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
484: {
485:   BurgersCtx        *user;
486:   PetscErrorCode    ierr;
487:   RiemannFunction   r;
488:   PetscFunctionList rlist      = 0;
489:   char              rname[256] = "exact";

492:   PetscNew(&user);

494:   ctx->physics.sample         = PhysicsSample_Burgers;
495:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
496:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
497:   ctx->physics.user           = user;
498:   ctx->physics.dof            = 1;

500:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
501:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Burgers_Exact);
502:   RiemannListAdd(&rlist,"roe",    PhysicsRiemann_Burgers_Roe);
503:   RiemannListAdd(&rlist,"lxf",    PhysicsRiemann_Burgers_LxF);
504:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Burgers_Rusanov);
505:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
506:   {
507:     PetscOptionsFList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
508:   }
509:   PetscOptionsEnd();
510:   RiemannListFind(rlist,rname,&r);
511:   PetscFunctionListDestroy(&rlist);
512:   ctx->physics.riemann = r;

514:   /* *
515:   * Hack to deal with LxF in semi-discrete form
516:   * max speed is 1 for the basic initial conditions (where |u| <= 1)
517:   * */
518:   if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
519:   return(0);
520: }



524: /* --------------------------------- Traffic ----------------------------------- */

526: typedef struct {
527:   PetscReal lxf_speed;
528:   PetscReal a;
529: } TrafficCtx;

531: PETSC_STATIC_INLINE PetscScalar TrafficFlux(PetscScalar a,PetscScalar u) { return a*u*(1-u); }

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

540:   if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
541:   switch (initial) {
542:     case 0:
543:       u[0] = (-a*t < x) ? 2 : 0; break;
544:     case 1:
545:       if      (x < PetscMin(2*a*t,0.5+a*t)) u[0] = -1;
546:       else if (x < 1)                       u[0] = 0;
547:       else                                  u[0] = 1;
548:       break;
549:     case 2:
550:       if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
551:       u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin)));
552:       break;
553:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
554:   }
555:   return(0);
556: }

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

565:   if (uL[0] < uR[0]) {
566:     flux[0] = PetscMin(TrafficFlux(a,uL[0]),
567:                        TrafficFlux(a,uR[0]));
568:   } else {
569:     flux[0] = (uR[0] < 0.5 && 0.5 < uL[0])
570:       ? TrafficFlux(a,0.5)
571:       : PetscMax(TrafficFlux(a,uL[0]),
572:                  TrafficFlux(a,uR[0]));
573:   }
574:   *maxspeed = a*MaxAbs(1-2*uL[0],1-2*uR[0]);
575:   return(0);
576: }

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

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

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

601:   speed     = a*(1 - (uL[0] + uR[0]));
602:   flux[0]   = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]);
603:   *maxspeed = speed;
604:   return(0);
605: }

609: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
610: {
611:   PetscReal a = ((TrafficCtx*)vctx)->a;
612:   PetscReal speed;

615:   speed     = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0]));
616:   flux[0]   = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]);
617:   *maxspeed = speed;
618:   return(0);
619: }

623: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
624: {
625:   PetscErrorCode    ierr;
626:   TrafficCtx        *user;
627:   RiemannFunction   r;
628:   PetscFunctionList rlist      = 0;
629:   char              rname[256] = "exact";

632:   PetscNew(&user);
633:   ctx->physics.sample         = PhysicsSample_Traffic;
634:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
635:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
636:   ctx->physics.user           = user;
637:   ctx->physics.dof            = 1;

639:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
640:   user->a = 0.5;
641:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Traffic_Exact);
642:   RiemannListAdd(&rlist,"roe",    PhysicsRiemann_Traffic_Roe);
643:   RiemannListAdd(&rlist,"lxf",    PhysicsRiemann_Traffic_LxF);
644:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Traffic_Rusanov);
645:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");
646:   {
647:     PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,NULL);
648:     PetscOptionsFList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
649:   }
650:   PetscOptionsEnd();

652:   RiemannListFind(rlist,rname,&r);
653:   PetscFunctionListDestroy(&rlist);

655:   ctx->physics.riemann = r;

657:   /* *
658:   * Hack to deal with LxF in semi-discrete form
659:   * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
660:   * */
661:   if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3*user->a;
662:   return(0);
663: }

665: /* --------------------------------- Linear Acoustics ----------------------------------- */

667: /* Flux: u_t + (A u)_x
668:  * z = sqrt(rho*bulk), c = sqrt(rho/bulk)
669:  * Spectral decomposition: A = R * D * Rinv
670:  * [    cz] = [-z   z] [-c    ] [-1/2z  1/2]
671:  * [c/z   ] = [ 1   1] [     c] [ 1/2z  1/2]
672:  *
673:  * We decompose this into the left-traveling waves Al = R * D^- Rinv
674:  * and the right-traveling waves Ar = R * D^+ * Rinv
675:  * Multiplying out these expressions produces the following two matrices
676:  */

678: typedef struct {
679:   PetscReal c;                  /* speed of sound: c = sqrt(bulk/rho) */
680:   PetscReal z;                  /* impedence: z = sqrt(rho*bulk) */
681: } AcousticsCtx;

683: PETSC_UNUSED PETSC_STATIC_INLINE void AcousticsFlux(AcousticsCtx *ctx,const PetscScalar *u,PetscScalar *f)
684: {
685:   f[0] = ctx->c*ctx->z*u[1];
686:   f[1] = ctx->c/ctx->z*u[0];
687: }

691: static PetscErrorCode PhysicsCharacteristic_Acoustics(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
692: {
693:   AcousticsCtx *phys = (AcousticsCtx*)vctx;
694:   PetscReal    z     = phys->z,c = phys->c;

697:   X[0*2+0]  = -z;
698:   X[0*2+1]  = z;
699:   X[1*2+0]  = 1;
700:   X[1*2+1]  = 1;
701:   Xi[0*2+0] = -1./(2*z);
702:   Xi[0*2+1] = 1./2;
703:   Xi[1*2+0] = 1./(2*z);
704:   Xi[1*2+1] = 1./2;
705:   speeds[0] = -c;
706:   speeds[1] = c;
707:   return(0);
708: }

712: static PetscErrorCode PhysicsSample_Acoustics_Initial(AcousticsCtx *phys,PetscInt initial,PetscReal xmin,PetscReal xmax,PetscReal x,PetscReal *u)
713: {
715:   switch (initial) {
716:   case 0:
717:     u[0] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.2) < 0.1) ? 1 : 0.5;
718:     u[1] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.7) < 0.1) ? 1 : -0.5;
719:     break;
720:   case 1:
721:     u[0] = PetscCosReal(3 * 2*PETSC_PI*x/(xmax-xmin));
722:     u[1] = PetscExpReal(-PetscSqr(x - (xmax + xmin)/2) / (2*PetscSqr(0.2*(xmax - xmin)))) - 0.5;
723:     break;
724:   default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
725:   }
726:   return(0);
727: }

731: static PetscErrorCode PhysicsSample_Acoustics(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
732: {
733:   AcousticsCtx   *phys = (AcousticsCtx*)vctx;
734:   PetscReal      c     = phys->c;
735:   PetscReal      x0a,x0b,u0a[2],u0b[2],tmp[2];
736:   PetscReal      X[2][2],Xi[2][2],dummy[2];

740:   switch (bctype) {
741:   case FVBC_OUTFLOW:
742:     x0a = x+c*t;
743:     x0b = x-c*t;
744:     break;
745:   case FVBC_PERIODIC:
746:     x0a = RangeMod(x+c*t,xmin,xmax);
747:     x0b = RangeMod(x-c*t,xmin,xmax);
748:     break;
749:   default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
750:   }
751:   PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0a,u0a);
752:   PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0b,u0b);
753:   PhysicsCharacteristic_Acoustics(vctx,2,u,&X[0][0],&Xi[0][0],dummy);
754:   tmp[0] = Xi[0][0]*u0a[0] + Xi[0][1]*u0a[1];
755:   tmp[1] = Xi[1][0]*u0b[0] + Xi[1][1]*u0b[1];
756:   u[0]   = X[0][0]*tmp[0] + X[0][1]*tmp[1];
757:   u[1]   = X[1][0]*tmp[0] + X[1][1]*tmp[1];
758:   return(0);
759: }

763: static PetscErrorCode PhysicsRiemann_Acoustics_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
764: {
765:   AcousticsCtx *phys = (AcousticsCtx*)vctx;
766:   PetscReal    c     = phys->c,z = phys->z;
767:   PetscReal
768:     Al[2][2] = {{-c/2     , c*z/2  },
769:                 {c/(2*z)  , -c/2   }}, /* Left traveling waves */
770:     Ar[2][2] = {{c/2      , c*z/2  },
771:                 {c/(2*z)  , c/2    }}; /* Right traveling waves */

774:   flux[0]   = Al[0][0]*uR[0] + Al[0][1]*uR[1] + Ar[0][0]*uL[0] + Ar[0][1]*uL[1];
775:   flux[1]   = Al[1][0]*uR[0] + Al[1][1]*uR[1] + Ar[1][0]*uL[0] + Ar[1][1]*uL[1];
776:   *maxspeed = c;
777:   return(0);
778: }

782: static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx)
783: {
784:   PetscErrorCode    ierr;
785:   AcousticsCtx      *user;
786:   PetscFunctionList rlist      = 0,rclist = 0;
787:   char              rname[256] = "exact",rcname[256] = "characteristic";

790:   PetscNew(&user);
791:   ctx->physics.sample         = PhysicsSample_Acoustics;
792:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
793:   ctx->physics.user           = user;
794:   ctx->physics.dof            = 2;

796:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
797:   PetscStrallocpy("v",&ctx->physics.fieldname[1]);

799:   user->c = 1;
800:   user->z = 1;

802:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Acoustics_Exact);
803:   ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Acoustics);
804:   ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
805:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for linear Acoustics","");
806:   {
807:     PetscOptionsReal("-physics_acoustics_c","c = sqrt(bulk/rho)","",user->c,&user->c,NULL);
808:     PetscOptionsReal("-physics_acoustics_z","z = sqrt(bulk*rho)","",user->z,&user->z,NULL);
809:     PetscOptionsFList("-physics_acoustics_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
810:     PetscOptionsFList("-physics_acoustics_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
811:   }
812:   PetscOptionsEnd();
813:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
814:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
815:   PetscFunctionListDestroy(&rlist);
816:   PetscFunctionListDestroy(&rclist);
817:   return(0);
818: }

820: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */

822: typedef struct {
823:   PetscReal acoustic_speed;
824: } IsoGasCtx;

826: PETSC_STATIC_INLINE void IsoGasFlux(PetscReal c,const PetscScalar *u,PetscScalar *f)
827: {
828:   f[0] = u[1];
829:   f[1] = PetscSqr(u[1])/u[0] + c*c*u[0];
830: }

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

838:   if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solutions not implemented for t > 0");
839:   switch (initial) {
840:     case 0:
841:       u[0] = (x < 0) ? 1 : 0.5;
842:       u[1] = (x < 0) ? 1 : 0.7;
843:       break;
844:     case 1:
845:       u[0] = 1+0.5*PetscSinReal(2*PETSC_PI*x);
846:       u[1] = 1*u[0];
847:       break;
848:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
849:   }
850:   return(0);
851: }

855: static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
856: {
857:   IsoGasCtx   *phys = (IsoGasCtx*)vctx;
858:   PetscReal   c     = phys->acoustic_speed;
859:   PetscScalar ubar,du[2],a[2],fL[2],fR[2],lam[2],ustar[2],R[2][2];
860:   PetscInt    i;

863:   ubar = (uL[1]/PetscSqrtScalar(uL[0]) + uR[1]/PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0]));
864:   /* write fluxuations in characteristic basis */
865:   du[0] = uR[0] - uL[0];
866:   du[1] = uR[1] - uL[1];
867:   a[0]  = (1/(2*c)) * ((ubar + c)*du[0] - du[1]);
868:   a[1]  = (1/(2*c)) * ((-ubar + c)*du[0] + du[1]);
869:   /* wave speeds */
870:   lam[0] = ubar - c;
871:   lam[1] = ubar + c;
872:   /* Right eigenvectors */
873:   R[0][0] = 1; R[0][1] = ubar-c;
874:   R[1][0] = 1; R[1][1] = ubar+c;
875:   /* Compute state in star region (between the 1-wave and 2-wave) */
876:   for (i=0; i<2; i++) ustar[i] = uL[i] + a[0]*R[0][i];
877:   if (uL[1]/uL[0] < c && c < ustar[1]/ustar[0]) { /* 1-wave is sonic rarefaction */
878:     PetscScalar ufan[2];
879:     ufan[0] = uL[0]*PetscExpScalar(uL[1]/(uL[0]*c) - 1);
880:     ufan[1] = c*ufan[0];
881:     IsoGasFlux(c,ufan,flux);
882:   } else if (ustar[1]/ustar[0] < -c && -c < uR[1]/uR[0]) { /* 2-wave is sonic rarefaction */
883:     PetscScalar ufan[2];
884:     ufan[0] = uR[0]*PetscExpScalar(-uR[1]/(uR[0]*c) - 1);
885:     ufan[1] = -c*ufan[0];
886:     IsoGasFlux(c,ufan,flux);
887:   } else {                      /* Centered form */
888:     IsoGasFlux(c,uL,fL);
889:     IsoGasFlux(c,uR,fR);
890:     for (i=0; i<2; i++) {
891:       PetscScalar absdu = PetscAbsScalar(lam[0])*a[0]*R[0][i] + PetscAbsScalar(lam[1])*a[1]*R[1][i];
892:       flux[i] = 0.5*(fL[i]+fR[i]) - 0.5*absdu;
893:     }
894:   }
895:   *maxspeed = MaxAbs(lam[0],lam[1]);
896:   return(0);
897: }

901: static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
902: {
903:   IsoGasCtx                   *phys = (IsoGasCtx*)vctx;
904:   PetscReal                   c     = phys->acoustic_speed;
905:   PetscScalar                 ustar[2];
906:   struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
907:   PetscInt                    i;
908:   PetscErrorCode              ierr;

911:   if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
912:   {
913:     /* Solve for star state */
914:     PetscScalar res,tmp,rho = 0.5*(L.rho + R.rho); /* initial guess */
915:     for (i=0; i<20; i++) {
916:       PetscScalar fr,fl,dfr,dfl;
917:       fl = (L.rho < rho)
918:         ? (rho-L.rho)/PetscSqrtScalar(L.rho*rho)       /* shock */
919:         : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
920:       fr = (R.rho < rho)
921:         ? (rho-R.rho)/PetscSqrtScalar(R.rho*rho)       /* shock */
922:         : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
923:       res = R.u-L.u + c*(fr+fl);
924:       PetscIsInfOrNanScalar(res);
925:       if (PetscAbsScalar(res) < 1e-10) {
926:         star.rho = rho;
927:         star.u   = L.u - c*fl;
928:         goto converged;
929:       }
930:       dfl = (L.rho < rho)
931:         ? 1/PetscSqrtScalar(L.rho*rho)*(1 - 0.5*(rho-L.rho)/rho)
932:         : 1/rho;
933:       dfr = (R.rho < rho)
934:         ? 1/PetscSqrtScalar(R.rho*rho)*(1 - 0.5*(rho-R.rho)/rho)
935:         : 1/rho;
936:       tmp = rho - res/(c*(dfr+dfl));
937:       if (tmp <= 0) rho /= 2;   /* Guard against Newton shooting off to a negative density */
938:       else rho = tmp;
939:       if (!((rho > 0) && PetscIsNormalScalar(rho))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate rho=%g",(double)PetscRealPart(rho));
940:     }
941:     SETERRQ1(PETSC_COMM_SELF,1,"Newton iteration for star.rho diverged after %D iterations",i);
942:   }
943: converged:
944:   if (L.u-c < 0 && 0 < star.u-c) { /* 1-wave is sonic rarefaction */
945:     PetscScalar ufan[2];
946:     ufan[0] = L.rho*PetscExpScalar(L.u/c - 1);
947:     ufan[1] = c*ufan[0];
948:     IsoGasFlux(c,ufan,flux);
949:   } else if (star.u+c < 0 && 0 < R.u+c) { /* 2-wave is sonic rarefaction */
950:     PetscScalar ufan[2];
951:     ufan[0] = R.rho*PetscExpScalar(-R.u/c - 1);
952:     ufan[1] = -c*ufan[0];
953:     IsoGasFlux(c,ufan,flux);
954:   } else if ((L.rho >= star.rho && L.u-c >= 0)
955:              || (L.rho < star.rho && (star.rho*star.u-L.rho*L.u)/(star.rho-L.rho) > 0)) {
956:     /* 1-wave is supersonic rarefaction, or supersonic shock */
957:     IsoGasFlux(c,uL,flux);
958:   } else if ((star.rho <= R.rho && R.u+c <= 0)
959:              || (star.rho > R.rho && (R.rho*R.u-star.rho*star.u)/(R.rho-star.rho) < 0)) {
960:     /* 2-wave is supersonic rarefaction or supersonic shock */
961:     IsoGasFlux(c,uR,flux);
962:   } else {
963:     ustar[0] = star.rho;
964:     ustar[1] = star.rho*star.u;
965:     IsoGasFlux(c,ustar,flux);
966:   }
967:   *maxspeed = MaxAbs(MaxAbs(star.u-c,star.u+c),MaxAbs(L.u-c,R.u+c));
968:   return(0);
969: }

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

980:   if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
981:   IsoGasFlux(c,uL,fL);
982:   IsoGasFlux(c,uR,fR);
983:   s         = PetscMax(PetscAbs(L.u),PetscAbs(R.u))+c;
984:   flux[0]   = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
985:   flux[1]   = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
986:   *maxspeed = s;
987:   return(0);
988: }

992: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
993: {
994:   IsoGasCtx      *phys = (IsoGasCtx*)vctx;
995:   PetscReal      c     = phys->acoustic_speed;

999:   speeds[0] = u[1]/u[0] - c;
1000:   speeds[1] = u[1]/u[0] + c;
1001:   X[0*2+0]  = 1;
1002:   X[0*2+1]  = speeds[0];
1003:   X[1*2+0]  = 1;
1004:   X[1*2+1]  = speeds[1];

1006:   PetscMemcpy(Xi,X,4*sizeof(X[0]));
1007:   PetscKernel_A_gets_inverse_A_2(Xi,0);
1008:   return(0);
1009: }

1013: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
1014: {
1015:   PetscErrorCode    ierr;
1016:   IsoGasCtx         *user;
1017:   PetscFunctionList rlist = 0,rclist = 0;
1018:   char              rname[256] = "exact",rcname[256] = "characteristic";

1021:   PetscNew(&user);
1022:   ctx->physics.sample         = PhysicsSample_IsoGas;
1023:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
1024:   ctx->physics.user           = user;
1025:   ctx->physics.dof            = 2;

1027:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1028:   PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);

1030:   user->acoustic_speed = 1;

1032:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_IsoGas_Exact);
1033:   RiemannListAdd(&rlist,"roe",    PhysicsRiemann_IsoGas_Roe);
1034:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_IsoGas_Rusanov);
1035:   ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_IsoGas);
1036:   ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1037:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");
1038:   {
1039:     PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,NULL);
1040:     PetscOptionsFList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
1041:     PetscOptionsFList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
1042:   }
1043:   PetscOptionsEnd();
1044:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
1045:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1046:   PetscFunctionListDestroy(&rlist);
1047:   PetscFunctionListDestroy(&rclist);
1048:   return(0);
1049: }



1053: /* --------------------------------- Shallow Water ----------------------------------- */

1055: typedef struct {
1056:   PetscReal gravity;
1057: } ShallowCtx;

1059: PETSC_STATIC_INLINE void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
1060: {
1061:   f[0] = u[1];
1062:   f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
1063: }

1067: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
1068: {
1069:   ShallowCtx                *phys = (ShallowCtx*)vctx;
1070:   PetscScalar               g    = phys->gravity,ustar[2],cL,cR,c,cstar;
1071:   struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
1072:   PetscInt                  i;
1073:   PetscErrorCode            ierr;

1076:   if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed thickness is negative");
1077:   cL = PetscSqrtScalar(g*L.h);
1078:   cR = PetscSqrtScalar(g*R.h);
1079:   c  = PetscMax(cL,cR);
1080:   {
1081:     /* Solve for star state */
1082:     const PetscInt maxits = 50;
1083:     PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */
1084:     h0 = h;
1085:     for (i=0; i<maxits; i++) {
1086:       PetscScalar fr,fl,dfr,dfl;
1087:       fl = (L.h < h)
1088:         ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */
1089:         : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h);   /* rarefaction */
1090:       fr = (R.h < h)
1091:         ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */
1092:         : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h);   /* rarefaction */
1093:       res = R.u - L.u + fr + fl;
1094:       PetscIsInfOrNanScalar(res);
1095:       if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h-h0) < 1e-8)) {
1096:         star.h = h;
1097:         star.u = L.u - fl;
1098:         goto converged;
1099:       } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) {        /* Line search */
1100:         h = 0.8*h0 + 0.2*h;
1101:         continue;
1102:       }
1103:       /* Accept the last step and take another */
1104:       res0 = res;
1105:       h0 = h;
1106:       dfl = (L.h < h)
1107:         ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h)
1108:         : PetscSqrtScalar(g/h);
1109:       dfr = (R.h < h)
1110:         ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h)
1111:         : PetscSqrtScalar(g/h);
1112:       tmp = h - res/(dfr+dfl);
1113:       if (tmp <= 0) h /= 2;   /* Guard against Newton shooting off to a negative thickness */
1114:       else h = tmp;
1115:       if (!((h > 0) && PetscIsNormalScalar(h))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate h=%g",(double)h);
1116:     }
1117:     SETERRQ1(PETSC_COMM_SELF,1,"Newton iteration for star.h diverged after %D iterations",i);
1118:   }
1119: converged:
1120:   cstar = PetscSqrtScalar(g*star.h);
1121:   if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */
1122:     PetscScalar ufan[2];
1123:     ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL);
1124:     ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0];
1125:     ShallowFlux(phys,ufan,flux);
1126:   } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */
1127:     PetscScalar ufan[2];
1128:     ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR);
1129:     ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0];
1130:     ShallowFlux(phys,ufan,flux);
1131:   } else if ((L.h >= star.h && L.u-c >= 0)
1132:              || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) {
1133:     /* 1-wave is right-travelling shock (supersonic) */
1134:     ShallowFlux(phys,uL,flux);
1135:   } else if ((star.h <= R.h && R.u+c <= 0)
1136:              || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) {
1137:     /* 2-wave is left-travelling shock (supersonic) */
1138:     ShallowFlux(phys,uR,flux);
1139:   } else {
1140:     ustar[0] = star.h;
1141:     ustar[1] = star.h*star.u;
1142:     ShallowFlux(phys,ustar,flux);
1143:   }
1144:   *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR));
1145:   return(0);
1146: }

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

1157:   if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed thickness is negative");
1158:   ShallowFlux(phys,uL,fL);
1159:   ShallowFlux(phys,uR,fR);
1160:   s         = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h));
1161:   flux[0]   = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
1162:   flux[1]   = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
1163:   *maxspeed = s;
1164:   return(0);
1165: }

1169: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
1170: {
1171:   ShallowCtx     *phys = (ShallowCtx*)vctx;
1172:   PetscReal      c;

1176:   c         = PetscSqrtScalar(u[0]*phys->gravity);
1177:   speeds[0] = u[1]/u[0] - c;
1178:   speeds[1] = u[1]/u[0] + c;
1179:   X[0*2+0]  = 1;
1180:   X[0*2+1]  = speeds[0];
1181:   X[1*2+0]  = 1;
1182:   X[1*2+1]  = speeds[1];

1184:   PetscMemcpy(Xi,X,4*sizeof(X[0]));
1185:   PetscKernel_A_gets_inverse_A_2(Xi,0);
1186:   return(0);
1187: }

1191: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
1192: {
1193:   PetscErrorCode    ierr;
1194:   ShallowCtx        *user;
1195:   PetscFunctionList rlist = 0,rclist = 0;
1196:   char              rname[256] = "exact",rcname[256] = "characteristic";

1199:   PetscNew(&user);
1200:   /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
1201:   ctx->physics.sample         = PhysicsSample_IsoGas;
1202:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
1203:   ctx->physics.user           = user;
1204:   ctx->physics.dof            = 2;

1206:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1207:   PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);

1209:   user->gravity = 1;

1211:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Shallow_Exact);
1212:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov);
1213:   ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Shallow);
1214:   ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1215:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");
1216:   {
1217:     PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,NULL);
1218:     PetscOptionsFList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
1219:     PetscOptionsFList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
1220:   }
1221:   PetscOptionsEnd();
1222:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
1223:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1224:   PetscFunctionListDestroy(&rlist);
1225:   PetscFunctionListDestroy(&rclist);
1226:   return(0);
1227: }



1231: /* --------------------------------- Finite Volume Solver ----------------------------------- */

1235: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
1236: {
1237:   FVCtx          *ctx = (FVCtx*)vctx;
1239:   PetscInt       i,j,k,Mx,dof,xs,xm;
1240:   PetscReal      hx,cfl_idt = 0;
1241:   PetscScalar    *x,*f,*slope;
1242:   Vec            Xloc;
1243:   DM             da;

1246:   TSGetDM(ts,&da);
1247:   DMGetLocalVector(da,&Xloc);
1248:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1249:   hx   = (ctx->xmax - ctx->xmin)/Mx;
1250:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1251:   DMGlobalToLocalEnd  (da,X,INSERT_VALUES,Xloc);

1253:   VecZeroEntries(F);

1255:   DMDAVecGetArray(da,Xloc,&x);
1256:   DMDAVecGetArray(da,F,&f);
1257:   DMDAGetArray(da,PETSC_TRUE,&slope);

1259:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

1261:   if (ctx->bctype == FVBC_OUTFLOW) {
1262:     for (i=xs-2; i<0; i++) {
1263:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
1264:     }
1265:     for (i=Mx; i<xs+xm+2; i++) {
1266:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
1267:     }
1268:   }
1269:   for (i=xs-1; i<xs+xm+1; i++) {
1270:     struct _LimitInfo info;
1271:     PetscScalar       *cjmpL,*cjmpR;
1272:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
1273:     (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1274:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
1275:     PetscMemzero(ctx->cjmpLR,2*dof*sizeof(ctx->cjmpLR[0]));
1276:     cjmpL = &ctx->cjmpLR[0];
1277:     cjmpR = &ctx->cjmpLR[dof];
1278:     for (j=0; j<dof; j++) {
1279:       PetscScalar jmpL,jmpR;
1280:       jmpL = x[(i+0)*dof+j] - x[(i-1)*dof+j];
1281:       jmpR = x[(i+1)*dof+j] - x[(i+0)*dof+j];
1282:       for (k=0; k<dof; k++) {
1283:         cjmpL[k] += ctx->Rinv[k+j*dof] * jmpL;
1284:         cjmpR[k] += ctx->Rinv[k+j*dof] * jmpR;
1285:       }
1286:     }
1287:     /* Apply limiter to the left and right characteristic jumps */
1288:     info.m  = dof;
1289:     info.hx = hx;
1290:     (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
1291:     for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
1292:     for (j=0; j<dof; j++) {
1293:       PetscScalar tmp = 0;
1294:       for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof] * ctx->cslope[k];
1295:       slope[i*dof+j] = tmp;
1296:     }
1297:   }

1299:   for (i=xs; i<xs+xm+1; i++) {
1300:     PetscReal   maxspeed;
1301:     PetscScalar *uL,*uR;
1302:     uL = &ctx->uLR[0];
1303:     uR = &ctx->uLR[dof];
1304:     for (j=0; j<dof; j++) {
1305:       uL[j] = x[(i-1)*dof+j] + slope[(i-1)*dof+j]*hx/2;
1306:       uR[j] = x[(i-0)*dof+j] - slope[(i-0)*dof+j]*hx/2;
1307:     }
1308:     (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed);
1309:     cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */

1311:     if (i > xs) {
1312:       for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
1313:     }
1314:     if (i < xs+xm) {
1315:       for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
1316:     }
1317:   }

1319:   DMDAVecRestoreArray(da,Xloc,&x);
1320:   DMDAVecRestoreArray(da,F,&f);
1321:   DMDARestoreArray(da,PETSC_TRUE,&slope);
1322:   DMRestoreLocalVector(da,&Xloc);

1324:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));
1325:   if (0) {
1326:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
1327:     PetscReal dt,tnow;
1328:     TSGetTimeStep(ts,&dt);
1329:     TSGetTime(ts,&tnow);
1330:     if (dt > 0.5/ctx->cfl_idt) {
1331:       if (1) {
1332:         PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
1333:       } else SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt));
1334:     }
1335:   }
1336:   return(0);
1337: }

1341: static PetscErrorCode SmallMatMultADB(PetscScalar *C,PetscInt bs,const PetscScalar *A,const PetscReal *D,const PetscScalar *B)
1342: {
1343:   PetscInt i,j,k;

1346:   for (i=0; i<bs; i++) {
1347:     for (j=0; j<bs; j++) {
1348:       PetscScalar tmp = 0;
1349:       for (k=0; k<bs; k++) tmp += A[i*bs+k] * D[k] * B[k*bs+j];
1350:       C[i*bs+j] = tmp;
1351:     }
1352:   }
1353:   return(0);
1354: }


1359: static PetscErrorCode FVIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void *vctx)
1360: {
1361:   FVCtx          *ctx = (FVCtx*)vctx;
1363:   PetscInt       i,j,dof = ctx->physics.dof;
1364:   PetscScalar    *x,*J;
1365:   PetscReal      hx;
1366:   DM             da;
1367:   DMDALocalInfo  dainfo;

1370:   TSGetDM(ts,&da);
1371:   DMDAVecGetArray(da,X,&x);
1372:   DMDAGetLocalInfo(da,&dainfo);
1373:   hx   = (ctx->xmax - ctx->xmin)/dainfo.mx;
1374:   PetscMalloc1(dof*dof,&J);
1375:   for (i=dainfo.xs; i<dainfo.xs+dainfo.xm; i++) {
1376:     (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1377:     for (j=0; j<dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]);
1378:     SmallMatMultADB(J,dof,ctx->R,ctx->speeds,ctx->Rinv);
1379:     for (j=0; j<dof*dof; j++) J[j] = J[j]/hx + shift*(j/dof == j%dof);
1380:     MatSetValuesBlocked(B,1,&i,1,&i,J,INSERT_VALUES);
1381:   }
1382:   PetscFree(J);
1383:   DMDAVecRestoreArray(da,X,&x);

1385:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1386:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1387:   if (A != B) {
1388:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1389:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1390:   }
1391:   return(0);
1392: }

1396: static PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
1397: {
1399:   PetscScalar    *u,*uj;
1400:   PetscInt       i,j,k,dof,xs,xm,Mx;

1403:   if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,1,"Physics has not provided a sampling function");
1404:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1405:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1406:   DMDAVecGetArray(da,U,&u);
1407:   PetscMalloc1(dof,&uj);
1408:   for (i=xs; i<xs+xm; i++) {
1409:     const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h;
1410:     const PetscInt  N = 200;
1411:     /* Integrate over cell i using trapezoid rule with N points. */
1412:     for (k=0; k<dof; k++) u[i*dof+k] = 0;
1413:     for (j=0; j<N+1; j++) {
1414:       PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N;
1415:       (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
1416:       for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
1417:     }
1418:   }
1419:   DMDAVecRestoreArray(da,U,&u);
1420:   PetscFree(uj);
1421:   return(0);
1422: }

1426: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
1427: {
1429:   PetscReal      xmin,xmax;
1430:   PetscScalar    sum,*x,tvsum,tvgsum;
1431:   PetscInt       imin,imax,Mx,i,j,xs,xm,dof;
1432:   Vec            Xloc;
1433:   PetscBool      iascii;

1436:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1437:   if (iascii) {
1438:     /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
1439:     DMGetLocalVector(da,&Xloc);
1440:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1441:     DMGlobalToLocalEnd  (da,X,INSERT_VALUES,Xloc);
1442:     DMDAVecGetArray(da,Xloc,&x);
1443:     DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1444:     DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1445:     tvsum = 0;
1446:     for (i=xs; i<xs+xm; i++) {
1447:       for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j] - x[(i-1)*dof+j]);
1448:     }
1449:     MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));
1450:     DMDAVecRestoreArray(da,Xloc,&x);
1451:     DMRestoreLocalVector(da,&Xloc);

1453:     VecMin(X,&imin,&xmin);
1454:     VecMax(X,&imax,&xmax);
1455:     VecSum(X,&sum);
1456:     PetscViewerASCIIPrintf(viewer,"Solution range [%8.5f,%8.5f] with extrema at %D and %D, mean %8.5f, ||x||_TV %8.5f\n",(double)xmin,(double)xmax,imin,imax,(double)(sum/Mx),(double)(tvgsum/Mx));
1457:   } else SETERRQ(PETSC_COMM_SELF,1,"Viewer type not supported");
1458:   return(0);
1459: }

1463: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1,PetscReal *nrmsup)
1464: {
1466:   Vec            Y;
1467:   PetscInt       Mx;

1470:   VecGetSize(X,&Mx);
1471:   VecDuplicate(X,&Y);
1472:   FVSample(ctx,da,t,Y);
1473:   VecAYPX(Y,-1,X);
1474:   VecNorm(Y,NORM_1,nrm1);
1475:   VecNorm(Y,NORM_INFINITY,nrmsup);
1476:   *nrm1 /= Mx;
1477:   VecDestroy(&Y);
1478:   return(0);
1479: }

1483: int main(int argc,char *argv[])
1484: {
1485:   char              lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
1486:   PetscFunctionList limiters   = 0,physics = 0;
1487:   MPI_Comm          comm;
1488:   TS                ts;
1489:   DM                da;
1490:   Vec               X,X0,R;
1491:   Mat               B;
1492:   FVCtx             ctx;
1493:   PetscInt          i,dof,xs,xm,Mx,draw = 0;
1494:   PetscBool         view_final = PETSC_FALSE;
1495:   PetscReal         ptime;
1496:   PetscErrorCode    ierr;

1498:   PetscInitialize(&argc,&argv,0,help);
1499:   comm = PETSC_COMM_WORLD;
1500:   PetscMemzero(&ctx,sizeof(ctx));

1502:   /* Register limiters to be available on the command line */
1503:   PetscFunctionListAdd(&limiters,"upwind"              ,Limit_Upwind);
1504:   PetscFunctionListAdd(&limiters,"lax-wendroff"        ,Limit_LaxWendroff);
1505:   PetscFunctionListAdd(&limiters,"beam-warming"        ,Limit_BeamWarming);
1506:   PetscFunctionListAdd(&limiters,"fromm"               ,Limit_Fromm);
1507:   PetscFunctionListAdd(&limiters,"minmod"              ,Limit_Minmod);
1508:   PetscFunctionListAdd(&limiters,"superbee"            ,Limit_Superbee);
1509:   PetscFunctionListAdd(&limiters,"mc"                  ,Limit_MC);
1510:   PetscFunctionListAdd(&limiters,"vanleer"             ,Limit_VanLeer);
1511:   PetscFunctionListAdd(&limiters,"vanalbada"           ,Limit_VanAlbada);
1512:   PetscFunctionListAdd(&limiters,"vanalbadatvd"        ,Limit_VanAlbadaTVD);
1513:   PetscFunctionListAdd(&limiters,"koren"               ,Limit_Koren);
1514:   PetscFunctionListAdd(&limiters,"korensym"            ,Limit_KorenSym);
1515:   PetscFunctionListAdd(&limiters,"koren3"              ,Limit_Koren3);
1516:   PetscFunctionListAdd(&limiters,"cada-torrilhon2"     ,Limit_CadaTorrilhon2);
1517:   PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1);
1518:   PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1"  ,Limit_CadaTorrilhon3R1);
1519:   PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10);
1520:   PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100);

1522:   /* Register physical models to be available on the command line */
1523:   PetscFunctionListAdd(&physics,"advect"          ,PhysicsCreate_Advect);
1524:   PetscFunctionListAdd(&physics,"burgers"         ,PhysicsCreate_Burgers);
1525:   PetscFunctionListAdd(&physics,"traffic"         ,PhysicsCreate_Traffic);
1526:   PetscFunctionListAdd(&physics,"acoustics"       ,PhysicsCreate_Acoustics);
1527:   PetscFunctionListAdd(&physics,"isogas"          ,PhysicsCreate_IsoGas);
1528:   PetscFunctionListAdd(&physics,"shallow"         ,PhysicsCreate_Shallow);

1530:   ctx.comm = comm;
1531:   ctx.cfl  = 0.9; ctx.bctype = FVBC_PERIODIC;
1532:   ctx.xmin = -1; ctx.xmax = 1;
1533:   PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
1534:   {
1535:     PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
1536:     PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
1537:     PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);
1538:     PetscOptionsFList("-physics","Name of physics (Riemann solver and characteristics) to use","",physics,physname,physname,sizeof(physname),NULL);
1539:     PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
1540:     PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
1541:     PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
1542:     PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
1543:     PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
1544:     PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
1545:   }
1546:   PetscOptionsEnd();

1548:   /* Choose the limiter from the list of registered limiters */
1549:   PetscFunctionListFind(limiters,lname,&ctx.limit);
1550:   if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);

1552:   /* Choose the physics from the list of registered models */
1553:   {
1554:     PetscErrorCode (*r)(FVCtx*);
1555:     PetscFunctionListFind(physics,physname,&r);
1556:     if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
1557:     /* Create the physics, will set the number of fields and their names */
1558:     (*r)(&ctx);
1559:   }

1561:   /* Create a DMDA to manage the parallel grid */
1562:   DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,-50,ctx.physics.dof,2,NULL,&da);
1563:   /* Inform the DMDA of the field names provided by the physics. */
1564:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1565:   for (i=0; i<ctx.physics.dof; i++) {
1566:     DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
1567:   }
1568:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1569:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

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

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

1578:   /* Create a vector to store the solution and to save the initial state */
1579:   DMCreateGlobalVector(da,&X);
1580:   VecDuplicate(X,&X0);
1581:   VecDuplicate(X,&R);

1583:   DMCreateMatrix(da,&B);

1585:   /* Create a time-stepping object */
1586:   TSCreate(comm,&ts);
1587:   TSSetDM(ts,da);
1588:   TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
1589:   TSSetIJacobian(ts,B,B,FVIJacobian,&ctx);
1590:   TSSetType(ts,TSSSP);
1591:   TSSetDuration(ts,1000,10);

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

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

1601:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);

1603:   {
1604:     PetscReal nrm1,nrmsup;
1605:     PetscInt  steps;

1607:     TSSolve(ts,X);
1608:     TSGetSolveTime(ts,&ptime);
1609:     TSGetTimeStepNumber(ts,&steps);

1611:     PetscPrintf(comm,"Final time %8.5f, steps %D\n",(double)ptime,steps);
1612:     if (ctx.exact) {
1613:       SolutionErrorNorms(&ctx,da,ptime,X,&nrm1,&nrmsup);
1614:       PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e  ||x-x_e||_sup %8.4e\n",(double)nrm1,(double)nrmsup);
1615:     }
1616:   }

1618:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);

1620:   if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
1621:   if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
1622:   if (draw & 0x4) {
1623:     Vec Y;
1624:     VecDuplicate(X,&Y);
1625:     FVSample(&ctx,da,ptime,Y);
1626:     VecAYPX(Y,-1,X);
1627:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1628:     VecDestroy(&Y);
1629:   }

1631:   if (view_final) {
1632:     PetscViewer viewer;
1633:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1634:     PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1635:     VecView(X,viewer);
1636:     PetscViewerDestroy(&viewer);
1637:   }

1639:   /* Clean up */
1640:   (*ctx.physics.destroy)(ctx.physics.user);
1641:   for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
1642:   PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
1643:   PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
1644:   VecDestroy(&X);
1645:   VecDestroy(&X0);
1646:   VecDestroy(&R);
1647:   MatDestroy(&B);
1648:   DMDestroy(&da);
1649:   TSDestroy(&ts);
1650:   PetscFunctionListDestroy(&limiters);
1651:   PetscFunctionListDestroy(&physics);
1652:   PetscFinalize();
1653:   return 0;
1654: }