Skip to content

Commit

Permalink
address 2d spherical transverse correction (#2961)
Browse files Browse the repository at this point in the history
This adds the necessary pressure correction during transverse correction since our usual flux doesn't contain pressure in spherical2d
  • Loading branch information
zhichen3 authored Sep 19, 2024
1 parent f39525f commit 7e9f5e0
Showing 1 changed file with 16 additions and 2 deletions.
18 changes: 16 additions & 2 deletions Source/hydro/trans.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,8 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member

#if AMREX_SPACEDIM == 2
int coord = geom.Coord();
const auto dx = geom.CellSizeArray();
const auto problo = geom.ProbLoArray();
#endif

bool reset_density = transverse_reset_density;
Expand Down Expand Up @@ -292,16 +294,28 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member
//
// in cylindrical coords -- note that the p term is not
// in a divergence for UMX in the x-direction, so there
// are no area factors. For this geometry, we do not
// are no area factors.
//
// Similarly for Spherical2D geometry, we have:
// d(rho u)/dt + d(rho u v)/(rdtheta) = -1/r^2 d(r^2 rho u u)/dr - dp/dr
// d(rho v)/dt + d(rho u v)/dr = -1/(r sin(theta)) d(sin(theta) rho v v)/dtheta - 1/r dp/dtheta
//
// For these non-cartesian geometries, we do not
// include p in our definition of the flux in the
// x-direction, for we need to fix this now.
// x-direction for Cylindrical2D or both x- and y-direction for spherical 2D
// So we need to fix this now.

Real runewn = run - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMX) -
area_t(il,jl,kl) * flux_t(il,jl,kl,UMX)) * volinv;
if (idir_t == 0 && !mom_flux_has_p(0, idir_t, coord)) {
runewn = runewn - cdtdx * (pgp - pgm);
}
Real rvnewn = rvn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMY) -
area_t(il,jl,kl) * flux_t(il,jl,kl,UMY)) * volinv;
if (idir_t == 1 && !mom_flux_has_p(1, idir_t, coord)) {
Real r = problo[0] + static_cast<Real>(il + 0.5_rt) * dx[0];
rvnewn = rvnewn - cdtdx / r * (pgp - pgm);
}
Real rwnewn = rwn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMZ) -
area_t(il,jl,kl) * flux_t(il,jl,kl,UMZ)) * volinv;
Real renewn = ren - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UEDEN) -
Expand Down

0 comments on commit 7e9f5e0

Please sign in to comment.