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) PetscPrintf(PETSC_COMM_SELF, " | %g |\n", flux[j]);
224: #endif
225: }
227: #ifdef PETSC_HAVE_LIBCEED
228: CEED_QFUNCTION(PhysicsRiemann_SW_Rusanov_CEED)(void *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
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];
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
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;
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: }
333: static PetscErrorCode Pressure_PG(const PetscReal gamma, const EulerNode *x, PetscReal *p)
334: {
335: PetscReal ru2;
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: }
343: static PetscErrorCode SpeedOfSound_PG(const PetscReal gamma, const EulerNode *x, PetscReal *c)
344: {
345: PetscReal p;
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: }
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: */
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;
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: }
377: /* Godunov fluxs */
378: static PetscScalar cvmgp_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
379: {
380: /* System generated locals */
381: PetscScalar ret_val;
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_ */
391: static PetscScalar cvmgm_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
392: {
393: /* System generated locals */
394: PetscScalar ret_val;
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_ */
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 */
408: static PetscScalar smallp = 1e-8;
410: /* System generated locals */
411: int i__1;
412: PetscScalar d__1, d__2;
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;
423: /* gascl1 = *gaml - 1.; */
424: /* gascl2 = (*gaml + 1.) * .5; */
425: /* gascl3 = gascl2 / *gaml; */
426: gascl4 = 1. / (*gaml - 1.);
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: }
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: }
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: /* -------------------------------------------------------------------- */
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;
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;
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;
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);
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: }
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;
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;
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;
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;
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);
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_ */
783: /* PetscReal* => EulerNode* conversion */
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;
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: }
839: #ifdef PETSC_HAVE_LIBCEED
840: CEED_QFUNCTION(PhysicsRiemann_Euler_Godunov_CEED)(void *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;
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];
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