Actual source code: ex9.c
petsc-3.6.1 2015-08-06
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> /* For the Kernel_*_gets_* stuff for BAIJ */
40: PETSC_STATIC_INLINE PetscReal Sgn(PetscReal a) { return (a<0) ? -1 : 1; }
41: PETSC_STATIC_INLINE PetscReal Abs(PetscReal a) { return (a<0) ? 0 : a; }
42: PETSC_STATIC_INLINE PetscReal Sqr(PetscReal a) { return a*a; }
43: PETSC_STATIC_INLINE PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; }
44: PETSC_UNUSED PETSC_STATIC_INLINE PetscReal MinAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) < PetscAbs(b)) ? a : b; }
45: PETSC_STATIC_INLINE PetscReal MinMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscAbs(b)); }
46: PETSC_STATIC_INLINE PetscReal MaxMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMax(PetscAbs(a),PetscAbs(b)); }
47: PETSC_STATIC_INLINE PetscReal MinMod3(PetscReal a,PetscReal b,PetscReal c) {return (a*b<0 || a*c<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscMin(PetscAbs(b),PetscAbs(c))); }
49: PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin + fmod(range+fmod(a,range),range); }
52: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
53: typedef struct _LimitInfo {
54: PetscReal hx;
55: PetscInt m;
56: } *LimitInfo;
57: static void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
58: {
59: PetscInt i;
60: for (i=0; i<info->m; i++) lmt[i] = 0;
61: }
62: static void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
63: {
64: PetscInt i;
65: for (i=0; i<info->m; i++) lmt[i] = jR[i];
66: }
67: static void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
68: {
69: PetscInt i;
70: for (i=0; i<info->m; i++) lmt[i] = jL[i];
71: }
72: static void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
73: {
74: PetscInt i;
75: for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i] + jR[i]);
76: }
77: static void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
78: {
79: PetscInt i;
80: for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]);
81: }
82: static void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
83: {
84: PetscInt i;
85: for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]));
86: }
87: static void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
88: {
89: PetscInt i;
90: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]);
91: }
92: static void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
93: { /* phi = (t + abs(t)) / (1 + abs(t)) */
94: PetscInt i;
95: 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);
96: }
97: static void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
98: { /* phi = (t + t^2) / (1 + t^2) */
99: PetscInt i;
100: 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);
101: }
102: static void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
103: { /* phi = (t + t^2) / (1 + t^2) */
104: PetscInt i;
105: for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0
106: : (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
107: }
108: static void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
109: { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
110: PetscInt i;
111: for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i]) + 2*Sqr(jL[i])*jR[i])
112: / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
113: }
114: static void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
115: { /* Symmetric version of above */
116: PetscInt i;
117: for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i])
118: / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
119: }
120: static void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
121: { /* Eq 11 of Cada-Torrilhon 2009 */
122: PetscInt i;
123: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]);
124: }
126: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)
127: {
128: return PetscMax(0,PetscMin((L+2*R)/3,
129: PetscMax(-0.5*L,
130: PetscMin(2*L,
131: PetscMin((L+2*R)/3,1.6*R)))));
132: }
133: static void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
134: { /* Cada-Torrilhon 2009, Eq 13 */
135: PetscInt i;
136: for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]);
137: }
138: static void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
139: { /* Cada-Torrilhon 2009, Eq 22 */
140: /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
141: const PetscReal eps = 1e-7,hx = info->hx;
142: PetscInt i;
143: for (i=0; i<info->m; i++) {
144: const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r*hx);
145: lmt[i] = ((eta < 1-eps)
146: ? (jL[i] + 2*jR[i]) / 3
147: : ((eta > 1+eps)
148: ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i])
149: : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3
150: + (1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]))));
151: }
152: }
153: static void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
154: { Limit_CadaTorrilhon3R(0.1,info,jL,jR,lmt); }
155: static void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
156: { Limit_CadaTorrilhon3R(1,info,jL,jR,lmt); }
157: static void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
158: { Limit_CadaTorrilhon3R(10,info,jL,jR,lmt); }
159: static void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
160: { Limit_CadaTorrilhon3R(100,info,jL,jR,lmt); }
163: /* --------------------------------- Finite Volume data structures ----------------------------------- */
165: typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType;
166: static const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0};
167: typedef PetscErrorCode (*RiemannFunction)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*);
168: typedef PetscErrorCode (*ReconstructFunction)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*,PetscReal*);
170: typedef struct {
171: PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*);
172: RiemannFunction riemann;
173: ReconstructFunction characteristic;
174: PetscErrorCode (*destroy)(void*);
175: void *user;
176: PetscInt dof;
177: char *fieldname[16];
178: } PhysicsCtx;
180: typedef struct {
181: void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*);
182: PhysicsCtx physics;
184: MPI_Comm comm;
185: char prefix[256];
187: /* Local work arrays */
188: PetscScalar *R,*Rinv; /* Characteristic basis, and it's inverse. COLUMN-MAJOR */
189: PetscScalar *cjmpLR; /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
190: PetscScalar *cslope; /* Limited slope, written in characteristic basis */
191: PetscScalar *uLR; /* Solution at left and right of interface, conservative variables, len=2*dof */
192: PetscScalar *flux; /* Flux across interface */
193: PetscReal *speeds; /* Speeds of each wave */
195: PetscReal cfl_idt; /* Max allowable value of 1/Delta t */
196: PetscReal cfl;
197: PetscReal xmin,xmax;
198: PetscInt initial;
199: PetscBool exact;
200: FVBCType bctype;
201: } FVCtx;
204: /* Utility */
208: PetscErrorCode RiemannListAdd(PetscFunctionList *flist,const char *name,RiemannFunction rsolve)
209: {
213: PetscFunctionListAdd(flist,name,rsolve);
214: return(0);
215: }
219: PetscErrorCode RiemannListFind(PetscFunctionList flist,const char *name,RiemannFunction *rsolve)
220: {
224: PetscFunctionListFind(flist,name,rsolve);
225: if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,1,"Riemann solver \"%s\" could not be found",name);
226: return(0);
227: }
231: PetscErrorCode ReconstructListAdd(PetscFunctionList *flist,const char *name,ReconstructFunction r)
232: {
236: PetscFunctionListAdd(flist,name,r);
237: return(0);
238: }
242: PetscErrorCode ReconstructListFind(PetscFunctionList flist,const char *name,ReconstructFunction *r)
243: {
247: PetscFunctionListFind(flist,name,r);
248: if (!*r) SETERRQ1(PETSC_COMM_SELF,1,"Reconstruction \"%s\" could not be found",name);
249: return(0);
250: }
253: /* --------------------------------- Physics ----------------------------------- */
254: /**
255: * Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction. These
256: * are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the
257: * number of fields and their names, and a function to deallocate private storage.
258: **/
260: /* First a few functions useful to several different physics */
263: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
264: {
265: PetscInt i,j;
268: for (i=0; i<m; i++) {
269: for (j=0; j<m; j++) Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j);
270: speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
271: }
272: return(0);
273: }
277: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
278: {
282: PetscFree(vctx);
283: return(0);
284: }
288: /* --------------------------------- Advection ----------------------------------- */
290: typedef struct {
291: PetscReal a; /* advective velocity */
292: } AdvectCtx;
296: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
297: {
298: AdvectCtx *ctx = (AdvectCtx*)vctx;
299: PetscReal speed;
302: speed = ctx->a;
303: flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
304: *maxspeed = speed;
305: return(0);
306: }
310: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
311: {
312: AdvectCtx *ctx = (AdvectCtx*)vctx;
315: X[0] = 1.;
316: Xi[0] = 1.;
317: speeds[0] = ctx->a;
318: return(0);
319: }
323: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
324: {
325: AdvectCtx *ctx = (AdvectCtx*)vctx;
326: PetscReal a = ctx->a,x0;
329: switch (bctype) {
330: case FVBC_OUTFLOW: x0 = x-a*t; break;
331: case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
332: default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
333: }
334: switch (initial) {
335: case 0: u[0] = (x0 < 0) ? 1 : -1; break;
336: case 1: u[0] = (x0 < 0) ? -1 : 1; break;
337: case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
338: case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break;
339: case 4: u[0] = PetscAbs(x0); break;
340: case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break;
341: case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
342: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
343: }
344: return(0);
345: }
349: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
350: {
352: AdvectCtx *user;
355: PetscNew(&user);
356: ctx->physics.sample = PhysicsSample_Advect;
357: ctx->physics.riemann = PhysicsRiemann_Advect;
358: ctx->physics.characteristic = PhysicsCharacteristic_Advect;
359: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
360: ctx->physics.user = user;
361: ctx->physics.dof = 1;
362: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
363: user->a = 1;
364: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
365: {
366: PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
367: }
368: PetscOptionsEnd();
369: return(0);
370: }
374: /* --------------------------------- Burgers ----------------------------------- */
376: typedef struct {
377: PetscReal lxf_speed;
378: } BurgersCtx;
382: static PetscErrorCode PhysicsSample_Burgers(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
383: {
386: if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
387: switch (initial) {
388: case 0: u[0] = (x < 0) ? 1 : -1; break;
389: case 1:
390: if (x < -t) u[0] = -1;
391: else if (x < t) u[0] = x/t;
392: else u[0] = 1;
393: break;
394: case 2:
395: if (x < 0) u[0] = 0;
396: else if (x <= t) u[0] = x/t;
397: else if (x < 1+0.5*t) u[0] = 1;
398: else u[0] = 0;
399: break;
400: case 3:
401: if (x < 0.2*t) u[0] = 0.2;
402: else if (x < t) u[0] = x/t;
403: else u[0] = 1;
404: break;
405: case 4:
406: if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
407: u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin)));
408: break;
409: case 5: /* Pure shock solution */
410: if (x < 0.5*t) u[0] = 1;
411: else u[0] = 0;
412: break;
413: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
414: }
415: return(0);
416: }
420: static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
421: {
424: if (uL[0] < uR[0]) { /* rarefaction */
425: flux[0] = (uL[0]*uR[0] < 0)
426: ? 0 /* sonic rarefaction */
427: : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0]));
428: } else { /* shock */
429: flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0]));
430: }
431: *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
432: return(0);
433: }
437: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
438: {
439: PetscReal speed;
442: speed = 0.5*(uL[0] + uR[0]);
443: flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
444: if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
445: *maxspeed = speed;
446: return(0);
447: }
451: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
452: {
453: PetscReal c;
454: PetscScalar fL,fR;
457: c = ((BurgersCtx*)vctx)->lxf_speed;
458: fL = 0.5*PetscSqr(uL[0]);
459: fR = 0.5*PetscSqr(uR[0]);
460: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
461: *maxspeed = c;
462: return(0);
463: }
467: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
468: {
469: PetscReal c;
470: PetscScalar fL,fR;
473: c = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0]));
474: fL = 0.5*PetscSqr(uL[0]);
475: fR = 0.5*PetscSqr(uR[0]);
476: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
477: *maxspeed = c;
478: return(0);
479: }
483: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
484: {
485: BurgersCtx *user;
486: PetscErrorCode ierr;
487: RiemannFunction r;
488: PetscFunctionList rlist = 0;
489: char rname[256] = "exact";
492: PetscNew(&user);
494: ctx->physics.sample = PhysicsSample_Burgers;
495: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
496: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
497: ctx->physics.user = user;
498: ctx->physics.dof = 1;
500: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
501: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Burgers_Exact);
502: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Burgers_Roe);
503: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Burgers_LxF);
504: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Burgers_Rusanov);
505: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
506: {
507: PetscOptionsFList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
508: }
509: PetscOptionsEnd();
510: RiemannListFind(rlist,rname,&r);
511: PetscFunctionListDestroy(&rlist);
512: ctx->physics.riemann = r;
514: /* *
515: * Hack to deal with LxF in semi-discrete form
516: * max speed is 1 for the basic initial conditions (where |u| <= 1)
517: * */
518: if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
519: return(0);
520: }
524: /* --------------------------------- Traffic ----------------------------------- */
526: typedef struct {
527: PetscReal lxf_speed;
528: PetscReal a;
529: } TrafficCtx;
531: PETSC_STATIC_INLINE PetscScalar TrafficFlux(PetscScalar a,PetscScalar u) { return a*u*(1-u); }
535: static PetscErrorCode PhysicsSample_Traffic(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
536: {
537: PetscReal a = ((TrafficCtx*)vctx)->a;
540: if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
541: switch (initial) {
542: case 0:
543: u[0] = (-a*t < x) ? 2 : 0; break;
544: case 1:
545: if (x < PetscMin(2*a*t,0.5+a*t)) u[0] = -1;
546: else if (x < 1) u[0] = 0;
547: else u[0] = 1;
548: break;
549: case 2:
550: if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
551: u[0] = 0.7 + 0.3*PetscSinReal(2*PETSC_PI*((x-xmin)/(xmax-xmin)));
552: break;
553: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
554: }
555: return(0);
556: }
560: static PetscErrorCode PhysicsRiemann_Traffic_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
561: {
562: PetscReal a = ((TrafficCtx*)vctx)->a;
565: if (uL[0] < uR[0]) {
566: flux[0] = PetscMin(TrafficFlux(a,uL[0]),
567: TrafficFlux(a,uR[0]));
568: } else {
569: flux[0] = (uR[0] < 0.5 && 0.5 < uL[0])
570: ? TrafficFlux(a,0.5)
571: : PetscMax(TrafficFlux(a,uL[0]),
572: TrafficFlux(a,uR[0]));
573: }
574: *maxspeed = a*MaxAbs(1-2*uL[0],1-2*uR[0]);
575: return(0);
576: }
580: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
581: {
582: PetscReal a = ((TrafficCtx*)vctx)->a;
583: PetscReal speed;
586: speed = a*(1 - (uL[0] + uR[0]));
587: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
588: *maxspeed = speed;
589: return(0);
590: }
594: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
595: {
596: TrafficCtx *phys = (TrafficCtx*)vctx;
597: PetscReal a = phys->a;
598: PetscReal speed;
601: speed = a*(1 - (uL[0] + uR[0]));
602: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]);
603: *maxspeed = speed;
604: return(0);
605: }
609: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
610: {
611: PetscReal a = ((TrafficCtx*)vctx)->a;
612: PetscReal speed;
615: speed = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0]));
616: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]);
617: *maxspeed = speed;
618: return(0);
619: }
623: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
624: {
625: PetscErrorCode ierr;
626: TrafficCtx *user;
627: RiemannFunction r;
628: PetscFunctionList rlist = 0;
629: char rname[256] = "exact";
632: PetscNew(&user);
633: ctx->physics.sample = PhysicsSample_Traffic;
634: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
635: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
636: ctx->physics.user = user;
637: ctx->physics.dof = 1;
639: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
640: user->a = 0.5;
641: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Traffic_Exact);
642: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Traffic_Roe);
643: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Traffic_LxF);
644: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Traffic_Rusanov);
645: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");
646: {
647: PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,NULL);
648: PetscOptionsFList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
649: }
650: PetscOptionsEnd();
652: RiemannListFind(rlist,rname,&r);
653: PetscFunctionListDestroy(&rlist);
655: ctx->physics.riemann = r;
657: /* *
658: * Hack to deal with LxF in semi-discrete form
659: * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
660: * */
661: if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3*user->a;
662: return(0);
663: }
665: /* --------------------------------- Linear Acoustics ----------------------------------- */
667: /* Flux: u_t + (A u)_x
668: * z = sqrt(rho*bulk), c = sqrt(rho/bulk)
669: * Spectral decomposition: A = R * D * Rinv
670: * [ cz] = [-z z] [-c ] [-1/2z 1/2]
671: * [c/z ] = [ 1 1] [ c] [ 1/2z 1/2]
672: *
673: * We decompose this into the left-traveling waves Al = R * D^- Rinv
674: * and the right-traveling waves Ar = R * D^+ * Rinv
675: * Multiplying out these expressions produces the following two matrices
676: */
678: typedef struct {
679: PetscReal c; /* speed of sound: c = sqrt(bulk/rho) */
680: PetscReal z; /* impedence: z = sqrt(rho*bulk) */
681: } AcousticsCtx;
683: PETSC_UNUSED PETSC_STATIC_INLINE void AcousticsFlux(AcousticsCtx *ctx,const PetscScalar *u,PetscScalar *f)
684: {
685: f[0] = ctx->c*ctx->z*u[1];
686: f[1] = ctx->c/ctx->z*u[0];
687: }
691: static PetscErrorCode PhysicsCharacteristic_Acoustics(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
692: {
693: AcousticsCtx *phys = (AcousticsCtx*)vctx;
694: PetscReal z = phys->z,c = phys->c;
697: X[0*2+0] = -z;
698: X[0*2+1] = z;
699: X[1*2+0] = 1;
700: X[1*2+1] = 1;
701: Xi[0*2+0] = -1./(2*z);
702: Xi[0*2+1] = 1./2;
703: Xi[1*2+0] = 1./(2*z);
704: Xi[1*2+1] = 1./2;
705: speeds[0] = -c;
706: speeds[1] = c;
707: return(0);
708: }
712: static PetscErrorCode PhysicsSample_Acoustics_Initial(AcousticsCtx *phys,PetscInt initial,PetscReal xmin,PetscReal xmax,PetscReal x,PetscReal *u)
713: {
715: switch (initial) {
716: case 0:
717: u[0] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.2) < 0.1) ? 1 : 0.5;
718: u[1] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.7) < 0.1) ? 1 : -0.5;
719: break;
720: case 1:
721: u[0] = PetscCosReal(3 * 2*PETSC_PI*x/(xmax-xmin));
722: u[1] = PetscExpReal(-PetscSqr(x - (xmax + xmin)/2) / (2*PetscSqr(0.2*(xmax - xmin)))) - 0.5;
723: break;
724: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
725: }
726: return(0);
727: }
731: static PetscErrorCode PhysicsSample_Acoustics(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
732: {
733: AcousticsCtx *phys = (AcousticsCtx*)vctx;
734: PetscReal c = phys->c;
735: PetscReal x0a,x0b,u0a[2],u0b[2],tmp[2];
736: PetscReal X[2][2],Xi[2][2],dummy[2];
740: switch (bctype) {
741: case FVBC_OUTFLOW:
742: x0a = x+c*t;
743: x0b = x-c*t;
744: break;
745: case FVBC_PERIODIC:
746: x0a = RangeMod(x+c*t,xmin,xmax);
747: x0b = RangeMod(x-c*t,xmin,xmax);
748: break;
749: default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
750: }
751: PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0a,u0a);
752: PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0b,u0b);
753: PhysicsCharacteristic_Acoustics(vctx,2,u,&X[0][0],&Xi[0][0],dummy);
754: tmp[0] = Xi[0][0]*u0a[0] + Xi[0][1]*u0a[1];
755: tmp[1] = Xi[1][0]*u0b[0] + Xi[1][1]*u0b[1];
756: u[0] = X[0][0]*tmp[0] + X[0][1]*tmp[1];
757: u[1] = X[1][0]*tmp[0] + X[1][1]*tmp[1];
758: return(0);
759: }
763: static PetscErrorCode PhysicsRiemann_Acoustics_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
764: {
765: AcousticsCtx *phys = (AcousticsCtx*)vctx;
766: PetscReal c = phys->c,z = phys->z;
767: PetscReal
768: Al[2][2] = {{-c/2 , c*z/2 },
769: {c/(2*z) , -c/2 }}, /* Left traveling waves */
770: Ar[2][2] = {{c/2 , c*z/2 },
771: {c/(2*z) , c/2 }}; /* Right traveling waves */
774: flux[0] = Al[0][0]*uR[0] + Al[0][1]*uR[1] + Ar[0][0]*uL[0] + Ar[0][1]*uL[1];
775: flux[1] = Al[1][0]*uR[0] + Al[1][1]*uR[1] + Ar[1][0]*uL[0] + Ar[1][1]*uL[1];
776: *maxspeed = c;
777: return(0);
778: }
782: static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx)
783: {
784: PetscErrorCode ierr;
785: AcousticsCtx *user;
786: PetscFunctionList rlist = 0,rclist = 0;
787: char rname[256] = "exact",rcname[256] = "characteristic";
790: PetscNew(&user);
791: ctx->physics.sample = PhysicsSample_Acoustics;
792: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
793: ctx->physics.user = user;
794: ctx->physics.dof = 2;
796: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
797: PetscStrallocpy("v",&ctx->physics.fieldname[1]);
799: user->c = 1;
800: 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,NULL);
808: PetscOptionsReal("-physics_acoustics_z","z = sqrt(bulk*rho)","",user->z,&user->z,NULL);
809: PetscOptionsFList("-physics_acoustics_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
810: PetscOptionsFList("-physics_acoustics_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
811: }
812: PetscOptionsEnd();
813: RiemannListFind(rlist,rname,&ctx->physics.riemann);
814: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
815: PetscFunctionListDestroy(&rlist);
816: PetscFunctionListDestroy(&rclist);
817: return(0);
818: }
820: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */
822: typedef struct {
823: PetscReal acoustic_speed;
824: } IsoGasCtx;
826: PETSC_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*PetscSinReal(2*PETSC_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;
908: PetscErrorCode ierr;
911: if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
912: {
913: /* Solve for star state */
914: PetscScalar res,tmp,rho = 0.5*(L.rho + R.rho); /* initial guess */
915: for (i=0; i<20; i++) {
916: PetscScalar fr,fl,dfr,dfl;
917: fl = (L.rho < rho)
918: ? (rho-L.rho)/PetscSqrtScalar(L.rho*rho) /* shock */
919: : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
920: fr = (R.rho < rho)
921: ? (rho-R.rho)/PetscSqrtScalar(R.rho*rho) /* shock */
922: : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
923: res = R.u-L.u + c*(fr+fl);
924: PetscIsInfOrNanScalar(res);
925: if (PetscAbsScalar(res) < 1e-10) {
926: star.rho = rho;
927: star.u = L.u - c*fl;
928: goto converged;
929: }
930: dfl = (L.rho < rho)
931: ? 1/PetscSqrtScalar(L.rho*rho)*(1 - 0.5*(rho-L.rho)/rho)
932: : 1/rho;
933: dfr = (R.rho < rho)
934: ? 1/PetscSqrtScalar(R.rho*rho)*(1 - 0.5*(rho-R.rho)/rho)
935: : 1/rho;
936: tmp = rho - res/(c*(dfr+dfl));
937: if (tmp <= 0) rho /= 2; /* Guard against Newton shooting off to a negative density */
938: else rho = tmp;
939: if (!((rho > 0) && PetscIsNormalScalar(rho))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate rho=%g",(double)PetscRealPart(rho));
940: }
941: SETERRQ1(PETSC_COMM_SELF,1,"Newton iteration for star.rho diverged after %D iterations",i);
942: }
943: converged:
944: if (L.u-c < 0 && 0 < star.u-c) { /* 1-wave is sonic rarefaction */
945: PetscScalar ufan[2];
946: ufan[0] = L.rho*PetscExpScalar(L.u/c - 1);
947: ufan[1] = c*ufan[0];
948: IsoGasFlux(c,ufan,flux);
949: } else if (star.u+c < 0 && 0 < R.u+c) { /* 2-wave is sonic rarefaction */
950: PetscScalar ufan[2];
951: ufan[0] = R.rho*PetscExpScalar(-R.u/c - 1);
952: ufan[1] = -c*ufan[0];
953: IsoGasFlux(c,ufan,flux);
954: } else if ((L.rho >= star.rho && L.u-c >= 0)
955: || (L.rho < star.rho && (star.rho*star.u-L.rho*L.u)/(star.rho-L.rho) > 0)) {
956: /* 1-wave is supersonic rarefaction, or supersonic shock */
957: IsoGasFlux(c,uL,flux);
958: } else if ((star.rho <= R.rho && R.u+c <= 0)
959: || (star.rho > R.rho && (R.rho*R.u-star.rho*star.u)/(R.rho-star.rho) < 0)) {
960: /* 2-wave is supersonic rarefaction or supersonic shock */
961: IsoGasFlux(c,uR,flux);
962: } else {
963: ustar[0] = star.rho;
964: ustar[1] = star.rho*star.u;
965: IsoGasFlux(c,ustar,flux);
966: }
967: *maxspeed = MaxAbs(MaxAbs(star.u-c,star.u+c),MaxAbs(L.u-c,R.u+c));
968: return(0);
969: }
973: static PetscErrorCode PhysicsRiemann_IsoGas_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
974: {
975: IsoGasCtx *phys = (IsoGasCtx*)vctx;
976: PetscScalar c = phys->acoustic_speed,fL[2],fR[2],s;
977: struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
980: if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
981: IsoGasFlux(c,uL,fL);
982: IsoGasFlux(c,uR,fR);
983: s = PetscMax(PetscAbs(L.u),PetscAbs(R.u))+c;
984: flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
985: flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
986: *maxspeed = s;
987: return(0);
988: }
992: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
993: {
994: IsoGasCtx *phys = (IsoGasCtx*)vctx;
995: PetscReal c = phys->acoustic_speed;
999: speeds[0] = u[1]/u[0] - c;
1000: speeds[1] = u[1]/u[0] + c;
1001: X[0*2+0] = 1;
1002: X[0*2+1] = speeds[0];
1003: X[1*2+0] = 1;
1004: X[1*2+1] = speeds[1];
1006: PetscMemcpy(Xi,X,4*sizeof(X[0]));
1007: PetscKernel_A_gets_inverse_A_2(Xi,0);
1008: return(0);
1009: }
1013: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
1014: {
1015: PetscErrorCode ierr;
1016: IsoGasCtx *user;
1017: PetscFunctionList rlist = 0,rclist = 0;
1018: char rname[256] = "exact",rcname[256] = "characteristic";
1021: PetscNew(&user);
1022: ctx->physics.sample = PhysicsSample_IsoGas;
1023: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
1024: ctx->physics.user = user;
1025: ctx->physics.dof = 2;
1027: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1028: PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1030: user->acoustic_speed = 1;
1032: RiemannListAdd(&rlist,"exact", PhysicsRiemann_IsoGas_Exact);
1033: RiemannListAdd(&rlist,"roe", PhysicsRiemann_IsoGas_Roe);
1034: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_IsoGas_Rusanov);
1035: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_IsoGas);
1036: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1037: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");
1038: {
1039: PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,NULL);
1040: PetscOptionsFList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
1041: PetscOptionsFList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
1042: }
1043: PetscOptionsEnd();
1044: RiemannListFind(rlist,rname,&ctx->physics.riemann);
1045: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1046: PetscFunctionListDestroy(&rlist);
1047: PetscFunctionListDestroy(&rclist);
1048: return(0);
1049: }
1053: /* --------------------------------- Shallow Water ----------------------------------- */
1055: typedef struct {
1056: PetscReal gravity;
1057: } ShallowCtx;
1059: PETSC_STATIC_INLINE void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
1060: {
1061: f[0] = u[1];
1062: f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
1063: }
1067: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
1068: {
1069: ShallowCtx *phys = (ShallowCtx*)vctx;
1070: PetscScalar g = phys->gravity,ustar[2],cL,cR,c,cstar;
1071: struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
1072: PetscInt i;
1073: PetscErrorCode ierr;
1076: if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed thickness is negative");
1077: cL = PetscSqrtScalar(g*L.h);
1078: cR = PetscSqrtScalar(g*R.h);
1079: c = PetscMax(cL,cR);
1080: {
1081: /* Solve for star state */
1082: const PetscInt maxits = 50;
1083: PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */
1084: h0 = h;
1085: for (i=0; i<maxits; i++) {
1086: PetscScalar fr,fl,dfr,dfl;
1087: fl = (L.h < h)
1088: ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */
1089: : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h); /* rarefaction */
1090: fr = (R.h < h)
1091: ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */
1092: : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h); /* rarefaction */
1093: res = R.u - L.u + fr + fl;
1094: PetscIsInfOrNanScalar(res);
1095: if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h-h0) < 1e-8)) {
1096: star.h = h;
1097: star.u = L.u - fl;
1098: goto converged;
1099: } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */
1100: h = 0.8*h0 + 0.2*h;
1101: continue;
1102: }
1103: /* Accept the last step and take another */
1104: res0 = res;
1105: h0 = h;
1106: dfl = (L.h < h)
1107: ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h)
1108: : PetscSqrtScalar(g/h);
1109: dfr = (R.h < h)
1110: ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h)
1111: : PetscSqrtScalar(g/h);
1112: tmp = h - res/(dfr+dfl);
1113: if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */
1114: else h = tmp;
1115: if (!((h > 0) && PetscIsNormalScalar(h))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate h=%g",(double)h);
1116: }
1117: SETERRQ1(PETSC_COMM_SELF,1,"Newton iteration for star.h diverged after %D iterations",i);
1118: }
1119: converged:
1120: cstar = PetscSqrtScalar(g*star.h);
1121: if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */
1122: PetscScalar ufan[2];
1123: ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL);
1124: ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0];
1125: ShallowFlux(phys,ufan,flux);
1126: } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */
1127: PetscScalar ufan[2];
1128: ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR);
1129: ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0];
1130: ShallowFlux(phys,ufan,flux);
1131: } else if ((L.h >= star.h && L.u-c >= 0)
1132: || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) {
1133: /* 1-wave is right-travelling shock (supersonic) */
1134: ShallowFlux(phys,uL,flux);
1135: } else if ((star.h <= R.h && R.u+c <= 0)
1136: || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) {
1137: /* 2-wave is left-travelling shock (supersonic) */
1138: ShallowFlux(phys,uR,flux);
1139: } else {
1140: ustar[0] = star.h;
1141: ustar[1] = star.h*star.u;
1142: ShallowFlux(phys,ustar,flux);
1143: }
1144: *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR));
1145: return(0);
1146: }
1150: static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
1151: {
1152: ShallowCtx *phys = (ShallowCtx*)vctx;
1153: PetscScalar g = phys->gravity,fL[2],fR[2],s;
1154: struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
1157: if (!(L.h > 0 && R.h > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed thickness is negative");
1158: ShallowFlux(phys,uL,fL);
1159: ShallowFlux(phys,uR,fR);
1160: s = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h));
1161: flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
1162: flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
1163: *maxspeed = s;
1164: return(0);
1165: }
1169: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
1170: {
1171: ShallowCtx *phys = (ShallowCtx*)vctx;
1172: PetscReal c;
1176: c = PetscSqrtScalar(u[0]*phys->gravity);
1177: speeds[0] = u[1]/u[0] - c;
1178: speeds[1] = u[1]/u[0] + c;
1179: X[0*2+0] = 1;
1180: X[0*2+1] = speeds[0];
1181: X[1*2+0] = 1;
1182: X[1*2+1] = speeds[1];
1184: PetscMemcpy(Xi,X,4*sizeof(X[0]));
1185: PetscKernel_A_gets_inverse_A_2(Xi,0);
1186: return(0);
1187: }
1191: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
1192: {
1193: PetscErrorCode ierr;
1194: ShallowCtx *user;
1195: PetscFunctionList rlist = 0,rclist = 0;
1196: char rname[256] = "exact",rcname[256] = "characteristic";
1199: PetscNew(&user);
1200: /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
1201: ctx->physics.sample = PhysicsSample_IsoGas;
1202: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
1203: ctx->physics.user = user;
1204: ctx->physics.dof = 2;
1206: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1207: PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1209: user->gravity = 1;
1211: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Shallow_Exact);
1212: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov);
1213: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Shallow);
1214: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1215: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow","");
1216: {
1217: PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,NULL);
1218: PetscOptionsFList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
1219: PetscOptionsFList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
1220: }
1221: PetscOptionsEnd();
1222: RiemannListFind(rlist,rname,&ctx->physics.riemann);
1223: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1224: PetscFunctionListDestroy(&rlist);
1225: PetscFunctionListDestroy(&rclist);
1226: return(0);
1227: }
1231: /* --------------------------------- Finite Volume Solver ----------------------------------- */
1235: static PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
1236: {
1237: FVCtx *ctx = (FVCtx*)vctx;
1238: PetscErrorCode ierr;
1239: PetscInt i,j,k,Mx,dof,xs,xm;
1240: PetscReal hx,cfl_idt = 0;
1241: PetscScalar *x,*f,*slope;
1242: Vec Xloc;
1243: DM da;
1246: TSGetDM(ts,&da);
1247: DMGetLocalVector(da,&Xloc);
1248: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1249: hx = (ctx->xmax - ctx->xmin)/Mx;
1250: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1251: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
1253: VecZeroEntries(F);
1255: DMDAVecGetArray(da,Xloc,&x);
1256: DMDAVecGetArray(da,F,&f);
1257: DMDAGetArray(da,PETSC_TRUE,&slope);
1259: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1261: if (ctx->bctype == FVBC_OUTFLOW) {
1262: for (i=xs-2; i<0; i++) {
1263: for (j=0; j<dof; j++) x[i*dof+j] = x[j];
1264: }
1265: for (i=Mx; i<xs+xm+2; i++) {
1266: for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
1267: }
1268: }
1269: for (i=xs-1; i<xs+xm+1; i++) {
1270: struct _LimitInfo info;
1271: PetscScalar *cjmpL,*cjmpR;
1272: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
1273: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1274: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
1275: PetscMemzero(ctx->cjmpLR,2*dof*sizeof(ctx->cjmpLR[0]));
1276: cjmpL = &ctx->cjmpLR[0];
1277: cjmpR = &ctx->cjmpLR[dof];
1278: for (j=0; j<dof; j++) {
1279: PetscScalar jmpL,jmpR;
1280: jmpL = x[(i+0)*dof+j] - x[(i-1)*dof+j];
1281: jmpR = x[(i+1)*dof+j] - x[(i+0)*dof+j];
1282: for (k=0; k<dof; k++) {
1283: cjmpL[k] += ctx->Rinv[k+j*dof] * jmpL;
1284: cjmpR[k] += ctx->Rinv[k+j*dof] * jmpR;
1285: }
1286: }
1287: /* Apply limiter to the left and right characteristic jumps */
1288: info.m = dof;
1289: info.hx = hx;
1290: (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
1291: for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
1292: for (j=0; j<dof; j++) {
1293: PetscScalar tmp = 0;
1294: for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof] * ctx->cslope[k];
1295: slope[i*dof+j] = tmp;
1296: }
1297: }
1299: for (i=xs; i<xs+xm+1; i++) {
1300: PetscReal maxspeed;
1301: PetscScalar *uL,*uR;
1302: uL = &ctx->uLR[0];
1303: uR = &ctx->uLR[dof];
1304: for (j=0; j<dof; j++) {
1305: uL[j] = x[(i-1)*dof+j] + slope[(i-1)*dof+j]*hx/2;
1306: uR[j] = x[(i-0)*dof+j] - slope[(i-0)*dof+j]*hx/2;
1307: }
1308: (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed);
1309: cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
1311: if (i > xs) {
1312: for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
1313: }
1314: if (i < xs+xm) {
1315: for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
1316: }
1317: }
1319: DMDAVecRestoreArray(da,Xloc,&x);
1320: DMDAVecRestoreArray(da,F,&f);
1321: DMDARestoreArray(da,PETSC_TRUE,&slope);
1322: DMRestoreLocalVector(da,&Xloc);
1324: MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));
1325: if (0) {
1326: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
1327: PetscReal dt,tnow;
1328: TSGetTimeStep(ts,&dt);
1329: TSGetTime(ts,&tnow);
1330: if (dt > 0.5/ctx->cfl_idt) {
1331: if (1) {
1332: PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));
1333: } else SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt));
1334: }
1335: }
1336: return(0);
1337: }
1341: static PetscErrorCode SmallMatMultADB(PetscScalar *C,PetscInt bs,const PetscScalar *A,const PetscReal *D,const PetscScalar *B)
1342: {
1343: PetscInt i,j,k;
1346: for (i=0; i<bs; i++) {
1347: for (j=0; j<bs; j++) {
1348: PetscScalar tmp = 0;
1349: for (k=0; k<bs; k++) tmp += A[i*bs+k] * D[k] * B[k*bs+j];
1350: C[i*bs+j] = tmp;
1351: }
1352: }
1353: return(0);
1354: }
1359: static PetscErrorCode FVIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void *vctx)
1360: {
1361: FVCtx *ctx = (FVCtx*)vctx;
1362: PetscErrorCode ierr;
1363: PetscInt i,j,dof = ctx->physics.dof;
1364: PetscScalar *J;
1365: const PetscScalar *x;
1366: PetscReal hx;
1367: DM da;
1368: DMDALocalInfo dainfo;
1371: TSGetDM(ts,&da);
1372: DMDAVecGetArrayRead(da,X,(void*)&x);
1373: DMDAGetLocalInfo(da,&dainfo);
1374: hx = (ctx->xmax - ctx->xmin)/dainfo.mx;
1375: PetscMalloc1(dof*dof,&J);
1376: for (i=dainfo.xs; i<dainfo.xs+dainfo.xm; i++) {
1377: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1378: for (j=0; j<dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]);
1379: SmallMatMultADB(J,dof,ctx->R,ctx->speeds,ctx->Rinv);
1380: for (j=0; j<dof*dof; j++) J[j] = J[j]/hx + shift*(j/dof == j%dof);
1381: MatSetValuesBlocked(B,1,&i,1,&i,J,INSERT_VALUES);
1382: }
1383: PetscFree(J);
1384: DMDAVecRestoreArrayRead(da,X,(void*)&x);
1386: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1387: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1388: if (A != B) {
1389: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1390: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1391: }
1392: return(0);
1393: }
1397: static PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
1398: {
1400: PetscScalar *u,*uj;
1401: PetscInt i,j,k,dof,xs,xm,Mx;
1404: if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,1,"Physics has not provided a sampling function");
1405: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1406: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1407: DMDAVecGetArray(da,U,&u);
1408: PetscMalloc1(dof,&uj);
1409: for (i=xs; i<xs+xm; i++) {
1410: const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h;
1411: const PetscInt N = 200;
1412: /* Integrate over cell i using trapezoid rule with N points. */
1413: for (k=0; k<dof; k++) u[i*dof+k] = 0;
1414: for (j=0; j<N+1; j++) {
1415: PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N;
1416: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
1417: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
1418: }
1419: }
1420: DMDAVecRestoreArray(da,U,&u);
1421: PetscFree(uj);
1422: return(0);
1423: }
1427: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
1428: {
1429: PetscErrorCode ierr;
1430: PetscReal xmin,xmax;
1431: PetscScalar sum,tvsum,tvgsum;
1432: const PetscScalar *x;
1433: PetscInt imin,imax,Mx,i,j,xs,xm,dof;
1434: Vec Xloc;
1435: PetscBool iascii;
1438: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1439: if (iascii) {
1440: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
1441: DMGetLocalVector(da,&Xloc);
1442: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1443: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
1444: DMDAVecGetArrayRead(da,Xloc,(void*)&x);
1445: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1446: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1447: tvsum = 0;
1448: for (i=xs; i<xs+xm; i++) {
1449: for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j] - x[(i-1)*dof+j]);
1450: }
1451: MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));
1452: DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);
1453: DMRestoreLocalVector(da,&Xloc);
1455: VecMin(X,&imin,&xmin);
1456: VecMax(X,&imax,&xmax);
1457: VecSum(X,&sum);
1458: 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));
1459: } else SETERRQ(PETSC_COMM_SELF,1,"Viewer type not supported");
1460: return(0);
1461: }
1465: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1,PetscReal *nrmsup)
1466: {
1468: Vec Y;
1469: PetscInt Mx;
1472: VecGetSize(X,&Mx);
1473: VecDuplicate(X,&Y);
1474: FVSample(ctx,da,t,Y);
1475: VecAYPX(Y,-1,X);
1476: VecNorm(Y,NORM_1,nrm1);
1477: VecNorm(Y,NORM_INFINITY,nrmsup);
1478: *nrm1 /= Mx;
1479: VecDestroy(&Y);
1480: return(0);
1481: }
1485: int main(int argc,char *argv[])
1486: {
1487: char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
1488: PetscFunctionList limiters = 0,physics = 0;
1489: MPI_Comm comm;
1490: TS ts;
1491: DM da;
1492: Vec X,X0,R;
1493: Mat B;
1494: FVCtx ctx;
1495: PetscInt i,dof,xs,xm,Mx,draw = 0;
1496: PetscBool view_final = PETSC_FALSE;
1497: PetscReal ptime;
1498: PetscErrorCode ierr;
1500: PetscInitialize(&argc,&argv,0,help);
1501: comm = PETSC_COMM_WORLD;
1502: PetscMemzero(&ctx,sizeof(ctx));
1504: /* Register limiters to be available on the command line */
1505: PetscFunctionListAdd(&limiters,"upwind" ,Limit_Upwind);
1506: PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit_LaxWendroff);
1507: PetscFunctionListAdd(&limiters,"beam-warming" ,Limit_BeamWarming);
1508: PetscFunctionListAdd(&limiters,"fromm" ,Limit_Fromm);
1509: PetscFunctionListAdd(&limiters,"minmod" ,Limit_Minmod);
1510: PetscFunctionListAdd(&limiters,"superbee" ,Limit_Superbee);
1511: PetscFunctionListAdd(&limiters,"mc" ,Limit_MC);
1512: PetscFunctionListAdd(&limiters,"vanleer" ,Limit_VanLeer);
1513: PetscFunctionListAdd(&limiters,"vanalbada" ,Limit_VanAlbada);
1514: PetscFunctionListAdd(&limiters,"vanalbadatvd" ,Limit_VanAlbadaTVD);
1515: PetscFunctionListAdd(&limiters,"koren" ,Limit_Koren);
1516: PetscFunctionListAdd(&limiters,"korensym" ,Limit_KorenSym);
1517: PetscFunctionListAdd(&limiters,"koren3" ,Limit_Koren3);
1518: PetscFunctionListAdd(&limiters,"cada-torrilhon2" ,Limit_CadaTorrilhon2);
1519: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1);
1520: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1" ,Limit_CadaTorrilhon3R1);
1521: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10);
1522: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100);
1524: /* Register physical models to be available on the command line */
1525: PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect);
1526: PetscFunctionListAdd(&physics,"burgers" ,PhysicsCreate_Burgers);
1527: PetscFunctionListAdd(&physics,"traffic" ,PhysicsCreate_Traffic);
1528: PetscFunctionListAdd(&physics,"acoustics" ,PhysicsCreate_Acoustics);
1529: PetscFunctionListAdd(&physics,"isogas" ,PhysicsCreate_IsoGas);
1530: PetscFunctionListAdd(&physics,"shallow" ,PhysicsCreate_Shallow);
1532: ctx.comm = comm;
1533: ctx.cfl = 0.9; ctx.bctype = FVBC_PERIODIC;
1534: ctx.xmin = -1; ctx.xmax = 1;
1535: PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
1536: {
1537: PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
1538: PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
1539: PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);
1540: PetscOptionsFList("-physics","Name of physics (Riemann solver and characteristics) to use","",physics,physname,physname,sizeof(physname),NULL);
1541: PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
1542: PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
1543: PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
1544: PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
1545: PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
1546: PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
1547: }
1548: PetscOptionsEnd();
1550: /* Choose the limiter from the list of registered limiters */
1551: PetscFunctionListFind(limiters,lname,&ctx.limit);
1552: if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);
1554: /* Choose the physics from the list of registered models */
1555: {
1556: PetscErrorCode (*r)(FVCtx*);
1557: PetscFunctionListFind(physics,physname,&r);
1558: if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
1559: /* Create the physics, will set the number of fields and their names */
1560: (*r)(&ctx);
1561: }
1563: /* Create a DMDA to manage the parallel grid */
1564: DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,-50,ctx.physics.dof,2,NULL,&da);
1565: /* Inform the DMDA of the field names provided by the physics. */
1566: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1567: for (i=0; i<ctx.physics.dof; i++) {
1568: DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
1569: }
1570: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1571: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1573: /* Set coordinates of cell centers */
1574: DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);
1576: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1577: PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope);
1578: PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds);
1580: /* Create a vector to store the solution and to save the initial state */
1581: DMCreateGlobalVector(da,&X);
1582: VecDuplicate(X,&X0);
1583: VecDuplicate(X,&R);
1585: DMCreateMatrix(da,&B);
1587: /* Create a time-stepping object */
1588: TSCreate(comm,&ts);
1589: TSSetDM(ts,da);
1590: TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
1591: TSSetIJacobian(ts,B,B,FVIJacobian,&ctx);
1592: TSSetType(ts,TSSSP);
1593: TSSetDuration(ts,1000,10);
1595: /* Compute initial conditions and starting time step */
1596: FVSample(&ctx,da,0,X0);
1597: FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
1598: VecCopy(X0,X); /* The function value was not used so we set X=X0 again */
1599: TSSetInitialTimeStep(ts,0,ctx.cfl/ctx.cfl_idt);
1601: TSSetFromOptions(ts); /* Take runtime options */
1603: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1605: {
1606: PetscReal nrm1,nrmsup;
1607: PetscInt steps;
1609: TSSolve(ts,X);
1610: TSGetSolveTime(ts,&ptime);
1611: TSGetTimeStepNumber(ts,&steps);
1613: PetscPrintf(comm,"Final time %8.5f, steps %D\n",(double)ptime,steps);
1614: if (ctx.exact) {
1615: SolutionErrorNorms(&ctx,da,ptime,X,&nrm1,&nrmsup);
1616: PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e ||x-x_e||_sup %8.4e\n",(double)nrm1,(double)nrmsup);
1617: }
1618: }
1620: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1622: if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
1623: if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
1624: if (draw & 0x4) {
1625: Vec Y;
1626: VecDuplicate(X,&Y);
1627: FVSample(&ctx,da,ptime,Y);
1628: VecAYPX(Y,-1,X);
1629: VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1630: VecDestroy(&Y);
1631: }
1633: if (view_final) {
1634: PetscViewer viewer;
1635: PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1636: PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1637: VecView(X,viewer);
1638: PetscViewerDestroy(&viewer);
1639: }
1641: /* Clean up */
1642: (*ctx.physics.destroy)(ctx.physics.user);
1643: for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
1644: PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
1645: PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
1646: VecDestroy(&X);
1647: VecDestroy(&X0);
1648: VecDestroy(&R);
1649: MatDestroy(&B);
1650: DMDestroy(&da);
1651: TSDestroy(&ts);
1652: PetscFunctionListDestroy(&limiters);
1653: PetscFunctionListDestroy(&physics);
1654: PetscFinalize();
1655: return 0;
1656: }