Actual source code: ex9.c

```  1: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
2:   "Solves scalar and vector problems, choose the physical model with -physics\n"
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"
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>

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 +PetscFmodReal(range+PetscFmodReal(a,range),range); }

51: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
52: typedef struct _LimitInfo {
53:   PetscReal hx;
54:   PetscInt  m;
55: } *LimitInfo;
56: static void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
57: {
58:   PetscInt i;
59:   for (i=0; i<info->m; i++) lmt[i] = 0;
60: }
61: static void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
62: {
63:   PetscInt i;
64:   for (i=0; i<info->m; i++) lmt[i] = jR[i];
65: }
66: static void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
67: {
68:   PetscInt i;
69:   for (i=0; i<info->m; i++) lmt[i] = jL[i];
70: }
71: static void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
72: {
73:   PetscInt i;
74:   for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i] + jR[i]);
75: }
76: static void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
77: {
78:   PetscInt i;
79:   for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]);
80: }
81: static void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
82: {
83:   PetscInt i;
84:   for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]));
85: }
86: static void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
87: {
88:   PetscInt i;
89:   for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]);
90: }
91: static void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
92: { /* phi = (t + abs(t)) / (1 + abs(t)) */
93:   PetscInt i;
94:   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);
95: }
96: static void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
97: { /* phi = (t + t^2) / (1 + t^2) */
98:   PetscInt i;
99:   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);
100: }
101: static void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
102: { /* phi = (t + t^2) / (1 + t^2) */
103:   PetscInt i;
104:   for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0 : (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
105: }
106: static void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
107: { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
108:   PetscInt i;
109:   for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i]) + 2*Sqr(jL[i])*jR[i])/(2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
110: }
111: static void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
112: { /* Symmetric version of above */
113:   PetscInt i;
114:   for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i])/(2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
115: }
116: static void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
117: { /* Eq 11 of Cada-Torrilhon 2009 */
118:   PetscInt i;
119:   for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]);
120: }
121: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)
122: {
123:   return PetscMax(0,PetscMin((L+2*R)/3,PetscMax(-0.5*L,PetscMin(2*L,PetscMin((L+2*R)/3,1.6*R)))));
124: }
125: static void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
126: { /* Cada-Torrilhon 2009, Eq 13 */
127:   PetscInt i;
128:   for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]);
129: }
130: static void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
131: { /* Cada-Torrilhon 2009, Eq 22 */
132:   /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
133:   const PetscReal eps = 1e-7,hx = info->hx;
134:   PetscInt        i;
135:   for (i=0; i<info->m; i++) {
136:     const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r*hx);
137:     lmt[i] = ((eta < 1-eps) ? (jL[i] + 2*jR[i]) / 3 : ((eta > 1+eps) ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]) : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3 + (1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]))));
138:   }
139: }
140: static void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
141: {
143: }
144: static void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
145: {
147: }
148: static void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
149: {
151: }
152: static void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
153: {
155: }

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

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

164: typedef struct {
165:   PetscErrorCode      (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*);
166:   RiemannFunction     riemann;
167:   ReconstructFunction characteristic;
168:   PetscErrorCode      (*destroy)(void*);
169:   void                *user;
170:   PetscInt            dof;
171:   char                *fieldname[16];
172: } PhysicsCtx;

174: typedef struct {
175:   void        (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*);
176:   PhysicsCtx  physics;
177:   MPI_Comm    comm;
178:   char        prefix[256];

180:   /* Local work arrays */
181:   PetscScalar *R,*Rinv;         /* Characteristic basis, and it's inverse.  COLUMN-MAJOR */
182:   PetscScalar *cjmpLR;          /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
183:   PetscScalar *cslope;          /* Limited slope, written in characteristic basis */
184:   PetscScalar *uLR;             /* Solution at left and right of interface, conservative variables, len=2*dof */
185:   PetscScalar *flux;            /* Flux across interface */
186:   PetscReal   *speeds;          /* Speeds of each wave */

188:   PetscReal   cfl_idt;            /* Max allowable value of 1/Delta t */
189:   PetscReal   cfl;
190:   PetscReal   xmin,xmax;
191:   PetscInt    initial;
192:   PetscBool   exact;
193:   FVBCType    bctype;
194: } FVCtx;

196: PetscErrorCode RiemannListAdd(PetscFunctionList *flist,const char *name,RiemannFunction rsolve)
197: {

202:   return(0);
203: }

205: PetscErrorCode RiemannListFind(PetscFunctionList flist,const char *name,RiemannFunction *rsolve)
206: {

210:   PetscFunctionListFind(flist,name,rsolve);
211:   if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Riemann solver \"%s\" could not be found",name);
212:   return(0);
213: }

215: PetscErrorCode ReconstructListAdd(PetscFunctionList *flist,const char *name,ReconstructFunction r)
216: {

221:   return(0);
222: }

224: PetscErrorCode ReconstructListFind(PetscFunctionList flist,const char *name,ReconstructFunction *r)
225: {

229:   PetscFunctionListFind(flist,name,r);
230:   if (!*r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Reconstruction \"%s\" could not be found",name);
231:   return(0);
232: }

234: /* --------------------------------- Physics ----------------------------------- */
235: /*
236:   Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction.  These
237:   are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the
238:   number of fields and their names, and a function to deallocate private storage.
239: */

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

247:   for (i=0; i<m; i++) {
248:     for (j=0; j<m; j++) Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j);
249:     speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
250:   }
251:   return(0);
252: }

254: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
255: {

259:   PetscFree(vctx);
260:   return(0);
261: }

263: /* --------------------------------- Advection ----------------------------------- */

265: typedef struct {
266:   PetscReal a;                  /* advective velocity */

269: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
270: {
272:   PetscReal speed;

275:   speed     = ctx->a;
276:   flux[0]   = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
277:   *maxspeed = speed;
278:   return(0);
279: }

281: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
282: {

286:   X[0]      = 1.;
287:   Xi[0]     = 1.;
288:   speeds[0] = ctx->a;
289:   return(0);
290: }

292: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
293: {
295:   PetscReal a    = ctx->a,x0;

298:   switch (bctype) {
299:     case FVBC_OUTFLOW:   x0 = x-a*t; break;
300:     case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
301:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
302:   }
303:   switch (initial) {
304:     case 0: u[0] = (x0 < 0) ? 1 : -1; break;
305:     case 1: u[0] = (x0 < 0) ? -1 : 1; break;
306:     case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
307:     case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
308:     case 4: u[0] = PetscAbs(x0); break;
309:     case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
310:     case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
311:     case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
312:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
313:   }
314:   return(0);
315: }

318: {

323:   PetscNew(&user);
327:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
328:   ctx->physics.user           = user;
329:   ctx->physics.dof            = 1;
330:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
331:   user->a = 1;
333:   {
335:   }
336:   PetscOptionsEnd();
337:   return(0);
338: }

340: /* --------------------------------- Burgers ----------------------------------- */

342: typedef struct {
343:   PetscReal lxf_speed;
344: } BurgersCtx;

346: static PetscErrorCode PhysicsSample_Burgers(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
347: {
349:   if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exact solution not implemented for periodic");
350:   switch (initial) {
351:     case 0: u[0] = (x < 0) ? 1 : -1; break;
352:     case 1:
353:       if       (x < -t) u[0] = -1;
354:       else if  (x < t)  u[0] = x/t;
355:       else              u[0] = 1;
356:       break;
357:     case 2:
358:       if      (x < 0)       u[0] = 0;
359:       else if (x <= t)      u[0] = x/t;
360:       else if (x < 1+0.5*t) u[0] = 1;
361:       else                  u[0] = 0;
362:       break;
363:     case 3:
364:       if       (x < 0.2*t) u[0] = 0.2;
365:       else if  (x < t) u[0] = x/t;
366:       else             u[0] = 1;
367:       break;
368:     case 4:
369:       if (t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only initial condition available");
370:       u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin)));
371:       break;
372:     case 5:                     /* Pure shock solution */
373:       if (x < 0.5*t) u[0] = 1;
374:       else u[0] = 0;
375:       break;
376:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
377:   }
378:   return(0);
379: }

381: static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
382: {
384:   if (uL[0] < uR[0]) {          /* rarefaction */
385:     flux[0] = (uL[0]*uR[0] < 0)
386:       ? 0                       /* sonic rarefaction */
387:       : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0]));
388:   } else {                      /* shock */
389:     flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0]));
390:   }
391:   *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
392:   return(0);
393: }

395: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
396: {
397:   PetscReal speed;

400:   speed   = 0.5*(uL[0] + uR[0]);
401:   flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
402:   if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
403:   *maxspeed = speed;
404:   return(0);
405: }

407: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
408: {
409:   PetscReal   c;
410:   PetscScalar fL,fR;

413:   c         = ((BurgersCtx*)vctx)->lxf_speed;
414:   fL        = 0.5*PetscSqr(uL[0]);
415:   fR        = 0.5*PetscSqr(uR[0]);
416:   flux[0]   = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
417:   *maxspeed = c;
418:   return(0);
419: }

421: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
422: {
423:   PetscReal   c;
424:   PetscScalar fL,fR;

427:   c         = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0]));
428:   fL        = 0.5*PetscSqr(uL[0]);
429:   fR        = 0.5*PetscSqr(uR[0]);
430:   flux[0]   = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
431:   *maxspeed = c;
432:   return(0);
433: }

435: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
436: {
437:   BurgersCtx        *user;
438:   PetscErrorCode    ierr;
439:   RiemannFunction   r;
440:   PetscFunctionList rlist      = 0;
441:   char              rname[256] = "exact";

444:   PetscNew(&user);

446:   ctx->physics.sample         = PhysicsSample_Burgers;
447:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
448:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
449:   ctx->physics.user           = user;
450:   ctx->physics.dof            = 1;

452:   PetscStrallocpy("u",&ctx->physics.fieldname[0]);
458:   {
459:     PetscOptionsFList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
460:   }
461:   PetscOptionsEnd();
462:   RiemannListFind(rlist,rname,&r);
463:   PetscFunctionListDestroy(&rlist);
464:   ctx->physics.riemann = r;

466:   /* *
467:   * Hack to deal with LxF in semi-discrete form
468:   * max speed is 1 for the basic initial conditions (where |u| <= 1)
469:   * */
470:   if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
471:   return(0);
472: }

474: /* --------------------------------- Traffic ----------------------------------- */

476: typedef struct {
477:   PetscReal lxf_speed;
478:   PetscReal a;
479: } TrafficCtx;

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

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

488:   if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exact solution not implemented for periodic");
489:   switch (initial) {
490:     case 0:
491:       u[0] = (-a*t < x) ? 2 : 0; break;
492:     case 1:
493:       if      (x < PetscMin(2*a*t,0.5+a*t)) u[0] = -1;
494:       else if (x < 1)                       u[0] = 0;
495:       else                                  u[0] = 1;
496:       break;
497:     case 2:
498:       if (t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only initial condition available");
499:       u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin)));
500:       break;
501:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
502:   }
503:   return(0);
504: }

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

511:   if (uL[0] < uR[0]) {
512:     flux[0] = PetscMin(TrafficFlux(a,uL[0]),TrafficFlux(a,uR[0]));
513:   } else {
514:     flux[0] = (uR[0] < 0.5 && 0.5 < uL[0]) ? TrafficFlux(a,0.5) : PetscMax(TrafficFlux(a,uL[0]),TrafficFlux(a,uR[0]));
515:   }
516:   *maxspeed = a*MaxAbs(1-2*uL[0],1-2*uR[0]);
517:   return(0);
518: }

520: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
521: {
522:   PetscReal a = ((TrafficCtx*)vctx)->a;
523:   PetscReal speed;

526:   speed = a*(1 - (uL[0] + uR[0]));
527:   flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
528:   *maxspeed = speed;
529:   return(0);
530: }

532: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
533: {
534:   TrafficCtx *phys = (TrafficCtx*)vctx;
535:   PetscReal  a     = phys->a;
536:   PetscReal  speed;

539:   speed     = a*(1 - (uL[0] + uR[0]));
540:   flux[0]   = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]);
541:   *maxspeed = speed;
542:   return(0);
543: }

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

551:   speed     = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0]));
552:   flux[0]   = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]);
553:   *maxspeed = speed;
554:   return(0);
555: }

557: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
558: {
559:   PetscErrorCode    ierr;
560:   TrafficCtx        *user;
561:   RiemannFunction   r;
562:   PetscFunctionList rlist      = 0;
563:   char              rname[256] = "exact";

566:   PetscNew(&user);
567:   ctx->physics.sample         = PhysicsSample_Traffic;
568:   ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
569:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
570:   ctx->physics.user           = user;
571:   ctx->physics.dof            = 1;

573:   PetscStrallocpy("density",&ctx->physics.fieldname[0]);
574:   user->a = 0.5;
579:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");
580:     PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,NULL);
581:     PetscOptionsFList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
582:   PetscOptionsEnd();

584:   RiemannListFind(rlist,rname,&r);
585:   PetscFunctionListDestroy(&rlist);

587:   ctx->physics.riemann = r;

589:   /* *
590:   * Hack to deal with LxF in semi-discrete form
591:   * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
592:   * */
593:   if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3*user->a;
594:   return(0);
595: }

597: /* --------------------------------- Linear Acoustics ----------------------------------- */

599: /* Flux: u_t + (A u)_x
600:  * z = sqrt(rho*bulk), c = sqrt(rho/bulk)
601:  * Spectral decomposition: A = R * D * Rinv
602:  * [    cz] = [-z   z] [-c    ] [-1/2z  1/2]
603:  * [c/z   ] = [ 1   1] [     c] [ 1/2z  1/2]
604:  *
605:  * We decompose this into the left-traveling waves Al = R * D^- Rinv
606:  * and the right-traveling waves Ar = R * D^+ * Rinv
607:  * Multiplying out these expressions produces the following two matrices
608:  */

610: typedef struct {
611:   PetscReal c;                  /* speed of sound: c = sqrt(bulk/rho) */
612:   PetscReal z;                  /* impedence: z = sqrt(rho*bulk) */
613: } AcousticsCtx;

615: PETSC_UNUSED PETSC_STATIC_INLINE void AcousticsFlux(AcousticsCtx *ctx,const PetscScalar *u,PetscScalar *f)
616: {
617:   f[0] = ctx->c*ctx->z*u[1];
618:   f[1] = ctx->c/ctx->z*u[0];
619: }

621: static PetscErrorCode PhysicsCharacteristic_Acoustics(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
622: {
623:   AcousticsCtx *phys = (AcousticsCtx*)vctx;
624:   PetscReal    z     = phys->z,c = phys->c;

627:   X[0*2+0]  = -z;
628:   X[0*2+1]  = z;
629:   X[1*2+0]  = 1;
630:   X[1*2+1]  = 1;
631:   Xi[0*2+0] = -1./(2*z);
632:   Xi[0*2+1] = 1./2;
633:   Xi[1*2+0] = 1./(2*z);
634:   Xi[1*2+1] = 1./2;
635:   speeds[0] = -c;
636:   speeds[1] = c;
637:   return(0);
638: }

640: static PetscErrorCode PhysicsSample_Acoustics_Initial(AcousticsCtx *phys,PetscInt initial,PetscReal xmin,PetscReal xmax,PetscReal x,PetscReal *u)
641: {
643:   switch (initial) {
644:   case 0:
645:     u[0] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.2) < 0.1) ? 1 : 0.5;
646:     u[1] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.7) < 0.1) ? 1 : -0.5;
647:     break;
648:   case 1:
649:     u[0] = PetscCosReal(3 * 2*PETSC_PI*x/(xmax-xmin));
650:     u[1] = PetscExpReal(-PetscSqr(x - (xmax + xmin)/2) / (2*PetscSqr(0.2*(xmax - xmin)))) - 0.5;
651:     break;
652:   default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
653:   }
654:   return(0);
655: }

657: static PetscErrorCode PhysicsSample_Acoustics(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
658: {
659:   AcousticsCtx   *phys = (AcousticsCtx*)vctx;
660:   PetscReal      c     = phys->c;
661:   PetscReal      x0a,x0b,u0a[2],u0b[2],tmp[2];
662:   PetscReal      X[2][2],Xi[2][2],dummy[2];

666:   switch (bctype) {
667:   case FVBC_OUTFLOW:
668:     x0a = x+c*t;
669:     x0b = x-c*t;
670:     break;
671:   case FVBC_PERIODIC:
672:     x0a = RangeMod(x+c*t,xmin,xmax);
673:     x0b = RangeMod(x-c*t,xmin,xmax);
674:     break;
675:   default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
676:   }
677:   PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0a,u0a);
678:   PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0b,u0b);
679:   PhysicsCharacteristic_Acoustics(vctx,2,u,&X[0][0],&Xi[0][0],dummy);
680:   tmp[0] = Xi[0][0]*u0a[0] + Xi[0][1]*u0a[1];
681:   tmp[1] = Xi[1][0]*u0b[0] + Xi[1][1]*u0b[1];
682:   u[0]   = X[0][0]*tmp[0] + X[0][1]*tmp[1];
683:   u[1]   = X[1][0]*tmp[0] + X[1][1]*tmp[1];
684:   return(0);
685: }

687: static PetscErrorCode PhysicsRiemann_Acoustics_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
688: {
689:   AcousticsCtx *phys = (AcousticsCtx*)vctx;
690:   PetscReal    c     = phys->c,z = phys->z;
691:   PetscReal
692:     Al[2][2] = {{-c/2     , c*z/2  },
693:                 {c/(2*z)  , -c/2   }}, /* Left traveling waves */
694:     Ar[2][2] = {{c/2      , c*z/2  },
695:                 {c/(2*z)  , c/2    }}; /* Right traveling waves */

698:   flux[0]   = Al[0][0]*uR[0] + Al[0][1]*uR[1] + Ar[0][0]*uL[0] + Ar[0][1]*uL[1];
699:   flux[1]   = Al[1][0]*uR[0] + Al[1][1]*uR[1] + Ar[1][0]*uL[0] + Ar[1][1]*uL[1];
700:   *maxspeed = c;
701:   return(0);
702: }

704: static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx)
705: {
706:   PetscErrorCode    ierr;
707:   AcousticsCtx      *user;
708:   PetscFunctionList rlist      = 0,rclist = 0;
709:   char              rname[256] = "exact",rcname[256] = "characteristic";

712:   PetscNew(&user);
713:   ctx->physics.sample         = PhysicsSample_Acoustics;
714:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
715:   ctx->physics.user           = user;
716:   ctx->physics.dof            = 2;

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

721:   user->c = 1;
722:   user->z = 1;

727:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for linear Acoustics","");
728:   {
729:     PetscOptionsReal("-physics_acoustics_c","c = sqrt(bulk/rho)","",user->c,&user->c,NULL);
730:     PetscOptionsReal("-physics_acoustics_z","z = sqrt(bulk*rho)","",user->z,&user->z,NULL);
731:     PetscOptionsFList("-physics_acoustics_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
732:     PetscOptionsFList("-physics_acoustics_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
733:   }
734:   PetscOptionsEnd();
735:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
736:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
737:   PetscFunctionListDestroy(&rlist);
738:   PetscFunctionListDestroy(&rclist);
739:   return(0);
740: }

742: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */

744: typedef struct {
745:   PetscReal acoustic_speed;
746: } IsoGasCtx;

748: PETSC_STATIC_INLINE void IsoGasFlux(PetscReal c,const PetscScalar *u,PetscScalar *f)
749: {
750:   f[0] = u[1];
751:   f[1] = PetscSqr(u[1])/u[0] + c*c*u[0];
752: }

754: static PetscErrorCode PhysicsSample_IsoGas(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
755: {
757:   if (t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exact solutions not implemented for t > 0");
758:   switch (initial) {
759:     case 0:
760:       u[0] = (x < 0) ? 1 : 0.5;
761:       u[1] = (x < 0) ? 1 : 0.7;
762:       break;
763:     case 1:
764:       u[0] = 1+0.5*PetscSinReal(2*PETSC_PI*x);
765:       u[1] = 1*u[0];
766:       break;
767:     default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
768:   }
769:   return(0);
770: }

772: static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
773: {
774:   IsoGasCtx   *phys = (IsoGasCtx*)vctx;
775:   PetscReal   c     = phys->acoustic_speed;
776:   PetscScalar ubar,du[2],a[2],fL[2],fR[2],lam[2],ustar[2],R[2][2];
777:   PetscInt    i;

780:   ubar = (uL[1]/PetscSqrtScalar(uL[0]) + uR[1]/PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0]));
781:   /* write fluxuations in characteristic basis */
782:   du[0] = uR[0] - uL[0];
783:   du[1] = uR[1] - uL[1];
784:   a[0]  = (1/(2*c)) * ((ubar + c)*du[0] - du[1]);
785:   a[1]  = (1/(2*c)) * ((-ubar + c)*du[0] + du[1]);
786:   /* wave speeds */
787:   lam[0] = ubar - c;
788:   lam[1] = ubar + c;
789:   /* Right eigenvectors */
790:   R[0][0] = 1; R[0][1] = ubar-c;
791:   R[1][0] = 1; R[1][1] = ubar+c;
792:   /* Compute state in star region (between the 1-wave and 2-wave) */
793:   for (i=0; i<2; i++) ustar[i] = uL[i] + a[0]*R[0][i];
794:   if (uL[1]/uL[0] < c && c < ustar[1]/ustar[0]) { /* 1-wave is sonic rarefaction */
795:     PetscScalar ufan[2];
796:     ufan[0] = uL[0]*PetscExpScalar(uL[1]/(uL[0]*c) - 1);
797:     ufan[1] = c*ufan[0];
798:     IsoGasFlux(c,ufan,flux);
799:   } else if (ustar[1]/ustar[0] < -c && -c < uR[1]/uR[0]) { /* 2-wave is sonic rarefaction */
800:     PetscScalar ufan[2];
801:     ufan[0] = uR[0]*PetscExpScalar(-uR[1]/(uR[0]*c) - 1);
802:     ufan[1] = -c*ufan[0];
803:     IsoGasFlux(c,ufan,flux);
804:   } else {                      /* Centered form */
805:     IsoGasFlux(c,uL,fL);
806:     IsoGasFlux(c,uR,fR);
807:     for (i=0; i<2; i++) {
808:       PetscScalar absdu = PetscAbsScalar(lam[0])*a[0]*R[0][i] + PetscAbsScalar(lam[1])*a[1]*R[1][i];
809:       flux[i] = 0.5*(fL[i]+fR[i]) - 0.5*absdu;
810:     }
811:   }
812:   *maxspeed = MaxAbs(lam[0],lam[1]);
813:   return(0);
814: }

816: static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
817: {
818:   IsoGasCtx                   *phys = (IsoGasCtx*)vctx;
819:   PetscReal                   c     = phys->acoustic_speed;
820:   PetscScalar                 ustar[2];
821:   struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
822:   PetscInt                    i;

825:   if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed density is negative");
826:   {
827:     /* Solve for star state */
828:     PetscScalar res,tmp,rho = 0.5*(L.rho + R.rho); /* initial guess */
829:     for (i=0; i<20; i++) {
830:       PetscScalar fr,fl,dfr,dfl;
831:       fl = (L.rho < rho)
832:         ? (rho-L.rho)/PetscSqrtScalar(L.rho*rho)       /* shock */
833:         : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
834:       fr = (R.rho < rho)
835:         ? (rho-R.rho)/PetscSqrtScalar(R.rho*rho)       /* shock */
836:         : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
837:       res = R.u-L.u + c*(fr+fl);
838:       if (PetscIsInfOrNanScalar(res)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinity or Not-a-Number generated in computation");
839:       if (PetscAbsScalar(res) < 1e-10) {
840:         star.rho = rho;
841:         star.u   = L.u - c*fl;
842:         goto converged;
843:       }
844:       dfl = (L.rho < rho) ? 1/PetscSqrtScalar(L.rho*rho)*(1 - 0.5*(rho-L.rho)/rho) : 1/rho;
845:       dfr = (R.rho < rho) ? 1/PetscSqrtScalar(R.rho*rho)*(1 - 0.5*(rho-R.rho)/rho) : 1/rho;
846:       tmp = rho - res/(c*(dfr+dfl));
847:       if (tmp <= 0) rho /= 2;   /* Guard against Newton shooting off to a negative density */
848:       else rho = tmp;
849:       if (!((rho > 0) && PetscIsNormalScalar(rho))) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP,"non-normal iterate rho=%g",(double)PetscRealPart(rho));
850:     }
851:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Newton iteration for star.rho diverged after %D iterations",i);
852:   }
853: converged:
854:   if (L.u-c < 0 && 0 < star.u-c) { /* 1-wave is sonic rarefaction */
855:     PetscScalar ufan[2];
856:     ufan[0] = L.rho*PetscExpScalar(L.u/c - 1);
857:     ufan[1] = c*ufan[0];
858:     IsoGasFlux(c,ufan,flux);
859:   } else if (star.u+c < 0 && 0 < R.u+c) { /* 2-wave is sonic rarefaction */
860:     PetscScalar ufan[2];
861:     ufan[0] = R.rho*PetscExpScalar(-R.u/c - 1);
862:     ufan[1] = -c*ufan[0];
863:     IsoGasFlux(c,ufan,flux);
864:   } else if ((L.rho >= star.rho && L.u-c >= 0) || (L.rho < star.rho && (star.rho*star.u-L.rho*L.u)/(star.rho-L.rho) > 0)) {
865:     /* 1-wave is supersonic rarefaction, or supersonic shock */
866:     IsoGasFlux(c,uL,flux);
867:   } else if ((star.rho <= R.rho && R.u+c <= 0) || (star.rho > R.rho && (R.rho*R.u-star.rho*star.u)/(R.rho-star.rho) < 0)) {
868:     /* 2-wave is supersonic rarefaction or supersonic shock */
869:     IsoGasFlux(c,uR,flux);
870:   } else {
871:     ustar[0] = star.rho;
872:     ustar[1] = star.rho*star.u;
873:     IsoGasFlux(c,ustar,flux);
874:   }
875:   *maxspeed = MaxAbs(MaxAbs(star.u-c,star.u+c),MaxAbs(L.u-c,R.u+c));
876:   return(0);
877: }

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

886:   if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed density is negative");
887:   IsoGasFlux(c,uL,fL);
888:   IsoGasFlux(c,uR,fR);
889:   s         = PetscMax(PetscAbs(L.u),PetscAbs(R.u))+c;
890:   flux[0]   = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
891:   flux[1]   = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
892:   *maxspeed = s;
893:   return(0);
894: }

896: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
897: {
898:   IsoGasCtx      *phys = (IsoGasCtx*)vctx;
899:   PetscReal      c     = phys->acoustic_speed;

903:   speeds[0] = u[1]/u[0] - c;
904:   speeds[1] = u[1]/u[0] + c;
905:   X[0*2+0]  = 1;
906:   X[0*2+1]  = speeds[0];
907:   X[1*2+0]  = 1;
908:   X[1*2+1]  = speeds[1];
909:   PetscArraycpy(Xi,X,4);
910:   PetscKernel_A_gets_inverse_A_2(Xi,0,PETSC_FALSE,NULL);
911:   return(0);
912: }

914: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
915: {
916:   PetscErrorCode    ierr;
917:   IsoGasCtx         *user;
918:   PetscFunctionList rlist = 0,rclist = 0;
919:   char              rname[256] = "exact",rcname[256] = "characteristic";

922:   PetscNew(&user);
923:   ctx->physics.sample         = PhysicsSample_IsoGas;
924:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
925:   ctx->physics.user           = user;
926:   ctx->physics.dof            = 2;

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

931:   user->acoustic_speed = 1;

938:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");
939:     PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,NULL);
940:     PetscOptionsFList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
941:     PetscOptionsFList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
942:   PetscOptionsEnd();
943:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
944:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
945:   PetscFunctionListDestroy(&rlist);
946:   PetscFunctionListDestroy(&rclist);
947:   return(0);
948: }

950: /* --------------------------------- Shallow Water ----------------------------------- */
951: typedef struct {
952:   PetscReal gravity;
953: } ShallowCtx;

955: PETSC_STATIC_INLINE void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
956: {
957:   f[0] = u[1];
958:   f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
959: }

961: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
962: {
963:   ShallowCtx                *phys = (ShallowCtx*)vctx;
964:   PetscScalar               g    = phys->gravity,ustar[2],cL,cR,c,cstar;
965:   struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
966:   PetscInt                  i;

969:   if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative");
970:   cL = PetscSqrtScalar(g*L.h);
971:   cR = PetscSqrtScalar(g*R.h);
972:   c  = PetscMax(cL,cR);
973:   {
974:     /* Solve for star state */
975:     const PetscInt maxits = 50;
976:     PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */
977:     h0 = h;
978:     for (i=0; i<maxits; i++) {
979:       PetscScalar fr,fl,dfr,dfl;
980:       fl = (L.h < h)
981:         ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */
982:         : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h);   /* rarefaction */
983:       fr = (R.h < h)
984:         ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */
985:         : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h);   /* rarefaction */
986:       res = R.u - L.u + fr + fl;
987:       if (PetscIsInfOrNanScalar(res)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinity or Not-a-Number generated in computation");
988:       if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h-h0) < 1e-8)) {
989:         star.h = h;
990:         star.u = L.u - fl;
991:         goto converged;
992:       } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) {        /* Line search */
993:         h = 0.8*h0 + 0.2*h;
994:         continue;
995:       }
996:       /* Accept the last step and take another */
997:       res0 = res;
998:       h0 = h;
999:       dfl = (L.h < h) ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h) : PetscSqrtScalar(g/h);
1000:       dfr = (R.h < h) ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h) : PetscSqrtScalar(g/h);
1001:       tmp = h - res/(dfr+dfl);
1002:       if (tmp <= 0) h /= 2;   /* Guard against Newton shooting off to a negative thickness */
1003:       else h = tmp;
1004:       if (!((h > 0) && PetscIsNormalScalar(h))) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP,"non-normal iterate h=%g",(double)h);
1005:     }
1006:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Newton iteration for star.h diverged after %D iterations",i);
1007:   }
1008: converged:
1009:   cstar = PetscSqrtScalar(g*star.h);
1010:   if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */
1011:     PetscScalar ufan[2];
1012:     ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL);
1013:     ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0];
1014:     ShallowFlux(phys,ufan,flux);
1015:   } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */
1016:     PetscScalar ufan[2];
1017:     ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR);
1018:     ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0];
1019:     ShallowFlux(phys,ufan,flux);
1020:   } else if ((L.h >= star.h && L.u-c >= 0) || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) {
1021:     /* 1-wave is right-travelling shock (supersonic) */
1022:     ShallowFlux(phys,uL,flux);
1023:   } else if ((star.h <= R.h && R.u+c <= 0) || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) {
1024:     /* 2-wave is left-travelling shock (supersonic) */
1025:     ShallowFlux(phys,uR,flux);
1026:   } else {
1027:     ustar[0] = star.h;
1028:     ustar[1] = star.h*star.u;
1029:     ShallowFlux(phys,ustar,flux);
1030:   }
1031:   *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR));
1032:   return(0);
1033: }

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

1042:   if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative");
1043:   ShallowFlux(phys,uL,fL);
1044:   ShallowFlux(phys,uR,fR);
1045:   s         = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h));
1046:   flux[0]   = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
1047:   flux[1]   = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
1048:   *maxspeed = s;
1049:   return(0);
1050: }

1052: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
1053: {
1054:   ShallowCtx     *phys = (ShallowCtx*)vctx;
1055:   PetscReal      c;

1059:   c         = PetscSqrtScalar(u[0]*phys->gravity);
1060:   speeds[0] = u[1]/u[0] - c;
1061:   speeds[1] = u[1]/u[0] + c;
1062:   X[0*2+0]  = 1;
1063:   X[0*2+1]  = speeds[0];
1064:   X[1*2+0]  = 1;
1065:   X[1*2+1]  = speeds[1];
1066:   PetscArraycpy(Xi,X,4);
1067:   PetscKernel_A_gets_inverse_A_2(Xi,0,PETSC_FALSE,NULL);
1068:   return(0);
1069: }

1071: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
1072: {
1073:   PetscErrorCode    ierr;
1074:   ShallowCtx        *user;
1075:   PetscFunctionList rlist = 0,rclist = 0;
1076:   char              rname[256] = "exact",rcname[256] = "characteristic";

1079:   PetscNew(&user);
1080:   /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
1081:   ctx->physics.sample         = PhysicsSample_IsoGas;
1082:   ctx->physics.destroy        = PhysicsDestroy_SimpleFree;
1083:   ctx->physics.user           = user;
1084:   ctx->physics.dof            = 2;

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

1089:   user->gravity = 1;

1095:   PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");
1096:     PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,NULL);
1097:     PetscOptionsFList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
1098:     PetscOptionsFList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
1099:   PetscOptionsEnd();
1100:   RiemannListFind(rlist,rname,&ctx->physics.riemann);
1101:   ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1102:   PetscFunctionListDestroy(&rlist);
1103:   PetscFunctionListDestroy(&rclist);
1104:   return(0);
1105: }

1107: /* --------------------------------- Finite Volume Solver ----------------------------------- */

1109: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
1110: {
1111:   FVCtx             *ctx = (FVCtx*)vctx;
1112:   PetscErrorCode    ierr;
1113:   PetscInt          i,j,k,Mx,dof,xs,xm;
1114:   PetscReal         hx,cfl_idt = 0;
1115:   PetscScalar       *x,*f,*slope;
1116:   Vec               Xloc;
1117:   DM                da;

1120:   TSGetDM(ts,&da);
1121:   DMGetLocalVector(da,&Xloc);
1122:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1123:   hx   = (ctx->xmax - ctx->xmin)/Mx;
1124:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1125:   DMGlobalToLocalEnd  (da,X,INSERT_VALUES,Xloc);

1127:   VecZeroEntries(F);

1129:   DMDAVecGetArray(da,Xloc,&x);
1130:   DMDAVecGetArray(da,F,&f);
1131:   DMDAGetArray(da,PETSC_TRUE,&slope);

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

1135:   if (ctx->bctype == FVBC_OUTFLOW) {
1136:     for (i=xs-2; i<0; i++) {
1137:       for (j=0; j<dof; j++) x[i*dof+j] = x[j];
1138:     }
1139:     for (i=Mx; i<xs+xm+2; i++) {
1140:       for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
1141:     }
1142:   }
1143:   for (i=xs-1; i<xs+xm+1; i++) {
1144:     struct _LimitInfo info;
1145:     PetscScalar       *cjmpL,*cjmpR;
1146:     /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
1147:     (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1148:     /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
1149:     PetscArrayzero(ctx->cjmpLR,2*dof);
1150:     cjmpL = &ctx->cjmpLR[0];
1151:     cjmpR = &ctx->cjmpLR[dof];
1152:     for (j=0; j<dof; j++) {
1153:       PetscScalar jmpL,jmpR;
1154:       jmpL = x[(i+0)*dof+j] - x[(i-1)*dof+j];
1155:       jmpR = x[(i+1)*dof+j] - x[(i+0)*dof+j];
1156:       for (k=0; k<dof; k++) {
1157:         cjmpL[k] += ctx->Rinv[k+j*dof] * jmpL;
1158:         cjmpR[k] += ctx->Rinv[k+j*dof] * jmpR;
1159:       }
1160:     }
1161:     /* Apply limiter to the left and right characteristic jumps */
1162:     info.m  = dof;
1163:     info.hx = hx;
1164:     (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
1165:     for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
1166:     for (j=0; j<dof; j++) {
1167:       PetscScalar tmp = 0;
1168:       for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof] * ctx->cslope[k];
1169:       slope[i*dof+j] = tmp;
1170:     }
1171:   }

1173:   for (i=xs; i<xs+xm+1; i++) {
1174:     PetscReal   maxspeed;
1175:     PetscScalar *uL,*uR;
1176:     uL = &ctx->uLR[0];
1177:     uR = &ctx->uLR[dof];
1178:     for (j=0; j<dof; j++) {
1179:       uL[j] = x[(i-1)*dof+j] + slope[(i-1)*dof+j]*hx/2;
1180:       uR[j] = x[(i-0)*dof+j] - slope[(i-0)*dof+j]*hx/2;
1181:     }
1182:     (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed);
1183:     cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */

1185:     if (i > xs) {
1186:       for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
1187:     }
1188:     if (i < xs+xm) {
1189:       for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
1190:     }
1191:   }

1193:   DMDAVecRestoreArray(da,Xloc,&x);
1194:   DMDAVecRestoreArray(da,F,&f);
1195:   DMDARestoreArray(da,PETSC_TRUE,&slope);
1196:   DMRestoreLocalVector(da,&Xloc);

1198:   MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));
1199:   if (0) {
1200:     /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
1201:     PetscReal dt,tnow;
1202:     TSGetTimeStep(ts,&dt);
1203:     TSGetTime(ts,&tnow);
1204:     if (dt > 0.5/ctx->cfl_idt) {
1205:       PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
1206:     }
1207:   }
1208:   return(0);
1209: }

1211: static PetscErrorCode SmallMatMultADB(PetscScalar *C,PetscInt bs,const PetscScalar *A,const PetscReal *D,const PetscScalar *B)
1212: {
1213:   PetscInt i,j,k;

1216:   for (i=0; i<bs; i++) {
1217:     for (j=0; j<bs; j++) {
1218:       PetscScalar tmp = 0;
1219:       for (k=0; k<bs; k++) tmp += A[i*bs+k] * D[k] * B[k*bs+j];
1220:       C[i*bs+j] = tmp;
1221:     }
1222:   }
1223:   return(0);
1224: }

1226: static PetscErrorCode FVIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void *vctx)
1227: {
1228:   FVCtx             *ctx = (FVCtx*)vctx;
1229:   PetscErrorCode    ierr;
1230:   PetscInt          i,j,dof = ctx->physics.dof;
1231:   PetscScalar       *J;
1232:   const PetscScalar *x;
1233:   PetscReal         hx;
1234:   DM                da;
1235:   DMDALocalInfo     dainfo;

1238:   TSGetDM(ts,&da);
1240:   DMDAGetLocalInfo(da,&dainfo);
1241:   hx   = (ctx->xmax - ctx->xmin)/dainfo.mx;
1242:   PetscMalloc1(dof*dof,&J);
1243:   for (i=dainfo.xs; i<dainfo.xs+dainfo.xm; i++) {
1244:     (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1245:     for (j=0; j<dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]);
1247:     for (j=0; j<dof*dof; j++) J[j] = J[j]/hx + shift*(j/dof == j%dof);
1248:     MatSetValuesBlocked(B,1,&i,1,&i,J,INSERT_VALUES);
1249:   }
1250:   PetscFree(J);

1253:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1254:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1255:   if (A != B) {
1256:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1257:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1258:   }
1259:   return(0);
1260: }

1262: static PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
1263: {
1265:   PetscScalar    *u,*uj;
1266:   PetscInt       i,j,k,dof,xs,xm,Mx;

1269:   if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function");
1270:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1271:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1272:   DMDAVecGetArray(da,U,&u);
1273:   PetscMalloc1(dof,&uj);
1274:   for (i=xs; i<xs+xm; i++) {
1275:     const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h;
1276:     const PetscInt  N = 200;
1277:     /* Integrate over cell i using trapezoid rule with N points. */
1278:     for (k=0; k<dof; k++) u[i*dof+k] = 0;
1279:     for (j=0; j<N+1; j++) {
1280:       PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N;
1281:       (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
1282:       for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
1283:     }
1284:   }
1285:   DMDAVecRestoreArray(da,U,&u);
1286:   PetscFree(uj);
1287:   return(0);
1288: }

1290: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
1291: {
1292:   PetscErrorCode    ierr;
1293:   PetscReal         xmin,xmax;
1294:   PetscScalar       sum,tvsum,tvgsum;
1295:   const PetscScalar *x;
1296:   PetscInt          imin,imax,Mx,i,j,xs,xm,dof;
1297:   Vec               Xloc;
1298:   PetscBool         iascii;

1301:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1302:   if (iascii) {
1303:     /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
1304:     DMGetLocalVector(da,&Xloc);
1305:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1306:     DMGlobalToLocalEnd  (da,X,INSERT_VALUES,Xloc);
1308:     DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1309:     DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1310:     tvsum = 0;
1311:     for (i=xs; i<xs+xm; i++) {
1312:       for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j] - x[(i-1)*dof+j]);
1313:     }
1314:     MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)da));
1316:     DMRestoreLocalVector(da,&Xloc);

1318:     VecMin(X,&imin,&xmin);
1319:     VecMax(X,&imax,&xmax);
1320:     VecSum(X,&sum);
1321:     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));
1322:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported");
1323:   return(0);
1324: }

1326: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1,PetscReal *nrmsup)
1327: {
1329:   Vec            Y;
1330:   PetscInt       Mx;

1333:   VecGetSize(X,&Mx);
1334:   VecDuplicate(X,&Y);
1335:   FVSample(ctx,da,t,Y);
1336:   VecAYPX(Y,-1,X);
1337:   VecNorm(Y,NORM_1,nrm1);
1338:   VecNorm(Y,NORM_INFINITY,nrmsup);
1339:   *nrm1 /= Mx;
1340:   VecDestroy(&Y);
1341:   return(0);
1342: }

1344: int main(int argc,char *argv[])
1345: {
1346:   char              lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
1347:   PetscFunctionList limiters   = 0,physics = 0;
1348:   MPI_Comm          comm;
1349:   TS                ts;
1350:   DM                da;
1351:   Vec               X,X0,R;
1352:   Mat               B;
1353:   FVCtx             ctx;
1354:   PetscInt          i,dof,xs,xm,Mx,draw = 0;
1355:   PetscBool         view_final = PETSC_FALSE;
1356:   PetscReal         ptime;
1357:   PetscErrorCode    ierr;

1359:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
1360:   comm = PETSC_COMM_WORLD;
1361:   PetscMemzero(&ctx,sizeof(ctx));

1363:   /* Register limiters to be available on the command line */

1383:   /* Register physical models to be available on the command line */

1391:   ctx.comm = comm;
1392:   ctx.cfl  = 0.9; ctx.bctype = FVBC_PERIODIC;
1393:   ctx.xmin = -1; ctx.xmax = 1;
1394:   PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
1395:     PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
1396:     PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
1397:     PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);
1398:     PetscOptionsFList("-physics","Name of physics (Riemann solver and characteristics) to use","",physics,physname,physname,sizeof(physname),NULL);
1399:     PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
1400:     PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
1401:     PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
1402:     PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
1403:     PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
1404:     PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
1405:   PetscOptionsEnd();

1407:   /* Choose the limiter from the list of registered limiters */
1408:   PetscFunctionListFind(limiters,lname,&ctx.limit);

1411:   /* Choose the physics from the list of registered models */
1412:   {
1413:     PetscErrorCode (*r)(FVCtx*);
1414:     PetscFunctionListFind(physics,physname,&r);
1416:     /* Create the physics, will set the number of fields and their names */
1417:     (*r)(&ctx);
1418:   }

1420:   /* Create a DMDA to manage the parallel grid */
1421:   DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);
1422:   DMSetFromOptions(da);
1423:   DMSetUp(da);
1424:   /* Inform the DMDA of the field names provided by the physics. */
1425:   /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1426:   for (i=0; i<ctx.physics.dof; i++) {
1427:     DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
1428:   }
1429:   DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1430:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

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

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

1439:   /* Create a vector to store the solution and to save the initial state */
1440:   DMCreateGlobalVector(da,&X);
1441:   VecDuplicate(X,&X0);
1442:   VecDuplicate(X,&R);

1444:   DMCreateMatrix(da,&B);

1446:   /* Create a time-stepping object */
1447:   TSCreate(comm,&ts);
1448:   TSSetDM(ts,da);
1449:   TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
1450:   TSSetIJacobian(ts,B,B,FVIJacobian,&ctx);
1451:   TSSetType(ts,TSSSP);
1452:   TSSetMaxTime(ts,10);
1453:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

1455:   /* Compute initial conditions and starting time step */
1456:   FVSample(&ctx,da,0,X0);
1457:   FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
1458:   VecCopy(X0,X);                        /* The function value was not used so we set X=X0 again */
1459:   TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
1460:   TSSetFromOptions(ts); /* Take runtime options */
1461:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1462:   {
1463:     PetscReal nrm1,nrmsup;
1464:     PetscInt  steps;

1466:     TSSolve(ts,X);
1467:     TSGetSolveTime(ts,&ptime);
1468:     TSGetStepNumber(ts,&steps);

1470:     PetscPrintf(comm,"Final time %8.5f, steps %D\n",(double)ptime,steps);
1471:     if (ctx.exact) {
1472:       SolutionErrorNorms(&ctx,da,ptime,X,&nrm1,&nrmsup);
1473:       PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e  ||x-x_e||_sup %8.4e\n",(double)nrm1,(double)nrmsup);
1474:     }
1475:   }

1477:   SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1478:   if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
1479:   if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
1480:   if (draw & 0x4) {
1481:     Vec Y;
1482:     VecDuplicate(X,&Y);
1483:     FVSample(&ctx,da,ptime,Y);
1484:     VecAYPX(Y,-1,X);
1485:     VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1486:     VecDestroy(&Y);
1487:   }

1489:   if (view_final) {
1490:     PetscViewer viewer;
1491:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1492:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1493:     VecView(X,viewer);
1494:     PetscViewerPopFormat(viewer);
1495:     PetscViewerDestroy(&viewer);
1496:   }

1498:   /* Clean up */
1499:   (*ctx.physics.destroy)(ctx.physics.user);
1500:   for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
1501:   PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
1502:   PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
1503:   VecDestroy(&X);
1504:   VecDestroy(&X0);
1505:   VecDestroy(&R);
1506:   MatDestroy(&B);
1507:   DMDestroy(&da);
1508:   TSDestroy(&ts);
1509:   PetscFunctionListDestroy(&limiters);
1510:   PetscFunctionListDestroy(&physics);
1511:   PetscFinalize();
1512:   return ierr;
1513: }

1515: /*TEST

1517:     build:
1518:       requires: !complex

1520:     test:
1521:       args: -da_grid_x 100 -initial 1 -xmin -2 -xmax 5 -exact -limit mc
1522:       requires: !complex !single

1524:     test:
1525:       suffix: 2
1526:       args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1
1527:       filter:  sed "s/at 48/at 0/g"
1528:       requires: !complex !single

1530:     test:
1531:       suffix: 3
1532:       args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1
1533:       nsize: 3
1534:       filter:  sed "s/at 48/at 0/g"
1535:       requires: !complex !single

1537: TEST*/
```