Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add eos information to floors, misc UserWorkBeforeOutput fix #210

Merged
merged 14 commits into from
Mar 20, 2024
Merged

Add eos information to floors, misc UserWorkBeforeOutput fix #210

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
55be7e9
feat: threaded eos bounds into floors/ceilings. tweak ye logic
AstroBarker Mar 14, 2024
7e78834
leftover print
AstroBarker Mar 14, 2024
4e1391a
feat: pull ideal gas eos bounds from <eos> not from fixup
AstroBarker Mar 14, 2024
003d75c
bounds in moments_source
AstroBarker Mar 14, 2024
ed064de
feat: call SetEOSBnds
AstroBarker Mar 19, 2024
21a1753
fix: proper copy by value for Bounds
AstroBarker Mar 20, 2024
1162e41
feat: remove eos from Ceilings
AstroBarker Mar 20, 2024
f1af909
fix: pull package outside of kernel in UserWorkBeforeOutput
AstroBarker Mar 20, 2024
72a317e
fix: missed some -> to .
AstroBarker Mar 20, 2024
26b5971
fix: some missed bounds copies
AstroBarker Mar 20, 2024
ad10726
fix/feat: copy bounds into FluxState for riemann solve
AstroBarker Mar 20, 2024
ecd49fa
fix: remove unused vars
AstroBarker Mar 20, 2024
e445b09
fix: revert mystery unintended changes....
AstroBarker Mar 20, 2024
e6d5b4a
merge main fluid.hpp ....
AstroBarker Mar 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 31 additions & 21 deletions src/fixup/fixup.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ namespace fixup {
std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
auto fix = std::make_shared<StateDescriptor>("fixup");
Params &params = fix->AllParams();
auto mutable_param = parthenon::Params::Mutability::Mutable;

const bool enable_mhd = pin->GetOrAddBoolean("fluid", "mhd", false);
const bool enable_rad = pin->GetOrAddBoolean("physics", "rad", false);
Expand Down Expand Up @@ -273,7 +274,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
Bounds(params.Get<Floors>("floor"), params.Get<Ceilings>("ceiling"),
params.Get<MHDCeilings>("mhd_ceiling"),
params.Get<RadiationFloors>("rad_floor"),
params.Get<RadiationCeilings>("rad_ceiling")));
params.Get<RadiationCeilings>("rad_ceiling")),
mutable_param);

return fix;
}
Expand Down Expand Up @@ -307,6 +309,9 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {

if (!enable_floors) return TaskStatus::complete;

const Real ye_min = eos_pkg->Param<Real>("ye_min");
const Real ye_max = eos_pkg->Param<Real>("ye_max");

const std::vector<std::string> vars(
{p::density::name(), c::density::name(), p::velocity::name(), c::momentum::name(),
p::energy::name(), c::energy::name(), p::bfield::name(), p::ye::name(),
Expand Down Expand Up @@ -344,7 +349,10 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {

auto eos = eos_pkg->Param<EOS>("d.EOS");
auto geom = Geometry::GetCoordinateSystem(rc);
auto bounds = fix_pkg->Param<Bounds>("bounds");
Bounds *bounds = fix_pkg->MutableParam<Bounds>("bounds");

// BLB: Setting EOS bnds for Ceilings/Floors here.
bounds->SetEOSBnds(eos_pkg);
AstroBarker marked this conversation as resolved.
Show resolved Hide resolved

const Real c2p_tol = fluid_pkg->Param<Real>("c2p_tol");
const int c2p_max_iter = fluid_pkg->Param<int>("c2p_max_iter");
Expand All @@ -367,21 +375,21 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
eos_lambda[1] = std::log10(v(b, tmp, k, j, i)); // use last temp as initial guess

double rho_floor, sie_floor;
bounds.GetFloors(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), rho_floor, sie_floor);
bounds->GetFloors(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), rho_floor, sie_floor);
AstroBarker marked this conversation as resolved.
Show resolved Hide resolved
double gamma_max, e_max;
bounds.GetCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), gamma_max, e_max);
bounds->GetCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), gamma_max, e_max);
Real bsqorho_max, bsqou_max;
bounds.GetMHDCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), bsqorho_max, bsqou_max);
bounds->GetMHDCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), bsqorho_max, bsqou_max);
Real J_floor;
bounds.GetRadiationFloors(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), J_floor);
bounds->GetRadiationFloors(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), J_floor);
Real xi_max;
Real garbage;
bounds.GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), xi_max, garbage);
bounds->GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), xi_max, garbage);

