Actual source code: ex9.c
petsc-3.3-p7 2013-05-11
1: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
2: "Solves scalar and vector problems, choose the physical model with -physics\n"
3: " advection - Constant coefficient scalar advection\n"
4: " u_t + (a*u)_x = 0\n"
5: " burgers - Burgers equation\n"
6: " u_t + (u^2/2)_x = 0\n"
7: " traffic - Traffic equation\n"
8: " u_t + (u*(1-u))_x = 0\n"
9: " acoustics - Acoustic wave propagation\n"
10: " u_t + (c*z*v)_x = 0\n"
11: " v_t + (c/z*u)_x = 0\n"
12: " isogas - Isothermal gas dynamics\n"
13: " rho_t + (rho*u)_x = 0\n"
14: " (rho*u)_t + (rho*u^2 + c^2*rho)_x = 0\n"
15: " shallow - Shallow water equations\n"
16: " h_t + (h*u)_x = 0\n"
17: " (h*u)_t + (h*u^2 + g*h^2/2)_x = 0\n"
18: "Some of these physical models have multiple Riemann solvers, select these with -physics_xxx_riemann\n"
19: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
20: " the states across shocks and rarefactions\n"
21: " roe - Linearized scheme, usually with an entropy fix inside sonic rarefactions\n"
22: "The systems provide a choice of reconstructions with -physics_xxx_reconstruct\n"
23: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
24: " conservative - Limit the conservative variables directly, can cause undesired interaction of waves\n\n"
25: "A variety of limiters for high-resolution TVD limiters are available with -limit\n"
26: " upwind,minmod,superbee,mc,vanleer,vanalbada,koren,cada-torillhon (last two are nominally third order)\n"
27: " and non-TVD schemes lax-wendroff,beam-warming,fromm\n\n"
28: "To preserve the TVD property, one should time step with a strong stability preserving method.\n"
29: "The optimal high order explicit Runge-Kutta methods in TSSSP are recommended for non-stiff problems.\n\n"
30: "Several initial conditions can be chosen with -initial N\n\n"
31: "The problem size should be set with -da_grid_x M\n\n";
33: /* To get isfinite in math.h */
34: #define _XOPEN_SOURCE 600
36: #include <petscts.h>
37: #include <petscdmda.h>
39: #include <../src/mat/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */
41: static inline PetscReal Sgn(PetscReal a) { return (a<0) ? -1 : 1; }
42: static inline PetscReal Abs(PetscReal a) { return (a<0) ? 0 : a; }
43: static inline PetscReal Sqr(PetscReal a) { return a*a; }
44: static inline PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; }
45: static inline PetscReal MinAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) < PetscAbs(b)) ? a : b; }
46: static inline PetscReal MinMod2(PetscReal a,PetscReal b)
47: { return (a*b<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscAbs(b)); }
48: static inline PetscReal MaxMod2(PetscReal a,PetscReal b)
49: { return (a*b<0) ? 0 : Sgn(a)*PetscMax(PetscAbs(a),PetscAbs(b)); }
50: static inline PetscReal MinMod3(PetscReal a,PetscReal b,PetscReal c)
51: {return (a*b<0 || a*c<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscMin(PetscAbs(b),PetscAbs(c))); }
53: static inline PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax)
54: { PetscReal range = xmax-xmin; return xmin + fmod(range+fmod(a,range),range); }
57: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
58: typedef struct _LimitInfo {
59: PetscReal hx;
60: PetscInt m;
61: } *LimitInfo;
62: static void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
63: {
64: PetscInt i;
65: for (i=0; i<info->m; i++) lmt[i] = 0;
66: }
67: static void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
68: {
69: PetscInt i;
70: for (i=0; i<info->m; i++) lmt[i] = jR[i];
71: }
72: static void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
73: {
74: PetscInt i;
75: for (i=0; i<info->m; i++) lmt[i] = jL[i];
76: }
77: static void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
78: {
79: PetscInt i;
80: for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i] + jR[i]);
81: }
82: static void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
83: {
84: PetscInt i;
85: for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]);
86: }
87: static void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
88: {
89: PetscInt i;
90: for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]));
91: }
92: static void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
93: {
94: PetscInt i;
95: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]);
96: }
97: static void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
98: { /* phi = (t + abs(t)) / (1 + abs(t)) */
99: PetscInt i;
100: for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Abs(jR[i]) + Abs(jL[i])*jR[i]) / (Abs(jL[i]) + Abs(jR[i]) + 1e-15);
101: }
102: static void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
103: { /* phi = (t + t^2) / (1 + t^2) */
104: PetscInt i;
105: for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
106: }
107: static void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
108: { /* phi = (t + t^2) / (1 + t^2) */
109: PetscInt i;
110: for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0
111: : (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
112: }
113: static void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
114: { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
115: PetscInt i;
116: for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i]) + 2*Sqr(jL[i])*jR[i])
117: / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
118: }
119: static void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
120: { /* Symmetric version of above */
121: PetscInt i;
122: for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i])
123: / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
124: }
125: static void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
126: { /* Eq 11 of Cada-Torrilhon 2009 */
127: PetscInt i;
128: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]);
129: }
131: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)
132: { return PetscMax(0,PetscMin((L+2*R)/3,
133: PetscMax(-0.5*L,
134: PetscMin(2*L,
135: PetscMin((L+2*R)/3,1.6*R)))));
136: }
137: static void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
138: { /* Cada-Torrilhon 2009, Eq 13 */
139: PetscInt i;
140: for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]);
141: }
142: static void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
143: { /* Cada-Torrilhon 2009, Eq 22 */
144: /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
145: const PetscReal eps = 1e-7,hx = info->hx;
146: PetscInt i;
147: for (i=0; i<info->m; i++) {
148: const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r*hx);
149: lmt[i] = ((eta < 1-eps)
150: ? (jL[i] + 2*jR[i]) / 3
151: : ((eta > 1+eps)
152: ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i])
153: : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3
154: + (1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]))));
155: }
156: }
157: static void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
158: { Limit_CadaTorrilhon3R(0.1,info,jL,jR,lmt); }
159: static void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
160: { Limit_CadaTorrilhon3R(1,info,jL,jR,lmt); }
161: static void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
162: { Limit_CadaTorrilhon3R(10,info,jL,jR,lmt); }
163: static void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
164: { Limit_CadaTorrilhon3R(100,info,jL,jR,lmt); }
167: /* --------------------------------- Finite Volume data structures ----------------------------------- */
169: typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType;
170: static const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0};
171: typedef PetscErrorCode (*RiemannFunction)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*);
172: typedef PetscErrorCode (*ReconstructFunction)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*,PetscReal*);
174: typedef struct {
175: PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*);
176: RiemannFunction riemann;
177: ReconstructFunction characteristic;
178: PetscErrorCode (*destroy)(void*);
179: void *user;
180: PetscInt dof;
181: char *fieldname[16];
182: } PhysicsCtx;
184: typedef struct {
185: void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*);
186: PhysicsCtx physics;
188: MPI_Comm comm;
189: char prefix[256];
191: /* Local work arrays */
192: PetscScalar *R,*Rinv; /* Characteristic basis, and it's inverse. COLUMN-MAJOR */
193: PetscScalar *cjmpLR; /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
194: PetscScalar *cslope; /* Limited slope, written in characteristic basis */
195: PetscScalar *uLR; /* Solution at left and right of interface, conservative variables, len=2*dof */
196: PetscScalar *flux; /* Flux across interface */
197: PetscReal *speeds; /* Speeds of each wave */
199: PetscReal cfl_idt; /* Max allowable value of 1/Delta t */
200: PetscReal cfl;
201: PetscReal xmin,xmax;
202: PetscInt initial;
203: PetscBool exact;
204: FVBCType bctype;
205: } FVCtx;
208: /* Utility */
212: PetscErrorCode RiemannListAdd(PetscFList *flist,const char *name,RiemannFunction rsolve)
213: {
217: PetscFListAdd(flist,name,"",(void(*)(void))rsolve);
218: return(0);
219: }
223: PetscErrorCode RiemannListFind(PetscFList flist,const char *name,RiemannFunction *rsolve)
224: {
228: PetscFListFind(flist,PETSC_COMM_WORLD,name,PETSC_FALSE,(void(**)(void))rsolve);
229: if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,1,"Riemann solver \"%s\" could not be found",name);
230: return(0);
231: }
235: PetscErrorCode ReconstructListAdd(PetscFList *flist,const char *name,ReconstructFunction r)
236: {
240: PetscFListAdd(flist,name,"",(void(*)(void))r);
241: return(0);
242: }
246: PetscErrorCode ReconstructListFind(PetscFList flist,const char *name,ReconstructFunction *r)
247: {
251: PetscFListFind(flist,PETSC_COMM_WORLD,name,PETSC_FALSE,(void(**)(void))r);
252: if (!*r) SETERRQ1(PETSC_COMM_SELF,1,"Reconstruction \"%s\" could not be found",name);
253: return(0);
254: }
257: /* --------------------------------- Physics ----------------------------------- */
258: /**
259: * Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction. These
260: * are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the
261: * number of fields and their names, and a function to deallocate private storage.
262: **/
264: /* First a few functions useful to several different physics */
267: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
268: {
269: PetscInt i,j;
272: for (i=0; i<m; i++) {
273: for (j=0; j<m; j++) {
274: Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j);
275: }
276: speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
277: }
278: return(0);
279: }
283: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
284: {
288: PetscFree(vctx);
289: return(0);
290: }
294: /* --------------------------------- Advection ----------------------------------- */
296: typedef struct {
297: PetscReal a; /* advective velocity */
298: } AdvectCtx;
302: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
303: {
304: AdvectCtx *ctx = (AdvectCtx*)vctx;
305: PetscReal speed;
308: speed = ctx->a;
309: flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
310: *maxspeed = speed;
311: return(0);
312: }
316: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
317: {
318: AdvectCtx *ctx = (AdvectCtx*)vctx;
321: X[0] = 1.;
322: Xi[0] = 1.;
323: speeds[0] = ctx->a;
324: return(0);
325: }
329: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
330: {
331: AdvectCtx *ctx = (AdvectCtx*)vctx;
332: PetscReal a = ctx->a,x0;
335: switch (bctype) {
336: case FVBC_OUTFLOW: x0 = x-a*t; break;
337: case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
338: default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
339: }
340: switch (initial) {
341: case 0: u[0] = (x0 < 0) ? 1 : -1; break;
342: case 1: u[0] = (x0 < 0) ? -1 : 1; break;
343: case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
344: case 3: u[0] = sin(2*M_PI*x0); break;
345: case 4: u[0] = PetscAbs(x0); break;
346: case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(sin(2*M_PI*x0)); break;
347: case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
348: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
349: }
350: return(0);
351: }
355: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
356: {
358: AdvectCtx *user;
361: PetscNew(*user,&user);
362: ctx->physics.sample = PhysicsSample_Advect;
363: ctx->physics.riemann = PhysicsRiemann_Advect;
364: ctx->physics.characteristic = PhysicsCharacteristic_Advect;
365: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
366: ctx->physics.user = user;
367: ctx->physics.dof = 1;
368: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
369: user->a = 1;
370: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
371: {
372: PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,PETSC_NULL);
373: }
374: PetscOptionsEnd();
375: return(0);
376: }
380: /* --------------------------------- Burgers ----------------------------------- */
382: typedef struct {
383: PetscReal lxf_speed;
384: } BurgersCtx;
388: static PetscErrorCode PhysicsSample_Burgers(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
389: {
392: if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
393: switch (initial) {
394: case 0: u[0] = (x < 0) ? 1 : -1; break;
395: case 1:
396: if (x < -t) u[0] = -1;
397: else if (x < t) u[0] = x/t;
398: else u[0] = 1;
399: break;
400: case 2:
401: if (x < 0) u[0] = 0;
402: else if (x <= t) u[0] = x/t;
403: else if (x < 1+0.5*t) u[0] = 1;
404: else u[0] = 0;
405: break;
406: case 3:
407: if (x < 0.2*t) u[0] = 0.2;
408: else if (x < t) u[0] = x/t;
409: else u[0] = 1;
410: break;
411: case 4:
412: if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
413: u[0] = 0.7 + 0.3*sin(2*M_PI*((x-xmin)/(xmax-xmin)));
414: break;
415: case 5: /* Pure shock solution */
416: if (x < 0.5*t) u[0] = 1;
417: else u[0] = 0;
418: break;
419: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
420: }
421: return(0);
422: }
426: static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
427: {
430: if (uL[0] < uR[0]) { /* rarefaction */
431: flux[0] = (uL[0]*uR[0] < 0)
432: ? 0 /* sonic rarefaction */
433: : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0]));
434: } else { /* shock */
435: flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0]));
436: }
437: *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
438: return(0);
439: }
443: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
444: {
445: PetscReal speed;
448: speed = 0.5*(uL[0] + uR[0]);
449: flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
450: if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
451: *maxspeed = speed;
452: return(0);
453: }
457: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
458: {
459: PetscReal c;
460: PetscScalar fL,fR;
463: c = ((BurgersCtx*)vctx)->lxf_speed;
464: fL = 0.5*PetscSqr(uL[0]);
465: fR = 0.5*PetscSqr(uR[0]);
466: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
467: *maxspeed = c;
468: return(0);
469: }
473: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
474: {
475: PetscReal c;
476: PetscScalar fL,fR;
479: c = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0]));
480: fL = 0.5*PetscSqr(uL[0]);
481: fR = 0.5*PetscSqr(uR[0]);
482: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
483: *maxspeed = c;
484: return(0);
485: }
489: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
490: {
491: BurgersCtx *user;
493: RiemannFunction r;
494: PetscFList rlist = 0;
495: char rname[256] = "exact";
498: PetscNew(*user,&user);
499: ctx->physics.sample = PhysicsSample_Burgers;
500: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
501: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
502: ctx->physics.user = user;
503: ctx->physics.dof = 1;
504: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
505: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Burgers_Exact);
506: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Burgers_Roe);
507: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Burgers_LxF);
508: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Burgers_Rusanov);
509: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
510: {
511: PetscOptionsList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
512: }
513: PetscOptionsEnd();
514: RiemannListFind(rlist,rname,&r);
515: PetscFListDestroy(&rlist);
516: ctx->physics.riemann = r;
518: /* *
519: * Hack to deal with LxF in semi-discrete form
520: * max speed is 1 for the basic initial conditions (where |u| <= 1)
521: * */
522: if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
523: return(0);
524: }
528: /* --------------------------------- Traffic ----------------------------------- */
530: typedef struct {
531: PetscReal lxf_speed;
532: PetscReal a;
533: } TrafficCtx;
535: static inline PetscScalar TrafficFlux(PetscScalar a,PetscScalar u) { return a*u*(1-u); }
539: static PetscErrorCode PhysicsSample_Traffic(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
540: {
541: PetscReal a = ((TrafficCtx*)vctx)->a;
544: if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
545: switch (initial) {
546: case 0:
547: u[0] = (-a*t < x) ? 2 : 0; break;
548: case 1:
549: if (x < PetscMin(2*a*t,0.5+a*t)) u[0] = -1;
550: else if (x < 1) u[0] = 0;
551: else u[0] = 1;
552: break;
553: case 2:
554: if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
555: u[0] = 0.7 + 0.3*sin(2*M_PI*((x-xmin)/(xmax-xmin)));
556: break;
557: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
558: }
559: return(0);
560: }
564: static PetscErrorCode PhysicsRiemann_Traffic_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
565: {
566: PetscReal a = ((TrafficCtx*)vctx)->a;
569: if (uL[0] < uR[0]) {
570: flux[0] = PetscMin(TrafficFlux(a,uL[0]),
571: TrafficFlux(a,uR[0]));
572: } else {
573: flux[0] = (uR[0] < 0.5 && 0.5 < uL[0])
574: ? TrafficFlux(a,0.5)
575: : PetscMax(TrafficFlux(a,uL[0]),
576: TrafficFlux(a,uR[0]));
577: }
578: *maxspeed = a*MaxAbs(1-2*uL[0],1-2*uR[0]);
579: return(0);
580: }
584: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
585: {
586: PetscReal a = ((TrafficCtx*)vctx)->a;
587: PetscReal speed;
590: speed = a*(1 - (uL[0] + uR[0]));
591: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
592: *maxspeed = speed;
593: return(0);
594: }
598: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
599: {
600: TrafficCtx *phys = (TrafficCtx*)vctx;
601: PetscReal a = phys->a;
602: PetscReal speed;
605: speed = a*(1 - (uL[0] + uR[0]));
606: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]);
607: *maxspeed = speed;
608: return(0);
609: }
613: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
614: {
615: PetscReal a = ((TrafficCtx*)vctx)->a;
616: PetscReal speed;
619: speed = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0]));
620: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]);
621: *maxspeed = speed;
622: return(0);
623: }
627: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
628: {
630: TrafficCtx *user;
631: RiemannFunction r;
632: PetscFList rlist = 0;
633: char rname[256] = "exact";
636: PetscNew(*user,&user);
637: ctx->physics.sample = PhysicsSample_Traffic;
638: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
639: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
640: ctx->physics.user = user;
641: ctx->physics.dof = 1;
642: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
643: user->a = 0.5;
644: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Traffic_Exact);
645: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Traffic_Roe);
646: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Traffic_LxF);
647: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Traffic_Rusanov);
648: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");
649: {
650: PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,PETSC_NULL);
651: PetscOptionsList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
652: }
653: PetscOptionsEnd();
655: RiemannListFind(rlist,rname,&r);
656: PetscFListDestroy(&rlist);
657: ctx->physics.riemann = r;
659: /* *
660: * Hack to deal with LxF in semi-discrete form
661: * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
662: * */
663: if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3*user->a;
665: return(0);
666: }
668: /* --------------------------------- Linear Acoustics ----------------------------------- */
670: /* Flux: u_t + (A u)_x
671: * z = sqrt(rho*bulk), c = sqrt(rho/bulk)
672: * Spectral decomposition: A = R * D * Rinv
673: * [ cz] = [-z z] [-c ] [-1/2z 1/2]
674: * [c/z ] = [ 1 1] [ c] [ 1/2z 1/2]
675: *
676: * We decompose this into the left-traveling waves Al = R * D^- Rinv
677: * and the right-traveling waves Ar = R * D^+ * Rinv
678: * Multiplying out these expressions produces the following two matrices
679: */
681: typedef struct {
682: PetscReal c; /* speed of sound: c = sqrt(bulk/rho) */
683: PetscReal z; /* impedence: z = sqrt(rho*bulk) */
684: } AcousticsCtx;
686: static inline void AcousticsFlux(AcousticsCtx *ctx,const PetscScalar *u,PetscScalar *f)
687: {
688: f[0] = ctx->c*ctx->z*u[1];
689: f[1] = ctx->c/ctx->z*u[0];
690: }
694: static PetscErrorCode PhysicsCharacteristic_Acoustics(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
695: {
696: AcousticsCtx *phys = (AcousticsCtx*)vctx;
697: PetscReal z = phys->z,c = phys->c;
700: X[0*2+0] = -z;
701: X[0*2+1] = z;
702: X[1*2+0] = 1;
703: X[1*2+1] = 1;
704: Xi[0*2+0] = -1./(2*z);
705: Xi[0*2+1] = 1./2;
706: Xi[1*2+0] = 1./(2*z);
707: Xi[1*2+1] = 1./2;
708: speeds[0] = -c;
709: speeds[1] = c;
710: return(0);
711: }
715: static PetscErrorCode PhysicsSample_Acoustics_Initial(AcousticsCtx *phys,PetscInt initial,PetscReal xmin,PetscReal xmax,PetscReal x,PetscReal *u)
716: {
718: switch (initial) {
719: case 0:
720: u[0] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.2) < 0.1) ? 1 : 0.5;
721: u[1] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.7) < 0.1) ? 1 : -0.5;
722: break;
723: case 1:
724: u[0] = cos(3 * 2*PETSC_PI*x/(xmax-xmin));
725: u[1] = exp(-PetscSqr(x - (xmax + xmin)/2) / (2*PetscSqr(0.2*(xmax - xmin)))) - 0.5;
726: break;
727: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
728: }
729: return(0);
730: }
734: static PetscErrorCode PhysicsSample_Acoustics(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
735: {
736: AcousticsCtx *phys = (AcousticsCtx*)vctx;
737: PetscReal c = phys->c;
738: PetscReal x0a,x0b,u0a[2],u0b[2],tmp[2];
739: PetscReal X[2][2],Xi[2][2],dummy[2];
743: switch (bctype) {
744: case FVBC_OUTFLOW:
745: x0a = x+c*t;
746: x0b = x-c*t;
747: break;
748: case FVBC_PERIODIC:
749: x0a = RangeMod(x+c*t,xmin,xmax);
750: x0b = RangeMod(x-c*t,xmin,xmax);
751: break;
752: default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
753: }
754: PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0a,u0a);
755: PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0b,u0b);
756: PhysicsCharacteristic_Acoustics(vctx,2,u,&X[0][0],&Xi[0][0],dummy);
757: tmp[0] = Xi[0][0]*u0a[0] + Xi[0][1]*u0a[1];
758: tmp[1] = Xi[1][0]*u0b[0] + Xi[1][1]*u0b[1];
759: u[0] = X[0][0]*tmp[0] + X[0][1]*tmp[1];
760: u[1] = X[1][0]*tmp[0] + X[1][1]*tmp[1];
761: return(0);
762: }
766: static PetscErrorCode PhysicsRiemann_Acoustics_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
767: {
768: AcousticsCtx *phys = (AcousticsCtx*)vctx;
769: PetscReal c = phys->c,z = phys->z;
770: PetscReal
771: Al[2][2] = {{-c/2 , c*z/2 },
772: {c/(2*z) , -c/2 }}, /* Left traveling waves */
773: Ar[2][2] = {{c/2 , c*z/2 },
774: {c/(2*z) , c/2 }}; /* Right traveling waves */
777: flux[0] = Al[0][0]*uR[0] + Al[0][1]*uR[1] + Ar[0][0]*uL[0] + Ar[0][1]*uL[1];
778: flux[1] = Al[1][0]*uR[0] + Al[1][1]*uR[1] + Ar[1][0]*uL[0] + Ar[1][1]*uL[1];
779: *maxspeed = c;
780: return(0);
781: }
785: static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx)
786: {
788: AcousticsCtx *user;
789: PetscFList rlist = 0,rclist = 0;
790: char rname[256] = "exact",rcname[256] = "characteristic";
793: PetscNew(*user,&user);
794: ctx->physics.sample = PhysicsSample_Acoustics;
795: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
796: ctx->physics.user = user;
797: ctx->physics.dof = 2;
798: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
799: PetscStrallocpy("v",&ctx->physics.fieldname[1]);
800: user->c = 1;
801: user->z = 1;
802: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Acoustics_Exact);
803: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Acoustics);
804: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
805: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for linear Acoustics","");
806: {
807: PetscOptionsReal("-physics_acoustics_c","c = sqrt(bulk/rho)","",user->c,&user->c,PETSC_NULL);
808: PetscOptionsReal("-physics_acoustics_z","z = sqrt(bulk*rho)","",user->z,&user->z,PETSC_NULL);
809: PetscOptionsList("-physics_acoustics_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
810: PetscOptionsList("-physics_acoustics_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
811: }
812: PetscOptionsEnd();
813: RiemannListFind(rlist,rname,&ctx->physics.riemann);
814: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
815: PetscFListDestroy(&rlist);
816: PetscFListDestroy(&rclist);
817: return(0);
818: }
820: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */
822: typedef struct {
823: PetscReal acoustic_speed;
824: } IsoGasCtx;
826: static inline void IsoGasFlux(PetscReal c,const PetscScalar *u,PetscScalar *f)
827: {
828: f[0] = u[1];
829: f[1] = PetscSqr(u[1])/u[0] + c*c*u[0];
830: }
834: static PetscErrorCode PhysicsSample_IsoGas(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
835: {
838: if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solutions not implemented for t > 0");
839: switch (initial) {
840: case 0:
841: u[0] = (x < 0) ? 1 : 0.5;
842: u[1] = (x < 0) ? 1 : 0.7;
843: break;
844: case 1:
845: u[0] = 1+0.5*sin(2*M_PI*x);
846: u[1] = 1*u[0];
847: break;
848: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
849: }
850: return(0);
851: }
855: static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
856: {
857: IsoGasCtx *phys = (IsoGasCtx*)vctx;
858: PetscReal c = phys->acoustic_speed;
859: PetscScalar ubar,du[2],a[2],fL[2],fR[2],lam[2],ustar[2],R[2][2];
860: PetscInt i;
863: ubar = (uL[1]/PetscSqrtScalar(uL[0]) + uR[1]/PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0]));
864: /* write fluxuations in characteristic basis */
865: du[0] = uR[0] - uL[0];
866: du[1] = uR[1] - uL[1];
867: a[0] = (1/(2*c)) * ((ubar + c)*du[0] - du[1]);
868: a[1] = (1/(2*c)) * ((-ubar + c)*du[0] + du[1]);
869: /* wave speeds */
870: lam[0] = ubar - c;
871: lam[1] = ubar + c;
872: /* Right eigenvectors */
873: R[0][0] = 1; R[0][1] = ubar-c;
874: R[1][0] = 1; R[1][1] = ubar+c;
875: /* Compute state in star region (between the 1-wave and 2-wave) */
876: for (i=0; i<2; i++) ustar[i] = uL[i] + a[0]*R[0][i];
877: if (uL[1]/uL[0] < c && c < ustar[1]/ustar[0]) { /* 1-wave is sonic rarefaction */
878: PetscScalar ufan[2];
879: ufan[0] = uL[0]*PetscExpScalar(uL[1]/(uL[0]*c) - 1);
880: ufan[1] = c*ufan[0];
881: IsoGasFlux(c,ufan,flux);
882: } else if (ustar[1]/ustar[0] < -c && -c < uR[1]/uR[0]) { /* 2-wave is sonic rarefaction */
883: PetscScalar ufan[2];
884: ufan[0] = uR[0]*PetscExpScalar(-uR[1]/(uR[0]*c) - 1);
885: ufan[1] = -c*ufan[0];
886: IsoGasFlux(c,ufan,flux);
887: } else { /* Centered form */
888: IsoGasFlux(c,uL,fL);
889: IsoGasFlux(c,uR,fR);
890: for (i=0; i<2; i++) {
891: PetscScalar absdu = PetscAbsScalar(lam[0])*a[0]*R[0][i] + PetscAbsScalar(lam[1])*a[1]*R[1][i];
892: flux[i] = 0.5*(fL[i]+fR[i]) - 0.5*absdu;
893: }
894: }
895: *maxspeed = MaxAbs(lam[0],lam[1]);
896: return(0);
897: }
901: static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
902: {
903: IsoGasCtx *phys = (IsoGasCtx*)vctx;
904: PetscReal c = phys->acoustic_speed;
905: PetscScalar ustar[2];
906: struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
907: PetscInt i;
910: if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
911: {
912: /* Solve for star state */
913: PetscScalar res,tmp,rho = 0.5*(L.rho + R.rho); /* initial guess */
914: for (i=0; i<20; i++) {
915: PetscScalar fr,fl,dfr,dfl;
916: fl = (L.rho < rho)
917: ? (rho-L.rho)/PetscSqrtScalar(L.rho*rho) /* shock */
918: : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
919: fr = (R.rho < rho)
920: ? (rho-R.rho)/PetscSqrtScalar(R.rho*rho) /* shock */
921: : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
922: res = R.u-L.u + c*(fr+fl);
923: if (!isfinite(res)) SETERRQ1(PETSC_COMM_SELF,1,"non-finite residual=%g",res);
924: if (PetscAbsScalar(res) < 1e-10) {
925: star.rho = rho;
926: star.u = L.u - c*fl;
927: goto converged;
928: }
929: dfl = (L.rho < rho)
930: ? 1/PetscSqrtScalar(L.rho*rho)*(1 - 0.5*(rho-L.rho)/rho)
931: : 1/rho;
932: dfr = (R.rho < rho)
933: ? 1/PetscSqrtScalar(R.rho*rho)*(1 - 0.5*(rho-R.rho)/rho)
934: : 1/rho;
935: tmp = rho - res/(c*(dfr+dfl));
936: if (tmp <= 0) rho /= 2; /* Guard against Newton shooting off to a negative density */
937: else rho = tmp;
938: if (!((rho > 0) && isnormal(rho))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate rho=%g",rho);
939: }
940: SETERRQ1(PETSC_COMM_SELF,1,"Newton iteration for star.rho diverged after %d iterations",i);
941: }
942: converged:
943: if (L.u-c < 0 && 0 < star.u-c) { /* 1-wave is sonic rarefaction */
944: PetscScalar ufan[2];
945: ufan[0] = L.rho*PetscExpScalar(L.u/c - 1);
946: ufan[1] = c*ufan[0];
947: IsoGasFlux(c,ufan,flux);
948: } else if (star.u+c < 0 && 0 < R.u+c) { /* 2-wave is sonic rarefaction */
949: PetscScalar ufan[2];
950: ufan[0] = R.rho*PetscExpScalar(-R.u/c - 1);
951: ufan[1] = -c*ufan[0];
952: IsoGasFlux(c,ufan,flux);
953: } else if ((L.rho >= star.rho && L.u-c >= 0)
954: || (L.rho < star.rho && (star.rho*star.u-L.rho*L.u)/(star.rho-L.rho) > 0)) {
955: /* 1-wave is supersonic rarefaction, or supersonic shock */
956: IsoGasFlux(c,uL,flux);
957: } else if ((star.rho <= R.rho && R.u+c <= 0)
958: || (star.rho > R.rho && (R.rho*R.u-star.rho*star.u)/(R.rho-star.rho) < 0)) {
959: /* 2-wave is supersonic rarefaction or supersonic shock */
960: IsoGasFlux(c,uR,flux);
961: } else {
962: ustar[0] = star.rho;
963: ustar[1] = star.rho*star.u;
964: IsoGasFlux(c,ustar,flux);
965: }
966: *maxspeed = MaxAbs(MaxAbs(star.u-c,star.u+c),MaxAbs(L.u-c,R.u+c));
967: return(0);
968: }
972: static PetscErrorCode PhysicsRiemann_IsoGas_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
973: {
974: IsoGasCtx *phys = (IsoGasCtx*)vctx;
975: PetscScalar c = phys->acoustic_speed,fL[2],fR[2],s;
976: struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
979: if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
980: IsoGasFlux(c,uL,fL);
981: IsoGasFlux(c,uR,fR);
982: s = PetscMax(PetscAbs(L.u),PetscAbs(R.u))+c;
983: flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
984: flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
985: *maxspeed = s;
986: return(0);
987: }
991: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
992: {
993: IsoGasCtx *phys = (IsoGasCtx*)vctx;
994: PetscReal c = phys->acoustic_speed;
998: speeds[0] = u[1]/u[0] - c;
999: speeds[1] = u[1]/u[0] + c;
1000: X[0*2+0] = 1;
1001: X[0*2+1] = speeds[0];
1002: X[1*2+0] = 1;
1003: X[1*2+1] = speeds[1];
1004: PetscMemcpy(Xi,X,4*sizeof(X[0]));
1005: PetscKernel_A_gets_inverse_A_2(Xi,0);
1006: return(0);
1007: }
1011: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
1012: {
1014: IsoGasCtx *user;
1015: PetscFList rlist = 0,rclist = 0;
1016: char rname[256] = "exact",rcname[256] = "characteristic";
1019: PetscNew(*user,&user);
1020: ctx->physics.sample = PhysicsSample_IsoGas;
1021: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
1022: ctx->physics.user = user;
1023: ctx->physics.dof = 2;
1024: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1025: PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1026: user->acoustic_speed = 1;
1027: RiemannListAdd(&rlist,"exact", PhysicsRiemann_IsoGas_Exact);
1028: RiemannListAdd(&rlist,"roe", PhysicsRiemann_IsoGas_Roe);
1029: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_IsoGas_Rusanov);
1030: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_IsoGas);
1031: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1032: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");
1033: {
1034: PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,PETSC_NULL);
1035: PetscOptionsList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
1036: PetscOptionsList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
1037: }
1038: PetscOptionsEnd();
1039: RiemannListFind(rlist,rname,&ctx->physics.riemann);
1040: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1041: PetscFListDestroy(&rlist);
1042: PetscFListDestroy(&rclist);
1043: return(0);
1044: }
1048: /* --------------------------------- Shallow Water ----------------------------------- */
1050: typedef struct {
1051: PetscReal gravity;
1052: } ShallowCtx;
1054: static inline void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
1055: {
1056: f[0] = u[1];
1057: f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
1058: }
1062: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
1063: {
1064: ShallowCtx *phys = (ShallowCtx*)vctx;
1065: PetscScalar g = phys->gravity,ustar[2],cL,cR,c,cstar;
1066: struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
1067: PetscInt i;
1070: if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed thickness is negative");
1071: cL = PetscSqrtScalar(g*L.h);
1072: cR = PetscSqrtScalar(g*R.h);
1073: c = PetscMax(cL,cR);
1074: {
1075: /* Solve for star state */
1076: const PetscInt maxits = 50;
1077: PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */
1078: h0 = h;
1079: for (i=0; i<maxits; i++) {
1080: PetscScalar fr,fl,dfr,dfl;
1081: fl = (L.h < h)
1082: ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */
1083: : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h); /* rarefaction */
1084: fr = (R.h < h)
1085: ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */
1086: : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h); /* rarefaction */
1087: res = R.u - L.u + fr + fl;
1088: if (!isfinite(res)) SETERRQ1(PETSC_COMM_SELF,1,"non-finite residual=%g",res);
1089: if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h-h0) < 1e-8)) {
1090: star.h = h;
1091: star.u = L.u - fl;
1092: goto converged;
1093: } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */
1094: h = 0.8*h0 + 0.2*h;
1095: continue;
1096: }
1097: /* Accept the last step and take another */
1098: res0 = res;
1099: h0 = h;
1100: dfl = (L.h < h)
1101: ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h)
1102: : PetscSqrtScalar(g/h);
1103: dfr = (R.h < h)
1104: ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h)
1105: : PetscSqrtScalar(g/h);
1106: tmp = h - res/(dfr+dfl);
1107: if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */
1108: else h = tmp;
1109: if (!((h > 0) && isnormal(h))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate h=%g",h);
1110: }
1111: SETERRQ1(PETSC_COMM_SELF,1,"Newton iteration for star.h diverged after %d iterations",i);
1112: }
1113: converged:
1114: cstar = PetscSqrtScalar(g*star.h);
1115: if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */
1116: PetscScalar ufan[2];
1117: ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL);
1118: ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0];
1119: ShallowFlux(phys,ufan,flux);
1120: } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */
1121: PetscScalar ufan[2];
1122: ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR);
1123: ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0];
1124: ShallowFlux(phys,ufan,flux);
1125: } else if ((L.h >= star.h && L.u-c >= 0)
1126: || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) {
1127: /* 1-wave is right-travelling shock (supersonic) */
1128: ShallowFlux(phys,uL,flux);
1129: } else if ((star.h <= R.h && R.u+c <= 0)
1130: || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) {
1131: /* 2-wave is left-travelling shock (supersonic) */
1132: ShallowFlux(phys,uR,flux);
1133: } else {
1134: ustar[0] = star.h;
1135: ustar[1] = star.h*star.u;
1136: ShallowFlux(phys,ustar,flux);
1137: }
1138: *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR));
1139: return(0);
1140: }
1144: static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
1145: {
1146: ShallowCtx *phys = (ShallowCtx*)vctx;
1147: PetscScalar g = phys->gravity,fL[2],fR[2],s;
1148: struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
1151: if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed thickness is negative");
1152: ShallowFlux(phys,uL,fL);
1153: ShallowFlux(phys,uR,fR);
1154: s = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h));
1155: flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
1156: flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
1157: *maxspeed = s;
1158: return(0);
1159: }
1163: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
1164: {
1165: ShallowCtx *phys = (ShallowCtx*)vctx;
1166: PetscReal c;
1170: c = PetscSqrtScalar(u[0]*phys->gravity);
1171: speeds[0] = u[1]/u[0] - c;
1172: speeds[1] = u[1]/u[0] + c;
1173: X[0*2+0] = 1;
1174: X[0*2+1] = speeds[0];
1175: X[1*2+0] = 1;
1176: X[1*2+1] = speeds[1];
1177: PetscMemcpy(Xi,X,4*sizeof(X[0]));
1178: PetscKernel_A_gets_inverse_A_2(Xi,0);
1179: return(0);
1180: }
1184: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
1185: {
1187: ShallowCtx *user;
1188: PetscFList rlist = 0,rclist = 0;
1189: char rname[256] = "exact",rcname[256] = "characteristic";
1192: PetscNew(*user,&user);
1193: /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
1194: ctx->physics.sample = PhysicsSample_IsoGas;
1195: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
1196: ctx->physics.user = user;
1197: ctx->physics.dof = 2;
1198: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1199: PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1200: user->gravity = 1;
1201: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Shallow_Exact);
1202: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov);
1203: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Shallow);
1204: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1205: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");
1206: {
1207: PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,PETSC_NULL);
1208: PetscOptionsList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof rname,PETSC_NULL);
1209: PetscOptionsList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof rcname,PETSC_NULL);
1210: }
1211: PetscOptionsEnd();
1212: RiemannListFind(rlist,rname,&ctx->physics.riemann);
1213: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1214: PetscFListDestroy(&rlist);
1215: PetscFListDestroy(&rclist);
1216: return(0);
1217: }
1221: /* --------------------------------- Finite Volume Solver ----------------------------------- */
1225: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
1226: {
1227: FVCtx *ctx = (FVCtx*)vctx;
1228: PetscErrorCode ierr;
1229: PetscInt i,j,k,Mx,dof,xs,xm;
1230: PetscReal hx,cfl_idt = 0;
1231: PetscScalar *x,*f,*slope;
1232: Vec Xloc;
1233: DM da;
1236: TSGetDM(ts,&da);
1237: DMGetLocalVector(da,&Xloc);
1238: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1239: hx = (ctx->xmax - ctx->xmin)/Mx;
1240: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1241: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
1243: VecZeroEntries(F);
1245: DMDAVecGetArray(da,Xloc,&x);
1246: DMDAVecGetArray(da,F,&f);
1247: DMDAGetArray(da,PETSC_TRUE,&slope);
1249: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1251: if (ctx->bctype == FVBC_OUTFLOW) {
1252: for (i=xs-2; i<0; i++) {
1253: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
1254: }
1255: for (i=Mx; i<xs+xm+2; i++) {
1256: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
1257: }
1258: }
1259: for (i=xs-1; i<xs+xm+1; i++) {
1260: struct _LimitInfo info;
1261: PetscScalar *cjmpL,*cjmpR;
1262: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
1263: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1264: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
1265: PetscMemzero(ctx->cjmpLR,2*dof*sizeof(ctx->cjmpLR[0]));
1266: cjmpL = &ctx->cjmpLR[0];
1267: cjmpR = &ctx->cjmpLR[dof];
1268: for (j=0; j<dof; j++) {
1269: PetscScalar jmpL,jmpR;
1270: jmpL = x[(i+0)*dof+j] - x[(i-1)*dof+j];
1271: jmpR = x[(i+1)*dof+j] - x[(i+0)*dof+j];
1272: for (k=0; k<dof; k++) {
1273: cjmpL[k] += ctx->Rinv[k+j*dof] * jmpL;
1274: cjmpR[k] += ctx->Rinv[k+j*dof] * jmpR;
1275: }
1276: }
1277: /* Apply limiter to the left and right characteristic jumps */
1278: info.m = dof;
1279: info.hx = hx;
1280: (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
1281: for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
1282: for (j=0; j<dof; j++) {
1283: PetscScalar tmp = 0;
1284: for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof] * ctx->cslope[k];
1285: slope[i*dof+j] = tmp;
1286: }
1287: }
1289: for (i=xs; i<xs+xm+1; i++) {
1290: PetscReal maxspeed;
1291: PetscScalar *uL,*uR;
1292: uL = &ctx->uLR[0];
1293: uR = &ctx->uLR[dof];
1294: for (j=0; j<dof; j++) {
1295: uL[j] = x[(i-1)*dof+j] + slope[(i-1)*dof+j]*hx/2;
1296: uR[j] = x[(i-0)*dof+j] - slope[(i-0)*dof+j]*hx/2;
1297: }
1298: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed);
1299: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
1301: if (i > xs) {
1302: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
1303: }
1304: if (i < xs+xm) {
1305: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
1306: }
1307: }
1309: DMDAVecRestoreArray(da,Xloc,&x);
1310: DMDAVecRestoreArray(da,F,&f);
1311: DMDARestoreArray(da,PETSC_TRUE,&slope);
1312: DMRestoreLocalVector(da,&Xloc);
1314: MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,((PetscObject)da)->comm);
1315: if (0) {
1316: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
1317: PetscReal dt,tnow;
1318: TSGetTimeStep(ts,&dt);
1319: TSGetTime(ts,&tnow);
1320: if (dt > 0.5/ctx->cfl_idt) {
1321: if (1) {
1322: PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",tnow,dt,0.5/ctx->cfl_idt);
1323: } else {
1324: SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",dt,ctx->cfl/ctx->cfl_idt);
1325: }
1326: }
1327: }
1328: return(0);
1329: }
1333: static PetscErrorCode SmallMatMultADB(PetscScalar *C,PetscInt bs,const PetscScalar *A,const PetscReal *D,const PetscScalar *B)
1334: {
1335: PetscInt i,j,k;
1338: for (i=0; i<bs; i++) {
1339: for (j=0; j<bs; j++) {
1340: PetscScalar tmp = 0;
1341: for (k=0; k<bs; k++) {
1342: tmp += A[i*bs+k] * D[k] * B[k*bs+j];
1343: }
1344: C[i*bs+j] = tmp;
1345: }
1346: }
1347: return(0);
1348: }
1353: static PetscErrorCode FVIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *vctx)
1354: {
1355: FVCtx *ctx = (FVCtx*)vctx;
1357: PetscInt i,j,dof = ctx->physics.dof;
1358: PetscScalar *x,*J;
1359: PetscReal hx;
1360: DM da;
1361: DMDALocalInfo dainfo;
1364: TSGetDM(ts,&da);
1365: DMDAVecGetArray(da,X,&x);
1366: DMDAGetLocalInfo(da,&dainfo);
1367: hx = (ctx->xmax - ctx->xmin)/dainfo.mx;
1368: PetscMalloc(dof*dof*sizeof(PetscScalar),&J);
1369: for (i=dainfo.xs; i<dainfo.xs+dainfo.xm; i++) {
1370: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1371: for (j=0; j<dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]);
1372: SmallMatMultADB(J,dof,ctx->R,ctx->speeds,ctx->Rinv);
1373: for (j=0; j<dof*dof; j++) J[j] = J[j]/hx + shift*(j/dof == j%dof);
1374: MatSetValuesBlocked(*B,1,&i,1,&i,J,INSERT_VALUES);
1375: }
1376: PetscFree(J);
1377: DMDAVecRestoreArray(da,X,&x);
1379: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1380: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1381: if (*A != *B) {
1382: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1383: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1384: }
1385: return(0);
1386: }
1390: static PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
1391: {
1393: PetscScalar *u,*uj;
1394: PetscInt i,j,k,dof,xs,xm,Mx;
1397: if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,1,"Physics has not provided a sampling function");
1398: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1399: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1400: DMDAVecGetArray(da,U,&u);
1401: PetscMalloc(dof*sizeof uj[0],&uj);
1402: for (i=xs; i<xs+xm; i++) {
1403: const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h;
1404: const PetscInt N = 200;
1405: /* Integrate over cell i using trapezoid rule with N points. */
1406: for (k=0; k<dof; k++) u[i*dof+k] = 0;
1407: for (j=0; j<N+1; j++) {
1408: PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N;
1409: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
1410: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N)?0.5:1.0)*uj[k]/N;
1411: }
1412: }
1413: DMDAVecRestoreArray(da,U,&u);
1414: PetscFree(uj);
1415: return(0);
1416: }
1420: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
1421: {
1423: PetscReal xmin,xmax;
1424: PetscScalar sum,*x,tvsum,tvgsum;
1425: PetscInt imin,imax,Mx,i,j,xs,xm,dof;
1426: Vec Xloc;
1427: PetscBool iascii;
1430: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1431: if (iascii) {
1432: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
1433: DMGetLocalVector(da,&Xloc);
1434: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1435: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
1436: DMDAVecGetArray(da,Xloc,&x);
1437: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1438: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1439: tvsum = 0;
1440: for (i=xs; i<xs+xm; i++) {
1441: for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j] - x[(i-1)*dof+j]);
1442: }
1443: MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_MAX,((PetscObject)da)->comm);
1444: DMDAVecRestoreArray(da,Xloc,&x);
1445: DMRestoreLocalVector(da,&Xloc);
1447: VecMin(X,&imin,&xmin);
1448: VecMax(X,&imax,&xmax);
1449: VecSum(X,&sum);
1450: PetscViewerASCIIPrintf(viewer,"Solution range [%8.5f,%8.5f] with extrema at %d and %d, mean %8.5f, ||x||_TV %8.5f\n",xmin,xmax,imin,imax,sum/Mx,tvgsum/Mx);
1451: } else {
1452: SETERRQ(PETSC_COMM_SELF,1,"Viewer type not supported");
1453: }
1454: return(0);
1455: }
1459: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1,PetscReal *nrmsup)
1460: {
1462: Vec Y;
1463: PetscInt Mx;
1466: VecGetSize(X,&Mx);
1467: VecDuplicate(X,&Y);
1468: FVSample(ctx,da,t,Y);
1469: VecAYPX(Y,-1,X);
1470: VecNorm(Y,NORM_1,nrm1);
1471: VecNorm(Y,NORM_INFINITY,nrmsup);
1472: *nrm1 /= Mx;
1473: VecDestroy(&Y);
1474: return(0);
1475: }
1479: int main(int argc,char *argv[])
1480: {
1481: char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
1482: PetscFList limiters = 0,physics = 0;
1483: MPI_Comm comm;
1484: TS ts;
1485: DM da;
1486: Vec X,X0,R;
1487: Mat B;
1488: FVCtx ctx;
1489: PetscInt i,dof,xs,xm,Mx,draw = 0;
1490: PetscBool view_final = PETSC_FALSE;
1491: PetscReal ptime;
1494: PetscInitialize(&argc,&argv,0,help);
1495: comm = PETSC_COMM_WORLD;
1496: PetscMemzero(&ctx,sizeof(ctx));
1498: /* Register limiters to be available on the command line */
1499: PetscFListAdd(&limiters,"upwind" ,"",(void(*)(void))Limit_Upwind);
1500: PetscFListAdd(&limiters,"lax-wendroff" ,"",(void(*)(void))Limit_LaxWendroff);
1501: PetscFListAdd(&limiters,"beam-warming" ,"",(void(*)(void))Limit_BeamWarming);
1502: PetscFListAdd(&limiters,"fromm" ,"",(void(*)(void))Limit_Fromm);
1503: PetscFListAdd(&limiters,"minmod" ,"",(void(*)(void))Limit_Minmod);
1504: PetscFListAdd(&limiters,"superbee" ,"",(void(*)(void))Limit_Superbee);
1505: PetscFListAdd(&limiters,"mc" ,"",(void(*)(void))Limit_MC);
1506: PetscFListAdd(&limiters,"vanleer" ,"",(void(*)(void))Limit_VanLeer);
1507: PetscFListAdd(&limiters,"vanalbada" ,"",(void(*)(void))Limit_VanAlbada);
1508: PetscFListAdd(&limiters,"vanalbadatvd" ,"",(void(*)(void))Limit_VanAlbadaTVD);
1509: PetscFListAdd(&limiters,"koren" ,"",(void(*)(void))Limit_Koren);
1510: PetscFListAdd(&limiters,"korensym" ,"",(void(*)(void))Limit_KorenSym);
1511: PetscFListAdd(&limiters,"koren3" ,"",(void(*)(void))Limit_Koren3);
1512: PetscFListAdd(&limiters,"cada-torrilhon2" ,"",(void(*)(void))Limit_CadaTorrilhon2);
1513: PetscFListAdd(&limiters,"cada-torrilhon3-r0p1","",(void(*)(void))Limit_CadaTorrilhon3R0p1);
1514: PetscFListAdd(&limiters,"cada-torrilhon3-r1" ,"",(void(*)(void))Limit_CadaTorrilhon3R1);
1515: PetscFListAdd(&limiters,"cada-torrilhon3-r10" ,"",(void(*)(void))Limit_CadaTorrilhon3R10);
1516: PetscFListAdd(&limiters,"cada-torrilhon3-r100","",(void(*)(void))Limit_CadaTorrilhon3R100);
1518: /* Register physical models to be available on the command line */
1519: PetscFListAdd(&physics,"advect" ,"",(void(*)(void))PhysicsCreate_Advect);
1520: PetscFListAdd(&physics,"burgers" ,"",(void(*)(void))PhysicsCreate_Burgers);
1521: PetscFListAdd(&physics,"traffic" ,"",(void(*)(void))PhysicsCreate_Traffic);
1522: PetscFListAdd(&physics,"acoustics" ,"",(void(*)(void))PhysicsCreate_Acoustics);
1523: PetscFListAdd(&physics,"isogas" ,"",(void(*)(void))PhysicsCreate_IsoGas);
1524: PetscFListAdd(&physics,"shallow" ,"",(void(*)(void))PhysicsCreate_Shallow);
1526: ctx.comm = comm;
1527: ctx.cfl = 0.9; ctx.bctype = FVBC_PERIODIC;
1528: ctx.xmin = -1; ctx.xmax = 1;
1529: PetscOptionsBegin(comm,PETSC_NULL,"Finite Volume solver options","");
1530: {
1531: PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,PETSC_NULL);
1532: PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,PETSC_NULL);
1533: PetscOptionsList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),PETSC_NULL);
1534: PetscOptionsList("-physics","Name of physics (Riemann solver and characteristics) to use","",physics,physname,physname,sizeof(physname),PETSC_NULL);
1535: PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,PETSC_NULL);
1536: PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof final_fname,&view_final);
1537: PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,PETSC_NULL);
1538: PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,PETSC_NULL);
1539: PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,PETSC_NULL);
1540: PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,PETSC_NULL);
1541: }
1542: PetscOptionsEnd();
1544: /* Choose the limiter from the list of registered limiters */
1545: PetscFListFind(limiters,comm,lname,PETSC_FALSE,(void(**)(void))&ctx.limit);
1546: if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);
1548: /* Choose the physics from the list of registered models */
1549: {
1550: PetscErrorCode (*r)(FVCtx*);
1551: PetscFListFind(physics,comm,physname,PETSC_FALSE,(void(**)(void))&r);
1552: if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
1553: /* Create the physics, will set the number of fields and their names */
1554: (*r)(&ctx);
1555: }
1557: /* Create a DMDA to manage the parallel grid */
1558: DMDACreate1d(comm,DMDA_BOUNDARY_PERIODIC,-50,ctx.physics.dof,2,PETSC_NULL,&da);
1559: /* Inform the DMDA of the field names provided by the physics. */
1560: /* The names will be shown in the title bars when run with -ts_monitor_solution */
1561: for (i=0; i<ctx.physics.dof; i++) {
1562: DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
1563: }
1564: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1565: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1567: /* Set coordinates of cell centers */
1568: DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);
1570: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1571: PetscMalloc4(dof*dof,PetscScalar,&ctx.R,dof*dof,PetscScalar,&ctx.Rinv,2*dof,PetscScalar,&ctx.cjmpLR,1*dof,PetscScalar,&ctx.cslope);
1572: PetscMalloc3(2*dof,PetscScalar,&ctx.uLR,dof,PetscScalar,&ctx.flux,dof,PetscReal,&ctx.speeds);
1574: /* Create a vector to store the solution and to save the initial state */
1575: DMCreateGlobalVector(da,&X);
1576: VecDuplicate(X,&X0);
1577: VecDuplicate(X,&R);
1579: DMCreateMatrix(da,PETSC_NULL,&B);
1581: /* Create a time-stepping object */
1582: TSCreate(comm,&ts);
1583: TSSetDM(ts,da);
1584: TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
1585: TSSetIJacobian(ts,B,B,FVIJacobian,&ctx);
1586: TSSetType(ts,TSSSP);
1587: TSSetDuration(ts,1000,10);
1589: /* Compute initial conditions and starting time step */
1590: FVSample(&ctx,da,0,X0);
1591: FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
1592: VecCopy(X0,X); /* The function value was not used so we set X=X0 again */
1593: TSSetInitialTimeStep(ts,0,ctx.cfl/ctx.cfl_idt);
1595: TSSetFromOptions(ts); /* Take runtime options */
1597: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1599: {
1600: PetscReal nrm1,nrmsup;
1601: PetscInt steps;
1603: TSSolve(ts,X,&ptime);
1604: TSGetTimeStepNumber(ts,&steps);
1606: PetscPrintf(comm,"Final time %8.5f, steps %d\n",ptime,steps);
1607: if (ctx.exact) {
1608: SolutionErrorNorms(&ctx,da,ptime,X,&nrm1,&nrmsup);
1609: PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e ||x-x_e||_sup %8.4e\n",nrm1,nrmsup);
1610: }
1611: }
1613: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1615: if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
1616: if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
1617: if (draw & 0x4) {
1618: Vec Y;
1619: VecDuplicate(X,&Y);
1620: FVSample(&ctx,da,ptime,Y);
1621: VecAYPX(Y,-1,X);
1622: VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1623: VecDestroy(&Y);
1624: }
1626: if (view_final) {
1627: PetscViewer viewer;
1628: PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1629: PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1630: VecView(X,viewer);
1631: PetscViewerDestroy(&viewer);
1632: }
1634: /* Clean up */
1635: (*ctx.physics.destroy)(ctx.physics.user);
1636: for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
1637: PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
1638: PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
1639: VecDestroy(&X);
1640: VecDestroy(&X0);
1641: VecDestroy(&R);
1642: MatDestroy(&B);
1643: DMDestroy(&da);
1644: TSDestroy(&ts);
1645: PetscFinalize();
1646: return 0;
1647: }