From 5348c4debbc6759e5714a834d0493a9bb43c6e60 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 7 Jul 2026 12:10:38 -0600 Subject: [PATCH 1/9] Prepare for multigroup radiation in artemis. + Move gas opacity initialization inside radiation initialization. + Update jaybenne, singularity submodules. + Add scattering powerlaw to artemis opacity variants. --- external/jaybenne | 2 +- external/singularity-eos | 2 +- external/singularity-opac | 2 +- src/CMakeLists.txt | 2 + src/artemis.cpp | 2 +- src/artemis.hpp | 3 + src/gas/gas.cpp | 98 --------------------------- src/radiation/gas_opacity.cpp | 123 ++++++++++++++++++++++++++++++++++ src/radiation/gas_opacity.hpp | 29 ++++++++ src/radiation/radiation.cpp | 28 +++++++- src/radiation/radiation.hpp | 6 +- src/utils/opacity/opacity.hpp | 1 + 12 files changed, 192 insertions(+), 106 deletions(-) create mode 100644 src/radiation/gas_opacity.cpp create mode 100644 src/radiation/gas_opacity.hpp diff --git a/external/jaybenne b/external/jaybenne index 1320212e..bebe963e 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit 1320212e6de773a9a964e680b59e488d0d193baf +Subproject commit bebe963e757015994053140e7e9c84ea8fed02b8 diff --git a/external/singularity-eos b/external/singularity-eos index 44a61202..271afd6e 160000 --- a/external/singularity-eos +++ b/external/singularity-eos @@ -1 +1 @@ -Subproject commit 44a612021ba8fa8dd771d638c8543a671fa8947e +Subproject commit 271afd6e52106333ea551462a3b256f982547c03 diff --git a/external/singularity-opac b/external/singularity-opac index cdd365ff..74a791f0 160000 --- a/external/singularity-opac +++ b/external/singularity-opac @@ -1 +1 @@ -Subproject commit cdd365ffd56ccdab1ae4b4410d276d0a46a4a314 +Subproject commit 74a791f0f4dd285c304a77cd38f3136639c25637 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ac8db0ee..deb753bb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -79,6 +79,8 @@ set (SRC_LIST pgen/strat.hpp pgen/thermalization.hpp + radiation/gas_opacity.cpp + radiation/gas_opacity.hpp radiation/radiation.cpp radiation/radiation.hpp radiation/imc/imc_driver.cpp diff --git a/src/artemis.cpp b/src/artemis.cpp index 7139c290..0de3ac5c 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -185,7 +185,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { // Operator split radiation if (do_radiation) { // Top-level radiation package - packages.Add(Radiation::Initialize(pin.get(), constants, do_imc)); + packages.Add(Radiation::Initialize(pin.get(), units, constants, do_imc)); // Select between Jaybenne IMC or Moments if (do_imc) { auto eos_h = packages.Get("gas")->Param("eos_h"); diff --git a/src/artemis.hpp b/src/artemis.hpp index c0fe684e..c6ead3fb 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -183,6 +183,9 @@ enum class ArtemisBC { // Tensor indexing (currently used in radiation moments) enum TensIdx { X11 = 0, X22 = 1, X33 = 2, X23 = 3, X13 = 4, X12 = 5 }; +// Radiation solver frequency type enum +enum class FrequencyType { gray, multigroup }; + // Floating point limits template KOKKOS_FORCEINLINE_FUNCTION constexpr auto Big() { diff --git a/src/gas/gas.cpp b/src/gas/gas.cpp index f643b472..b16fe7f9 100644 --- a/src/gas/gas.cpp +++ b/src/gas/gas.cpp @@ -28,7 +28,6 @@ #include "utils/fluxes/fluid_fluxes.hpp" #include "utils/history.hpp" #include "utils/integrators/artemis_integrator.hpp" -#include "utils/opacity/opacity.hpp" #include "utils/refinement/amr_criteria.hpp" #include "utils/units.hpp" @@ -43,7 +42,6 @@ std::shared_ptr Initialize(ParameterInput *pin, ArtemisUtils::Units &units, ArtemisUtils::Constants &constants, Packages_t &packages) { - using namespace singularity::photons; auto gas = std::make_shared("gas"); Params ¶ms = gas->AllParams(); @@ -207,102 +205,6 @@ std::shared_ptr Initialize(ParameterInput *pin, } params.Add("rsolver", riemann_solver); - // Opacity models - const Real time = units.GetTimeCodeToPhysical(); - const Real mass = units.GetMassCodeToPhysical(); - const Real length = units.GetLengthCodeToPhysical(); - const Real temp = units.GetTemperatureCodeToPhysical(); - - // Absorption opacity model - std::string opacity_model_name = - pin->GetOrAddString("gas/opacity/absorption", "opacity_model", "constant"); - - // Mean absorption opacity (either read from table or uses model - ArtemisUtils::MeanOpacity opacity; - if (opacity_model_name == "table") { - std::string table_filename = - pin->GetString("gas/opacity/absorption", "opacity_table"); - opacity = - singularity::photons::MeanNonCGSUnits( - singularity::photons::MeanOpacityBase(table_filename), time, mass, length, - temp); - } else { - // Instantiate mean absorption opacity object (i.e., table) - const Real lRhoMin_a = pin->GetOrAddReal("gas/opacity/absorption", "lRhoMin", -1.0); - const Real lRhoMax_a = pin->GetOrAddReal("gas/opacity/absorption", "lRhoMax", 1.0); - const int NRho_a = pin->GetOrAddInteger("gas/opacity/absorption", "NRho", 2); - const Real lTMin_a = pin->GetOrAddReal("gas/opacity/absorption", "lTMin", -1.0); - const Real lTMax_a = pin->GetOrAddReal("gas/opacity/absorption", "lTMax", 1.0); - const int NT_a = pin->GetOrAddInteger("gas/opacity/absorption", "NT", 2); - - if (opacity_model_name == "none") { - auto model = Gray(0.0); - opacity = - singularity::photons::MeanNonCGSUnits( - singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, - lTMin_a, lTMax_a, NT_a), - time, mass, length, temp); - } else if (opacity_model_name == "constant") { - const Real kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "kappa_a", 0.0); - auto model = Gray(kappa_a); - opacity = - singularity::photons::MeanNonCGSUnits( - singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, - lTMin_a, lTMax_a, NT_a), - time, mass, length, temp); - } else if (opacity_model_name == "powerlaw") { - const Real coef_kappa_a = - pin->GetOrAddReal("gas/opacity/absorption", "coef_kappa_a", 0.0); - const Real rho_exp = pin->GetOrAddReal("gas/opacity/absorption", "rho_exp", 0.0); - const Real temp_exp = pin->GetOrAddReal("gas/opacity/absorption", "temp_exp", 0.0); - auto model = PowerLaw(coef_kappa_a, rho_exp, temp_exp); - opacity = - singularity::photons::MeanNonCGSUnits( - singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, - lTMin_a, lTMax_a, NT_a), - time, mass, length, temp); - } else { - PARTHENON_FAIL("Opacity model not recognized!"); - } - } - params.Add("opacity_h", opacity); - params.Add("opacity_d", opacity.GetOnDevice()); - - // Scattering opacity model - std::string scattering_model_name = - pin->GetOrAddString("gas/opacity/scattering", "scattering_model", "none"); - - // Instantiate mean scattering opacity object (i.e., table) - const Real lRhoMin_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMin", -1.0); - const Real lRhoMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMax", 1.0); - const int NRho_s = pin->GetOrAddInteger("gas/opacity/scattering", "NRho", 2); - const Real lTMin_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMin", -1.0); - const Real lTMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMax", 1.0); - const int NT_s = pin->GetOrAddInteger("gas/opacity/scattering", "NT", 2); - - ArtemisUtils::MeanScattering scattering; - if (scattering_model_name == "none") { - auto smodel = GrayS(0.0, 1.0); - scattering = - singularity::photons::MeanNonCGSUnitsS( - singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s, - lTMin_s, lTMax_s, NT_s), - time, mass, length, temp); - } else if (scattering_model_name == "constant") { - const Real kappa_s = pin->GetOrAddReal("gas/opacity/scattering", "kappa_s", 0.0); - auto smodel = GrayS(kappa_s, 1.0); - scattering = - singularity::photons::MeanNonCGSUnitsS( - singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s, - lTMin_s, lTMax_s, NT_s), - time, mass, length, temp); - } else { - PARTHENON_FAIL("Scattering model not recognized!"); - } - - params.Add("scattering_h", scattering); - params.Add("scattering_d", scattering.GetOnDevice()); - // Dual energy switch // When internal > de_switch * total we use the total // The default turns off the switch diff --git a/src/radiation/gas_opacity.cpp b/src/radiation/gas_opacity.cpp new file mode 100644 index 00000000..6e20b8a4 --- /dev/null +++ b/src/radiation/gas_opacity.cpp @@ -0,0 +1,123 @@ +//======================================================================================== +// (C) (or copyright) 2026. 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. +//======================================================================================== +#include "gas_opacity.hpp" +#include "utils/opacity/opacity.hpp" + +namespace Gas { +void InitGasOpacity(ParameterInput *pin, + const ArtemisUtils::Units &units, + Params ¶ms) { + using namespace singularity::photons; + + // // get the gas state here + // auto gas = std::make_shared("gas"); + // Params ¶ms = gas->AllParams(); + + // Opacity models + const Real time = units.GetTimeCodeToPhysical(); + const Real mass = units.GetMassCodeToPhysical(); + const Real length = units.GetLengthCodeToPhysical(); + const Real temp = units.GetTemperatureCodeToPhysical(); + + // Absorption opacity model + std::string opacity_model_name = + pin->GetOrAddString("gas/opacity/absorption", "opacity_model", "constant"); + + // Mean absorption opacity (either read from table or uses model + ArtemisUtils::MeanOpacity opacity; + if (opacity_model_name == "table") { + std::string table_filename = + pin->GetString("gas/opacity/absorption", "opacity_table"); + opacity = + singularity::photons::MeanNonCGSUnits( + singularity::photons::MeanOpacityBase(table_filename), time, mass, length, + temp); + } else { + // Instantiate mean absorption opacity object (i.e., table) + const Real lRhoMin_a = pin->GetOrAddReal("gas/opacity/absorption", "lRhoMin", -1.0); + const Real lRhoMax_a = pin->GetOrAddReal("gas/opacity/absorption", "lRhoMax", 1.0); + const int NRho_a = pin->GetOrAddInteger("gas/opacity/absorption", "NRho", 2); + const Real lTMin_a = pin->GetOrAddReal("gas/opacity/absorption", "lTMin", -1.0); + const Real lTMax_a = pin->GetOrAddReal("gas/opacity/absorption", "lTMax", 1.0); + const int NT_a = pin->GetOrAddInteger("gas/opacity/absorption", "NT", 2); + + if (opacity_model_name == "none") { + auto model = Gray(0.0); + opacity = + singularity::photons::MeanNonCGSUnits( + singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, + lTMin_a, lTMax_a, NT_a), + time, mass, length, temp); + } else if (opacity_model_name == "constant") { + const Real kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "kappa_a", 0.0); + auto model = Gray(kappa_a); + opacity = + singularity::photons::MeanNonCGSUnits( + singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, + lTMin_a, lTMax_a, NT_a), + time, mass, length, temp); + } else if (opacity_model_name == "powerlaw") { + const Real coef_kappa_a = + pin->GetOrAddReal("gas/opacity/absorption", "coef_kappa_a", 0.0); + const Real rho_exp = pin->GetOrAddReal("gas/opacity/absorption", "rho_exp", 0.0); + const Real temp_exp = pin->GetOrAddReal("gas/opacity/absorption", "temp_exp", 0.0); + auto model = PowerLaw(coef_kappa_a, rho_exp, temp_exp); + opacity = + singularity::photons::MeanNonCGSUnits( + singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, + lTMin_a, lTMax_a, NT_a), + time, mass, length, temp); + } else { + PARTHENON_FAIL("Opacity model not recognized!"); + } + } + params.Add("opacity_h", opacity); + params.Add("opacity_d", opacity.GetOnDevice()); + + // Scattering opacity model + std::string scattering_model_name = + pin->GetOrAddString("gas/opacity/scattering", "scattering_model", "none"); + + // Instantiate mean scattering opacity object (i.e., table) + const Real lRhoMin_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMin", -1.0); + const Real lRhoMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMax", 1.0); + const int NRho_s = pin->GetOrAddInteger("gas/opacity/scattering", "NRho", 2); + const Real lTMin_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMin", -1.0); + const Real lTMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMax", 1.0); + const int NT_s = pin->GetOrAddInteger("gas/opacity/scattering", "NT", 2); + + ArtemisUtils::MeanScattering scattering; + if (scattering_model_name == "none") { + auto smodel = GrayS(0.0, 1.0); + scattering = + singularity::photons::MeanNonCGSUnitsS( + singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s, + lTMin_s, lTMax_s, NT_s), + time, mass, length, temp); + } else if (scattering_model_name == "constant") { + const Real kappa_s = pin->GetOrAddReal("gas/opacity/scattering", "kappa_s", 0.0); + auto smodel = GrayS(kappa_s, 1.0); + scattering = + singularity::photons::MeanNonCGSUnitsS( + singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s, + lTMin_s, lTMax_s, NT_s), + time, mass, length, temp); + } else { + PARTHENON_FAIL("Scattering model not recognized!"); + } + + params.Add("scattering_h", scattering); + params.Add("scattering_d", scattering.GetOnDevice()); +} + +} // namespace Gas diff --git a/src/radiation/gas_opacity.hpp b/src/radiation/gas_opacity.hpp new file mode 100644 index 00000000..65928ae3 --- /dev/null +++ b/src/radiation/gas_opacity.hpp @@ -0,0 +1,29 @@ +//======================================================================================== +// (C) (or copyright) 2026. 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. +//======================================================================================== +#ifndef RADIATION_GAS_OPACITY_HPP_ +#define RADIATION_GAS_OPACITY_HPP_ + +// Artemis includes +#include "utils/units.hpp" + +// Parthenon includes +#include +#include + +namespace Gas { +void InitGasOpacity(ParameterInput *pin, + const ArtemisUtils::Units &units, + Params ¶ms); +} // namespace Gas + +#endif // RADIATION_GAS_OPACITY_HPP_ diff --git a/src/radiation/radiation.cpp b/src/radiation/radiation.cpp index 9a10d34d..699bc75b 100644 --- a/src/radiation/radiation.cpp +++ b/src/radiation/radiation.cpp @@ -14,6 +14,7 @@ // Artemis includes #include "radiation.hpp" #include "artemis.hpp" +#include "gas_opacity.hpp" #include "geometry/geometry.hpp" #include "utils/artemis_utils.hpp" #include "utils/eos/eos.hpp" @@ -29,8 +30,10 @@ 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) { +std::shared_ptr Initialize(ParameterInput *pin, + ArtemisUtils::Units &units, + ArtemisUtils::Constants &constants, + const bool do_imc) { auto radiation = std::make_shared("radiation"); Params ¶ms = radiation->AllParams(); @@ -55,8 +58,20 @@ Initialize(ParameterInput *pin, ArtemisUtils::Constants &constants, const bool d params.Add("chat", light); } + // frequency type (determined below) + FrequencyType frequency_type; + // Add derived radiation fields expected by Jaybenne if (do_imc) { + // Get multigroup indicator + std::string frequency_type_name = pin->GetString("radiation/imc", "frequency_type"); + if (frequency_type_name == "gray") { + frequency_type = FrequencyType::gray; + } else if (frequency_type_name == "multigroup") { + frequency_type = FrequencyType::multigroup; + } else { + PARTHENON_FAIL("\"mcblock/frequency_type\" not recognized!"); + } // Number of radiation species (i.e., groups) const int nspecies = pin->GetOrAddInteger("radiation/imc", "nspecies", 1); params.Add("nspecies", nspecies); @@ -72,8 +87,17 @@ Initialize(ParameterInput *pin, ArtemisUtils::Constants &constants, const bool d MetadataRadiation, MetadataOperatorSplit}); radiation->AddField(m); radiation->AddField(m); + } else { + // TODO: extend MG frequency_type option to moments + frequency_type = FrequencyType::gray; } + // incorporate frequency type for gas opacity initialization + params.Add("frequency_type", frequency_type); + + // Initialize gas opacity + Gas::InitGasOpacity(pin, units, params); + // Enroll in tstart/tstop machinery ArtemisUtils::AddPackageTimeParams( params, (do_imc) ? "radiation/imc" : "radiation/moment", pin); diff --git a/src/radiation/radiation.hpp b/src/radiation/radiation.hpp index 57ffa419..ed1991f3 100644 --- a/src/radiation/radiation.hpp +++ b/src/radiation/radiation.hpp @@ -18,8 +18,10 @@ namespace Radiation { -std::shared_ptr -Initialize(ParameterInput *pin, ArtemisUtils::Constants &constants, const bool do_imc); +std::shared_ptr Initialize(ParameterInput *pin, + ArtemisUtils::Units &units, + ArtemisUtils::Constants &constants, + const bool do_imc); TaskStatus SetOpacities(MeshData *md); TaskCollection UpdateRadiationFields(Mesh *pmesh); diff --git a/src/utils/opacity/opacity.hpp b/src/utils/opacity/opacity.hpp index 981d5691..d2049c0e 100644 --- a/src/utils/opacity/opacity.hpp +++ b/src/utils/opacity/opacity.hpp @@ -30,6 +30,7 @@ using Opacity = singularity::photons::impl::Variant< // Reduced scattering variant for this codebase using Scattering = singularity::photons::impl::S_Variant< singularity::photons::NonCGSUnitsS, + singularity::photons::NonCGSUnitsS, singularity::photons::NonCGSUnitsS>; // Reduced variant for mean absorption opacities From 89685f6e6b2c064236ce5f30bdc2b37cfb3c131e Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 7 Jul 2026 17:33:39 -0600 Subject: [PATCH 2/9] Correct which package gas opacity is being accessed from (radiation). --- src/artemis.cpp | 21 ++++++++--- src/radiation/gas_opacity.cpp | 43 +++++++++++++++++++---- src/radiation/gas_opacity.hpp | 3 +- src/radiation/moments/matter_coupling.hpp | 10 +++--- src/radiation/moments/moments.cpp | 3 +- src/radiation/radiation.cpp | 11 +++--- 6 files changed, 68 insertions(+), 23 deletions(-) diff --git a/src/artemis.cpp b/src/artemis.cpp index 0de3ac5c..0a1a9eeb 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -189,10 +189,23 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { // 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(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h, - "radiation/imc")); + const auto frequency_type = + packages.Get("artemis")->Param("frequency_type"); + if (frequency_type == FrequencyType::gray) { + auto opacity_h = packages.Get("radiation")->Param("opacity_h"); + auto scattering_h = + packages.Get("radiation")->Param("scattering_h"); + packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h, + "radiation/imc")); + } else { + PARTHENON_REQUIRE(frequency_type == FrequencyType::multigroup, + "Invalid frequency_type!"); + auto opacity_h = packages.Get("radiation")->Param("mg_opacity_h"); + auto scattering_h = + packages.Get("radiation")->Param("mg_scattering_h"); + 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) { diff --git a/src/radiation/gas_opacity.cpp b/src/radiation/gas_opacity.cpp index 6e20b8a4..fe043841 100644 --- a/src/radiation/gas_opacity.cpp +++ b/src/radiation/gas_opacity.cpp @@ -14,28 +14,29 @@ #include "utils/opacity/opacity.hpp" namespace Gas { -void InitGasOpacity(ParameterInput *pin, - const ArtemisUtils::Units &units, +void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, Params ¶ms) { using namespace singularity::photons; - // // get the gas state here - // auto gas = std::make_shared("gas"); - // Params ¶ms = gas->AllParams(); - // Opacity models const Real time = units.GetTimeCodeToPhysical(); const Real mass = units.GetMassCodeToPhysical(); const Real length = units.GetLengthCodeToPhysical(); const Real temp = units.GetTemperatureCodeToPhysical(); + // Get frequency type (it should already be set in radiation Initialization) + const auto frequency_type = params.Get("frequency_type"); + // Absorption opacity model std::string opacity_model_name = pin->GetOrAddString("gas/opacity/absorption", "opacity_model", "constant"); // Mean absorption opacity (either read from table or uses model + ArtemisUtils::Opacity mg_opacity; ArtemisUtils::MeanOpacity opacity; if (opacity_model_name == "table") { + PARTHENON_REQUIRE(frequency_type == FrequencyType::gray, + "Only gray table opacity permitted, for now."); std::string table_filename = pin->GetString("gas/opacity/absorption", "opacity_table"); opacity = @@ -58,6 +59,10 @@ void InitGasOpacity(ParameterInput *pin, singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, lTMin_a, lTMax_a, NT_a), time, mass, length, temp); + if (frequency_type == FrequencyType::multigroup) { + mg_opacity = singularity::photons::NonCGSUnits( + std::move(model), time, mass, length, temp); + } } else if (opacity_model_name == "constant") { const Real kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "kappa_a", 0.0); auto model = Gray(kappa_a); @@ -66,6 +71,10 @@ void InitGasOpacity(ParameterInput *pin, singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, lTMin_a, lTMax_a, NT_a), time, mass, length, temp); + if (frequency_type == FrequencyType::multigroup) { + mg_opacity = singularity::photons::NonCGSUnits( + std::move(model), time, mass, length, temp); + } } else if (opacity_model_name == "powerlaw") { const Real coef_kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "coef_kappa_a", 0.0); @@ -77,12 +86,21 @@ void InitGasOpacity(ParameterInput *pin, singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, lTMin_a, lTMax_a, NT_a), time, mass, length, temp); + if (frequency_type == FrequencyType::multigroup) { + mg_opacity = singularity::photons::NonCGSUnits( + std::move(model), time, mass, length, temp); + } } else { PARTHENON_FAIL("Opacity model not recognized!"); } } + params.Add("opacity_h", opacity); params.Add("opacity_d", opacity.GetOnDevice()); + if (frequency_type == FrequencyType::multigroup) { + params.Add("mg_opacity_h", mg_opacity); + params.Add("mg_opacity_d", mg_opacity.GetOnDevice()); + } // Scattering opacity model std::string scattering_model_name = @@ -96,6 +114,7 @@ void InitGasOpacity(ParameterInput *pin, const Real lTMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMax", 1.0); const int NT_s = pin->GetOrAddInteger("gas/opacity/scattering", "NT", 2); + ArtemisUtils::Scattering mg_scattering; ArtemisUtils::MeanScattering scattering; if (scattering_model_name == "none") { auto smodel = GrayS(0.0, 1.0); @@ -104,6 +123,10 @@ void InitGasOpacity(ParameterInput *pin, singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s, lTMin_s, lTMax_s, NT_s), time, mass, length, temp); + if (frequency_type == FrequencyType::multigroup) { + mg_scattering = singularity::photons::NonCGSUnitsS( + std::move(smodel), time, mass, length, temp); + } } else if (scattering_model_name == "constant") { const Real kappa_s = pin->GetOrAddReal("gas/opacity/scattering", "kappa_s", 0.0); auto smodel = GrayS(kappa_s, 1.0); @@ -112,12 +135,20 @@ void InitGasOpacity(ParameterInput *pin, singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s, lTMin_s, lTMax_s, NT_s), time, mass, length, temp); + if (frequency_type == FrequencyType::multigroup) { + mg_scattering = singularity::photons::NonCGSUnitsS( + std::move(smodel), time, mass, length, temp); + } } else { PARTHENON_FAIL("Scattering model not recognized!"); } params.Add("scattering_h", scattering); params.Add("scattering_d", scattering.GetOnDevice()); + if (frequency_type == FrequencyType::multigroup) { + params.Add("mg_scattering_h", mg_scattering); + params.Add("mg_scattering_d", mg_scattering.GetOnDevice()); + } } } // namespace Gas diff --git a/src/radiation/gas_opacity.hpp b/src/radiation/gas_opacity.hpp index 65928ae3..d4a24b77 100644 --- a/src/radiation/gas_opacity.hpp +++ b/src/radiation/gas_opacity.hpp @@ -21,8 +21,7 @@ #include namespace Gas { -void InitGasOpacity(ParameterInput *pin, - const ArtemisUtils::Units &units, +void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, Params ¶ms); } // namespace Gas diff --git a/src/radiation/moments/matter_coupling.hpp b/src/radiation/moments/matter_coupling.hpp index c471f3a4..cdde2029 100644 --- a/src/radiation/moments/matter_coupling.hpp +++ b/src/radiation/moments/matter_coupling.hpp @@ -41,9 +41,10 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { // Extract gas package and params auto &gas_pkg = pm->packages.Get("gas"); + auto &rad_pkg = pm->packages.Get("radiation"); auto eos_d = gas_pkg->template Param("eos_d"); - auto opac_d = gas_pkg->template Param("opacity_d"); - auto scat_d = gas_pkg->template Param("scattering_d"); + auto opac_d = rad_pkg->template Param("opacity_d"); + auto scat_d = rad_pkg->template Param("scattering_d"); auto dflr = gas_pkg->template Param("dfloor"); auto de_switch = gas_pkg->template Param("de_switch"); @@ -210,9 +211,10 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // Extract gas package and params auto &gas_pkg = pm->packages.Get("gas"); + auto &rad_pkg = pm->packages.Get("radiation"); auto eos_d = gas_pkg->template Param("eos_d"); - auto opac_d = gas_pkg->template Param("opacity_d"); - auto scat_d = gas_pkg->template Param("scattering_d"); + auto opac_d = rad_pkg->template Param("opacity_d"); + auto scat_d = rad_pkg->template Param("scattering_d"); auto dflr = gas_pkg->template Param("dfloor"); auto de_switch = gas_pkg->template Param("de_switch"); diff --git a/src/radiation/moments/moments.cpp b/src/radiation/moments/moments.cpp index 4cb1ca85..d6fb6ca4 100644 --- a/src/radiation/moments/moments.cpp +++ b/src/radiation/moments/moments.cpp @@ -386,6 +386,7 @@ void InitMesh(parthenon::Mesh *pmesh) { PARTHENON_INSTRUMENT auto &moments_pkg = pmesh->packages.Get("moments"); auto &gas_pkg = pmesh->packages.Get("gas"); + auto &rad_pkg = pmesh->packages.Get("radiation"); const Real arad = moments_pkg->Param("arad"); const bool use_opac = moments_pkg->Param("use_opac"); @@ -410,7 +411,7 @@ void InitMesh(parthenon::Mesh *pmesh) { "coord_params"); if (use_opac) { - const auto &opac_d = gas_pkg->Param("opacity_d"); + const auto &opac_d = rad_pkg->Param("opacity_d"); const bool multi_d = pmesh->ndim >= 2; const bool three_d = pmesh->ndim == 3; parthenon::par_for( diff --git a/src/radiation/radiation.cpp b/src/radiation/radiation.cpp index 699bc75b..1ecebf8b 100644 --- a/src/radiation/radiation.cpp +++ b/src/radiation/radiation.cpp @@ -58,19 +58,18 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("chat", light); } - // frequency type (determined below) - FrequencyType frequency_type; + // frequency type (assume gray unless set below) + FrequencyType frequency_type = FrequencyType::gray; // Add derived radiation fields expected by Jaybenne if (do_imc) { // Get multigroup indicator std::string frequency_type_name = pin->GetString("radiation/imc", "frequency_type"); - if (frequency_type_name == "gray") { - frequency_type = FrequencyType::gray; - } else if (frequency_type_name == "multigroup") { + if (frequency_type_name == "multigroup") { frequency_type = FrequencyType::multigroup; } else { - PARTHENON_FAIL("\"mcblock/frequency_type\" not recognized!"); + PARTHENON_REQUIRE(frequency_type_name == "gray", + "Supported frequency_type are gray or multigroup!"); } // Number of radiation species (i.e., groups) const int nspecies = pin->GetOrAddInteger("radiation/imc", "nspecies", 1); From 6309b03505d21856eaf20a015ef58fb7d35c6a07 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Wed, 8 Jul 2026 10:44:06 -0600 Subject: [PATCH 3/9] + Fix more affected package paramater accesses. + Add frequency_type to radiation/imc input nodes. --- inputs/radiation/crooked_pipe.in | 1 + inputs/radiation/rad_shock.in | 1 + inputs/radiation/thermalization_imc.in | 1 + src/artemis.cpp | 2 +- src/radiation/radiation.cpp | 5 +++-- 5 files changed, 7 insertions(+), 3 deletions(-) diff --git a/inputs/radiation/crooked_pipe.in b/inputs/radiation/crooked_pipe.in index db51fc94..c5c80dff 100644 --- a/inputs/radiation/crooked_pipe.in +++ b/inputs/radiation/crooked_pipe.in @@ -114,3 +114,4 @@ use_ddmc = false # use DDMC? cutoff = 1.0e-6 # default is 1.0e-6 emit_temp_threshold = 1.2e4 # from ryan, 1.1e4 is background #min_swarm_occupancy = 0.5 +frequency_type = gray \ No newline at end of file diff --git a/inputs/radiation/rad_shock.in b/inputs/radiation/rad_shock.in index 7f852cff..91b0aeee 100644 --- a/inputs/radiation/rad_shock.in +++ b/inputs/radiation/rad_shock.in @@ -108,6 +108,7 @@ temp_exp = 0.0 # temp exponent for opacity powerlaw num_particles = 100000 # particle resolution use_ddmc = false # use DDMC? +frequency_type = gray rhol = 5.69 # density (left) diff --git a/inputs/radiation/thermalization_imc.in b/inputs/radiation/thermalization_imc.in index cd022f48..66f31ee7 100644 --- a/inputs/radiation/thermalization_imc.in +++ b/inputs/radiation/thermalization_imc.in @@ -85,6 +85,7 @@ kappa_a = 2.0 num_particles = 200000 use_ddmc = true +frequency_type = gray # # cfl = 0.3 diff --git a/src/artemis.cpp b/src/artemis.cpp index 0a1a9eeb..01eb09c0 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -190,7 +190,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { if (do_imc) { auto eos_h = packages.Get("gas")->Param("eos_h"); const auto frequency_type = - packages.Get("artemis")->Param("frequency_type"); + packages.Get("radiation")->Param("frequency_type"); if (frequency_type == FrequencyType::gray) { auto opacity_h = packages.Get("radiation")->Param("opacity_h"); auto scattering_h = diff --git a/src/radiation/radiation.cpp b/src/radiation/radiation.cpp index 1ecebf8b..d127ba9c 100644 --- a/src/radiation/radiation.cpp +++ b/src/radiation/radiation.cpp @@ -112,10 +112,11 @@ TaskStatus SetOpacities(MeshData *md) { auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; auto &gas_pkg = pm->packages.Get("gas"); + auto &rad_pkg = pm->packages.Get("radiation"); 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"); + MeanOpacity opacity_d = rad_pkg->template Param("opacity_d"); + MeanScattering scattering_d = rad_pkg->template Param("scattering_d"); // Packing and indexing // TODO(): Will eventually incorporate other fluids From 0fc437cd55566635c2c25b754b37de807faa6cfa Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Fri, 31 Jul 2026 18:16:12 -0600 Subject: [PATCH 4/9] Bump jaybenne to develop. --- external/jaybenne | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/jaybenne b/external/jaybenne index bebe963e..ab27f67a 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit bebe963e757015994053140e7e9c84ea8fed02b8 +Subproject commit ab27f67ade9374083a5e82cb407b760cd347a30a From db5fd89dbb98356a4c1bc3e06134c0b5563f34dd Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Fri, 31 Jul 2026 18:19:10 -0600 Subject: [PATCH 5/9] Bump singularity-opac. --- external/singularity-opac | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/singularity-opac b/external/singularity-opac index 74a791f0..31ec436b 160000 --- a/external/singularity-opac +++ b/external/singularity-opac @@ -1 +1 @@ -Subproject commit 74a791f0f4dd285c304a77cd38f3136639c25637 +Subproject commit 31ec436bb466674c91ea9fa22f6b03b1703bba51 From fa909c31bc551be885f8f152a6ebb6a928ef66e0 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Mon, 3 Aug 2026 15:52:13 -0600 Subject: [PATCH 6/9] Make minimal opacity API updates. + Bump singularity-opac. + Use AbsorptionCoefficient and ScatteringCoefficient indexed at 0. + Hard-code frequency bounds for gray singularity-opac models. --- external/jaybenne | 2 +- external/singularity-opac | 2 +- src/radiation/gas_opacity.cpp | 27 +++++++++++++++-------- src/radiation/moments/matter_coupling.hpp | 18 +++++++-------- src/radiation/moments/moments.cpp | 2 +- src/radiation/radiation.cpp | 6 ++--- src/utils/opacity/opacity.hpp | 2 +- 7 files changed, 33 insertions(+), 26 deletions(-) diff --git a/external/jaybenne b/external/jaybenne index ab27f67a..69b27fdf 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit ab27f67ade9374083a5e82cb407b760cd347a30a +Subproject commit 69b27fdf2d02448ee844df51bc60a292c7e3d57a diff --git a/external/singularity-opac b/external/singularity-opac index 31ec436b..b97515f6 160000 --- a/external/singularity-opac +++ b/external/singularity-opac @@ -1 +1 @@ -Subproject commit 31ec436bb466674c91ea9fa22f6b03b1703bba51 +Subproject commit b97515f695f96fabeca28852b5483bc0cd49528e diff --git a/src/radiation/gas_opacity.cpp b/src/radiation/gas_opacity.cpp index fe043841..73852422 100644 --- a/src/radiation/gas_opacity.cpp +++ b/src/radiation/gas_opacity.cpp @@ -27,6 +27,10 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, // Get frequency type (it should already be set in radiation Initialization) const auto frequency_type = params.Get("frequency_type"); + // TODO: read this (and consolidate gray/multigroup modes?) + const std::vector gray_bounds = {1.e12, 3.e20}; + const int NG = static_cast(gray_bounds.size()) - 1; + // Absorption opacity model std::string opacity_model_name = pin->GetOrAddString("gas/opacity/absorption", "opacity_model", "constant"); @@ -57,7 +61,8 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, opacity = singularity::photons::MeanNonCGSUnits( singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, - lTMin_a, lTMax_a, NT_a), + lTMin_a, lTMax_a, NT_a, gray_bounds, + NG), time, mass, length, temp); if (frequency_type == FrequencyType::multigroup) { mg_opacity = singularity::photons::NonCGSUnits( @@ -69,7 +74,8 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, opacity = singularity::photons::MeanNonCGSUnits( singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, - lTMin_a, lTMax_a, NT_a), + lTMin_a, lTMax_a, NT_a, gray_bounds, + NG), time, mass, length, temp); if (frequency_type == FrequencyType::multigroup) { mg_opacity = singularity::photons::NonCGSUnits( @@ -84,7 +90,8 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, opacity = singularity::photons::MeanNonCGSUnits( singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, - lTMin_a, lTMax_a, NT_a), + lTMin_a, lTMax_a, NT_a, gray_bounds, + NG), time, mass, length, temp); if (frequency_type == FrequencyType::multigroup) { mg_opacity = singularity::photons::NonCGSUnits( @@ -119,9 +126,10 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, if (scattering_model_name == "none") { auto smodel = GrayS(0.0, 1.0); scattering = - singularity::photons::MeanNonCGSUnitsS( - singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s, - lTMin_s, lTMax_s, NT_s), + singularity::photons::MeanNonCGSUnitsS( + singularity::photons::MeanSOpacityBase(smodel, lRhoMin_s, lRhoMax_s, NRho_s, + lTMin_s, lTMax_s, NT_s, gray_bounds, + NG), time, mass, length, temp); if (frequency_type == FrequencyType::multigroup) { mg_scattering = singularity::photons::NonCGSUnitsS( @@ -131,9 +139,10 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, const Real kappa_s = pin->GetOrAddReal("gas/opacity/scattering", "kappa_s", 0.0); auto smodel = GrayS(kappa_s, 1.0); scattering = - singularity::photons::MeanNonCGSUnitsS( - singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s, - lTMin_s, lTMax_s, NT_s), + singularity::photons::MeanNonCGSUnitsS( + singularity::photons::MeanSOpacityBase(smodel, lRhoMin_s, lRhoMax_s, NRho_s, + lTMin_s, lTMax_s, NT_s, gray_bounds, + NG), time, mass, length, temp); if (frequency_type == FrequencyType::multigroup) { mg_scattering = singularity::photons::NonCGSUnitsS( diff --git a/src/radiation/moments/matter_coupling.hpp b/src/radiation/moments/matter_coupling.hpp index cdde2029..12540431 100644 --- a/src/radiation/moments/matter_coupling.hpp +++ b/src/radiation/moments/matter_coupling.hpp @@ -143,7 +143,7 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { T = std::pow(B / arad, 0.25); e = eos_d.InternalEnergyFromDensityTemperature(dens, T) * dens; const Real Cv = dens * eos_d.SpecificHeatFromDensityTemperature(dens, T); - const Real a = chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); + const Real a = chat * dt * opac_d.PlanckGroupAbsorptionCoefficient(dens, T, 0); const Real fleck = FleckFactor(arad, T, Cv); const Real Ri = a * (E - B); @@ -172,8 +172,8 @@ TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { e = eos_d.InternalEnergyFromDensityTemperature(dens, T) * dens; const Real dEg = e - e0; Real a = chat * dt * - (opac_d.RosselandMeanAbsorptionCoefficient(dens, T) + - scat_d.RosselandMeanTotalScatteringCoefficient(dens, T)); + (opac_d.AbsorptionCoefficient(dens, T, 0) + + scat_d.ScatteringCoefficient(dens, T, 0)); std::array dF{-a / (1. + a) * Fr0[0], -a / (1. + a) * Fr0[1], -a / (1. + a) * Fr0[2]}; const Real icc = -1. / (c * chat * dens); @@ -366,9 +366,9 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const Real Cv = dens * eos_d.SpecificHeatFromDensityTemperature(dens, T); const Real fleck = FleckFactor(arad, T, Cv); - const Real sigp = chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); - const Real sigs = - chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); + const Real sigp = + chat * dt * opac_d.PlanckGroupAbsorptionCoefficient(dens, T, 0); + const Real sigs = chat * dt * scat_d.ScatteringCoefficient(dens, T, 0); const Real sigf = sigp + sigs; const Real ca = g * (sigf - g2 * sigs * (1. + bdbdp)); @@ -409,10 +409,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real eg = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / eref; dEg = eg - eg0; - const Real sigp = - chat * dt * opac_d.RosselandMeanAbsorptionCoefficient(dens, T); - const Real sigs = - chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); + const Real sigp = chat * dt * opac_d.AbsorptionCoefficient(dens, T, 0); + const Real sigs = chat * dt * scat_d.ScatteringCoefficient(dens, T, 0); const Real sigf = sigp + sigs; const Real a = g * sigf; diff --git a/src/radiation/moments/moments.cpp b/src/radiation/moments/moments.cpp index d6fb6ca4..20ccb046 100644 --- a/src/radiation/moments/moments.cpp +++ b/src/radiation/moments/moments.cpp @@ -427,7 +427,7 @@ void InitMesh(parthenon::Mesh *pmesh) { if (multi_d) dx_min = std::min(dx_min, dx[1]); if (three_d) dx_min = std::min(dx_min, dx[2]); const Real tau = - std::min(1.0, dx_min * opac_d.RosselandMeanAbsorptionCoefficient(rho, T)); + std::min(1.0, dx_min * opac_d.AbsorptionCoefficient(rho, T, 0)); const Real Erad = tau * arad * SQR(SQR(T)); vmesh(b, rad::cons::energy(0), k, j, i) = Erad; vmesh(b, rad::prim::energy(0), k, j, i) = Erad; diff --git a/src/radiation/radiation.cpp b/src/radiation/radiation.cpp index d127ba9c..6873885e 100644 --- a/src/radiation/radiation.cpp +++ b/src/radiation/radiation.cpp @@ -94,7 +94,7 @@ std::shared_ptr Initialize(ParameterInput *pin, // incorporate frequency type for gas opacity initialization params.Add("frequency_type", frequency_type); - // Initialize gas opacity + // Initialize gas opacity models Gas::InitGasOpacity(pin, units, params); // Enroll in tstart/tstop machinery @@ -139,8 +139,8 @@ TaskStatus SetOpacities(MeshData *md) { 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); + aa = opacity_d.AbsorptionCoefficient(rho, temp, 0); + ss = scattering_d.ScatteringCoefficient(rho, temp, 0); }); return TaskStatus::complete; diff --git a/src/utils/opacity/opacity.hpp b/src/utils/opacity/opacity.hpp index d2049c0e..4bc63452 100644 --- a/src/utils/opacity/opacity.hpp +++ b/src/utils/opacity/opacity.hpp @@ -39,7 +39,7 @@ using MeanOpacity = // Reduced variant for mean scattering opacities using MeanScattering = - singularity::photons::MeanNonCGSUnitsS; + singularity::photons::MeanNonCGSUnitsS; } // namespace ArtemisUtils From cf946dbc1d742d209f85621a97348d28f3b040e2 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 4 Aug 2026 11:00:01 -0600 Subject: [PATCH 7/9] Attempt fix for non-cgs regression. + Scale opacity bounds. + Bump jaybenne. --- external/jaybenne | 2 +- src/radiation/gas_opacity.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/external/jaybenne b/external/jaybenne index 69b27fdf..cc83bf1d 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit 69b27fdf2d02448ee844df51bc60a292c7e3d57a +Subproject commit cc83bf1d8bd4728cd2444eac1fc8b4fb95b773f1 diff --git a/src/radiation/gas_opacity.cpp b/src/radiation/gas_opacity.cpp index 73852422..38119543 100644 --- a/src/radiation/gas_opacity.cpp +++ b/src/radiation/gas_opacity.cpp @@ -28,7 +28,8 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, const auto frequency_type = params.Get("frequency_type"); // TODO: read this (and consolidate gray/multigroup modes?) - const std::vector gray_bounds = {1.e12, 3.e20}; + // hard-coded numbers in Hz + const std::vector gray_bounds = {time * 1.e12, time * 3.e20}; const int NG = static_cast(gray_bounds.size()) - 1; // Absorption opacity model From a46d96e675d90ad62e67dac7c5d7ffaa4a96784c Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Fri, 7 Aug 2026 10:59:49 -0600 Subject: [PATCH 8/9] Update singularity-opac, jaybenne, extend gray frequency bounds. --- external/jaybenne | 2 +- external/singularity-opac | 2 +- src/radiation/gas_opacity.cpp | 3 +-- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/external/jaybenne b/external/jaybenne index cc83bf1d..29b0b154 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit cc83bf1d8bd4728cd2444eac1fc8b4fb95b773f1 +Subproject commit 29b0b1544a76115806fea35797088d87089151ce diff --git a/external/singularity-opac b/external/singularity-opac index b97515f6..d57886d1 160000 --- a/external/singularity-opac +++ b/external/singularity-opac @@ -1 +1 @@ -Subproject commit b97515f695f96fabeca28852b5483bc0cd49528e +Subproject commit d57886d10f4068de119502b93a30af7114a542de diff --git a/src/radiation/gas_opacity.cpp b/src/radiation/gas_opacity.cpp index 38119543..891e83c9 100644 --- a/src/radiation/gas_opacity.cpp +++ b/src/radiation/gas_opacity.cpp @@ -28,8 +28,7 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, const auto frequency_type = params.Get("frequency_type"); // TODO: read this (and consolidate gray/multigroup modes?) - // hard-coded numbers in Hz - const std::vector gray_bounds = {time * 1.e12, time * 3.e20}; + const std::vector gray_bounds = {0.0, std::numeric_limits::infinity()}; const int NG = static_cast(gray_bounds.size()) - 1; // Absorption opacity model From 5f7e69cb7c0da89493fbc9ca39834b86d77a819d Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Fri, 14 Aug 2026 17:25:17 -0600 Subject: [PATCH 9/9] Use MeanOpacity and MeanScattering in multigroup mode. + Remove stored monochromatic opacity models. + Parse frequency grid from radblock_name (input: imc or moments). + Add use_opac_grps to use pre-tabulated absorption frequency grid. - Applied as both radiation and scattering grid. + Remove imc's monochromatic overload of Initialization. + Update jaybenne and singularity-opac. --- external/jaybenne | 2 +- external/singularity-opac | 2 +- src/artemis.cpp | 22 +++------- src/radiation/gas_opacity.cpp | 80 +++++++++++++++++------------------ src/radiation/gas_opacity.hpp | 4 +- src/radiation/radiation.cpp | 13 +++--- 6 files changed, 56 insertions(+), 67 deletions(-) diff --git a/external/jaybenne b/external/jaybenne index 29b0b154..9f80d223 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit 29b0b1544a76115806fea35797088d87089151ce +Subproject commit 9f80d223a811626cf2d9b4194ea5cc7f1ca07ab5 diff --git a/external/singularity-opac b/external/singularity-opac index d57886d1..7d72f0f1 160000 --- a/external/singularity-opac +++ b/external/singularity-opac @@ -1 +1 @@ -Subproject commit d57886d10f4068de119502b93a30af7114a542de +Subproject commit 7d72f0f16376e4e710911cb56ae14b49f9072848 diff --git a/src/artemis.cpp b/src/artemis.cpp index 01eb09c0..fafbef76 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -189,23 +189,11 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { // Select between Jaybenne IMC or Moments if (do_imc) { auto eos_h = packages.Get("gas")->Param("eos_h"); - const auto frequency_type = - packages.Get("radiation")->Param("frequency_type"); - if (frequency_type == FrequencyType::gray) { - auto opacity_h = packages.Get("radiation")->Param("opacity_h"); - auto scattering_h = - packages.Get("radiation")->Param("scattering_h"); - packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h, - "radiation/imc")); - } else { - PARTHENON_REQUIRE(frequency_type == FrequencyType::multigroup, - "Invalid frequency_type!"); - auto opacity_h = packages.Get("radiation")->Param("mg_opacity_h"); - auto scattering_h = - packages.Get("radiation")->Param("mg_scattering_h"); - packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h, - "radiation/imc")); - } + auto opacity_h = packages.Get("radiation")->Param("opacity_h"); + auto scattering_h = + packages.Get("radiation")->Param("scattering_h"); + 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) { diff --git a/src/radiation/gas_opacity.cpp b/src/radiation/gas_opacity.cpp index 891e83c9..8e27922a 100644 --- a/src/radiation/gas_opacity.cpp +++ b/src/radiation/gas_opacity.cpp @@ -14,8 +14,8 @@ #include "utils/opacity/opacity.hpp" namespace Gas { -void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, - Params ¶ms) { +void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, Params ¶ms, + const std::string &radblock_name) { using namespace singularity::photons; // Opacity models @@ -27,16 +27,36 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, // Get frequency type (it should already be set in radiation Initialization) const auto frequency_type = params.Get("frequency_type"); - // TODO: read this (and consolidate gray/multigroup modes?) - const std::vector gray_bounds = {0.0, std::numeric_limits::infinity()}; - const int NG = static_cast(gray_bounds.size()) - 1; - // Absorption opacity model std::string opacity_model_name = pin->GetOrAddString("gas/opacity/absorption", "opacity_model", "constant"); + // check if using groups from parsed-in opacity table + const bool use_opac_grps = + pin->GetOrAddBoolean(radblock_name, "use_opac_groups", false); + PARTHENON_REQUIRE( + use_opac_grps ? opacity_model_name == "table" : true, + "Opacity group bounds can only be used with opacity_model_name=table!"); + + // set opacity group bounds: gray goes from 0 to infty + std::vector opac_grp_bnds = {0.0, std::numeric_limits::infinity()}; + // reset to parsed input (for non-tabular opacity models) for multigroup + if (frequency_type == FrequencyType::multigroup && !use_opac_grps) { + const Real numin = pin->GetReal(radblock_name, "numin"); // in Hz + const Real numax = pin->GetReal(radblock_name, "numax"); // in Hz + const int n_nubins = pin->GetInteger(radblock_name, "n_nubins"); + // reset opacity group bounds + // NOTE: these are group edges, not interior points + opac_grp_bnds.assign(n_nubins + 1, 0.0); + // assume uniform log-spacing, grid is midpoints in log-space + const Real dlnu = (std::log(numax) - std::log(numin)) / n_nubins; + for (int n = 0; n < n_nubins + 1; ++n) { + opac_grp_bnds[n] = numin * std::exp(n * dlnu); + } + } + int NG = static_cast(opac_grp_bnds.size()) - 1; + // Mean absorption opacity (either read from table or uses model - ArtemisUtils::Opacity mg_opacity; ArtemisUtils::MeanOpacity opacity; if (opacity_model_name == "table") { PARTHENON_REQUIRE(frequency_type == FrequencyType::gray, @@ -61,26 +81,18 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, opacity = singularity::photons::MeanNonCGSUnits( singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, - lTMin_a, lTMax_a, NT_a, gray_bounds, + lTMin_a, lTMax_a, NT_a, opac_grp_bnds, NG), time, mass, length, temp); - if (frequency_type == FrequencyType::multigroup) { - mg_opacity = singularity::photons::NonCGSUnits( - std::move(model), time, mass, length, temp); - } } else if (opacity_model_name == "constant") { const Real kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "kappa_a", 0.0); auto model = Gray(kappa_a); opacity = singularity::photons::MeanNonCGSUnits( singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, - lTMin_a, lTMax_a, NT_a, gray_bounds, + lTMin_a, lTMax_a, NT_a, opac_grp_bnds, NG), time, mass, length, temp); - if (frequency_type == FrequencyType::multigroup) { - mg_opacity = singularity::photons::NonCGSUnits( - std::move(model), time, mass, length, temp); - } } else if (opacity_model_name == "powerlaw") { const Real coef_kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "coef_kappa_a", 0.0); @@ -90,13 +102,9 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, opacity = singularity::photons::MeanNonCGSUnits( singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a, - lTMin_a, lTMax_a, NT_a, gray_bounds, + lTMin_a, lTMax_a, NT_a, opac_grp_bnds, NG), time, mass, length, temp); - if (frequency_type == FrequencyType::multigroup) { - mg_opacity = singularity::photons::NonCGSUnits( - std::move(model), time, mass, length, temp); - } } else { PARTHENON_FAIL("Opacity model not recognized!"); } @@ -104,15 +112,18 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, params.Add("opacity_h", opacity); params.Add("opacity_d", opacity.GetOnDevice()); - if (frequency_type == FrequencyType::multigroup) { - params.Add("mg_opacity_h", mg_opacity); - params.Add("mg_opacity_d", mg_opacity.GetOnDevice()); - } // Scattering opacity model std::string scattering_model_name = pin->GetOrAddString("gas/opacity/scattering", "scattering_model", "none"); + // ensure analytic scattering uses tabular absorption group bounds + // TODO: table scattering opacity + if (use_opac_grps && opacity_model_name == "table") { + opac_grp_bnds = opacity.GetGroupBounds(); + NG = opacity.ngroups(); + } + // Instantiate mean scattering opacity object (i.e., table) const Real lRhoMin_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMin", -1.0); const Real lRhoMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMax", 1.0); @@ -121,43 +132,30 @@ void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, const Real lTMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMax", 1.0); const int NT_s = pin->GetOrAddInteger("gas/opacity/scattering", "NT", 2); - ArtemisUtils::Scattering mg_scattering; ArtemisUtils::MeanScattering scattering; if (scattering_model_name == "none") { auto smodel = GrayS(0.0, 1.0); scattering = singularity::photons::MeanNonCGSUnitsS( singularity::photons::MeanSOpacityBase(smodel, lRhoMin_s, lRhoMax_s, NRho_s, - lTMin_s, lTMax_s, NT_s, gray_bounds, + lTMin_s, lTMax_s, NT_s, opac_grp_bnds, NG), time, mass, length, temp); - if (frequency_type == FrequencyType::multigroup) { - mg_scattering = singularity::photons::NonCGSUnitsS( - std::move(smodel), time, mass, length, temp); - } } else if (scattering_model_name == "constant") { const Real kappa_s = pin->GetOrAddReal("gas/opacity/scattering", "kappa_s", 0.0); auto smodel = GrayS(kappa_s, 1.0); scattering = singularity::photons::MeanNonCGSUnitsS( singularity::photons::MeanSOpacityBase(smodel, lRhoMin_s, lRhoMax_s, NRho_s, - lTMin_s, lTMax_s, NT_s, gray_bounds, + lTMin_s, lTMax_s, NT_s, opac_grp_bnds, NG), time, mass, length, temp); - if (frequency_type == FrequencyType::multigroup) { - mg_scattering = singularity::photons::NonCGSUnitsS( - std::move(smodel), time, mass, length, temp); - } } else { PARTHENON_FAIL("Scattering model not recognized!"); } params.Add("scattering_h", scattering); params.Add("scattering_d", scattering.GetOnDevice()); - if (frequency_type == FrequencyType::multigroup) { - params.Add("mg_scattering_h", mg_scattering); - params.Add("mg_scattering_d", mg_scattering.GetOnDevice()); - } } } // namespace Gas diff --git a/src/radiation/gas_opacity.hpp b/src/radiation/gas_opacity.hpp index d4a24b77..b84ee142 100644 --- a/src/radiation/gas_opacity.hpp +++ b/src/radiation/gas_opacity.hpp @@ -21,8 +21,8 @@ #include namespace Gas { -void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, - Params ¶ms); +void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units, Params ¶ms, + const std::string &radblock_name); } // namespace Gas #endif // RADIATION_GAS_OPACITY_HPP_ diff --git a/src/radiation/radiation.cpp b/src/radiation/radiation.cpp index 6873885e..5407ca84 100644 --- a/src/radiation/radiation.cpp +++ b/src/radiation/radiation.cpp @@ -79,8 +79,6 @@ std::shared_ptr Initialize(ParameterInput *pin, for (int n = 0; n < nspecies; ++n) fluidids.push_back(n); - // Control field for sparse gas fields - // Absorption and scattering opacity Metadata m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy, MetadataRadiation, MetadataOperatorSplit}); @@ -94,12 +92,13 @@ std::shared_ptr Initialize(ParameterInput *pin, // incorporate frequency type for gas opacity initialization params.Add("frequency_type", frequency_type); + std::string radblock_name = (do_imc ? "radiation/imc" : "radiation/moment"); + // Initialize gas opacity models - Gas::InitGasOpacity(pin, units, params); + Gas::InitGasOpacity(pin, units, params, radblock_name); // Enroll in tstart/tstop machinery - ArtemisUtils::AddPackageTimeParams( - params, (do_imc) ? "radiation/imc" : "radiation/moment", pin); + ArtemisUtils::AddPackageTimeParams(params, radblock_name, pin); return radiation; } @@ -114,6 +113,10 @@ TaskStatus SetOpacities(MeshData *md) { auto &gas_pkg = pm->packages.Get("gas"); auto &rad_pkg = pm->packages.Get("radiation"); + // do not populate gray fields in MG + const auto frequency_type = rad_pkg->Param("frequency_type"); + if (frequency_type == FrequencyType::multigroup) return TaskStatus::complete; + EOS eos_d = gas_pkg->template Param("eos_d"); MeanOpacity opacity_d = rad_pkg->template Param("opacity_d"); MeanScattering scattering_d = rad_pkg->template Param("scattering_d");