From 5809482c2bad967025fbc15043fde6c93f231cba Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Wed, 18 Jun 2025 21:58:52 -0600 Subject: [PATCH 01/10] Precompute gray opacity per cell. --- CMakeLists.txt | 2 ++ src/artemis.hpp | 4 ++++ src/derived/fill_derived.cpp | 24 ++++++++++++++++++++++++ src/gas/gas.cpp | 7 +++++++ 4 files changed, 37 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index d00b1513..64174856 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -228,6 +228,8 @@ set(JAYBENNE_HOST_VARIABLE_HEADER "src/artemis.hpp") set(JAYBENNE_HOST_DENSITY_VARIABLE "gas::prim::density") set(JAYBENNE_HOST_SPECIFIC_INTERNAL_ENERGY_VARIABLE "gas::prim::sie") set(JAYBENNE_HOST_UPDATE_ENERGY_VARIABLE "gas::cons::internal_energy") +set(JAYBENNE_HOST_ABSORPTION_OPACITY_VARIABLE "gas::opac::absorption") +set(JAYBENNE_HOST_SCATTERING_OPACITY_VARIABLE "gas::opac::scattering") # Add jaybenne message(STATUS "Adding jaybenne and dependencies") diff --git a/src/artemis.hpp b/src/artemis.hpp index 384d2099..0033fa0b 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -56,6 +56,10 @@ ARTEMIS_VARIABLE(gas.diff, energy); namespace face { ARTEMIS_VARIABLE(gas.face, velocity); } // namespace face +namespace opac { +ARTEMIS_VARIABLE(gas.opac, absorption); +ARTEMIS_VARIABLE(gas.opac, scattering); +} // namespace opac } // namespace gas namespace dust { diff --git a/src/derived/fill_derived.cpp b/src/derived/fill_derived.cpp index a670c8a0..890814ac 100644 --- a/src/derived/fill_derived.cpp +++ b/src/derived/fill_derived.cpp @@ -18,9 +18,12 @@ #include "radiation/moments/moments.hpp" #include "utils/artemis_utils.hpp" #include "utils/eos/eos.hpp" +#include "utils/opacity/opacity.hpp" using ArtemisUtils::EOS; using ArtemisUtils::VI; +using ArtemisUtils::MeanOpacity; +using ArtemisUtils::MeanScattering; namespace ArtemisDerived { //---------------------------------------------------------------------------------------- @@ -74,6 +77,9 @@ TaskStatus SetAuxillaryFields(MeshData *md) { u_u = (ufloor)*utmp + (!ufloor) * uflr; } }); + + + return TaskStatus::complete; } @@ -220,16 +226,24 @@ void PrimToCons(T *md) { const bool do_gas = artemis_pkg->template Param("do_gas"); const bool do_dust = artemis_pkg->template Param("do_dust"); const bool do_rad = artemis_pkg->template Param("do_moment"); + const bool do_imc = artemis_pkg->template Param("do_imc"); // Extract gas parameters Real dflr_gas = Null(); Real sieflr_gas = Null(); EOS eos_d; + MeanOpacity mopacity_d; + MeanScattering mscattering_d; if (do_gas) { auto &gas_pkg = pm->packages.Get("gas"); dflr_gas = gas_pkg->template Param("dfloor"); sieflr_gas = gas_pkg->template Param("siefloor"); eos_d = gas_pkg->template Param("eos_d"); + // opacity types + if (do_imc) { + mopacity_d = gas_pkg->template Param("mopacity_d"); + mscattering_d = gas_pkg->template Param("mscattering_d"); + } } // Extract dust parameters @@ -252,6 +266,7 @@ void PrimToCons(T *md) { MakePackDescriptor( @@ -304,6 +319,15 @@ void PrimToCons(T *md) { const Real ke = 0.5 * w_d * (SQR(vel1) + SQR(vel2) + SQR(vel3)); Real &u_e = vmesh(b, gas::cons::total_energy(n), k, j, i); u_e = u_u + ke; + + // Sync opacity (TODO: do_rad) + if (do_imc) { + Real &aa = vmesh(b, gas::opac::absorption(), k, j, i); + Real &ss = vmesh(b, gas::opac::scattering(), k, j, i); + const Real temp = eos_d.TemperatureFromDensityInternalEnergy(w_d, w_s); + aa = mopacity_d.AbsorptionCoefficient(w_d, temp); + ss = mscattering_d.RosselandMeanTotalScatteringCoefficient(w_d, temp); + } } } diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index 3aa86f68..511fb947 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -329,6 +329,13 @@ std::shared_ptr Initialize(ParameterInput *pin, m.SetSparseThresholds(0.0, 0.0, 0.0); gas->AddSparsePool(m, control_field, fluidids); + // Absorption and scattering opacity + m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}); + ArtemisUtils::EnrollArtemisRefinementOps(m, coords); + m.SetSparseThresholds(0.0, 0.0, 0.0); + gas->AddSparsePool(m, control_field, fluidids); + gas->AddSparsePool(m, control_field, fluidids); + // Normal face Velocity for PdV evaluation of internal energy m = Metadata({Metadata::Face, Metadata::Derived, Metadata::OneCopy, Metadata::Sparse}); m.SetSparseThresholds(0.0, 0.0, 0.0); From 242c54023a007a44953bd362995031cf10d809b4 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Thu, 19 Jun 2025 14:59:38 -0600 Subject: [PATCH 02/10] Move opacity precomputation to dedicated task (collection). + Add radiation.hpp and radiation.cpp for task (collection). + Add InitializeRadDataFields and UpdateRadDataFields tasks here. + Call these tasts only for IMC for now. + Update jaybenne submodule to commit with host opacity fields. --- CMakeLists.txt | 4 +- external/jaybenne | 2 +- src/CMakeLists.txt | 2 + src/artemis.cpp | 2 + src/artemis.hpp | 8 +-- src/derived/fill_derived.cpp | 26 -------- src/gas/gas.cpp | 7 -- src/radiation/imc/imc_driver.cpp | 4 +- src/radiation/radiation.cpp | 106 +++++++++++++++++++++++++++++++ src/radiation/radiation.hpp | 18 ++++++ 10 files changed, 138 insertions(+), 41 deletions(-) create mode 100644 src/radiation/radiation.cpp create mode 100644 src/radiation/radiation.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 64174856..410f5e84 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -228,8 +228,8 @@ set(JAYBENNE_HOST_VARIABLE_HEADER "src/artemis.hpp") set(JAYBENNE_HOST_DENSITY_VARIABLE "gas::prim::density") set(JAYBENNE_HOST_SPECIFIC_INTERNAL_ENERGY_VARIABLE "gas::prim::sie") set(JAYBENNE_HOST_UPDATE_ENERGY_VARIABLE "gas::cons::internal_energy") -set(JAYBENNE_HOST_ABSORPTION_OPACITY_VARIABLE "gas::opac::absorption") -set(JAYBENNE_HOST_SCATTERING_OPACITY_VARIABLE "gas::opac::scattering") +set(JAYBENNE_HOST_ABSORPTION_OPACITY_VARIABLE "rad::opac::absorption") +set(JAYBENNE_HOST_SCATTERING_OPACITY_VARIABLE "rad::opac::scattering") # Add jaybenne message(STATUS "Adding jaybenne and dependencies") diff --git a/external/jaybenne b/external/jaybenne index 14d1ab6b..6b19009f 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit 14d1ab6bc4d476ab4bd3abd009eaff68eff94b73 +Subproject commit 6b19009f8878b2843b31feb2b066e27065acdc3d diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 97a75395..a661fa15 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -71,6 +71,8 @@ set (SRC_LIST pgen/strat.hpp pgen/thermalization.hpp + radiation/radiation.hpp + radiation/radiation.cpp radiation/imc/imc_driver.cpp radiation/imc/imc.hpp radiation/moments/moments.cpp diff --git a/src/artemis.cpp b/src/artemis.cpp index 175e4cb2..5f3b9a7d 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -22,6 +22,7 @@ #include "gravity/gravity.hpp" #include "nbody/nbody.hpp" #include "radiation/moments/moments.hpp" +#include "radiation/radiation.hpp" #include "rotating_frame/rotating_frame.hpp" #include "utils/artemis_utils.hpp" #include "utils/history.hpp" @@ -132,6 +133,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { auto eos_h = packages.Get("gas")->Param("eos_h"); auto opacity_h = packages.Get("gas")->Param("opacity_h"); auto scattering_h = packages.Get("gas")->Param("scattering_h"); + packages.Add(rad::InitializeRadDataFields(pin.get())); packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h, "radiation/imc")); PARTHENON_REQUIRE(coords == Coordinates::cartesian, diff --git a/src/artemis.hpp b/src/artemis.hpp index 0033fa0b..df0105d6 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -56,10 +56,6 @@ ARTEMIS_VARIABLE(gas.diff, energy); namespace face { ARTEMIS_VARIABLE(gas.face, velocity); } // namespace face -namespace opac { -ARTEMIS_VARIABLE(gas.opac, absorption); -ARTEMIS_VARIABLE(gas.opac, scattering); -} // namespace opac } // namespace gas namespace dust { @@ -83,6 +79,10 @@ ARTEMIS_VARIABLE(rad.prim, energy); ARTEMIS_VARIABLE(rad.prim, pressure); ARTEMIS_VARIABLE(rad.prim, flux); } // namespace prim +namespace opac { +ARTEMIS_VARIABLE(rad.opac, absorption); +ARTEMIS_VARIABLE(rad.opac, scattering); +} // namespace opac } // namespace rad #undef ARTEMIS_VARIABLE diff --git a/src/derived/fill_derived.cpp b/src/derived/fill_derived.cpp index 890814ac..937b91fc 100644 --- a/src/derived/fill_derived.cpp +++ b/src/derived/fill_derived.cpp @@ -18,12 +18,9 @@ #include "radiation/moments/moments.hpp" #include "utils/artemis_utils.hpp" #include "utils/eos/eos.hpp" -#include "utils/opacity/opacity.hpp" using ArtemisUtils::EOS; using ArtemisUtils::VI; -using ArtemisUtils::MeanOpacity; -using ArtemisUtils::MeanScattering; namespace ArtemisDerived { //---------------------------------------------------------------------------------------- @@ -33,7 +30,6 @@ namespace ArtemisDerived { template TaskStatus SetAuxillaryFields(MeshData *md) { using parthenon::MakePackDescriptor; - using TE = parthenon::TopologicalElement; auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; @@ -77,9 +73,6 @@ TaskStatus SetAuxillaryFields(MeshData *md) { u_u = (ufloor)*utmp + (!ufloor) * uflr; } }); - - - return TaskStatus::complete; } @@ -90,7 +83,6 @@ TaskStatus SetAuxillaryFields(MeshData *md) { template void ConsToPrim(MeshData *md) { using parthenon::MakePackDescriptor; - using TE = parthenon::TopologicalElement; auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; @@ -217,7 +209,6 @@ void ConsToPrim(MeshData *md) { template void PrimToCons(T *md) { using parthenon::MakePackDescriptor; - using TE = parthenon::TopologicalElement; auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; @@ -232,18 +223,11 @@ void PrimToCons(T *md) { Real dflr_gas = Null(); Real sieflr_gas = Null(); EOS eos_d; - MeanOpacity mopacity_d; - MeanScattering mscattering_d; if (do_gas) { auto &gas_pkg = pm->packages.Get("gas"); dflr_gas = gas_pkg->template Param("dfloor"); sieflr_gas = gas_pkg->template Param("siefloor"); eos_d = gas_pkg->template Param("eos_d"); - // opacity types - if (do_imc) { - mopacity_d = gas_pkg->template Param("mopacity_d"); - mscattering_d = gas_pkg->template Param("mscattering_d"); - } } // Extract dust parameters @@ -266,7 +250,6 @@ void PrimToCons(T *md) { MakePackDescriptor( @@ -319,15 +302,6 @@ void PrimToCons(T *md) { const Real ke = 0.5 * w_d * (SQR(vel1) + SQR(vel2) + SQR(vel3)); Real &u_e = vmesh(b, gas::cons::total_energy(n), k, j, i); u_e = u_u + ke; - - // Sync opacity (TODO: do_rad) - if (do_imc) { - Real &aa = vmesh(b, gas::opac::absorption(), k, j, i); - Real &ss = vmesh(b, gas::opac::scattering(), k, j, i); - const Real temp = eos_d.TemperatureFromDensityInternalEnergy(w_d, w_s); - aa = mopacity_d.AbsorptionCoefficient(w_d, temp); - ss = mscattering_d.RosselandMeanTotalScatteringCoefficient(w_d, temp); - } } } diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index 511fb947..3aa86f68 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -329,13 +329,6 @@ std::shared_ptr Initialize(ParameterInput *pin, m.SetSparseThresholds(0.0, 0.0, 0.0); gas->AddSparsePool(m, control_field, fluidids); - // Absorption and scattering opacity - m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}); - ArtemisUtils::EnrollArtemisRefinementOps(m, coords); - m.SetSparseThresholds(0.0, 0.0, 0.0); - gas->AddSparsePool(m, control_field, fluidids); - gas->AddSparsePool(m, control_field, fluidids); - // Normal face Velocity for PdV evaluation of internal energy m = Metadata({Metadata::Face, Metadata::Derived, Metadata::OneCopy, Metadata::Sparse}); m.SetSparseThresholds(0.0, 0.0, 0.0); diff --git a/src/radiation/imc/imc_driver.cpp b/src/radiation/imc/imc_driver.cpp index 3388baaa..f680ccaa 100644 --- a/src/radiation/imc/imc_driver.cpp +++ b/src/radiation/imc/imc_driver.cpp @@ -15,6 +15,7 @@ #include "artemis.hpp" #include "derived/fill_derived.hpp" #include "radiation/imc/imc.hpp" +#include "radiation/radiation.hpp" // Jaybenne includes #include "jaybenne.hpp" @@ -28,7 +29,8 @@ namespace IMC { //! \brief Executes thermal IMC transport (Jaybenne) and syncs updated fields template TaskListStatus JaybenneIMC(Mesh *pmesh, const Real time, const Real dt) { - auto status = jaybenne::RadiationStep(pmesh, time, dt).Execute(); + auto status = rad::UpdateRadDataFields(pmesh).Execute(); + status = jaybenne::RadiationStep(pmesh, time, dt).Execute(); if (status != TaskListStatus::complete) return status; status = ArtemisDerived::SyncFields(pmesh, time, dt).Execute(); return status; diff --git a/src/radiation/radiation.cpp b/src/radiation/radiation.cpp new file mode 100644 index 00000000..791dc759 --- /dev/null +++ b/src/radiation/radiation.cpp @@ -0,0 +1,106 @@ +//======================================================================================== +// (C) (or copyright) 2025. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +// Artemis includes +#include "artemis.hpp" +#include "geometry/geometry.hpp" +#include "utils/artemis_utils.hpp" +#include "utils/eos/eos.hpp" +#include "utils/opacity/opacity.hpp" + +using ArtemisUtils::EOS; +using ArtemisUtils::MeanOpacity; +using ArtemisUtils::MeanScattering; + +namespace rad { + +std::shared_ptr InitializeRadDataFields(ParameterInput *pin) { + + // instantiate a state descriptor for the fields + auto rad_fs = std::make_shared("rad_fields"); + + // Only one photon species + std::vector radids = {0}; + + // Control field for sparse gas fields + const std::string control_field = "rad_fields"; + + // const int ndim = ProblemDimension(pin); + // std::string sys = pin->GetOrAddString("artemis", "coordinates", "cartesian"); + // Coordinates coords = geometry::CoordSelect(sys, ndim); + + // Absorption and scattering opacity + Metadata m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy, Metadata::Sparse}); + //ArtemisUtils::EnrollArtemisRefinementOps(m, coords); + //m.SetSparseThresholds(0.0, 0.0, 0.0); + rad_fs->AddSparsePool(m, control_field, radids); + rad_fs->AddSparsePool(m, control_field, radids); + + // return rad_fields state descriptor + return rad_fs; +} + +TaskStatus UpRadDataFields(MeshData *md) { + using parthenon::MakePackDescriptor; + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + auto &gas_pkg = pm->packages.Get("gas"); + + EOS eos_d = gas_pkg->template Param("eos_d"); + MeanOpacity opacity_d = gas_pkg->template Param("opacity_d"); + MeanScattering scattering_d = gas_pkg->template Param("scattering_d"); + + // Packing and indexing (TODO: use dust sie, density) + static auto desc = MakePackDescriptor< + gas::prim::density, gas::prim::sie, + rad::opac::absorption, rad::opac::scattering>(resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + const int nblocks = md->NumBlocks(); + IndexRange ib = md->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBoundsK(IndexDomain::interior); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "ConsToPrim", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + + const Real &rho = vmesh(b, gas::prim::density(), k, j, i); + const Real &sie = vmesh(b, gas::prim::sie(), k, j, i); + const Real temp = eos_d.TemperatureFromDensityInternalEnergy(rho, sie); + Real &aa = vmesh(b, rad::opac::absorption(), k, j, i); + Real &ss = vmesh(b, rad::opac::scattering(), k, j, i); + + aa = opacity_d.AbsorptionCoefficient(rho, temp); + ss = scattering_d.RosselandMeanTotalScatteringCoefficient(rho, temp); + }); + + return TaskStatus::complete; +} + +// task collection for updating data fields +TaskCollection UpdateRadDataFields(Mesh *pmesh) { + TaskCollection tc; + TaskID none(0); + const int num_partitions = pmesh->DefaultNumPartitions(); + auto ® = tc.AddRegion(num_partitions); + for (int i = 0; i < num_partitions; i++) { + auto &tl = reg[i]; + auto &base = pmesh->mesh_data.GetOrAdd("base", i); + auto upradf = tl.AddTask(none, UpRadDataFields, base.get()); + } + + return tc; +} + +} // namespace rad diff --git a/src/radiation/radiation.hpp b/src/radiation/radiation.hpp new file mode 100644 index 00000000..e275b2ec --- /dev/null +++ b/src/radiation/radiation.hpp @@ -0,0 +1,18 @@ +//======================================================================================== +// (C) (or copyright) 2025. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== + +namespace rad { +std::shared_ptr InitializeRadDataFields(ParameterInput *pin); +TaskStatus UpRadDataFields(MeshData *md); +TaskCollection UpdateRadDataFields(Mesh *pmesh); +} // namespace rad From a4cec016e6f3cd763163d9fb167de429eb1af25e Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Fri, 20 Jun 2025 10:07:30 -0400 Subject: [PATCH 03/10] @pdmullen pass at A--J API change --- src/artemis.cpp | 7 +- src/artemis_driver.cpp | 2 +- src/derived/fill_derived.cpp | 4 +- src/radiation/imc/imc_driver.cpp | 3 +- src/radiation/moments/matter_coupling.hpp | 48 +++++----- src/radiation/moments/moments.cpp | 75 ++++++++-------- src/radiation/moments/moments.hpp | 24 ++--- src/radiation/moments/moments_driver.cpp | 13 ++- src/radiation/params.yaml | 5 -- src/radiation/radiation.cpp | 104 ++++++++++++++-------- src/radiation/radiation.hpp | 21 +++-- src/utils/fluxes/fluid_fluxes.hpp | 2 +- src/utils/fluxes/riemann/hlle.hpp | 8 +- src/utils/fluxes/riemann/llf.hpp | 8 +- 14 files changed, 179 insertions(+), 145 deletions(-) diff --git a/src/artemis.cpp b/src/artemis.cpp index 5f3b9a7d..31ceb4cc 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -128,18 +128,19 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { if (do_cooling) packages.Add(Gas::Cooling::Initialize(pin.get())); if (do_drag) packages.Add(Drag::Initialize(pin.get())); if (do_radiation) { - // swap between native artemis radiation and jaybenne imc + // Top-level radiation package + packages.Add(Radiation::Initialize(pin.get(), constants, do_imc)); + // Select between Jaybenne IMC or Moments if (do_imc) { auto eos_h = packages.Get("gas")->Param("eos_h"); auto opacity_h = packages.Get("gas")->Param("opacity_h"); auto scattering_h = packages.Get("gas")->Param("scattering_h"); - packages.Add(rad::InitializeRadDataFields(pin.get())); packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h, "radiation/imc")); PARTHENON_REQUIRE(coords == Coordinates::cartesian, "Jaybenne currently supports only Cartesian coordinates!"); } else if (do_moment) { - packages.Add(Radiation::Initialize(pin.get(), constants)); + packages.Add(Moments::Initialize(pin.get(), constants)); } else { PARTHENON_FAIL("Unknown radiation model!"); } diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index d2d4ed76..ffac5e4b 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -128,7 +128,7 @@ TaskListStatus ArtemisDriver::Step() { if (status != TaskListStatus::complete) return status; // Operator split, moments subcyling (M1 or P1) - if (do_moment) status = Radiation::MomentsDriver(pmesh, tm, rad_integrator.get()); + if (do_moment) status = Moments::MomentsDriver(pmesh, tm, rad_integrator.get()); if (status != TaskListStatus::complete) return status; // Compute new dt, (de)refine, and handle sparse (if enabled) diff --git a/src/derived/fill_derived.cpp b/src/derived/fill_derived.cpp index 937b91fc..945098ec 100644 --- a/src/derived/fill_derived.cpp +++ b/src/derived/fill_derived.cpp @@ -191,7 +191,7 @@ void ConsToPrim(MeshData *md) { const Real hfx1 = vmesh(b, rad::cons::flux(VI(n, 0)), k, j, i) / conv[0]; const Real hfx2 = vmesh(b, rad::cons::flux(VI(n, 1)), k, j, i) / conv[1]; const Real hfx3 = vmesh(b, rad::cons::flux(VI(n, 2)), k, j, i) / conv[2]; - const auto fx = Radiation::NormalizeFlux(hfx1, hfx2, hfx3); + const auto fx = Moments::NormalizeFlux(hfx1, hfx2, hfx3); vmesh(b, rad::cons::flux(VI(n, 0)), k, j, i) = fx[0] * conv[0]; vmesh(b, rad::cons::flux(VI(n, 1)), k, j, i) = fx[1] * conv[1]; vmesh(b, rad::cons::flux(VI(n, 2)), k, j, i) = fx[2] * conv[2]; @@ -342,7 +342,7 @@ void PrimToCons(T *md) { const Real fx1 = vmesh(b, rad::prim::flux(VI(n, 0)), k, j, i); const Real fx2 = vmesh(b, rad::prim::flux(VI(n, 1)), k, j, i); const Real fx3 = vmesh(b, rad::prim::flux(VI(n, 2)), k, j, i); - const auto fx = Radiation::NormalizeFlux(fx1, fx2, fx3); + const auto fx = Moments::NormalizeFlux(fx1, fx2, fx3); vmesh(b, rad::cons::flux(VI(n, 0)), k, j, i) = fx[0] * conv[0]; vmesh(b, rad::cons::flux(VI(n, 1)), k, j, i) = fx[1] * conv[1]; vmesh(b, rad::cons::flux(VI(n, 2)), k, j, i) = fx[2] * conv[2]; diff --git a/src/radiation/imc/imc_driver.cpp b/src/radiation/imc/imc_driver.cpp index f680ccaa..f8bdfdd7 100644 --- a/src/radiation/imc/imc_driver.cpp +++ b/src/radiation/imc/imc_driver.cpp @@ -29,7 +29,8 @@ namespace IMC { //! \brief Executes thermal IMC transport (Jaybenne) and syncs updated fields template TaskListStatus JaybenneIMC(Mesh *pmesh, const Real time, const Real dt) { - auto status = rad::UpdateRadDataFields(pmesh).Execute(); + auto status = Radiation::UpdateRadiationFields(pmesh).Execute(); + if (status != TaskListStatus::complete) return status; status = jaybenne::RadiationStep(pmesh, time, dt).Execute(); if (status != TaskListStatus::complete) return status; status = ArtemisDerived::SyncFields(pmesh, time, dt).Execute(); diff --git a/src/radiation/moments/matter_coupling.hpp b/src/radiation/moments/matter_coupling.hpp index 93d50136..78744937 100644 --- a/src/radiation/moments/matter_coupling.hpp +++ b/src/radiation/moments/matter_coupling.hpp @@ -26,10 +26,10 @@ using ArtemisUtils::MeanOpacity; using ArtemisUtils::MeanScattering; using ArtemisUtils::VI; -namespace Radiation { +namespace Moments { //---------------------------------------------------------------------------------------- -//! \fn TaskStatus Radiation::MatterCouplingSimpleImpl +//! \fn TaskStatus Moments::MatterCouplingSimpleImpl //! \brief Implementation for simple radiation-matter coupling source template TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { @@ -47,15 +47,15 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { auto dflr = gas_pkg->template Param("dfloor"); auto de_switch = gas_pkg->template Param("de_switch"); - // Extract radiation package and params - auto &rad_pkg = pm->packages.Get("moments"); - const auto chat = rad_pkg->template Param("chat"); - const auto c = rad_pkg->template Param("c"); - const auto arad = rad_pkg->template Param("arad"); - const auto outer_max = rad_pkg->template Param("outer_iteration_max"); - const auto inner_max = rad_pkg->template Param("inner_iteration_max"); - const auto outer_tol = rad_pkg->template Param("outer_iteration_tol"); - const auto inner_tol = rad_pkg->template Param("inner_iteration_tol"); + // Extract radiation and moments package and params + auto &moments_pkg = pm->packages.Get("moments"); + const auto chat = moments_pkg->template Param("chat"); + const auto c = moments_pkg->template Param("c"); + const auto arad = moments_pkg->template Param("arad"); + const auto outer_max = moments_pkg->template Param("outer_iteration_max"); + const auto inner_max = moments_pkg->template Param("inner_iteration_max"); + const auto outer_tol = moments_pkg->template Param("outer_iteration_tol"); + const auto inner_tol = moments_pkg->template Param("inner_iteration_tol"); // Extract rotating frame quantities Real om0 = 0.0; @@ -80,7 +80,7 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { // Prepare scratch pad memory // const int ncells1 = ib.e - ib.s + 1 + 2 * parthenon::Globals::nghost; // int scr_size = ScratchPad1D::shmem_size(ncells1) * 12; - // const int scr_level = rad_pkg->template Param("scr_level"); + // const int scr_level = moments_pkg->template Param("scr_level"); parthenon::par_for( DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -182,7 +182,7 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { } //---------------------------------------------------------------------------------------- -//! \fn TaskStatus Radiation::MatterCouplingSimpleImpl +//! \fn TaskStatus Moments::MatterCouplingSimpleImpl //! \brief Implementation for "full" radiation-matter coupling source template TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { @@ -201,14 +201,14 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { auto de_switch = gas_pkg->template Param("de_switch"); // Extract radiation package and params - auto &rad_pkg = pm->packages.Get("moments"); - const auto chat = rad_pkg->template Param("chat"); - const auto c = rad_pkg->template Param("c"); - const auto arad = rad_pkg->template Param("arad"); - const auto outer_max = rad_pkg->template Param("outer_iteration_max"); - const auto inner_max = rad_pkg->template Param("inner_iteration_max"); - const auto outer_tol = rad_pkg->template Param("outer_iteration_tol"); - const auto inner_tol = rad_pkg->template Param("inner_iteration_tol"); + auto &moments_pkg = pm->packages.Get("moments"); + const auto chat = moments_pkg->template Param("chat"); + const auto c = moments_pkg->template Param("c"); + const auto arad = moments_pkg->template Param("arad"); + const auto outer_max = moments_pkg->template Param("outer_iteration_max"); + const auto inner_max = moments_pkg->template Param("inner_iteration_max"); + const auto outer_tol = moments_pkg->template Param("outer_iteration_tol"); + const auto inner_tol = moments_pkg->template Param("inner_iteration_tol"); // Extract rotating frame quantities Real om0 = 0.0; @@ -234,7 +234,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // Prepare scratch pad memory // const int ncells1 = ib.e - ib.s + 1 + 2 * parthenon::Globals::nghost; // int scr_size = ScratchPad1D::shmem_size(ncells1) * 12; - // const int scr_level = rad_pkg->template Param("scr_level"); + // const int scr_level = moments_pkg->template Param("scr_level"); parthenon::par_for( DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -421,6 +421,6 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { return TaskStatus::complete; } -} // namespace Radiation +} // namespace Moments -#endif // RADIATION_MOMENTS_MATTER_COUPLING_HPP_ +#endif // RADIATION_MOMENTS_MATTER_COUPLING_HPP_ diff --git a/src/radiation/moments/moments.cpp b/src/radiation/moments/moments.cpp index 402af1c8..3e546ac6 100644 --- a/src/radiation/moments/moments.cpp +++ b/src/radiation/moments/moments.cpp @@ -31,17 +31,17 @@ using ArtemisUtils::MeanOpacity; using ArtemisUtils::MeanScattering; using ArtemisUtils::VI; -namespace Radiation { +namespace Moments { //---------------------------------------------------------------------------------------- -//! \fn StateDescriptor Radiation::Initialize -//! \brief Adds intialization function for radiation hydrodynamics package +//! \fn StateDescriptor Moments::Initialize +//! \brief Adds intialization function for moments package std::shared_ptr Initialize(ParameterInput *pin, ArtemisUtils::Constants &constants) { - auto radiation = std::make_shared("moments"); - Params ¶ms = radiation->AllParams(); + auto moments = std::make_shared("moments"); + Params ¶ms = moments->AllParams(); // Metadata flags - auto MetadataMoments = radiation->GetMetadataFlag(); + auto MetadataMoments = moments->GetMetadataFlag(); auto MetadataOperatorSplit = Metadata::GetUserFlag("OperatorSplit"); // Closure type @@ -106,13 +106,14 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("full_coupling", pin->GetOrAddBoolean("radiation/moment", "full_coupling", true)); - // We stuff some constants into params so that can be used in post-processing + // Radiation constants (including chat for Moments) + // NOTE(@pdmullen): These are also stored in top level radiation package... const Real light = constants.GetCCode(); params.Add("c", light); - const Real creduc = pin->GetOrAddReal("radiation/moment", "creduc", 1.0); - params.Add("chat", light / creduc); const Real arad = constants.GetARCode(); params.Add("arad", arad); + const Real creduc = pin->GetOrAddReal("radiation/moment", "creduc", 1.0); + params.Add("chat", light / creduc); // Floors const Real efloor = pin->GetOrAddReal("radiation/moment", "efloor", 1.0e-20); @@ -151,7 +152,7 @@ std::shared_ptr Initialize(ParameterInput *pin, MetadataOperatorSplit}); ArtemisUtils::EnrollArtemisRefinementOps(m, coords); m.SetSparseThresholds(0.0, 0.0, 0.0); - radiation->AddSparsePool(m, control_field, fluidids); + moments->AddSparsePool(m, control_field, fluidids); // Conserved Flux m = Metadata({Metadata::Cell, Metadata::Vector, Metadata::Conserved, @@ -160,7 +161,7 @@ std::shared_ptr Initialize(ParameterInput *pin, std::vector({3})); ArtemisUtils::EnrollArtemisRefinementOps(m, coords); m.SetSparseThresholds(0.0, 0.0, 0.0); - radiation->AddSparsePool(m, control_field, fluidids); + moments->AddSparsePool(m, control_field, fluidids); // Primitive Energy Density m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, Metadata::OneCopy, @@ -168,7 +169,7 @@ std::shared_ptr Initialize(ParameterInput *pin, MetadataOperatorSplit}); ArtemisUtils::EnrollArtemisRefinementOps(m, coords); m.SetSparseThresholds(0.0, 0.0, 0.0); - radiation->AddSparsePool(m, control_field, fluidids); + moments->AddSparsePool(m, control_field, fluidids); // Primitive Pressure (and associated Riemann pressures) m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, Metadata::OneCopy, @@ -176,7 +177,7 @@ std::shared_ptr Initialize(ParameterInput *pin, MetadataOperatorSplit}); ArtemisUtils::EnrollArtemisRefinementOps(m, coords); m.SetSparseThresholds(0.0, 0.0, 0.0); - radiation->AddSparsePool(m, control_field, fluidids); + moments->AddSparsePool(m, control_field, fluidids); // Primitive Reduced Flux m = Metadata({Metadata::Cell, Metadata::Vector, Metadata::Derived, Metadata::Intensive, @@ -185,7 +186,7 @@ std::shared_ptr Initialize(ParameterInput *pin, std::vector({3})); ArtemisUtils::EnrollArtemisRefinementOps(m, coords); m.SetSparseThresholds(0.0, 0.0, 0.0); - radiation->AddSparsePool(m, control_field, fluidids); + moments->AddSparsePool(m, control_field, fluidids); // Radiation refinement criterion const std::string refine_field = @@ -217,42 +218,42 @@ std::shared_ptr Initialize(ParameterInput *pin, // Cartesian if (coords == G::cartesian) { if (ref_dens) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } else if (ref_pres) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } // Spherical } else if (coords == G::spherical1D) { if (ref_dens) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } else if (ref_pres) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } } else if (coords == G::spherical2D) { if (ref_dens) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } else if (ref_pres) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } } else if (coords == G::spherical3D) { if (ref_dens) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } else if (ref_pres) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } // Cylindrical } else if (coords == G::cylindrical) { if (ref_dens) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } else if (ref_pres) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } // Axisymmetric } else if (coords == G::axisymmetric) { if (ref_dens) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } else if (ref_pres) { - radiation->CheckRefinementBlock = ScalarFirstDerivative; + moments->CheckRefinementBlock = ScalarFirstDerivative; } } } else if (ref_mag) { @@ -262,19 +263,19 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("refine_thr", rthr); params.Add("deref_thr", dthr); if (ref_dens) { - radiation->CheckRefinementBlock = ScalarMagnitude; + moments->CheckRefinementBlock = ScalarMagnitude; } else if (ref_pres) { - radiation->CheckRefinementBlock = ScalarMagnitude; + moments->CheckRefinementBlock = ScalarMagnitude; } } } - return radiation; + return moments; } //---------------------------------------------------------------------------------------- -//! \fn TaskStatus Radiation::CalculateFluxes -//! \brief Evaluates advective fluxes for radiation evolution +//! \fn TaskStatus Moments::CalculateFluxes +//! \brief Evaluates advective fluxes for moments evolution TaskStatus CalculateFluxes(MeshData *md) { auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; @@ -305,8 +306,8 @@ TaskStatus CalculateFluxes(MeshData *md) { } //---------------------------------------------------------------------------------------- -//! \fn TaskStatus Radiation::FluxSource -//! \brief Evaluates coordinate terms from advective fluxes for radiation evolution +//! \fn TaskStatus Moments::FluxSource +//! \brief Evaluates coordinate terms from advective fluxes for moments evolution TaskStatus FluxSource(MeshData *md, const Real dt) { auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; @@ -348,9 +349,9 @@ TaskStatus MatterCoupling(MeshData *u0, const Real dt) { if (!(do_gas)) return TaskStatus::complete; // Extract moments package and params - auto &radiation_pkg = pm->packages.Get("moments"); - auto closure_type = radiation_pkg->template Param("closure_type"); - auto full_coupling = radiation_pkg->template Param("full_coupling"); + auto &moments_pkg = pm->packages.Get("moments"); + auto closure_type = moments_pkg->template Param("closure_type"); + auto full_coupling = moments_pkg->template Param("full_coupling"); // Call MatterCoupling with appropriate GEOM, Fluid, and Closure type given coupling if (closure_type == Closure::m1) { @@ -380,4 +381,4 @@ template TaskStatus MatterCoupling(MD *u0, const Real dt); template TaskStatus MatterCoupling(MD *u0, const Real dt); template TaskStatus MatterCoupling(MD *u0, const Real dt); -} // namespace Radiation +} // namespace Moments diff --git a/src/radiation/moments/moments.hpp b/src/radiation/moments/moments.hpp index f6be1f71..eb8418f4 100644 --- a/src/radiation/moments/moments.hpp +++ b/src/radiation/moments/moments.hpp @@ -18,7 +18,7 @@ #include "utils/integrators/artemis_integrator.hpp" #include "utils/units.hpp" -namespace Radiation { +namespace Moments { //---------------------------------------------------------------------------------------- std::shared_ptr Initialize(ParameterInput *pin, @@ -40,12 +40,12 @@ TaskCollection MomentsTasks(Mesh *pmesh, const SimTime &tm, parthenon::LowStorageIntegrator *integrator); //---------------------------------------------------------------------------------------- -//! \fn Real Radiation::EstimateTimeStep +//! \fn Real Moments::EstimateTimeStep //! \brief Not enrolled in parthenon's determination for global dt template Real EstimateTimeStep(parthenon::Mesh *pmesh) { - auto &radiation_pkg = pmesh->packages.Get("moments"); - auto ¶ms = radiation_pkg->AllParams(); + auto &moments_pkg = pmesh->packages.Get("moments"); + auto ¶ms = moments_pkg->AllParams(); Real dxmin = Big(); if constexpr (geometry::is_cartesian()) { @@ -71,7 +71,7 @@ Real EstimateTimeStep(parthenon::Mesh *pmesh) { // Compute minimum dx Real min_dx = Big(); parthenon::par_reduce( - parthenon::loop_pattern_mdrange_tag, "Radiation::EstimateTimestepMesh", + parthenon::loop_pattern_mdrange_tag, "Moments::EstimateTimestepMesh", DevExecSpace(), 0, md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &ldx_m) { // Extract coordinates @@ -97,7 +97,7 @@ Real EstimateTimeStep(parthenon::Mesh *pmesh) { } //---------------------------------------------------------------------------------------- -//! \fn Real Radiation::EddingtonFactor +//! \fn Real Moments::EddingtonFactor //! \brief Computes Eddington factor given closure model template KOKKOS_INLINE_FUNCTION Real EddingtonFactor(const Real f) { @@ -113,7 +113,7 @@ KOKKOS_INLINE_FUNCTION Real EddingtonFactor(const Real f) { } //---------------------------------------------------------------------------------------- -//! \fn std::array Radiation::EddingtonTensor +//! \fn std::array Moments::EddingtonTensor //! \brief Computes entries of Eddington tensor given closure model template KOKKOS_INLINE_FUNCTION std::array @@ -139,7 +139,7 @@ EddingtonTensor(const std::array fred) { } //---------------------------------------------------------------------------------------- -//! \fn std::tuple Radiation::WaveSpeed +//! \fn std::tuple Moments::WaveSpeed //! \brief Computes wavespeed given closure model template KOKKOS_INLINE_FUNCTION std::tuple WaveSpeed(const Real mu, const Real f) { @@ -160,7 +160,7 @@ KOKKOS_INLINE_FUNCTION std::tuple WaveSpeed(const Real mu, const Rea } //---------------------------------------------------------------------------------------- -//! \fn std::array Radiation::NormalizeFlux +//! \fn std::array Moments::NormalizeFlux //! \brief Normalize radiation flux KOKKOS_INLINE_FUNCTION std::array NormalizeFlux(const Real fx1, const Real fx2, const Real fx3) { @@ -174,7 +174,7 @@ std::array NormalizeFlux(const Real fx1, const Real fx2, const Real fx3 } //---------------------------------------------------------------------------------------- -//! \fn Real Radiation::FleckFactor +//! \fn Real Moments::FleckFactor //! \brief Returns Fleck factor dB/dE KOKKOS_INLINE_FUNCTION Real FleckFactor(const Real ar, const Real T, const Real cv) { @@ -182,7 +182,7 @@ Real FleckFactor(const Real ar, const Real T, const Real cv) { } //---------------------------------------------------------------------------------------- -//! \fn std::array Radiation::SolveRadFlux +//! \fn std::array Moments::SolveRadFlux //! \brief //! //! Invert this matrix: @@ -211,6 +211,6 @@ std::array SolveRadFlux(const Real a, const Real b, idet}; } -} // namespace Radiation +} // namespace Moments #endif // RADIATION_MOMENTS_MOMENTS_HPP_ diff --git a/src/radiation/moments/moments_driver.cpp b/src/radiation/moments/moments_driver.cpp index 71be2d17..0fb2676a 100644 --- a/src/radiation/moments/moments_driver.cpp +++ b/src/radiation/moments/moments_driver.cpp @@ -18,7 +18,7 @@ #include "utils/integrators/artemis_integrator.hpp" #include "utils/units.hpp" -namespace Radiation { +namespace Moments { //---------------------------------------------------------------------------------------- //! \fn TaskListStatus MomentsDriver @@ -27,7 +27,7 @@ template TaskListStatus MomentsDriver(Mesh *pmesh, const SimTime &tm, parthenon::LowStorageIntegrator *integrator) { // Craft a series of **equal** substeps that sum to the unsplit step - const Real dtlimit = Radiation::EstimateTimeStep(pmesh); + const Real dtlimit = Moments::EstimateTimeStep(pmesh); const int nsteps = static_cast(std::ceil(integrator->dt / dtlimit)); integrator->dt = integrator->dt / nsteps; @@ -95,7 +95,7 @@ TaskCollection MomentsTasks(Mesh *pmesh, const SimTime &tm, auto start_flx_recv = tl.AddTask(none, parthenon::StartReceiveFluxCorrections, u0m); // Compute radiation fluxes - auto rad_flx = tl.AddTask(none, Radiation::CalculateFluxes, u0m.get()); + auto rad_flx = tl.AddTask(none, Moments::CalculateFluxes, u0m.get()); // Communicate and set fluxes auto send_flx = tl.AddTask( @@ -110,12 +110,11 @@ TaskCollection MomentsTasks(Mesh *pmesh, const SimTime &tm, u1.get(), g0, g1, 0.0); // Apply "coordinate source terms" - auto coord_src = - tl.AddTask(rupdate | cupdate, Radiation::FluxSource, u0m.get(), bdt); + auto coord_src = tl.AddTask(rupdate | cupdate, Moments::FluxSource, u0m.get(), bdt); // Apply matter-coupling step auto coupling = - tl.AddTask(coord_src, Radiation::MatterCoupling, u0c.get(), bdt); + tl.AddTask(coord_src, Moments::MatterCoupling, u0c.get(), bdt); // Set auxillary fields auto set_aux = @@ -155,4 +154,4 @@ template TaskCollection MomentsTasks(M *pm, const ST &t, LSI *ii template TaskCollection MomentsTasks(M *pm, const ST &t, LSI *ii); template TaskCollection MomentsTasks(M *pm, const ST &t, LSI *ii); -} // namespace Radiation +} // namespace Moments diff --git a/src/radiation/params.yaml b/src/radiation/params.yaml index d4ac03b9..826e96b5 100644 --- a/src/radiation/params.yaml +++ b/src/radiation/params.yaml @@ -149,8 +149,3 @@ radiation: _type: bool _description: "Radiation affects host fluid fields." _default: true - - - - - diff --git a/src/radiation/radiation.cpp b/src/radiation/radiation.cpp index 791dc759..304505fc 100644 --- a/src/radiation/radiation.cpp +++ b/src/radiation/radiation.cpp @@ -12,45 +12,69 @@ //======================================================================================== // Artemis includes +#include "radiation.hpp" #include "artemis.hpp" #include "geometry/geometry.hpp" #include "utils/artemis_utils.hpp" #include "utils/eos/eos.hpp" #include "utils/opacity/opacity.hpp" +#include "utils/units.hpp" using ArtemisUtils::EOS; using ArtemisUtils::MeanOpacity; using ArtemisUtils::MeanScattering; -namespace rad { - -std::shared_ptr InitializeRadDataFields(ParameterInput *pin) { - - // instantiate a state descriptor for the fields - auto rad_fs = std::make_shared("rad_fields"); - - // Only one photon species - std::vector radids = {0}; - - // Control field for sparse gas fields - const std::string control_field = "rad_fields"; - - // const int ndim = ProblemDimension(pin); - // std::string sys = pin->GetOrAddString("artemis", "coordinates", "cartesian"); - // Coordinates coords = geometry::CoordSelect(sys, ndim); - - // Absorption and scattering opacity - Metadata m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy, Metadata::Sparse}); - //ArtemisUtils::EnrollArtemisRefinementOps(m, coords); - //m.SetSparseThresholds(0.0, 0.0, 0.0); - rad_fs->AddSparsePool(m, control_field, radids); - rad_fs->AddSparsePool(m, control_field, radids); +namespace Radiation { +//---------------------------------------------------------------------------------------- +//! \fn StateDescriptor Radiation::Initialize +//! \brief Adds intialization function for radiation package +//! NOTE(@pdmullen): ...to become a top-level package for radiation utils commmon to impl +std::shared_ptr +Initialize(ParameterInput *pin, ArtemisUtils::Constants &constants, const bool do_imc) { + auto radiation = std::make_shared("radiation"); + Params ¶ms = radiation->AllParams(); + + // Metadata flags + auto MetadataRadiation = radiation->GetMetadataFlag(); + auto MetadataOperatorSplit = Metadata::GetUserFlag("OperatorSplit"); + + // Radiation constants (including chat for Moments) + // NOTE(@pdmullen): We require constants again in moments sector (when + // moments enabled) due to our anonymous flux machinery + const Real light = constants.GetCCode(); + params.Add("c", light); + const Real arad = constants.GetARCode(); + params.Add("arad", arad); + const Real creduc = pin->GetOrAddReal("radiation/moment", "creduc", 1.0); + params.Add("chat", light / creduc); + + // Add derived radiation fields expected by Jaybenne + if (do_imc) { + // Number of radiation species (i.e., groups) + const int nspecies = pin->GetOrAddInteger("radiation/imc", "nspecies", 1); + params.Add("nspecies", nspecies); + PARTHENON_REQUIRE(nspecies == 1, "Jaybenne IMC only works with nspecies=1!"); + std::vector fluidids; + for (int n = 0; n < nspecies; ++n) + fluidids.push_back(n); + + // Control field for sparse gas fields + const std::string control_field = rad::opac::absorption::name(); + + // Absorption and scattering opacity + Metadata m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy, + Metadata::Sparse, MetadataRadiation, MetadataOperatorSplit}); + radiation->AddSparsePool(m, control_field, fluidids); + radiation->AddSparsePool(m, control_field, fluidids); + } - // return rad_fields state descriptor - return rad_fs; + return radiation; } -TaskStatus UpRadDataFields(MeshData *md) { +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus Radiation::SetOpacities +//! \brief Routine to set opacitiy fields (when required, e.g., for Jaybenne IMC) +TaskStatus SetOpacities(MeshData *md) { using parthenon::MakePackDescriptor; auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; @@ -60,21 +84,21 @@ TaskStatus UpRadDataFields(MeshData *md) { MeanOpacity opacity_d = gas_pkg->template Param("opacity_d"); MeanScattering scattering_d = gas_pkg->template Param("scattering_d"); - // Packing and indexing (TODO: use dust sie, density) - static auto desc = MakePackDescriptor< - gas::prim::density, gas::prim::sie, - rad::opac::absorption, rad::opac::scattering>(resolved_pkgs.get()); + // Packing and indexing + // TODO(): Will eventually incorporate other fluids + static auto desc = + MakePackDescriptor(resolved_pkgs.get()); auto vmesh = desc.GetPack(md); - const int nblocks = md->NumBlocks(); - IndexRange ib = md->GetBoundsI(IndexDomain::interior); - IndexRange jb = md->GetBoundsJ(IndexDomain::interior); - IndexRange kb = md->GetBoundsK(IndexDomain::interior); + IndexRange ib = md->GetBoundsI(IndexDomain::entire); + IndexRange jb = md->GetBoundsJ(IndexDomain::entire); + IndexRange kb = md->GetBoundsK(IndexDomain::entire); + // Set opacities parthenon::par_for( DEFAULT_LOOP_PATTERN, "ConsToPrim", parthenon::DevExecSpace(), 0, md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { - const Real &rho = vmesh(b, gas::prim::density(), k, j, i); const Real &sie = vmesh(b, gas::prim::sie(), k, j, i); const Real temp = eos_d.TemperatureFromDensityInternalEnergy(rho, sie); @@ -88,8 +112,10 @@ TaskStatus UpRadDataFields(MeshData *md) { return TaskStatus::complete; } -// task collection for updating data fields -TaskCollection UpdateRadDataFields(Mesh *pmesh) { +//---------------------------------------------------------------------------------------- +//! \fn TaskCollection Radiation::UpdateRadiationFields +//! \brief TaskCollection to set radiation fields (when required, e.g., for Jaybenne IMC) +TaskCollection UpdateRadiationFields(Mesh *pmesh) { TaskCollection tc; TaskID none(0); const int num_partitions = pmesh->DefaultNumPartitions(); @@ -97,10 +123,10 @@ TaskCollection UpdateRadDataFields(Mesh *pmesh) { for (int i = 0; i < num_partitions; i++) { auto &tl = reg[i]; auto &base = pmesh->mesh_data.GetOrAdd("base", i); - auto upradf = tl.AddTask(none, UpRadDataFields, base.get()); + auto set_opac = tl.AddTask(none, SetOpacities, base.get()); } return tc; } -} // namespace rad +} // namespace Radiation diff --git a/src/radiation/radiation.hpp b/src/radiation/radiation.hpp index e275b2ec..6d652d70 100644 --- a/src/radiation/radiation.hpp +++ b/src/radiation/radiation.hpp @@ -10,9 +10,20 @@ // license in this material to reproduce, prepare derivative works, distribute copies to // the public, perform publicly and display publicly, and to permit others to do so. //======================================================================================== +#ifndef RADIATION_RADIATION_HPP_ +#define RADIATION_RADIATION_HPP_ -namespace rad { -std::shared_ptr InitializeRadDataFields(ParameterInput *pin); -TaskStatus UpRadDataFields(MeshData *md); -TaskCollection UpdateRadDataFields(Mesh *pmesh); -} // namespace rad +#include "artemis.hpp" +#include "utils/units.hpp" + +namespace Radiation { + +std::shared_ptr +Initialize(ParameterInput *pin, ArtemisUtils::Constants &constants, const bool do_imc); + +TaskStatus SetOpacities(MeshData *md); +TaskCollection UpdateRadiationFields(Mesh *pmesh); + +} // namespace Radiation + +#endif // RADIATION_RADIATION_HPP_ \ No newline at end of file diff --git a/src/utils/fluxes/fluid_fluxes.hpp b/src/utils/fluxes/fluid_fluxes.hpp index f51b0318..2fab7a50 100644 --- a/src/utils/fluxes/fluid_fluxes.hpp +++ b/src/utils/fluxes/fluid_fluxes.hpp @@ -377,7 +377,7 @@ TaskStatus FluxSourceImpl(MeshData *md, PKG &pkg, PRIM vp, CONS vcons, FAC const Real &fy = vp_(b, IVY, k, j, i); const Real &fz = vp_(b, IVZ, k, j, i); const Real ff = std::sqrt(SQR(fx) + SQR(fy) + SQR(fz)); - const Real chi = Radiation::EddingtonFactor(ff); + const Real chi = Moments::EddingtonFactor(ff); wdt *= (3.0 * chi - 1.0) * hcchat_ / (ff + Fuzz()); } diff --git a/src/utils/fluxes/riemann/hlle.hpp b/src/utils/fluxes/riemann/hlle.hpp index 7d5826f5..2deeab44 100644 --- a/src/utils/fluxes/riemann/hlle.hpp +++ b/src/utils/fluxes/riemann/hlle.hpp @@ -273,10 +273,10 @@ struct RiemannSolver(fl); - const Real chir = Radiation::EddingtonFactor(fr); - const auto [sla, slb] = Radiation::WaveSpeed(nlx, fl); - const auto [sra, srb] = Radiation::WaveSpeed(nrx, fr); + const Real chil = Moments::EddingtonFactor(fl); + const Real chir = Moments::EddingtonFactor(fr); + const auto [sla, slb] = Moments::WaveSpeed(nlx, fl); + const auto [sra, srb] = Moments::WaveSpeed(nrx, fr); const Real sl = std::min(sla, slb); const Real sr = std::max(sra, srb); diff --git a/src/utils/fluxes/riemann/llf.hpp b/src/utils/fluxes/riemann/llf.hpp index 49810689..f30abe6b 100644 --- a/src/utils/fluxes/riemann/llf.hpp +++ b/src/utils/fluxes/riemann/llf.hpp @@ -223,10 +223,10 @@ struct RiemannSolver(fl); - const Real chir = Radiation::EddingtonFactor(fr); - const auto [sla, slb] = Radiation::WaveSpeed(nlx, fl); - const auto [sra, srb] = Radiation::WaveSpeed(nrx, fr); + const Real chil = Moments::EddingtonFactor(fl); + const Real chir = Moments::EddingtonFactor(fr); + const auto [sla, slb] = Moments::WaveSpeed(nlx, fl); + const auto [sra, srb] = Moments::WaveSpeed(nrx, fr); const Real sl = std::min(sla, slb); const Real sr = std::max(sra, srb); From b2e2ec69983b71605c9720be22cdecd767ef90fc Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Fri, 20 Jun 2025 10:49:52 -0400 Subject: [PATCH 04/10] Opt in to compiler timings --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 410f5e84..de3559f4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,7 +23,7 @@ option(ARTEMIS_ENABLE_HIP "Enable hip for artemis and all dependencies" OFF) option(ARTEMIS_ENABLE_HDF5 "Enable HDF5 for artemis and all dependencies" ON) option(ARTEMIS_ENABLE_MPI "Enable MPI for artemis and all dependencies" ON) option(ARTEMIS_ENABLE_OPENMP "Enable OpenMP for artemis and parthenon" OFF) -option(ARTEMIS_ENABLE_COMPILE_TIMING "Enable timing of compilation of artemis" ON) +option(ARTEMIS_ENABLE_COMPILE_TIMING "Enable timing of compilation of artemis" OFF) option(ARTEMIS_ENABLE_ASAN "Enable AddressSanitizer to detect memory errors" OFF) option(ARTEMIS_PORTABLE_RESTART "Enable portable data types for REBOUND restarts" ON) From d64fa64156f53a0b60c258533dbe64c16ffcfc80 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Fri, 20 Jun 2025 10:50:07 -0400 Subject: [PATCH 05/10] Consistency in ordering --- src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a661fa15..29a2d9de 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -71,8 +71,8 @@ set (SRC_LIST pgen/strat.hpp pgen/thermalization.hpp - radiation/radiation.hpp radiation/radiation.cpp + radiation/radiation.hpp radiation/imc/imc_driver.cpp radiation/imc/imc.hpp radiation/moments/moments.cpp From 14399b48e4d98463dce5d5ae9be93053daf7cb32 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Fri, 20 Jun 2025 10:50:19 -0400 Subject: [PATCH 06/10] Bump jaybenne --- external/jaybenne | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/jaybenne b/external/jaybenne index 6b19009f..36fbd743 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit 6b19009f8878b2843b31feb2b066e27065acdc3d +Subproject commit 36fbd74363aa9356f2e46f956368bc4b9c2ce8df From d73a7f4f62774d8b4b1bd2acea789086d1ff57d9 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Fri, 20 Jun 2025 11:00:36 -0400 Subject: [PATCH 07/10] Extract params from radiation pkg --- tst/scripts/radiation/rad_shock_cgs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tst/scripts/radiation/rad_shock_cgs.py b/tst/scripts/radiation/rad_shock_cgs.py index efc11b14..440a3160 100644 --- a/tst/scripts/radiation/rad_shock_cgs.py +++ b/tst/scripts/radiation/rad_shock_cgs.py @@ -73,7 +73,7 @@ def analyze(): os.path.join(artemis.get_data_dir(), "{}.out1.final.phdf".format(_file_id)), "r" ) as f: cv = f["Params"].attrs["gas/cv"] - ar = f["Params"].attrs["moments/arad"] + ar = f["Params"].attrs["radiation/arad"] xm = f["Locations/x"][...].ravel() xc = 0.5 * (xm[:-1] + xm[1:]) sie = f["gas.prim.sie_0"][...].ravel() From 5e71ed0cbdb2d92c60d0f8bdf7e2fa28fe46d7aa Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Wed, 25 Jun 2025 10:48:52 -0400 Subject: [PATCH 08/10] Bump jaybenne --- external/jaybenne | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/jaybenne b/external/jaybenne index 36fbd743..d3ecb7a1 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit 36fbd74363aa9356f2e46f956368bc4b9c2ce8df +Subproject commit d3ecb7a1db53a8d8abb96fb40ef02a4ef549a1ac From 88388f40206c9fa6f189f2e28a9c4d162a7f300a Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Wed, 25 Jun 2025 18:23:56 -0400 Subject: [PATCH 09/10] Bump jaybenne again --- external/jaybenne | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/jaybenne b/external/jaybenne index d3ecb7a1..13294529 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit d3ecb7a1db53a8d8abb96fb40ef02a4ef549a1ac +Subproject commit 13294529210664fb39e3d47ec5e9bdac92a13fdf From b8e8d64e5b46da87a6cf5de73bba8235efa3086c Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Wed, 25 Jun 2025 19:31:22 -0400 Subject: [PATCH 10/10] Bump jaybenne --- external/jaybenne | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/jaybenne b/external/jaybenne index 13294529..9bd9e13b 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit 13294529210664fb39e3d47ec5e9bdac92a13fdf +Subproject commit 9bd9e13b5bf21e8d70f8c424ba3d9948d9cf99f5