Real rho_floor_max = rho_floor;
Real u_floor_max = rho_floor * sie_floor;
Expand Down Expand Up @@ -437,16 +445,18 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {

if (floor_applied) {
Real vp_normalobs[3] = {0}; // Inject floors at rest in normal observer frame
Real ye_prim_default = 0.5;
eos_lambda[0] = ye_prim_default;
Real ye = 0.5;
if (pye > 0) {
ye = v(b, cye, k, j, i);
}
eos_lambda[0] = ye;
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
Real dprs =
eos.PressureFromDensityInternalEnergy(drho, ratio(du, drho), eos_lambda);
Real dgm1 = ratio(
eos.BulkModulusFromDensityInternalEnergy(drho, ratio(du, drho), eos_lambda),
dprs);
prim2con::p2c(drho, vp_normalobs, bp, du, ye_prim_default, dprs, dgm1, gcov,
gammacon, betacon, alpha, sdetgam, dcrho, dS, dBcons, dtau,
dyecons);
prim2con::p2c(drho, vp_normalobs, bp, du, ye, dprs, dgm1, gcov, gammacon,
betacon, alpha, sdetgam, dcrho, dS, dBcons, dtau, dyecons);

// Update cons vars (not B field)
v(b, crho, k, j, i) += dcrho;
Expand All @@ -464,7 +474,7 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
SPACELOOP(ii) { v(b, pvel_lo + ii, k, j, i) = vp_normalobs[ii]; }
v(b, peng, k, j, i) = du;
if (pye > 0) {
v(b, pye, k, j, i) = ye_prim_default;
v(b, pye, k, j, i) = ye_min;
}

// Update auxiliary primitives
Expand All @@ -476,9 +486,9 @@ TaskStatus ApplyFloorsImpl(T *rc, IndexDomain domain = IndexDomain::entire) {
eos_lambda);

// Update cons vars (not B field)
prim2con::p2c(drho, vp_normalobs, bp, du, ye_prim_default, dprs, dgm1, gcov,
gammacon, betacon, alpha, sdetgam, v(b, crho, k, j, i), dS,
dBcons, v(b, ceng, k, j, i), dyecons);
prim2con::p2c(drho, vp_normalobs, bp, du, ye, dprs, dgm1, gcov, gammacon,
betacon, alpha, sdetgam, v(b, crho, k, j, i), dS, dBcons,
v(b, ceng, k, j, i), dyecons);
SPACELOOP(ii) { v(b, cmom_lo + ii, k, j, i) = dS[ii]; }
if (pye > 0) {
v(b, cye, k, j, i) = dyecons;
Expand Down
54 changes: 47 additions & 7 deletions src/fixup/fixup.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,37 +91,60 @@ class Floors {
KOKKOS_INLINE_FUNCTION
void GetFloors(const Real x1, const Real x2, const Real x3, Real &rflr,
Real &sflr) const {
if (!eos_bnds_set_) {
PARTHENON_FAIL("EOS bounds not set in floors.");
}
AstroBarker marked this conversation as resolved.
Show resolved Hide resolved

switch (floor_flag_) {
case FloorFlag::ConstantRhoSie:
rflr = r0_;
sflr = s0_;
rflr = std::max(rflr, rho_min_eos_);
sflr = std::max(sflr, sie_min_eos_);
break;
case FloorFlag::ExpX1RhoSie:
rflr = r0_ * std::exp(ralpha_ * x1);
sflr = s0_ * std::exp(salpha_ * x1);
rflr = std::max(rflr, rho_min_eos_);
sflr = std::max(sflr, sie_min_eos_);
break;
case FloorFlag::ExpX1RhoU: {
Real scratch = r0_ * std::exp(ralpha_ * x1);
sflr = s0_ * std::exp(salpha_ * x1) / std::max(rflr, scratch);
rflr = scratch;
rflr = std::max(rflr, rho_min_eos_);
sflr = std::max(sflr, sie_min_eos_);
} break;
case FloorFlag::X1RhoSie:
rflr = r0_ * std::min(1.0, std::pow(x1, ralpha_));
sflr = s0_ * std::min(1.0, std::pow(x1, salpha_));
rflr = std::max(rflr, rho_min_eos_);
sflr = std::max(sflr, sie_min_eos_);
break;
case FloorFlag::RRhoSie: {
Real r = std::sqrt(x1 * x1 + x2 * x2 + x3 * x3);
rflr = r0_ * std::min(1.0, std::pow(r, ralpha_));
sflr = s0_ * std::min(1.0, std::pow(r, salpha_));
rflr = std::max(rflr, rho_min_eos_);
sflr = std::max(sflr, sie_min_eos_);
} break;
default:
PARTHENON_FAIL("No valid floor set.");
}
}

void SetEOSBnds(StateDescriptor *eos_pkg) {
if (!eos_bnds_set_) {
rho_min_eos_ = eos_pkg->Param<Real>("rho_min");
sie_min_eos_ = eos_pkg->Param<Real>("sie_min");
eos_bnds_set_ = true;
}
}

private:
Real r0_, s0_, ralpha_, salpha_;
Real r0_, s0_, ralpha_, salpha_, rho_min_eos_, sie_min_eos_;
const FloorFlag floor_flag_;
bool eos_bnds_set_;
};

static struct ConstantJFloor {
Expand Down Expand Up @@ -173,19 +196,30 @@ class Ceilings {
KOKKOS_INLINE_FUNCTION
void GetCeilings(const Real x1, const Real x2, const Real x3, Real &gmax,
Real &smax) const {
if (!eos_bnds_set_) {
PARTHENON_FAIL("EOS bounds not set in ceilings.");
}
AstroBarker marked this conversation as resolved.
Show resolved Hide resolved
switch (ceiling_flag_) {
case 1:
gmax = g0_;
smax = s0_;
smax = std::min(s0_, sie_max_eos_);
break;
default:
PARTHENON_FAIL("No valid ceiling set.");
}
}

void SetEOSBnds(StateDescriptor *eos_pkg) {
if (!eos_bnds_set_) {
sie_max_eos_ = eos_pkg->Param<Real>("sie_max");
eos_bnds_set_ = true;
}
}

private:
Real g0_, s0_;
Real g0_, s0_, sie_max_eos_;
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
const int ceiling_flag_;
bool eos_bnds_set_;
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
};

static struct ConstantBsqRatCeiling {
Expand Down Expand Up @@ -257,13 +291,19 @@ class Bounds {
const RadiationFloors &rfl, const RadiationCeilings &rcl)
: floors_(fl), ceilings_(cl), mhd_ceilings_(mcl), radiation_floors_(rfl),
radiation_ceilings_(rcl) {}
explicit Bounds(const Floors &fl)
explicit Bounds(Floors &fl)
: floors_(fl), ceilings_(Ceilings()), mhd_ceilings_(MHDCeilings()),
radiation_floors_(RadiationFloors()), radiation_ceilings_(RadiationCeilings()) {}
explicit Bounds(const Ceilings &cl)
explicit Bounds(Ceilings &cl)
: floors_(Floors()), ceilings_(cl), mhd_ceilings_(MHDCeilings()),
radiation_floors_(RadiationFloors()), radiation_ceilings_(RadiationCeilings()) {}

template <class... Args>
KOKKOS_INLINE_FUNCTION void SetEOSBnds(Args &&...args) {
floors_.SetEOSBnds(std::forward<Args>(args)...);
ceilings_.SetEOSBnds(std::forward<Args>(args)...);
}
AstroBarker marked this conversation as resolved.
Show resolved Hide resolved

template <class... Args>
KOKKOS_INLINE_FUNCTION void GetFloors(Args &&...args) const {
floors_.GetFloors(std::forward<Args>(args)...);
Expand All @@ -290,8 +330,8 @@ class Bounds {
}

private:
const Floors floors_;
const Ceilings ceilings_;
Floors floors_;
Ceilings ceilings_;
const MHDCeilings mhd_ceilings_;
const RadiationFloors radiation_floors_;
const RadiationCeilings radiation_ceilings_;
Expand Down
6 changes: 3 additions & 3 deletions src/fixup/fixup_c2p.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) {

auto eos = eos_pkg->Param<Microphysics::EOS::EOS>("d.EOS");
auto geom = Geometry::GetCoordinateSystem(rc);
auto bounds = fix_pkg->Param<Bounds>("bounds");
Bounds *bounds = fix_pkg->MutableParam<Bounds>("bounds");

Coordinates_t coords = rc->GetParentPointer()->coords;

Expand All @@ -176,8 +176,8 @@ TaskStatus ConservedToPrimitiveFixupImpl(T *rc) {
eos_lambda[1] = std::log10(v(b, tmp, k, j, i));

Real gamma_max, e_max;
bounds.GetCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), gamma_max, e_max);
bounds->GetCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), gamma_max, e_max);

if (c2p_failure_force_fixup_both && rad_active) {
if (v(b, ifail, k, j, i) == con2prim_robust::FailFlags::fail ||
Expand Down
6 changes: 3 additions & 3 deletions src/fixup/fixup_radc2p.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ TaskStatus RadConservedToPrimitiveFixupImpl(T *rc) {
}

auto geom = Geometry::GetCoordinateSystem(rc);
auto bounds = fix_pkg->Param<Bounds>("bounds");
Bounds *bounds = fix_pkg->MutableParam<Bounds>("bounds");

Coordinates_t coords = rc->GetParentPointer()->coords;

Expand All @@ -129,8 +129,8 @@ TaskStatus RadConservedToPrimitiveFixupImpl(T *rc) {
KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
Real xi_max;
Real garbage;
bounds.GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), xi_max, garbage);
bounds->GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), xi_max, garbage);

// It is assumed that the fluid is already fixed up
auto fail = [&](const int k, const int j, const int i) {
Expand Down
10 changes: 5 additions & 5 deletions src/fixup/fixup_src.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ TaskStatus SourceFixupImpl(T *rc) {
}

auto eos = eos_pkg->Param<EOS>("d.EOS");
auto bounds = fix_pkg->Param<Bounds>("bounds");
Bounds *bounds = fix_pkg->MutableParam<Bounds>("bounds");

const std::vector<std::string> vars(
{p::density::name(), c::density::name(), p::velocity::name(), c::momentum::name(),
Expand Down Expand Up @@ -151,12 +151,12 @@ TaskStatus SourceFixupImpl(T *rc) {
kb.e, jb.s, jb.e, ib.s, ib.e,
KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
double gamma_max, e_max;
bounds.GetCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), gamma_max, e_max);
bounds->GetCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), gamma_max, e_max);
Real xi_max;
Real garbage;
bounds.GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), xi_max, garbage);
bounds->GetRadiationCeilings(coords.Xc<1>(k, j, i), coords.Xc<2>(k, j, i),
coords.Xc<3>(k, j, i), xi_max, garbage);

