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