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