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)(PetscCtx ctx, CeedInt Q, const CeedScalar *const in[], CeedScalar *const out[])
    229: {
    230:   const CeedScalar *xL = in[0], *xR = in[1], *geom = in[2];
    231:   CeedScalar       *cL = out[0], *cR = out[1];
    232:   const Physics_SW *sw = (Physics_SW *)ctx;
    233:   struct _n_Physics phys;
    234:   #if 0
    235:   const CeedScalar *info = in[3];
    236:   #endif

    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)(PetscCtx ctx, CeedInt Q, const CeedScalar *const in[], CeedScalar *const out[])
    841: {
    842:   const CeedScalar    *xL = in[0], *xR = in[1], *geom = in[2];
    843:   CeedScalar          *cL = out[0], *cR = out[1];
    844:   const Physics_Euler *eu = (Physics_Euler *)ctx;
    845:   struct _n_Physics    phys;

    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