double eos_lambda[2]; // used for stellarcollapse eos and other EOS's that require
// root finding.
Expand Down
20 changes: 10 additions & 10 deletions src/fluid/con2prim_robust.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,15 +48,15 @@ class Residual {
KOKKOS_FUNCTION
Residual(const Real D, const Real q, const Real bsq, const Real bsq_rpsq,
const Real rsq, const Real rbsq, const Real v0sq, const Real Ye,
const Microphysics::EOS::EOS &eos, const fixup::Bounds &bnds, const Real x1,
const Microphysics::EOS::EOS &eos, const fixup::Bounds *bnds, const Real x1,
const Real x2, const Real x3, const Real floor_scale_fac)
: D_(D), q_(q), bsq_(bsq), bsq_rpsq_(bsq_rpsq), rsq_(rsq), rbsq_(rbsq), v0sq_(v0sq),
eos_(eos), bounds_(bnds), x1_(x1), x2_(x2), x3_(x3),
floor_scale_fac_(floor_scale_fac) {
lambda_[0] = Ye;
Real garbage = 0.0;
bounds_.GetFloors(x1_, x2_, x3_, rho_floor_, garbage);
bounds_.GetCeilings(x1_, x2_, x3_, gam_max_, e_max_);
bounds_->GetFloors(x1_, x2_, x3_, rho_floor_, garbage);
bounds_->GetCeilings(x1_, x2_, x3_, gam_max_, e_max_);

rho_floor_ *= floor_scale_fac_;
}
Expand Down Expand Up @@ -100,7 +100,7 @@ class Residual {
Real ehat_mu(const Real mu, const Real qbar, const Real rbarsq, const Real vhatsq,
const Real What) {
Real rho = rhohat_mu(1.0 / What);
bounds_.GetFloors(x1_, x2_, x3_, rho, e_floor_);
bounds_->GetFloors(x1_, x2_, x3_, rho, e_floor_);
e_floor_ *= floor_scale_fac_;
const Real ehat_trial =
What * (qbar - mu * rbarsq) + vhatsq * What * What / (1.0 + What);
Expand Down Expand Up @@ -164,7 +164,7 @@ class Residual {
private:
const Real D_, q_, bsq_, bsq_rpsq_, rsq_, rbsq_, v0sq_;
const Microphysics::EOS::EOS &eos_;
const fixup::Bounds &bounds_;
const fixup::Bounds *bounds_;
const Real x1_, x2_, x3_;
const Real floor_scale_fac_;
Real lambda_[2];
Expand Down Expand Up @@ -224,7 +224,7 @@ struct CellGeom {
template <typename Data_t, typename T>
class ConToPrim {
public:
ConToPrim(Data_t *rc, fixup::Bounds bnds, const Real tol, const int max_iterations,
ConToPrim(Data_t *rc, fixup::Bounds *bnds, const Real tol, const int max_iterations,
const Real floor_scale_fac, const bool fail_on_floors,
const bool fail_on_ceilings)
: bounds(bnds), var(rc->PackVariables(Vars(), imap)),
Expand Down Expand Up @@ -278,7 +278,7 @@ class ConToPrim {
int NumBlocks() { return var.GetDim(5); }

private:
fixup::Bounds bounds;
fixup::Bounds *bounds;
PackIndexMap imap;
const T var;
const int prho, crho;
Expand Down Expand Up @@ -306,11 +306,11 @@ class ConToPrim {

Real rhoflr = 0.0;
Real epsflr;
bounds.GetFloors(x1, x2, x3, rhoflr, epsflr);
bounds->GetFloors(x1, x2, x3, rhoflr, epsflr);
rhoflr *= floor_scale_fac_;
epsflr *= floor_scale_fac_;
Real gam_max, eps_max;
bounds.GetCeilings(x1, x2, x3, gam_max, eps_max);
bounds->GetCeilings(x1, x2, x3, gam_max, eps_max);

const Real D = v(crho) * igdet;
#if USE_VALENCIA
Expand Down Expand Up @@ -456,7 +456,7 @@ class ConToPrim {
using C2P_Block_t = ConToPrim<MeshBlockData<Real>, VariablePack<Real>>;
using C2P_Mesh_t = ConToPrim<MeshData<Real>, MeshBlockPack<Real>>;

inline C2P_Block_t ConToPrimSetup(MeshBlockData<Real> *rc, fixup::Bounds bounds,
inline C2P_Block_t ConToPrimSetup(MeshBlockData<Real> *rc, fixup::Bounds *bounds,
const Real tol, const int max_iter,
const Real c2p_floor_scale_fac,
const bool fail_on_floors,
Expand Down
Loading
Loading