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

WIP: Number Moments #168

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open

WIP: Number Moments #168

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
98594b9
Update to newest Parthenon
lroberts36 Apr 19, 2023
dc7c40b
Part way through adding number evolution options
lroberts36 Apr 20, 2023
774b7cf
Update all explicit tasks to include number evolution
lroberts36 Apr 20, 2023
8f0d6b3
Move radiation variables
lroberts36 Apr 20, 2023
99c48e2
Do fixup for number as well
lroberts36 Apr 20, 2023
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
2 changes: 1 addition & 1 deletion external/parthenon
Open in desktop
Submodule parthenon updated from eb16b5 to 01cd49
40 changes: 22 additions & 18 deletions src/fixup/fixup_radc2p.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,14 +41,11 @@ using robust::ratio;

namespace fixup {

template <typename T, class CLOSURE>
template <typename T, class CLOSURE, radiation::RadEnergyMoment EMOMENT>
TaskStatus RadConservedToPrimitiveFixupImpl(T *rc) {
namespace p = fluid_prim;
namespace c = fluid_cons;
namespace impl = internal_variables;
namespace ir = radmoment_internal;
namespace pr = radmoment_prim;
namespace cr = radmoment_cons;

auto *pmb = rc->GetParentPointer().get();
IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
Expand All @@ -62,23 +59,24 @@ TaskStatus RadConservedToPrimitiveFixupImpl(T *rc) {
bool enable_c2p_fixup = fix_pkg->Param<bool>("enable_c2p_fixup");
bool update_rad = rad_pkg->Param<bool>("active");
if (!enable_c2p_fixup || !update_rad) return TaskStatus::complete;

const std::vector<std::string> vars({p::velocity, p::ye, c::ye, pr::J, pr::H, cr::E,
cr::F, ir::tilPi, ir::c2pfail, impl::fail});

auto var = radiation::RadiationVariableNames::GetByMomentType(EMOMENT);
const std::vector<std::string> vars({p::velocity, p::ye, c::ye, var.J, var.H, var.E,
var.F, var.tilPi, var.c2pfail, impl::fail});

PackIndexMap imap;
auto v = rc->PackVariables(vars, imap);

auto idx_pvel = imap.GetFlatIdx(p::velocity);
int pye = imap[p::ye].second; // negative if not present
int cye = imap[c::ye].second;
auto idx_J = imap.GetFlatIdx(pr::J, false);
auto idx_H = imap.GetFlatIdx(pr::H, false);
auto idx_E = imap.GetFlatIdx(cr::E, false);
auto idx_F = imap.GetFlatIdx(cr::F, false);
auto idx_J = imap.GetFlatIdx(var.J, false);
auto idx_H = imap.GetFlatIdx(var.H, false);
auto idx_E = imap.GetFlatIdx(var.E, false);
auto idx_F = imap.GetFlatIdx(var.F, false);
int ifluidfail = imap[impl::fail].first;
int iradfail = imap[ir::c2pfail].first;
auto iTilPi = imap.GetFlatIdx(ir::tilPi, false);
int iradfail = imap[var.c2pfail].first;
auto iTilPi = imap.GetFlatIdx(var.tilPi, false);

bool report_c2p_fails = fix_pkg->Param<bool>("report_c2p_fails");
if (report_c2p_fails) {
Expand Down Expand Up @@ -227,22 +225,28 @@ TaskStatus RadConservedToPrimitiveFixup(T *rc) {
StateDescriptor *rad_pkg = pm->packages.Get("radiation").get();
StateDescriptor *fix_pkg = pm->packages.Get("fixup").get();
const bool enable_rad_floors = fix_pkg->Param<bool>("enable_rad_floors");
auto do_num = rad_pkg->Param<bool>("do_number_evolution");
std::string method;
if (enable_rad_floors) {
method = rad_pkg->Param<std::string>("method");
} else {
return TaskStatus::complete;
}

using namespace radiation;
// TODO(BRR) share these settings somewhere else. Set at configure time?
using settings =
using settings_N =
ClosureSettings<ClosureEquation::number_conserve, ClosureVerbosity::quiet>;
using settings_E =
ClosureSettings<ClosureEquation::energy_conserve, ClosureVerbosity::quiet>;
if (method == "moment_m1") {
return RadConservedToPrimitiveFixupImpl<T, radiation::ClosureM1<settings>>(rc);
if (do_num) RadConservedToPrimitiveFixupImpl<T, ClosureM1<settings_N>, RadEnergyMoment::Number>(rc);
return RadConservedToPrimitiveFixupImpl<T, ClosureM1<settings_E>, RadEnergyMoment::Energy>(rc);
} else if (method == "moment_eddington") {
return RadConservedToPrimitiveFixupImpl<T, radiation::ClosureEdd<settings>>(rc);
if (do_num) RadConservedToPrimitiveFixupImpl<T, ClosureEdd<settings_N>, RadEnergyMoment::Number>(rc);
return RadConservedToPrimitiveFixupImpl<T, ClosureEdd<settings_E>, RadEnergyMoment::Energy>(rc);
} else if (method == "mocmc") {
return RadConservedToPrimitiveFixupImpl<T, radiation::ClosureMOCMC<settings>>(rc);
if (do_num) RadConservedToPrimitiveFixupImpl<T, ClosureMOCMC<settings_N>, RadEnergyMoment::Number>(rc);
return RadConservedToPrimitiveFixupImpl<T, ClosureMOCMC<settings_E>, RadEnergyMoment::Energy>(rc);
}
return TaskStatus::fail;
}
Expand Down
28 changes: 21 additions & 7 deletions src/phoebus_driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ namespace phoebus {
PhoebusDriver::PhoebusDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm,
const bool is_restart)
: EvolutionDriver(pin, app_in, pm),
integrator(std::make_unique<StagedIntegrator>(pin)), is_restart_(is_restart) {
integrator(std::make_unique<parthenon::LowStorageIntegrator>(pin)), is_restart_(is_restart) {

// fail if these are not specified in the input file
pin->CheckRequired("parthenon/mesh", "ix1_bc");
Expand Down Expand Up @@ -192,7 +192,12 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) {
bool rad_mocmc_active = false;
if (rad_active) {
rad_mocmc_active = (rad->Param<std::string>("method") == "mocmc");
}
bool rad_moment_number_evolve = false;
if (rad_moments_active) {
rad_moment_number_evolve = rad->Param<bool>("do_number_evolution");
}

const auto monopole_enabled = monopole->Param<bool>("enable_monopole_gr");
// Force static here means monopole only called at initialization.
// and source terms are disabled
Expand Down Expand Up @@ -312,13 +317,22 @@ TaskCollection PhoebusDriver::RungeKuttaStage(const int stage) {

if (rad_moments_active) {
using MDT = std::remove_pointer<decltype(sc0.get())>::type;
auto moment_recon =
tl.AddTask(none, radiation::ReconstructEdgeStates<MDT>, sc0.get());
// TODO(BRR) Remove from task list if not using due to MOCMC
auto get_opacities =
tl.AddTask(moment_recon, radiation::MomentCalculateOpacities<MDT>, sc0.get());
auto moment_flux =
tl.AddTask(get_opacities, radiation::CalculateFluxes<MDT>, sc0.get());
tl.AddTask(none, radiation::MomentCalculateOpacities<MDT>, sc0.get());

auto moment_recon_E =
tl.AddTask(get_opacities, radiation::ReconstructEdgeStates<MDT>, sc0.get());
// TODO(BRR) Remove from task list if not using due to MOCMC
auto moment_flux_E =
tl.AddTask(moment_recon_E, radiation::CalculateFluxes<MDT>, sc0.get());
auto moment_recon = moment_recon_E;
auto moment_flux = moment_flux_E;
if (rad_moment_number_evolve) {
moment_recon =
tl.AddTask(moment_flux_E, radiation::ReconstructEdgeStates<MDT, radiation::RadEnergyMoment::Number>, sc0.get());
moment_flux =
tl.AddTask(moment_recon, radiation::CalculateFluxes<MDT, radiation::RadEnergyMoment::Number>, sc0.get());
}
auto moment_geom_src = tl.AddTask(none, radiation::CalculateGeometricSource<MDT>,
sc0.get(), gsrc.get());
sndrcv_flux_depend = sndrcv_flux_depend | moment_flux;
Expand Down
2 changes: 1 addition & 1 deletion src/phoebus_driver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ class PhoebusDriver : public EvolutionDriver {
TaskListStatus Step();

private:
std::unique_ptr<StagedIntegrator> integrator;
std::unique_ptr<parthenon::LowStorageIntegrator> integrator;
const bool is_restart_;
Real dt_init, dt_init_fact;

Expand Down
18 changes: 18 additions & 0 deletions src/phoebus_utils/variables.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,11 +36,19 @@ constexpr char ye[] = "c.ye";
namespace radmoment_prim {
constexpr char J[] = "r.p.J";
constexpr char H[] = "r.p.H";
namespace num {
constexpr char J[] = "r.p.n.J";
constexpr char H[] = "r.p.n.H";
} // namespace num
} // namespace radmoment_prim

namespace radmoment_cons {
constexpr char E[] = "r.c.E";
constexpr char F[] = "r.c.F";
namespace num {
constexpr char E[] = "r.c.n.E";
constexpr char F[] = "r.c.n.F";
} // namespace radmoment_cons
} // namespace radmoment_cons

namespace radmoment_internal {
Expand All @@ -58,6 +66,16 @@ constexpr char tilPi[] = "r.i.tilPi";
constexpr char kappaH_mean[] = "r.i.kappaH_mean";
constexpr char c2pfail[] = "r.i.c2p_fail";
constexpr char srcfail[] = "r.i.src_fail";
namespace num {
constexpr char xi[] = "r.i.n.xi";
constexpr char phi[] = "r.i.n.phi";
constexpr char dJ[] = "r.i.n.dJ";
constexpr char kappaJ[] = "r.i.n.kappaJ";
constexpr char kappaH[] = "r.i.n.kappaH";
constexpr char kappaH_mean[] = "r.i.n.kappaH_mean";
constexpr char JBB[] = "r.i.n.JBB";
constexpr char tilPi[] = "r.i.n.tilPi";
}
} // namespace radmoment_internal

namespace mocmc_internal {
Expand Down
6 changes: 3 additions & 3 deletions src/radiation/closure.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,8 @@ enum class ClosureVerbosity { quiet, v1, v2 };
template <ClosureEquation EQ = ClosureEquation::energy_conserve,
ClosureVerbosity VB = ClosureVerbosity::quiet>
struct ClosureSettings {
static const ClosureEquation eqn_type = EQ;
static const ClosureVerbosity verbosity = VB;
static constexpr ClosureEquation eqn_type = EQ;
static constexpr ClosureVerbosity verbosity = VB;
};

enum class ClosureCon2PrimStrategy {
Expand All @@ -66,7 +66,7 @@ template <class SET = ClosureSettings<>>
class ClosureEdd {
public:
using LocalGeometryType = LocalThreeGeometry;

static constexpr ClosureEquation eqn_type = SET::eqn_type;
//-------------------------------------------------------------------------------------
/// Constructor just calculates the inverse 3-metric, covariant three-velocity, and the
/// Lorentz factor for the given background state.
Expand Down
Loading