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