From 41c4a4527c37ba3e55aeeb10b8779723e70f39b5 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 18 Dec 2024 08:34:24 -0700 Subject: [PATCH 01/16] ryan's mean singularity-opac --- external/singularity-opac | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/singularity-opac b/external/singularity-opac index b4542ed..b61be52 160000 --- a/external/singularity-opac +++ b/external/singularity-opac @@ -1 +1 @@ -Subproject commit b4542ed7cd618e6ad0d9331205f7ed3d0c2e3a0d +Subproject commit b61be5226fee505d05ffeb6419408740792b51c7 From acb2cfdf04ee3be2f65332a4f97d70d0ff38a6d6 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 18 Dec 2024 08:46:53 -0700 Subject: [PATCH 02/16] Multiple initialize commands --- src/jaybenne/jaybenne.cpp | 66 +++++++++++++++++++++++++++++---------- src/jaybenne/jaybenne.hpp | 5 +++ src/mcblock/opacity.hpp | 6 ++++ 3 files changed, 61 insertions(+), 16 deletions(-) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 0a167cf..a3345c5 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -151,12 +151,11 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) { } //---------------------------------------------------------------------------------------- -//! \fn StateDescriptor Jaybenne::Initialize +//! \fn StateDescriptor Jaybenne::Initialize_impl //! \brief Initialize the Jaybenne physics package. This function defines and sets the //! parameters associated with Jaybenne, and enrolls the data variables associated with -//! this physics package. -std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacity, - Scattering &scattering, EOS &eos) { +//! this physics package, for everything not related to frequency discretization. +std::shared_ptr Initialize_impl(ParameterInput *pin, EOS &eos) { auto pkg = std::make_shared("jaybenne"); // Total number of particles @@ -172,12 +171,6 @@ std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacit "Minimum allowable swarm occupancy must be >= 0 and less than 1"); pkg->AddParam<>("min_swarm_occupancy", min_swarm_occupancy); - // Frequency range - Real numin = pin->GetOrAddReal("jaybenne", "numin", std::numeric_limits::min()); - pkg->AddParam<>("numin", numin); - Real numax = pin->GetOrAddReal("jaybenne", "numax", std::numeric_limits::max()); - pkg->AddParam<>("numax", numax); - // Physical constants const auto units = opacity.GetRuntimePhysicalConstants(); pkg->AddParam<>("speed_of_light", units.c); @@ -226,12 +219,6 @@ std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacit // Equation of state model pkg->AddParam<>("eos_d", eos.GetOnDevice()); - // Opacity model - pkg->AddParam<>("opacity_d", opacity.GetOnDevice()); - - // Scattering model - pkg->AddParam<>("scattering_d", scattering.GetOnDevice()); - // Swarm and swarm variables Metadata swarm_metadata({Metadata::Provides, Metadata::None}); pkg->AddSwarm(photons_swarm_name, swarm_metadata); @@ -265,6 +252,53 @@ std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacit return pkg; } +//---------------------------------------------------------------------------------------- +//! \fn StateDescriptor Jaybenne::Initialize +//! \brief Initialize the Jaybenne physics package. This function defines and sets the +//! parameters associated with Jaybenne, and enrolls the data variables associated with +//! this physics package, specific to the gray frequency discretization. +std::shared_ptr Initialize(ParameterInput *pin, MeanOpacity &mopacity, + MeanScattering &mscattering, EOS &eos) { + auto pkg = Initialize_impl(pin, eos); + + pkg->AddParam<>("frequency_discretization", FrequencyDiscretization::multigroup); + + // Opacity model + pkg->AddParam<>("mopacity_d", mopacity.GetOnDevice()); + + // Scattering model + pkg->AddParam<>("mscattering_d", mscattering.GetOnDevice()); + + return pkg; +} + +//---------------------------------------------------------------------------------------- +//! \fn StateDescriptor Jaybenne::Initialize +//! \brief Initialize the Jaybenne physics package. This function defines and sets the +//! parameters associated with Jaybenne, and enrolls the data variables associated with +//! this physics package, specific to the multigroup frequency discretization. +std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacity, + Scattering &scattering, EOS &eos) { + + auto pkg = Initialize_impl(pin, eos); + + // Frequency range + Real numin = pin->GetOrAddReal("jaybenne", "numin", std::numeric_limits::min()); + pkg->AddParam<>("numin", numin); + Real numax = pin->GetOrAddReal("jaybenne", "numax", std::numeric_limits::max()); + pkg->AddParam<>("numax", numax); + + pkg->AddParam<>("frequency_discretization", FrequencyDiscretization::multigroup); + + // Opacity model + pkg->AddParam<>("opacity_d", opacity.GetOnDevice()); + + // Scattering model + pkg->AddParam<>("scattering_d", scattering.GetOnDevice()); + + return pkg; +} + //---------------------------------------------------------------------------------------- //! \fn Real Jaybenne::EstimateTimestepMesh //! \brief Compute radiation timestep diff --git a/src/jaybenne/jaybenne.hpp b/src/jaybenne/jaybenne.hpp index 0a5f88f..203a9e2 100644 --- a/src/jaybenne/jaybenne.hpp +++ b/src/jaybenne/jaybenne.hpp @@ -50,10 +50,15 @@ namespace jaybenne { std::shared_ptr Initialize(parthenon::ParameterInput *pin, Opacity &opacity, Scattering &scattering, EOS &eos); +std::shared_ptr Initialize(parthenon::ParameterInput *pin, + MeanOpacity &mopacity, + MeanScattering &mscattering, + EOS &eos); // Model enums enum class SourceStrategy { uniform, energy }; enum class SourceType { thermal, emission }; +enum class FrequencyDiscretization { gray, multigroup }; // Tasks TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real dt); diff --git a/src/mcblock/opacity.hpp b/src/mcblock/opacity.hpp index 580d992..e32767a 100644 --- a/src/mcblock/opacity.hpp +++ b/src/mcblock/opacity.hpp @@ -14,6 +14,8 @@ #define MCBLOCK_OPACITY_HPP_ // Singularity-opac includes +#include +#include #include #include @@ -24,11 +26,15 @@ using Opacity = singularity::photons::impl::Variant< singularity::photons::NonCGSUnits, singularity::photons::NonCGSUnits>; +using MeanOpacity = singularity::photons::MeanOpacity; + // Reduced scattering variant just for jaybenne using Scattering = singularity::photons::impl::S_Variant< singularity::photons::NonCGSUnitsS, singularity::photons::NonCGSUnitsS>; +using MeanScattering = singularity::photons::MeanSOpacity; + } // namespace mcblock #endif // MCBLOCK_OPACITY_HPP_ From 86a9111cbda6a79668b29ac5e68367190f7ad1de Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 18 Dec 2024 08:54:51 -0700 Subject: [PATCH 03/16] fix some issues --- src/CMakeLists.txt | 6 ++++++ src/jaybenne/jaybenne.cpp | 13 +++++++++---- src/jaybenne/jaybenne_config.hpp.in | 2 ++ 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d7cb917..e45dd48 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -31,6 +31,12 @@ if(JAYBENNE_STANDALONE_MODE) set(JAYBENNE_SCATTERING_OPACITY_TYPE "mcblock::Scattering" CACHE STRING "Namespace for scattering opacity Singularity type" FORCE) + set(JAYBENNE_ABSORPTION_MEAN_OPACITY_TYPE + "mcblock::MeanOpacity" + CACHE STRING "Namespace for absorption mean opacity Singularity type" FORCE) + set(JAYBENNE_SCATTERING_MEAN_OPACITY_TYPE + "mcblock::MeanScattering" + CACHE STRING "Namespace for scattering mean opacity Singularity type" FORCE) set(JAYBENNE_HOST_VARIABLE_HEADER "../mcblock/mcblock_variables.hpp" CACHE STRING "Path to host code variables header file" FORCE) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index a3345c5..b7bc31a 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -155,7 +155,9 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) { //! \brief Initialize the Jaybenne physics package. This function defines and sets the //! parameters associated with Jaybenne, and enrolls the data variables associated with //! this physics package, for everything not related to frequency discretization. -std::shared_ptr Initialize_impl(ParameterInput *pin, EOS &eos) { +std::shared_ptr +Initialize_impl(ParameterInput *pin, EOS &eos, + singularity::RuntimePhysicalConstants units) { auto pkg = std::make_shared("jaybenne"); // Total number of particles @@ -172,7 +174,7 @@ std::shared_ptr Initialize_impl(ParameterInput *pin, EOS &eos) pkg->AddParam<>("min_swarm_occupancy", min_swarm_occupancy); // Physical constants - const auto units = opacity.GetRuntimePhysicalConstants(); + // const auto units = opacity.GetRuntimePhysicalConstants(); pkg->AddParam<>("speed_of_light", units.c); pkg->AddParam<>("stefan_boltzmann", units.sb); @@ -259,7 +261,9 @@ std::shared_ptr Initialize_impl(ParameterInput *pin, EOS &eos) //! this physics package, specific to the gray frequency discretization. std::shared_ptr Initialize(ParameterInput *pin, MeanOpacity &mopacity, MeanScattering &mscattering, EOS &eos) { - auto pkg = Initialize_impl(pin, eos); + const auto units = mopacity.GetRuntimePhysicalConstants(); + + auto pkg = Initialize_impl(pin, eos, units); pkg->AddParam<>("frequency_discretization", FrequencyDiscretization::multigroup); @@ -279,8 +283,9 @@ std::shared_ptr Initialize(ParameterInput *pin, MeanOpacity &mo //! this physics package, specific to the multigroup frequency discretization. std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacity, Scattering &scattering, EOS &eos) { + const auto units = opacity.GetRuntimePhysicalConstants(); - auto pkg = Initialize_impl(pin, eos); + auto pkg = Initialize_impl(pin, eos, units); // Frequency range Real numin = pin->GetOrAddReal("jaybenne", "numin", std::numeric_limits::min()); diff --git a/src/jaybenne/jaybenne_config.hpp.in b/src/jaybenne/jaybenne_config.hpp.in index 6eb4518..662fdfa 100644 --- a/src/jaybenne/jaybenne_config.hpp.in +++ b/src/jaybenne/jaybenne_config.hpp.in @@ -22,6 +22,8 @@ using @JAYBENNE_EOS_TYPE@; using @JAYBENNE_ABSORPTION_OPACITY_TYPE@; using @JAYBENNE_SCATTERING_OPACITY_TYPE@; +using @JAYBENNE_ABSORPTION_MEAN_OPACITY_TYPE@; +using @JAYBENNE_SCATTERING_MEAN_OPACITY_TYPE@; #include "@JAYBENNE_HOST_VARIABLE_HEADER@" From a72bded234b682e2adfd0f6af769ea1b390b15e9 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 18 Dec 2024 11:07:53 -0700 Subject: [PATCH 04/16] Sample MG CDF --- src/jaybenne/jaybenne.cpp | 52 ++++++++++++---- src/jaybenne/jaybenne.hpp | 4 +- src/jaybenne/jaybenne_variables.hpp | 1 + src/jaybenne/sourcing.cpp | 93 +++++++++++++++++++++++++---- 4 files changed, 125 insertions(+), 25 deletions(-) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index b7bc31a..1a39c55 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -71,6 +71,7 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) { auto &jb_pkg = pmesh->packages.Get("jaybenne"); const auto &max_transport_iterations = jb_pkg->Param("max_transport_iterations"); const bool &use_ddmc = jb_pkg->Param("use_ddmc"); + const auto &fd = jb_pkg->Param("frequency_type"); // MeshData subsets auto ddmc_field_names = std::vector{fj::ddmc_face_prob::name()}; @@ -102,9 +103,20 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) { // prepare for iterative transport loop auto derived = tl.AddTask(none, UpdateDerivedTransportFields, base.get(), dt); - auto source = tl.AddTask( - derived, jaybenne::SourcePhotons, jaybenne::SourceType::emission>, - base.get(), t_start, dt); + TaskID source = derived; + if (fd == FrequencyType::gray) { + auto source = tl.AddTask( + derived, + jaybenne::SourcePhotons, jaybenne::SourceType::emission, + FrequencyType::gray>, + base.get(), t_start, dt); + } else if (fd == FrequencyType::multigroup) { + auto source = tl.AddTask( + derived, + jaybenne::SourcePhotons, jaybenne::SourceType::emission, + FrequencyType::multigroup>, + base.get(), t_start, dt); + } auto bcs = use_ddmc ? parthenon::AddBoundaryExchangeTasks(source, tl, md_ddmc, pmesh->multilevel) : source; @@ -177,6 +189,7 @@ Initialize_impl(ParameterInput *pin, EOS &eos, // const auto units = opacity.GetRuntimePhysicalConstants(); pkg->AddParam<>("speed_of_light", units.c); pkg->AddParam<>("stefan_boltzmann", units.sb); + pkg->AddParam<>("planck_constant", units.h); // RNG bool unique_rank_seeds = pin->GetOrAddBoolean("jaybenne", "unique_rank_seeds", true); @@ -265,7 +278,7 @@ std::shared_ptr Initialize(ParameterInput *pin, MeanOpacity &mo auto pkg = Initialize_impl(pin, eos, units); - pkg->AddParam<>("frequency_discretization", FrequencyDiscretization::multigroup); + pkg->AddParam<>("frequency_type", FrequencyType::multigroup); // Opacity model pkg->AddParam<>("mopacity_d", mopacity.GetOnDevice()); @@ -287,13 +300,20 @@ std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacit auto pkg = Initialize_impl(pin, eos, units); - // Frequency range - Real numin = pin->GetOrAddReal("jaybenne", "numin", std::numeric_limits::min()); - pkg->AddParam<>("numin", numin); - Real numax = pin->GetOrAddReal("jaybenne", "numax", std::numeric_limits::max()); - pkg->AddParam<>("numax", numax); + // Frequency discretization + auto time = units.time; + Real numin = pin->GetReal("jaybenne", "numin"); // in Hz + pkg->AddParam<>("numin", numin * time); // in code units + Real numax = pin->GetReal("jaybenne", "numax"); // in Hz + pkg->AddParam<>("numax", numax * time); // in code units + int n_nubins = pin->GetInteger("jaybenne", "n_nubins"); + pkg->AddParam<>("n_nubins", n_nubins); + + pkg->AddParam<>("frequency_type", FrequencyType::multigroup); - pkg->AddParam<>("frequency_discretization", FrequencyDiscretization::multigroup); + // Emission CDF + Metadata m_onecopy({Metadata::Cell, Metadata::OneCopy}, std::vector({n_nubins})); + pkg->AddField(field::jaybenne::emission_cdf::name(), m_onecopy); // Opacity model pkg->AddParam<>("opacity_d", opacity.GetOnDevice()); @@ -595,9 +615,17 @@ TaskStatus EvaluateRadiationEnergy(T *md) { //! \brief Initialize radiation based on material temperature and either thermal or //! zero initial radiation. void InitializeRadiation(MeshBlockData *mbd, const bool is_thermal) { + auto &jb_pkg = mbd->GetBlockPointer()->packages.Get("jaybenne"); + const auto &fd = jb_pkg->Param("frequency_type"); + if (is_thermal) { - jaybenne::SourcePhotons, jaybenne::SourceType::thermal>(mbd, 0.0, - 0.0); + if (fd == FrequencyType::gray) { + jaybenne::SourcePhotons, jaybenne::SourceType::thermal, + FrequencyType::gray>(mbd, 0.0, 0.0); + } else if (fd == FrequencyType::multigroup) { + jaybenne::SourcePhotons, jaybenne::SourceType::thermal, + FrequencyType::multigroup>(mbd, 0.0, 0.0); + } } // Call this so radiation field variables are properly initialized for outputs diff --git a/src/jaybenne/jaybenne.hpp b/src/jaybenne/jaybenne.hpp index 203a9e2..60dd5ab 100644 --- a/src/jaybenne/jaybenne.hpp +++ b/src/jaybenne/jaybenne.hpp @@ -58,14 +58,14 @@ std::shared_ptr Initialize(parthenon::ParameterInput // Model enums enum class SourceStrategy { uniform, energy }; enum class SourceType { thermal, emission }; -enum class FrequencyDiscretization { gray, multigroup }; +enum class FrequencyType { gray, multigroup }; // Tasks TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real dt); TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const Real dt); TaskStatus SampleDDMCBlockFace(MeshData *md); TaskStatus CheckCompletion(MeshData *md, const Real t_end); -template +template TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt); TaskStatus DefragParticles(MeshBlock *pmb); TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt); diff --git a/src/jaybenne/jaybenne_variables.hpp b/src/jaybenne/jaybenne_variables.hpp index 644ebe2..fe5e100 100644 --- a/src/jaybenne/jaybenne_variables.hpp +++ b/src/jaybenne/jaybenne_variables.hpp @@ -37,6 +37,7 @@ JAYBENNE_FIELD_VARIABLE(field.jaybenne, fleck_factor); JAYBENNE_FIELD_VARIABLE(field.jaybenne, ddmc_face_prob); JAYBENNE_FIELD_VARIABLE(field.jaybenne, source_ew_per_cell); JAYBENNE_FIELD_VARIABLE(field.jaybenne, source_num_per_cell); +JAYBENNE_FIELD_VARIABLE(field.jaybenne, emission_cdf); JAYBENNE_FIELD_VARIABLE(field.jaybenne, energy_delta); namespace host { typedef HOST_DENSITY density; diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 2b82f04..a07b597 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -22,7 +22,7 @@ namespace jaybenne { //! \fn TaskStatus SourcePhotons //! \brief Create photons, either during initialization or during a timestep. //! TODO(BRR) modify interface so we don't need t_start, dt for initialization -template +template TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { namespace fj = field::jaybenne; namespace fjh = field::jaybenne::host; @@ -31,8 +31,21 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; auto &jb_pkg = pm->packages.Get("jaybenne"); + const Real h = jb_pkg->template Param("planck_constant"); auto &eos = jb_pkg->template Param("eos_d"); - auto &opacity = jb_pkg->template Param("opacity_d"); + const MeanOpacity *mopacity = nullptr; + const Opacity *opacity = nullptr; + int n_nubins = -1; + Real numin = -1.; + Real numax = -1.; + if constexpr (FT == FrequencyType::gray) { + mopacity = &(jb_pkg->template Param("mopacity_d")); + } else if constexpr (FT == FrequencyType::multigroup) { + opacity = &(jb_pkg->template Param("opacity_d")); + n_nubins = jb_pkg->template Param("n_nubins"); + numin = jb_pkg->template Param("numin"); + numax = jb_pkg->template Param("numax"); + } auto &do_emission = jb_pkg->template Param("do_emission"); auto &source_strategy = jb_pkg->template Param("source_strategy"); PARTHENON_REQUIRE(source_strategy != SourceStrategy::energy, @@ -51,7 +64,8 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { // Create pack static auto desc = MakePackDescriptor(resolved_pkgs.get()); + fj::source_num_per_cell, fj::emission_cdf, fj::energy_delta>( + resolved_pkgs.get()); auto vmesh = desc.GetPack(md); // Indexing and dimensionality @@ -87,12 +101,37 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { [[maybe_unused]] const Real &sbd = sb; [[maybe_unused]] const Real &vvd = vv; [[maybe_unused]] const Real &dtd = dt; - [[maybe_unused]] auto &opac = opacity; + [[maybe_unused]] auto *opac = opacity; + [[maybe_unused]] auto *mopac = mopacity; + [[maybe_unused]] const auto numind = numin; + [[maybe_unused]] const auto numaxd = numax; + [[maybe_unused]] const auto n_nubinsd = n_nubins; Real erad = 0.0; if constexpr (ST == SourceType::thermal) { erad = (4.0 * sbd / vvd) * std::pow(temp, 4.0) * dv; } else if constexpr (ST == SourceType::emission) { - const Real emis = opac.Emissivity(rho, temp); + Real emis = -1.; + if constexpr (FT == FrequencyType::gray) { + emis = opac->Emissivity(rho, temp); + } else if constexpr (FT == FrequencyType::multigroup) { + // Construct emission CDF + const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; + vmesh(b, fj::emission_cdf(0), k, j, i) = opac->EmissivityPerNu( + rho, temp, std::exp(std::log(numind) + 0.5 * dlnu)); + for (int n = 1; n < n_nubinsd; n++) { + const Real nu = std::exp(std::log(numind) + (n + 0.5) * dlnu); + vmesh(b, fj::emission_cdf(n), k, j, i) = + opac->EmissivityPerNu(rho, temp, nu) + + vmesh(b, fj::emission_cdf(n - 1), k, j, i); + } + for (int n = 0; n < n_nubinsd; n++) { + // Get total emissivity + emis += vmesh(b, fj::emission_cdf(n), k, j, i); + // Normalize emission CDF + vmesh(b, fj::emission_cdf(n), k, j, i) /= + vmesh(b, fj::emission_cdf(n_nubins - 1), k, j, i); + } + } erad = vmesh(b, fj::fleck_factor(), k, j, i) * emis * dv * dtd; } // Set source_num_per_cell @@ -156,6 +195,10 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { auto rng_gen = rng_pool.get_state(); [[maybe_unused]] const Real &dtd = dt; [[maybe_unused]] const Real &t_startd = t_start; + [[maybe_unused]] const Real hd = h; + [[maybe_unused]] const Real numaxd = numax; + [[maybe_unused]] const Real numind = numin; + [[maybe_unused]] const int n_nubinsd = n_nubins; // Starting index and length of particles in this cell const int &pstart_idx = prefix_sum(b, cell_idx_1d); @@ -189,8 +232,23 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { const Real &rho = vmesh(b, fjh::density(), k, j, i); const Real &sie = vmesh(b, fjh::sie(), k, j, i); const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie); - ppack_r(b, ph::energy(), n) = sample_Planck_energy(rng_gen, sb, temp); - ppack_r(b, ph::weight(), n) = vmesh(b, fj::source_ew_per_cell(), k, j, i); + if constexpr (FT == FrequencyType::gray) { + ppack_r(b, ph::energy(), n) = sample_Planck_energy(rng_gen, sb, temp); + } else if constexpr (FT == FrequencyType::multigroup) { + // Sample energy from CDF + const Real rand = rng_gen.drand(); + int n; + for (n = 0; n < n_nubinsd; n++) { + if (vmesh(b, fj::emission_cdf(n), k, j, i) >= rand) { + break; + } + } + const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; + const Real nu = std::exp(std::log(numind) + (n + 0.5) * dlnu); + ppack_r(b, ph::energy(), n) = hd * nu; + } + ppack_r(b, ph::weight(), n) = + vmesh(b, fj::source_ew_per_cell(), k, j, i) / ppack_r(b, ph::energy(), n); if constexpr (ST == SourceType::emission) { dejbn -= ppack_r(b, ph::weight(), n); @@ -212,9 +270,22 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { typedef MeshBlockData BD; typedef MeshData D; typedef SourceType ST; -template TaskStatus SourcePhotons(BD *md, const Real t0, const Real dt); -template TaskStatus SourcePhotons(BD *md, const Real t0, const Real dt); -template TaskStatus SourcePhotons(D *md, const Real t0, const Real dt); -template TaskStatus SourcePhotons(D *md, const Real t0, const Real dt); +typedef FrequencyType FT; +template TaskStatus SourcePhotons(BD *md, const Real t0, + const Real dt); +template TaskStatus SourcePhotons(BD *md, const Real t0, + const Real dt); +template TaskStatus SourcePhotons(D *md, const Real t0, + const Real dt); +template TaskStatus SourcePhotons(D *md, const Real t0, + const Real dt); +template TaskStatus SourcePhotons(BD *md, const Real t0, + const Real dt); +template TaskStatus SourcePhotons(BD *md, const Real t0, + const Real dt); +template TaskStatus SourcePhotons(D *md, const Real t0, + const Real dt); +template TaskStatus SourcePhotons(D *md, const Real t0, + const Real dt); } // namespace jaybenne From da7c6ec109c998877025def10b510e1a65145517 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Wed, 18 Dec 2024 13:46:39 -0700 Subject: [PATCH 05/16] n_nubinsd, emis poison fix, work on transporters --- src/jaybenne/jaybenne.cpp | 20 ++++++++++--- src/jaybenne/jaybenne.hpp | 2 ++ src/jaybenne/sourcing.cpp | 3 +- src/jaybenne/transport.cpp | 51 ++++++++++++++++++++++++--------- src/jaybenne/transport_ddmc.cpp | 45 +++++++++++++++++------------ 5 files changed, 85 insertions(+), 36 deletions(-) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 1a39c55..875627c 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -103,7 +103,7 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) { // prepare for iterative transport loop auto derived = tl.AddTask(none, UpdateDerivedTransportFields, base.get(), dt); - TaskID source = derived; + auto source = derived; if (fd == FrequencyType::gray) { auto source = tl.AddTask( derived, @@ -127,9 +127,21 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) { Kokkos::Profiling::pushRegion("Jaybenne::TransportLoop"); return TaskStatus::complete; }); - auto transport = - use_ddmc ? itl.AddTask(time_start, TransportPhotons_DDMC, base.get(), t_start, dt) - : itl.AddTask(time_start, TransportPhotons, base.get(), t_start, dt); + auto transport = time_start; + if (fd == FrequencyType::gray) { + transport = + use_ddmc ? itl.AddTask(time_start, TransportPhotons_DDMC, + base.get(), t_start, dt) + : itl.AddTask(time_start, TransportPhotons, + base.get(), t_start, dt); + } else if (fd == FrequencyType::multigroup) { + transport = + use_ddmc + ? itl.AddTask(time_start, TransportPhotons_DDMC, + base.get(), t_start, dt) + : itl.AddTask(time_start, TransportPhotons, + base.get(), t_start, dt); + } auto reset_comms = itl.AddTask(transport, MeshResetCommunication, base.get()); auto send = itl.AddTask(reset_comms, MeshSend, base.get()); auto receive = itl.AddTask(transport | send, MeshReceive, base.get()); diff --git a/src/jaybenne/jaybenne.hpp b/src/jaybenne/jaybenne.hpp index 60dd5ab..a2b8ff1 100644 --- a/src/jaybenne/jaybenne.hpp +++ b/src/jaybenne/jaybenne.hpp @@ -61,7 +61,9 @@ enum class SourceType { thermal, emission }; enum class FrequencyType { gray, multigroup }; // Tasks +template TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real dt); +template TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const Real dt); TaskStatus SampleDDMCBlockFace(MeshData *md); TaskStatus CheckCompletion(MeshData *md, const Real t_end); diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index a07b597..3406e37 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -124,12 +124,13 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { opac->EmissivityPerNu(rho, temp, nu) + vmesh(b, fj::emission_cdf(n - 1), k, j, i); } + emis = 0.; for (int n = 0; n < n_nubinsd; n++) { // Get total emissivity emis += vmesh(b, fj::emission_cdf(n), k, j, i); // Normalize emission CDF vmesh(b, fj::emission_cdf(n), k, j, i) /= - vmesh(b, fj::emission_cdf(n_nubins - 1), k, j, i); + vmesh(b, fj::emission_cdf(n_nubinsd - 1), k, j, i); } } erad = vmesh(b, fj::fleck_factor(), k, j, i) * emis * dv * dtd; diff --git a/src/jaybenne/transport.cpp b/src/jaybenne/transport.cpp index 1e7c292..7d649b4 100644 --- a/src/jaybenne/transport.cpp +++ b/src/jaybenne/transport.cpp @@ -25,6 +25,7 @@ namespace jaybenne { //---------------------------------------------------------------------------------------- //! \fn TaskStatus TransportPhotons //! \brief +template TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real dt) { namespace fj = field::jaybenne; namespace fjh = field::jaybenne::host; @@ -34,11 +35,26 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; auto &jb_pkg = pm->packages.Get("jaybenne"); - auto &eos = jb_pkg->Param("eos_d"); - auto &opacity = jb_pkg->Param("opacity_d"); - auto &scattering = jb_pkg->Param("scattering_d"); - auto &rng_pool = jb_pkg->Param("rng_pool"); - const Real vv = jb_pkg->Param("speed_of_light"); + auto &eos = jb_pkg->template Param("eos_d"); + const Opacity *opacity = nullptr; + const Scattering *scattering = nullptr; + const MeanOpacity *mopacity = nullptr; + const MeanScattering *mscattering = nullptr; + int n_nubins = -1; + Real numin = -1.; + Real numax = -1.; + if constexpr (FT == FrequencyType::gray) { + mopacity = &(jb_pkg->template Param("mopacity_d")); + mscattering = &(jb_pkg->template Param("mscattering_d")); + } else if constexpr (FT == FrequencyType::multigroup) { + opacity = &(jb_pkg->template Param("opacity_d")); + scattering = &(jb_pkg->template Param("scattering_d")); + n_nubins = jb_pkg->template Param("n_nubins"); + numin = jb_pkg->template Param("numin"); + numax = jb_pkg->template Param("numax"); + } + auto &rng_pool = jb_pkg->template Param("rng_pool"); + const Real vv = jb_pkg->template Param("speed_of_light"); // Create SparsePack static auto desc = @@ -111,20 +127,20 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d "Particle initially outside X3 logical bnds!"); // Calculate cell bounds - const Real xl = coords.Xc(ip) - 0.5 * dx_i; - const Real xu = coords.Xc(ip) + 0.5 * dx_i; - const Real yl = coords.Xc(jp) - 0.5 * dx_j; - const Real yu = coords.Xc(jp) + 0.5 * dx_j; - const Real zl = coords.Xc(kp) - 0.5 * dx_k; - const Real zu = coords.Xc(kp) + 0.5 * dx_k; + const Real xl = coords.template Xc(ip) - 0.5 * dx_i; + const Real xu = coords.template Xc(ip) + 0.5 * dx_i; + const Real yl = coords.template Xc(jp) - 0.5 * dx_j; + const Real yu = coords.template Xc(jp) + 0.5 * dx_j; + const Real zl = coords.template Xc(kp) - 0.5 * dx_k; + const Real zu = coords.template Xc(kp) + 0.5 * dx_k; // Extract physical quantities const Real &rho = vmesh(b, fjh::density(), kp, jp, ip); const Real &sie = vmesh(b, fjh::sie(), kp, jp, ip); const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie); const Real &ff = vmesh(b, fj::fleck_factor(), kp, jp, ip); - const Real ss = scattering.TotalScatteringCoefficient(rho, temp, ee); - const Real aa = opacity.AbsorptionCoefficient(rho, temp, ee); + const Real ss = scattering->TotalScatteringCoefficient(rho, temp, ee); + const Real aa = opacity->AbsorptionCoefficient(rho, temp, ee); // reset collision indicators bool is_absorbed = false; @@ -215,4 +231,13 @@ TaskStatus CheckCompletion(MeshData *md, const Real t_end) { } } +//---------------------------------------------------------------------------------------- +//! template instantiations +template TaskStatus TransportPhotons(MeshData *md, + const Real t_start, + const Real dt); +template TaskStatus TransportPhotons(MeshData *md, + const Real t_start, + const Real dt); + } // namespace jaybenne diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index 6c626b5..6959e92 100644 --- a/src/jaybenne/transport_ddmc.cpp +++ b/src/jaybenne/transport_ddmc.cpp @@ -25,6 +25,7 @@ namespace jaybenne { //---------------------------------------------------------------------------------------- //! \fn TaskStatus TransportPhotons_DDMC //! \brief +template TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const Real dt) { namespace fj = field::jaybenne; namespace fjh = field::jaybenne::host; @@ -35,12 +36,12 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; auto &jb_pkg = pm->packages.Get("jaybenne"); - auto &eos = jb_pkg->Param("eos_d"); - auto &opacity = jb_pkg->Param("opacity_d"); - auto &scattering = jb_pkg->Param("scattering_d"); - auto &rng_pool = jb_pkg->Param("rng_pool"); - const Real vv = jb_pkg->Param("speed_of_light"); - const Real &tau_ddmc = jb_pkg->Param("tau_ddmc"); + auto &eos = jb_pkg->template Param("eos_d"); + auto &opacity = jb_pkg->template Param("opacity_d"); + auto &scattering = jb_pkg->template Param("scattering_d"); + auto &rng_pool = jb_pkg->template Param("rng_pool"); + const Real vv = jb_pkg->template Param("speed_of_light"); + const Real &tau_ddmc = jb_pkg->template Param("tau_ddmc"); // Create SparsePack static auto desc = @@ -113,12 +114,12 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R "Particle initially outside X3 logical bnds!"); // calculate cell bounds - const Real xl = coords.Xc(ip) - 0.5 * dx_i; - const Real xu = coords.Xc(ip) + 0.5 * dx_i; - const Real yl = coords.Xc(jp) - 0.5 * dx_j; - const Real yu = coords.Xc(jp) + 0.5 * dx_j; - const Real zl = coords.Xc(kp) - 0.5 * dx_k; - const Real zu = coords.Xc(kp) + 0.5 * dx_k; + const Real xl = coords.template Xc(ip) - 0.5 * dx_i; + const Real xu = coords.template Xc(ip) + 0.5 * dx_i; + const Real yl = coords.template Xc(jp) - 0.5 * dx_j; + const Real yu = coords.template Xc(jp) + 0.5 * dx_j; + const Real zl = coords.template Xc(kp) - 0.5 * dx_k; + const Real zu = coords.template Xc(kp) + 0.5 * dx_k; // Extract physical quantities const Real &rho = vmesh(b, fjh::density(), kp, jp, ip); @@ -139,12 +140,12 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R swarm_d.Xtoijk(x, y, z, ip, jp, kp); // calculate cell bounds - const Real xl = coords.Xc(ip) - 0.5 * dx_i; - const Real xu = coords.Xc(ip) + 0.5 * dx_i; - const Real yl = coords.Xc(jp) - 0.5 * dx_j; - const Real yu = coords.Xc(jp) + 0.5 * dx_j; - const Real zl = coords.Xc(kp) - 0.5 * dx_k; - const Real zu = coords.Xc(kp) + 0.5 * dx_k; + const Real xl = coords.template Xc(ip) - 0.5 * dx_i; + const Real xu = coords.template Xc(ip) + 0.5 * dx_i; + const Real yl = coords.template Xc(jp) - 0.5 * dx_j; + const Real yu = coords.template Xc(jp) + 0.5 * dx_j; + const Real zl = coords.template Xc(kp) - 0.5 * dx_k; + const Real zu = coords.template Xc(kp) + 0.5 * dx_k; // get face probabilities const Real &Px_l = vmesh(b, TE::F1, fj::ddmc_face_prob(), kp, jp, ip); @@ -235,5 +236,13 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R return TaskStatus::complete; } +//---------------------------------------------------------------------------------------- +//! template instantiations +template TaskStatus TransportPhotons_DDMC(MeshData *md, + const Real t_start, + const Real dt); +template TaskStatus TransportPhotons_DDMC(MeshData *md, + const Real t_start, + const Real dt); } // namespace jaybenne From a4b0491f430e5713a4ea438b74e3adbe6fe8261e Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Thu, 19 Dec 2024 15:41:03 -0700 Subject: [PATCH 06/16] Fixing little errors --- inputs/inf.in | 1 + src/mcblock/mcblock.cpp | 142 +++++++++++++++++++++++++++++----------- src/mcblock/mcblock.hpp | 3 +- src/mcblock/opacity.hpp | 6 +- 4 files changed, 109 insertions(+), 43 deletions(-) diff --git a/inputs/inf.in b/inputs/inf.in index 66a780d..367ceb8 100644 --- a/inputs/inf.in +++ b/inputs/inf.in @@ -58,6 +58,7 @@ source_strategy = uniform seed = 349857 +frequency_type = gray opacity_model = constant opacity_constant_value = 1.0 # cm^2/g scattering_model = constant diff --git a/src/mcblock/mcblock.cpp b/src/mcblock/mcblock.cpp index d5496bd..689e6e1 100644 --- a/src/mcblock/mcblock.cpp +++ b/src/mcblock/mcblock.cpp @@ -91,58 +91,115 @@ std::shared_ptr Initialize(ParameterInput *pin) { pkg->AddParam<>("length_scale", length_scale); pkg->AddParam<>("temperature_scale", temperature_scale); + // Frequency discretization + std::string frequency_type_name = pin->GetString("mcblock", "frequency_type"); + FrequencyType 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!"); + } + // Absorption opacity model - OpacityModel opacity_model; Opacity opacity; - std::string opacity_model_name = pin->GetString("mcblock", "opacity_model"); - if (opacity_model_name == "none") { - opacity_model = OpacityModel::none; - opacity = singularity::photons::NonCGSUnits( - singularity::photons::Gray(0.0), time_scale, mass_scale, length_scale, - temperature_scale); - } else if (opacity_model_name == "constant") { - opacity_model = OpacityModel::constant; - Real kappa = pin->GetReal("mcblock", "opacity_constant_value"); - opacity = singularity::photons::NonCGSUnits( - singularity::photons::Gray(kappa), time_scale, mass_scale, length_scale, - temperature_scale); - } else if (opacity_model_name == "ep_bremss") { - opacity_model = OpacityModel::epbremss; - opacity = singularity::photons::NonCGSUnits( - singularity::photons::EPBremss(), time_scale, mass_scale, length_scale, - temperature_scale); - } else { - // nothing else supported for now - PARTHENON_FAIL("Only none or constant opacity models supported!"); + MeanOpacity mopacity; + std::string opacity_model = pin->GetString("mcblock", "opacity_model"); + if (frequency_type == FrequencyType::gray) { + if (opacity_model == "none") { + auto opac = singularity::photons::Gray(1.e-100); + mopacity = + singularity::photons::MeanNonCGSUnits( + singularity::photons::MeanOpacityBase(opac, 1.e-10, 1.e10, 1, 1.e-10, 1.e10, + 1), + time_scale, mass_scale, length_scale, 1.); + } else if (opacity_model == "constant") { + Real kappa = pin->GetReal("mcblock", "opacity_constant_value"); + auto opac = singularity::photons::Gray(kappa); + mopacity = + singularity::photons::MeanNonCGSUnits( + singularity::photons::MeanOpacityBase(opac, 1.e-10, 1.e10, 1, 1.e-10, 1.e10, + 1), + time_scale, mass_scale, length_scale, 1.); + } else if (opacity_model == "table") { + std::string table_filename = pin->GetString("mcblock", "opacity_table"); + mopacity = + singularity::photons::MeanNonCGSUnits( + singularity::photons::MeanOpacityBase(table_filename), time_scale, + mass_scale, length_scale, 1.); + } else { + PARTHENON_FAIL("Only none, constant, or table opacity models supported!"); + } + } else if (frequency_type == FrequencyType::multigroup) { + if (opacity_model == "none") { + opacity = singularity::photons::NonCGSUnits( + singularity::photons::Gray(0.0), time_scale, mass_scale, length_scale, + temperature_scale); + } else if (opacity_model == "constant") { + Real kappa = pin->GetReal("mcblock", "opacity_constant_value"); + opacity = singularity::photons::NonCGSUnits( + singularity::photons::Gray(kappa), time_scale, mass_scale, length_scale, + temperature_scale); + } else if (opacity_model == "ep_bremss") { + opacity = singularity::photons::NonCGSUnits( + singularity::photons::EPBremss(), time_scale, mass_scale, length_scale, + temperature_scale); + } else { + // nothing else supported for now + PARTHENON_FAIL("Only none or constant opacity models supported!"); + } } + pkg->AddParam<>("frequency_type", frequency_type); pkg->AddParam<>("opacity_model", opacity_model); pkg->AddParam<>("opacity_h", opacity); + pkg->AddParam<>("mopacity_h", mopacity); // Scattering opacity model // TODO(BRR) Remove apm with switch in singularity-opac to cm^2/g opacities? const Real apm = pin->GetOrAddReal("mcblock", "apm", 1.); // Average particle mass (code units) pkg->AddParam<>("apm", apm); - ScatteringModel scattering_model; Scattering scattering; - std::string scattering_model_name = + MeanScattering mscattering; + std::string scattering_model = pin->GetOrAddString("mcblock", "scattering_model", "none"); - if (scattering_model_name == "none") { - scattering_model = ScatteringModel::none; - scattering = singularity::photons::NonCGSUnitsS( - singularity::photons::GrayS(0.0, apm), time_scale, mass_scale, length_scale, - temperature_scale); - } else if (scattering_model_name == "constant") { - scattering_model = ScatteringModel::constant; - Real kappa_s = pin->GetReal("mcblock", "scattering_constant_value"); - scattering = singularity::photons::NonCGSUnitsS( - singularity::photons::GrayS(kappa_s, apm), time_scale, mass_scale, length_scale, - temperature_scale); - } else { - PARTHENON_FAIL("Only none or constant scattering models supported!"); + if (frequency_type == FrequencyType::gray) { + if (scattering_model == "none") { + auto sopac = singularity::photons::GrayS(0.0, apm); + mscattering = + singularity::photons::MeanNonCGSUnitsS( + singularity::photons::MeanSOpacityCGS(sopac, 1.e-10, 1.e10, 1, 1.e-10, + 1.e10, 1), + time_scale, mass_scale, length_scale, 1.); + } else if (scattering_model == "constant") { + Real kappa_s = pin->GetReal("mcblock", "scattering_constant_value"); + auto sopac = singularity::photons::GrayS(kappa_s, apm); + mscattering = + singularity::photons::MeanNonCGSUnitsS( + singularity::photons::MeanSOpacityCGS(sopac, 1.e-10, 1.e10, 1, 1.e-10, + 1.e10, 1), + time_scale, mass_scale, length_scale, 1.); + } else { + PARTHENON_FAIL("Only none or constant scattering models supported!"); + } + } else if (frequency_type == FrequencyType::multigroup) { + if (scattering_model == "none") { + scattering = singularity::photons::NonCGSUnitsS( + singularity::photons::GrayS(0.0, apm), time_scale, mass_scale, length_scale, + temperature_scale); + } else if (scattering_model == "constant") { + Real kappa_s = pin->GetReal("mcblock", "scattering_constant_value"); + scattering = singularity::photons::NonCGSUnitsS( + singularity::photons::GrayS(kappa_s, apm), time_scale, mass_scale, length_scale, + temperature_scale); + } else { + PARTHENON_FAIL("Only none or constant scattering models supported!"); + } } pkg->AddParam<>("scattering_model", scattering_model); pkg->AddParam<>("scattering_h", scattering); + pkg->AddParam<>("mscattering_h", mscattering); pkg->PreFillDerivedMesh = UpdateDerived; @@ -290,9 +347,16 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { packages.Add(mcblock::Initialize(pin.get())); auto &mcblock = packages.Get("mcblock"); auto eos_h = mcblock->Param("eos_h"); - auto opacity_h = mcblock->Param("opacity_h"); - auto scattering_h = mcblock->Param("scattering_h"); - packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h)); + auto frequency_type = mcblock->Param("frequency_type"); + if (frequency_type == FrequencyType::gray) { + auto mopacity_h = mcblock->Param("mopacity_h"); + auto mscattering_h = mcblock->Param("mscattering_h"); + packages.Add(jaybenne::Initialize(pin.get(), mopacity_h, mscattering_h, eos_h)); + } else if (frequency_type == FrequencyType::multigroup) { + auto opacity_h = mcblock->Param("opacity_h"); + auto scattering_h = mcblock->Param("scattering_h"); + packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h)); + } return packages; } diff --git a/src/mcblock/mcblock.hpp b/src/mcblock/mcblock.hpp index a186a92..faeb23d 100644 --- a/src/mcblock/mcblock.hpp +++ b/src/mcblock/mcblock.hpp @@ -29,8 +29,7 @@ namespace mcblock { // Model enums enum class InitialRadiation { none, thermal }; -enum class OpacityModel { none, constant, epbremss }; -enum class ScatteringModel { none, constant }; +enum class FrequencyType { gray, multigroup }; std::shared_ptr Initialize(ParameterInput *pin); void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin); diff --git a/src/mcblock/opacity.hpp b/src/mcblock/opacity.hpp index e32767a..4396d09 100644 --- a/src/mcblock/opacity.hpp +++ b/src/mcblock/opacity.hpp @@ -26,14 +26,16 @@ using Opacity = singularity::photons::impl::Variant< singularity::photons::NonCGSUnits, singularity::photons::NonCGSUnits>; -using MeanOpacity = singularity::photons::MeanOpacity; +using MeanOpacity = + singularity::photons::MeanNonCGSUnits; // Reduced scattering variant just for jaybenne using Scattering = singularity::photons::impl::S_Variant< singularity::photons::NonCGSUnitsS, singularity::photons::NonCGSUnitsS>; -using MeanScattering = singularity::photons::MeanSOpacity; +using MeanScattering = + singularity::photons::MeanNonCGSUnitsS; } // namespace mcblock From 3ded5a3a091be65a64f5242d34b4fcf211e1a76e Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 20 Dec 2024 14:37:07 -0700 Subject: [PATCH 07/16] Need to fix UpdateDerivedTransportFields --- src/jaybenne/jaybenne.cpp | 2 +- src/mcblock/mcblock.cpp | 16 ++++++---------- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 875627c..6b43dad 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -290,7 +290,7 @@ std::shared_ptr Initialize(ParameterInput *pin, MeanOpacity &mo auto pkg = Initialize_impl(pin, eos, units); - pkg->AddParam<>("frequency_type", FrequencyType::multigroup); + pkg->AddParam<>("frequency_type", FrequencyType::gray); // Opacity model pkg->AddParam<>("mopacity_d", mopacity.GetOnDevice()); diff --git a/src/mcblock/mcblock.cpp b/src/mcblock/mcblock.cpp index 689e6e1..b9ae449 100644 --- a/src/mcblock/mcblock.cpp +++ b/src/mcblock/mcblock.cpp @@ -111,17 +111,15 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto opac = singularity::photons::Gray(1.e-100); mopacity = singularity::photons::MeanNonCGSUnits( - singularity::photons::MeanOpacityBase(opac, 1.e-10, 1.e10, 1, 1.e-10, 1.e10, - 1), - time_scale, mass_scale, length_scale, 1.); + singularity::photons::MeanOpacityBase(opac, -1, 1, 1, -1, 1, 1), time_scale, + mass_scale, length_scale, 1.); } else if (opacity_model == "constant") { Real kappa = pin->GetReal("mcblock", "opacity_constant_value"); auto opac = singularity::photons::Gray(kappa); mopacity = singularity::photons::MeanNonCGSUnits( - singularity::photons::MeanOpacityBase(opac, 1.e-10, 1.e10, 1, 1.e-10, 1.e10, - 1), - time_scale, mass_scale, length_scale, 1.); + singularity::photons::MeanOpacityBase(opac, -1, 1, 2, -1, 1, 2), time_scale, + mass_scale, length_scale, 1.); } else if (opacity_model == "table") { std::string table_filename = pin->GetString("mcblock", "opacity_table"); mopacity = @@ -169,16 +167,14 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto sopac = singularity::photons::GrayS(0.0, apm); mscattering = singularity::photons::MeanNonCGSUnitsS( - singularity::photons::MeanSOpacityCGS(sopac, 1.e-10, 1.e10, 1, 1.e-10, - 1.e10, 1), + singularity::photons::MeanSOpacityCGS(sopac, -1., 1., 2, -1., 1., 2), time_scale, mass_scale, length_scale, 1.); } else if (scattering_model == "constant") { Real kappa_s = pin->GetReal("mcblock", "scattering_constant_value"); auto sopac = singularity::photons::GrayS(kappa_s, apm); mscattering = singularity::photons::MeanNonCGSUnitsS( - singularity::photons::MeanSOpacityCGS(sopac, 1.e-10, 1.e10, 1, 1.e-10, - 1.e10, 1), + singularity::photons::MeanSOpacityCGS(sopac, -1., 1., 2, -1., 1., 2), time_scale, mass_scale, length_scale, 1.); } else { PARTHENON_FAIL("Only none or constant scattering models supported!"); From f9f2a711efff9905664f311875928b2b6677c37d Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 20 Dec 2024 15:32:26 -0700 Subject: [PATCH 08/16] multiply Jnu by dnu --- src/jaybenne/jaybenne.cpp | 112 ++++++++++++++++++++++++++++++-------- src/jaybenne/sourcing.cpp | 15 ++--- 2 files changed, 98 insertions(+), 29 deletions(-) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 6b43dad..b5ed2d6 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -346,14 +346,15 @@ Real EstimateTimestepMesh(MeshData *md) { } //---------------------------------------------------------------------------------------- -//! \fn TaskStatus UpdateDerivedTransportFields +//! \fn TaskStatus UpdateDerivedTransportFieldsImpl //! \brief fields set: Fleck factor, DDMC face probabilities //! Fleck factor formula: //! betaf = 4aT^3 / (rho * cv) //! f = 1 / (1 + betaf * opacP * c * dt) //! NOTE: if J = opacP * c * aR * T^4, //! then f = 1 / (1 + 4 * J * dt / (rho * cv * T)) -TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt) { +template +TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { namespace fj = field::jaybenne; namespace fjh = field::jaybenne::host; @@ -361,7 +362,17 @@ TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt) { auto &resolved_pkgs = pm->resolved_packages; auto &jbn = pm->packages.Get("jaybenne"); auto &eos = jbn->Param("eos_d"); - auto &opacity = jbn->Param("opacity_d"); + Opacity opacity; + Scattering scattering; + MeanOpacity mopacity; + MeanScattering mscattering; + if constexpr (FT == FrequencyType::gray) { + mopacity = jbn->Param("mopacity_d"); + mscattering = jbn->Param("mscattering_d"); + } else if constexpr (FT == FrequencyType::multigroup) { + opacity = jbn->Param("opacity_d"); + scattering = jbn->Param("scattering_d"); + } const auto &ib = md->GetBoundsI(IndexDomain::interior); const auto &jb = md->GetBoundsJ(IndexDomain::interior); @@ -381,7 +392,13 @@ TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt) { const Real &sie = vmesh(b, fjh::sie(), k, j, i); const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie); const Real cv = eos.SpecificHeatFromDensityInternalEnergy(rho, sie); - const Real emis = opacity.Emissivity(rho, temp); + Real emis; + if constexpr (FT == FrequencyType::gray) { + emis = mopacity.Emissivity(rho, temp); + + } else if constexpr (FT == FrequencyType::multigroup) { + emis = opacity.Emissivity(rho, temp); + } vmesh(b, fj::fleck_factor(), k, j, i) = 1.0 / (1.0 + (4.0 * emis / (rho * cv * temp)) * dt); }); @@ -399,9 +416,6 @@ TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt) { // get DDMC cell optical thickness threshold const Real tau_ddmc = jbn->Param("tau_ddmc"); - // get scattering to calculate DDMC threshold - auto &scattering = jbn->Param("scattering_d"); - // calculate DDMC face probabilities in X1 direction const int iu = ib.e + 1; parthenon::par_for( @@ -432,11 +446,22 @@ TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt) { const Real &rho_u = vmesh(b, fjh::density(), k, j, i); const Real &sie_u = vmesh(b, fjh::sie(), k, j, i); const Real temp_u = eos.TemperatureFromDensityInternalEnergy(rho_u, sie_u); - // TODO: replace 3rd argument when this routine operates in multigroup - const Real ss_l = scattering.TotalScatteringCoefficient(rho_l, temp_l, 1.0); - const Real aa_l = opacity.AbsorptionCoefficient(rho_l, temp_l, 1.0); - const Real ss_u = scattering.TotalScatteringCoefficient(rho_u, temp_u, 1.0); - const Real aa_u = opacity.AbsorptionCoefficient(rho_u, temp_u, 1.0); + Real ss_l; + Real aa_l; + Real ss_u; + Real aa_u; + if constexpr (FT == FrequencyType::gray) { + ss_l = mscattering.RosselandMeanTotalScatteringCoefficient(rho_l, temp_l); + aa_l = mopacity.RosselandMeanAbsorptionCoefficient(rho_l, temp_l); + ss_u = mscattering.RosselandMeanTotalScatteringCoefficient(rho_u, temp_u); + aa_u = mopacity.RosselandMeanAbsorptionCoefficient(rho_u, temp_u); + } else if constexpr (FT == FrequencyType::multigroup) { + // TODO: replace 3rd argument when this routine operates in multigroup + ss_l = scattering.TotalScatteringCoefficient(rho_l, temp_l, 1.0); + aa_l = opacity.AbsorptionCoefficient(rho_l, temp_l, 1.0); + ss_u = scattering.TotalScatteringCoefficient(rho_u, temp_u, 1.0); + aa_u = opacity.AbsorptionCoefficient(rho_u, temp_u, 1.0); + } // calculate optical thicknesses from lower and upper cell Real tau_l = dx_lx * (ss_l + aa_l); @@ -479,11 +504,22 @@ TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt) { const Real &rho_u = vmesh(b, fjh::density(), k, j, i); const Real &sie_u = vmesh(b, fjh::sie(), k, j, i); const Real temp_u = eos.TemperatureFromDensityInternalEnergy(rho_u, sie_u); - // TODO: replace 3rd argument when this routine operates in multigroup - const Real ss_l = scattering.TotalScatteringCoefficient(rho_l, temp_l, 1.0); - const Real aa_l = opacity.AbsorptionCoefficient(rho_l, temp_l, 1.0); - const Real ss_u = scattering.TotalScatteringCoefficient(rho_u, temp_u, 1.0); - const Real aa_u = opacity.AbsorptionCoefficient(rho_u, temp_u, 1.0); + Real ss_l; + Real aa_l; + Real ss_u; + Real aa_u; + if constexpr (FT == FrequencyType::gray) { + ss_l = mscattering.RosselandMeanTotalScatteringCoefficient(rho_l, temp_l); + aa_l = mopacity.RosselandMeanAbsorptionCoefficient(rho_l, temp_l); + ss_u = mscattering.RosselandMeanTotalScatteringCoefficient(rho_u, temp_u); + aa_u = mopacity.RosselandMeanAbsorptionCoefficient(rho_u, temp_u); + } else if constexpr (FT == FrequencyType::multigroup) { + // TODO: replace 3rd argument when this routine operates in multigroup + ss_l = scattering.TotalScatteringCoefficient(rho_l, temp_l, 1.0); + aa_l = opacity.AbsorptionCoefficient(rho_l, temp_l, 1.0); + ss_u = scattering.TotalScatteringCoefficient(rho_u, temp_u, 1.0); + aa_u = opacity.AbsorptionCoefficient(rho_u, temp_u, 1.0); + } // calculate optical thicknesses from lower and upper cell Real tau_l = dx_ly * (ss_l + aa_l); @@ -528,11 +564,22 @@ TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt) { const Real &rho_u = vmesh(b, fjh::density(), k, j, i); const Real &sie_u = vmesh(b, fjh::sie(), k, j, i); const Real temp_u = eos.TemperatureFromDensityInternalEnergy(rho_u, sie_u); - // TODO: replace 3rd argument when this routine operates in multigroup - const Real ss_l = scattering.TotalScatteringCoefficient(rho_l, temp_l, 1.0); - const Real aa_l = opacity.AbsorptionCoefficient(rho_l, temp_l, 1.0); - const Real ss_u = scattering.TotalScatteringCoefficient(rho_u, temp_u, 1.0); - const Real aa_u = opacity.AbsorptionCoefficient(rho_u, temp_u, 1.0); + Real ss_l; + Real aa_l; + Real ss_u; + Real aa_u; + if constexpr (FT == FrequencyType::gray) { + ss_l = mscattering.RosselandMeanTotalScatteringCoefficient(rho_l, temp_l); + aa_l = mopacity.RosselandMeanAbsorptionCoefficient(rho_l, temp_l); + ss_u = mscattering.RosselandMeanTotalScatteringCoefficient(rho_u, temp_u); + aa_u = mopacity.RosselandMeanAbsorptionCoefficient(rho_u, temp_u); + } else if constexpr (FT == FrequencyType::multigroup) { + // TODO: replace 3rd argument when this routine operates in multigroup + ss_l = scattering.TotalScatteringCoefficient(rho_l, temp_l, 1.0); + aa_l = opacity.AbsorptionCoefficient(rho_l, temp_l, 1.0); + ss_u = scattering.TotalScatteringCoefficient(rho_u, temp_u, 1.0); + aa_u = opacity.AbsorptionCoefficient(rho_u, temp_u, 1.0); + } // calculate optical thicknesses from lower and upper cell Real tau_l = dx_lz * (ss_l + aa_l); @@ -550,6 +597,27 @@ TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt) { return TaskStatus::complete; } +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus UpdateDerivedTransportFieldsImpl +//! \brief fields set: Fleck factor, DDMC face probabilities +//! Fleck factor formula: +//! betaf = 4aT^3 / (rho * cv) +//! f = 1 / (1 + betaf * opacP * c * dt) +//! NOTE: if J = opacP * c * aR * T^4, +//! then f = 1 / (1 + 4 * J * dt / (rho * cv * T)) +TaskStatus UpdateDerivedTransportFields(MeshData *md, const Real dt) { + auto pm = md->GetParentPointer(); + auto &jbn = pm->packages.Get("jaybenne"); + auto &frequency_type = jbn->Param("frequency_type"); + if (frequency_type == FrequencyType::gray) { + return UpdateDerivedTransportFieldsImpl(md, dt); + } else if (frequency_type == FrequencyType::multigroup) { + return UpdateDerivedTransportFieldsImpl(md, dt); + } else { + PARTHENON_FAIL("frequency_type not recognized!"); + } +} + //---------------------------------------------------------------------------------------- //! \fn TaskStatus TaskStatus DefragParticles(MeshBlock *pmb) //! \brief NOTE(PDM): currently unused??? diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 3406e37..82ab721 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -33,15 +33,15 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { auto &jb_pkg = pm->packages.Get("jaybenne"); const Real h = jb_pkg->template Param("planck_constant"); auto &eos = jb_pkg->template Param("eos_d"); - const MeanOpacity *mopacity = nullptr; - const Opacity *opacity = nullptr; int n_nubins = -1; Real numin = -1.; Real numax = -1.; + MeanOpacity mopacity; + Opacity opacity; if constexpr (FT == FrequencyType::gray) { - mopacity = &(jb_pkg->template Param("mopacity_d")); + mopacity = jb_pkg->template Param("mopacity_d"); } else if constexpr (FT == FrequencyType::multigroup) { - opacity = &(jb_pkg->template Param("opacity_d")); + opacity = jb_pkg->template Param("opacity_d"); n_nubins = jb_pkg->template Param("n_nubins"); numin = jb_pkg->template Param("numin"); numax = jb_pkg->template Param("numax"); @@ -112,7 +112,7 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { } else if constexpr (ST == SourceType::emission) { Real emis = -1.; if constexpr (FT == FrequencyType::gray) { - emis = opac->Emissivity(rho, temp); + emis = mopac.Emissivity(rho, temp); } else if constexpr (FT == FrequencyType::multigroup) { // Construct emission CDF const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; @@ -121,13 +121,14 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { for (int n = 1; n < n_nubinsd; n++) { const Real nu = std::exp(std::log(numind) + (n + 0.5) * dlnu); vmesh(b, fj::emission_cdf(n), k, j, i) = - opac->EmissivityPerNu(rho, temp, nu) + + opac.EmissivityPerNu(rho, temp, nu) + vmesh(b, fj::emission_cdf(n - 1), k, j, i); } emis = 0.; for (int n = 0; n < n_nubinsd; n++) { + const Real dnu = dlnu * std::exp(std::log(numind) + (n + 0.5) * dlnu); // Get total emissivity - emis += vmesh(b, fj::emission_cdf(n), k, j, i); + emis += vmesh(b, fj::emission_cdf(n), k, j, i) * dnu; // Normalize emission CDF vmesh(b, fj::emission_cdf(n), k, j, i) /= vmesh(b, fj::emission_cdf(n_nubinsd - 1), k, j, i); From 0fc595b75ba3e452b81a497c01e4857fde1088bd Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 20 Dec 2024 15:36:14 -0700 Subject: [PATCH 09/16] fix cdf integration --- src/jaybenne/sourcing.cpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 82ab721..7d06682 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -101,8 +101,8 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { [[maybe_unused]] const Real &sbd = sb; [[maybe_unused]] const Real &vvd = vv; [[maybe_unused]] const Real &dtd = dt; - [[maybe_unused]] auto *opac = opacity; - [[maybe_unused]] auto *mopac = mopacity; + [[maybe_unused]] auto opac = opacity; + [[maybe_unused]] auto mopac = mopacity; [[maybe_unused]] const auto numind = numin; [[maybe_unused]] const auto numaxd = numax; [[maybe_unused]] const auto n_nubinsd = n_nubins; @@ -116,19 +116,23 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { } else if constexpr (FT == FrequencyType::multigroup) { // Construct emission CDF const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; - vmesh(b, fj::emission_cdf(0), k, j, i) = opac->EmissivityPerNu( - rho, temp, std::exp(std::log(numind) + 0.5 * dlnu)); + Real nu = std::exp(std::log(numind) + 0.5 * dlnu) Real dnu = nu * dlnu; + vmesh(b, fj::emission_cdf(0), k, j, i) = + opac.EmissivityPerNu(rho, temp, + std::exp(std::log(numind) + 0.5 * dlnu)) * + dnu; for (int n = 1; n < n_nubinsd; n++) { - const Real nu = std::exp(std::log(numind) + (n + 0.5) * dlnu); + nu = std::exp(std::log(numind) + (n + 0.5) * dlnu); + dnu = dlnu * nu; vmesh(b, fj::emission_cdf(n), k, j, i) = - opac.EmissivityPerNu(rho, temp, nu) + + opac.EmissivityPerNu(rho, temp, nu) * dnu + vmesh(b, fj::emission_cdf(n - 1), k, j, i); } emis = 0.; for (int n = 0; n < n_nubinsd; n++) { const Real dnu = dlnu * std::exp(std::log(numind) + (n + 0.5) * dlnu); // Get total emissivity - emis += vmesh(b, fj::emission_cdf(n), k, j, i) * dnu; + emis += vmesh(b, fj::emission_cdf(n), k, j, i); // Normalize emission CDF vmesh(b, fj::emission_cdf(n), k, j, i) /= vmesh(b, fj::emission_cdf(n_nubinsd - 1), k, j, i); From 4d38f84469db72123342f7526553bec2472e8c49 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 20 Dec 2024 15:46:12 -0700 Subject: [PATCH 10/16] inf example cycles --- src/jaybenne/jaybenne.cpp | 3 +++ src/jaybenne/sourcing.cpp | 3 ++- src/jaybenne/transport.cpp | 27 +++++++++++++++++---------- 3 files changed, 22 insertions(+), 11 deletions(-) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index b5ed2d6..4ce6a9c 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -312,6 +312,9 @@ std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacit auto pkg = Initialize_impl(pin, eos, units); + PARTHENON_REQUIRE(pkg->Param("use_ddmc") == false, + "DDMC not supported for multigroup currently!"); + // Frequency discretization auto time = units.time; Real numin = pin->GetReal("jaybenne", "numin"); // in Hz diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 7d06682..5a27cb8 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -116,7 +116,8 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { } else if constexpr (FT == FrequencyType::multigroup) { // Construct emission CDF const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; - Real nu = std::exp(std::log(numind) + 0.5 * dlnu) Real dnu = nu * dlnu; + Real nu = std::exp(std::log(numind) + 0.5 * dlnu); + Real dnu = nu * dlnu; vmesh(b, fj::emission_cdf(0), k, j, i) = opac.EmissivityPerNu(rho, temp, std::exp(std::log(numind) + 0.5 * dlnu)) * diff --git a/src/jaybenne/transport.cpp b/src/jaybenne/transport.cpp index 7d649b4..8b59777 100644 --- a/src/jaybenne/transport.cpp +++ b/src/jaybenne/transport.cpp @@ -36,19 +36,19 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d auto &resolved_pkgs = pm->resolved_packages; auto &jb_pkg = pm->packages.Get("jaybenne"); auto &eos = jb_pkg->template Param("eos_d"); - const Opacity *opacity = nullptr; - const Scattering *scattering = nullptr; - const MeanOpacity *mopacity = nullptr; - const MeanScattering *mscattering = nullptr; + Opacity opacity; + Scattering scattering; + MeanOpacity mopacity; + MeanScattering mscattering; int n_nubins = -1; Real numin = -1.; Real numax = -1.; if constexpr (FT == FrequencyType::gray) { - mopacity = &(jb_pkg->template Param("mopacity_d")); - mscattering = &(jb_pkg->template Param("mscattering_d")); + mopacity = jb_pkg->template Param("mopacity_d"); + mscattering = jb_pkg->template Param("mscattering_d"); } else if constexpr (FT == FrequencyType::multigroup) { - opacity = &(jb_pkg->template Param("opacity_d")); - scattering = &(jb_pkg->template Param("scattering_d")); + opacity = jb_pkg->template Param("opacity_d"); + scattering = jb_pkg->template Param("scattering_d"); n_nubins = jb_pkg->template Param("n_nubins"); numin = jb_pkg->template Param("numin"); numax = jb_pkg->template Param("numax"); @@ -139,8 +139,15 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d const Real &sie = vmesh(b, fjh::sie(), kp, jp, ip); const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie); const Real &ff = vmesh(b, fj::fleck_factor(), kp, jp, ip); - const Real ss = scattering->TotalScatteringCoefficient(rho, temp, ee); - const Real aa = opacity->AbsorptionCoefficient(rho, temp, ee); + Real ss; + Real aa; + if constexpr (FT == FrequencyType::gray) { + ss = mscattering.RosselandMeanTotalScatteringCoefficient(rho, temp); + aa = mopacity.RosselandMeanAbsorptionCoefficient(rho, temp); + } else if constexpr (FT == FrequencyType::multigroup) { + ss = scattering.TotalScatteringCoefficient(rho, temp, ee); + aa = opacity.AbsorptionCoefficient(rho, temp, ee); + } // reset collision indicators bool is_absorbed = false; From 3643bc9a086990615d179a2ac9ab36f43b893912 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 20 Dec 2024 15:48:45 -0700 Subject: [PATCH 11/16] Modified inputs... but stepdiff is nan... --- inputs/inf_stiff.in | 1 + inputs/stepdiff.in | 1 + inputs/stepdiff_ddmc.in | 1 + inputs/stepdiff_smr.in | 1 + inputs/stepdiff_smr_ddmc.in | 1 + inputs/stepdiff_smr_hybrid.in | 1 + 6 files changed, 6 insertions(+) diff --git a/inputs/inf_stiff.in b/inputs/inf_stiff.in index 520e0f5..7d72f27 100644 --- a/inputs/inf_stiff.in +++ b/inputs/inf_stiff.in @@ -53,6 +53,7 @@ source_strategy = uniform seed = 349856 +frequency_type = gray opacity_model = constant opacity_constant_value = 1000.0 # cm^2/g scattering_model = constant diff --git a/inputs/stepdiff.in b/inputs/stepdiff.in index 4a69394..83beb5d 100644 --- a/inputs/stepdiff.in +++ b/inputs/stepdiff.in @@ -65,6 +65,7 @@ source_strategy = uniform seed = 349857 +frequency_type = gray opacity_model = none scattering_model = constant scattering_constant_value = 1.0e3 diff --git a/inputs/stepdiff_ddmc.in b/inputs/stepdiff_ddmc.in index fdc0734..f86ce3c 100644 --- a/inputs/stepdiff_ddmc.in +++ b/inputs/stepdiff_ddmc.in @@ -66,6 +66,7 @@ source_strategy = uniform seed = 349857 +frequency_type = gray opacity_model = none scattering_model = constant scattering_constant_value = 1.0e3 diff --git a/inputs/stepdiff_smr.in b/inputs/stepdiff_smr.in index 6c82284..9985f2e 100644 --- a/inputs/stepdiff_smr.in +++ b/inputs/stepdiff_smr.in @@ -75,6 +75,7 @@ source_strategy = uniform seed = 349857 +frequency_type = gray opacity_model = none scattering_model = constant scattering_constant_value = 1.0e3 diff --git a/inputs/stepdiff_smr_ddmc.in b/inputs/stepdiff_smr_ddmc.in index ba45948..57142ca 100644 --- a/inputs/stepdiff_smr_ddmc.in +++ b/inputs/stepdiff_smr_ddmc.in @@ -77,6 +77,7 @@ source_strategy = uniform seed = 349857 +frequency_type = gray opacity_model = none scattering_model = constant scattering_constant_value = 1.0e3 diff --git a/inputs/stepdiff_smr_hybrid.in b/inputs/stepdiff_smr_hybrid.in index e19ff7a..d7928cc 100644 --- a/inputs/stepdiff_smr_hybrid.in +++ b/inputs/stepdiff_smr_hybrid.in @@ -76,6 +76,7 @@ source_strategy = uniform seed = 349857 +frequency_type = gray opacity_model = none scattering_model = constant scattering_constant_value = 1.0e3 From d264ed015d08c53b8cee8ab3633e832419d75af9 Mon Sep 17 00:00:00 2001 From: Ben Ryan Date: Fri, 20 Dec 2024 15:52:39 -0700 Subject: [PATCH 12/16] OK stepdiff runs now but the test fails --- src/mcblock/mcblock.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mcblock/mcblock.cpp b/src/mcblock/mcblock.cpp index b9ae449..dd48c4b 100644 --- a/src/mcblock/mcblock.cpp +++ b/src/mcblock/mcblock.cpp @@ -111,7 +111,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { auto opac = singularity::photons::Gray(1.e-100); mopacity = singularity::photons::MeanNonCGSUnits( - singularity::photons::MeanOpacityBase(opac, -1, 1, 1, -1, 1, 1), time_scale, + singularity::photons::MeanOpacityBase(opac, -1, 1, 2, -1, 1, 2), time_scale, mass_scale, length_scale, 1.); } else if (opacity_model == "constant") { Real kappa = pin->GetReal("mcblock", "opacity_constant_value"); From f026b921f21aba497a72b10a978a969ea876c40b Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 21 Jan 2025 14:29:09 -0700 Subject: [PATCH 13/16] Apply suggestions from code review + Cancel log and exponential where possible. + Add comment explaining nu formula is log-space group midpoint. --- src/jaybenne/sourcing.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 5a27cb8..fb7d71c 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -116,14 +116,16 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { } else if constexpr (FT == FrequencyType::multigroup) { // Construct emission CDF const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; - Real nu = std::exp(std::log(numind) + 0.5 * dlnu); + // this is the mid-point of the 1st group in log-space: nu=exp(log(numin) + 0.5*dlnu) + Real nu = numind * std::exp(0.5 * dlnu); Real dnu = nu * dlnu; vmesh(b, fj::emission_cdf(0), k, j, i) = opac.EmissivityPerNu(rho, temp, - std::exp(std::log(numind) + 0.5 * dlnu)) * + numind * std::exp(0.5 * dlnu)) * dnu; for (int n = 1; n < n_nubinsd; n++) { - nu = std::exp(std::log(numind) + (n + 0.5) * dlnu); + // this is the mid-point of group n in log-space: nu=exp(log(numin)+(n+0.5)*dlnu) + nu = numind * std::exp((n + 0.5) * dlnu); dnu = dlnu * nu; vmesh(b, fj::emission_cdf(n), k, j, i) = opac.EmissivityPerNu(rho, temp, nu) * dnu + From 8ed31613d2450233b94055d1454340679ba3ce02 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 21 Jan 2025 17:14:35 -0700 Subject: [PATCH 14/16] Make a few minor updates. + Fix the formatting using style/format.sh. + Use new AbsorptionCoefficient API for mean opacity. + Update some parts of the README.md file. Note: similar to MeanOpacity's Emissivity, AbsorptionCoefficient defaults to Rosseland-weighted grey opacity. --- README.md | 5 ++--- src/jaybenne/sourcing.cpp | 9 +++++---- src/jaybenne/transport.cpp | 3 ++- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index b786871..0e71cc1 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,6 @@ Currently supported computers/partitions are: ## Darwin - power9-rhel7 volta-x86 skylake-gold @@ -92,5 +91,5 @@ of `jaybenne`. # Run driver executable - cd build/src - mpiexec -n 1 ./mcblock -i ../../prob/jbinput.stepdiff + cd build + mpiexec -n 1 ./mcblock -i ../inputs/stepdiff.in diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index fb7d71c..0801c1d 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -116,15 +116,16 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { } else if constexpr (FT == FrequencyType::multigroup) { // Construct emission CDF const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; - // this is the mid-point of the 1st group in log-space: nu=exp(log(numin) + 0.5*dlnu) + // this is the mid-point of the 1st group in log-space: + // nu=exp(log(numin) + 0.5*dlnu) Real nu = numind * std::exp(0.5 * dlnu); Real dnu = nu * dlnu; vmesh(b, fj::emission_cdf(0), k, j, i) = - opac.EmissivityPerNu(rho, temp, - numind * std::exp(0.5 * dlnu)) * + opac.EmissivityPerNu(rho, temp, numind * std::exp(0.5 * dlnu)) * dnu; for (int n = 1; n < n_nubinsd; n++) { - // this is the mid-point of group n in log-space: nu=exp(log(numin)+(n+0.5)*dlnu) + // this is the mid-point of group n in log-space: + // nu=exp(log(numin)+(n+0.5)*dlnu) nu = numind * std::exp((n + 0.5) * dlnu); dnu = dlnu * nu; vmesh(b, fj::emission_cdf(n), k, j, i) = diff --git a/src/jaybenne/transport.cpp b/src/jaybenne/transport.cpp index 8b59777..dcb6651 100644 --- a/src/jaybenne/transport.cpp +++ b/src/jaybenne/transport.cpp @@ -142,8 +142,9 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d Real ss; Real aa; if constexpr (FT == FrequencyType::gray) { + // TODO: use TotalScatteringCoefficient(rho, temp), when available ss = mscattering.RosselandMeanTotalScatteringCoefficient(rho, temp); - aa = mopacity.RosselandMeanAbsorptionCoefficient(rho, temp); + aa = mopacity.AbsorptionCoefficient(rho, temp); } else if constexpr (FT == FrequencyType::multigroup) { ss = scattering.TotalScatteringCoefficient(rho, temp, ee); aa = opacity.AbsorptionCoefficient(rho, temp, ee); From 33d68b95cc75a6d420bf75724effd22e3a3210ec Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 21 Jan 2025 17:57:37 -0700 Subject: [PATCH 15/16] Cancel more exp(log()) occurrences. --- src/jaybenne/sourcing.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 0801c1d..3921bf5 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -134,7 +134,7 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { } emis = 0.; for (int n = 0; n < n_nubinsd; n++) { - const Real dnu = dlnu * std::exp(std::log(numind) + (n + 0.5) * dlnu); + const Real dnu = dlnu * numind * std::exp((n + 0.5) * dlnu); // Get total emissivity emis += vmesh(b, fj::emission_cdf(n), k, j, i); // Normalize emission CDF @@ -254,7 +254,7 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { } } const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; - const Real nu = std::exp(std::log(numind) + (n + 0.5) * dlnu); + const Real nu = numind * std::exp((n + 0.5) * dlnu); ppack_r(b, ph::energy(), n) = hd * nu; } ppack_r(b, ph::weight(), n) = From 5bffda07d5d8b08a18202e5e7fa5d8004b099a67 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Wed, 29 Jan 2025 12:28:36 -0700 Subject: [PATCH 16/16] Do not divide particle energy weight by frequency. + Remove division by frequency in source particle creation. + Update singularity-opac commit to correct mean opacity accesses. + Have regressions fail if mean errors are NaN. --- external/singularity-opac | 2 +- src/jaybenne/sourcing.cpp | 5 +++-- tst/regression_test.py | 4 ++-- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/external/singularity-opac b/external/singularity-opac index b61be52..b05c6e1 160000 --- a/external/singularity-opac +++ b/external/singularity-opac @@ -1 +1 @@ -Subproject commit b61be5226fee505d05ffeb6419408740792b51c7 +Subproject commit b05c6e12f003edae4cea38f6aa279c395aa11098 diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 3921bf5..733c4f3 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -223,6 +223,9 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { ppack_i(b, ph::ijk(1), n) = j; ppack_i(b, ph::ijk(2), n) = k; + // Set energy weight + ppack_r(b, ph::weight(), n) = vmesh(b, fj::source_ew_per_cell(), k, j, i); + // Sample position uniformly in space over cell // TODO(BRR) only valid for Cartesian ppack_r(b, swarm_position::x(), n) = xi + dx_i * (rng_gen.drand() - 0.5); @@ -257,8 +260,6 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { const Real nu = numind * std::exp((n + 0.5) * dlnu); ppack_r(b, ph::energy(), n) = hd * nu; } - ppack_r(b, ph::weight(), n) = - vmesh(b, fj::source_ew_per_cell(), k, j, i) / ppack_r(b, ph::energy(), n); if constexpr (ST == SourceType::emission) { dejbn -= ppack_r(b, ph::weight(), n); diff --git a/tst/regression_test.py b/tst/regression_test.py index 47aae00..981757d 100644 --- a/tst/regression_test.py +++ b/tst/regression_test.py @@ -412,10 +412,10 @@ def analytic_comparison( print(f"Max fractional error: {max_frac_error:.2e}") if args.comparison == "mean": - if mean_frac_error > tolerance: + if mean_frac_error > tolerance or np.isnan(mean_frac_error): success = False elif args.comparison == "weighted_mean": - if mean_frac_error_weighted > tolerance: + if mean_frac_error_weighted > tolerance or np.isnan(mean_frac_error): success = False if args.visualize: