Actual source code: ex9.c

petsc-3.3-p7 2013-05-11
  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: /* To get isfinite in math.h */
 34: #define _XOPEN_SOURCE 600

 36: #include <petscts.h>
 37: #include <petscdmda.h>

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

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

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


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

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


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

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

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

184: typedef struct {
185:   void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*);
186:   PhysicsCtx physics;

188:   MPI_Comm comm;
189:   char prefix[256];

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

199:   PetscReal cfl_idt;            /* Max allowable value of 1/Delta t */
200:   PetscReal cfl;
201:   PetscReal xmin,xmax;
202:   PetscInt initial;
203:   PetscBool  exact;
204:   FVBCType bctype;
205: } FVCtx;


208: /* Utility */

212: PetscErrorCode RiemannListAdd(PetscFList *flist,const char *name,RiemannFunction rsolve)
213: {

217:   PetscFListAdd(flist,name,"",(void(*)(void))rsolve);
218:   return(0);
219: }

223: PetscErrorCode RiemannListFind(PetscFList flist,const char *name,RiemannFunction *rsolve)
224: {

228:   PetscFListFind(flist,PETSC_COMM_WORLD,name,PETSC_FALSE,(void(**)(void))rsolve);
229:   if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,1,"Riemann solver \"%s\" could not be found",name);
230:   return(0);
231: }

235: PetscErrorCode ReconstructListAdd(PetscFList *flist,const char *name,ReconstructFunction r)
236: {

240:   PetscFListAdd(flist,name,"",(void(*)(void))r);
241:   return(0);
242: }

246: PetscErrorCode ReconstructListFind(PetscFList flist,const char *name,ReconstructFunction *r)
247: {

251:   PetscFListFind(flist,PETSC_COMM_WORLD,name,PETSC_FALSE,(void(**)(void))r);
252:   if (!*r) SETERRQ1(PETSC_COMM_SELF,1,"Reconstruction \"%s\" could not be found",name);
253:   return(0);
254: }


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

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

272:   for (i=0; i<m; i++) {
273:     for (j=0; j<m; j++) {
274:       Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j);
275:     }
276:     speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
277:   }
278:   return(0);
279: }

283: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
284: {

288:   PetscFree(vctx);
289:   return(0);
290: }



294: /* --------------------------------- Advection ----------------------------------- */

296: typedef struct {
297:   PetscReal a;                  /* advective velocity */
298: } AdvectCtx;

302: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
303: {
304:   AdvectCtx *ctx = (AdvectCtx*)vctx;
305:   PetscReal speed;

308:   speed = ctx->a;
309:   flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
310:   *maxspeed = speed;
311:   return(0);
312: }

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

321:   X[0] = 1.;
322:   Xi[0] = 1.;
323:   speeds[0] = ctx->a;
324:   return(0);
325: }

329: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
330: {
331:   AdvectCtx *ctx = (AdvectCtx*)vctx;
332:   PetscReal a = ctx->a,x0;

335:   switch (bctype) {
336:     case FVBC_OUTFLOW: x0 = x-a*t; break;
337:     case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
338:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
339:   }
340:   switch (initial) {
341:     case 0: u[0] = (x0 < 0) ? 1 : -1; break;
342:     case 1: u[0] = (x0 < 0) ? -1 : 1; break;
343:     case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
344:     case 3: u[0] = sin(2*M_PI*x0); break;
345:     case 4: u[0] = PetscAbs(x0); break;
346:     case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(sin(2*M_PI*x0)); break;
347:     case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
348:     default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
349:   }
350:   return(0);
351: }

355: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
356: {
358:   AdvectCtx *user;

361:   PetscNew(*user,&user);
362:   ctx->physics.sample         = PhysicsSample_Advect;
363:   ctx->physics.riemann        = PhysicsRiemann_Advect;
364:   ctx->physics.characteristic = PhysicsCharacteristic_Advect;
365:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
366:   ctx->physics.user           = user;
367:   ctx->physics.dof            = 1;
368:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
369:   user->a = 1;
370:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
371:   {
372:     PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,PETSC_NULL);
373:   }
374:   PetscOptionsEnd();
375:   return(0);
376: }



380: /* --------------------------------- Burgers ----------------------------------- */

382: typedef struct {
383:   PetscReal lxf_speed;
384: } BurgersCtx;

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

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

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

430:   if (uL[0] < uR[0]) {          /* rarefaction */
431:     flux[0] = (uL[0]*uR[0] < 0)
432:       ? 0                       /* sonic rarefaction */
433:       : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0]));
434:   } else {                      /* shock */
435:     flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0]));
436:   }
437:   *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
438:   return(0);
439: }

443: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
444: {
445:   PetscReal speed;

448:   speed = 0.5*(uL[0] + uR[0]);
449:   flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
450:   if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
451:   *maxspeed = speed;
452:   return(0);
453: }

457: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
458: {
459:   PetscReal c;
460:   PetscScalar fL,fR;

463:   c = ((BurgersCtx*)vctx)->lxf_speed;
464:   fL = 0.5*PetscSqr(uL[0]);
465:   fR = 0.5*PetscSqr(uR[0]);
466:   flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
467:   *maxspeed = c;
468:   return(0);
469: }

473: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
474: {
475:   PetscReal c;
476:   PetscScalar fL,fR;

479:   c = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0]));
480:   fL = 0.5*PetscSqr(uL[0]);
481:   fR = 0.5*PetscSqr(uR[0]);
482:   flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
483:   *maxspeed = c;
484:   return(0);
485: }

489: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
490: {
491:   BurgersCtx *user;
493:   RiemannFunction r;
494:   PetscFList rlist = 0;
495:   char rname[256] = "exact";

498:   PetscNew(*user,&user);
499:   ctx->physics.sample         = PhysicsSample_Burgers;
500:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
501:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
502:   ctx->physics.user           = user;
503:   ctx->physics.dof            = 1;
504:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
505:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Burgers_Exact);
506:   RiemannListAdd(&rlist,"roe",    PhysicsRiemann_Burgers_Roe);
507:   RiemannListAdd(&rlist,"lxf",    PhysicsRiemann_Burgers_LxF);
508:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Burgers_Rusanov);
509:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
510:   {
511:     PetscOptionsList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
512:   }
513:   PetscOptionsEnd();
514:   RiemannListFind(rlist,rname,&r);
515:   PetscFListDestroy(&rlist);
516:   ctx->physics.riemann = r;

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



528: /* --------------------------------- Traffic ----------------------------------- */

530: typedef struct {
531:   PetscReal lxf_speed;
532:   PetscReal a;
533: } TrafficCtx;

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

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

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

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

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

584: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
585: {
586:   PetscReal a = ((TrafficCtx*)vctx)->a;
587:   PetscReal speed;

590:   speed = a*(1 - (uL[0] + uR[0]));
591:   flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
592:   *maxspeed = speed;
593:   return(0);
594: }

598: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
599: {
600:   TrafficCtx *phys = (TrafficCtx*)vctx;
601:   PetscReal a = phys->a;
602:   PetscReal speed;

605:   speed = a*(1 - (uL[0] + uR[0]));
606:   flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]);
607:   *maxspeed = speed;
608:   return(0);
609: }

613: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
614: {
615:   PetscReal a = ((TrafficCtx*)vctx)->a;
616:   PetscReal speed;

619:   speed = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0]));
620:   flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]);
621:   *maxspeed = speed;
622:   return(0);
623: }

627: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
628: {
630:   TrafficCtx *user;
631:   RiemannFunction r;
632:   PetscFList rlist = 0;
633:   char rname[256] = "exact";

636:   PetscNew(*user,&user);
637:   ctx->physics.sample         = PhysicsSample_Traffic;
638:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
639:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
640:   ctx->physics.user           = user;
641:   ctx->physics.dof            = 1;
642:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
643:   user->a = 0.5;
644:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Traffic_Exact);
645:   RiemannListAdd(&rlist,"roe",    PhysicsRiemann_Traffic_Roe);
646:   RiemannListAdd(&rlist,"lxf",    PhysicsRiemann_Traffic_LxF);
647:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Traffic_Rusanov);
648:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");
649:   {
650:     PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,PETSC_NULL);
651:     PetscOptionsList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
652:   }
653:   PetscOptionsEnd();

655:   RiemannListFind(rlist,rname,&r);
656:   PetscFListDestroy(&rlist);
657:   ctx->physics.riemann = r;

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

665:   return(0);
666: }

668: /* --------------------------------- Linear Acoustics ----------------------------------- */

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

681: typedef struct {
682:   PetscReal c;                  /* speed of sound: c = sqrt(bulk/rho) */
683:   PetscReal z;                  /* impedence: z = sqrt(rho*bulk) */
684: } AcousticsCtx;

686: static inline void AcousticsFlux(AcousticsCtx *ctx,const PetscScalar *u,PetscScalar *f)
687: {
688:   f[0] = ctx->c*ctx->z*u[1];
689:   f[1] = ctx->c/ctx->z*u[0];
690: }

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

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

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

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

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

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

777:   flux[0] = Al[0][0]*uR[0] + Al[0][1]*uR[1] + Ar[0][0]*uL[0] + Ar[0][1]*uL[1];
778:   flux[1] = Al[1][0]*uR[0] + Al[1][1]*uR[1] + Ar[1][0]*uL[0] + Ar[1][1]*uL[1];
779:   *maxspeed = c;
780:   return(0);
781: }

785: static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx)
786: {
788:   AcousticsCtx *user;
789:   PetscFList rlist = 0,rclist = 0;
790:   char rname[256] = "exact",rcname[256] = "characteristic";

793:   PetscNew(*user,&user);
794:   ctx->physics.sample         = PhysicsSample_Acoustics;
795:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
796:   ctx->physics.user           = user;
797:   ctx->physics.dof            = 2;
798:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
799:   PetscStrallocpy("v",&ctx->physics.fieldname[1]);
800:   user->c = 1;
801:   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,PETSC_NULL);
808:     PetscOptionsReal("-physics_acoustics_z","z = sqrt(bulk*rho)","",user->z,&user->z,PETSC_NULL);
809:     PetscOptionsList("-physics_acoustics_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
810:     PetscOptionsList("-physics_acoustics_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
811:   }
812:   PetscOptionsEnd();
813:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
814:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
815:   PetscFListDestroy(&rlist);
816:   PetscFListDestroy(&rclist);
817:   return(0);
818: }

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

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

826: 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*sin(2*M_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;

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

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

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

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

998:   speeds[0] = u[1]/u[0] - c;
999:   speeds[1] = u[1]/u[0] + c;
1000:   X[0*2+0] = 1;
1001:   X[0*2+1] = speeds[0];
1002:   X[1*2+0] = 1;
1003:   X[1*2+1] = speeds[1];
1004:   PetscMemcpy(Xi,X,4*sizeof(X[0]));
1005:   PetscKernel_A_gets_inverse_A_2(Xi,0);
1006:   return(0);
1007: }

1011: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
1012: {
1014:   IsoGasCtx *user;
1015:   PetscFList rlist = 0,rclist = 0;
1016:   char rname[256] = "exact",rcname[256] = "characteristic";

1019:   PetscNew(*user,&user);
1020:   ctx->physics.sample         = PhysicsSample_IsoGas;
1021:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
1022:   ctx->physics.user           = user;
1023:   ctx->physics.dof            = 2;
1024:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1025:   PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1026:   user->acoustic_speed = 1;
1027:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_IsoGas_Exact);
1028:   RiemannListAdd(&rlist,"roe",    PhysicsRiemann_IsoGas_Roe);
1029:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_IsoGas_Rusanov);
1030:   ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_IsoGas);
1031:   ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1032:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");
1033:   {
1034:     PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,PETSC_NULL);
1035:     PetscOptionsList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
1036:     PetscOptionsList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
1037:   }
1038:   PetscOptionsEnd();
1039:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
1040:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1041:   PetscFListDestroy(&rlist);
1042:   PetscFListDestroy(&rclist);
1043:   return(0);
1044: }



1048: /* --------------------------------- Shallow Water ----------------------------------- */

1050: typedef struct {
1051:   PetscReal gravity;
1052: } ShallowCtx;

1054: static inline void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
1055: {
1056:   f[0] = u[1];
1057:   f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
1058: }

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

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

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

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

1163: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
1164: {
1165:   ShallowCtx *phys = (ShallowCtx*)vctx;
1166:   PetscReal c;

1170:   c = PetscSqrtScalar(u[0]*phys->gravity);
1171:   speeds[0] = u[1]/u[0] - c;
1172:   speeds[1] = u[1]/u[0] + c;
1173:   X[0*2+0] = 1;
1174:   X[0*2+1] = speeds[0];
1175:   X[1*2+0] = 1;
1176:   X[1*2+1] = speeds[1];
1177:   PetscMemcpy(Xi,X,4*sizeof(X[0]));
1178:   PetscKernel_A_gets_inverse_A_2(Xi,0);
1179:   return(0);
1180: }

