diff --git a/external/jaybenne b/external/jaybenne index 1320212e..9f80d223 160000 --- a/external/jaybenne +++ b/external/jaybenne @@ -1 +1 @@ -Subproject commit 1320212e6de773a9a964e680b59e488d0d193baf +Subproject commit 9f80d223a811626cf2d9b4194ea5cc7f1ca07ab5 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..7d72f0f1 160000 --- a/external/singularity-opac +++ b/external/singularity-opac @@ -1 +1 @@ -Subproject commit cdd365ffd56ccdab1ae4b4410d276d0a46a4a314 +Subproject commit 7d72f0f16376e4e710911cb56ae14b49f9072848 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 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..fafbef76 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -185,12 +185,13 @@ 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"); + 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, 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..8e27922a --- /dev/null +++ b/src/radiation/gas_opacity.cpp @@ -0,0 +1,161 @@ +//======================================================================================== +// (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, + const std::string &radblock_name) { + 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"); + + // 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::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, opac_grp_bnds, + NG), + 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, opac_grp_bnds, + NG), + 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, opac_grp_bnds, + NG), + 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"); + + // 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); + 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::MeanSOpacityBase(smodel, lRhoMin_s, lRhoMax_s, NRho_s, + lTMin_s, lTMax_s, NT_s, opac_grp_bnds, + NG), + 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, opac_grp_bnds, + NG), + 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..b84ee142 --- /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, + const std::string &radblock_name); +} // 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..12540431 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"); @@ -142,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); @@ -171,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); @@ -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"); @@ -364,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)); @@ -407,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 4cb1ca85..20ccb046 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( @@ -426,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 9a10d34d..5407ca84 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); @@ -65,18 +79,26 @@ Initialize(ParameterInput *pin, ArtemisUtils::Constants &constants, const bool d 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}); 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); + + std::string radblock_name = (do_imc ? "radiation/imc" : "radiation/moment"); + + // Initialize gas opacity models + 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; } @@ -89,10 +111,15 @@ 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"); + + // 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 = 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 @@ -115,8 +142,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/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..4bc63452 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 @@ -38,7 +39,7 @@ using MeanOpacity = // Reduced variant for mean scattering opacities using MeanScattering = - singularity::photons::MeanNonCGSUnitsS; + singularity::photons::MeanNonCGSUnitsS; } // namespace ArtemisUtils