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/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/CMakeLists.txt b/src/CMakeLists.txt index 51386b94..7199ce3b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -77,6 +77,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..01eb09c0 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -185,14 +185,27 @@ 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"); - 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("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")); + } PARTHENON_REQUIRE(coords == Coordinates::cartesian, "Jaybenne currently supports only Cartesian coordinates!"); } else if (do_moment) { diff --git a/src/artemis.hpp b/src/artemis.hpp index 20e13e93..5e593429 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -185,6 +185,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..fe043841 --- /dev/null +++ b/src/radiation/gas_opacity.cpp @@ -0,0 +1,154 @@ +//======================================================================================== +// (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; + + // 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 = + 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); + 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), + 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); + 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); + 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 = + 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::Scattering mg_scattering; + 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); + 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::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 new file mode 100644 index 00000000..d4a24b77 --- /dev/null +++ b/src/radiation/gas_opacity.hpp @@ -0,0 +1,28 @@ +//======================================================================================== +// (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/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 9a10d34d..d127ba9c 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,19 @@ Initialize(ParameterInput *pin, ArtemisUtils::Constants &constants, const bool d params.Add("chat", light); } + // 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 == "multigroup") { + frequency_type = FrequencyType::multigroup; + } else { + 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); params.Add("nspecies", nspecies); @@ -72,8 +86,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); @@ -89,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 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