1 #include <petscdm.h>
2 #include <petscdmceed.h>
3
4 #ifdef __CUDACC_RTC__
5 #define PETSC_HAVE_LIBCEED
6 // Define PETSc types to be equal to Ceed types
7 typedef CeedInt PetscInt;
8 typedef CeedScalar PetscReal;
9 typedef CeedScalar PetscScalar;
10 typedef CeedInt PetscErrorCode;
11 // Define things we are missing from PETSc headers
12 #undef PETSC_SUCCESS
13 #define PETSC_SUCCESS 0
14 #define PETSC_COMM_SELF MPI_COMM_SELF
15 #undef PetscFunctionBeginUser
16 #define PetscFunctionBeginUser
17 #undef PetscFunctionReturn
18 #define PetscFunctionReturn(x) return x
19 #undef PetscCall
20 #define PetscCall(a) a
21 #define PetscFunctionReturnVoid() return
22 // Math definitions
23 #undef PetscSqrtReal
24 #define PetscSqrtReal(x) sqrt(x)
25 #undef PetscSqrtScalar
26 #define PetscSqrtScalar(x) sqrt(x)
27 #undef PetscSqr
28 #define PetscSqr(x) (x * x)
29 #define PetscSqrReal(x) (x * x)
30 #define PetscAbsReal(x) abs(x)
31 #define PetscAbsScalar(x) abs(x)
32 #define PetscMax(x, y) x > y ? x : y
33 #define PetscMin(x, y) x < y ? x : y
34 #define PetscRealPart(a) a
35 #define PetscPowScalar(a, b) pow(a, b)
36 #endif
37
38 #define DIM 2 /* Geometric dimension */
39
40 /* Represents continuum physical equations. */
41 typedef struct _n_Physics *Physics;
42
43 /* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
44 * discretization-independent, but its members depend on the scenario being solved. */
45 typedef struct _n_Model *Model;
46
47 struct FieldDescription {
48 const char *name;
49 PetscInt dof;
50 };
51
52 struct _n_Physics {
53 void (*riemann)(PetscInt, PetscInt, const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscInt, const PetscScalar[], PetscScalar[], void *);
54 PetscInt dof; /* number of degrees of freedom per cell */
55 PetscReal maxspeed; /* kludge to pick initial time step, need to add monitoring and step control */
56 void *data;
57 PetscInt nfields;
58 const struct FieldDescription *field_desc;
59 };
60
61 typedef struct {
62 PetscReal gravity;
63 struct {
64 PetscInt Height;
65 PetscInt Speed;
66 PetscInt Energy;
67 } functional;
68 } Physics_SW;
69
70 typedef struct {
71 PetscReal h;
72 PetscReal uh[DIM];
73 } SWNode;
74 typedef union
75 {
76 SWNode swnode;
77 PetscReal vals[DIM + 1];
78 } SWNodeUnion;
79
80 typedef enum {
81 EULER_IV_SHOCK,
82 EULER_SS_SHOCK,
83 EULER_SHOCK_TUBE,
84 EULER_LINEAR_WAVE
85 } EulerType;
86
87 typedef struct {
88 PetscReal gamma;
89 PetscReal rhoR;
90 PetscReal amach;
91 PetscReal itana;
92 EulerType type;
93 struct {
94 PetscInt Density;
95 PetscInt Momentum;
96 PetscInt Energy;
97 PetscInt Pressure;
98 PetscInt Speed;
99 } monitor;
100 } Physics_Euler;
101
102 typedef struct {
103 PetscReal r;
104 PetscReal ru[DIM];
105 PetscReal E;
106 } EulerNode;
107 typedef union
108 {
109 EulerNode eulernode;
110 PetscReal vals[DIM + 2];
111 } EulerNodeUnion;
112
Dot2Real(const PetscReal * x,const PetscReal * y)113 static inline PetscReal Dot2Real(const PetscReal *x, const PetscReal *y)
114 {
115 return x[0] * y[0] + x[1] * y[1];
116 }
Norm2Real(const PetscReal * x)117 static inline PetscReal Norm2Real(const PetscReal *x)
118 {
119 return PetscSqrtReal(PetscAbsReal(Dot2Real(x, x)));
120 }
Normalize2Real(PetscReal * x)121 static inline void Normalize2Real(PetscReal *x)
122 {
123 PetscReal a = 1. / Norm2Real(x);
124 x[0] *= a;
125 x[1] *= a;
126 }
Scale2Real(PetscReal a,const PetscReal * x,PetscReal * y)127 static inline void Scale2Real(PetscReal a, const PetscReal *x, PetscReal *y)
128 {
129 y[0] = a * x[0];
130 y[1] = a * x[1];
131 }
132
DotDIMReal(const PetscReal * x,const PetscReal * y)133 static inline PetscReal DotDIMReal(const PetscReal *x, const PetscReal *y)
134 {
135 PetscInt i;
136 PetscReal prod = 0.0;
137
138 for (i = 0; i < DIM; i++) prod += x[i] * y[i];
139 return prod;
140 }
NormDIM(const PetscReal * x)141 static inline PetscReal NormDIM(const PetscReal *x)
142 {
143 return PetscSqrtReal(PetscAbsReal(DotDIMReal(x, x)));
144 }
Waxpy2Real(PetscReal a,const PetscReal * x,const PetscReal * y,PetscReal * w)145 static inline void Waxpy2Real(PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
146 {
147 w[0] = a * x[0] + y[0];
148 w[1] = a * x[1] + y[1];
149 }
150
151 /*
152 * h_t + div(uh) = 0
153 * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
154 *
155 * */
SWFlux(Physics phys,const PetscReal * n,const SWNode * x,SWNode * f)156 static PetscErrorCode SWFlux(Physics phys, const PetscReal *n, const SWNode *x, SWNode *f)
157 {
158 Physics_SW *sw = (Physics_SW *)phys->data;
159 PetscReal uhn, u[DIM];
160 PetscInt i;
161
162 PetscFunctionBeginUser;
163 Scale2Real(1. / x->h, x->uh, u);
164 uhn = x->uh[0] * n[0] + x->uh[1] * n[1];
165 f->h = uhn;
166 for (i = 0; i < DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
167 PetscFunctionReturn(PETSC_SUCCESS);
168 }
169
PhysicsRiemann_SW_Rusanov(PetscInt dim,PetscInt Nf,const PetscReal * qp,const PetscReal * n,const PetscScalar * xL,const PetscScalar * xR,PetscInt numConstants,const PetscScalar constants[],PetscScalar * flux,Physics phys)170 static void PhysicsRiemann_SW_Rusanov(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys)
171 {
172 Physics_SW *sw = (Physics_SW *)phys->data;
173 PetscReal cL, cR, speed;
174 PetscReal nn[DIM];
175 #if !defined(PETSC_USE_COMPLEX)
176 const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
177 #else
178 SWNodeUnion uLreal, uRreal;
179 const SWNode *uL = &uLreal.swnode;
180 const SWNode *uR = &uRreal.swnode;
181 #endif
182 SWNodeUnion fL, fR;
183 PetscInt i;
184 PetscReal zero = 0.;
185 PetscErrorCode ierr;
186
187 #if defined(PETSC_USE_COMPLEX)
188 uLreal.swnode.h = 0;
189 uRreal.swnode.h = 0;
190 for (i = 0; i < 1 + dim; i++) uLreal.vals[i] = PetscRealPart(xL[i]);
191 for (i = 0; i < 1 + dim; i++) uRreal.vals[i] = PetscRealPart(xR[i]);
192 #endif
193
194 if (uL->h < 0 || uR->h < 0) {
195 // reconstructed thickness is negative
196 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
197 for (i = 0; i < 1 + dim; ++i) flux[i] = zero / zero;
198 PetscCallVoid(PetscFPTrapPop());
199 return;
200 }
201
202 nn[0] = n[0];
203 nn[1] = n[1];
204 Normalize2Real(nn);
205 ierr = SWFlux(phys, nn, uL, &fL.swnode);
206 if (ierr) {
207 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
208 for (i = 0; i < 1 + dim; ++i) fL.vals[i] = zero / zero;
209 PetscCallVoid(PetscFPTrapPop());
210 }
211 ierr = SWFlux(phys, nn, uR, &fR.swnode);
212 if (ierr) {
213 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
214 for (i = 0; i < 1 + dim; ++i) fR.vals[i] = zero / zero;
215 PetscCallVoid(PetscFPTrapPop());
216 }
217 cL = PetscSqrtReal(sw->gravity * uL->h);
218 cR = PetscSqrtReal(sw->gravity * uR->h); /* gravity wave speed */
219 speed = PetscMax(PetscAbsReal(Dot2Real(uL->uh, nn) / uL->h) + cL, PetscAbsReal(Dot2Real(uR->uh, nn) / uR->h) + cR);
220 for (i = 0; i < 1 + dim; i++) flux[i] = (0.5 * (fL.vals[i] + fR.vals[i]) + 0.5 * speed * (xL[i] - xR[i])) * Norm2Real(n);
221 #if 0
222 PetscPrintf(PETSC_COMM_SELF, "Rusanov Flux (%g)\n", sw->gravity);
223 for (PetscInt j = 0; j < 3; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g |\n", flux[j]);
224 #endif
225 }
226
227 #ifdef PETSC_HAVE_LIBCEED
PhysicsRiemann_SW_Rusanov_CEED(PetscCtx ctx,CeedInt Q,const CeedScalar * const in[],CeedScalar * const out[])228 CEED_QFUNCTION(PhysicsRiemann_SW_Rusanov_CEED)(PetscCtx ctx, CeedInt Q, const CeedScalar *const in[], CeedScalar *const out[])
229 {
230 const CeedScalar *xL = in[0], *xR = in[1], *geom = in[2];
231 CeedScalar *cL = out[0], *cR = out[1];
232 const Physics_SW *sw = (Physics_SW *)ctx;
233 struct _n_Physics phys;
234 #if 0
235 const CeedScalar *info = in[3];
236 #endif
237
238 phys.data = (void *)sw;
239 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
240 {
241 const CeedScalar qL[3] = {xL[i + Q * 0], xL[i + Q * 1], xL[i + Q * 2]};
242 const CeedScalar qR[3] = {xR[i + Q * 0], xR[i + Q * 1], xR[i + Q * 2]};
243 const CeedScalar n[2] = {geom[i + Q * 0], geom[i + Q * 1]};
244 CeedScalar flux[3];
245
246 #if 0
247 PetscPrintf(PETSC_COMM_SELF, "Face %d Normal\n", (int)info[i + Q * 0]);
248 for (CeedInt j = 0; j < DIM; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g |\n", n[j]);
249 PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: left state\n", (int)info[i + Q * 1]);
250 for (CeedInt j = 0; j < DIM + 1; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g |\n", qL[j]);
251 PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: right state\n", (int)info[i + Q * 2]);
252 for (CeedInt j = 0; j < DIM + 1; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g |\n", qR[j]);
253 #endif
254 PhysicsRiemann_SW_Rusanov(DIM, DIM + 1, NULL, n, qL, qR, 0, NULL, flux, &phys);
255 for (CeedInt j = 0; j < 3; ++j) {
256 cL[i + Q * j] = -flux[j] / geom[i + Q * 2];
257 cR[i + Q * j] = flux[j] / geom[i + Q * 3];
258 }
259 #if 0
260 PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: left flux\n", (int)info[i + Q * 1]);
261 for (CeedInt j = 0; j < DIM + 1; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g | (%g)\n", cL[i + Q * j], geom[i + Q * 2]);
262 PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: right flux\n", (int)info[i + Q * 2]);
263 for (CeedInt j = 0; j < DIM + 1; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g | (%g)\n", cR[i + Q * j], geom[i + Q * 3]);
264 #endif
265 }
266 return CEED_ERROR_SUCCESS;
267 }
268 #endif
269
PhysicsRiemann_SW_HLL(PetscInt dim,PetscInt Nf,const PetscReal * qp,const PetscReal * n,const PetscScalar * xL,const PetscScalar * xR,PetscInt numConstants,const PetscScalar constants[],PetscScalar * flux,Physics phys)270 static void PhysicsRiemann_SW_HLL(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys)
271 {
272 Physics_SW *sw = (Physics_SW *)phys->data;
273 PetscReal aL, aR;
274 PetscReal nn[DIM];
275 #if !defined(PETSC_USE_COMPLEX)
276 const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
277 #else
278 SWNodeUnion uLreal, uRreal;
279 const SWNode *uL = &uLreal.swnode;
280 const SWNode *uR = &uRreal.swnode;
281 #endif
282 SWNodeUnion fL, fR;
283 PetscInt i;
284 PetscReal zero = 0.;
285 PetscErrorCode ierr;
286
287 #if defined(PETSC_USE_COMPLEX)
288 uLreal.swnode.h = 0;
289 uRreal.swnode.h = 0;
290 for (i = 0; i < 1 + dim; i++) uLreal.vals[i] = PetscRealPart(xL[i]);
291 for (i = 0; i < 1 + dim; i++) uRreal.vals[i] = PetscRealPart(xR[i]);
292 #endif
293 if (uL->h <= 0 || uR->h <= 0) {
294 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
295 for (i = 0; i < 1 + dim; i++) flux[i] = zero;
296 PetscCallVoid(PetscFPTrapPop());
297 return;
298 } /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative"); */
299 nn[0] = n[0];
300 nn[1] = n[1];
301 Normalize2Real(nn);
302 ierr = SWFlux(phys, nn, uL, &fL.swnode);
303 if (ierr) {
304 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
305 for (i = 0; i < 1 + dim; ++i) fL.vals[i] = zero / zero;
306 PetscCallVoid(PetscFPTrapPop());
307 }
308 ierr = SWFlux(phys, nn, uR, &fR.swnode);
309 if (ierr) {
310 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
311 for (i = 0; i < 1 + dim; ++i) fR.vals[i] = zero / zero;
312 PetscCallVoid(PetscFPTrapPop());
313 }
314 /* gravity wave speed */
315 aL = PetscSqrtReal(sw->gravity * uL->h);
316 aR = PetscSqrtReal(sw->gravity * uR->h);
317 // Defining u_tilda and v_tilda as u and v
318 PetscReal u_L, u_R;
319 u_L = Dot2Real(uL->uh, nn) / uL->h;
320 u_R = Dot2Real(uR->uh, nn) / uR->h;
321 PetscReal sL, sR;
322 sL = PetscMin(u_L - aL, u_R - aR);
323 sR = PetscMax(u_L + aL, u_R + aR);
324 if (sL > zero) {
325 for (i = 0; i < dim + 1; i++) flux[i] = fL.vals[i] * Norm2Real(n);
326 } else if (sR < zero) {
327 for (i = 0; i < dim + 1; i++) flux[i] = fR.vals[i] * Norm2Real(n);
328 } else {
329 for (i = 0; i < dim + 1; i++) flux[i] = ((sR * fL.vals[i] - sL * fR.vals[i] + sR * sL * (xR[i] - xL[i])) / (sR - sL)) * Norm2Real(n);
330 }
331 }
332
Pressure_PG(const PetscReal gamma,const EulerNode * x,PetscReal * p)333 static PetscErrorCode Pressure_PG(const PetscReal gamma, const EulerNode *x, PetscReal *p)
334 {
335 PetscReal ru2;
336
337 PetscFunctionBeginUser;
338 ru2 = DotDIMReal(x->ru, x->ru);
339 (*p) = (x->E - 0.5 * ru2 / x->r) * (gamma - 1.0); /* (E - rho V^2/2)(gamma-1) = e rho (gamma-1) */
340 PetscFunctionReturn(PETSC_SUCCESS);
341 }
342
SpeedOfSound_PG(const PetscReal gamma,const EulerNode * x,PetscReal * c)343 static PetscErrorCode SpeedOfSound_PG(const PetscReal gamma, const EulerNode *x, PetscReal *c)
344 {
345 PetscReal p;
346
347 PetscFunctionBeginUser;
348 PetscCall(Pressure_PG(gamma, x, &p));
349 /* gamma = heat capacity ratio */
350 (*c) = PetscSqrtReal(gamma * p / x->r);
351 PetscFunctionReturn(PETSC_SUCCESS);
352 }
353
354 /*
355 * x = (rho,rho*(u_1),...,rho*e)^T
356 * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
357 *
358 * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
359 *
360 */
EulerFlux(Physics phys,const PetscReal * n,const EulerNode * x,EulerNode * f)361 static PetscErrorCode EulerFlux(Physics phys, const PetscReal *n, const EulerNode *x, EulerNode *f)
362 {
363 Physics_Euler *eu = (Physics_Euler *)phys->data;
364 PetscReal nu, p;
365 PetscInt i;
366
367 PetscFunctionBeginUser;
368 PetscCall(Pressure_PG(eu->gamma, x, &p));
369 nu = DotDIMReal(x->ru, n);
370 f->r = nu; /* A rho u */
371 nu /= x->r; /* A u */
372 for (i = 0; i < DIM; i++) f->ru[i] = nu * x->ru[i] + n[i] * p; /* r u^2 + p */
373 f->E = nu * (x->E + p); /* u(e+p) */
374 PetscFunctionReturn(PETSC_SUCCESS);
375 }
376
377 /* Godunov fluxs */
cvmgp_(PetscScalar * a,PetscScalar * b,PetscScalar * test)378 static PetscScalar cvmgp_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
379 {
380 /* System generated locals */
381 PetscScalar ret_val;
382
383 if (PetscRealPart(*test) > 0.) goto L10;
384 ret_val = *b;
385 return ret_val;
386 L10:
387 ret_val = *a;
388 return ret_val;
389 } /* cvmgp_ */
390
cvmgm_(PetscScalar * a,PetscScalar * b,PetscScalar * test)391 static PetscScalar cvmgm_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
392 {
393 /* System generated locals */
394 PetscScalar ret_val;
395
396 if (PetscRealPart(*test) < 0.) goto L10;
397 ret_val = *b;
398 return ret_val;
399 L10:
400 ret_val = *a;
401 return ret_val;
402 } /* cvmgm_ */
403
riem1mdt(PetscScalar * gaml,PetscScalar * gamr,PetscScalar * rl,PetscScalar * pl,PetscScalar * uxl,PetscScalar * rr,PetscScalar * pr,PetscScalar * uxr,PetscScalar * rstarl,PetscScalar * rstarr,PetscScalar * pstar,PetscScalar * ustar)404 static int riem1mdt(PetscScalar *gaml, PetscScalar *gamr, PetscScalar *rl, PetscScalar *pl, PetscScalar *uxl, PetscScalar *rr, PetscScalar *pr, PetscScalar *uxr, PetscScalar *rstarl, PetscScalar *rstarr, PetscScalar *pstar, PetscScalar *ustar)
405 {
406 /* Initialized data */
407
408 static PetscScalar smallp = 1e-8;
409
410 /* System generated locals */
411 int i__1;
412 PetscScalar d__1, d__2;
413
414 /* Local variables */
415 static int i0;
416 static PetscScalar cl, cr, wl, zl, wr, zr, pst, durl, skpr1, skpr2;
417 static int iwave;
418 static PetscScalar gascl4, gascr4, cstarl, dpstar, cstarr;
419 /* static PetscScalar csqrl, csqrr, gascl1, gascl2, gascl3, gascr1, gascr2, gascr3; */
420 static int iterno;
421 static PetscScalar ustarl, ustarr, rarepr1, rarepr2;
422
423 /* gascl1 = *gaml - 1.; */
424 /* gascl2 = (*gaml + 1.) * .5; */
425 /* gascl3 = gascl2 / *gaml; */
426 gascl4 = 1. / (*gaml - 1.);
427
428 /* gascr1 = *gamr - 1.; */
429 /* gascr2 = (*gamr + 1.) * .5; */
430 /* gascr3 = gascr2 / *gamr; */
431 gascr4 = 1. / (*gamr - 1.);
432 iterno = 10;
433 /* find pstar: */
434 cl = PetscSqrtScalar(*gaml * *pl / *rl);
435 cr = PetscSqrtScalar(*gamr * *pr / *rr);
436 wl = *rl * cl;
437 wr = *rr * cr;
438 /* csqrl = wl * wl; */
439 /* csqrr = wr * wr; */
440 *pstar = (wl * *pr + wr * *pl) / (wl + wr);
441 *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
442 pst = *pl / *pr;
443 skpr1 = cr * (pst - 1.) * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
444 d__1 = (*gamr - 1.) / (*gamr * 2.);
445 rarepr2 = gascr4 * 2. * cr * (1. - PetscPowScalar(pst, d__1));
446 pst = *pr / *pl;
447 skpr2 = cl * (pst - 1.) * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
448 d__1 = (*gaml - 1.) / (*gaml * 2.);
449 rarepr1 = gascl4 * 2. * cl * (1. - PetscPowScalar(pst, d__1));
450 durl = *uxr - *uxl;
451 if (PetscRealPart(*pr) < PetscRealPart(*pl)) {
452 if (PetscRealPart(durl) >= PetscRealPart(rarepr1)) {
453 iwave = 100;
454 } else if (PetscRealPart(durl) <= PetscRealPart(-skpr1)) {
455 iwave = 300;
456 } else {
457 iwave = 400;
458 }
459 } else {
460 if (PetscRealPart(durl) >= PetscRealPart(rarepr2)) {
461 iwave = 100;
462 } else if (PetscRealPart(durl) <= PetscRealPart(-skpr2)) {
463 iwave = 300;
464 } else {
465 iwave = 200;
466 }
467 }
468 if (iwave == 100) {
469 /* 1-wave: rarefaction wave, 3-wave: rarefaction wave */
470 /* case (100) */
471 i__1 = iterno;
472 for (i0 = 1; i0 <= i__1; ++i0) {
473 d__1 = *pstar / *pl;
474 d__2 = 1. / *gaml;
475 *rstarl = *rl * PetscPowScalar(d__1, d__2);
476 cstarl = PetscSqrtScalar(*gaml * *pstar / *rstarl);
477 ustarl = *uxl - gascl4 * 2. * (cstarl - cl);
478 zl = *rstarl * cstarl;
479 d__1 = *pstar / *pr;
480 d__2 = 1. / *gamr;
481 *rstarr = *rr * PetscPowScalar(d__1, d__2);
482 cstarr = PetscSqrtScalar(*gamr * *pstar / *rstarr);
483 ustarr = *uxr + gascr4 * 2. * (cstarr - cr);
484 zr = *rstarr * cstarr;
485 dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
486 *pstar -= dpstar;
487 *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
488 if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
489 #if 0
490 break;
491 #endif
492 }
493 }
494 /* 1-wave: shock wave, 3-wave: rarefaction wave */
495 } else if (iwave == 200) {
496 /* case (200) */
497 i__1 = iterno;
498 for (i0 = 1; i0 <= i__1; ++i0) {
499 pst = *pstar / *pl;
500 ustarl = *uxl - (pst - 1.) * cl * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
501 zl = *pl / cl * PetscSqrtScalar(*gaml * 2. * (*gaml - 1. + (*gaml + 1.) * pst)) * (*gaml - 1. + (*gaml + 1.) * pst) / (*gaml * 3. - 1. + (*gaml + 1.) * pst);
502 d__1 = *pstar / *pr;
503 d__2 = 1. / *gamr;
504 *rstarr = *rr * PetscPowScalar(d__1, d__2);
505 cstarr = PetscSqrtScalar(*gamr * *pstar / *rstarr);
506 zr = *rstarr * cstarr;
507 ustarr = *uxr + gascr4 * 2. * (cstarr - cr);
508 dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
509 *pstar -= dpstar;
510 *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
511 if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
512 #if 0
513 break;
514 #endif
515 }
516 }
517 /* 1-wave: shock wave, 3-wave: shock */
518 } else if (iwave == 300) {
519 /* case (300) */
520 i__1 = iterno;
521 for (i0 = 1; i0 <= i__1; ++i0) {
522 pst = *pstar / *pl;
523 ustarl = *uxl - (pst - 1.) * cl * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
524 zl = *pl / cl * PetscSqrtScalar(*gaml * 2. * (*gaml - 1. + (*gaml + 1.) * pst)) * (*gaml - 1. + (*gaml + 1.) * pst) / (*gaml * 3. - 1. + (*gaml + 1.) * pst);
525 pst = *pstar / *pr;
526 ustarr = *uxr + (pst - 1.) * cr * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
527 zr = *pr / cr * PetscSqrtScalar(*gamr * 2. * (*gamr - 1. + (*gamr + 1.) * pst)) * (*gamr - 1. + (*gamr + 1.) * pst) / (*gamr * 3. - 1. + (*gamr + 1.) * pst);
528 dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
529 *pstar -= dpstar;
530 *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
531 if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
532 #if 0
533 break;
534 #endif
535 }
536 }
537 /* 1-wave: rarefaction wave, 3-wave: shock */
538 } else if (iwave == 400) {
539 /* case (400) */
540 i__1 = iterno;
541 for (i0 = 1; i0 <= i__1; ++i0) {
542 d__1 = *pstar / *pl;
543 d__2 = 1. / *gaml;
544 *rstarl = *rl * PetscPowScalar(d__1, d__2);
545 cstarl = PetscSqrtScalar(*gaml * *pstar / *rstarl);
546 ustarl = *uxl - gascl4 * 2. * (cstarl - cl);
547 zl = *rstarl * cstarl;
548 pst = *pstar / *pr;
549 ustarr = *uxr + (pst - 1.) * cr * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
550 zr = *pr / cr * PetscSqrtScalar(*gamr * 2. * (*gamr - 1. + (*gamr + 1.) * pst)) * (*gamr - 1. + (*gamr + 1.) * pst) / (*gamr * 3. - 1. + (*gamr + 1.) * pst);
551 dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
552 *pstar -= dpstar;
553 *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
554 if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
555 #if 0
556 break;
557 #endif
558 }
559 }
560 }
561
562 *ustar = (zl * ustarr + zr * ustarl) / (zl + zr);
563 if (PetscRealPart(*pstar) > PetscRealPart(*pl)) {
564 pst = *pstar / *pl;
565 *rstarl = ((*gaml + 1.) * pst + *gaml - 1.) / ((*gaml - 1.) * pst + *gaml + 1.) * *rl;
566 }
567 if (PetscRealPart(*pstar) > PetscRealPart(*pr)) {
568 pst = *pstar / *pr;
569 *rstarr = ((*gamr + 1.) * pst + *gamr - 1.) / ((*gamr - 1.) * pst + *gamr + 1.) * *rr;
570 }
571 return iwave;
572 }
573
sign(PetscScalar x)574 static PetscScalar sign(PetscScalar x)
575 {
576 if (PetscRealPart(x) > 0) return 1.0;
577 if (PetscRealPart(x) < 0) return -1.0;
578 return 0.0;
579 }
580 /* Riemann Solver */
581 /* -------------------------------------------------------------------- */
riemannsolver(PetscScalar * xcen,PetscScalar * xp,PetscScalar * dtt,PetscScalar * rl,PetscScalar * uxl,PetscScalar * pl,PetscScalar * utl,PetscScalar * ubl,PetscScalar * gaml,PetscScalar * rho1l,PetscScalar * rr,PetscScalar * uxr,PetscScalar * pr,PetscScalar * utr,PetscScalar * ubr,PetscScalar * gamr,PetscScalar * rho1r,PetscScalar * rx,PetscScalar * uxm,PetscScalar * px,PetscScalar * utx,PetscScalar * ubx,PetscScalar * gam,PetscScalar * rho1)582 static int riemannsolver(PetscScalar *xcen, PetscScalar *xp, PetscScalar *dtt, PetscScalar *rl, PetscScalar *uxl, PetscScalar *pl, PetscScalar *utl, PetscScalar *ubl, PetscScalar *gaml, PetscScalar *rho1l, PetscScalar *rr, PetscScalar *uxr, PetscScalar *pr, PetscScalar *utr, PetscScalar *ubr, PetscScalar *gamr, PetscScalar *rho1r, PetscScalar *rx, PetscScalar *uxm, PetscScalar *px, PetscScalar *utx, PetscScalar *ubx, PetscScalar *gam, PetscScalar *rho1)
583 {
584 /* System generated locals */
585 PetscScalar d__1, d__2;
586
587 /* Local variables */
588 static PetscScalar s, c0, p0, r0, u0, w0, x0, x2, ri, cx, sgn0, wsp0, gasc1, gasc2, gasc3, gasc4;
589 static PetscScalar cstar, pstar, rstar, ustar, xstar, wspst, ushock, streng, rstarl, rstarr, rstars;
590 int iwave;
591
592 if (*rl == *rr && *pr == *pl && *uxl == *uxr && *gaml == *gamr) {
593 *rx = *rl;
594 *px = *pl;
595 *uxm = *uxl;
596 *gam = *gaml;
597 x2 = *xcen + *uxm * *dtt;
598
599 if (PetscRealPart(*xp) >= PetscRealPart(x2)) {
600 *utx = *utr;
601 *ubx = *ubr;
602 *rho1 = *rho1r;
603 } else {
604 *utx = *utl;
605 *ubx = *ubl;
606 *rho1 = *rho1l;
607 }
608 return 0;
609 }
610 iwave = riem1mdt(gaml, gamr, rl, pl, uxl, rr, pr, uxr, &rstarl, &rstarr, &pstar, &ustar);
611
612 x2 = *xcen + ustar * *dtt;
613 d__1 = *xp - x2;
614 sgn0 = sign(d__1);
615 /* x is in 3-wave if sgn0 = 1 */
616 /* x is in 1-wave if sgn0 = -1 */
617 r0 = cvmgm_(rl, rr, &sgn0);
618 p0 = cvmgm_(pl, pr, &sgn0);
619 u0 = cvmgm_(uxl, uxr, &sgn0);
620 *gam = cvmgm_(gaml, gamr, &sgn0);
621 gasc1 = *gam - 1.;
622 gasc2 = (*gam + 1.) * .5;
623 gasc3 = gasc2 / *gam;
624 gasc4 = 1. / (*gam - 1.);
625 c0 = PetscSqrtScalar(*gam * p0 / r0);
626 streng = pstar - p0;
627 w0 = *gam * r0 * p0 * (gasc3 * streng / p0 + 1.);
628 rstars = r0 / (1. - r0 * streng / w0);
629 d__1 = p0 / pstar;
630 d__2 = -1. / *gam;
631 rstarr = r0 * PetscPowScalar(d__1, d__2);
632 rstar = cvmgm_(&rstarr, &rstars, &streng);
633 w0 = PetscSqrtScalar(w0);
634 cstar = PetscSqrtScalar(*gam * pstar / rstar);
635 wsp0 = u0 + sgn0 * c0;
636 wspst = ustar + sgn0 * cstar;
637 ushock = ustar + sgn0 * w0 / rstar;
638 wspst = cvmgp_(&ushock, &wspst, &streng);
639 wsp0 = cvmgp_(&ushock, &wsp0, &streng);
640 x0 = *xcen + wsp0 * *dtt;
641 xstar = *xcen + wspst * *dtt;
642 /* using gas formula to evaluate rarefaction wave */
643 /* ri : reiman invariant */
644 ri = u0 - sgn0 * 2. * gasc4 * c0;
645 cx = sgn0 * .5 * gasc1 / gasc2 * ((*xp - *xcen) / *dtt - ri);
646 *uxm = ri + sgn0 * 2. * gasc4 * cx;
647 s = p0 / PetscPowScalar(r0, *gam);
648 d__1 = cx * cx / (*gam * s);
649 *rx = PetscPowScalar(d__1, gasc4);
650 *px = cx * cx * *rx / *gam;
651 d__1 = sgn0 * (x0 - *xp);
652 *rx = cvmgp_(rx, &r0, &d__1);
653 d__1 = sgn0 * (x0 - *xp);
654 *px = cvmgp_(px, &p0, &d__1);
655 d__1 = sgn0 * (x0 - *xp);
656 *uxm = cvmgp_(uxm, &u0, &d__1);
657 d__1 = sgn0 * (xstar - *xp);
658 *rx = cvmgm_(rx, &rstar, &d__1);
659 d__1 = sgn0 * (xstar - *xp);
660 *px = cvmgm_(px, &pstar, &d__1);
661 d__1 = sgn0 * (xstar - *xp);
662 *uxm = cvmgm_(uxm, &ustar, &d__1);
663 if (PetscRealPart(*xp) >= PetscRealPart(x2)) {
664 *utx = *utr;
665 *ubx = *ubr;
666 *rho1 = *rho1r;
667 } else {
668 *utx = *utl;
669 *ubx = *ubl;
670 *rho1 = *rho1l;
671 }
672 return iwave;
673 }
674
godunovflux(const PetscScalar * ul,const PetscScalar * ur,PetscScalar * flux,const PetscReal * nn,int ndim,PetscReal gamma)675 static int godunovflux(const PetscScalar *ul, const PetscScalar *ur, PetscScalar *flux, const PetscReal *nn, int ndim, PetscReal gamma)
676 {
677 /* System generated locals */
678 int i__1, iwave;
679 PetscScalar d__1, d__2, d__3;
680
681 /* Local variables */
682 static int k;
683 static PetscScalar bn[3], fn, ft, tg[3], pl, rl, pm, pr, rr, xp, ubl, ubm, ubr, dtt, unm, tmp, utl, utm, uxl, utr, uxr, gaml, gamm, gamr, xcen, rhom, rho1l, rho1m, rho1r;
684
685 /* Function Body */
686 xcen = 0.;
687 xp = 0.;
688 i__1 = ndim;
689 for (k = 1; k <= i__1; ++k) {
690 tg[k - 1] = 0.;
691 bn[k - 1] = 0.;
692 }
693 dtt = 1.;
694 if (ndim == 3) {
695 if (nn[0] == 0. && nn[1] == 0.) {
696 tg[0] = 1.;
697 } else {
698 tg[0] = -nn[1];
699 tg[1] = nn[0];
700 }
701 /* tmp=dsqrt(tg(1)**2+tg(2)**2) */
702 /* tg=tg/tmp */
703 bn[0] = -nn[2] * tg[1];
704 bn[1] = nn[2] * tg[0];
705 bn[2] = nn[0] * tg[1] - nn[1] * tg[0];
706 /* Computing 2nd power */
707 d__1 = bn[0];
708 /* Computing 2nd power */
709 d__2 = bn[1];
710 /* Computing 2nd power */
711 d__3 = bn[2];
712 tmp = PetscSqrtScalar(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
713 i__1 = ndim;
714 for (k = 1; k <= i__1; ++k) bn[k - 1] /= tmp;
715 } else if (ndim == 2) {
716 tg[0] = -nn[1];
717 tg[1] = nn[0];
718 /* tmp=dsqrt(tg(1)**2+tg(2)**2) */
719 /* tg=tg/tmp */
720 bn[0] = 0.;
721 bn[1] = 0.;
722 bn[2] = 1.;
723 }
724 rl = ul[0];
725 rr = ur[0];
726 uxl = 0.;
727 uxr = 0.;
728 utl = 0.;
729 utr = 0.;
730 ubl = 0.;
731 ubr = 0.;
732 i__1 = ndim;
733 for (k = 1; k <= i__1; ++k) {
734 uxl += ul[k] * nn[k - 1];
735 uxr += ur[k] * nn[k - 1];
736 utl += ul[k] * tg[k - 1];
737 utr += ur[k] * tg[k - 1];
738 ubl += ul[k] * bn[k - 1];
739 ubr += ur[k] * bn[k - 1];
740 }
741 uxl /= rl;
742 uxr /= rr;
743 utl /= rl;
744 utr /= rr;
745 ubl /= rl;
746 ubr /= rr;
747
748 gaml = gamma;
749 gamr = gamma;
750 /* Computing 2nd power */
751 d__1 = uxl;
752 /* Computing 2nd power */
753 d__2 = utl;
754 /* Computing 2nd power */
755 d__3 = ubl;
756 pl = (gamma - 1.) * (ul[ndim + 1] - rl * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
757 /* Computing 2nd power */
758 d__1 = uxr;
759 /* Computing 2nd power */
760 d__2 = utr;
761 /* Computing 2nd power */
762 d__3 = ubr;
763 pr = (gamma - 1.) * (ur[ndim + 1] - rr * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
764 rho1l = rl;
765 rho1r = rr;
766
767 iwave = riemannsolver(&xcen, &xp, &dtt, &rl, &uxl, &pl, &utl, &ubl, &gaml, &rho1l, &rr, &uxr, &pr, &utr, &ubr, &gamr, &rho1r, &rhom, &unm, &pm, &utm, &ubm, &gamm, &rho1m);
768
769 flux[0] = rhom * unm;
770 fn = rhom * unm * unm + pm;
771 ft = rhom * unm * utm;
772 /* flux(2)=fn*nn(1)+ft*nn(2) */
773 /* flux(3)=fn*tg(1)+ft*tg(2) */
774 flux[1] = fn * nn[0] + ft * tg[0];
775 flux[2] = fn * nn[1] + ft * tg[1];
776 /* flux(2)=rhom*unm*(unm)+pm */
777 /* flux(3)=rhom*(unm)*utm */
778 if (ndim == 3) flux[3] = rhom * unm * ubm;
779 flux[ndim + 1] = (rhom * .5 * (unm * unm + utm * utm + ubm * ubm) + gamm / (gamm - 1.) * pm) * unm;
780 return iwave;
781 } /* godunovflux_ */
782
783 /* PetscReal* => EulerNode* conversion */
PhysicsRiemann_Euler_Godunov(PetscInt dim,PetscInt Nf,const PetscReal * qp,const PetscReal * n,const PetscScalar * xL,const PetscScalar * xR,PetscInt numConstants,const PetscScalar constants[],PetscScalar * flux,Physics phys)784 static void PhysicsRiemann_Euler_Godunov(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys)
785 {
786 Physics_Euler *eu = (Physics_Euler *)phys->data;
787 const PetscReal gamma = eu->gamma;
788 PetscReal zero = 0.;
789 PetscReal cL, cR, speed, velL, velR, nn[DIM], s2;
790 PetscInt i;
791 PetscErrorCode ierr;
792
793 PetscFunctionBeginUser;
794 for (i = 0, s2 = 0.; i < DIM; i++) {
795 nn[i] = n[i];
796 s2 += nn[i] * nn[i];
797 }
798 s2 = PetscSqrtReal(s2); /* |n|_2 = sum(n^2)^1/2 */
799 for (i = 0.; i < DIM; i++) nn[i] /= s2;
800 if (0) { /* Rusanov */
801 const EulerNode *uL = (const EulerNode *)xL, *uR = (const EulerNode *)xR;
802 EulerNodeUnion fL, fR;
803 ierr = EulerFlux(phys, nn, uL, &fL.eulernode);
804 if (ierr) {
805 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
806 for (i = 0; i < 2 + dim; i++) fL.vals[i] = zero / zero;
807 PetscCallVoid(PetscFPTrapPop());
808 }
809 ierr = EulerFlux(phys, nn, uR, &fR.eulernode);
810 if (ierr) {
811 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
812 for (i = 0; i < 2 + dim; i++) fR.vals[i] = zero / zero;
813 PetscCallVoid(PetscFPTrapPop());
814 }
815 ierr = SpeedOfSound_PG(gamma, uL, &cL);
816 if (ierr) {
817 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
818 cL = zero / zero;
819 PetscCallVoid(PetscFPTrapPop());
820 }
821 ierr = SpeedOfSound_PG(gamma, uR, &cR);
822 if (ierr) {
823 PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
824 cR = zero / zero;
825 PetscCallVoid(PetscFPTrapPop());
826 }
827 velL = DotDIMReal(uL->ru, nn) / uL->r;
828 velR = DotDIMReal(uR->ru, nn) / uR->r;
829 speed = PetscMax(velR + cR, velL + cL);
830 for (i = 0; i < 2 + dim; i++) flux[i] = 0.5 * ((fL.vals[i] + fR.vals[i]) + speed * (xL[i] - xR[i])) * s2;
831 } else {
832 /* int iwave = */
833 godunovflux(xL, xR, flux, nn, DIM, gamma);
834 for (i = 0; i < 2 + dim; i++) flux[i] *= s2;
835 }
836 PetscFunctionReturnVoid();
837 }
838
839 #ifdef PETSC_HAVE_LIBCEED
PhysicsRiemann_Euler_Godunov_CEED(PetscCtx ctx,CeedInt Q,const CeedScalar * const in[],CeedScalar * const out[])840 CEED_QFUNCTION(PhysicsRiemann_Euler_Godunov_CEED)(PetscCtx ctx, CeedInt Q, const CeedScalar *const in[], CeedScalar *const out[])
841 {
842 const CeedScalar *xL = in[0], *xR = in[1], *geom = in[2];
843 CeedScalar *cL = out[0], *cR = out[1];
844 const Physics_Euler *eu = (Physics_Euler *)ctx;
845 struct _n_Physics phys;
846
847 phys.data = (void *)eu;
848 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
849 {
850 const CeedScalar qL[DIM + 2] = {xL[i + Q * 0], xL[i + Q * 1], xL[i + Q * 2], xL[i + Q * 3]};
851 const CeedScalar qR[DIM + 2] = {xR[i + Q * 0], xR[i + Q * 1], xR[i + Q * 2], xR[i + Q * 3]};
852 const CeedScalar n[DIM] = {geom[i + Q * 0], geom[i + Q * 1]};
853 CeedScalar flux[DIM + 2];
854
855 #if 0
856 PetscPrintf(PETSC_COMM_SELF, "Cell %d Normal\n", 0);
857 for (CeedInt j = 0; j < DIM; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g |\n", n[j]);
858 PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: left state\n", 0);
859 for (CeedInt j = 0; j < DIM + 2; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g |\n", qL[j]);
860 PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: right state\n", 0);
861 for (CeedInt j = 0; j < DIM + 2; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g |\n", qR[j]);
862 #endif
863 PhysicsRiemann_Euler_Godunov(DIM, DIM + 2, NULL, n, qL, qR, 0, NULL, flux, &phys);
864 for (CeedInt j = 0; j < DIM + 2; ++j) {
865 cL[i + Q * j] = -flux[j] / geom[i + Q * 2];
866 cR[i + Q * j] = flux[j] / geom[i + Q * 3];
867 }
868 #if 0
869 PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: left flux\n", 0);
870 for (CeedInt j = 0; j < DIM + 2; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g | (%g)\n", cL[i + Q * j], geom[i + Q * 2]);
871 PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: right flux\n", 0);
872 for (CeedInt j = 0; j < DIM + 2; ++j) PetscPrintf(PETSC_COMM_SELF, " | %g | (%g)\n", cR[i + Q * j], geom[i + Q * 3]);
873 #endif
874 }
875 return CEED_ERROR_SUCCESS;
876 }
877 #endif
878