1184: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
1185: {
1187:   ShallowCtx *user;
1188:   PetscFList rlist = 0,rclist = 0;
1189:   char rname[256] = "exact",rcname[256] = "characteristic";

1192:   PetscNew(*user,&user);
1193:   /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
1194:   ctx->physics.sample         = PhysicsSample_IsoGas;
1195:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
1196:   ctx->physics.user           = user;
1197:   ctx->physics.dof            = 2;
1198:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1199:   PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1200:   user->gravity = 1;
1201:   RiemannListAdd(&rlist,"exact",  PhysicsRiemann_Shallow_Exact);
1202:   RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov);
1203:   ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Shallow);
1204:   ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1205:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");
1206:   {
1207:     PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,PETSC_NULL);
1208:     PetscOptionsList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
1209:     PetscOptionsList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
1210:   }
1211:   PetscOptionsEnd();
1212:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
1213:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1214:   PetscFListDestroy(&rlist);
1215:   PetscFListDestroy(&rclist);
1216:   return(0);
1217: }



1221: /* --------------------------------- Finite Volume Solver ----------------------------------- */

1225: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
1226: {
1227:   FVCtx          *ctx = (FVCtx*)vctx;
1228:   PetscErrorCode  ierr;
1229:   PetscInt        i,j,k,Mx,dof,xs,xm;
1230:   PetscReal       hx,cfl_idt = 0;
1231:   PetscScalar    *x,*f,*slope;
1232:   Vec             Xloc;
1233:   DM              da;

1236:   TSGetDM(ts,&da);
1237:   DMGetLocalVector(da,&Xloc);
1238:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1239:   hx = (ctx->xmax - ctx->xmin)/Mx;
1240:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1241:   DMGlobalToLocalEnd  (da,X,INSERT_VALUES,Xloc);

1243:   VecZeroEntries(F);

1245:   DMDAVecGetArray(da,Xloc,&x);
1246:   DMDAVecGetArray(da,F,&f);
1247:   DMDAGetArray(da,PETSC_TRUE,&slope);

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

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

1289:   for (i=xs; i<xs+xm+1; i++) {
1290:     PetscReal maxspeed;
1291:     PetscScalar *uL,*uR;
1292:     uL = &ctx->uLR[0];
1293:     uR = &ctx->uLR[dof];
1294:     for (j=0; j<dof; j++) {
1295:       uL[j] = x[(i-1)*dof+j] + slope[(i-1)*dof+j]*hx/2;
1296:       uR[j] = x[(i-0)*dof+j] - slope[(i-0)*dof+j]*hx/2;
1297:     }
1298:     (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed);
1299:     cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */

1301:     if (i > xs) {
1302:       for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
1303:     }
1304:     if (i < xs+xm) {
1305:       for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
1306:     }
1307:   }

1309:   DMDAVecRestoreArray(da,Xloc,&x);
1310:   DMDAVecRestoreArray(da,F,&f);
1311:   DMDARestoreArray(da,PETSC_TRUE,&slope);
1312:   DMRestoreLocalVector(da,&Xloc);

1314:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,((PetscObject)da)->comm);
1315:   if (0) {
1316:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
1317:     PetscReal dt,tnow;
1318:     TSGetTimeStep(ts,&dt);
1319:     TSGetTime(ts,&tnow);
1320:     if (dt > 0.5/ctx->cfl_idt) {
1321:       if (1) {
1322:         PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",tnow,dt,0.5/ctx->cfl_idt);
1323:       } else {
1324:         SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",dt,ctx->cfl/ctx->cfl_idt);
1325:       }
1326:     }
1327:   }
1328:   return(0);
1329: }

1333: static PetscErrorCode SmallMatMultADB(PetscScalar *C,PetscInt bs,const PetscScalar *A,const PetscReal *D,const PetscScalar *B)
1334: {
1335:   PetscInt i,j,k;

1338:   for (i=0; i<bs; i++) {
1339:     for (j=0; j<bs; j++) {
1340:       PetscScalar tmp = 0;
1341:       for (k=0; k<bs; k++) {
1342:         tmp += A[i*bs+k] * D[k] * B[k*bs+j];
1343:       }
1344:       C[i*bs+j] = tmp;
1345:     }
1346:   }
1347:   return(0);
1348: }


