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