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