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: " acoustics - Acoustic wave propagation\n"
10: " u_t + (c*z*v)_x = 0\n"
11: " v_t + (c/z*u)_x = 0\n"
12: " isogas - Isothermal gas dynamics\n"
13: " rho_t + (rho*u)_x = 0\n"
14: " (rho*u)_t + (rho*u^2 + c^2*rho)_x = 0\n"
15: " shallow - Shallow water equations\n"
16: " h_t + (h*u)_x = 0\n"
17: " (h*u)_t + (h*u^2 + g*h^2/2)_x = 0\n"
18: "Some of these physical models have multiple Riemann solvers, select these with -physics_xxx_riemann\n"
19: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
20: " the states across shocks and rarefactions\n"
21: " roe - Linearized scheme, usually with an entropy fix inside sonic rarefactions\n"
22: "The systems provide a choice of reconstructions with -physics_xxx_reconstruct\n"
23: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
24: " conservative - Limit the conservative variables directly, can cause undesired interaction of waves\n\n"
25: "A variety of limiters for high-resolution TVD limiters are available with -limit\n"
26: " upwind,minmod,superbee,mc,vanleer,vanalbada,koren,cada-torillhon (last two are nominally third order)\n"
27: " and non-TVD schemes lax-wendroff,beam-warming,fromm\n\n"
28: "To preserve the TVD property, one should time step with a strong stability preserving method.\n"
29: "The optimal high order explicit Runge-Kutta methods in TSSSP are recommended for non-stiff problems.\n\n"
30: "Several initial conditions can be chosen with -initial N\n\n"
31: "The problem size should be set with -da_grid_x M\n\n";
33: #include <petscts.h>
34: #include <petscdm.h>
35: #include <petscdmda.h>
36: #include <petscdraw.h>
38: #include <petsc/private/kernels/blockinvert.h>
40: static inline PetscReal Sgn(PetscReal a) { return (a<0) ? -1 : 1; }
41: static inline PetscReal Abs(PetscReal a) { return (a<0) ? 0 : a; }
42: static inline PetscReal Sqr(PetscReal a) { return a*a; }
43: static inline PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; }
44: PETSC_UNUSED static inline PetscReal MinAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) < PetscAbs(b)) ? a : b; }
45: static inline PetscReal MinMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscAbs(b)); }
46: static inline PetscReal MaxMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMax(PetscAbs(a),PetscAbs(b)); }
47: static inline PetscReal MinMod3(PetscReal a,PetscReal b,PetscReal c) {return (a*b<0 || a*c<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscMin(PetscAbs(b),PetscAbs(c))); }
49: static inline PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }
51: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
52: typedef struct _LimitInfo {
53: PetscReal hx;
54: PetscInt m;
55: } *LimitInfo;
56: static void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
57: {
58: PetscInt i;
59: for (i=0; i<info->m; i++) lmt[i] = 0;
60: }
61: static void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
62: {
63: PetscInt i;
64: for (i=0; i<info->m; i++) lmt[i] = jR[i];
65: }
66: static void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
67: {
68: PetscInt i;
69: for (i=0; i<info->m; i++) lmt[i] = jL[i];
70: }
71: static void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
72: {
73: PetscInt i;
74: for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i] + jR[i]);
75: }
76: static void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
77: {
78: PetscInt i;
79: for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]);
80: }
81: static void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
82: {
83: PetscInt i;
84: for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]));
85: }
86: static void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
87: {
88: PetscInt i;
89: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]);
90: }
91: static void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
92: { /* phi = (t + abs(t)) / (1 + abs(t)) */
93: PetscInt i;
94: for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Abs(jR[i]) + Abs(jL[i])*jR[i]) / (Abs(jL[i]) + Abs(jR[i]) + 1e-15);
95: }
96: static void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
97: { /* phi = (t + t^2) / (1 + t^2) */
98: PetscInt i;
99: for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
100: }
101: static void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
102: { /* phi = (t + t^2) / (1 + t^2) */
103: PetscInt i;
104: for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0 : (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
105: }
106: static void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
107: { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
108: PetscInt i;
109: for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i]) + 2*Sqr(jL[i])*jR[i])/(2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
110: }
111: static void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
112: { /* Symmetric version of above */
113: PetscInt i;
114: for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i])/(2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
115: }
116: static void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
117: { /* Eq 11 of Cada-Torrilhon 2009 */
118: PetscInt i;
119: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]);
120: }
121: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)
122: {
123: return PetscMax(0,PetscMin((L+2*R)/3,PetscMax(-0.5*L,PetscMin(2*L,PetscMin((L+2*R)/3,1.6*R)))));
124: }
125: static void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
126: { /* Cada-Torrilhon 2009, Eq 13 */
127: PetscInt i;
128: for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]);
129: }
130: static void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
131: { /* Cada-Torrilhon 2009, Eq 22 */
132: /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
133: const PetscReal eps = 1e-7,hx = info->hx;
134: PetscInt i;
135: for (i=0; i<info->m; i++) {
136: const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r*hx);
137: lmt[i] = ((eta < 1-eps) ? (jL[i] + 2*jR[i]) / 3 : ((eta > 1+eps) ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]) : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3 + (1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]))));
138: }
139: }
140: static void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
141: {
142: Limit_CadaTorrilhon3R(0.1,info,jL,jR,lmt);
143: }
144: static void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
145: {
146: Limit_CadaTorrilhon3R(1,info,jL,jR,lmt);
147: }
148: static void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
149: {
150: Limit_CadaTorrilhon3R(10,info,jL,jR,lmt);
151: }
152: static void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
153: {
154: Limit_CadaTorrilhon3R(100,info,jL,jR,lmt);
155: }
157: /* --------------------------------- Finite Volume data structures ----------------------------------- */
159: typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType;
160: static const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0};
161: typedef PetscErrorCode (*RiemannFunction)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*);
162: typedef PetscErrorCode (*ReconstructFunction)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*,PetscReal*);
164: typedef struct {
165: PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*);
166: RiemannFunction riemann;
167: ReconstructFunction characteristic;
168: PetscErrorCode (*destroy)(void*);
169: void *user;
170: PetscInt dof;
171: char *fieldname[16];
172: } PhysicsCtx;
174: typedef struct {
175: void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*);
176: PhysicsCtx physics;
177: MPI_Comm comm;
178: char prefix[256];
180: /* Local work arrays */
181: PetscScalar *R,*Rinv; /* Characteristic basis, and it's inverse. COLUMN-MAJOR */
182: PetscScalar *cjmpLR; /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
183: PetscScalar *cslope; /* Limited slope, written in characteristic basis */
184: PetscScalar *uLR; /* Solution at left and right of interface, conservative variables, len=2*dof */
185: PetscScalar *flux; /* Flux across interface */
186: PetscReal *speeds; /* Speeds of each wave */
188: PetscReal cfl_idt; /* Max allowable value of 1/Delta t */
189: PetscReal cfl;
190: PetscReal xmin,xmax;
191: PetscInt initial;
192: PetscBool exact;
193: FVBCType bctype;
194: } FVCtx;
196: PetscErrorCode RiemannListAdd(PetscFunctionList *flist,const char *name,RiemannFunction rsolve)
197: {
199: PetscFunctionListAdd(flist,name,rsolve);
200: return 0;
201: }
203: PetscErrorCode RiemannListFind(PetscFunctionList flist,const char *name,RiemannFunction *rsolve)
204: {
206: PetscFunctionListFind(flist,name,rsolve);
208: return 0;
209: }
211: PetscErrorCode ReconstructListAdd(PetscFunctionList *flist,const char *name,ReconstructFunction r)
212: {
214: PetscFunctionListAdd(flist,name,r);
215: return 0;
216: }
218: PetscErrorCode ReconstructListFind(PetscFunctionList flist,const char *name,ReconstructFunction *r)
219: {
221: PetscFunctionListFind(flist,name,r);
223: return 0;
224: }
226: /* --------------------------------- Physics ----------------------------------- */
227: /*
228: Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction. These
229: are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the
230: number of fields and their names, and a function to deallocate private storage.
231: */
233: /* First a few functions useful to several different physics */
234: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
235: {
236: PetscInt i,j;
239: for (i=0; i<m; i++) {
240: for (j=0; j<m; j++) Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j);
241: speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
242: }
243: return 0;
244: }
246: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
247: {
249: PetscFree(vctx);
250: return 0;
251: }
253: /* --------------------------------- Advection ----------------------------------- */
255: typedef struct {
256: PetscReal a; /* advective velocity */
257: } AdvectCtx;
259: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
260: {
261: AdvectCtx *ctx = (AdvectCtx*)vctx;
262: PetscReal speed;
265: speed = ctx->a;
266: flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
267: *maxspeed = speed;
268: return 0;
269: }
271: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
272: {
273: AdvectCtx *ctx = (AdvectCtx*)vctx;
276: X[0] = 1.;
277: Xi[0] = 1.;
278: speeds[0] = ctx->a;
279: return 0;
280: }
282: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
283: {
284: AdvectCtx *ctx = (AdvectCtx*)vctx;
285: PetscReal a = ctx->a,x0;
288: switch (bctype) {
289: case FVBC_OUTFLOW: x0 = x-a*t; break;
290: case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
291: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
292: }
293: switch (initial) {
294: case 0: u[0] = (x0 < 0) ? 1 : -1; break;
295: case 1: u[0] = (x0 < 0) ? -1 : 1; break;
296: case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
297: case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
298: case 4: u[0] = PetscAbs(x0); break;
299: case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
300: case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
301: case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break;
302: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
303: }
304: return 0;
305: }
307: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
308: {
310: AdvectCtx *user;
313: PetscNew(&user);
314: ctx->physics.sample = PhysicsSample_Advect;
315: ctx->physics.riemann = PhysicsRiemann_Advect;
316: ctx->physics.characteristic = PhysicsCharacteristic_Advect;
317: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
318: ctx->physics.user = user;
319: ctx->physics.dof = 1;
320: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
321: user->a = 1;
322: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
323: {
324: PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
325: }
326: PetscOptionsEnd();
327: return 0;
328: }
330: /* --------------------------------- Burgers ----------------------------------- */
332: typedef struct {
333: PetscReal lxf_speed;
334: } BurgersCtx;
336: static PetscErrorCode PhysicsSample_Burgers(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
337: {
340: switch (initial) {
341: case 0: u[0] = (x < 0) ? 1 : -1; break;
342: case 1:
343: if (x < -t) u[0] = -1;
344: else if (x < t) u[0] = x/t;
345: else u[0] = 1;
346: break;
347: case 2:
348: if (x <= 0) u[0] = 0;
349: else if (x < t) u[0] = x/t;
350: else if (x < 1+0.5*t) u[0] = 1;
351: else u[0] = 0;
352: break;
353: case 3:
354: if (x < 0.2*t) u[0] = 0.2;
355: else if (x < t) u[0] = x/t;
356: else u[0] = 1;
357: break;
358: case 4:
360: u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin)));
361: break;
362: case 5: /* Pure shock solution */
363: if (x < 0.5*t) u[0] = 1;
364: else u[0] = 0;
365: break;
366: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
367: }
368: return 0;
369: }
371: static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
372: {
374: if (uL[0] < uR[0]) { /* rarefaction */
375: flux[0] = (uL[0]*uR[0] < 0)
376: ? 0 /* sonic rarefaction */
377: : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0]));
378: } else { /* shock */
379: flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0]));
380: }
381: *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
382: return 0;
383: }
385: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
386: {
387: PetscReal speed;
390: speed = 0.5*(uL[0] + uR[0]);
391: flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
392: if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
393: *maxspeed = speed;
394: return 0;
395: }
397: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
398: {
399: PetscReal c;
400: PetscScalar fL,fR;
403: c = ((BurgersCtx*)vctx)->lxf_speed;
404: fL = 0.5*PetscSqr(uL[0]);
405: fR = 0.5*PetscSqr(uR[0]);
406: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
407: *maxspeed = c;
408: return 0;
409: }
411: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
412: {
413: PetscReal c;
414: PetscScalar fL,fR;
417: c = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0]));
418: fL = 0.5*PetscSqr(uL[0]);
419: fR = 0.5*PetscSqr(uR[0]);
420: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
421: *maxspeed = c;
422: return 0;
423: }
425: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
426: {
427: BurgersCtx *user;
428: PetscErrorCode ierr;
429: RiemannFunction r;
430: PetscFunctionList rlist = 0;
431: char rname[256] = "exact";
434: PetscNew(&user);
436: ctx->physics.sample = PhysicsSample_Burgers;
437: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
438: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
439: ctx->physics.user = user;
440: ctx->physics.dof = 1;
442: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
443: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Burgers_Exact);
444: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Burgers_Roe);
445: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Burgers_LxF);
446: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Burgers_Rusanov);
447: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
448: {
449: PetscOptionsFList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
450: }
451: PetscOptionsEnd();
452: RiemannListFind(rlist,rname,&r);
453: PetscFunctionListDestroy(&rlist);
454: ctx->physics.riemann = r;
456: /* *
457: * Hack to deal with LxF in semi-discrete form
458: * max speed is 1 for the basic initial conditions (where |u| <= 1)
459: * */
460: if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
461: return 0;
462: }
464: /* --------------------------------- Traffic ----------------------------------- */
466: typedef struct {
467: PetscReal lxf_speed;
468: PetscReal a;
469: } TrafficCtx;
471: static inline PetscScalar TrafficFlux(PetscScalar a,PetscScalar u) { return a*u*(1-u); }
473: static PetscErrorCode PhysicsSample_Traffic(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
474: {
475: PetscReal a = ((TrafficCtx*)vctx)->a;
479: switch (initial) {
480: case 0:
481: u[0] = (-a*t < x) ? 2 : 0; break;
482: case 1:
483: if (x < PetscMin(2*a*t,0.5+a*t)) u[0] = -1;
484: else if (x < 1) u[0] = 0;
485: else u[0] = 1;
486: break;
487: case 2:
489: u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin)));
490: break;
491: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
492: }
493: return 0;
494: }
496: static PetscErrorCode PhysicsRiemann_Traffic_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
497: {
498: PetscReal a = ((TrafficCtx*)vctx)->a;
501: if (uL[0] < uR[0]) {
502: flux[0] = PetscMin(TrafficFlux(a,uL[0]),TrafficFlux(a,uR[0]));
503: } else {
504: flux[0] = (uR[0] < 0.5 && 0.5 < uL[0]) ? TrafficFlux(a,0.5) : PetscMax(TrafficFlux(a,uL[0]),TrafficFlux(a,uR[0]));
505: }
506: *maxspeed = a*MaxAbs(1-2*uL[0],1-2*uR[0]);
507: return 0;
508: }
510: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
511: {
512: PetscReal a = ((TrafficCtx*)vctx)->a;
513: PetscReal speed;
516: speed = a*(1 - (uL[0] + uR[0]));
517: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
518: *maxspeed = speed;
519: return 0;
520: }
522: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
523: {
524: TrafficCtx *phys = (TrafficCtx*)vctx;
525: PetscReal a = phys->a;
526: PetscReal speed;
529: speed = a*(1 - (uL[0] + uR[0]));
530: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]);
531: *maxspeed = speed;
532: return 0;
533: }
535: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
536: {
537: PetscReal a = ((TrafficCtx*)vctx)->a;
538: PetscReal speed;
541: speed = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0]));
542: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]);
543: *maxspeed = speed;
544: return 0;
545: }
547: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
548: {
549: PetscErrorCode ierr;
550: TrafficCtx *user;
551: RiemannFunction r;
552: PetscFunctionList rlist = 0;
553: char rname[256] = "exact";
556: PetscNew(&user);
557: ctx->physics.sample = PhysicsSample_Traffic;
558: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
559: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
560: ctx->physics.user = user;
561: ctx->physics.dof = 1;
563: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
564: user->a = 0.5;
565: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Traffic_Exact);
566: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Traffic_Roe);
567: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Traffic_LxF);
568: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Traffic_Rusanov);
569: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");
570: PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,NULL);
571: PetscOptionsFList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
572: PetscOptionsEnd();
574: RiemannListFind(rlist,rname,&r);
575: PetscFunctionListDestroy(&rlist);
577: ctx->physics.riemann = r;
579: /* *
580: * Hack to deal with LxF in semi-discrete form
581: * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
582: * */
583: if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3*user->a;
584: return 0;
585: }
587: /* --------------------------------- Linear Acoustics ----------------------------------- */
589: /* Flux: u_t + (A u)_x
590: * z = sqrt(rho*bulk), c = sqrt(rho/bulk)
591: * Spectral decomposition: A = R * D * Rinv
592: * [ cz] = [-z z] [-c ] [-1/2z 1/2]
593: * [c/z ] = [ 1 1] [ c] [ 1/2z 1/2]
594: *
595: * We decompose this into the left-traveling waves Al = R * D^- Rinv
596: * and the right-traveling waves Ar = R * D^+ * Rinv
597: * Multiplying out these expressions produces the following two matrices
598: */
600: typedef struct {
601: PetscReal c; /* speed of sound: c = sqrt(bulk/rho) */
602: PetscReal z; /* impedence: z = sqrt(rho*bulk) */
603: } AcousticsCtx;
605: PETSC_UNUSED static inline void AcousticsFlux(AcousticsCtx *ctx,const PetscScalar *u,PetscScalar *f)
606: {
607: f[0] = ctx->c*ctx->z*u[1];
608: f[1] = ctx->c/ctx->z*u[0];
609: }
611: static PetscErrorCode PhysicsCharacteristic_Acoustics(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
612: {
613: AcousticsCtx *phys = (AcousticsCtx*)vctx;
614: PetscReal z = phys->z,c = phys->c;
617: X[0*2+0] = -z;
618: X[0*2+1] = z;
619: X[1*2+0] = 1;
620: X[1*2+1] = 1;
621: Xi[0*2+0] = -1./(2*z);
622: Xi[0*2+1] = 1./2;
623: Xi[1*2+0] = 1./(2*z);
624: Xi[1*2+1] = 1./2;
625: speeds[0] = -c;
626: speeds[1] = c;
627: return 0;
628: }
630: static PetscErrorCode PhysicsSample_Acoustics_Initial(AcousticsCtx *phys,PetscInt initial,PetscReal xmin,PetscReal xmax,PetscReal x,PetscReal *u)
631: {
633: switch (initial) {
634: case 0:
635: u[0] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.2) < 0.1) ? 1 : 0.5;
636: u[1] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.7) < 0.1) ? 1 : -0.5;
637: break;
638: case 1:
639: u[0] = PetscCosReal(3 * 2*PETSC_PI*x/(xmax-xmin));
640: u[1] = PetscExpReal(-PetscSqr(x - (xmax + xmin)/2) / (2*PetscSqr(0.2*(xmax - xmin)))) - 0.5;
641: break;
642: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
643: }
644: return 0;
645: }
647: static PetscErrorCode PhysicsSample_Acoustics(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
648: {
649: AcousticsCtx *phys = (AcousticsCtx*)vctx;
650: PetscReal c = phys->c;
651: PetscReal x0a,x0b,u0a[2],u0b[2],tmp[2];
652: PetscReal X[2][2],Xi[2][2],dummy[2];
655: switch (bctype) {
656: case FVBC_OUTFLOW:
657: x0a = x+c*t;
658: x0b = x-c*t;
659: break;
660: case FVBC_PERIODIC:
661: x0a = RangeMod(x+c*t,xmin,xmax);
662: x0b = RangeMod(x-c*t,xmin,xmax);
663: break;
664: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType");
665: }
666: PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0a,u0a);
667: PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0b,u0b);
668: PhysicsCharacteristic_Acoustics(vctx,2,u,&X[0][0],&Xi[0][0],dummy);
669: tmp[0] = Xi[0][0]*u0a[0] + Xi[0][1]*u0a[1];
670: tmp[1] = Xi[1][0]*u0b[0] + Xi[1][1]*u0b[1];
671: u[0] = X[0][0]*tmp[0] + X[0][1]*tmp[1];
672: u[1] = X[1][0]*tmp[0] + X[1][1]*tmp[1];
673: return 0;
674: }
676: static PetscErrorCode PhysicsRiemann_Acoustics_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
677: {
678: AcousticsCtx *phys = (AcousticsCtx*)vctx;
679: PetscReal c = phys->c,z = phys->z;
680: PetscReal
681: Al[2][2] = {{-c/2 , c*z/2 },
682: {c/(2*z) , -c/2 }}, /* Left traveling waves */
683: Ar[2][2] = {{c/2 , c*z/2 },
684: {c/(2*z) , c/2 }}; /* Right traveling waves */
687: flux[0] = Al[0][0]*uR[0] + Al[0][1]*uR[1] + Ar[0][0]*uL[0] + Ar[0][1]*uL[1];
688: flux[1] = Al[1][0]*uR[0] + Al[1][1]*uR[1] + Ar[1][0]*uL[0] + Ar[1][1]*uL[1];
689: *maxspeed = c;
690: return 0;
691: }
693: static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx)
694: {
695: PetscErrorCode ierr;
696: AcousticsCtx *user;
697: PetscFunctionList rlist = 0,rclist = 0;
698: char rname[256] = "exact",rcname[256] = "characteristic";
701: PetscNew(&user);
702: ctx->physics.sample = PhysicsSample_Acoustics;
703: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
704: ctx->physics.user = user;
705: ctx->physics.dof = 2;
707: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
708: PetscStrallocpy("v",&ctx->physics.fieldname[1]);
710: user->c = 1;
711: user->z = 1;
713: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Acoustics_Exact);
714: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Acoustics);
715: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
716: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for linear Acoustics","");
717: {
718: PetscOptionsReal("-physics_acoustics_c","c = sqrt(bulk/rho)","",user->c,&user->c,NULL);
719: PetscOptionsReal("-physics_acoustics_z","z = sqrt(bulk*rho)","",user->z,&user->z,NULL);
720: PetscOptionsFList("-physics_acoustics_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
721: PetscOptionsFList("-physics_acoustics_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
722: }
723: PetscOptionsEnd();
724: RiemannListFind(rlist,rname,&ctx->physics.riemann);
725: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
726: PetscFunctionListDestroy(&rlist);
727: PetscFunctionListDestroy(&rclist);
728: return 0;
729: }
731: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */
733: typedef struct {
734: PetscReal acoustic_speed;
735: } IsoGasCtx;
737: static inline void IsoGasFlux(PetscReal c,const PetscScalar *u,PetscScalar *f)
738: {
739: f[0] = u[1];
740: f[1] = PetscSqr(u[1])/u[0] + c*c*u[0];
741: }
743: static PetscErrorCode PhysicsSample_IsoGas(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
744: {
747: switch (initial) {
748: case 0:
749: u[0] = (x < 0) ? 1 : 0.5;
750: u[1] = (x < 0) ? 1 : 0.7;
751: break;
752: case 1:
753: u[0] = 1+0.5*PetscSinReal(2*PETSC_PI*x);
754: u[1] = 1*u[0];
755: break;
756: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition");
757: }
758: return 0;
759: }
761: static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
762: {
763: IsoGasCtx *phys = (IsoGasCtx*)vctx;
764: PetscReal c = phys->acoustic_speed;
765: PetscScalar ubar,du[2],a[2],fL[2],fR[2],lam[2],ustar[2],R[2][2];
766: PetscInt i;
769: ubar = (uL[1]/PetscSqrtScalar(uL[0]) + uR[1]/PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0]));
770: /* write fluxuations in characteristic basis */
771: du[0] = uR[0] - uL[0];
772: du[1] = uR[1] - uL[1];
773: a[0] = (1/(2*c)) * ((ubar + c)*du[0] - du[1]);
774: a[1] = (1/(2*c)) * ((-ubar + c)*du[0] + du[1]);
775: /* wave speeds */
776: lam[0] = ubar - c;
777: lam[1] = ubar + c;
778: /* Right eigenvectors */
779: R[0][0] = 1; R[0][1] = ubar-c;
780: R[1][0] = 1; R[1][1] = ubar+c;
781: /* Compute state in star region (between the 1-wave and 2-wave) */
782: for (i=0; i<2; i++) ustar[i] = uL[i] + a[0]*R[0][i];
783: if (uL[1]/uL[0] < c && c < ustar[1]/ustar[0]) { /* 1-wave is sonic rarefaction */
784: PetscScalar ufan[2];
785: ufan[0] = uL[0]*PetscExpScalar(uL[1]/(uL[0]*c) - 1);
786: ufan[1] = c*ufan[0];
787: IsoGasFlux(c,ufan,flux);
788: } else if (ustar[1]/ustar[0] < -c && -c < uR[1]/uR[0]) { /* 2-wave is sonic rarefaction */
789: PetscScalar ufan[2];
790: ufan[0] = uR[0]*PetscExpScalar(-uR[1]/(uR[0]*c) - 1);
791: ufan[1] = -c*ufan[0];
792: IsoGasFlux(c,ufan,flux);
793: } else { /* Centered form */
794: IsoGasFlux(c,uL,fL);
795: IsoGasFlux(c,uR,fR);
796: for (i=0; i<2; i++) {
797: PetscScalar absdu = PetscAbsScalar(lam[0])*a[0]*R[0][i] + PetscAbsScalar(lam[1])*a[1]*R[1][i];
798: flux[i] = 0.5*(fL[i]+fR[i]) - 0.5*absdu;
799: }
800: }
801: *maxspeed = MaxAbs(lam[0],lam[1]);
802: return 0;
803: }
805: static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
806: {
807: IsoGasCtx *phys = (IsoGasCtx*)vctx;
808: PetscReal c = phys->acoustic_speed;
809: PetscScalar ustar[2];
810: struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
811: PetscInt i;
815: {
816: /* Solve for star state */
817: PetscScalar res,tmp,rho = 0.5*(L.rho + R.rho); /* initial guess */
818: for (i=0; i<20; i++) {
819: PetscScalar fr,fl,dfr,dfl;
820: fl = (L.rho < rho)
821: ? (rho-L.rho)/PetscSqrtScalar(L.rho*rho) /* shock */
822: : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
823: fr = (R.rho < rho)
824: ? (rho-R.rho)/PetscSqrtScalar(R.rho*rho) /* shock */
825: : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
826: res = R.u-L.u + c*(fr+fl);
828: if (PetscAbsScalar(res) < 1e-10) {
829: star.rho = rho;
830: star.u = L.u - c*fl;
831: goto converged;
832: }
833: dfl = (L.rho < rho) ? 1/PetscSqrtScalar(L.rho*rho)*(1 - 0.5*(rho-L.rho)/rho) : 1/rho;
834: dfr = (R.rho < rho) ? 1/PetscSqrtScalar(R.rho*rho)*(1 - 0.5*(rho-R.rho)/rho) : 1/rho;
835: tmp = rho - res/(c*(dfr+dfl));
836: if (tmp <= 0) rho /= 2; /* Guard against Newton shooting off to a negative density */
837: else rho = tmp;
839: }
840: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Newton iteration for star.rho diverged after %D iterations",i);
841: }
842: converged:
843: if (L.u-c < 0 && 0 < star.u-c) { /* 1-wave is sonic rarefaction */
844: PetscScalar ufan[2];
845: ufan[0] = L.rho*PetscExpScalar(L.u/c - 1);
846: ufan[1] = c*ufan[0];
847: IsoGasFlux(c,ufan,flux);
848: } else if (star.u+c < 0 && 0 < R.u+c) { /* 2-wave is sonic rarefaction */
849: PetscScalar ufan[2];
850: ufan[0] = R.rho*PetscExpScalar(-R.u/c - 1);
851: ufan[1] = -c*ufan[0];
852: IsoGasFlux(c,ufan,flux);
853: } else if ((L.rho >= star.rho && L.u-c >= 0) || (L.rho < star.rho && (star.rho*star.u-L.rho*L.u)/(star.rho-L.rho) > 0)) {
854: /* 1-wave is supersonic rarefaction, or supersonic shock */
855: IsoGasFlux(c,uL,flux);
856: } else if ((star.rho <= R.rho && R.u+c <= 0) || (star.rho > R.rho && (R.rho*R.u-star.rho*star.u)/(R.rho-star.rho) < 0)) {
857: /* 2-wave is supersonic rarefaction or supersonic shock */
858: IsoGasFlux(c,uR,flux);
859: } else {
860: ustar[0] = star.rho;
861: ustar[1] = star.rho*star.u;
862: IsoGasFlux(c,ustar,flux);
863: }
864: *maxspeed = MaxAbs(MaxAbs(star.u-c,star.u+c),MaxAbs(L.u-c,R.u+c));
865: return 0;
866: }
868: static PetscErrorCode PhysicsRiemann_IsoGas_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
869: {
870: IsoGasCtx *phys = (IsoGasCtx*)vctx;
871: PetscScalar c = phys->acoustic_speed,fL[2],fR[2],s;
872: struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
876: IsoGasFlux(c,uL,fL);
877: IsoGasFlux(c,uR,fR);
878: s = PetscMax(PetscAbs(L.u),PetscAbs(R.u))+c;
879: flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
880: flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
881: *maxspeed = s;
882: return 0;
883: }
885: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
886: {
887: IsoGasCtx *phys = (IsoGasCtx*)vctx;
888: PetscReal c = phys->acoustic_speed;
891: speeds[0] = u[1]/u[0] - c;
892: speeds[1] = u[1]/u[0] + c;
893: X[0*2+0] = 1;
894: X[0*2+1] = speeds[0];
895: X[1*2+0] = 1;
896: X[1*2+1] = speeds[1];
897: PetscArraycpy(Xi,X,4);
898: PetscKernel_A_gets_inverse_A_2(Xi,0,PETSC_FALSE,NULL);
899: return 0;
900: }
902: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
903: {
904: PetscErrorCode ierr;
905: IsoGasCtx *user;
906: PetscFunctionList rlist = 0,rclist = 0;
907: char rname[256] = "exact",rcname[256] = "characteristic";
910: PetscNew(&user);
911: ctx->physics.sample = PhysicsSample_IsoGas;
912: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
913: ctx->physics.user = user;
914: ctx->physics.dof = 2;
916: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
917: PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
919: user->acoustic_speed = 1;
921: RiemannListAdd(&rlist,"exact", PhysicsRiemann_IsoGas_Exact);
922: RiemannListAdd(&rlist,"roe", PhysicsRiemann_IsoGas_Roe);
923: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_IsoGas_Rusanov);
924: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_IsoGas);
925: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
926: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");
927: PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,NULL);
928: PetscOptionsFList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
929: PetscOptionsFList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
930: PetscOptionsEnd();
931: RiemannListFind(rlist,rname,&ctx->physics.riemann);
932: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
933: PetscFunctionListDestroy(&rlist);
934: PetscFunctionListDestroy(&rclist);
935: return 0;
936: }
938: /* --------------------------------- Shallow Water ----------------------------------- */
939: typedef struct {
940: PetscReal gravity;
941: } ShallowCtx;
943: static inline void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
944: {
945: f[0] = u[1];
946: f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
947: }
949: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
950: {
951: ShallowCtx *phys = (ShallowCtx*)vctx;
952: PetscScalar g = phys->gravity,ustar[2],cL,cR,c,cstar;
953: struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
954: PetscInt i;
958: cL = PetscSqrtScalar(g*L.h);
959: cR = PetscSqrtScalar(g*R.h);
960: c = PetscMax(cL,cR);
961: {
962: /* Solve for star state */
963: const PetscInt maxits = 50;
964: PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */
965: h0 = h;
966: for (i=0; i<maxits; i++) {
967: PetscScalar fr,fl,dfr,dfl;
968: fl = (L.h < h)
969: ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */
970: : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h); /* rarefaction */
971: fr = (R.h < h)
972: ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */
973: : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h); /* rarefaction */
974: res = R.u - L.u + fr + fl;
976: if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h-h0) < 1e-8)) {
977: star.h = h;
978: star.u = L.u - fl;
979: goto converged;
980: } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */
981: h = 0.8*h0 + 0.2*h;
982: continue;
983: }
984: /* Accept the last step and take another */
985: res0 = res;
986: h0 = h;
987: dfl = (L.h < h) ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h) : PetscSqrtScalar(g/h);
988: dfr = (R.h < h) ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h) : PetscSqrtScalar(g/h);
989: tmp = h - res/(dfr+dfl);
990: if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */
991: else h = tmp;
993: }
994: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Newton iteration for star.h diverged after %D iterations",i);
995: }
996: converged:
997: cstar = PetscSqrtScalar(g*star.h);
998: if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */
999: PetscScalar ufan[2];
1000: ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL);
1001: ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0];
1002: ShallowFlux(phys,ufan,flux);
1003: } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */
1004: PetscScalar ufan[2];
1005: ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR);
1006: ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0];
1007: ShallowFlux(phys,ufan,flux);
1008: } else if ((L.h >= star.h && L.u-c >= 0) || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) {
1009: /* 1-wave is right-travelling shock (supersonic) */
1010: ShallowFlux(phys,uL,flux);
1011: } else if ((star.h <= R.h && R.u+c <= 0) || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) {
1012: /* 2-wave is left-travelling shock (supersonic) */
1013: ShallowFlux(phys,uR,flux);
1014: } else {
1015: ustar[0] = star.h;
1016: ustar[1] = star.h*star.u;
1017: ShallowFlux(phys,ustar,flux);
1018: }
1019: *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR));
1020: return 0;
1021: }
1023: static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
1024: {
1025: ShallowCtx *phys = (ShallowCtx*)vctx;
1026: PetscScalar g = phys->gravity,fL[2],fR[2],s;
1027: struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
1031: ShallowFlux(phys,uL,fL);
1032: ShallowFlux(phys,uR,fR);
1033: s = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h));
1034: flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
1035: flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
1036: *maxspeed = s;
1037: return 0;
1038: }
1040: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
1041: {
1042: ShallowCtx *phys = (ShallowCtx*)vctx;
1043: PetscReal c;
1046: c = PetscSqrtScalar(u[0]*phys->gravity);
1047: speeds[0] = u[1]/u[0] - c;
1048: speeds[1] = u[1]/u[0] + c;
1049: X[0*2+0] = 1;
1050: X[0*2+1] = speeds[0];
1051: X[1*2+0] = 1;
1052: X[1*2+1] = speeds[1];
1053: PetscArraycpy(Xi,X,4);
1054: PetscKernel_A_gets_inverse_A_2(Xi,0,PETSC_FALSE,NULL);
1055: return 0;
1056: }
1058: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
1059: {
1060: PetscErrorCode ierr;
1061: ShallowCtx *user;
1062: PetscFunctionList rlist = 0,rclist = 0;
1063: char rname[256] = "exact",rcname[256] = "characteristic";
1066: PetscNew(&user);
1067: /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
1068: ctx->physics.sample = PhysicsSample_IsoGas;
1069: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
1070: ctx->physics.user = user;
1071: ctx->physics.dof = 2;
1073: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1074: PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1076: user->gravity = 1;
1078: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Shallow_Exact);
1079: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov);
1080: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Shallow);
1081: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1082: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");
1083: PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,NULL);
1084: PetscOptionsFList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
1085: PetscOptionsFList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
1086: PetscOptionsEnd();
1087: RiemannListFind(rlist,rname,&ctx->physics.riemann);
1088: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1089: PetscFunctionListDestroy(&rlist);
1090: PetscFunctionListDestroy(&rclist);
1091: return 0;
1092: }
1094: /* --------------------------------- Finite Volume Solver ----------------------------------- */
1096: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
1097: {
1098: FVCtx *ctx = (FVCtx*)vctx;
1099: PetscInt i,j,k,Mx,dof,xs,xm;
1100: PetscReal hx,cfl_idt = 0;
1101: PetscScalar *x,*f,*slope;
1102: Vec Xloc;
1103: DM da;
1106: TSGetDM(ts,&da);
1107: DMGetLocalVector(da,&Xloc);
1108: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1109: hx = (ctx->xmax - ctx->xmin)/Mx;
1110: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1111: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
1113: VecZeroEntries(F);
1115: DMDAVecGetArray(da,Xloc,&x);
1116: DMDAVecGetArray(da,F,&f);
1117: DMDAGetArray(da,PETSC_TRUE,&slope);
1119: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1121: if (ctx->bctype == FVBC_OUTFLOW) {
1122: for (i=xs-2; i<0; i++) {
1123: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
1124: }
1125: for (i=Mx; i<xs+xm+2; i++) {
1126: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
1127: }
1128: }
1129: for (i=xs-1; i<xs+xm+1; i++) {
1130: struct _LimitInfo info;
1131: PetscScalar *cjmpL,*cjmpR;
1132: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
1133: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1134: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
1135: PetscArrayzero(ctx->cjmpLR,2*dof);
1136: cjmpL = &ctx->cjmpLR[0];
1137: cjmpR = &ctx->cjmpLR[dof];
1138: for (j=0; j<dof; j++) {
1139: PetscScalar jmpL,jmpR;
1140: jmpL = x[(i+0)*dof+j] - x[(i-1)*dof+j];
1141: jmpR = x[(i+1)*dof+j] - x[(i+0)*dof+j];
1142: for (k=0; k<dof; k++) {
1143: cjmpL[k] += ctx->Rinv[k+j*dof] * jmpL;
1144: cjmpR[k] += ctx->Rinv[k+j*dof] * jmpR;
1145: }
1146: }
1147: /* Apply limiter to the left and right characteristic jumps */
1148: info.m = dof;
1149: info.hx = hx;
1150: (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
1151: for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
1152: for (j=0; j<dof; j++) {
1153: PetscScalar tmp = 0;
1154: for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof] * ctx->cslope[k];
1155: slope[i*dof+j] = tmp;
1156: }
1157: }
1159: for (i=xs; i<xs+xm+1; i++) {
1160: PetscReal maxspeed;
1161: PetscScalar *uL,*uR;
1162: uL = &ctx->uLR[0];
1163: uR = &ctx->uLR[dof];
1164: for (j=0; j<dof; j++) {
1165: uL[j] = x[(i-1)*dof+j] + slope[(i-1)*dof+j]*hx/2;
1166: uR[j] = x[(i-0)*dof+j] - slope[(i-0)*dof+j]*hx/2;
1167: }
1168: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed);
1169: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
1171: if (i > xs) {
1172: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
1173: }
1174: if (i < xs+xm) {
1175: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
1176: }
1177: }
1179: DMDAVecRestoreArray(da,Xloc,&x);
1180: DMDAVecRestoreArray(da,F,&f);
1181: DMDARestoreArray(da,PETSC_TRUE,&slope);
1182: DMRestoreLocalVector(da,&Xloc);
1184: MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));
1185: if (0) {
1186: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
1187: PetscReal dt,tnow;
1188: TSGetTimeStep(ts,&dt);
1189: TSGetTime(ts,&tnow);
1190: if (dt > 0.5/ctx->cfl_idt) {
1191: PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
1192: }
1193: }
1194: return 0;
1195: }
1197: static PetscErrorCode SmallMatMultADB(PetscScalar *C,PetscInt bs,const PetscScalar *A,const PetscReal *D,const PetscScalar *B)
1198: {
1199: PetscInt i,j,k;
1202: for (i=0; i<bs; i++) {
1203: for (j=0; j<bs; j++) {
1204: PetscScalar tmp = 0;
1205: for (k=0; k<bs; k++) tmp += A[i*bs+k] * D[k] * B[k*bs+j];
1206: C[i*bs+j] = tmp;
1207: }
1208: }
1209: return 0;
1210: }
1212: static PetscErrorCode FVIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void *vctx)
1213: {
1214: FVCtx *ctx = (FVCtx*)vctx;
1215: PetscInt i,j,dof = ctx->physics.dof;
1216: PetscScalar *J;
1217: const PetscScalar *x;
1218: PetscReal hx;
1219: DM da;
1220: DMDALocalInfo dainfo;
1223: TSGetDM(ts,&da);
1224: DMDAVecGetArrayRead(da,X,(void*)&x);
1225: DMDAGetLocalInfo(da,&dainfo);
1226: hx = (ctx->xmax - ctx->xmin)/dainfo.mx;
1227: PetscMalloc1(dof*dof,&J);
1228: for (i=dainfo.xs; i<dainfo.xs+dainfo.xm; i++) {
1229: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1230: for (j=0; j<dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]);
1231: SmallMatMultADB(J,dof,ctx->R,ctx->speeds,ctx->Rinv);
1232: for (j=0; j<dof*dof; j++) J[j] = J[j]/hx + shift*(j/dof == j%dof);
1233: MatSetValuesBlocked(B,1,&i,1,&i,J,INSERT_VALUES);
1234: }
1235: PetscFree(J);
1236: DMDAVecRestoreArrayRead(da,X,(void*)&x);
1238: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1239: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1240: if (A != B) {
1241: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1242: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1243: }
1244: return 0;
1245: }
1247: static PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
1248: {
1249: PetscScalar *u,*uj;
1250: PetscInt i,j,k,dof,xs,xm,Mx;
1254: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1255: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1256: DMDAVecGetArray(da,U,&u);
1257: PetscMalloc1(dof,&uj);
1258: for (i=xs; i<xs+xm; i++) {
1259: const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h;
1260: const PetscInt N = 200;
1261: /* Integrate over cell i using trapezoid rule with N points. */
1262: for (k=0; k<dof; k++) u[i*dof+k] = 0;
1263: for (j=0; j<N+1; j++) {
1264: PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N;
1265: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
1266: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
1267: }
1268: }
1269: DMDAVecRestoreArray(da,U,&u);
1270: PetscFree(uj);
1271: return 0;
1272: }
1274: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
1275: {
1276: PetscReal xmin,xmax;
1277: PetscScalar sum,tvsum,tvgsum;
1278: const PetscScalar *x;
1279: PetscInt imin,imax,Mx,i,j,xs,xm,dof;
1280: Vec Xloc;
1281: PetscBool iascii;
1284: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1285: if (iascii) {
1286: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
1287: DMGetLocalVector(da,&Xloc);
1288: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1289: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
1290: DMDAVecGetArrayRead(da,Xloc,(void*)&x);
1291: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1292: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1293: tvsum = 0;
1294: for (i=xs; i<xs+xm; i++) {
1295: for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j] - x[(i-1)*dof+j]);
1296: }
1297: MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)da));
1298: DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);
1299: DMRestoreLocalVector(da,&Xloc);
1301: VecMin(X,&imin,&xmin);
1302: VecMax(X,&imax,&xmax);
1303: VecSum(X,&sum);
1304: PetscViewerASCIIPrintf(viewer,"Solution range [%8.5f,%8.5f] with extrema at %D and %D, mean %8.5f, ||x||_TV %8.5f\n",(double)xmin,(double)xmax,imin,imax,(double)(sum/Mx),(double)(tvgsum/Mx));
1305: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported");
1306: return 0;
1307: }
1309: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1,PetscReal *nrmsup)
1310: {
1311: Vec Y;
1312: PetscInt Mx;
1315: VecGetSize(X,&Mx);
1316: VecDuplicate(X,&Y);
1317: FVSample(ctx,da,t,Y);
1318: VecAYPX(Y,-1,X);
1319: VecNorm(Y,NORM_1,nrm1);
1320: VecNorm(Y,NORM_INFINITY,nrmsup);
1321: *nrm1 /= Mx;
1322: VecDestroy(&Y);
1323: return 0;
1324: }
1326: int main(int argc,char *argv[])
1327: {
1328: char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
1329: PetscFunctionList limiters = 0,physics = 0;
1330: MPI_Comm comm;
1331: TS ts;
1332: DM da;
1333: Vec X,X0,R;
1334: Mat B;
1335: FVCtx ctx;
1336: PetscInt i,dof,xs,xm,Mx,draw = 0;
1337: PetscBool view_final = PETSC_FALSE;
1338: PetscReal ptime;
1339: PetscErrorCode ierr;
1341: PetscInitialize(&argc,&argv,0,help);
1342: comm = PETSC_COMM_WORLD;
1343: PetscMemzero(&ctx,sizeof(ctx));
1345: /* Register limiters to be available on the command line */
1346: PetscFunctionListAdd(&limiters,"upwind" ,Limit_Upwind);
1347: PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit_LaxWendroff);
1348: PetscFunctionListAdd(&limiters,"beam-warming" ,Limit_BeamWarming);
1349: PetscFunctionListAdd(&limiters,"fromm" ,Limit_Fromm);
1350: PetscFunctionListAdd(&limiters,"minmod" ,Limit_Minmod);
1351: PetscFunctionListAdd(&limiters,"superbee" ,Limit_Superbee);
1352: PetscFunctionListAdd(&limiters,"mc" ,Limit_MC);
1353: PetscFunctionListAdd(&limiters,"vanleer" ,Limit_VanLeer);
1354: PetscFunctionListAdd(&limiters,"vanalbada" ,Limit_VanAlbada);
1355: PetscFunctionListAdd(&limiters,"vanalbadatvd" ,Limit_VanAlbadaTVD);
1356: PetscFunctionListAdd(&limiters,"koren" ,Limit_Koren);
1357: PetscFunctionListAdd(&limiters,"korensym" ,Limit_KorenSym);
1358: PetscFunctionListAdd(&limiters,"koren3" ,Limit_Koren3);
1359: PetscFunctionListAdd(&limiters,"cada-torrilhon2" ,Limit_CadaTorrilhon2);
1360: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1);
1361: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1" ,Limit_CadaTorrilhon3R1);
1362: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10);
1363: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100);
1365: /* Register physical models to be available on the command line */
1366: PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect);
1367: PetscFunctionListAdd(&physics,"burgers" ,PhysicsCreate_Burgers);
1368: PetscFunctionListAdd(&physics,"traffic" ,PhysicsCreate_Traffic);
1369: PetscFunctionListAdd(&physics,"acoustics" ,PhysicsCreate_Acoustics);
1370: PetscFunctionListAdd(&physics,"isogas" ,PhysicsCreate_IsoGas);
1371: PetscFunctionListAdd(&physics,"shallow" ,PhysicsCreate_Shallow);
1373: ctx.comm = comm;
1374: ctx.cfl = 0.9; ctx.bctype = FVBC_PERIODIC;
1375: ctx.xmin = -1; ctx.xmax = 1;
1376: PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
1377: PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
1378: PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
1379: PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);
1380: PetscOptionsFList("-physics","Name of physics (Riemann solver and characteristics) to use","",physics,physname,physname,sizeof(physname),NULL);
1381: PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
1382: PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
1383: PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
1384: PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
1385: PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
1386: PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
1387: PetscOptionsEnd();
1389: /* Choose the limiter from the list of registered limiters */
1390: PetscFunctionListFind(limiters,lname,&ctx.limit);
1393: /* Choose the physics from the list of registered models */
1394: {
1395: PetscErrorCode (*r)(FVCtx*);
1396: PetscFunctionListFind(physics,physname,&r);
1398: /* Create the physics, will set the number of fields and their names */
1399: (*r)(&ctx);
1400: }
1402: /* Create a DMDA to manage the parallel grid */
1403: DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);
1404: DMSetFromOptions(da);
1405: DMSetUp(da);
1406: /* Inform the DMDA of the field names provided by the physics. */
1407: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1408: for (i=0; i<ctx.physics.dof; i++) {
1409: DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
1410: }
1411: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1412: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1414: /* Set coordinates of cell centers */
1415: DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);
1417: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1418: PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope);
1419: PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds);
1421: /* Create a vector to store the solution and to save the initial state */
1422: DMCreateGlobalVector(da,&X);
1423: VecDuplicate(X,&X0);
1424: VecDuplicate(X,&R);
1426: DMCreateMatrix(da,&B);
1428: /* Create a time-stepping object */
1429: TSCreate(comm,&ts);
1430: TSSetDM(ts,da);
1431: TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
1432: TSSetIJacobian(ts,B,B,FVIJacobian,&ctx);
1433: TSSetType(ts,TSSSP);
1434: TSSetMaxTime(ts,10);
1435: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
1437: /* Compute initial conditions and starting time step */
1438: FVSample(&ctx,da,0,X0);
1439: FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
1440: VecCopy(X0,X); /* The function value was not used so we set X=X0 again */
1441: TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);
1442: TSSetFromOptions(ts); /* Take runtime options */
1443: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1444: {
1445: PetscReal nrm1,nrmsup;
1446: PetscInt steps;
1448: TSSolve(ts,X);
1449: TSGetSolveTime(ts,&ptime);
1450: TSGetStepNumber(ts,&steps);
1452: PetscPrintf(comm,"Final time %8.5f, steps %D\n",(double)ptime,steps);
1453: if (ctx.exact) {
1454: SolutionErrorNorms(&ctx,da,ptime,X,&nrm1,&nrmsup);
1455: PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e ||x-x_e||_sup %8.4e\n",(double)nrm1,(double)nrmsup);
1456: }
1457: }
1459: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1460: if (draw & 0x1) VecView(X0,PETSC_VIEWER_DRAW_WORLD);
1461: if (draw & 0x2) VecView(X,PETSC_VIEWER_DRAW_WORLD);
1462: if (draw & 0x4) {
1463: Vec Y;
1464: VecDuplicate(X,&Y);
1465: FVSample(&ctx,da,ptime,Y);
1466: VecAYPX(Y,-1,X);
1467: VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1468: VecDestroy(&Y);
1469: }
1471: if (view_final) {
1472: PetscViewer viewer;
1473: PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1474: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1475: VecView(X,viewer);
1476: PetscViewerPopFormat(viewer);
1477: PetscViewerDestroy(&viewer);
1478: }
1480: /* Clean up */
1481: (*ctx.physics.destroy)(ctx.physics.user);
1482: for (i=0; i<ctx.physics.dof; i++) PetscFree(ctx.physics.fieldname[i]);
1483: PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
1484: PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
1485: VecDestroy(&X);
1486: VecDestroy(&X0);
1487: VecDestroy(&R);
1488: MatDestroy(&B);
1489: DMDestroy(&da);
1490: TSDestroy(&ts);
1491: PetscFunctionListDestroy(&limiters);
1492: PetscFunctionListDestroy(&physics);
1493: PetscFinalize();
1494: return 0;
1495: }
1497: /*TEST
1499: build:
1500: requires: !complex
1502: test:
1503: args: -da_grid_x 100 -initial 1 -xmin -2 -xmax 5 -exact -limit mc
1504: requires: !complex !single
1506: test:
1507: suffix: 2
1508: args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1
1509: filter: sed "s/at 48/at 0/g"
1510: requires: !complex !single
1512: test:
1513: suffix: 3
1514: args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1
1515: nsize: 3
1516: filter: sed "s/at 48/at 0/g"
1517: requires: !complex !single
1519: TEST*/