Actual source code: ex9.c
petsc-3.4.5 2014-06-29
1: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
2: "Solves scalar and vector problems, choose the physical model with -physics\n"
3: " advection - Constant coefficient scalar advection\n"
4: " u_t + (a*u)_x = 0\n"
5: " burgers - Burgers equation\n"
6: " u_t + (u^2/2)_x = 0\n"
7: " traffic - Traffic equation\n"
8: " u_t + (u*(1-u))_x = 0\n"
9: " acoustics - Acoustic wave propagation\n"
10: " u_t + (c*z*v)_x = 0\n"
11: " v_t + (c/z*u)_x = 0\n"
12: " isogas - Isothermal gas dynamics\n"
13: " rho_t + (rho*u)_x = 0\n"
14: " (rho*u)_t + (rho*u^2 + c^2*rho)_x = 0\n"
15: " shallow - Shallow water equations\n"
16: " h_t + (h*u)_x = 0\n"
17: " (h*u)_t + (h*u^2 + g*h^2/2)_x = 0\n"
18: "Some of these physical models have multiple Riemann solvers, select these with -physics_xxx_riemann\n"
19: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
20: " the states across shocks and rarefactions\n"
21: " roe - Linearized scheme, usually with an entropy fix inside sonic rarefactions\n"
22: "The systems provide a choice of reconstructions with -physics_xxx_reconstruct\n"
23: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
24: " conservative - Limit the conservative variables directly, can cause undesired interaction of waves\n\n"
25: "A variety of limiters for high-resolution TVD limiters are available with -limit\n"
26: " upwind,minmod,superbee,mc,vanleer,vanalbada,koren,cada-torillhon (last two are nominally third order)\n"
27: " and non-TVD schemes lax-wendroff,beam-warming,fromm\n\n"
28: "To preserve the TVD property, one should time step with a strong stability preserving method.\n"
29: "The optimal high order explicit Runge-Kutta methods in TSSSP are recommended for non-stiff problems.\n\n"
30: "Several initial conditions can be chosen with -initial N\n\n"
31: "The problem size should be set with -da_grid_x M\n\n";
33: /* To get isfinite in math.h */
34: #define _XOPEN_SOURCE 600
36: #include <petscts.h>
37: #include <petscdmda.h>
38: #include <petscdraw.h>
40: #include <petsc-private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */
42: PETSC_STATIC_INLINE PetscReal Sgn(PetscReal a) { return (a<0) ? -1 : 1; }
43: PETSC_STATIC_INLINE PetscReal Abs(PetscReal a) { return (a<0) ? 0 : a; }
44: PETSC_STATIC_INLINE PetscReal Sqr(PetscReal a) { return a*a; }
45: PETSC_STATIC_INLINE PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; }
46: PETSC_STATIC_INLINE PetscReal MinAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) < PetscAbs(b)) ? a : b; }
47: PETSC_STATIC_INLINE PetscReal MinMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscAbs(b)); }
48: PETSC_STATIC_INLINE PetscReal MaxMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMax(PetscAbs(a),PetscAbs(b)); }
49: 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))); }
51: PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin + fmod(range+fmod(a,range),range); }
54: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
55: typedef struct _LimitInfo {
56: PetscReal hx;
57: PetscInt m;
58: } *LimitInfo;
59: static void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
60: {
61: PetscInt i;
62: for (i=0; i<info->m; i++) lmt[i] = 0;
63: }
64: static void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
65: {
66: PetscInt i;
67: for (i=0; i<info->m; i++) lmt[i] = jR[i];
68: }
69: static void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
70: {
71: PetscInt i;
72: for (i=0; i<info->m; i++) lmt[i] = jL[i];
73: }
74: static void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
75: {
76: PetscInt i;
77: for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i] + jR[i]);
78: }
79: static void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
80: {
81: PetscInt i;
82: for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]);
83: }
84: static void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
85: {
86: PetscInt i;
87: for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]));
88: }
89: static void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
90: {
91: PetscInt i;
92: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]);
93: }
94: static void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
95: { /* phi = (t + abs(t)) / (1 + abs(t)) */
96: PetscInt i;
97: 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);
98: }
99: static void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
100: { /* phi = (t + t^2) / (1 + t^2) */
101: PetscInt i;
102: 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);
103: }
104: static void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
105: { /* phi = (t + t^2) / (1 + t^2) */
106: PetscInt i;
107: for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0
108: : (jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
109: }
110: static void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
111: { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
112: PetscInt i;
113: for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i]) + 2*Sqr(jL[i])*jR[i])
114: / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
115: }
116: static void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
117: { /* Symmetric version of above */
118: PetscInt i;
119: for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i]) + Sqr(jL[i])*jR[i])
120: / (2*Sqr(jL[i]) - jL[i]*jR[i] + 2*Sqr(jR[i]) + 1e-15));
121: }
122: static void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
123: { /* Eq 11 of Cada-Torrilhon 2009 */
124: PetscInt i;
125: for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]);
126: }
128: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)
129: {
130: return PetscMax(0,PetscMin((L+2*R)/3,
131: PetscMax(-0.5*L,
132: PetscMin(2*L,
133: PetscMin((L+2*R)/3,1.6*R)))));
134: }
135: static void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
136: { /* Cada-Torrilhon 2009, Eq 13 */
137: PetscInt i;
138: for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]);
139: }
140: static void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
141: { /* Cada-Torrilhon 2009, Eq 22 */
142: /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
143: const PetscReal eps = 1e-7,hx = info->hx;
144: PetscInt i;
145: for (i=0; i<info->m; i++) {
146: const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r*hx);
147: lmt[i] = ((eta < 1-eps)
148: ? (jL[i] + 2*jR[i]) / 3
149: : ((eta > 1+eps)
150: ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i])
151: : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3
152: + (1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]))));
153: }
154: }
155: static void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
156: { Limit_CadaTorrilhon3R(0.1,info,jL,jR,lmt); }
157: static void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
158: { Limit_CadaTorrilhon3R(1,info,jL,jR,lmt); }
159: static void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
160: { Limit_CadaTorrilhon3R(10,info,jL,jR,lmt); }
161: static void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
162: { Limit_CadaTorrilhon3R(100,info,jL,jR,lmt); }
165: /* --------------------------------- Finite Volume data structures ----------------------------------- */
167: typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType;
168: static const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0};
169: typedef PetscErrorCode (*RiemannFunction)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*);
170: typedef PetscErrorCode (*ReconstructFunction)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*,PetscReal*);
172: typedef struct {
173: PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*);
174: RiemannFunction riemann;
175: ReconstructFunction characteristic;
176: PetscErrorCode (*destroy)(void*);
177: void *user;
178: PetscInt dof;
179: char *fieldname[16];
180: } PhysicsCtx;
182: typedef struct {
183: void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*);
184: PhysicsCtx physics;
186: MPI_Comm comm;
187: char prefix[256];
189: /* Local work arrays */
190: PetscScalar *R,*Rinv; /* Characteristic basis, and it's inverse. COLUMN-MAJOR */
191: PetscScalar *cjmpLR; /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
192: PetscScalar *cslope; /* Limited slope, written in characteristic basis */
193: PetscScalar *uLR; /* Solution at left and right of interface, conservative variables, len=2*dof */
194: PetscScalar *flux; /* Flux across interface */
195: PetscReal *speeds; /* Speeds of each wave */
197: PetscReal cfl_idt; /* Max allowable value of 1/Delta t */
198: PetscReal cfl;
199: PetscReal xmin,xmax;
200: PetscInt initial;
201: PetscBool exact;
202: FVBCType bctype;
203: } FVCtx;
206: /* Utility */
210: PetscErrorCode RiemannListAdd(PetscFunctionList *flist,const char *name,RiemannFunction rsolve)
211: {
215: PetscFunctionListAdd(flist,name,rsolve);
216: return(0);
217: }
221: PetscErrorCode RiemannListFind(PetscFunctionList flist,const char *name,RiemannFunction *rsolve)
222: {
226: PetscFunctionListFind(flist,name,rsolve);
227: if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,1,"Riemann solver \"%s\" could not be found",name);
228: return(0);
229: }
233: PetscErrorCode ReconstructListAdd(PetscFunctionList *flist,const char *name,ReconstructFunction r)
234: {
238: PetscFunctionListAdd(flist,name,r);
239: return(0);
240: }
244: PetscErrorCode ReconstructListFind(PetscFunctionList flist,const char *name,ReconstructFunction *r)
245: {
249: PetscFunctionListFind(flist,name,r);
250: if (!*r) SETERRQ1(PETSC_COMM_SELF,1,"Reconstruction \"%s\" could not be found",name);
251: return(0);
252: }
255: /* --------------------------------- Physics ----------------------------------- */
256: /**
257: * Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction. These
258: * are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the
259: * number of fields and their names, and a function to deallocate private storage.
260: **/
262: /* First a few functions useful to several different physics */
265: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
266: {
267: PetscInt i,j;
270: for (i=0; i<m; i++) {
271: for (j=0; j<m; j++) Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j);
272: speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
273: }
274: return(0);
275: }
279: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
280: {
284: PetscFree(vctx);
285: return(0);
286: }
290: /* --------------------------------- Advection ----------------------------------- */
292: typedef struct {
293: PetscReal a; /* advective velocity */
294: } AdvectCtx;
298: static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
299: {
300: AdvectCtx *ctx = (AdvectCtx*)vctx;
301: PetscReal speed;
304: speed = ctx->a;
305: flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0];
306: *maxspeed = speed;
307: return(0);
308: }
312: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
313: {
314: AdvectCtx *ctx = (AdvectCtx*)vctx;
317: X[0] = 1.;
318: Xi[0] = 1.;
319: speeds[0] = ctx->a;
320: return(0);
321: }
325: static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
326: {
327: AdvectCtx *ctx = (AdvectCtx*)vctx;
328: PetscReal a = ctx->a,x0;
331: switch (bctype) {
332: case FVBC_OUTFLOW: x0 = x-a*t; break;
333: case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break;
334: default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
335: }
336: switch (initial) {
337: case 0: u[0] = (x0 < 0) ? 1 : -1; break;
338: case 1: u[0] = (x0 < 0) ? -1 : 1; break;
339: case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break;
340: case 3: u[0] = sin(2*M_PI*x0); break;
341: case 4: u[0] = PetscAbs(x0); break;
342: case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(sin(2*M_PI*x0)); break;
343: case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break;
344: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
345: }
346: return(0);
347: }
351: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
352: {
354: AdvectCtx *user;
357: PetscNew(*user,&user);
358: ctx->physics.sample = PhysicsSample_Advect;
359: ctx->physics.riemann = PhysicsRiemann_Advect;
360: ctx->physics.characteristic = PhysicsCharacteristic_Advect;
361: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
362: ctx->physics.user = user;
363: ctx->physics.dof = 1;
364: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
365: user->a = 1;
366: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
367: {
368: PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL);
369: }
370: PetscOptionsEnd();
371: return(0);
372: }
376: /* --------------------------------- Burgers ----------------------------------- */
378: typedef struct {
379: PetscReal lxf_speed;
380: } BurgersCtx;
384: static PetscErrorCode PhysicsSample_Burgers(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
385: {
388: if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
389: switch (initial) {
390: case 0: u[0] = (x < 0) ? 1 : -1; break;
391: case 1:
392: if (x < -t) u[0] = -1;
393: else if (x < t) u[0] = x/t;
394: else u[0] = 1;
395: break;
396: case 2:
397: if (x < 0) u[0] = 0;
398: else if (x <= t) u[0] = x/t;
399: else if (x < 1+0.5*t) u[0] = 1;
400: else u[0] = 0;
401: break;
402: case 3:
403: if (x < 0.2*t) u[0] = 0.2;
404: else if (x < t) u[0] = x/t;
405: else u[0] = 1;
406: break;
407: case 4:
408: if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
409: u[0] = 0.7 + 0.3*sin(2*M_PI*((x-xmin)/(xmax-xmin)));
410: break;
411: case 5: /* Pure shock solution */
412: if (x < 0.5*t) u[0] = 1;
413: else u[0] = 0;
414: break;
415: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
416: }
417: return(0);
418: }
422: static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
423: {
426: if (uL[0] < uR[0]) { /* rarefaction */
427: flux[0] = (uL[0]*uR[0] < 0)
428: ? 0 /* sonic rarefaction */
429: : 0.5*PetscMin(PetscSqr(uL[0]),PetscSqr(uR[0]));
430: } else { /* shock */
431: flux[0] = 0.5*PetscMax(PetscSqr(uL[0]),PetscSqr(uR[0]));
432: }
433: *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
434: return(0);
435: }
439: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
440: {
441: PetscReal speed;
444: speed = 0.5*(uL[0] + uR[0]);
445: flux[0] = 0.25*(PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
446: if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
447: *maxspeed = speed;
448: return(0);
449: }
453: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
454: {
455: PetscReal c;
456: PetscScalar fL,fR;
459: c = ((BurgersCtx*)vctx)->lxf_speed;
460: fL = 0.5*PetscSqr(uL[0]);
461: fR = 0.5*PetscSqr(uR[0]);
462: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
463: *maxspeed = c;
464: return(0);
465: }
469: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
470: {
471: PetscReal c;
472: PetscScalar fL,fR;
475: c = PetscMax(PetscAbs(uL[0]),PetscAbs(uR[0]));
476: fL = 0.5*PetscSqr(uL[0]);
477: fR = 0.5*PetscSqr(uR[0]);
478: flux[0] = 0.5*(fL + fR) - 0.5*c*(uR[0] - uL[0]);
479: *maxspeed = c;
480: return(0);
481: }
485: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
486: {
487: BurgersCtx *user;
488: PetscErrorCode ierr;
489: RiemannFunction r;
490: PetscFunctionList rlist = 0;
491: char rname[256] = "exact";
494: PetscNew(*user,&user);
496: ctx->physics.sample = PhysicsSample_Burgers;
497: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
498: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
499: ctx->physics.user = user;
500: ctx->physics.dof = 1;
502: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
503: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Burgers_Exact);
504: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Burgers_Roe);
505: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Burgers_LxF);
506: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Burgers_Rusanov);
507: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");
508: {
509: PetscOptionsList("-physics_burgers_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
510: }
511: PetscOptionsEnd();
512: RiemannListFind(rlist,rname,&r);
513: PetscFunctionListDestroy(&rlist);
514: ctx->physics.riemann = r;
516: /* *
517: * Hack to deal with LxF in semi-discrete form
518: * max speed is 1 for the basic initial conditions (where |u| <= 1)
519: * */
520: if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
521: return(0);
522: }
526: /* --------------------------------- Traffic ----------------------------------- */
528: typedef struct {
529: PetscReal lxf_speed;
530: PetscReal a;
531: } TrafficCtx;
533: PETSC_STATIC_INLINE PetscScalar TrafficFlux(PetscScalar a,PetscScalar u) { return a*u*(1-u); }
537: static PetscErrorCode PhysicsSample_Traffic(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
538: {
539: PetscReal a = ((TrafficCtx*)vctx)->a;
542: if (bctype == FVBC_PERIODIC && t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solution not implemented for periodic");
543: switch (initial) {
544: case 0:
545: u[0] = (-a*t < x) ? 2 : 0; break;
546: case 1:
547: if (x < PetscMin(2*a*t,0.5+a*t)) u[0] = -1;
548: else if (x < 1) u[0] = 0;
549: else u[0] = 1;
550: break;
551: case 2:
552: if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Only initial condition available");
553: u[0] = 0.7 + 0.3*sin(2*M_PI*((x-xmin)/(xmax-xmin)));
554: break;
555: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
556: }
557: return(0);
558: }
562: static PetscErrorCode PhysicsRiemann_Traffic_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
563: {
564: PetscReal a = ((TrafficCtx*)vctx)->a;
567: if (uL[0] < uR[0]) {
568: flux[0] = PetscMin(TrafficFlux(a,uL[0]),
569: TrafficFlux(a,uR[0]));
570: } else {
571: flux[0] = (uR[0] < 0.5 && 0.5 < uL[0])
572: ? TrafficFlux(a,0.5)
573: : PetscMax(TrafficFlux(a,uL[0]),
574: TrafficFlux(a,uR[0]));
575: }
576: *maxspeed = a*MaxAbs(1-2*uL[0],1-2*uR[0]);
577: return(0);
578: }
582: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
583: {
584: PetscReal a = ((TrafficCtx*)vctx)->a;
585: PetscReal speed;
588: speed = a*(1 - (uL[0] + uR[0]));
589: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*PetscAbs(speed)*(uR[0]-uL[0]);
590: *maxspeed = speed;
591: return(0);
592: }
596: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
597: {
598: TrafficCtx *phys = (TrafficCtx*)vctx;
599: PetscReal a = phys->a;
600: PetscReal speed;
603: speed = a*(1 - (uL[0] + uR[0]));
604: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*phys->lxf_speed*(uR[0]-uL[0]);
605: *maxspeed = speed;
606: return(0);
607: }
611: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
612: {
613: PetscReal a = ((TrafficCtx*)vctx)->a;
614: PetscReal speed;
617: speed = a*PetscMax(PetscAbs(1-2*uL[0]),PetscAbs(1-2*uR[0]));
618: flux[0] = 0.5*(TrafficFlux(a,uL[0]) + TrafficFlux(a,uR[0])) - 0.5*speed*(uR[0]-uL[0]);
619: *maxspeed = speed;
620: return(0);
621: }
625: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
626: {
627: PetscErrorCode ierr;
628: TrafficCtx *user;
629: RiemannFunction r;
630: PetscFunctionList rlist = 0;
631: char rname[256] = "exact";
634: PetscNew(*user,&user);
635: ctx->physics.sample = PhysicsSample_Traffic;
636: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
637: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
638: ctx->physics.user = user;
639: ctx->physics.dof = 1;
641: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
642: user->a = 0.5;
643: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Traffic_Exact);
644: RiemannListAdd(&rlist,"roe", PhysicsRiemann_Traffic_Roe);
645: RiemannListAdd(&rlist,"lxf", PhysicsRiemann_Traffic_LxF);
646: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_Traffic_Rusanov);
647: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Traffic","");
648: {
649: PetscOptionsReal("-physics_traffic_a","Flux = a*u*(1-u)","",user->a,&user->a,NULL);
650: PetscOptionsList("-physics_traffic_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
651: }
652: PetscOptionsEnd();
654: RiemannListFind(rlist,rname,&r);
655: PetscFunctionListDestroy(&rlist);
657: ctx->physics.riemann = r;
659: /* *
660: * Hack to deal with LxF in semi-discrete form
661: * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
662: * */
663: if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3*user->a;
664: return(0);
665: }
667: /* --------------------------------- Linear Acoustics ----------------------------------- */
669: /* Flux: u_t + (A u)_x
670: * z = sqrt(rho*bulk), c = sqrt(rho/bulk)
671: * Spectral decomposition: A = R * D * Rinv
672: * [ cz] = [-z z] [-c ] [-1/2z 1/2]
673: * [c/z ] = [ 1 1] [ c] [ 1/2z 1/2]
674: *
675: * We decompose this into the left-traveling waves Al = R * D^- Rinv
676: * and the right-traveling waves Ar = R * D^+ * Rinv
677: * Multiplying out these expressions produces the following two matrices
678: */
680: typedef struct {
681: PetscReal c; /* speed of sound: c = sqrt(bulk/rho) */
682: PetscReal z; /* impedence: z = sqrt(rho*bulk) */
683: } AcousticsCtx;
685: PETSC_STATIC_INLINE void AcousticsFlux(AcousticsCtx *ctx,const PetscScalar *u,PetscScalar *f)
686: {
687: f[0] = ctx->c*ctx->z*u[1];
688: f[1] = ctx->c/ctx->z*u[0];
689: }
693: static PetscErrorCode PhysicsCharacteristic_Acoustics(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
694: {
695: AcousticsCtx *phys = (AcousticsCtx*)vctx;
696: PetscReal z = phys->z,c = phys->c;
699: X[0*2+0] = -z;
700: X[0*2+1] = z;
701: X[1*2+0] = 1;
702: X[1*2+1] = 1;
703: Xi[0*2+0] = -1./(2*z);
704: Xi[0*2+1] = 1./2;
705: Xi[1*2+0] = 1./(2*z);
706: Xi[1*2+1] = 1./2;
707: speeds[0] = -c;
708: speeds[1] = c;
709: return(0);
710: }
714: static PetscErrorCode PhysicsSample_Acoustics_Initial(AcousticsCtx *phys,PetscInt initial,PetscReal xmin,PetscReal xmax,PetscReal x,PetscReal *u)
715: {
717: switch (initial) {
718: case 0:
719: u[0] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.2) < 0.1) ? 1 : 0.5;
720: u[1] = (PetscAbs((x - xmin)/(xmax - xmin) - 0.7) < 0.1) ? 1 : -0.5;
721: break;
722: case 1:
723: u[0] = cos(3 * 2*PETSC_PI*x/(xmax-xmin));
724: u[1] = exp(-PetscSqr(x - (xmax + xmin)/2) / (2*PetscSqr(0.2*(xmax - xmin)))) - 0.5;
725: break;
726: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
727: }
728: return(0);
729: }
733: static PetscErrorCode PhysicsSample_Acoustics(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
734: {
735: AcousticsCtx *phys = (AcousticsCtx*)vctx;
736: PetscReal c = phys->c;
737: PetscReal x0a,x0b,u0a[2],u0b[2],tmp[2];
738: PetscReal X[2][2],Xi[2][2],dummy[2];
742: switch (bctype) {
743: case FVBC_OUTFLOW:
744: x0a = x+c*t;
745: x0b = x-c*t;
746: break;
747: case FVBC_PERIODIC:
748: x0a = RangeMod(x+c*t,xmin,xmax);
749: x0b = RangeMod(x-c*t,xmin,xmax);
750: break;
751: default: SETERRQ(PETSC_COMM_SELF,1,"unknown BCType");
752: }
753: PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0a,u0a);
754: PhysicsSample_Acoustics_Initial(phys,initial,xmin,xmax,x0b,u0b);
755: PhysicsCharacteristic_Acoustics(vctx,2,u,&X[0][0],&Xi[0][0],dummy);
756: tmp[0] = Xi[0][0]*u0a[0] + Xi[0][1]*u0a[1];
757: tmp[1] = Xi[1][0]*u0b[0] + Xi[1][1]*u0b[1];
758: u[0] = X[0][0]*tmp[0] + X[0][1]*tmp[1];
759: u[1] = X[1][0]*tmp[0] + X[1][1]*tmp[1];
760: return(0);
761: }
765: static PetscErrorCode PhysicsRiemann_Acoustics_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
766: {
767: AcousticsCtx *phys = (AcousticsCtx*)vctx;
768: PetscReal c = phys->c,z = phys->z;
769: PetscReal
770: Al[2][2] = {{-c/2 , c*z/2 },
771: {c/(2*z) , -c/2 }}, /* Left traveling waves */
772: Ar[2][2] = {{c/2 , c*z/2 },
773: {c/(2*z) , c/2 }}; /* Right traveling waves */
776: flux[0] = Al[0][0]*uR[0] + Al[0][1]*uR[1] + Ar[0][0]*uL[0] + Ar[0][1]*uL[1];
777: flux[1] = Al[1][0]*uR[0] + Al[1][1]*uR[1] + Ar[1][0]*uL[0] + Ar[1][1]*uL[1];
778: *maxspeed = c;
779: return(0);
780: }
784: static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx)
785: {
786: PetscErrorCode ierr;
787: AcousticsCtx *user;
788: PetscFunctionList rlist = 0,rclist = 0;
789: char rname[256] = "exact",rcname[256] = "characteristic";
792: PetscNew(*user,&user);
793: ctx->physics.sample = PhysicsSample_Acoustics;
794: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
795: ctx->physics.user = user;
796: ctx->physics.dof = 2;
798: PetscStrallocpy("u",&ctx->physics.fieldname[0]);
799: PetscStrallocpy("v",&ctx->physics.fieldname[1]);
801: user->c = 1;
802: user->z = 1;
804: RiemannListAdd(&rlist,"exact", PhysicsRiemann_Acoustics_Exact);
805: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_Acoustics);
806: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
807: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for linear Acoustics","");
808: {
809: PetscOptionsReal("-physics_acoustics_c","c = sqrt(bulk/rho)","",user->c,&user->c,NULL);
810: PetscOptionsReal("-physics_acoustics_z","z = sqrt(bulk*rho)","",user->z,&user->z,NULL);
811: PetscOptionsList("-physics_acoustics_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
812: PetscOptionsList("-physics_acoustics_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
813: }
814: PetscOptionsEnd();
815: RiemannListFind(rlist,rname,&ctx->physics.riemann);
816: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
817: PetscFunctionListDestroy(&rlist);
818: PetscFunctionListDestroy(&rclist);
819: return(0);
820: }
822: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */
824: typedef struct {
825: PetscReal acoustic_speed;
826: } IsoGasCtx;
828: PETSC_STATIC_INLINE void IsoGasFlux(PetscReal c,const PetscScalar *u,PetscScalar *f)
829: {
830: f[0] = u[1];
831: f[1] = PetscSqr(u[1])/u[0] + c*c*u[0];
832: }
836: static PetscErrorCode PhysicsSample_IsoGas(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u)
837: {
840: if (t > 0) SETERRQ(PETSC_COMM_SELF,1,"Exact solutions not implemented for t > 0");
841: switch (initial) {
842: case 0:
843: u[0] = (x < 0) ? 1 : 0.5;
844: u[1] = (x < 0) ? 1 : 0.7;
845: break;
846: case 1:
847: u[0] = 1+0.5*sin(2*M_PI*x);
848: u[1] = 1*u[0];
849: break;
850: default: SETERRQ(PETSC_COMM_SELF,1,"unknown initial condition");
851: }
852: return(0);
853: }
857: static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
858: {
859: IsoGasCtx *phys = (IsoGasCtx*)vctx;
860: PetscReal c = phys->acoustic_speed;
861: PetscScalar ubar,du[2],a[2],fL[2],fR[2],lam[2],ustar[2],R[2][2];
862: PetscInt i;
865: ubar = (uL[1]/PetscSqrtScalar(uL[0]) + uR[1]/PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0]));
866: /* write fluxuations in characteristic basis */
867: du[0] = uR[0] - uL[0];
868: du[1] = uR[1] - uL[1];
869: a[0] = (1/(2*c)) * ((ubar + c)*du[0] - du[1]);
870: a[1] = (1/(2*c)) * ((-ubar + c)*du[0] + du[1]);
871: /* wave speeds */
872: lam[0] = ubar - c;
873: lam[1] = ubar + c;
874: /* Right eigenvectors */
875: R[0][0] = 1; R[0][1] = ubar-c;
876: R[1][0] = 1; R[1][1] = ubar+c;
877: /* Compute state in star region (between the 1-wave and 2-wave) */
878: for (i=0; i<2; i++) ustar[i] = uL[i] + a[0]*R[0][i];
879: if (uL[1]/uL[0] < c && c < ustar[1]/ustar[0]) { /* 1-wave is sonic rarefaction */
880: PetscScalar ufan[2];
881: ufan[0] = uL[0]*PetscExpScalar(uL[1]/(uL[0]*c) - 1);
882: ufan[1] = c*ufan[0];
883: IsoGasFlux(c,ufan,flux);
884: } else if (ustar[1]/ustar[0] < -c && -c < uR[1]/uR[0]) { /* 2-wave is sonic rarefaction */
885: PetscScalar ufan[2];
886: ufan[0] = uR[0]*PetscExpScalar(-uR[1]/(uR[0]*c) - 1);
887: ufan[1] = -c*ufan[0];
888: IsoGasFlux(c,ufan,flux);
889: } else { /* Centered form */
890: IsoGasFlux(c,uL,fL);
891: IsoGasFlux(c,uR,fR);
892: for (i=0; i<2; i++) {
893: PetscScalar absdu = PetscAbsScalar(lam[0])*a[0]*R[0][i] + PetscAbsScalar(lam[1])*a[1]*R[1][i];
894: flux[i] = 0.5*(fL[i]+fR[i]) - 0.5*absdu;
895: }
896: }
897: *maxspeed = MaxAbs(lam[0],lam[1]);
898: return(0);
899: }
903: static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
904: {
905: IsoGasCtx *phys = (IsoGasCtx*)vctx;
906: PetscReal c = phys->acoustic_speed;
907: PetscScalar ustar[2];
908: struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
909: PetscInt i;
912: if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
913: {
914: /* Solve for star state */
915: PetscScalar res,tmp,rho = 0.5*(L.rho + R.rho); /* initial guess */
916: for (i=0; i<20; i++) {
917: PetscScalar fr,fl,dfr,dfl;
918: fl = (L.rho < rho)
919: ? (rho-L.rho)/PetscSqrtScalar(L.rho*rho) /* shock */
920: : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
921: fr = (R.rho < rho)
922: ? (rho-R.rho)/PetscSqrtScalar(R.rho*rho) /* shock */
923: : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
924: res = R.u-L.u + c*(fr+fl);
925: if (!isfinite(res)) SETERRQ1(PETSC_COMM_SELF,1,"non-finite residual=%g",res);
926: if (PetscAbsScalar(res) < 1e-10) {
927: star.rho = rho;
928: star.u = L.u - c*fl;
929: goto converged;
930: }
931: dfl = (L.rho < rho)
932: ? 1/PetscSqrtScalar(L.rho*rho)*(1 - 0.5*(rho-L.rho)/rho)
933: : 1/rho;
934: dfr = (R.rho < rho)
935: ? 1/PetscSqrtScalar(R.rho*rho)*(1 - 0.5*(rho-R.rho)/rho)
936: : 1/rho;
937: tmp = rho - res/(c*(dfr+dfl));
938: if (tmp <= 0) rho /= 2; /* Guard against Newton shooting off to a negative density */
939: else rho = tmp;
940: if (!((rho > 0) && isnormal(rho))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate rho=%g",rho);
941: }
942: SETERRQ1(PETSC_COMM_SELF,1,"Newton iteration for star.rho diverged after %d iterations",i);
943: }
944: converged:
945: if (L.u-c < 0 && 0 < star.u-c) { /* 1-wave is sonic rarefaction */
946: PetscScalar ufan[2];
947: ufan[0] = L.rho*PetscExpScalar(L.u/c - 1);
948: ufan[1] = c*ufan[0];
949: IsoGasFlux(c,ufan,flux);
950: } else if (star.u+c < 0 && 0 < R.u+c) { /* 2-wave is sonic rarefaction */
951: PetscScalar ufan[2];
952: ufan[0] = R.rho*PetscExpScalar(-R.u/c - 1);
953: ufan[1] = -c*ufan[0];
954: IsoGasFlux(c,ufan,flux);
955: } else if ((L.rho >= star.rho && L.u-c >= 0)
956: || (L.rho < star.rho && (star.rho*star.u-L.rho*L.u)/(star.rho-L.rho) > 0)) {
957: /* 1-wave is supersonic rarefaction, or supersonic shock */
958: IsoGasFlux(c,uL,flux);
959: } else if ((star.rho <= R.rho && R.u+c <= 0)
960: || (star.rho > R.rho && (R.rho*R.u-star.rho*star.u)/(R.rho-star.rho) < 0)) {
961: /* 2-wave is supersonic rarefaction or supersonic shock */
962: IsoGasFlux(c,uR,flux);
963: } else {
964: ustar[0] = star.rho;
965: ustar[1] = star.rho*star.u;
966: IsoGasFlux(c,ustar,flux);
967: }
968: *maxspeed = MaxAbs(MaxAbs(star.u-c,star.u+c),MaxAbs(L.u-c,R.u+c));
969: return(0);
970: }
974: static PetscErrorCode PhysicsRiemann_IsoGas_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
975: {
976: IsoGasCtx *phys = (IsoGasCtx*)vctx;
977: PetscScalar c = phys->acoustic_speed,fL[2],fR[2],s;
978: struct {PetscScalar rho,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]};
981: if (!(L.rho > 0 && R.rho > 0)) SETERRQ(PETSC_COMM_SELF,1,"Reconstructed density is negative");
982: IsoGasFlux(c,uL,fL);
983: IsoGasFlux(c,uR,fR);
984: s = PetscMax(PetscAbs(L.u),PetscAbs(R.u))+c;
985: flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(uL[0] - uR[0]);
986: flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]);
987: *maxspeed = s;
988: return(0);
989: }
993: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds)
994: {
995: IsoGasCtx *phys = (IsoGasCtx*)vctx;
996: PetscReal c = phys->acoustic_speed;
1000: speeds[0] = u[1]/u[0] - c;
1001: speeds[1] = u[1]/u[0] + c;
1002: X[0*2+0] = 1;
1003: X[0*2+1] = speeds[0];
1004: X[1*2+0] = 1;
1005: X[1*2+1] = speeds[1];
1007: PetscMemcpy(Xi,X,4*sizeof(X[0]));
1008: PetscKernel_A_gets_inverse_A_2(Xi,0);
1009: return(0);
1010: }
1014: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
1015: {
1016: PetscErrorCode ierr;
1017: IsoGasCtx *user;
1018: PetscFunctionList rlist = 0,rclist = 0;
1019: char rname[256] = "exact",rcname[256] = "characteristic";
1022: PetscNew(*user,&user);
1023: ctx->physics.sample = PhysicsSample_IsoGas;
1024: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
1025: ctx->physics.user = user;
1026: ctx->physics.dof = 2;
1028: PetscStrallocpy("density",&ctx->physics.fieldname[0]);
1029: PetscStrallocpy("momentum",&ctx->physics.fieldname[1]);
1031: user->acoustic_speed = 1;
1033: RiemannListAdd(&rlist,"exact", PhysicsRiemann_IsoGas_Exact);
1034: RiemannListAdd(&rlist,"roe", PhysicsRiemann_IsoGas_Roe);
1035: RiemannListAdd(&rlist,"rusanov",PhysicsRiemann_IsoGas_Rusanov);
1036: ReconstructListAdd(&rclist,"characteristic",PhysicsCharacteristic_IsoGas);
1037: ReconstructListAdd(&rclist,"conservative",PhysicsCharacteristic_Conservative);
1038: PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for IsoGas","");
1039: {
1040: PetscOptionsReal("-physics_isogas_acoustic_speed","Acoustic speed","",user->acoustic_speed,&user->acoustic_speed,NULL);
1041: PetscOptionsList("-physics_isogas_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
1042: PetscOptionsList("-physics_isogas_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL);
1043: }
1044: PetscOptionsEnd();
1045: RiemannListFind(rlist,rname,&ctx->physics.riemann);
1046: ReconstructListFind(rclist,rcname,&ctx->physics.characteristic);
1047: PetscFunctionListDestroy(&rlist);
1048: PetscFunctionListDestroy(&rclist);
1049: return(0);
1050: }
1054: /* --------------------------------- Shallow Water ----------------------------------- */
1056: typedef struct {
1057: PetscReal gravity;
1058: } ShallowCtx;
1060: PETSC_STATIC_INLINE void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f)
1061: {
1062: f[0] = u[1];
1063: f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]);
1064: }
1068: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed)
1069: {
1070: ShallowCtx *phys = (ShallowCtx*)vctx;
1071: PetscScalar g = phys->gravity,ustar[2],cL,cR,c,cstar;
1072: struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star;
1073: PetscInt i;
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: if (!isfinite(res)) SETERRQ1(PETSC_COMM_SELF,1,"non-finite residual=%g",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) && isnormal(h))) SETERRQ1(PETSC_COMM_SELF,1,"non-normal iterate h=%g",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,&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: PetscOptionsList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL);
1219: PetscOptionsList("-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;
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",tnow,dt,0.5/ctx->cfl_idt);
1333: } else SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",dt,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,MatStructure *flg,void *vctx)
1360: {
1361: FVCtx *ctx = (FVCtx*)vctx;
1363: PetscInt i,j,dof = ctx->physics.dof;
1364: PetscScalar *x,*J;
1365: PetscReal hx;
1366: DM da;
1367: DMDALocalInfo dainfo;
1370: TSGetDM(ts,&da);
1371: DMDAVecGetArray(da,X,&x);
1372: DMDAGetLocalInfo(da,&dainfo);
1373: hx = (ctx->xmax - ctx->xmin)/dainfo.mx;
1374: PetscMalloc(dof*dof*sizeof(PetscScalar),&J);
1375: for (i=dainfo.xs; i<dainfo.xs+dainfo.xm; i++) {
1376: (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds);
1377: for (j=0; j<dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]);
1378: SmallMatMultADB(J,dof,ctx->R,ctx->speeds,ctx->Rinv);
1379: for (j=0; j<dof*dof; j++) J[j] = J[j]/hx + shift*(j/dof == j%dof);
1380: MatSetValuesBlocked(*B,1,&i,1,&i,J,INSERT_VALUES);
1381: }
1382: PetscFree(J);
1383: DMDAVecRestoreArray(da,X,&x);
1385: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1386: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1387: if (*A != *B) {
1388: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1389: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1390: }
1391: return(0);
1392: }
1396: static PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
1397: {
1399: PetscScalar *u,*uj;
1400: PetscInt i,j,k,dof,xs,xm,Mx;
1403: if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,1,"Physics has not provided a sampling function");
1404: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1405: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1406: DMDAVecGetArray(da,U,&u);
1407: PetscMalloc(dof*sizeof(uj[0]),&uj);
1408: for (i=xs; i<xs+xm; i++) {
1409: const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h;
1410: const PetscInt N = 200;
1411: /* Integrate over cell i using trapezoid rule with N points. */
1412: for (k=0; k<dof; k++) u[i*dof+k] = 0;
1413: for (j=0; j<N+1; j++) {
1414: PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N;
1415: (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);
1416: for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
1417: }
1418: }
1419: DMDAVecRestoreArray(da,U,&u);
1420: PetscFree(uj);
1421: return(0);
1422: }
1426: static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
1427: {
1429: PetscReal xmin,xmax;
1430: PetscScalar sum,*x,tvsum,tvgsum;
1431: PetscInt imin,imax,Mx,i,j,xs,xm,dof;
1432: Vec Xloc;
1433: PetscBool iascii;
1436: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1437: if (iascii) {
1438: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
1439: DMGetLocalVector(da,&Xloc);
1440: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
1441: DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);
1442: DMDAVecGetArray(da,Xloc,&x);
1443: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1444: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1445: tvsum = 0;
1446: for (i=xs; i<xs+xm; i++) {
1447: for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j] - x[(i-1)*dof+j]);
1448: }
1449: MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));
1450: DMDAVecRestoreArray(da,Xloc,&x);
1451: DMRestoreLocalVector(da,&Xloc);
1453: VecMin(X,&imin,&xmin);
1454: VecMax(X,&imax,&xmax);
1455: VecSum(X,&sum);
1456: PetscViewerASCIIPrintf(viewer,"Solution range [%8.5f,%8.5f] with extrema at %d and %d, mean %8.5f, ||x||_TV %8.5f\n",xmin,xmax,imin,imax,sum/Mx,tvgsum/Mx);
1457: } else SETERRQ(PETSC_COMM_SELF,1,"Viewer type not supported");
1458: return(0);
1459: }
1463: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1,PetscReal *nrmsup)
1464: {
1466: Vec Y;
1467: PetscInt Mx;
1470: VecGetSize(X,&Mx);
1471: VecDuplicate(X,&Y);
1472: FVSample(ctx,da,t,Y);
1473: VecAYPX(Y,-1,X);
1474: VecNorm(Y,NORM_1,nrm1);
1475: VecNorm(Y,NORM_INFINITY,nrmsup);
1476: *nrm1 /= Mx;
1477: VecDestroy(&Y);
1478: return(0);
1479: }
1483: int main(int argc,char *argv[])
1484: {
1485: char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m";
1486: PetscFunctionList limiters = 0,physics = 0;
1487: MPI_Comm comm;
1488: TS ts;
1489: DM da;
1490: Vec X,X0,R;
1491: Mat B;
1492: FVCtx ctx;
1493: PetscInt i,dof,xs,xm,Mx,draw = 0;
1494: PetscBool view_final = PETSC_FALSE;
1495: PetscReal ptime;
1496: PetscErrorCode ierr;
1498: PetscInitialize(&argc,&argv,0,help);
1499: comm = PETSC_COMM_WORLD;
1500: PetscMemzero(&ctx,sizeof(ctx));
1502: /* Register limiters to be available on the command line */
1503: PetscFunctionListAdd(&limiters,"upwind" ,Limit_Upwind);
1504: PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit_LaxWendroff);
1505: PetscFunctionListAdd(&limiters,"beam-warming" ,Limit_BeamWarming);
1506: PetscFunctionListAdd(&limiters,"fromm" ,Limit_Fromm);
1507: PetscFunctionListAdd(&limiters,"minmod" ,Limit_Minmod);
1508: PetscFunctionListAdd(&limiters,"superbee" ,Limit_Superbee);
1509: PetscFunctionListAdd(&limiters,"mc" ,Limit_MC);
1510: PetscFunctionListAdd(&limiters,"vanleer" ,Limit_VanLeer);
1511: PetscFunctionListAdd(&limiters,"vanalbada" ,Limit_VanAlbada);
1512: PetscFunctionListAdd(&limiters,"vanalbadatvd" ,Limit_VanAlbadaTVD);
1513: PetscFunctionListAdd(&limiters,"koren" ,Limit_Koren);
1514: PetscFunctionListAdd(&limiters,"korensym" ,Limit_KorenSym);
1515: PetscFunctionListAdd(&limiters,"koren3" ,Limit_Koren3);
1516: PetscFunctionListAdd(&limiters,"cada-torrilhon2" ,Limit_CadaTorrilhon2);
1517: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1);
1518: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1" ,Limit_CadaTorrilhon3R1);
1519: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10);
1520: PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100);
1522: /* Register physical models to be available on the command line */
1523: PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect);
1524: PetscFunctionListAdd(&physics,"burgers" ,PhysicsCreate_Burgers);
1525: PetscFunctionListAdd(&physics,"traffic" ,PhysicsCreate_Traffic);
1526: PetscFunctionListAdd(&physics,"acoustics" ,PhysicsCreate_Acoustics);
1527: PetscFunctionListAdd(&physics,"isogas" ,PhysicsCreate_IsoGas);
1528: PetscFunctionListAdd(&physics,"shallow" ,PhysicsCreate_Shallow);
1530: ctx.comm = comm;
1531: ctx.cfl = 0.9; ctx.bctype = FVBC_PERIODIC;
1532: ctx.xmin = -1; ctx.xmax = 1;
1533: PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");
1534: {
1535: PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);
1536: PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);
1537: PetscOptionsList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL);
1538: PetscOptionsList("-physics","Name of physics (Riemann solver and characteristics) to use","",physics,physname,physname,sizeof(physname),NULL);
1539: PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);
1540: PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final);
1541: PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);
1542: PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);
1543: PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);
1544: PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);
1545: }
1546: PetscOptionsEnd();
1548: /* Choose the limiter from the list of registered limiters */
1549: PetscFunctionListFind(limiters,lname,&ctx.limit);
1550: if (!ctx.limit) SETERRQ1(PETSC_COMM_SELF,1,"Limiter '%s' not found",lname);
1552: /* Choose the physics from the list of registered models */
1553: {
1554: PetscErrorCode (*r)(FVCtx*);
1555: PetscFunctionListFind(physics,physname,&r);
1556: if (!r) SETERRQ1(PETSC_COMM_SELF,1,"Physics '%s' not found",physname);
1557: /* Create the physics, will set the number of fields and their names */
1558: (*r)(&ctx);
1559: }
1561: /* Create a DMDA to manage the parallel grid */
1562: DMDACreate1d(comm,DMDA_BOUNDARY_PERIODIC,-50,ctx.physics.dof,2,NULL,&da);
1563: /* Inform the DMDA of the field names provided by the physics. */
1564: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1565: for (i=0; i<ctx.physics.dof; i++) {
1566: DMDASetFieldName(da,i,ctx.physics.fieldname[i]);
1567: }
1568: DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);
1569: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
1571: /* Set coordinates of cell centers */
1572: DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0);
1574: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1575: PetscMalloc4(dof*dof,PetscScalar,&ctx.R,dof*dof,PetscScalar,&ctx.Rinv,2*dof,PetscScalar,&ctx.cjmpLR,1*dof,PetscScalar,&ctx.cslope);
1576: PetscMalloc3(2*dof,PetscScalar,&ctx.uLR,dof,PetscScalar,&ctx.flux,dof,PetscReal,&ctx.speeds);
1578: /* Create a vector to store the solution and to save the initial state */
1579: DMCreateGlobalVector(da,&X);
1580: VecDuplicate(X,&X0);
1581: VecDuplicate(X,&R);
1583: DMCreateMatrix(da,NULL,&B);
1585: /* Create a time-stepping object */
1586: TSCreate(comm,&ts);
1587: TSSetDM(ts,da);
1588: TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);
1589: TSSetIJacobian(ts,B,B,FVIJacobian,&ctx);
1590: TSSetType(ts,TSSSP);
1591: TSSetDuration(ts,1000,10);
1593: /* Compute initial conditions and starting time step */
1594: FVSample(&ctx,da,0,X0);
1595: FVRHSFunction(ts,0,X0,X,(void*)&ctx); /* Initial function evaluation, only used to determine max speed */
1596: VecCopy(X0,X); /* The function value was not used so we set X=X0 again */
1597: TSSetInitialTimeStep(ts,0,ctx.cfl/ctx.cfl_idt);
1599: TSSetFromOptions(ts); /* Take runtime options */
1601: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1603: {
1604: PetscReal nrm1,nrmsup;
1605: PetscInt steps;
1607: TSSolve(ts,X);
1608: TSGetSolveTime(ts,&ptime);
1609: TSGetTimeStepNumber(ts,&steps);
1611: PetscPrintf(comm,"Final time %8.5f, steps %d\n",ptime,steps);
1612: if (ctx.exact) {
1613: SolutionErrorNorms(&ctx,da,ptime,X,&nrm1,&nrmsup);
1614: PetscPrintf(comm,"Error ||x-x_e||_1 %8.4e ||x-x_e||_sup %8.4e\n",nrm1,nrmsup);
1615: }
1616: }
1618: SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);
1620: if (draw & 0x1) {VecView(X0,PETSC_VIEWER_DRAW_WORLD);}
1621: if (draw & 0x2) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
1622: if (draw & 0x4) {
1623: Vec Y;
1624: VecDuplicate(X,&Y);
1625: FVSample(&ctx,da,ptime,Y);
1626: VecAYPX(Y,-1,X);
1627: VecView(Y,PETSC_VIEWER_DRAW_WORLD);
1628: VecDestroy(&Y);
1629: }
1631: if (view_final) {
1632: PetscViewer viewer;
1633: PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);
1634: PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1635: VecView(X,viewer);
1636: PetscViewerDestroy(&viewer);
1637: }
1639: /* Clean up */
1640: (*ctx.physics.destroy)(ctx.physics.user);
1641: for (i=0; i<ctx.physics.dof; i++) {PetscFree(ctx.physics.fieldname[i]);}
1642: PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope);
1643: PetscFree3(ctx.uLR,ctx.flux,ctx.speeds);
1644: VecDestroy(&X);
1645: VecDestroy(&X0);
1646: VecDestroy(&R);
1647: MatDestroy(&B);
1648: DMDestroy(&da);
1649: TSDestroy(&ts);
1650: PetscFunctionListDestroy(&limiters);
1651: PetscFunctionListDestroy(&physics);
1652: PetscFinalize();
1653: return 0;
1654: }