Actual source code: ex11.h
1: #include <petscdm.h>
2: #include <petscdmceed.h>
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
38: #define DIM 2 /* Geometric dimension */
40: /* Represents continuum physical equations. */
41: typedef struct _n_Physics *Physics;
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;
47: struct FieldDescription {
48: const char *name;
49: PetscInt dof;
50: };
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: };
61: typedef struct {
62: PetscReal gravity;
63: struct {
64: PetscInt Height;
65: PetscInt Speed;
66: PetscInt Energy;
67: } functional;
68: } Physics_SW;
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;
80: typedef enum {
81: EULER_IV_SHOCK,
82: EULER_SS_SHOCK,
83: EULER_SHOCK_TUBE,
84: EULER_LINEAR_WAVE
85: } EulerType;
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;
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;
113: static inline PetscReal Dot2Real(const PetscReal *x, const PetscReal *y)
114: {
115: return x[0] * y[0] + x[1] * y[1];
116: }
117: static inline PetscReal Norm2Real(const PetscReal *x)
118: {
119: return PetscSqrtReal(PetscAbsReal(Dot2Real(x, x)));
120: }
121: static inline void Normalize2Real(PetscReal *x)
122: {
123: PetscReal a = 1. / Norm2Real(x);
124: x[0] *= a;
125: x[1] *= a;
126: }
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: }
133: static inline PetscReal DotDIMReal(const PetscReal *x, const PetscReal *y)
134: {
135: PetscInt i;
136: PetscReal prod = 0.0;
138: for (i = 0; i < DIM; i++) prod += x[i] * y[i];
139: return prod;
140: }
141: static inline PetscReal NormDIM(const PetscReal *x)
142: {
143: return PetscSqrtReal(PetscAbsReal(DotDIMReal(x, x)));
144: }
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: }
151: /*
152: * h_t + div(uh) = 0
153: * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
154: *
155: * */
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;
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: }
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;
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
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: }
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) {
224: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", flux[j]);
225: }
226: #endif
227: }
229: #ifdef PETSC_HAVE_LIBCEED
230: CEED_QFUNCTION(PhysicsRiemann_SW_Rusanov_CEED)(void *ctx, CeedInt Q, const CeedScalar *const in[], CeedScalar *const out[])
231: {
232: const CeedScalar *xL = in[0], *xR = in[1], *geom = in[2];
233: CeedScalar *cL = out[0], *cR = out[1];
234: const Physics_SW *sw = (Physics_SW *)ctx;
235: struct _n_Physics phys;
236: #if 0
237: const CeedScalar *info = in[3];
238: #endif
240: phys.data = (void *)sw;
241: CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
242: {
243: const CeedScalar qL[3] = {xL[i + Q * 0], xL[i + Q * 1], xL[i + Q * 2]};
244: const CeedScalar qR[3] = {xR[i + Q * 0], xR[i + Q * 1], xR[i + Q * 2]};
245: const CeedScalar n[2] = {geom[i + Q * 0], geom[i + Q * 1]};
246: CeedScalar flux[3];
248: #if 0
249: PetscPrintf(PETSC_COMM_SELF, "Face %d Normal\n", (int)info[i + Q * 0]);
250: for (CeedInt j = 0; j < DIM; ++j) {
251: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", n[j]);
252: }
253: PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: left state\n", (int)info[i + Q * 1]);
254: for (CeedInt j = 0; j < DIM + 1; ++j) {
255: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", qL[j]);
256: }
257: PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: right state\n", (int)info[i + Q * 2]);
258: for (CeedInt j = 0; j < DIM + 1; ++j) {
259: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", qR[j]);
260: }
261: #endif
262: PhysicsRiemann_SW_Rusanov(DIM, DIM + 1, NULL, n, qL, qR, 0, NULL, flux, &phys);
263: for (CeedInt j = 0; j < 3; ++j) {
264: cL[i + Q * j] = -flux[j] / geom[i + Q * 2];
265: cR[i + Q * j] = flux[j] / geom[i + Q * 3];
266: }
267: #if 0
268: PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: left flux\n", (int)info[i + Q * 1]);
269: for (CeedInt j = 0; j < DIM + 1; ++j) {
270: PetscPrintf(PETSC_COMM_SELF, " | %g | (%g)\n", cL[i + Q * j], geom[i + Q * 2]);
271: }
272: PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: right flux\n", (int)info[i + Q * 2]);
273: for (CeedInt j = 0; j < DIM + 1; ++j) {
274: PetscPrintf(PETSC_COMM_SELF, " | %g | (%g)\n", cR[i + Q * j], geom[i + Q * 3]);
275: }
276: #endif
277: }
278: return CEED_ERROR_SUCCESS;
279: }
280: #endif
282: 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)
283: {
284: Physics_SW *sw = (Physics_SW *)phys->data;
285: PetscReal aL, aR;
286: PetscReal nn[DIM];
287: #if !defined(PETSC_USE_COMPLEX)
288: const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
289: #else
290: SWNodeUnion uLreal, uRreal;
291: const SWNode *uL = &uLreal.swnode;
292: const SWNode *uR = &uRreal.swnode;
293: #endif
294: SWNodeUnion fL, fR;
295: PetscInt i;
296: PetscReal zero = 0.;
297: PetscErrorCode ierr;
299: #if defined(PETSC_USE_COMPLEX)
300: uLreal.swnode.h = 0;
301: uRreal.swnode.h = 0;
302: for (i = 0; i < 1 + dim; i++) uLreal.vals[i] = PetscRealPart(xL[i]);
303: for (i = 0; i < 1 + dim; i++) uRreal.vals[i] = PetscRealPart(xR[i]);
304: #endif
305: if (uL->h <= 0 || uR->h <= 0) {
306: PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
307: for (i = 0; i < 1 + dim; i++) flux[i] = zero;
308: PetscCallVoid(PetscFPTrapPop());
309: return;
310: } /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative"); */
311: nn[0] = n[0];
312: nn[1] = n[1];
313: Normalize2Real(nn);
314: ierr = SWFlux(phys, nn, uL, &fL.swnode);
315: if (ierr) {
316: PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
317: for (i = 0; i < 1 + dim; ++i) fL.vals[i] = zero / zero;
318: PetscCallVoid(PetscFPTrapPop());
319: }
320: ierr = SWFlux(phys, nn, uR, &fR.swnode);
321: if (ierr) {
322: PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
323: for (i = 0; i < 1 + dim; ++i) fR.vals[i] = zero / zero;
324: PetscCallVoid(PetscFPTrapPop());
325: }
326: /* gravity wave speed */
327: aL = PetscSqrtReal(sw->gravity * uL->h);
328: aR = PetscSqrtReal(sw->gravity * uR->h);
329: // Defining u_tilda and v_tilda as u and v
330: PetscReal u_L, u_R;
331: u_L = Dot2Real(uL->uh, nn) / uL->h;
332: u_R = Dot2Real(uR->uh, nn) / uR->h;
333: PetscReal sL, sR;
334: sL = PetscMin(u_L - aL, u_R - aR);
335: sR = PetscMax(u_L + aL, u_R + aR);
336: if (sL > zero) {
337: for (i = 0; i < dim + 1; i++) flux[i] = fL.vals[i] * Norm2Real(n);
338: } else if (sR < zero) {
339: for (i = 0; i < dim + 1; i++) flux[i] = fR.vals[i] * Norm2Real(n);
340: } else {
341: 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);
342: }
343: }
345: static PetscErrorCode Pressure_PG(const PetscReal gamma, const EulerNode *x, PetscReal *p)
346: {
347: PetscReal ru2;
349: PetscFunctionBeginUser;
350: ru2 = DotDIMReal(x->ru, x->ru);
351: (*p) = (x->E - 0.5 * ru2 / x->r) * (gamma - 1.0); /* (E - rho V^2/2)(gamma-1) = e rho (gamma-1) */
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode SpeedOfSound_PG(const PetscReal gamma, const EulerNode *x, PetscReal *c)
356: {
357: PetscReal p;
359: PetscFunctionBeginUser;
360: PetscCall(Pressure_PG(gamma, x, &p));
361: /* gamma = heat capacity ratio */
362: (*c) = PetscSqrtReal(gamma * p / x->r);
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: /*
367: * x = (rho,rho*(u_1),...,rho*e)^T
368: * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
369: *
370: * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
371: *
372: */
373: static PetscErrorCode EulerFlux(Physics phys, const PetscReal *n, const EulerNode *x, EulerNode *f)
374: {
375: Physics_Euler *eu = (Physics_Euler *)phys->data;
376: PetscReal nu, p;
377: PetscInt i;
379: PetscFunctionBeginUser;
380: PetscCall(Pressure_PG(eu->gamma, x, &p));
381: nu = DotDIMReal(x->ru, n);
382: f->r = nu; /* A rho u */
383: nu /= x->r; /* A u */
384: for (i = 0; i < DIM; i++) f->ru[i] = nu * x->ru[i] + n[i] * p; /* r u^2 + p */
385: f->E = nu * (x->E + p); /* u(e+p) */
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: /* Godunov fluxs */
390: static PetscScalar cvmgp_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
391: {
392: /* System generated locals */
393: PetscScalar ret_val;
395: if (PetscRealPart(*test) > 0.) goto L10;
396: ret_val = *b;
397: return ret_val;
398: L10:
399: ret_val = *a;
400: return ret_val;
401: } /* cvmgp_ */
403: static PetscScalar cvmgm_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
404: {
405: /* System generated locals */
406: PetscScalar ret_val;
408: if (PetscRealPart(*test) < 0.) goto L10;
409: ret_val = *b;
410: return ret_val;
411: L10:
412: ret_val = *a;
413: return ret_val;
414: } /* cvmgm_ */
416: 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)
417: {
418: /* Initialized data */
420: static PetscScalar smallp = 1e-8;
422: /* System generated locals */
423: int i__1;
424: PetscScalar d__1, d__2;
426: /* Local variables */
427: static int i0;
428: static PetscScalar cl, cr, wl, zl, wr, zr, pst, durl, skpr1, skpr2;
429: static int iwave;
430: static PetscScalar gascl4, gascr4, cstarl, dpstar, cstarr;
431: /* static PetscScalar csqrl, csqrr, gascl1, gascl2, gascl3, gascr1, gascr2, gascr3; */
432: static int iterno;
433: static PetscScalar ustarl, ustarr, rarepr1, rarepr2;
435: /* gascl1 = *gaml - 1.; */
436: /* gascl2 = (*gaml + 1.) * .5; */
437: /* gascl3 = gascl2 / *gaml; */
438: gascl4 = 1. / (*gaml - 1.);
440: /* gascr1 = *gamr - 1.; */
441: /* gascr2 = (*gamr + 1.) * .5; */
442: /* gascr3 = gascr2 / *gamr; */
443: gascr4 = 1. / (*gamr - 1.);
444: iterno = 10;
445: /* find pstar: */
446: cl = PetscSqrtScalar(*gaml * *pl / *rl);
447: cr = PetscSqrtScalar(*gamr * *pr / *rr);
448: wl = *rl * cl;
449: wr = *rr * cr;
450: /* csqrl = wl * wl; */
451: /* csqrr = wr * wr; */
452: *pstar = (wl * *pr + wr * *pl) / (wl + wr);
453: *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
454: pst = *pl / *pr;
455: skpr1 = cr * (pst - 1.) * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
456: d__1 = (*gamr - 1.) / (*gamr * 2.);
457: rarepr2 = gascr4 * 2. * cr * (1. - PetscPowScalar(pst, d__1));
458: pst = *pr / *pl;
459: skpr2 = cl * (pst - 1.) * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
460: d__1 = (*gaml - 1.) / (*gaml * 2.);
461: rarepr1 = gascl4 * 2. * cl * (1. - PetscPowScalar(pst, d__1));
462: durl = *uxr - *uxl;
463: if (PetscRealPart(*pr) < PetscRealPart(*pl)) {
464: if (PetscRealPart(durl) >= PetscRealPart(rarepr1)) {
465: iwave = 100;
466: } else if (PetscRealPart(durl) <= PetscRealPart(-skpr1)) {
467: iwave = 300;
468: } else {
469: iwave = 400;
470: }
471: } else {
472: if (PetscRealPart(durl) >= PetscRealPart(rarepr2)) {
473: iwave = 100;
474: } else if (PetscRealPart(durl) <= PetscRealPart(-skpr2)) {
475: iwave = 300;
476: } else {
477: iwave = 200;
478: }
479: }
480: if (iwave == 100) {
481: /* 1-wave: rarefaction wave, 3-wave: rarefaction wave */
482: /* case (100) */
483: i__1 = iterno;
484: for (i0 = 1; i0 <= i__1; ++i0) {
485: d__1 = *pstar / *pl;
486: d__2 = 1. / *gaml;
487: *rstarl = *rl * PetscPowScalar(d__1, d__2);
488: cstarl = PetscSqrtScalar(*gaml * *pstar / *rstarl);
489: ustarl = *uxl - gascl4 * 2. * (cstarl - cl);
490: zl = *rstarl * cstarl;
491: d__1 = *pstar / *pr;
492: d__2 = 1. / *gamr;
493: *rstarr = *rr * PetscPowScalar(d__1, d__2);
494: cstarr = PetscSqrtScalar(*gamr * *pstar / *rstarr);
495: ustarr = *uxr + gascr4 * 2. * (cstarr - cr);
496: zr = *rstarr * cstarr;
497: dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
498: *pstar -= dpstar;
499: *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
500: if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
501: #if 0
502: break;
503: #endif
504: }
505: }
506: /* 1-wave: shock wave, 3-wave: rarefaction wave */
507: } else if (iwave == 200) {
508: /* case (200) */
509: i__1 = iterno;
510: for (i0 = 1; i0 <= i__1; ++i0) {
511: pst = *pstar / *pl;
512: ustarl = *uxl - (pst - 1.) * cl * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
513: zl = *pl / cl * PetscSqrtScalar(*gaml * 2. * (*gaml - 1. + (*gaml + 1.) * pst)) * (*gaml - 1. + (*gaml + 1.) * pst) / (*gaml * 3. - 1. + (*gaml + 1.) * pst);
514: d__1 = *pstar / *pr;
515: d__2 = 1. / *gamr;
516: *rstarr = *rr * PetscPowScalar(d__1, d__2);
517: cstarr = PetscSqrtScalar(*gamr * *pstar / *rstarr);
518: zr = *rstarr * cstarr;
519: ustarr = *uxr + gascr4 * 2. * (cstarr - cr);
520: dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
521: *pstar -= dpstar;
522: *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
523: if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
524: #if 0
525: break;
526: #endif
527: }
528: }
529: /* 1-wave: shock wave, 3-wave: shock */
530: } else if (iwave == 300) {
531: /* case (300) */
532: i__1 = iterno;
533: for (i0 = 1; i0 <= i__1; ++i0) {
534: pst = *pstar / *pl;
535: ustarl = *uxl - (pst - 1.) * cl * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
536: zl = *pl / cl * PetscSqrtScalar(*gaml * 2. * (*gaml - 1. + (*gaml + 1.) * pst)) * (*gaml - 1. + (*gaml + 1.) * pst) / (*gaml * 3. - 1. + (*gaml + 1.) * pst);
537: pst = *pstar / *pr;
538: ustarr = *uxr + (pst - 1.) * cr * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
539: zr = *pr / cr * PetscSqrtScalar(*gamr * 2. * (*gamr - 1. + (*gamr + 1.) * pst)) * (*gamr - 1. + (*gamr + 1.) * pst) / (*gamr * 3. - 1. + (*gamr + 1.) * pst);
540: dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
541: *pstar -= dpstar;
542: *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
543: if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
544: #if 0
545: break;
546: #endif
547: }
548: }
549: /* 1-wave: rarefaction wave, 3-wave: shock */
550: } else if (iwave == 400) {
551: /* case (400) */
552: i__1 = iterno;
553: for (i0 = 1; i0 <= i__1; ++i0) {
554: d__1 = *pstar / *pl;
555: d__2 = 1. / *gaml;
556: *rstarl = *rl * PetscPowScalar(d__1, d__2);
557: cstarl = PetscSqrtScalar(*gaml * *pstar / *rstarl);
558: ustarl = *uxl - gascl4 * 2. * (cstarl - cl);
559: zl = *rstarl * cstarl;
560: pst = *pstar / *pr;
561: ustarr = *uxr + (pst - 1.) * cr * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
562: zr = *pr / cr * PetscSqrtScalar(*gamr * 2. * (*gamr - 1. + (*gamr + 1.) * pst)) * (*gamr - 1. + (*gamr + 1.) * pst) / (*gamr * 3. - 1. + (*gamr + 1.) * pst);
563: dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
564: *pstar -= dpstar;
565: *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
566: if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
567: #if 0
568: break;
569: #endif
570: }
571: }
572: }
574: *ustar = (zl * ustarr + zr * ustarl) / (zl + zr);
575: if (PetscRealPart(*pstar) > PetscRealPart(*pl)) {
576: pst = *pstar / *pl;
577: *rstarl = ((*gaml + 1.) * pst + *gaml - 1.) / ((*gaml - 1.) * pst + *gaml + 1.) * *rl;
578: }
579: if (PetscRealPart(*pstar) > PetscRealPart(*pr)) {
580: pst = *pstar / *pr;
581: *rstarr = ((*gamr + 1.) * pst + *gamr - 1.) / ((*gamr - 1.) * pst + *gamr + 1.) * *rr;
582: }
583: return iwave;
584: }
586: static PetscScalar sign(PetscScalar x)
587: {
588: if (PetscRealPart(x) > 0) return 1.0;
589: if (PetscRealPart(x) < 0) return -1.0;
590: return 0.0;
591: }
592: /* Riemann Solver */
593: /* -------------------------------------------------------------------- */
594: 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)
595: {
596: /* System generated locals */
597: PetscScalar d__1, d__2;
599: /* Local variables */
600: static PetscScalar s, c0, p0, r0, u0, w0, x0, x2, ri, cx, sgn0, wsp0, gasc1, gasc2, gasc3, gasc4;
601: static PetscScalar cstar, pstar, rstar, ustar, xstar, wspst, ushock, streng, rstarl, rstarr, rstars;
602: int iwave;
604: if (*rl == *rr && *pr == *pl && *uxl == *uxr && *gaml == *gamr) {
605: *rx = *rl;
606: *px = *pl;
607: *uxm = *uxl;
608: *gam = *gaml;
609: x2 = *xcen + *uxm * *dtt;
611: if (PetscRealPart(*xp) >= PetscRealPart(x2)) {
612: *utx = *utr;
613: *ubx = *ubr;
614: *rho1 = *rho1r;
615: } else {
616: *utx = *utl;
617: *ubx = *ubl;
618: *rho1 = *rho1l;
619: }
620: return 0;
621: }
622: iwave = riem1mdt(gaml, gamr, rl, pl, uxl, rr, pr, uxr, &rstarl, &rstarr, &pstar, &ustar);
624: x2 = *xcen + ustar * *dtt;
625: d__1 = *xp - x2;
626: sgn0 = sign(d__1);
627: /* x is in 3-wave if sgn0 = 1 */
628: /* x is in 1-wave if sgn0 = -1 */
629: r0 = cvmgm_(rl, rr, &sgn0);
630: p0 = cvmgm_(pl, pr, &sgn0);
631: u0 = cvmgm_(uxl, uxr, &sgn0);
632: *gam = cvmgm_(gaml, gamr, &sgn0);
633: gasc1 = *gam - 1.;
634: gasc2 = (*gam + 1.) * .5;
635: gasc3 = gasc2 / *gam;
636: gasc4 = 1. / (*gam - 1.);
637: c0 = PetscSqrtScalar(*gam * p0 / r0);
638: streng = pstar - p0;
639: w0 = *gam * r0 * p0 * (gasc3 * streng / p0 + 1.);
640: rstars = r0 / (1. - r0 * streng / w0);
641: d__1 = p0 / pstar;
642: d__2 = -1. / *gam;
643: rstarr = r0 * PetscPowScalar(d__1, d__2);
644: rstar = cvmgm_(&rstarr, &rstars, &streng);
645: w0 = PetscSqrtScalar(w0);
646: cstar = PetscSqrtScalar(*gam * pstar / rstar);
647: wsp0 = u0 + sgn0 * c0;
648: wspst = ustar + sgn0 * cstar;
649: ushock = ustar + sgn0 * w0 / rstar;
650: wspst = cvmgp_(&ushock, &wspst, &streng);
651: wsp0 = cvmgp_(&ushock, &wsp0, &streng);
652: x0 = *xcen + wsp0 * *dtt;
653: xstar = *xcen + wspst * *dtt;
654: /* using gas formula to evaluate rarefaction wave */
655: /* ri : reiman invariant */
656: ri = u0 - sgn0 * 2. * gasc4 * c0;
657: cx = sgn0 * .5 * gasc1 / gasc2 * ((*xp - *xcen) / *dtt - ri);
658: *uxm = ri + sgn0 * 2. * gasc4 * cx;
659: s = p0 / PetscPowScalar(r0, *gam);
660: d__1 = cx * cx / (*gam * s);
661: *rx = PetscPowScalar(d__1, gasc4);
662: *px = cx * cx * *rx / *gam;
663: d__1 = sgn0 * (x0 - *xp);
664: *rx = cvmgp_(rx, &r0, &d__1);
665: d__1 = sgn0 * (x0 - *xp);
666: *px = cvmgp_(px, &p0, &d__1);
667: d__1 = sgn0 * (x0 - *xp);
668: *uxm = cvmgp_(uxm, &u0, &d__1);
669: d__1 = sgn0 * (xstar - *xp);
670: *rx = cvmgm_(rx, &rstar, &d__1);
671: d__1 = sgn0 * (xstar - *xp);
672: *px = cvmgm_(px, &pstar, &d__1);
673: d__1 = sgn0 * (xstar - *xp);
674: *uxm = cvmgm_(uxm, &ustar, &d__1);
675: if (PetscRealPart(*xp) >= PetscRealPart(x2)) {
676: *utx = *utr;
677: *ubx = *ubr;
678: *rho1 = *rho1r;
679: } else {
680: *utx = *utl;
681: *ubx = *ubl;
682: *rho1 = *rho1l;
683: }
684: return iwave;
685: }
687: static int godunovflux(const PetscScalar *ul, const PetscScalar *ur, PetscScalar *flux, const PetscReal *nn, int ndim, PetscReal gamma)
688: {
689: /* System generated locals */
690: int i__1, iwave;
691: PetscScalar d__1, d__2, d__3;
693: /* Local variables */
694: static int k;
695: 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;
697: /* Function Body */
698: xcen = 0.;
699: xp = 0.;
700: i__1 = ndim;
701: for (k = 1; k <= i__1; ++k) {
702: tg[k - 1] = 0.;
703: bn[k - 1] = 0.;
704: }
705: dtt = 1.;
706: if (ndim == 3) {
707: if (nn[0] == 0. && nn[1] == 0.) {
708: tg[0] = 1.;
709: } else {
710: tg[0] = -nn[1];
711: tg[1] = nn[0];
712: }
713: /* tmp=dsqrt(tg(1)**2+tg(2)**2) */
714: /* tg=tg/tmp */
715: bn[0] = -nn[2] * tg[1];
716: bn[1] = nn[2] * tg[0];
717: bn[2] = nn[0] * tg[1] - nn[1] * tg[0];
718: /* Computing 2nd power */
719: d__1 = bn[0];
720: /* Computing 2nd power */
721: d__2 = bn[1];
722: /* Computing 2nd power */
723: d__3 = bn[2];
724: tmp = PetscSqrtScalar(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
725: i__1 = ndim;
726: for (k = 1; k <= i__1; ++k) bn[k - 1] /= tmp;
727: } else if (ndim == 2) {
728: tg[0] = -nn[1];
729: tg[1] = nn[0];
730: /* tmp=dsqrt(tg(1)**2+tg(2)**2) */
731: /* tg=tg/tmp */
732: bn[0] = 0.;
733: bn[1] = 0.;
734: bn[2] = 1.;
735: }
736: rl = ul[0];
737: rr = ur[0];
738: uxl = 0.;
739: uxr = 0.;
740: utl = 0.;
741: utr = 0.;
742: ubl = 0.;
743: ubr = 0.;
744: i__1 = ndim;
745: for (k = 1; k <= i__1; ++k) {
746: uxl += ul[k] * nn[k - 1];
747: uxr += ur[k] * nn[k - 1];
748: utl += ul[k] * tg[k - 1];
749: utr += ur[k] * tg[k - 1];
750: ubl += ul[k] * bn[k - 1];
751: ubr += ur[k] * bn[k - 1];
752: }
753: uxl /= rl;
754: uxr /= rr;
755: utl /= rl;
756: utr /= rr;
757: ubl /= rl;
758: ubr /= rr;
760: gaml = gamma;
761: gamr = gamma;
762: /* Computing 2nd power */
763: d__1 = uxl;
764: /* Computing 2nd power */
765: d__2 = utl;
766: /* Computing 2nd power */
767: d__3 = ubl;
768: pl = (gamma - 1.) * (ul[ndim + 1] - rl * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
769: /* Computing 2nd power */
770: d__1 = uxr;
771: /* Computing 2nd power */
772: d__2 = utr;
773: /* Computing 2nd power */
774: d__3 = ubr;
775: pr = (gamma - 1.) * (ur[ndim + 1] - rr * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
776: rho1l = rl;
777: rho1r = rr;
779: 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);
781: flux[0] = rhom * unm;
782: fn = rhom * unm * unm + pm;
783: ft = rhom * unm * utm;
784: /* flux(2)=fn*nn(1)+ft*nn(2) */
785: /* flux(3)=fn*tg(1)+ft*tg(2) */
786: flux[1] = fn * nn[0] + ft * tg[0];
787: flux[2] = fn * nn[1] + ft * tg[1];
788: /* flux(2)=rhom*unm*(unm)+pm */
789: /* flux(3)=rhom*(unm)*utm */
790: if (ndim == 3) flux[3] = rhom * unm * ubm;
791: flux[ndim + 1] = (rhom * .5 * (unm * unm + utm * utm + ubm * ubm) + gamm / (gamm - 1.) * pm) * unm;
792: return iwave;
793: } /* godunovflux_ */
795: /* PetscReal* => EulerNode* conversion */
796: 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)
797: {
798: Physics_Euler *eu = (Physics_Euler *)phys->data;
799: const PetscReal gamma = eu->gamma;
800: PetscReal zero = 0.;
801: PetscReal cL, cR, speed, velL, velR, nn[DIM], s2;
802: PetscInt i;
803: PetscErrorCode ierr;
805: PetscFunctionBeginUser;
806: for (i = 0, s2 = 0.; i < DIM; i++) {
807: nn[i] = n[i];
808: s2 += nn[i] * nn[i];
809: }
810: s2 = PetscSqrtReal(s2); /* |n|_2 = sum(n^2)^1/2 */
811: for (i = 0.; i < DIM; i++) nn[i] /= s2;
812: if (0) { /* Rusanov */
813: const EulerNode *uL = (const EulerNode *)xL, *uR = (const EulerNode *)xR;
814: EulerNodeUnion fL, fR;
815: ierr = EulerFlux(phys, nn, uL, &fL.eulernode);
816: if (ierr) {
817: PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
818: for (i = 0; i < 2 + dim; i++) fL.vals[i] = zero / zero;
819: PetscCallVoid(PetscFPTrapPop());
820: }
821: ierr = EulerFlux(phys, nn, uR, &fR.eulernode);
822: if (ierr) {
823: PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
824: for (i = 0; i < 2 + dim; i++) fR.vals[i] = zero / zero;
825: PetscCallVoid(PetscFPTrapPop());
826: }
827: ierr = SpeedOfSound_PG(gamma, uL, &cL);
828: if (ierr) {
829: PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
830: cL = zero / zero;
831: PetscCallVoid(PetscFPTrapPop());
832: }
833: ierr = SpeedOfSound_PG(gamma, uR, &cR);
834: if (ierr) {
835: PetscCallVoid(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
836: cR = zero / zero;
837: PetscCallVoid(PetscFPTrapPop());
838: }
839: velL = DotDIMReal(uL->ru, nn) / uL->r;
840: velR = DotDIMReal(uR->ru, nn) / uR->r;
841: speed = PetscMax(velR + cR, velL + cL);
842: for (i = 0; i < 2 + dim; i++) flux[i] = 0.5 * ((fL.vals[i] + fR.vals[i]) + speed * (xL[i] - xR[i])) * s2;
843: } else {
844: /* int iwave = */
845: godunovflux(xL, xR, flux, nn, DIM, gamma);
846: for (i = 0; i < 2 + dim; i++) flux[i] *= s2;
847: }
848: PetscFunctionReturnVoid();
849: }
851: #ifdef PETSC_HAVE_LIBCEED
852: CEED_QFUNCTION(PhysicsRiemann_Euler_Godunov_CEED)(void *ctx, CeedInt Q, const CeedScalar *const in[], CeedScalar *const out[])
853: {
854: const CeedScalar *xL = in[0], *xR = in[1], *geom = in[2];
855: CeedScalar *cL = out[0], *cR = out[1];
856: const Physics_Euler *eu = (Physics_Euler *)ctx;
857: struct _n_Physics phys;
859: phys.data = (void *)eu;
860: CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
861: {
862: const CeedScalar qL[DIM + 2] = {xL[i + Q * 0], xL[i + Q * 1], xL[i + Q * 2], xL[i + Q * 3]};
863: const CeedScalar qR[DIM + 2] = {xR[i + Q * 0], xR[i + Q * 1], xR[i + Q * 2], xR[i + Q * 3]};
864: const CeedScalar n[DIM] = {geom[i + Q * 0], geom[i + Q * 1]};
865: CeedScalar flux[DIM + 2];
867: #if 0
868: PetscPrintf(PETSC_COMM_SELF, "Cell %d Normal\n", 0);
869: for (CeedInt j = 0; j < DIM; ++j) {
870: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", n[j]);
871: }
872: PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: left state\n", 0);
873: for (CeedInt j = 0; j < DIM + 2; ++j) {
874: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", qL[j]);
875: }
876: PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: right state\n", 0);
877: for (CeedInt j = 0; j < DIM + 2; ++j) {
878: PetscPrintf(PETSC_COMM_SELF, " | %g |\n", qR[j]);
879: }
880: #endif
881: PhysicsRiemann_Euler_Godunov(DIM, DIM + 2, NULL, n, qL, qR, 0, NULL, flux, &phys);
882: for (CeedInt j = 0; j < DIM + 2; ++j) {
883: cL[i + Q * j] = -flux[j] / geom[i + Q * 2];
884: cR[i + Q * j] = flux[j] / geom[i + Q * 3];
885: }
886: #if 0
887: PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: left flux\n", 0);
888: for (CeedInt j = 0; j < DIM + 2; ++j) {
889: PetscPrintf(PETSC_COMM_SELF, " | %g | (%g)\n", cL[i + Q * j], geom[i + Q * 2]);
890: }
891: PetscPrintf(PETSC_COMM_SELF, "Cell %d Element Residual: right flux\n", 0);
892: for (CeedInt j = 0; j < DIM + 2; ++j) {
893: PetscPrintf(PETSC_COMM_SELF, " | %g | (%g)\n", cR[i + Q * j], geom[i + Q * 3]);
894: }
895: #endif
896: }
897: return CEED_ERROR_SUCCESS;
898: }
899: #endif