1353: static PetscErrorCode FVIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *vctx)
1354: {
1355:   FVCtx *ctx = (FVCtx*)vctx;
1357:   PetscInt i,j,dof = ctx->physics.dof;
1358:   PetscScalar *x,*J;
1359:   PetscReal hx;
1360:   DM da;
1361:   DMDALocalInfo dainfo;

1364:   TSGetDM(ts,&da);
1365:   DMDAVecGetArray(da,X,&x);
1366:   DMDAGetLocalInfo(da,&dainfo);
1367:   hx = (ctx->xmax - ctx->xmin)/dainfo.mx;
1368:   PetscMalloc(dof*dof*sizeof(PetscScalar),&J);
1369:   for (i=dainfo.xs; i<dainfo.xs+dainfo.xm; i++) {
1370:     (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1371:     for (j=0; j<dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]);
1372:     SmallMatMultADB(J,dof,ctx->R,ctx->speeds,ctx->Rinv);
1373:     for (j=0; j<dof*dof; j++) J[j] = J[j]/hx + shift*(j/dof == j%dof);
1374:     MatSetValuesBlocked(*B,1,&i,1,&i,J,INSERT_VALUES);
1375:   }
1376:   PetscFree(J);
1377:   DMDAVecRestoreArray(da,X,&x);

1379:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1380:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1381:   if (*A != *B) {
1382:     MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1383:     MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1384:   }
1385:   return(0);
1386: }

1390: static PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
1391: {
1393:   PetscScalar *u,*uj;
1394:   PetscInt i,j,k,dof,xs,xm,Mx;

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

1420: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
1421: {
1423:   PetscReal xmin,xmax;
1424:   PetscScalar sum,*x,tvsum,tvgsum;
1425:   PetscInt imin,imax,Mx,i,j,xs,xm,dof;
1426:   Vec Xloc;
1427:   PetscBool  iascii;

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

1447:     VecMin(X,&imin,&xmin);
1448:     VecMax(X,&imax,&xmax);
1449:     VecSum(X,&sum);
1450:     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);
1451:   } else {
1452:     SETERRQ(PETSC_COMM_SELF,1,"Viewer type not supported");
1453:   }
1454:   return(0);
1455: }

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

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

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

1494:   PetscInitialize(&argc,&argv,0,help);
1495:   comm = PETSC_COMM_WORLD;
1496:   PetscMemzero(&ctx,sizeof(ctx));

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

1518:   /* Register physical models to be available on the command line */
1519:   PetscFListAdd(&physics,"advect"          ,"",(void(*)(void))PhysicsCreate_Advect);
1520:   PetscFListAdd(&physics,"burgers"         ,"",(void(*)(void))PhysicsCreate_Burgers);
1521:   PetscFListAdd(&physics,"traffic"         ,"",(void(*)(void))PhysicsCreate_Traffic);
1522:   PetscFListAdd(&physics,"acoustics"       ,"",(void(*)(void))PhysicsCreate_Acoustics);
1523:   PetscFListAdd(&physics,"isogas"          ,"",(void(*)(void))PhysicsCreate_IsoGas);
1524:   PetscFListAdd(&physics,"shallow"         ,"",(void(*)(void))PhysicsCreate_Shallow);

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

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

1548:   /* Choose the physics from the list of registered models */
1549:   {
1550:     PetscErrorCode (*r)(FVCtx*);
1551:     PetscFListFind(physics,comm,physname,PETSC_FALSE,(void(**)(void))&r);
1552:     if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
1553:     /* Create the physics, will set the number of fields and their names */
1554:     (*r)(&ctx);
1555:   }

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

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

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

1574:   /* Create a vector to store the solution and to save the initial state */
1575:   DMCreateGlobalVector(da,&X);
1576:   VecDuplicate(X,&X0);
1577:   VecDuplicate(X,&R);

1579:   DMCreateMatrix(da,PETSC_NULL,&B);

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

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

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

1597:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);

1599:   {
1600:     PetscReal nrm1,nrmsup;
1601:     PetscInt steps;

1603:     TSSolve(ts,X,&ptime);
1604:     TSGetTimeStepNumber(ts,&steps);

1606:     PetscPrintf(comm,"Final time %8.5f, steps %d\n",ptime,steps);
1607:     if (ctx.exact) {
1608:       SolutionErrorNorms(&ctx,da,ptime,X,&nrm1,&nrmsup);
1609:       PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e  ||x-x_e||_sup %8.4e\n",nrm1,nrmsup);
1610:     }
1611:   }

1613:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);

1615:   if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
1616:   if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
1617:   if (draw & 0x4) {
1618:     Vec Y;
1619:     VecDuplicate(X,&Y);
1620:     FVSample(&ctx,da,ptime,Y);
1621:     VecAYPX(Y,-1,X);
1622:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1623:     VecDestroy(&Y);
1624:   }

1626:   if (view_final) {
1627:     PetscViewer viewer;
1628:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1629:     PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1630:     VecView(X,viewer);
1631:     PetscViewerDestroy(&viewer);
1632:   }

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