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