Actual source code: ex9.c
1: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
2: "Solves scalar and vector problems, choose the physical model with -physics\n"
3: " advection - Constant coefficient scalar advection\n"
4: " u_t + (a*u)_x = 0\n"
5: " burgers - Burgers equation\n"
6: " u_t + (u^2/2)_x = 0\n"
7: " traffic - Traffic equation\n"
8: " u_t + (u*(1-u))_x = 0\n"
9: " isogas - Isothermal gas dynamics\n"
10: " rho_t + (rho*u)_x = 0\n"
11: " (rho*u)_t + (rho*u^2 + c^2*rho)_x = 0\n"
12: " shallow - Shallow water equations\n"
13: " h_t + (h*u)_x = 0\n"
14: " (h*u)_t + (h*u^2 + g*h^2/2)_x = 0\n"
15: "Some of these physical models have multiple Riemann solvers, select these with -physics_xxx_riemann\n"
16: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
17: " the states across shocks and rarefactions\n"
18: " roe - Linearized scheme, usually with an entropy fix inside sonic rarefactions\n"
19: "The systems provide a choice of reconstructions with -physics_xxx_reconstruct\n"
20: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
21: " conservative - Limit the conservative variables directly, can cause undesired interaction of waves\n\n"
22: "A variety of limiters for high-resolution TVD limiters are available with -limit\n"
23: " upwind,minmod,superbee,mc,vanleer,vanalbada,koren,cada-torillhon (last two are nominally third order)\n"
24: " and non-TVD schemes lax-wendroff,beam-warming,fromm\n\n"
25: "To preserve the TVD property, one should time step with a strong stability preserving method.\n"
26: "The optimal high order explicit Runge-Kutta methods in TSSSP are recommended for non-stiff problems.\n\n"
27: "Several initial conditions can be chosen with -initial N\n\n"
28: "The problem size should be set with -da_grid_x M\n\n";
30: /* To get isfinite in math.h */
31: #define _XOPEN_SOURCE 600
33: #include <unistd.h> /* usleep */
34: #include petscts.h
35: #include petscda.h
37: #include "../src/mat/blockinvert.h" /* For the Kernel_*_gets_* stuff for BAIJ */
39: static inline PetscReal Sgn(PetscReal a) { return (a<0) ? -1 : 1; }
40: static inline PetscReal Abs(PetscReal a) { return (a<0) ? 0 : a; }
41: static inline PetscReal Sqr(PetscReal a) { return a*a; }
42: static inline PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; }
43: static inline PetscReal MinAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) < PetscAbs(b)) ? a : b; }
44: static inline PetscReal MinMod2(PetscReal a,PetscReal b)
45: { return (a*b<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscAbs(b)); }
46: static inline PetscReal MaxMod2(PetscReal a,PetscReal b)
47: { return (a*b<0) ? 0 : Sgn(a)*PetscMax(PetscAbs(a),PetscAbs(b)); }
48: static inline PetscReal MinMod3(PetscReal a,PetscReal b,PetscReal c)
49: {return (a*b<0 || a*c<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscMin(PetscAbs(b),PetscAbs(c))); }
51: static inline PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax)
52: { PetscReal range = xmax-xmin; return xmin + fmod(range+fmod(a,range),range); }
55: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
56: typedef struct _LimitInfo {
57: PetscReal hx;
58: PetscInt m;
59: } *LimitInfo;
60: static void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
61: {
62: PetscInt i;
63: for (i=0; i<info->m; i++) lmt[i] = 0;
64: }
65: static void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
66: {
67: PetscInt i;
68: for (i=0; i<info->m; i++) lmt[i] = jR[i];
69: }
70: static void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
71: {
72: PetscInt i;
73: for (i=0; i<info->m; i++) lmt[i] = jL[i];
74: }
75: static void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
76: {
77: PetscInt i;
78: for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i] + jR[i]);
79: }
80: static void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
81: {
82: PetscInt i;
83: for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]);
84: }
85: static void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
86: {
87: PetscInt i;
88: for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]));
89: }
90: static void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
91: {
92: PetscInt i;
93: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]);
94: }
95: static void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
96: { /* phi = (t + abs(t)) / (1 + abs(t)) */
97: PetscInt i;
98: 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);
99: }
100: static void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
101: { /* phi = (t + t^2) / (1 + t^2) */
102: PetscInt i;
103: 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);
104: }
105: static void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
106: { /* phi = (t + t^2) / (1 + t^2) */
107: PetscInt i;
108: for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0
109: : (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
110: }
111: static void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
112: { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
113: PetscInt i;
114: for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i]) + 2*Sqr(jL[i])*jR[i])
115: / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
116: }
117: static void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
118: { /* Symmetric version of above */
119: PetscInt i;
120: for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i])
121: / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
122: }
123: static void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
124: { /* Eq 11 of Cada-Torrilhon 2009 */
125: PetscInt i;
126: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]);
127: }
129: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)
130: { return PetscMax(0,PetscMin((L+2*R)/3,
131: PetscMax(-0.5*L,
132: PetscMin(2*L,
133: PetscMin((L+2*R)/3,1.6*R)))));
134: }
135: static void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
136: { /* Cada-Torrilhon 2009, Eq 13 */
137: PetscInt i;
138: for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]);
139: }
140: static void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
141: { /* Cada-Torrilhon 2009, Eq 22 */
142: /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
143: const PetscReal eps = 1e-7,hx = info->hx;
144: PetscInt i;
145: for (i=0; i<info->m; i++) {
146: const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r*hx);
147: lmt[i] = ((eta < 1-eps)
148: ? (jL[i] + 2*jR[i]) / 3
149: : ((eta > 1+eps)
150: ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i])
151: : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3
152: + (1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]))));
153: }
154: }
155: static void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
156: { Limit_CadaTorrilhon3R(0.1,info,jL,jR,lmt); }
157: static void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
158: { Limit_CadaTorrilhon3R(1,info,jL,jR,lmt); }
159: static void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
160: { Limit_CadaTorrilhon3R(10,info,jL,jR,lmt); }
161: static void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
162: { Limit_CadaTorrilhon3R(100,info,jL,jR,lmt); }
165: /* --------------------------------- Finite Volume data structures ----------------------------------- */
167: typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType;
168: static const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0};
169: typedef PetscErrorCode (*RiemannFunction)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*);
170: typedef PetscErrorCode (*ReconstructFunction)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*);
172: typedef struct {
173: PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*);
174: RiemannFunction riemann;
175: ReconstructFunction characteristic;
176: PetscErrorCode (*destroy)(void*);
177: void *user;
178: PetscInt dof;
179: char *fieldname[16];
180: } PhysicsCtx;
182: typedef struct {
183: void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*);
184: PhysicsCtx physics;
186: MPI_Comm comm;
187: char prefix[256];
188: DA da;
189: /* Local work arrays */
190: PetscScalar *R,*Rinv; /* Characteristic basis, and it's inverse. COLUMN-MAJOR */
191: PetscScalar *cjmpLR; /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
192: PetscScalar *cslope; /* Limited slope, written in characteristic basis */
193: PetscScalar *uLR; /* Solution at left and right of interface, conservative variables, len=2*dof */
194: PetscScalar *flux; /* Flux across interface */
196: PetscReal cfl_idt; /* Max allowable value of 1/Delta t */
197: PetscReal cfl;
198: PetscReal xmin,xmax;
199: PetscInt initial;
200: PetscTruth exact;
201: FVBCType bctype;
202: } FVCtx;
205: /* Utility */
209: PetscErrorCode RiemannListAdd(PetscFList *flist,const char *name,RiemannFunction rsolve)
210: {
214: PetscFListAdd(flist,name,"",(void(*)(void))rsolve);
215: return(0);
216: }
220: PetscErrorCode RiemannListFind(PetscFList flist,const char *name,RiemannFunction *rsolve)
221: {
225: PetscFListFind(flist,PETSC_COMM_WORLD,name,(void(**)(void))rsolve);
226: if (!*rsolve) SETERRQ1(1,"Riemann solver \"%s\" could not be found",name);
227: return(0);
228: }
232: PetscErrorCode ReconstructListAdd(PetscFList *flist,const char *name,ReconstructFunction r)
233: {
237: PetscFListAdd(flist,name,"",(void(*)(void))r);
238: return(0);
239: }
243: PetscErrorCode ReconstructListFind(PetscFList flist,const char *name,ReconstructFunction *r)
244: {
248: PetscFListFind(flist,PETSC_COMM_WORLD,name,(void(**)(void))r);
249: if (!*r) SETERRQ1(1,"Reconstruction \"%s\" could not be found",name);
250: return(0);
251: }
254: /* --------------------------------- Physics ----------------------------------- */
255: /**
256: * Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction. These
257: * are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the
258: * number of fields and their names, and a function to deallocate private storage.
259: **/
261: /* First a few functions useful to several different physics */
264: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi)
265: {
266: PetscInt i,j;
269: for (i=0; i<m; i++) {
270: for (j=0; j<m; j++) {
271: Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j);
272: }
273: }
274: return(0);
275: }
279: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
280: {
284: PetscFree(vctx);
285: return(0);
286: }
290: /* --------------------------------- Advection ----------------------------------- */
292: typedef struct {
293: PetscReal a; /* advective velocity */
294: } AdvectCtx;
298: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
299: {
300: AdvectCtx *ctx = (AdvectCtx*)vctx;
301: PetscReal speed;
304: speed = ctx->a;
305: flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
306: *maxspeed = speed;
307: return(0);
308: }
312: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
313: {
314: AdvectCtx *ctx = (AdvectCtx*)vctx;
315: PetscReal a = ctx->a,x0;
318: switch (bctype) {
319: case FVBC_OUTFLOW: x0 = x-a*t; break;
320: case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
321: default: SETERRQ(1,"unknown BCType");
322: }
323: switch (initial) {
324: case 0: u[0] = (x0 < 0) ? 1 : -1; break;
325: case 1: u[0] = (x0 < 0) ? -1 : 1; break;
326: case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
327: case 3: u[0] = sin(2*M_PI*x0); break;
328: case 4: u[0] = PetscAbs(x0); break;
329: default: SETERRQ(1,"unknown initial condition");
330: }
331: return(0);
332: }
336: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
337: {
339: AdvectCtx *user;
342: PetscNew(*user,&user);
343: ctx->physics.sample = PhysicsSample_Advect;
344: ctx->physics.riemann = PhysicsRiemann_Advect;
345: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
346: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
347: ctx->physics.user = user;
348: ctx->physics.dof = 1;
349: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
350: user->a = 1;
351: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
352: {
353: PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,PETSC_NULL);
354: }
355: PetscOptionsEnd();
356: return(0);
357: }
361: /* --------------------------------- Burgers ----------------------------------- */
363: typedef struct {
364: PetscReal lxf_speed;
365: } BurgersCtx;
369: static PetscErrorCode PhysicsSample_Burgers(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
370: {
373: if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(1,"Exact solution not implemented for periodic");
374: switch (initial) {
375: case 0: u[0] = (x < 0) ? 1 : -1; break;
376: case 1:
377: if (x < -t) u[0] = -1;
378: else if (x < t) u[0] = x/t;
379: else u[0] = 1;
380: break;
381: case 2:
382: if (x < 0) u[0] = 0;
383: else if (x <= t) u[0] = x/t;
384: else if (x < 1+0.5*t) u[0] = 1;
385: else u[0] = 0;
386: break;
387: case 3:
388: if (x < 0.2*t) u[0] = 0.2;
389: else if (x < t) u[0] = x/t;
390: else u[0] = 1;
391: break;
392: case 4:
393: if (t > 0) SETERRQ(1,"Only initial condition available");
394: u[0] = 0.7 + 0.3*sin(2*M_PI*((x-xmin)/(xmax-xmin)));
395: break;
396: default: SETERRQ(1,"unknown initial condition");
397: }
398: return(0);
399: }
403: static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
404: {
407: if (uL[0] < uR[0]) { /* rarefaction */
408: flux[0] = (uL[0]*uR[0] < 0)
409: ? 0 /* sonic rarefaction */
410: : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0]));
411: } else { /* shock */
412: flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0]));
413: }
414: *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
415: return(0);
416: }
420: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
421: {
422: PetscReal speed;
425: speed = 0.5*(uL[0] + uR[0]);
426: flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
427: if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
428: *maxspeed = speed;
429: return(0);
430: }
434: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
435: {
436: PetscReal c;
437: PetscScalar fL,fR;
440: c = ((BurgersCtx*)vctx)->lxf_speed;
441: fL = 0.5*PetscSqr(uL[0]);
442: fR = 0.5*PetscSqr(uR[0]);
443: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
444: *maxspeed = c;
445: return(0);
446: }
450: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
451: {
452: PetscReal c;
453: PetscScalar fL,fR;
456: c = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0]));
457: fL = 0.5*PetscSqr(uL[0]);
458: fR = 0.5*PetscSqr(uR[0]);
459: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
460: *maxspeed = c;
461: return(0);
462: }
466: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
467: {
468: BurgersCtx *user;
470: RiemannFunction r;
471: PetscFList rlist = 0;
472: char rname[256] = "exact";
475: PetscNew(*user,&user);
476: ctx->physics.sample = PhysicsSample_Burgers;
477: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
478: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
479: ctx->physics.user = user;
480: ctx->physics.dof = 1;
481: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
482: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Burgers_Exact);
483: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Burgers_Roe);
484: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Burgers_LxF);
485: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Burgers_Rusanov);
486: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
487: {
488: PetscOptionsList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
489: }
490: PetscOptionsEnd();
491: RiemannListFind(rlist,rname,&r);
492: PetscFListDestroy(&rlist);
493: ctx->physics.riemann = r;
495: /* *
496: * Hack to deal with LxF in semi-discrete form
497: * max speed is 1 for the basic initial conditions (where |u| <= 1)
498: * */
499: if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
500: return(0);
501: }
505: /* --------------------------------- Traffic ----------------------------------- */
507: typedef struct {
508: PetscReal lxf_speed;
509: PetscReal a;
510: } TrafficCtx;
512: static inline PetscScalar TrafficFlux(PetscScalar a,PetscScalar u) { return a*u*(1-u); }
516: static PetscErrorCode PhysicsSample_Traffic(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
517: {
518: PetscReal a = ((TrafficCtx*)vctx)->a;
521: if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(1,"Exact solution not implemented for periodic");
522: switch (initial) {
523: case 0:
524: u[0] = (-a*t < x) ? 2 : 0; break;
525: case 1:
526: if (x < PetscMin(2*a*t,0.5+a*t)) u[0] = -1;
527: else if (x < 1) u[0] = 0;
528: else u[0] = 1;
529: break;
530: case 2:
531: if (t > 0) SETERRQ(1,"Only initial condition available");
532: u[0] = 0.7 + 0.3*sin(2*M_PI*((x-xmin)/(xmax-xmin)));
533: break;
534: default: SETERRQ(1,"unknown initial condition");
535: }
536: return(0);
537: }
541: static PetscErrorCode PhysicsRiemann_Traffic_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
542: {
543: PetscReal a = ((TrafficCtx*)vctx)->a;
546: if (uL[0] < uR[0]) {
547: flux[0] = PetscMin(TrafficFlux(a,uL[0]),
548: TrafficFlux(a,uR[0]));
549: } else {
550: flux[0] = (uR[0] < 0.5 && 0.5 < uL[0])
551: ? TrafficFlux(a,0.5)
552: : PetscMax(TrafficFlux(a,uL[0]),
553: TrafficFlux(a,uR[0]));
554: }
555: *maxspeed = a*MaxAbs(1-2*uL[0],1-2*uR[0]);
556: return(0);
557: }
561: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
562: {
563: PetscReal a = ((TrafficCtx*)vctx)->a;
564: PetscReal speed;
567: speed = a*(1 - (uL[0] + uR[0]));
568: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
569: *maxspeed = speed;
570: return(0);
571: }
575: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
576: {
577: TrafficCtx *phys = (TrafficCtx*)vctx;
578: PetscReal a = phys->a;
579: PetscReal speed;
582: speed = a*(1 - (uL[0] + uR[0]));
583: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]);
584: *maxspeed = speed;
585: return(0);
586: }
590: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
591: {
592: PetscReal a = ((TrafficCtx*)vctx)->a;
593: PetscReal speed;
596: speed = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0]));
597: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]);
598: *maxspeed = speed;
599: return(0);
600: }
604: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
605: {
607: TrafficCtx *user;
608: RiemannFunction r;
609: PetscFList rlist = 0;
610: char rname[256] = "exact";
613: PetscNew(*user,&user);
614: ctx->physics.sample = PhysicsSample_Traffic;
615: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
616: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
617: ctx->physics.user = user;
618: ctx->physics.dof = 1;
619: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
620: user->a = 0.5;
621: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Traffic_Exact);
622: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Traffic_Roe);
623: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Traffic_LxF);
624: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Traffic_Rusanov);
625: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");
626: {
627: PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,PETSC_NULL);
628: PetscOptionsList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
629: }
630: PetscOptionsEnd();
632: RiemannListFind(rlist,rname,&r);
633: PetscFListDestroy(&rlist);
634: ctx->physics.riemann = r;
636: /* *
637: * Hack to deal with LxF in semi-discrete form
638: * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
639: * */
640: if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3*user->a;
642: return(0);
643: }
648: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */
650: typedef struct {
651: PetscReal acoustic_speed;
652: } IsoGasCtx;
654: static inline void IsoGasFlux(PetscReal c,const PetscScalar *u,PetscScalar *f)
655: {
656: f[0] = u[1];
657: f[1] = PetscSqr(u[1])/u[0] + c*c*u[0];
658: }
662: static PetscErrorCode PhysicsSample_IsoGas(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
663: {
666: if (t > 0) SETERRQ(1,"Exact solutions not implemented for t > 0");
667: switch (initial) {
668: case 0:
669: u[0] = (x < 0) ? 1 : 0.5;
670: u[1] = (x < 0) ? 1 : 0.7;
671: break;
672: case 1:
673: u[0] = 1+0.5*sin(2*M_PI*x);
674: u[1] = 1*u[0];
675: break;
676: default: SETERRQ(1,"unknown initial condition");
677: }
678: return(0);
679: }
683: static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
684: {
685: IsoGasCtx *phys = (IsoGasCtx*)vctx;
686: PetscReal c = phys->acoustic_speed;
687: PetscScalar ubar,du[2],a[2],fL[2],fR[2],lam[2],ustar[2],R[2][2];
688: PetscInt i;
691: ubar = (uL[1]/PetscSqrtScalar(uL[0]) + uR[1]/PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0]));
692: /* write fluxuations in characteristic basis */
693: du[0] = uR[0] - uL[0];
694: du[1] = uR[1] - uL[1];
695: a[0] = (1/(2*c)) * ((ubar + c)*du[0] - du[1]);
696: a[1] = (1/(2*c)) * ((-ubar + c)*du[0] + du[1]);
697: /* wave speeds */
698: lam[0] = ubar - c;
699: lam[1] = ubar + c;
700: /* Right eigenvectors */
701: R[0][0] = 1; R[0][1] = ubar-c;
702: R[1][0] = 1; R[1][1] = ubar+c;
703: /* Compute state in star region (between the 1-wave and 2-wave) */
704: for (i=0; i<2; i++) ustar[i] = uL[i] + a[0]*R[0][i];
705: if (uL[1]/uL[0] < c && c < ustar[1]/ustar[0]) { /* 1-wave is sonic rarefaction */
706: PetscScalar ufan[2];
707: ufan[0] = uL[0]*PetscExpScalar(uL[1]/(uL[0]*c) - 1);
708: ufan[1] = c*ufan[0];
709: IsoGasFlux(c,ufan,flux);
710: } else if (ustar[1]/ustar[0] < -c && -c < uR[1]/uR[0]) { /* 2-wave is sonic rarefaction */
711: PetscScalar ufan[2];
712: ufan[0] = uR[0]*PetscExpScalar(-uR[1]/(uR[0]*c) - 1);
713: ufan[1] = -c*ufan[0];
714: IsoGasFlux(c,ufan,flux);
715: } else { /* Centered form */
716: IsoGasFlux(c,uL,fL);
717: IsoGasFlux(c,uR,fR);
718: for (i=0; i<2; i++) {
719: PetscScalar absdu = PetscAbsScalar(lam[0])*a[0]*R[0][i] + PetscAbsScalar(lam[1])*a[1]*R[1][i];
720: flux[i] = 0.5*(fL[i]+fR[i]) - 0.5*absdu;
721: }
722: }
723: *maxspeed = MaxAbs(lam[0],lam[1]);
724: return(0);
725: }
729: static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
730: {
731: IsoGasCtx *phys = (IsoGasCtx*)vctx;
732: PetscReal c = phys->acoustic_speed;
733: PetscScalar ustar[2];
734: struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
735: PetscInt i;
738: if (!(L.rho > 0 && R.rho > 0)) SETERRQ(1,"Reconstructed density is negative");
739: {
740: /* Solve for star state */
741: PetscScalar res,tmp,rho = 0.5*(L.rho + R.rho); /* initial guess */
742: for (i=0; i<20; i++) {
743: PetscScalar fr,fl,dfr,dfl;
744: fl = (L.rho < rho)
745: ? (rho-L.rho)/PetscSqrtScalar(L.rho*rho) /* shock */
746: : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
747: fr = (R.rho < rho)
748: ? (rho-R.rho)/PetscSqrtScalar(R.rho*rho) /* shock */
749: : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
750: res = R.u-L.u + c*(fr+fl);
751: if (!isfinite(res)) SETERRQ1(1,"non-finite residual=%g",res);
752: if (PetscAbsScalar(res) < 1e-10) {
753: star.rho = rho;
754: star.u = L.u - c*fl;
755: goto converged;
756: }
757: dfl = (L.rho < rho)
758: ? 1/PetscSqrtScalar(L.rho*rho)*(1 - 0.5*(rho-L.rho)/rho)
759: : 1/rho;
760: dfr = (R.rho < rho)
761: ? 1/PetscSqrtScalar(R.rho*rho)*(1 - 0.5*(rho-R.rho)/rho)
762: : 1/rho;
763: tmp = rho - res/(c*(dfr+dfl));
764: if (tmp <= 0) rho /= 2; /* Guard against Newton shooting off to a negative density */
765: else rho = tmp;
766: if (!((rho > 0) && isnormal(rho))) SETERRQ1(1,"non-normal iterate rho=%g",rho);
767: }
768: SETERRQ1(1,"Newton iteration for star.rho diverged after %d iterations",i);
769: }
770: converged:
771: if (L.u-c < 0 && 0 < star.u-c) { /* 1-wave is sonic rarefaction */
772: PetscScalar ufan[2];
773: ufan[0] = L.rho*PetscExpScalar(L.u/c - 1);
774: ufan[1] = c*ufan[0];
775: IsoGasFlux(c,ufan,flux);
776: } else if (star.u+c < 0 && 0 < R.u+c) { /* 2-wave is sonic rarefaction */
777: PetscScalar ufan[2];
778: ufan[0] = R.rho*PetscExpScalar(-R.u/c - 1);
779: ufan[1] = -c*ufan[0];
780: IsoGasFlux(c,ufan,flux);
781: } else if ((L.rho >= star.rho && L.u-c >= 0)
782: || (L.rho < star.rho && (star.rho*star.u-L.rho*L.u)/(star.rho-L.rho) > 0)) {
783: /* 1-wave is supersonic rarefaction, or supersonic shock */
784: IsoGasFlux(c,uL,flux);
785: } else if ((star.rho <= R.rho && R.u+c <= 0)
786: || (star.rho > R.rho && (R.rho*R.u-star.rho*star.u)/(R.rho-star.rho) < 0)) {
787: /* 2-wave is supersonic rarefaction or supersonic shock */
788: IsoGasFlux(c,uR,flux);
789: } else {
790: ustar[0] = star.rho;
791: ustar[1] = star.rho*star.u;
792: IsoGasFlux(c,ustar,flux);
793: }
794: *maxspeed = MaxAbs(MaxAbs(star.u-c,star.u+c),MaxAbs(L.u-c,R.u+c));
795: return(0);
796: }
800: static PetscErrorCode PhysicsRiemann_IsoGas_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
801: {
802: IsoGasCtx *phys = (IsoGasCtx*)vctx;
803: PetscScalar c = phys->acoustic_speed,fL[2],fR[2],s;
804: struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
807: if (!(L.rho > 0 && R.rho > 0)) SETERRQ(1,"Reconstructed density is negative");
808: IsoGasFlux(c,uL,fL);
809: IsoGasFlux(c,uR,fR);
810: s = PetscMax(PetscAbs(L.u),PetscAbs(R.u))+c;
811: flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
812: flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
813: *maxspeed = s;
814: return(0);
815: }
819: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi)
820: {
821: IsoGasCtx *phys = (IsoGasCtx*)vctx;
822: PetscReal c = phys->acoustic_speed;
826: X[0*2+0] = 1;
827: X[0*2+1] = u[1]/u[0] - c;
828: X[1*2+0] = 1;
829: X[1*2+1] = u[1]/u[0] + c;
830: PetscMemcpy(Xi,X,4*sizeof(X[0]));
831: Kernel_A_gets_inverse_A_2(Xi,0);
832: return(0);
833: }
837: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
838: {
840: IsoGasCtx *user;
841: PetscFList rlist = 0,rclist = 0;
842: char rname[256] = "exact",rcname[256] = "characteristic";
845: PetscNew(*user,&user);
846: ctx->physics.sample = PhysicsSample_IsoGas;
847: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
848: ctx->physics.user = user;
849: ctx->physics.dof = 2;
850: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
851: PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
852: user->acoustic_speed = 1;
853: RiemannListAdd(&rlist,"exact", PhysicsRiemann_IsoGas_Exact);
854: RiemannListAdd(&rlist,"roe", PhysicsRiemann_IsoGas_Roe);
855: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_IsoGas_Rusanov);
856: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_IsoGas);
857: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
858: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");
859: {
860: PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,PETSC_NULL);
861: PetscOptionsList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
862: PetscOptionsList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
863: }
864: PetscOptionsEnd();
865: RiemannListFind(rlist,rname,&ctx->physics.riemann);
866: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
867: PetscFListDestroy(&rlist);
868: PetscFListDestroy(&rclist);
869: return(0);
870: }
874: /* --------------------------------- Shallow Water ----------------------------------- */
876: typedef struct {
877: PetscReal gravity;
878: } ShallowCtx;
880: static inline void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
881: {
882: f[0] = u[1];
883: f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
884: }
888: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
889: {
890: ShallowCtx *phys = (ShallowCtx*)vctx;
891: PetscScalar g = phys->gravity,ustar[2],cL,cR,c,cstar;
892: struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
893: PetscInt i;
896: if (!(L.h > 0 && R.h > 0)) SETERRQ(1,"Reconstructed thickness is negative");
897: cL = PetscSqrtScalar(g*L.h);
898: cR = PetscSqrtScalar(g*R.h);
899: c = PetscMax(cL,cR);
900: {
901: /* Solve for star state */
902: const PetscInt maxits = 50;
903: PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */
904: h0 = h;
905: for (i=0; i<maxits; i++) {
906: PetscScalar fr,fl,dfr,dfl;
907: fl = (L.h < h)
908: ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */
909: : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h); /* rarefaction */
910: fr = (R.h < h)
911: ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */
912: : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h); /* rarefaction */
913: res = R.u - L.u + fr + fl;
914: if (!isfinite(res)) SETERRQ1(1,"non-finite residual=%g",res);
915: //PetscPrintf(PETSC_COMM_WORLD,"h=%g, res[%d] = %g\n",h,i,res);
916: if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h-h0) < 1e-8)) {
917: star.h = h;
918: star.u = L.u - fl;
919: goto converged;
920: } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */
921: h = 0.8*h0 + 0.2*h;
922: continue;
923: }
924: /* Accept the last step and take another */
925: res0 = res;
926: h0 = h;
927: dfl = (L.h < h)
928: ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h)
929: : PetscSqrtScalar(g/h);
930: dfr = (R.h < h)
931: ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h)
932: : PetscSqrtScalar(g/h);
933: tmp = h - res/(dfr+dfl);
934: if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */
935: else h = tmp;
936: if (!((h > 0) && isnormal(h))) SETERRQ1(1,"non-normal iterate h=%g",h);
937: }
938: SETERRQ1(1,"Newton iteration for star.h diverged after %d iterations",i);
939: }
940: converged:
941: cstar = PetscSqrtScalar(g*star.h);
942: if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */
943: PetscScalar ufan[2];
944: ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL);
945: ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0];
946: ShallowFlux(phys,ufan,flux);
947: } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */
948: PetscScalar ufan[2];
949: ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR);
950: ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0];
951: ShallowFlux(phys,ufan,flux);
952: } else if ((L.h >= star.h && L.u-c >= 0)
953: || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) {
954: /* 1-wave is right-travelling shock (supersonic) */
955: ShallowFlux(phys,uL,flux);
956: } else if ((star.h <= R.h && R.u+c <= 0)
957: || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) {
958: /* 2-wave is left-travelling shock (supersonic) */
959: ShallowFlux(phys,uR,flux);
960: } else {
961: ustar[0] = star.h;
962: ustar[1] = star.h*star.u;
963: ShallowFlux(phys,ustar,flux);
964: }
965: *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR));
966: return(0);
967: }
971: static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
972: {
973: ShallowCtx *phys = (ShallowCtx*)vctx;
974: PetscScalar g = phys->gravity,fL[2],fR[2],s;
975: struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
978: if (!(L.h > 0 && R.h > 0)) SETERRQ(1,"Reconstructed thickness is negative");
979: ShallowFlux(phys,uL,fL);
980: ShallowFlux(phys,uR,fR);
981: s = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h));
982: flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
983: flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
984: *maxspeed = s;
985: return(0);
986: }
990: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi)
991: {
992: ShallowCtx *phys = (ShallowCtx*)vctx;
993: PetscReal c;
997: c = PetscSqrtScalar(u[0]*phys->gravity);
998: X[0*2+0] = 1;
999: X[0*2+1] = u[1]/u[0] - c;
1000: X[1*2+0] = 1;
1001: X[1*2+1] = u[1]/u[0] + c;
1002: PetscMemcpy(Xi,X,4*sizeof(X[0]));
1003: Kernel_A_gets_inverse_A_2(Xi,0);
1004: return(0);
1005: }
1009: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
1010: {
1012: ShallowCtx *user;
1013: PetscFList rlist = 0,rclist = 0;
1014: char rname[256] = "exact",rcname[256] = "characteristic";
1017: PetscNew(*user,&user);
1018: /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
1019: ctx->physics.sample = PhysicsSample_IsoGas;
1020: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
1021: ctx->physics.user = user;
1022: ctx->physics.dof = 2;
1023: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1024: PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1025: user->gravity = 1;
1026: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Shallow_Exact);
1027: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov);
1028: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Shallow);
1029: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1030: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");
1031: {
1032: PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,PETSC_NULL);
1033: PetscOptionsList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
1034: PetscOptionsList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
1035: }
1036: PetscOptionsEnd();
1037: RiemannListFind(rlist,rname,&ctx->physics.riemann);
1038: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1039: PetscFListDestroy(&rlist);
1040: PetscFListDestroy(&rclist);
1041: return(0);
1042: }
1046: /* --------------------------------- Finite Volume Solver ----------------------------------- */
1050: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
1051: {
1052: FVCtx *ctx = (FVCtx*)vctx;
1053: PetscErrorCode ierr;
1054: PetscInt i,j,k,Mx,dof,xs,xm;
1055: PetscReal hx,cfl_idt = 0;
1056: PetscScalar *x,*f,*slope;
1057: Vec Xloc;
1060: DAGetLocalVector(ctx->da,&Xloc);
1061: DAGetInfo(ctx->da,0, &Mx,0,0, 0,0,0, &dof,0,0,0);
1062: hx = (ctx->xmax - ctx->xmin)/Mx;
1063: DAGlobalToLocalBegin(ctx->da,X,INSERT_VALUES,Xloc);
1064: DAGlobalToLocalEnd (ctx->da,X,INSERT_VALUES,Xloc);
1066: VecZeroEntries(F);
1068: DAVecGetArray(ctx->da,Xloc,&x);
1069: DAVecGetArray(ctx->da,F,&f);
1070: DAGetArray(ctx->da,PETSC_TRUE,(void**)&slope);
1072: DAGetCorners(ctx->da,&xs,0,0,&xm,0,0);
1074: if (ctx->bctype == FVBC_OUTFLOW) {
1075: for (i=xs-2; i<0; i++) {
1076: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
1077: }
1078: for (i=Mx; i<xs+xm+2; i++) {
1079: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
1080: }
1081: }
1082: for (i=xs-1; i<xs+xm+1; i++) {
1083: struct _LimitInfo info;
1084: PetscScalar *cjmpL,*cjmpR;
1085: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
1086: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv);
1087: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
1088: PetscMemzero(ctx->cjmpLR,2*dof*sizeof(ctx->cjmpLR[0]));
1089: cjmpL = &ctx->cjmpLR[0];
1090: cjmpR = &ctx->cjmpLR[dof];
1091: for (j=0; j<dof; j++) {
1092: PetscScalar jmpL,jmpR;
1093: jmpL = x[(i+0)*dof+j] - x[(i-1)*dof+j];
1094: jmpR = x[(i+1)*dof+j] - x[(i+0)*dof+j];
1095: for (k=0; k<dof; k++) {
1096: cjmpL[k] += ctx->Rinv[k+j*dof] * jmpL;
1097: cjmpR[k] += ctx->Rinv[k+j*dof] * jmpR;
1098: }
1099: }
1100: /* Apply limiter to the left and right characteristic jumps */
1101: info.m = dof;
1102: info.hx = hx;
1103: (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
1104: for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
1105: for (j=0; j<dof; j++) {
1106: PetscScalar tmp = 0;
1107: for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof] * ctx->cslope[k];
1108: slope[i*dof+j] = tmp;
1109: }
1110: }
1112: for (i=xs; i<xs+xm+1; i++) {
1113: PetscReal maxspeed;
1114: PetscScalar *uL,*uR;
1115: uL = &ctx->uLR[0];
1116: uR = &ctx->uLR[dof];
1117: for (j=0; j<dof; j++) {
1118: uL[j] = x[(i-1)*dof+j] + slope[(i-1)*dof+j]*hx/2;
1119: uR[j] = x[(i-0)*dof+j] - slope[(i-0)*dof+j]*hx/2;
1120: }
1121: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed);
1122: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
1124: if (i > xs) {
1125: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
1126: }
1127: if (i < xs+xm) {
1128: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
1129: }
1130: }
1132: DAVecRestoreArray(ctx->da,Xloc,&x);
1133: DAVecRestoreArray(ctx->da,F,&f);
1134: DARestoreArray(ctx->da,PETSC_TRUE,(void**)&slope);
1135: DARestoreLocalVector(ctx->da,&Xloc);
1137: PetscGlobalMax(&cfl_idt,&ctx->cfl_idt,((PetscObject)ctx->da)->comm);
1138: if (0) {
1139: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
1140: PetscReal dt,tnow;
1141: TSGetTimeStep(ts,&dt);
1142: TSGetTime(ts,&tnow);
1143: if (dt > 0.5/ctx->cfl_idt) {
1144: if (1) {
1145: PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",tnow,dt,0.5/ctx->cfl_idt);
1146: } else {
1147: SETERRQ2(1,"Stability constraint exceeded, %g > %g",dt,ctx->cfl/ctx->cfl_idt);
1148: }
1149: }
1150: }
1151: return(0);
1152: }
1156: static PetscErrorCode FVSample(FVCtx *ctx,PetscReal time,Vec U)
1157: {
1159: PetscScalar *u,*uj;
1160: PetscInt i,j,k,dof,xs,xm,Mx;
1163: if (!ctx->physics.sample) SETERRQ(1,"Physics has not provided a sampling function");
1164: DAGetInfo(ctx->da,0, &Mx,0,0, 0,0,0, &dof,0,0,0);
1165: DAGetCorners(ctx->da,&xs,0,0,&xm,0,0);
1166: DAVecGetArray(ctx->da,U,&u);
1167: PetscMalloc(dof*sizeof uj[0],&uj);
1168: for (i=xs; i<xs+xm; i++) {
1169: const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h;
1170: const PetscInt N = 200;
1171: /* Integrate over cell i using trapezoid rule with N points. */
1172: for (k=0; k<dof; k++) u[i*dof+k] = 0;
1173: for (j=0; j<N+1; j++) {
1174: PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N;
1175: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
1176: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N)?0.5:1.0)*uj[k]/N;
1177: }
1178: }
1179: DAVecRestoreArray(ctx->da,U,&u);
1180: PetscFree(uj);
1181: return(0);
1182: }
1186: static PetscErrorCode SolutionStatsView(DA da,Vec X,PetscViewer viewer)
1187: {
1189: PetscReal xmin,xmax;
1190: PetscScalar sum,*x,tvsum,tvgsum;
1191: PetscInt imin,imax,Mx,i,j,xs,xm,dof;
1192: Vec Xloc;
1193: PetscTruth iascii;
1196: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
1197: if (iascii) {
1198: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
1199: DAGetLocalVector(da,&Xloc);
1200: DAGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1201: DAGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
1202: DAVecGetArray(da,Xloc,&x);
1203: DAGetCorners(da,&xs,0,0,&xm,0,0);
1204: DAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0);
1205: tvsum = 0;
1206: for (i=xs; i<xs+xm; i++) {
1207: for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j] - x[(i-1)*dof+j]);
1208: }
1209: PetscGlobalMax(&tvsum,&tvgsum,((PetscObject)da)->comm);
1210: DAVecRestoreArray(da,Xloc,&x);
1211: DARestoreLocalVector(da,&Xloc);
1213: VecMin(X,&imin,&xmin);
1214: VecMax(X,&imax,&xmax);
1215: VecSum(X,&sum);
1216: 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);
1217: } else {
1218: SETERRQ(1,"Viewer type not supported");
1219: }
1220: return(0);
1221: }
1225: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,PetscReal t,Vec X,PetscReal *nrm1,PetscReal *nrmsup)
1226: {
1228: Vec Y;
1229: PetscInt Mx;
1232: VecGetSize(X,&Mx);
1233: VecDuplicate(X,&Y);
1234: FVSample(ctx,t,Y);
1235: VecAYPX(Y,-1,X);
1236: VecNorm(Y,NORM_1,nrm1);
1237: VecNorm(Y,NORM_INFINITY,nrmsup);
1238: *nrm1 /= Mx;
1239: VecDestroy(Y);
1240: return(0);
1241: }
1246: int main(int argc,char *argv[])
1247: {
1248: char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
1249: PetscFList limiters = 0,physics = 0;
1250: MPI_Comm comm;
1251: TS ts;
1252: Vec X,X0;
1253: FVCtx ctx;
1254: PetscInt i,dof,xs,xm,Mx,draw = 0;
1255: PetscTruth view_final = PETSC_FALSE;
1256: PetscReal ptime;
1259: PetscInitialize(&argc,&argv,0,help);
1260: comm = PETSC_COMM_WORLD;
1261: PetscMemzero(&ctx,sizeof(ctx));
1263: /* Register limiters to be available on the command line */
1264: PetscFListAdd(&limiters,"upwind" ,"",(void(*)(void))Limit_Upwind);
1265: PetscFListAdd(&limiters,"lax-wendroff" ,"",(void(*)(void))Limit_LaxWendroff);
1266: PetscFListAdd(&limiters,"beam-warming" ,"",(void(*)(void))Limit_BeamWarming);
1267: PetscFListAdd(&limiters,"fromm" ,"",(void(*)(void))Limit_Fromm);
1268: PetscFListAdd(&limiters,"minmod" ,"",(void(*)(void))Limit_Minmod);
1269: PetscFListAdd(&limiters,"superbee" ,"",(void(*)(void))Limit_Superbee);
1270: PetscFListAdd(&limiters,"mc" ,"",(void(*)(void))Limit_MC);
1271: PetscFListAdd(&limiters,"vanleer" ,"",(void(*)(void))Limit_VanLeer);
1272: PetscFListAdd(&limiters,"vanalbada" ,"",(void(*)(void))Limit_VanAlbada);
1273: PetscFListAdd(&limiters,"vanalbadatvd" ,"",(void(*)(void))Limit_VanAlbadaTVD);
1274: PetscFListAdd(&limiters,"koren" ,"",(void(*)(void))Limit_Koren);
1275: PetscFListAdd(&limiters,"korensym" ,"",(void(*)(void))Limit_KorenSym);
1276: PetscFListAdd(&limiters,"koren3" ,"",(void(*)(void))Limit_Koren3);
1277: PetscFListAdd(&limiters,"cada-torrilhon2" ,"",(void(*)(void))Limit_CadaTorrilhon2);
1278: PetscFListAdd(&limiters,"cada-torrilhon3-r0p1","",(void(*)(void))Limit_CadaTorrilhon3R0p1);
1279: PetscFListAdd(&limiters,"cada-torrilhon3-r1" ,"",(void(*)(void))Limit_CadaTorrilhon3R1);
1280: PetscFListAdd(&limiters,"cada-torrilhon3-r10" ,"",(void(*)(void))Limit_CadaTorrilhon3R10);
1281: PetscFListAdd(&limiters,"cada-torrilhon3-r100","",(void(*)(void))Limit_CadaTorrilhon3R100);
1283: /* Register physical models to be available on the command line */
1284: PetscFListAdd(&physics,"advect" ,"",(void(*)(void))PhysicsCreate_Advect);
1285: PetscFListAdd(&physics,"burgers" ,"",(void(*)(void))PhysicsCreate_Burgers);
1286: PetscFListAdd(&physics,"traffic" ,"",(void(*)(void))PhysicsCreate_Traffic);
1287: PetscFListAdd(&physics,"isogas" ,"",(void(*)(void))PhysicsCreate_IsoGas);
1288: PetscFListAdd(&physics,"shallow" ,"",(void(*)(void))PhysicsCreate_Shallow);
1290: ctx.cfl = 0.9; ctx.bctype = FVBC_PERIODIC;
1291: ctx.xmin = 0; ctx.xmax = 1;
1292: PetscOptionsBegin(comm,PETSC_NULL,"Finite Volume solver options","");
1293: {
1294: PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,PETSC_NULL);
1295: PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,PETSC_NULL);
1296: PetscOptionsList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),PETSC_NULL);
1297: PetscOptionsList("-physics","Name of physics (Riemann solver and characteristics) to use","",physics,physname,physname,sizeof(physname),PETSC_NULL);
1298: PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,PETSC_NULL);
1299: PetscOptionsString("-view_final","Write final solution in ASCII Matlab format to given file name","",final_fname,final_fname,sizeof final_fname,&view_final);
1300: PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,PETSC_NULL);
1301: PetscOptionsTruth("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,PETSC_NULL);
1302: PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,PETSC_NULL);
1303: PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,ctx.bctype,(PetscEnum*)&ctx.bctype,PETSC_NULL);
1304: }
1305: PetscOptionsEnd();
1307: /* Choose the limiter from the list of registered limiters */
1308: PetscFListFind(limiters,comm,lname,(void(**)(void))&ctx.limit);
1309: if (!ctx.limit) SETERRQ1(1,"Limiter '%s' not found",lname);
1311: /* Choose the physics from the list of registered models */
1312: {
1313: PetscErrorCode (*r)(FVCtx*);
1314: PetscFListFind(physics,comm,physname,(void(**)(void))&r);
1315: if (!r) SETERRQ1(1,"Physics '%s' not found",physname);
1316: /* Create the physics, will set the number of fields and their names */
1317: (*r)(&ctx);
1318: }
1320: /* Create a DA to manage the parallel grid */
1321: DACreate1d(comm,DA_XPERIODIC,-50,ctx.physics.dof,2,PETSC_NULL,&ctx.da);
1322: /* Inform the DA of the field names provided by the physics. */
1323: /* The names will be shown in the title bars when run with -ts_monitor_solution */
1324: for (i=0; i<ctx.physics.dof; i++) {
1325: DASetFieldName(ctx.da,i,ctx.physics.fieldname[i]);
1326: }
1327: /* Allow customization of the DA at runtime, mostly to change problem size with -da_grid_x M */
1328: DASetFromOptions(ctx.da);
1329: DAGetInfo(ctx.da,0, &Mx,0,0, 0,0,0, &dof,0,0,0);
1330: DAGetCorners(ctx.da,&xs,0,0,&xm,0,0);
1332: /* Set coordinates of cell centers */
1333: DASetUniformCoordinates(ctx.da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);
1335: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1336: PetscMalloc4(dof*dof,PetscScalar,&ctx.R,dof*dof,PetscScalar,&ctx.Rinv,2*dof,PetscScalar,&ctx.cjmpLR,1*dof,PetscScalar,&ctx.cslope);
1337: PetscMalloc2(2*dof,PetscScalar,&ctx.uLR,dof,PetscScalar,&ctx.flux);
1339: /* Create a vector to store the solution and to save the initial state */
1340: DACreateGlobalVector(ctx.da,&X);
1341: VecDuplicate(X,&X0);
1343: /* Create a time-stepping object */
1344: TSCreate(comm,&ts);
1345: TSSetProblemType(ts,TS_NONLINEAR);
1346: TSSetRHSFunction(ts,FVRHSFunction,&ctx);
1347: TSSetType(ts,TSSSP);
1348: TSSetDuration(ts,1000,10);
1350: /* Compute initial conditions and starting time step */
1351: FVSample(&ctx,0,X0);
1352: FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
1353: VecCopy(X0,X); /* The function value was not used so we set X=X0 again */
1354: TSSetInitialTimeStep(ts,0,ctx.cfl/ctx.cfl_idt);
1356: TSSetSolution(ts,X); /* The TS will use X for the solution, starting with it's current value as initial condition */
1357: TSSetFromOptions(ts); /* Take runtime options */
1359: SolutionStatsView(ctx.da,X,PETSC_VIEWER_STDOUT_WORLD);
1361: {
1362: PetscReal nrm1,nrmsup;
1363: PetscInt steps;
1365: TSStep(ts,&steps,&ptime);
1367: PetscPrintf(comm,"Final time %8.5f, steps %d\n",ptime,steps);
1368: if (ctx.exact) {
1369: SolutionErrorNorms(&ctx,ptime,X,&nrm1,&nrmsup);
1370: PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e ||x-x_e||_sup %8.4e\n",nrm1,nrmsup);
1371: }
1372: }
1374: SolutionStatsView(ctx.da,X,PETSC_VIEWER_STDOUT_WORLD);
1376: if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
1377: if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
1378: if (draw & 0x4) {
1379: Vec Y;
1380: VecDuplicate(X,&Y);
1381: FVSample(&ctx,ptime,Y);
1382: VecAYPX(Y,-1,X);
1383: VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1384: VecDestroy(Y);
1385: }
1387: if (view_final) {
1388: PetscViewer viewer;
1389: PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1390: PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1391: VecView(X,viewer);
1392: PetscViewerDestroy(viewer);
1393: }
1395: /* Clean up */
1396: (*ctx.physics.destroy)(ctx.physics.user);
1397: for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
1398: PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
1399: PetscFree2(ctx.uLR,ctx.flux);
1400: VecDestroy(X);
1401: VecDestroy(X0);
1402: DADestroy(ctx.da);
1403: TSDestroy(ts);
1404: PetscFinalize();
1405: return 0;
1406: }