From 6baaf6fb37680a0ae324feb90441cc837177e71f Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Thu, 19 Mar 2026 16:16:55 -0600 Subject: [PATCH 01/15] Implement preliminary multigroup (MG) version of DDMC. + Add ddmc_mg_utils.hpp for MG-specific DDMC functionality. + Add MG DDMC integration and sampling routines and structs. + Add midpoint rule integral for non-dimensional Planck function. + Enable multigroup template in transport_ddmc. + Use DDMC->IMC out-scatter and DDMC leakage frequency sampling. + Use integrals at face over group for total leakage proabilities. + Replace ddmc_face_prob with ddmc_(lo|hi)_face_prob. Note: ddmc_face_prob is only symmetric at face for grey DDMC, so ddmc_(lo|hi)_face_prob for low and high side of face are needed. --- src/jaybenne/ddmc_mg_utils.hpp | 291 ++++++++++++++++++++++++++++ src/jaybenne/jaybenne.cpp | 210 +++++++++++++++----- src/jaybenne/jaybenne.hpp | 6 +- src/jaybenne/jaybenne_variables.hpp | 3 +- src/jaybenne/model_enums.hpp | 25 +++ src/jaybenne/planck.hpp | 14 ++ src/jaybenne/sample_ddmc_bface.cpp | 67 +++++-- src/jaybenne/scattering.hpp | 72 ++++++- src/jaybenne/transport.cpp | 44 ++--- src/jaybenne/transport_ddmc.cpp | 197 +++++++++++++++++-- src/jaybenne/transport_utils.hpp | 13 +- src/mcblock/mcblock.cpp | 1 - 12 files changed, 819 insertions(+), 124 deletions(-) create mode 100644 src/jaybenne/ddmc_mg_utils.hpp create mode 100644 src/jaybenne/model_enums.hpp diff --git a/src/jaybenne/ddmc_mg_utils.hpp b/src/jaybenne/ddmc_mg_utils.hpp new file mode 100644 index 0000000..fe2f9d7 --- /dev/null +++ b/src/jaybenne/ddmc_mg_utils.hpp @@ -0,0 +1,291 @@ +#ifndef JAYBENNE_DDMC_MG_UTILS_HPP_ +#define JAYBENNE_DDMC_MG_UTILS_HPP_ + +// Host configuration file +#include "planck.hpp" + +// ddmc_mg_utils +namespace jaybenne { + +//---------------------------------------------------------------------------------------- +// helper struct to encapsulate data needed to integrate DDMC probabilities over groups +struct ddmc_mg_leak_args { + const int &n_nubins; // number of frequency groups (or bins) + const Real &numind; // minimum frequency + const Real &dlnu; // log-spacing of frequency groups + const Real &hd; // Planck constant + const Real &sbd; // Stefan-Boltzmann constant + const Real &tau_ddmc; // DDMC cell optical-thickness threshold + const Real &dx_l; // lower cell length + const Real &dx_u; // upper cell length + const Real &rho_l; // lower cell density + const Real &rho_u; // upper cell density + const Real &temp_l; // lower cell temperature + const Real &temp_u; // upper cell temperature +}; + +// helper struct for in-cell MG DDMC processes +struct ddmc_mg_cell_args { + const int &n_nubins; // number of frequency groups (or bins) + const Real &numind; // minimum frequency + const Real &dlnu; // log-spacing of frequency groups + const Real &hd; // Planck constant + const Real &sbd; // Stefan-Boltzmann constant + const Real &tau_ddmc; // DDMC cell optical-thickness threshold + const Real &dx_min; // min cell length scale used to activate DDMC + const Real ρ // cell density + const Real &temp; // cell temperature +}; + +//---------------------------------------------------------------------------------------- +template +KOKKOS_FORCEINLINE_FUNCTION std::pair +calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args &dmg, + const bool &use_lo) { + + // define extrapolation distance (Habetler & Matkowski 1975) + constexpr Real lam_ext = 0.7104; + + // initialize unnnormalized Planck integral and probability + Real planck_sum = 0.0; + Real leak_sum = 0.0; + + // integrate + for (int n = 0; n < dmg.n_nubins; ++n) { + + // this is the mid-point of group n in log-space: + // nu=exp(log(numin)+(n+0.5)*dlnu) + const Real nu = dmg.numind * std::exp((n + 0.5) * dmg.dlnu); + const Real dnu = nu * dmg.dlnu; + + // evaluate face opacities at nu + const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, nu); + const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, nu); + const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, nu); + const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, nu); + + // calculate optical thicknesses from lower and upper cell + const Real tau_l = dmg.dx_l * (ss_l + aa_l); + const Real tau_u = dmg.dx_u * (ss_u + aa_u); + + // check if group is included + const bool use_grp = (use_lo ? tau_l > dmg.tau_ddmc : tau_u > dmg.tau_ddmc); + + if (use_grp) { + const Real mtau_l = tau_l > dmg.tau_ddmc ? tau_l : 2.0 * lam_ext; + const Real mtau_u = tau_u > dmg.tau_ddmc ? tau_u : 2.0 * lam_ext; + + // calculate per-nu-bin leakage + const Real Pg = 2.0 / (3.0 * (mtau_l + mtau_u)); + + // evaluate a face temperature (note this is typically a 4-norm average) + const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); + + // get an unnormalized, non-dimensional Planck integral over group + const Real bg = midpoint_Planck(dmg.sbd * temp_f, dmg.hd * nu, dmg.hd * dnu); + + // aggregate values + planck_sum += bg; + leak_sum += bg * Pg; + } + } + + PARTHENON_DEBUG_REQUIRE(planck_sum > 0.0, "Planck integral <= 0.0"); + return {leak_sum, planck_sum}; +} + +//---------------------------------------------------------------------------------------- +// integrate unnormalized leakage probability on high or low side of face +template +KOKKOS_FORCEINLINE_FUNCTION Real calc_ddmc_mg_leakprob(const OP &abs, const SC &sct, + const ddmc_mg_leak_args &dmg, + const bool &use_lo) { + const auto leak_sum_pair = calc_ddmc_mg_leak_numdenom(abs, sct, dmg, use_lo); + // normalize sum (so that leakage probability is an average) + return leak_sum_pair.first / leak_sum_pair.second; +} + +//---------------------------------------------------------------------------------------- +template +KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &sct, + const ddmc_mg_leak_args &dmg, + const bool &use_lo, + RngGen &rng_gen) { + + // first get CDF totals + const auto leak_sum_pair = calc_ddmc_mg_leak_numdenom(abs, sct, dmg, use_lo); + const Real &leak_tot_sum = leak_sum_pair.first; + const Real &planck_tot_sum = leak_sum_pair.second; + + // define extrapolation distance (Habetler & Matkowski 1975) + constexpr Real lam_ext = 0.7104; + + Real leak_sum = 0.0; + Real planck_sum = 0.0; + + // sample + const Real rand1 = leak_tot_sum * rng_gen.drand(); + Real nu_sampled = -1.0; // poisoned initialization + + // integrate + for (int n = 0; n < dmg.n_nubins; ++n) { + + // this is the mid-point of group n in log-space: + // nu=exp(log(numin)+(n+0.5)*dlnu) + const Real nu = dmg.numind * std::exp((n + 0.5) * dmg.dlnu); + const Real dnu = nu * dmg.dlnu; + + // evaluate face opacities at nu + const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, nu); + const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, nu); + const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, nu); + const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, nu); + + // calculate optical thicknesses from lower and upper cell + const Real tau_l = dmg.dx_l * (ss_l + aa_l); + const Real tau_u = dmg.dx_u * (ss_u + aa_u); + + // check if group is included + const bool use_grp = (use_lo ? tau_l > dmg.tau_ddmc : tau_u > dmg.tau_ddmc); + + if (use_grp) { + const Real mtau_l = tau_l > dmg.tau_ddmc ? tau_l : 2.0 * lam_ext; + const Real mtau_u = tau_u > dmg.tau_ddmc ? tau_u : 2.0 * lam_ext; + + // calculate per-nu-bin leakage + const Real Pg = 2.0 / (3.0 * (mtau_l + mtau_u)); + + // evaluate a face temperature (note this is typically a 4-norm average) + const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); + + // get an unnormalized, non-dimensional Planck integral over group + const Real bg = midpoint_Planck(dmg.sbd * temp_f, dmg.hd * nu, dmg.hd * nu); + + // aggregate values + planck_sum += bg; + leak_sum += bg * Pg; + if (leak_sum > rand1) { + nu_sampled = nu; + break; + } + } + } + + PARTHENON_DEBUG_REQUIRE(nu_sampled > 0.0, "nu_sampled <= 0.0"); + + // return sampled energy value (in units of frequency) + return dmg.hd * nu_sampled; +} + +//---------------------------------------------------------------------------------------- +// calculate: absorption, outscatter probability +template +KOKKOS_FORCEINLINE_FUNCTION std::pair +calc_ddmc_mg_probs(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc) { + + // initialize unnnormalized Planck integral and probability + Real planck_sum = 0.0; + Real abs_sum = 0.0; + Real abs_tot_sum = 0.0; + + // integrate + for (int n = 0; n < dmgc.n_nubins; ++n) { + + // this is the mid-point of group n in log-space: + // nu=exp(log(numin)+(n+0.5)*dlnu) + const Real nu = dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); + const Real dnu = nu * dmgc.dlnu; + + // evaluate face opacities at nu + const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, nu); + const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu); + + // get an unnormalized, non-dimensional Planck integral over group + const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, dmgc.hd * nu, dmgc.hd * dnu); + + // sum total (Planck) + abs_tot_sum += bg * aa; + + // check group inclusion + if (dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc) { + + // aggregate values + planck_sum += bg; + abs_sum += bg * aa; + } + } + + // 1st entry = DDMC absorption, 2nd = out-scatter probability + return {abs_sum / planck_sum, 1.0 - abs_sum / abs_tot_sum}; +} + +//---------------------------------------------------------------------------------------- +// sample out-scatter IMC group +// NOTE(MGDDMC): this routine is assuming Kirchhoff's Law for emissivity (LTE) +template +KOKKOS_FORCEINLINE_FUNCTION Real sample_ddmc2imc_outscatter(const OP &abs, const SC &sct, + const ddmc_mg_cell_args &dmgc, + RngGen &rng_gen) { + + Real scat_out_tot_sum = 0.0; + + // integrate total cdf value + for (int n = 0; n < dmgc.n_nubins; ++n) { + + // this is the mid-point of group n in log-space: + // nu=exp(log(numin)+(n+0.5)*dlnu) + const Real nu = dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); + const Real dnu = nu * dmgc.dlnu; + + // evaluate face opacities at nu + const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, nu); + const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu); + + // check group exclusion + if (!(dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc)) { + // get an unnormalized, non-dimensional Planck integral over group + const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, dmgc.hd * nu, dmgc.hd * dnu); + scat_out_tot_sum += bg * aa; + } + } + + // sample + const Real rand1 = scat_out_tot_sum * rng_gen.drand(); + Real nu_sampled = -1.0; // poisoned initialization + + Real abs_sum = 0.0; + + // find bin for sample + for (int n = 0; n < dmgc.n_nubins; ++n) { + + // this is the mid-point of group n in log-space: + // nu=exp(log(numin)+(n+0.5)*dlnu) + const Real nu = dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); + const Real dnu = nu * dmgc.dlnu; + + // evaluate face opacities at nu + const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, nu); + const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu); + + // check group exclusion + if (!(dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc)) { + + // get an unnormalized, non-dimensional Planck integral over group + const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, dmgc.hd * nu, dmgc.hd * dnu); + + // aggregate values + abs_sum += bg * aa; + + if (abs_sum > rand1) { + nu_sampled = nu; + break; + } + } + } + + return dmgc.hd * nu_sampled; +} + +} // namespace jaybenne + +#endif // JAYBENNE_DDMC_MG_UTILS_HPP_ diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 7ac2186..36feb4f 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -15,6 +15,7 @@ #include // Jaybenne includes +#include "ddmc_mg_utils.hpp" #include "jaybenne.hpp" #include "jaybenne_utils.hpp" #include @@ -87,7 +88,8 @@ TaskCollection RadiationStep(Mesh *pmesh, const SimTime &tm, const Real dt) { const auto &fd = jb_pkg->template Param("frequency_type"); // MeshData subsets - auto ddmc_field_names = std::vector{fj::ddmc_face_prob::name()}; + auto ddmc_field_names = std::vector{fj::ddmc_lo_face_prob::name(), + fj::ddmc_hi_face_prob::name()}; auto &ddmc_reg = pmesh->mesh_data.AddShallow("ddmc_reg", pmesh->mesh_data.Get(), ddmc_field_names); @@ -324,7 +326,8 @@ Initialize_impl(ParameterInput *pin, EOS &eos, // Face-based radiation fields Metadata mface({Metadata::Face, Metadata::Derived, Metadata::FillGhost}); - pkg->AddField(field::jaybenne::ddmc_face_prob::name(), mface); + pkg->AddField(field::jaybenne::ddmc_lo_face_prob::name(), mface); + pkg->AddField(field::jaybenne::ddmc_hi_face_prob::name(), mface); // Radiation timestep pkg->EstimateTimestepMesh = EstimateTimestepMesh; @@ -429,6 +432,12 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { Scattering scattering; MeanOpacity mopacity; MeanScattering mscattering; + int n_nubins = -1; + Real numin = -1.0; + Real numax = -1.0; + Real dlnu = -1.0; + Real h = -1.0; + Real sb = -1.0; // get opacity average indicators const auto &use_planck = jbn->template Param("use_planck"); @@ -444,6 +453,13 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { "Modified Fleck factor is not compatible with multigroup!"); opacity = jbn->template Param("opacity_d"); scattering = jbn->template Param("scattering_d"); + n_nubins = jbn->template Param("n_nubins"); + numin = jbn->template Param("numin"); + numax = jbn->template Param("numax"); + h = jbn->template Param("planck_constant"); + sb = jbn->template Param("stefan_boltzmann"); + // initialize (assumed) log spacing + dlnu = (std::log(numax) - std::log(numin)) / n_nubins; } const auto &ib = md->GetBoundsI(IndexDomain::interior); @@ -451,8 +467,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { const auto &kb = md->GetBoundsK(IndexDomain::interior); static auto desc = - MakePackDescriptor( - resolved_pkgs.get()); + MakePackDescriptor(resolved_pkgs.get()); auto vmesh = desc.GetPack(md); parthenon::par_for( @@ -497,8 +513,6 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { const bool use_ddmc = jbn->template Param("use_ddmc"); if (use_ddmc) { - PARTHENON_REQUIRE(FT == FrequencyType::gray, "DDMC only works in gray!"); - // use Planck for all-Planck mode in leakage coefficients too const OpacityAveraging gmode2 = (use_planck && !use_rosseland) ? Planck : Rosseland; @@ -551,27 +565,59 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] auto mscatter = mscattering; [[maybe_unused]] auto opac = opacity; [[maybe_unused]] auto scatter = scattering; + [[maybe_unused]] const auto numind = numin; + [[maybe_unused]] const auto dlnud = dlnu; + [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto hd = h; + [[maybe_unused]] const auto sbd = sb; if constexpr (FT == FrequencyType::gray) { ss_l = mscatter.RosselandMeanTotalScatteringCoefficient(rho_l, temp_l); aa_l = mopac.AbsorptionCoefficient(rho_l, temp_l, gmode2d); ss_u = mscatter.RosselandMeanTotalScatteringCoefficient(rho_u, temp_u); aa_u = mopac.AbsorptionCoefficient(rho_u, temp_u, gmode2d); - } else if constexpr (FT == FrequencyType::multigroup) { - // TODO: replace 3rd argument when this routine operates in multigroup - ss_l = scatter.TotalScatteringCoefficient(rho_l, temp_l, 1.0); - aa_l = opac.AbsorptionCoefficient(rho_l, temp_l, 1.0); - ss_u = scatter.TotalScatteringCoefficient(rho_u, temp_u, 1.0); - aa_u = opac.AbsorptionCoefficient(rho_u, temp_u, 1.0); - } - // calculate optical thicknesses from lower and upper cell - Real tau_l = dx_lx * (ss_l + aa_l); - Real tau_u = dx_ux * (ss_u + aa_u); - tau_l = tau_l > tau_ddmc ? tau_l : 2.0 * lam_ext; - tau_u = tau_u > tau_ddmc ? tau_u : 2.0 * lam_ext; + // calculate optical thicknesses from lower and upper cell + Real tau_l = dx_lx * (ss_l + aa_l); + Real tau_u = dx_ux * (ss_u + aa_u); + tau_l = tau_l > tau_ddmc ? tau_l : 2.0 * lam_ext; + tau_u = tau_u > tau_ddmc ? tau_u : 2.0 * lam_ext; + + // set probability (face DDMC albedo); for grey mode these are copies + vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), k, j, i) = + 2.0 / (3.0 * (tau_l + tau_u)); + vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), k, j, i) = + 2.0 / (3.0 * (tau_l + tau_u)); + + } else if constexpr (FT == FrequencyType::multigroup) { - // set probability (face DDMC albedo) - vmesh(b, TE::F1, fj::ddmc_face_prob(), k, j, i) = 2.0 / (3.0 * (tau_l + tau_u)); + // create DDMC MG leakage data helper argument (defined in ddmc_mg_utils.hpp) + // clang-format off + const ddmc_mg_leak_args dmg{n_nubinsd, + numind, + dlnud, + hd, + sbd, + tau_ddmc, + dx_lx, + dx_ux, + rho_l, + rho_u, + temp_l, + temp_u}; + // clang-format on + + // face side indicator + constexpr bool use_lo = true; + constexpr bool use_hi = false; + + // integrate lo-x leakage probability + vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), k, j, i) = + calc_ddmc_mg_leakprob(opac, scatter, dmg, use_lo); + + // integrate hi-x leakage probability + vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), k, j, i) = + calc_ddmc_mg_leakprob(opac, scatter, dmg, use_hi); + } }); // set face probabilities in X2 direction @@ -618,28 +664,60 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] auto mscatter = mscattering; [[maybe_unused]] auto opac = opacity; [[maybe_unused]] auto scatter = scattering; + [[maybe_unused]] const auto numind = numin; + [[maybe_unused]] const auto dlnud = dlnu; + [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto hd = h; + [[maybe_unused]] const auto sbd = sb; if constexpr (FT == FrequencyType::gray) { ss_l = mscatter.RosselandMeanTotalScatteringCoefficient(rho_l, temp_l); aa_l = mopac.AbsorptionCoefficient(rho_l, temp_l, gmode2d); ss_u = mscatter.RosselandMeanTotalScatteringCoefficient(rho_u, temp_u); aa_u = mopac.AbsorptionCoefficient(rho_u, temp_u, gmode2d); - } else if constexpr (FT == FrequencyType::multigroup) { - // TODO: replace 3rd argument when this routine operates in multigroup - ss_l = scatter.TotalScatteringCoefficient(rho_l, temp_l, 1.0); - aa_l = opac.AbsorptionCoefficient(rho_l, temp_l, 1.0); - ss_u = scatter.TotalScatteringCoefficient(rho_u, temp_u, 1.0); - aa_u = opac.AbsorptionCoefficient(rho_u, temp_u, 1.0); - } - // calculate optical thicknesses from lower and upper cell - Real tau_l = dx_ly * (ss_l + aa_l); - Real tau_u = dx_uy * (ss_u + aa_u); - tau_l = tau_l > tau_ddmc ? tau_l : 2.0 * lam_ext; - tau_u = tau_u > tau_ddmc ? tau_u : 2.0 * lam_ext; + // calculate optical thicknesses from lower and upper cell + Real tau_l = dx_ly * (ss_l + aa_l); + Real tau_u = dx_uy * (ss_u + aa_u); + tau_l = tau_l > tau_ddmc ? tau_l : 2.0 * lam_ext; + tau_u = tau_u > tau_ddmc ? tau_u : 2.0 * lam_ext; - // set probability (face DDMC albedo) - vmesh(b, TE::F2, fj::ddmc_face_prob(), k, j, i) = - 2.0 / (3.0 * (tau_l + tau_u)); + // set probability (face DDMC albedo); for grey mode these are copies + vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), k, j, i) = + 2.0 / (3.0 * (tau_l + tau_u)); + vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), k, j, i) = + 2.0 / (3.0 * (tau_l + tau_u)); + + } else if constexpr (FT == FrequencyType::multigroup) { + + // create DDMC MG leakage data helper argument (defined in + // ddmc_mg_utils.hpp) + // clang-format off + const ddmc_mg_leak_args dmg{n_nubins, + numind, + dlnud, + hd, + sbd, + tau_ddmc, + dx_ly, + dx_uy, + rho_l, + rho_u, + temp_l, + temp_u}; + // clang-format on + + // face side indicator + constexpr bool use_lo = true; + constexpr bool use_hi = false; + + // integrate lo-y leakage probability + vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), k, j, i) = + calc_ddmc_mg_leakprob(opac, scatter, dmg, use_lo); + + // integrate hi-y leakage probability + vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), k, j, i) = + calc_ddmc_mg_leakprob(opac, scatter, dmg, use_hi); + } }); } @@ -687,28 +765,60 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] auto mscatter = mscattering; [[maybe_unused]] auto opac = opacity; [[maybe_unused]] auto scatter = scattering; + [[maybe_unused]] const auto numind = numin; + [[maybe_unused]] const auto dlnud = dlnu; + [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto hd = h; + [[maybe_unused]] const auto sbd = sb; if constexpr (FT == FrequencyType::gray) { ss_l = mscatter.RosselandMeanTotalScatteringCoefficient(rho_l, temp_l); aa_l = mopac.AbsorptionCoefficient(rho_l, temp_l, gmode2d); ss_u = mscatter.RosselandMeanTotalScatteringCoefficient(rho_u, temp_u); aa_u = mopac.AbsorptionCoefficient(rho_u, temp_u, gmode2d); - } else if constexpr (FT == FrequencyType::multigroup) { - // TODO: replace 3rd argument when this routine operates in multigroup - ss_l = scatter.TotalScatteringCoefficient(rho_l, temp_l, 1.0); - aa_l = opac.AbsorptionCoefficient(rho_l, temp_l, 1.0); - ss_u = scatter.TotalScatteringCoefficient(rho_u, temp_u, 1.0); - aa_u = opac.AbsorptionCoefficient(rho_u, temp_u, 1.0); - } - // calculate optical thicknesses from lower and upper cell - Real tau_l = dx_lz * (ss_l + aa_l); - Real tau_u = dx_uz * (ss_u + aa_u); - tau_l = tau_l > tau_ddmc ? tau_l : 2.0 * lam_ext; - tau_u = tau_u > tau_ddmc ? tau_u : 2.0 * lam_ext; + // calculate optical thicknesses from lower and upper cell + Real tau_l = dx_lz * (ss_l + aa_l); + Real tau_u = dx_uz * (ss_u + aa_u); + tau_l = tau_l > tau_ddmc ? tau_l : 2.0 * lam_ext; + tau_u = tau_u > tau_ddmc ? tau_u : 2.0 * lam_ext; - // set probability (face DDMC albedo) - vmesh(b, TE::F3, fj::ddmc_face_prob(), k, j, i) = - 2.0 / (3.0 * (tau_l + tau_u)); + // set probability (face DDMC albedo); for grey mode these are copies + vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), k, j, i) = + 2.0 / (3.0 * (tau_l + tau_u)); + vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), k, j, i) = + 2.0 / (3.0 * (tau_l + tau_u)); + + } else if constexpr (FT == FrequencyType::multigroup) { + + // create DDMC MG leakage data helper argument (defined in + // ddmc_mg_utils.hpp) + // clang-format off + const ddmc_mg_leak_args dmg{n_nubins, + numind, + dlnud, + hd, + sbd, + tau_ddmc, + dx_lz, + dx_uz, + rho_l, + rho_u, + temp_l, + temp_u}; + // clang-format on + + // face side indicator + constexpr bool use_lo = true; + constexpr bool use_hi = false; + + // integrate lo-z leakage probability + vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), k, j, i) = + calc_ddmc_mg_leakprob(opac, scatter, dmg, use_lo); + + // integrate hi-z leakage probability + vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), k, j, i) = + calc_ddmc_mg_leakprob(opac, scatter, dmg, use_hi); + } }); } } diff --git a/src/jaybenne/jaybenne.hpp b/src/jaybenne/jaybenne.hpp index 1b79e12..fc3cadc 100644 --- a/src/jaybenne/jaybenne.hpp +++ b/src/jaybenne/jaybenne.hpp @@ -43,6 +43,7 @@ using namespace parthenon::package::prelude; // TODO(BRR) don't include these here? #include "boundaries.hpp" #include "jaybenne_variables.hpp" +#include "model_enums.hpp" #include "planck.hpp" #include "scattering.hpp" @@ -55,11 +56,6 @@ std::shared_ptr Initialize(parthenon::ParameterInput *pin, MeanOpacity &mopacity, MeanScattering &mscattering, EOS &eos, std::string block_name = "jaybenne"); -// Model enums -enum class SourceStrategy { uniform, energy }; -enum class SourceType { thermal, emission }; -enum class FrequencyType { gray, multigroup }; - // Initialization nulls template KOKKOS_FORCEINLINE_FUNCTION constexpr auto JaybenneNull() { diff --git a/src/jaybenne/jaybenne_variables.hpp b/src/jaybenne/jaybenne_variables.hpp index 2f7599c..172b104 100644 --- a/src/jaybenne/jaybenne_variables.hpp +++ b/src/jaybenne/jaybenne_variables.hpp @@ -34,7 +34,8 @@ namespace field { namespace jaybenne { JAYBENNE_FIELD_VARIABLE(field.jaybenne, energy_tally); JAYBENNE_FIELD_VARIABLE(field.jaybenne, fleck_factor); -JAYBENNE_FIELD_VARIABLE(field.jaybenne, ddmc_face_prob); +JAYBENNE_FIELD_VARIABLE(field.jaybenne, ddmc_lo_face_prob); +JAYBENNE_FIELD_VARIABLE(field.jaybenne, ddmc_hi_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, delta_num_per_cell); diff --git a/src/jaybenne/model_enums.hpp b/src/jaybenne/model_enums.hpp new file mode 100644 index 0000000..30e0e6c --- /dev/null +++ b/src/jaybenne/model_enums.hpp @@ -0,0 +1,25 @@ +//======================================================================================== +// (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 JAYBENNE_MODEL_ENUMS_HPP_ +#define JAYBENNE_MODEL_ENUMS_HPP_ + +namespace jaybenne { + +// Model enums +enum class SourceStrategy { uniform, energy }; +enum class SourceType { thermal, emission }; +enum class FrequencyType { gray, multigroup }; + +} // namespace jaybenne + +#endif // JAYBENNE_MODEL_ENUMS_HPP_ diff --git a/src/jaybenne/planck.hpp b/src/jaybenne/planck.hpp index dc1709e..4861bea 100644 --- a/src/jaybenne/planck.hpp +++ b/src/jaybenne/planck.hpp @@ -18,6 +18,20 @@ namespace jaybenne { +//---------------------------------------------------------------------------------------- +//! \fn Midpoint Rule Planck integral, returning unnormalized non-dimensional value +//! \brief Note the arguments must be in the same units (code temperature units). +//! The onus is on the calling routine(s) to normalize or dimensionalize if needed. +//! Since the midpoint is provided, this could be made 1st-order as well +//! (left/right) +KOKKOS_FORCEINLINE_FUNCTION +Real midpoint_Planck(const Real &temp, const Real &nu, const Real &dnu) { + const Real tinv = 1.0 / temp; + const Real x = nu * tinv; + const Real efac = std::exp(-x); + return (dnu * tinv) * (x * x * x) * efac / (1.0 - efac); +} + //---------------------------------------------------------------------------------------- //! \fn Real sample_Planck_energy //! \brief Efficiently samples the Planck distribution for particle energy diff --git a/src/jaybenne/sample_ddmc_bface.cpp b/src/jaybenne/sample_ddmc_bface.cpp index ffadb12..289a040 100644 --- a/src/jaybenne/sample_ddmc_bface.cpp +++ b/src/jaybenne/sample_ddmc_bface.cpp @@ -100,7 +100,8 @@ TaskStatus SampleDDMCBlockFace(MeshData *md) { const bool three_d = (pm->ndim > 2); // Create SparsePack - static auto desc = MakePackDescriptor(resolved_pkgs.get()); + static auto desc = MakePackDescriptor( + resolved_pkgs.get()); auto vmesh = desc.GetPack(md); // set tolerance for checking particle coordinate @@ -184,9 +185,13 @@ TaskStatus SampleDDMCBlockFace(MeshData *md) { // get face probabilities for bounding faces const Real &Px_uy = - vmesh(b, TE::F1, fj::ddmc_face_prob(), kp, jp_u, ip_b); + at_block_x_min + ? vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), kp, jp_u, ip_b) + : vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), kp, jp_u, ip_b); const Real &Px_ly = - vmesh(b, TE::F1, fj::ddmc_face_prob(), kp, jp_l, ip_b); + at_block_x_min + ? vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), kp, jp_l, ip_b) + : vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), kp, jp_l, ip_b); SampleFace2D(jp_l, dy_j, Px_ly, Px_uy, rng_gen, jp, y); } @@ -213,9 +218,13 @@ TaskStatus SampleDDMCBlockFace(MeshData *md) { // get face probabilities for bounding faces const Real &Py_ux = - vmesh(b, TE::F2, fj::ddmc_face_prob(), kp, jp_b, ip_u); + at_block_y_min + ? vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), kp, jp_b, ip_u) + : vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), kp, jp_b, ip_u); const Real &Py_lx = - vmesh(b, TE::F2, fj::ddmc_face_prob(), kp, jp_b, ip_l); + at_block_y_min + ? vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), kp, jp_b, ip_l) + : vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), kp, jp_b, ip_l); SampleFace2D(ip_l, dx_i, Py_lx, Py_ux, rng_gen, ip, x); } @@ -311,13 +320,21 @@ TaskStatus SampleDDMCBlockFace(MeshData *md) { // get face probabilities for bounding faces const Real &Px_uu = - vmesh(b, TE::F1, fj::ddmc_face_prob(), kp_u, jp_u, ip_b); + at_block_x_min + ? vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), kp_u, jp_u, ip_b) + : vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), kp_u, jp_u, ip_b); const Real &Px_ul = - vmesh(b, TE::F1, fj::ddmc_face_prob(), kp_u, jp_l, ip_b); + at_block_x_min + ? vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), kp_u, jp_l, ip_b) + : vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), kp_u, jp_l, ip_b); const Real &Px_lu = - vmesh(b, TE::F1, fj::ddmc_face_prob(), kp_l, jp_u, ip_b); + at_block_x_min + ? vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), kp_l, jp_u, ip_b) + : vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), kp_l, jp_u, ip_b); const Real &Px_ll = - vmesh(b, TE::F1, fj::ddmc_face_prob(), kp_l, jp_l, ip_b); + at_block_x_min + ? vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), kp_l, jp_l, ip_b) + : vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), kp_l, jp_l, ip_b); // sample a refined yz-face SampleFace3D(jp_l, kp_l, dy_j, dz_k, Px_ll, Px_lu, Px_ul, Px_uu, @@ -350,13 +367,21 @@ TaskStatus SampleDDMCBlockFace(MeshData *md) { // get face probabilities for bounding faces const Real &Py_uu = - vmesh(b, TE::F2, fj::ddmc_face_prob(), kp_u, jp_b, ip_u); + at_block_y_min + ? vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), kp_u, jp_b, ip_u) + : vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), kp_u, jp_b, ip_u); const Real &Py_ul = - vmesh(b, TE::F2, fj::ddmc_face_prob(), kp_u, jp_b, ip_l); + at_block_y_min + ? vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), kp_u, jp_b, ip_l) + : vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), kp_u, jp_b, ip_l); const Real &Py_lu = - vmesh(b, TE::F2, fj::ddmc_face_prob(), kp_l, jp_b, ip_u); + at_block_y_min + ? vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), kp_l, jp_b, ip_u) + : vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), kp_l, jp_b, ip_u); const Real &Py_ll = - vmesh(b, TE::F2, fj::ddmc_face_prob(), kp_l, jp_b, ip_l); + at_block_y_min + ? vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), kp_l, jp_b, ip_l) + : vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), kp_l, jp_b, ip_l); // sample a refined zx-face SampleFace3D(ip_l, kp_l, dx_i, dz_k, Py_ll, Py_lu, Py_ul, Py_uu, @@ -389,13 +414,21 @@ TaskStatus SampleDDMCBlockFace(MeshData *md) { // get face probabilities for bounding faces const Real &Pz_uu = - vmesh(b, TE::F3, fj::ddmc_face_prob(), kp_b, jp_u, ip_u); + at_block_z_min + ? vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), kp_b, jp_u, ip_u) + : vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), kp_b, jp_u, ip_u); const Real &Pz_ul = - vmesh(b, TE::F3, fj::ddmc_face_prob(), kp_b, jp_u, ip_l); + at_block_z_min + ? vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), kp_b, jp_u, ip_l) + : vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), kp_b, jp_u, ip_l); const Real &Pz_lu = - vmesh(b, TE::F3, fj::ddmc_face_prob(), kp_b, jp_l, ip_u); + at_block_z_min + ? vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), kp_b, jp_l, ip_u) + : vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), kp_b, jp_l, ip_u); const Real &Pz_ll = - vmesh(b, TE::F3, fj::ddmc_face_prob(), kp_b, jp_l, ip_l); + at_block_z_min + ? vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), kp_b, jp_l, ip_l) + : vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), kp_b, jp_l, ip_l); // sample a refined xy-face SampleFace3D(ip_l, jp_l, dx_i, dy_j, Pz_ll, Pz_lu, Pz_ul, Pz_uu, diff --git a/src/jaybenne/scattering.hpp b/src/jaybenne/scattering.hpp index 7d3be39..7358755 100644 --- a/src/jaybenne/scattering.hpp +++ b/src/jaybenne/scattering.hpp @@ -16,17 +16,71 @@ namespace jaybenne { //---------------------------------------------------------------------------------------- -//! \fn void scatter -//! \brief TODO(RTW): template on ScatteringModel, when pertinent +// scattering input data helper structs +struct cell_scat_args { + const int &b; // block index + const int &ip; // x/X1 block cell index + const int &jp; // y/X2 block cell index + const int &kp; // z/X3 block cell index + const Real &ff; // Fleck factor + const Real &aa; // absorption opacity (1/length) + const Real &ss; // scattering opacity (1/length) + // frequency data + const int &n_nubinsd; // number of frequency groups + const Real &numind; // mininmum frequency + const Real &numaxd; // maximum frequency + const Real &hd; // Planck constant +}; +struct ptcl_scat_args { + RngGen &rng_gen; + const Real &vv; + Real &vx; // particle x/X1-direction speed + Real &vy; // particle y/X2-direction speed + Real &vz; // particle z/X3-direction speed + Real ⅇ // particle frequency (in units of energy) +}; + +//---------------------------------------------------------------------------------------- +//! \fn void isotropic direction redistribution KOKKOS_FORCEINLINE_FUNCTION -void ScatterKernel(RngGen &rng_gen, const Real &vv, Real &vx_out, Real &vy_out, - Real &vz_out) { - const Real mu = 2.0 * rng_gen.drand() - 1.0; - const Real phi = 2.0 * M_PI * rng_gen.drand(); +void sample_vol_iso_dir(ptcl_scat_args sa) { + const Real mu = 2.0 * sa.rng_gen.drand() - 1.0; + const Real phi = 2.0 * M_PI * sa.rng_gen.drand(); const Real stheta = std::sqrt(1.0 - mu * mu); - vx_out = vv * stheta * std::cos(phi); - vy_out = vv * stheta * std::sin(phi); - vz_out = vv * mu; + sa.vx = sa.vv * stheta * std::cos(phi); + sa.vy = sa.vv * stheta * std::sin(phi); + sa.vz = sa.vv * mu; +} + +//---------------------------------------------------------------------------------------- +//! \fn void scatter kernel +//! TODO(BRR): template on scattering kernel type? +template +KOKKOS_FORCEINLINE_FUNCTION void scatter_kernel(const T &vmesh, cell_scat_args csa, + ptcl_scat_args psa) { + namespace fj = field::jaybenne; + + // sample direction isotropically (see TODO above about more general kernels) + sample_vol_iso_dir(psa); + + // sample whether effective scattering occurred + const Real rand1 = psa.rng_gen.drand(); + if (rand1 * ((1.0 - csa.ff) * csa.aa + csa.ss) < (1.0 - csa.ff) * csa.aa) { + + // Sample energy from CDF + const Real rand2 = psa.rng_gen.drand(); + int n; + for (n = 0; n < csa.n_nubinsd; n++) { + if (vmesh(csa.b, fj::emission_cdf(n), csa.kp, csa.jp, csa.ip) >= rand2) { + break; + } + } + const Real dlnu = (std::log(csa.numaxd) - std::log(csa.numind)) / csa.n_nubinsd; + const Real nu = csa.numind * std::exp((n + 0.5) * dlnu); + + // reset particle frequency + psa.ee = csa.hd * nu; + } } } // namespace jaybenne diff --git a/src/jaybenne/transport.cpp b/src/jaybenne/transport.cpp index 1925e85..c268ee0 100644 --- a/src/jaybenne/transport.cpp +++ b/src/jaybenne/transport.cpp @@ -150,7 +150,6 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d [[maybe_unused]] auto scatter = scattering; [[maybe_unused]] auto eost = eos; if constexpr (FT == FrequencyType::gray) { - // TODO: use TotalScatteringCoefficient(rho, temp), when available ss = vmesh(b, fjh::scattering_opacity(), kp, jp, ip); aa = vmesh(b, fjh::absorption_opacity(), kp, jp, ip); } else if constexpr (FT == FrequencyType::multigroup) { @@ -221,28 +220,27 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d } if (is_scattered) { - // process scattering - // TODO(BRR): template on scattering model - ScatterKernel(rng_gen, vv, vx, vy, vz); - - // if multigroup eff scatter, redistribute frequency - if constexpr (FT == FrequencyType::multigroup) { - // sample whether effective scattering occurred - const Real rand1 = rng_gen.drand(); - if (rand1 * ((1.0 - ff) * aa + ss) < (1.0 - ff) * aa) { - - // Sample energy from CDF - const Real rand2 = rng_gen.drand(); - int n; - for (n = 0; n < n_nubinsd; n++) { - if (vmesh(b, fj::emission_cdf(n), kp, jp, ip) >= rand2) { - break; - } - } - const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; - const Real nu = numind * std::exp((n + 0.5) * dlnu); - ee = hd * nu; - } + + // form particle scattering argument struct + ptcl_scat_args psa{rng_gen, vv, vx, vy, vz, ee}; + + if constexpr (FT == FrequencyType::gray) { + + // just do direction-sampling + sample_vol_iso_dir(psa); + + // if multigroup eff scatter, redistribute frequency + } else if constexpr (FT == FrequencyType::multigroup) { + + // form cell scattering argument struct + // clang-format off + cell_scat_args csa{b, ip, jp, kp, + ff, aa, ss, + n_nubinsd, numind, numaxd, hd}; + // clang-format on + + // invoke frequency-dependent scattering kernel + scatter_kernel(vmesh, csa, psa); } } diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index be28bb6..ac3ec39 100644 --- a/src/jaybenne/transport_ddmc.cpp +++ b/src/jaybenne/transport_ddmc.cpp @@ -15,6 +15,7 @@ #include // Jaybenne includes +#include "ddmc_mg_utils.hpp" #include "jaybenne.hpp" #include "jaybenne_utils.hpp" #include "scattering.hpp" @@ -34,20 +35,40 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R namespace ph = particle::photons; using TE = parthenon::TopologicalElement; - PARTHENON_REQUIRE(FT == FrequencyType::gray, "DDMC only works in gray!"); - auto pm = md->GetParentPointer(); auto &resolved_pkgs = pm->resolved_packages; auto &jb_pkg = pm->packages.Get("jaybenne"); + auto &eos = jb_pkg->template Param("eos_d"); auto &rng_pool = jb_pkg->template Param("rng_pool"); const Real vv = jb_pkg->template Param("speed_of_light"); const Real ske = 0.5 * SQR(vv); const Real &tau_ddmc = jb_pkg->template Param("tau_ddmc"); const Real cutoff = jb_pkg->template Param("cutoff"); + // data needed for multigroup frequency sampling + const Real h = jb_pkg->template Param("planck_constant"); + const Real sb = jb_pkg->template Param("stefan_boltzmann"); + int n_nubins = JaybenneNull(); + Real numin = JaybenneNull(); + Real numax = JaybenneNull(); + Real dlnu = JaybenneNull(); + Opacity opacity; + Scattering scattering; + if constexpr (FT == FrequencyType::multigroup) { + n_nubins = jb_pkg->template Param("n_nubins"); + numin = jb_pkg->template Param("numin"); + numax = jb_pkg->template Param("numax"); + // initialize (assumed) log spacing + dlnu = (std::log(numax) - std::log(numin)) / n_nubins; + // set opacity objects + opacity = jb_pkg->template Param("opacity_d"); + scattering = jb_pkg->template Param("scattering_d"); + } + // Create SparsePack static auto desc = - MakePackDescriptor( resolved_pkgs.get()); auto vmesh = desc.GetPack(md); @@ -81,6 +102,15 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R const auto &swarm_d = ppack_r.GetContext(b); if (swarm_d.IsActive(n)) { auto rng_gen = rng_pool.get_state(); + + // frequency data, needed for multigroup + [[maybe_unused]] const auto hd = h; + [[maybe_unused]] const auto sbd = sb; + [[maybe_unused]] const auto numind = numin; + [[maybe_unused]] const auto numaxd = numax; + [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto dlnud = dlnu; + auto &coords = vmesh.GetCoordinates(b); const Real &dx_i = coords.template Dxc(0, 0, 0); const Real &dx_j = coords.template Dxc(0, 0, 0); @@ -131,8 +161,23 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // Extract physical quantities const Real &ff = vmesh(b, fj::fleck_factor(), kp, jp, ip); - const Real &ss = vmesh(b, fjh::scattering_opacity(), kp, jp, ip); - const Real &aa = vmesh(b, fjh::absorption_opacity(), kp, jp, ip); + Real ss = JaybenneNull(); + Real aa = JaybenneNull(); + [[maybe_unused]] Real rho = JaybenneNull(); + [[maybe_unused]] Real temp = JaybenneNull(); + [[maybe_unused]] auto opac = opacity; + [[maybe_unused]] auto scatter = scattering; + [[maybe_unused]] auto eost = eos; + if constexpr (FT == FrequencyType::gray) { + ss = vmesh(b, fjh::scattering_opacity(), kp, jp, ip); + aa = vmesh(b, fjh::absorption_opacity(), kp, jp, ip); + } else if constexpr (FT == FrequencyType::multigroup) { + rho = vmesh(b, fjh::density(), kp, jp, ip); + const Real &sie = vmesh(b, fjh::sie(), kp, jp, ip); + temp = eost.TemperatureFromDensityInternalEnergy(rho, sie); + ss = scatter.TotalScatteringCoefficient(rho, temp, ee); + aa = opac.AbsorptionCoefficient(rho, temp, ee); + } // reset collision indicators bool is_absorbed = false; @@ -156,23 +201,51 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R 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); - const Real &Px_u = vmesh(b, TE::F1, fj::ddmc_face_prob(), kp, jp, ip + 1); + const Real &Px_l = vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), kp, jp, ip); + const Real &Px_u = + vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), kp, jp, ip + 1); const Real &Py_l = - multi_d ? vmesh(b, TE::F2, fj::ddmc_face_prob(), kp, jp, ip) : 0.0; + multi_d ? vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), kp, jp, ip) : 0.0; const Real &Py_u = - multi_d ? vmesh(b, TE::F2, fj::ddmc_face_prob(), kp, jp + 1, ip) : 0.0; + multi_d ? vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), kp, jp + 1, ip) + : 0.0; const Real &Pz_l = - three_d ? vmesh(b, TE::F3, fj::ddmc_face_prob(), kp, jp, ip) : 0.0; + three_d ? vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), kp, jp, ip) : 0.0; const Real &Pz_u = - three_d ? vmesh(b, TE::F3, fj::ddmc_face_prob(), kp + 1, jp, ip) : 0.0; + three_d ? vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), kp + 1, jp, ip) + : 0.0; + + // store old cell indices (this is only needed for multigroup) + const int ip_old = ip; + const int jp_old = jp; + const int kp_old = kp; + + // NOTE(MGDDMC): grey ss will be needed for inelastic scattering + Real aa_g = aa; + Real gm_g = 0.0; + if constexpr (FT == FrequencyType::multigroup) { + // clang-format off + ddmc_mg_cell_args dmgc{n_nubinsd, + numind, + dlnud, + hd, + sbd, + tau_ddmc, + dx_push, + rho, + temp}; + // clang-format on + const auto aagm = calc_ddmc_mg_probs(opac, scatter, dmgc); + aa_g = aagm.first; + gm_g = aagm.second; + } // create DDMC step argument list // clang-format off ddmc_step_args dia{ // constants rng_gen, t_start, dt, - ff, aa, ss, vv, + ff, aa_g, ss, gm_g, vv, multi_d, three_d, xl, yl, zl, xu, yu, zu, Px_l, Py_l, Pz_l, Px_u, Py_u, Pz_u, @@ -187,6 +260,64 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R if (!is_rejected) ptcl_ddmc_step(dia, cutoff); + if constexpr (FT == FrequencyType::multigroup) { + // NOTE: these ddmc_mg_leak_args do not use adjacent cell, so + // the face CDF does not sum here to ddmc_(hi|lo)_face_prob. + // This is hopefully a minor error. + + // particle must have leaked if nothing else + if (!(is_absorbed || is_scattered || is_census || is_rejected)) { + + // only one index should be +/-1 of the current index + const int ip_u = (ip_old == ip + 1 ? ip_old : ip); + const int ip_l = (ip_old == ip - 1 ? ip_old : ip); + const int jp_u = (jp_old == jp + 1 ? jp_old : jp); + const int jp_l = (jp_old == jp - 1 ? jp_old : jp); + const int kp_u = (kp_old == kp + 1 ? kp_old : kp); + const int kp_l = (kp_old == kp - 1 ? kp_old : kp); + PARTHENON_DEBUG_REQUIRE(ip_u + jp_u + kp_u - ip_l - jp_l - kp_l == 1, + "invalid index difference for DDMC leakage"); + + // select side of face + const bool use_lo_x = (ip_old == ip - 1); + const bool use_lo_y = (jp_old == jp - 1); + const bool use_lo_z = (kp_old == kp - 1); + // use_lo can only be false here if one of the old is the new index+1 + const bool use_lo = (use_lo_x || use_lo_y || use_lo_z); + + // TODO(MGDDMC): is this dx alone sufficient for nu-sampling at face? + const Real dx_f = (ip_old != ip ? dx_i : (jp_old != jp ? dx_j : dx_k)); + + // get rho and temperature + const Real &rho_l = vmesh(b, fjh::density(), kp_l, jp_l, ip_l); + const Real &sie_l = vmesh(b, fjh::sie(), kp_l, jp_l, ip_l); + const Real temp_l = + eos.TemperatureFromDensityInternalEnergy(rho_l, sie_l); + const Real &rho_u = vmesh(b, fjh::density(), kp_u, jp_u, ip_u); + const Real &sie_u = vmesh(b, fjh::sie(), kp_u, jp_u, ip_u); + const Real temp_u = + eos.TemperatureFromDensityInternalEnergy(rho_u, sie_u); + + // clang-format off + const ddmc_mg_leak_args dmg{n_nubinsd, + numind, + dlnud, + hd, + sbd, + tau_ddmc, + dx_f, + dx_f, + rho_l, + rho_u, + temp_l, + temp_u}; + // clang-format off + + // sample particle frequency (energy units) + ee = sample_leakage_group(opac, scatter, dmg, use_lo, rng_gen); + } + } + } else { // push particle @@ -253,10 +384,44 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R } if (is_scattered) { - // process scattering - // TODO(BRR): if eff scatter, redistribute frequency - // TODO(BRR): template on scattering model - ScatterKernel(rng_gen, vv, vx, vy, vz); + + // form particle scattering argument struct + ptcl_scat_args psa{rng_gen, vv, vx, vy, vz, ee}; + + // if multigroup eff scatter, redistribute frequency + if constexpr (FT == FrequencyType::gray) { + + // just do direction-sampling (only called by IMC, by construction) + sample_vol_iso_dir(psa); + + } else if constexpr (FT == FrequencyType::multigroup) { + + if (is_ddmc_step) { + // clang-format off + ddmc_mg_cell_args dmgc{n_nubinsd, + numind, + dlnud, + hd, + sbd, + tau_ddmc, + dx_push, + rho, + temp}; + // clang-format on + ee = sample_ddmc2imc_outscatter(opac, scatter, dmgc, rng_gen); + } else { + + // form cell scattering argument struct + // clang-format off + cell_scat_args csa{b, ip, jp, kp, + ff, aa, ss, + n_nubinsd, numind, numaxd, hd}; + // clang-format on + + // invoke frequency-dependent scattering kernel + scatter_kernel(vmesh, csa, psa); + } + } } if (is_census) { diff --git a/src/jaybenne/transport_utils.hpp b/src/jaybenne/transport_utils.hpp index cbc07ec..ddd50cf 100644 --- a/src/jaybenne/transport_utils.hpp +++ b/src/jaybenne/transport_utils.hpp @@ -85,6 +85,7 @@ struct ddmc_step_args { const Real &ff; // Fleck factor const Real &aa; // absorption opacity (1/length) const Real &ss; // scattering opacity (1/length) + const Real &gm; // out-scatter probability (unitless) const Real &vv; // particle speed (should be c) const bool &multi_d; // 2D or 3D const bool &three_d; // 3D @@ -255,8 +256,11 @@ void ptcl_ddmc_step(ddmc_step_args dia, const double cutoff) { // attenuate if fraction >= cutoff, analog absorb if fraction < cutoff const Real an_abs = (dia.fraction < cutoff) ? dia.ff * dia.aa : 0.0; + // effective out-scatter probability (gm=0 for grey DDMC) + const Real sct_out = dia.gm * (1.0 - dia.ff) * dia.aa; + // calculate time to DDMC event and compare to time to end of time step (census) - const Real cdf_ddmc = an_abs + leak_tot + rmin; + const Real cdf_ddmc = an_abs + sct_out + leak_tot + rmin; const Real dt_ddmc = -std::log(dia.rng_gen.drand()) / (dia.vv * cdf_ddmc); const Real dt_end = (dia.t_start + dia.dt) - dia.t; const bool is_ddmc_event = dt_ddmc < dt_end; @@ -283,7 +287,12 @@ void ptcl_ddmc_step(ddmc_step_args dia, const double cutoff) { // particle will be absorbed dia.is_absorbed = true; - } else if (xi < an_abs + leak_tot) { + } else if (xi < an_abs + sct_out) { + + // particle will be scattered into an IMC group + dia.is_scattered = true; + + } else if (xi < an_abs + sct_out + leak_tot) { // TODO(RTW): only sample direction if adjacent cell is below tau_ddmc diff --git a/src/mcblock/mcblock.cpp b/src/mcblock/mcblock.cpp index fc110c4..e87a108 100644 --- a/src/mcblock/mcblock.cpp +++ b/src/mcblock/mcblock.cpp @@ -72,7 +72,6 @@ std::shared_ptr Initialize(ParameterInput *pin) { // opacity fields m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}); - // TODO: maybe worth filling ghosts for expedited DDMC stencil pkg->AddField(field::material::absorption_opacity::name(), m); pkg->AddField(field::material::scattering_opacity::name(), m); From 7afe8133e98bd1fadbe927e9511bb92a53710c82 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 7 Apr 2026 12:01:48 -0600 Subject: [PATCH 02/15] Make DDMC threshold check consistent among direction-based stencils. --- src/jaybenne/ddmc_mg_utils.hpp | 18 ++++++--- src/jaybenne/jaybenne.cpp | 71 +++++++++++++++++++++++++++------ src/jaybenne/transport_ddmc.cpp | 2 + 3 files changed, 73 insertions(+), 18 deletions(-) diff --git a/src/jaybenne/ddmc_mg_utils.hpp b/src/jaybenne/ddmc_mg_utils.hpp index fe2f9d7..ff09ad0 100644 --- a/src/jaybenne/ddmc_mg_utils.hpp +++ b/src/jaybenne/ddmc_mg_utils.hpp @@ -16,6 +16,8 @@ struct ddmc_mg_leak_args { const Real &hd; // Planck constant const Real &sbd; // Stefan-Boltzmann constant const Real &tau_ddmc; // DDMC cell optical-thickness threshold + const Real &dx_lmin; // min cell length scale used to activate DDMC + const Real &dx_umin; // min cell length scale used to activate DDMC const Real &dx_l; // lower cell length const Real &dx_u; // upper cell length const Real &rho_l; // lower cell density @@ -65,15 +67,17 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, nu); // calculate optical thicknesses from lower and upper cell + const Real tau_lmin = dmg.dx_lmin * (ss_l + aa_l); + const Real tau_umin = dmg.dx_umin * (ss_u + aa_u); const Real tau_l = dmg.dx_l * (ss_l + aa_l); const Real tau_u = dmg.dx_u * (ss_u + aa_u); // check if group is included - const bool use_grp = (use_lo ? tau_l > dmg.tau_ddmc : tau_u > dmg.tau_ddmc); + const bool use_grp = (use_lo ? tau_lmin > dmg.tau_ddmc : tau_umin > dmg.tau_ddmc); if (use_grp) { - const Real mtau_l = tau_l > dmg.tau_ddmc ? tau_l : 2.0 * lam_ext; - const Real mtau_u = tau_u > dmg.tau_ddmc ? tau_u : 2.0 * lam_ext; + const Real mtau_l = tau_lmin > dmg.tau_ddmc ? tau_l : 2.0 * lam_ext; + const Real mtau_u = tau_umin > dmg.tau_ddmc ? tau_u : 2.0 * lam_ext; // calculate per-nu-bin leakage const Real Pg = 2.0 / (3.0 * (mtau_l + mtau_u)); @@ -142,15 +146,17 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &s const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, nu); // calculate optical thicknesses from lower and upper cell + const Real tau_lmin = dmg.dx_lmin * (ss_l + aa_l); + const Real tau_umin = dmg.dx_umin * (ss_u + aa_u); const Real tau_l = dmg.dx_l * (ss_l + aa_l); const Real tau_u = dmg.dx_u * (ss_u + aa_u); // check if group is included - const bool use_grp = (use_lo ? tau_l > dmg.tau_ddmc : tau_u > dmg.tau_ddmc); + const bool use_grp = (use_lo ? tau_lmin > dmg.tau_ddmc : tau_umin > dmg.tau_ddmc); if (use_grp) { - const Real mtau_l = tau_l > dmg.tau_ddmc ? tau_l : 2.0 * lam_ext; - const Real mtau_u = tau_u > dmg.tau_ddmc ? tau_u : 2.0 * lam_ext; + const Real mtau_l = tau_lmin > dmg.tau_ddmc ? tau_l : 2.0 * lam_ext; + const Real mtau_u = tau_umin > dmg.tau_ddmc ? tau_u : 2.0 * lam_ext; // calculate per-nu-bin leakage const Real Pg = 2.0 / (3.0 * (mtau_l + mtau_u)); diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 36feb4f..7fdf4c8 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -531,6 +531,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { // get coordinates of block auto &coords = vmesh.GetCoordinates(b); const Real &dx_i = coords.Dxc(0, 0, 0); + const Real &dx_j = coords.Dxc(0, 0, 0); + const Real &dx_k = coords.Dxc(0, 0, 0); // get current, lower, upper neighbor block levels in x-direction const Real rlev = static_cast(vmesh.GetLevel(b, 0, 0, 0)); @@ -541,9 +543,18 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { ? rlev : static_cast(vmesh.GetLevel(b, 0, 0, 1)); + // calculate neighbor refinement scaling factors + const Real scle_lx = i == ib.s ? std::pow(2.0, rlev - rlev_lx) : 1.0; + const Real scle_ux = i == iu ? std::pow(2.0, rlev - rlev_ux) : 1.0; + // calculate neighbor dx values - const Real dx_lx = i == ib.s ? std::pow(2.0, rlev - rlev_lx) * dx_i : dx_i; - const Real dx_ux = i == iu ? std::pow(2.0, rlev - rlev_ux) * dx_i : dx_i; + const Real dx_lx = scle_lx * dx_i; + const Real dx_ux = scle_ux * dx_i; + + // calculate neighbor min length for min optical thickness (for threshold check) + const Real dx_push = std::min(dx_i, std::min(dx_j, dx_k)); + const Real dx_lmin = scle_lx * dx_push; + const Real dx_umin = scle_ux * dx_push; // TODO: interpolate temperatures to evaluate face opacities? // (If opacity gradients are not large, maybe this is not needed) @@ -577,10 +588,12 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { aa_u = mopac.AbsorptionCoefficient(rho_u, temp_u, gmode2d); // calculate optical thicknesses from lower and upper cell + const Real tau_lmin = dx_lmin * (ss_l + aa_l); + const Real tau_umin = dx_umin * (ss_u + aa_u); Real tau_l = dx_lx * (ss_l + aa_l); Real tau_u = dx_ux * (ss_u + aa_u); - tau_l = tau_l > tau_ddmc ? tau_l : 2.0 * lam_ext; - tau_u = tau_u > tau_ddmc ? tau_u : 2.0 * lam_ext; + tau_l = tau_lmin > tau_ddmc ? tau_l : 2.0 * lam_ext; + tau_u = tau_umin > tau_ddmc ? tau_u : 2.0 * lam_ext; // set probability (face DDMC albedo); for grey mode these are copies vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), k, j, i) = @@ -598,6 +611,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { hd, sbd, tau_ddmc, + dx_lmin, + dx_umin, dx_lx, dx_ux, rho_l, @@ -630,7 +645,9 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { ib.e, KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { // get coordinates of block auto &coords = vmesh.GetCoordinates(b); + const Real &dx_i = coords.Dxc(0, 0, 0); const Real &dx_j = coords.Dxc(0, 0, 0); + const Real &dx_k = coords.Dxc(0, 0, 0); // get current, lower, upper neighbor block levels in x-direction const Real rlev = static_cast(vmesh.GetLevel(b, 0, 0, 0)); @@ -641,9 +658,19 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { ? rlev : static_cast(vmesh.GetLevel(b, 0, 1, 0)); + // calculate neighbor refinement scaling factors + const Real scle_ly = j == jb.s ? std::pow(2.0, rlev - rlev_ly) : 1.0; + const Real scle_uy = j == ju ? std::pow(2.0, rlev - rlev_uy) : 1.0; + // calculate neighbor dx values - const Real dx_ly = j == jb.s ? std::pow(2.0, rlev - rlev_ly) * dx_j : dx_j; - const Real dx_uy = j == ju ? std::pow(2.0, rlev - rlev_uy) * dx_j : dx_j; + const Real dx_ly = scle_ly * dx_j; + const Real dx_uy = scle_uy * dx_j; + + // calculate neighbor min length for min optical thickness (for threshold + // check) + const Real dx_push = std::min(dx_i, std::min(dx_j, dx_k)); + const Real dx_lmin = scle_ly * dx_push; + const Real dx_umin = scle_uy * dx_push; // TODO: interpolate temperatures to evaluate face opacities? // (If opacity gradients are not large, maybe this is not needed) @@ -676,10 +703,12 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { aa_u = mopac.AbsorptionCoefficient(rho_u, temp_u, gmode2d); // calculate optical thicknesses from lower and upper cell + const Real tau_lmin = dx_lmin * (ss_l + aa_l); + const Real tau_umin = dx_umin * (ss_u + aa_u); Real tau_l = dx_ly * (ss_l + aa_l); Real tau_u = dx_uy * (ss_u + aa_u); - tau_l = tau_l > tau_ddmc ? tau_l : 2.0 * lam_ext; - tau_u = tau_u > tau_ddmc ? tau_u : 2.0 * lam_ext; + tau_l = tau_lmin > tau_ddmc ? tau_l : 2.0 * lam_ext; + tau_u = tau_umin > tau_ddmc ? tau_u : 2.0 * lam_ext; // set probability (face DDMC albedo); for grey mode these are copies vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), k, j, i) = @@ -698,6 +727,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { hd, sbd, tau_ddmc, + dx_lmin, + dx_umin, dx_ly, dx_uy, rho_l, @@ -731,6 +762,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { ib.e, KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { // get coordinates of block auto &coords = vmesh.GetCoordinates(b); + const Real &dx_i = coords.Dxc(0, 0, 0); + const Real &dx_j = coords.Dxc(0, 0, 0); const Real &dx_k = coords.Dxc(0, 0, 0); // get current, lower, upper neighbor block levels in x-direction @@ -742,9 +775,19 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { ? rlev : static_cast(vmesh.GetLevel(b, 1, 0, 0)); + // calculate neighbor refinement scaling factors + const Real scle_lz = k == kb.s ? std::pow(2.0, rlev - rlev_lz) : 1.0; + const Real scle_uz = k == ku ? std::pow(2.0, rlev - rlev_uz) : 1.0; + // calculate neighbor dx values - const Real dx_lz = k == kb.s ? std::pow(2.0, rlev - rlev_lz) * dx_k : dx_k; - const Real dx_uz = k == ku ? std::pow(2.0, rlev - rlev_uz) * dx_k : dx_k; + const Real dx_lz = scle_lz * dx_k; + const Real dx_uz = scle_uz * dx_k; + + // calculate neighbor min length for min optical thickness (for threshold + // check) + const Real dx_push = std::min(dx_i, std::min(dx_j, dx_k)); + const Real dx_lmin = scle_lz * dx_push; + const Real dx_umin = scle_uz * dx_push; // TODO: interpolate temperatures to evaluate face opacities? // (If opacity gradients are not large, maybe this is not needed) @@ -777,10 +820,12 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { aa_u = mopac.AbsorptionCoefficient(rho_u, temp_u, gmode2d); // calculate optical thicknesses from lower and upper cell + const Real tau_lmin = dx_lmin * (ss_l + aa_l); + const Real tau_umin = dx_umin * (ss_u + aa_u); Real tau_l = dx_lz * (ss_l + aa_l); Real tau_u = dx_uz * (ss_u + aa_u); - tau_l = tau_l > tau_ddmc ? tau_l : 2.0 * lam_ext; - tau_u = tau_u > tau_ddmc ? tau_u : 2.0 * lam_ext; + tau_l = tau_lmin > tau_ddmc ? tau_l : 2.0 * lam_ext; + tau_u = tau_umin > tau_ddmc ? tau_u : 2.0 * lam_ext; // set probability (face DDMC albedo); for grey mode these are copies vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), k, j, i) = @@ -799,6 +844,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { hd, sbd, tau_ddmc, + dx_lmin, + dx_umin, dx_lz, dx_uz, rho_l, diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index ac3ec39..521eb9e 100644 --- a/src/jaybenne/transport_ddmc.cpp +++ b/src/jaybenne/transport_ddmc.cpp @@ -307,6 +307,8 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R tau_ddmc, dx_f, dx_f, + dx_f, + dx_f, rho_l, rho_u, temp_l, From 4a0d7756cbd630dbaa0939bc533b4600e8da6e1e Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Wed, 8 Apr 2026 17:34:26 -0600 Subject: [PATCH 03/15] Fix some bugs / issues in MG DDMC (caught so far). + Change Stefan-Boltzmann to Boltzmann constant in Planck function. + Use Rayleigh-Jeans approximation for nu/T<1e-4. + Return 0 for leakage probability if no groups contribute. + Sample isotropic direction after DDMC frequency out-scatter. --- src/jaybenne/ddmc_mg_utils.hpp | 6 +++++- src/jaybenne/jaybenne.cpp | 6 ++---- src/jaybenne/planck.hpp | 8 ++++++-- src/jaybenne/transport_ddmc.cpp | 6 +++++- 4 files changed, 18 insertions(+), 8 deletions(-) diff --git a/src/jaybenne/ddmc_mg_utils.hpp b/src/jaybenne/ddmc_mg_utils.hpp index ff09ad0..332bcf7 100644 --- a/src/jaybenne/ddmc_mg_utils.hpp +++ b/src/jaybenne/ddmc_mg_utils.hpp @@ -51,6 +51,7 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args // initialize unnnormalized Planck integral and probability Real planck_sum = 0.0; Real leak_sum = 0.0; + int nnubins_used = 0; // integrate for (int n = 0; n < dmg.n_nubins; ++n) { @@ -76,6 +77,8 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args const bool use_grp = (use_lo ? tau_lmin > dmg.tau_ddmc : tau_umin > dmg.tau_ddmc); if (use_grp) { + nnubins_used++; + const Real mtau_l = tau_lmin > dmg.tau_ddmc ? tau_l : 2.0 * lam_ext; const Real mtau_u = tau_umin > dmg.tau_ddmc ? tau_u : 2.0 * lam_ext; @@ -94,7 +97,7 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args } } - PARTHENON_DEBUG_REQUIRE(planck_sum > 0.0, "Planck integral <= 0.0"); + PARTHENON_REQUIRE(nnubins_used ? planck_sum > 0.0 : true, "Planck integral <= 0.0"); return {leak_sum, planck_sum}; } @@ -105,6 +108,7 @@ KOKKOS_FORCEINLINE_FUNCTION Real calc_ddmc_mg_leakprob(const OP &abs, const SC & const ddmc_mg_leak_args &dmg, const bool &use_lo) { const auto leak_sum_pair = calc_ddmc_mg_leak_numdenom(abs, sct, dmg, use_lo); + if (!(leak_sum_pair.second > 0.0)) return 0.0; // normalize sum (so that leakage probability is an average) return leak_sum_pair.first / leak_sum_pair.second; } diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 7fdf4c8..ade3c96 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -234,6 +234,7 @@ Initialize_impl(ParameterInput *pin, EOS &eos, pkg->AddParam<>("speed_of_light", units.c); pkg->AddParam<>("stefan_boltzmann", units.sb); pkg->AddParam<>("planck_constant", units.h); + pkg->AddParam<>("boltzmann", units.kb); // RNG bool unique_rank_seeds = pin->GetOrAddBoolean(block_name, "unique_rank_seeds", true); @@ -370,9 +371,6 @@ std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacit auto pkg = Initialize_impl(pin, eos, units, block_name); - PARTHENON_REQUIRE(pkg->template Param("use_ddmc") == false, - "DDMC not supported for multigroup currently!"); - // Frequency discretization auto time = units.time; Real numin = pin->GetReal(block_name, "numin"); // in Hz @@ -457,7 +455,7 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { numin = jbn->template Param("numin"); numax = jbn->template Param("numax"); h = jbn->template Param("planck_constant"); - sb = jbn->template Param("stefan_boltzmann"); + sb = jbn->template Param("boltzmann"); // initialize (assumed) log spacing dlnu = (std::log(numax) - std::log(numin)) / n_nubins; } diff --git a/src/jaybenne/planck.hpp b/src/jaybenne/planck.hpp index 4861bea..903f663 100644 --- a/src/jaybenne/planck.hpp +++ b/src/jaybenne/planck.hpp @@ -28,8 +28,12 @@ KOKKOS_FORCEINLINE_FUNCTION Real midpoint_Planck(const Real &temp, const Real &nu, const Real &dnu) { const Real tinv = 1.0 / temp; const Real x = nu * tinv; - const Real efac = std::exp(-x); - return (dnu * tinv) * (x * x * x) * efac / (1.0 - efac); + if (x > 1e-4) { + const Real efac = std::exp(-x); + return (dnu * tinv) * (x * x * x) * efac / (1.0 - efac); + } else { + return (dnu * tinv) * (x * x); + } } //---------------------------------------------------------------------------------------- diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index 521eb9e..423a90c 100644 --- a/src/jaybenne/transport_ddmc.cpp +++ b/src/jaybenne/transport_ddmc.cpp @@ -47,7 +47,7 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // data needed for multigroup frequency sampling const Real h = jb_pkg->template Param("planck_constant"); - const Real sb = jb_pkg->template Param("stefan_boltzmann"); + const Real sb = jb_pkg->template Param("boltzmann"); int n_nubins = JaybenneNull(); Real numin = JaybenneNull(); Real numax = JaybenneNull(); @@ -411,6 +411,10 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R temp}; // clang-format on ee = sample_ddmc2imc_outscatter(opac, scatter, dmgc, rng_gen); + + // resample particle direction + sample_vol_iso_dir(psa); + } else { // form cell scattering argument struct From 642aa676384229272ebc5aefcfb009f738ca6d84 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 21 Apr 2026 16:31:53 -0600 Subject: [PATCH 04/15] Fix several more bugs in multigroup DDMC. + Convert frequency to energy when calling opacity. + Fix nu->dnu (ee->dee) in call to midpoint_Planck in leakage. + Check for 0-absorption opacity in DDMC, if 0 return 0 outscatter. --- src/jaybenne/ddmc_mg_utils.hpp | 83 ++++++++++++++++++---------------- 1 file changed, 45 insertions(+), 38 deletions(-) diff --git a/src/jaybenne/ddmc_mg_utils.hpp b/src/jaybenne/ddmc_mg_utils.hpp index 332bcf7..1007bae 100644 --- a/src/jaybenne/ddmc_mg_utils.hpp +++ b/src/jaybenne/ddmc_mg_utils.hpp @@ -58,14 +58,14 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args // this is the mid-point of group n in log-space: // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real nu = dmg.numind * std::exp((n + 0.5) * dmg.dlnu); - const Real dnu = nu * dmg.dlnu; + const Real ee = dmg.hd * dmg.numind * std::exp((n + 0.5) * dmg.dlnu); + const Real dee = ee * dmg.dlnu; // evaluate face opacities at nu - const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, nu); - const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, nu); - const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, nu); - const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, nu); + const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, ee); + const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, ee); + const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, ee); + const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, ee); // calculate optical thicknesses from lower and upper cell const Real tau_lmin = dmg.dx_lmin * (ss_l + aa_l); @@ -89,7 +89,7 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); // get an unnormalized, non-dimensional Planck integral over group - const Real bg = midpoint_Planck(dmg.sbd * temp_f, dmg.hd * nu, dmg.hd * dnu); + const Real bg = midpoint_Planck(dmg.sbd * temp_f, ee, dee); // aggregate values planck_sum += bg; @@ -133,21 +133,21 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &s // sample const Real rand1 = leak_tot_sum * rng_gen.drand(); - Real nu_sampled = -1.0; // poisoned initialization + Real ee_sampled = -1.0; // poisoned initialization // integrate for (int n = 0; n < dmg.n_nubins; ++n) { // this is the mid-point of group n in log-space: // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real nu = dmg.numind * std::exp((n + 0.5) * dmg.dlnu); - const Real dnu = nu * dmg.dlnu; + const Real ee = dmg.hd * dmg.numind * std::exp((n + 0.5) * dmg.dlnu); + const Real dee = ee * dmg.dlnu; // evaluate face opacities at nu - const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, nu); - const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, nu); - const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, nu); - const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, nu); + const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, ee); + const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, ee); + const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, ee); + const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, ee); // calculate optical thicknesses from lower and upper cell const Real tau_lmin = dmg.dx_lmin * (ss_l + aa_l); @@ -169,22 +169,22 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &s const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); // get an unnormalized, non-dimensional Planck integral over group - const Real bg = midpoint_Planck(dmg.sbd * temp_f, dmg.hd * nu, dmg.hd * nu); + const Real bg = midpoint_Planck(dmg.sbd * temp_f, ee, dee); // aggregate values planck_sum += bg; leak_sum += bg * Pg; if (leak_sum > rand1) { - nu_sampled = nu; + ee_sampled = ee; break; } } } - PARTHENON_DEBUG_REQUIRE(nu_sampled > 0.0, "nu_sampled <= 0.0"); + PARTHENON_DEBUG_REQUIRE(ee_sampled > 0.0, "ee_sampled <= 0.0"); - // return sampled energy value (in units of frequency) - return dmg.hd * nu_sampled; + // return sampled energy value (in units of energy) + return ee_sampled; } //---------------------------------------------------------------------------------------- @@ -203,15 +203,15 @@ calc_ddmc_mg_probs(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc) // this is the mid-point of group n in log-space: // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real nu = dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); - const Real dnu = nu * dmgc.dlnu; + const Real ee = dmgc.hd * dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); + const Real dee = ee * dmgc.dlnu; // evaluate face opacities at nu - const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, nu); - const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu); + const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, ee); + const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, ee); // get an unnormalized, non-dimensional Planck integral over group - const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, dmgc.hd * nu, dmgc.hd * dnu); + const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, ee, dee); // sum total (Planck) abs_tot_sum += bg * aa; @@ -225,8 +225,15 @@ calc_ddmc_mg_probs(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc) } } + PARTHENON_DEBUG_REQUIRE(planck_sum > 0.0, "planck_sum = 0: no DDMC groups in DDMC."); + // 1st entry = DDMC absorption, 2nd = out-scatter probability - return {abs_sum / planck_sum, 1.0 - abs_sum / abs_tot_sum}; + if (abs_tot_sum > 0.0) { + return {abs_sum / planck_sum, 1.0 - abs_sum / abs_tot_sum}; + } else { + // no possible outscatter if there is no absorption/redistribution + return {abs_sum / planck_sum, 0.0}; + } } //---------------------------------------------------------------------------------------- @@ -244,24 +251,24 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_ddmc2imc_outscatter(const OP &abs, const // this is the mid-point of group n in log-space: // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real nu = dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); - const Real dnu = nu * dmgc.dlnu; + const Real ee = dmgc.hd * dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); + const Real dee = ee * dmgc.dlnu; // evaluate face opacities at nu - const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, nu); - const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu); + const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, ee); + const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, ee); // check group exclusion if (!(dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc)) { // get an unnormalized, non-dimensional Planck integral over group - const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, dmgc.hd * nu, dmgc.hd * dnu); + const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, ee, dee); scat_out_tot_sum += bg * aa; } } // sample const Real rand1 = scat_out_tot_sum * rng_gen.drand(); - Real nu_sampled = -1.0; // poisoned initialization + Real ee_sampled = -1.0; // poisoned initialization Real abs_sum = 0.0; @@ -270,30 +277,30 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_ddmc2imc_outscatter(const OP &abs, const // this is the mid-point of group n in log-space: // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real nu = dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); - const Real dnu = nu * dmgc.dlnu; + const Real ee = dmgc.hd * dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); + const Real dee = ee * dmgc.dlnu; // evaluate face opacities at nu - const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, nu); - const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu); + const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, ee); + const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, ee); // check group exclusion if (!(dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc)) { // get an unnormalized, non-dimensional Planck integral over group - const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, dmgc.hd * nu, dmgc.hd * dnu); + const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, ee, dee); // aggregate values abs_sum += bg * aa; if (abs_sum > rand1) { - nu_sampled = nu; + ee_sampled = ee; break; } } } - return dmgc.hd * nu_sampled; + return ee_sampled; } } // namespace jaybenne From 25d4497738541dcb9638ab9e4f759451d7341b98 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Wed, 22 Apr 2026 12:45:06 -0600 Subject: [PATCH 05/15] Fix GPU multigroup runs. + Update singularity-eos,-opac submodules. + Avoid some if constexpr capture build errors with maybe_unused. --- src/jaybenne/jaybenne.cpp | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index ade3c96..dc99f2e 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -579,6 +579,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto hd = h; [[maybe_unused]] const auto sbd = sb; + [[maybe_unused]] const auto tau_ddmcd = tau_ddmc; + [[maybe_unused]] const auto lam_extd = lam_ext; if constexpr (FT == FrequencyType::gray) { ss_l = mscatter.RosselandMeanTotalScatteringCoefficient(rho_l, temp_l); aa_l = mopac.AbsorptionCoefficient(rho_l, temp_l, gmode2d); @@ -590,8 +592,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { const Real tau_umin = dx_umin * (ss_u + aa_u); Real tau_l = dx_lx * (ss_l + aa_l); Real tau_u = dx_ux * (ss_u + aa_u); - tau_l = tau_lmin > tau_ddmc ? tau_l : 2.0 * lam_ext; - tau_u = tau_umin > tau_ddmc ? tau_u : 2.0 * lam_ext; + tau_l = tau_lmin > tau_ddmcd ? tau_l : 2.0 * lam_extd; + tau_u = tau_umin > tau_ddmcd ? tau_u : 2.0 * lam_extd; // set probability (face DDMC albedo); for grey mode these are copies vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), k, j, i) = @@ -608,7 +610,7 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { dlnud, hd, sbd, - tau_ddmc, + tau_ddmcd, dx_lmin, dx_umin, dx_lx, @@ -694,6 +696,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto hd = h; [[maybe_unused]] const auto sbd = sb; + [[maybe_unused]] const auto tau_ddmcd = tau_ddmc; + [[maybe_unused]] const auto lam_extd = lam_ext; if constexpr (FT == FrequencyType::gray) { ss_l = mscatter.RosselandMeanTotalScatteringCoefficient(rho_l, temp_l); aa_l = mopac.AbsorptionCoefficient(rho_l, temp_l, gmode2d); @@ -705,8 +709,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { const Real tau_umin = dx_umin * (ss_u + aa_u); Real tau_l = dx_ly * (ss_l + aa_l); Real tau_u = dx_uy * (ss_u + aa_u); - tau_l = tau_lmin > tau_ddmc ? tau_l : 2.0 * lam_ext; - tau_u = tau_umin > tau_ddmc ? tau_u : 2.0 * lam_ext; + tau_l = tau_lmin > tau_ddmcd ? tau_l : 2.0 * lam_extd; + tau_u = tau_umin > tau_ddmcd ? tau_u : 2.0 * lam_extd; // set probability (face DDMC albedo); for grey mode these are copies vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), k, j, i) = @@ -724,7 +728,7 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { dlnud, hd, sbd, - tau_ddmc, + tau_ddmcd, dx_lmin, dx_umin, dx_ly, @@ -811,6 +815,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto hd = h; [[maybe_unused]] const auto sbd = sb; + [[maybe_unused]] const auto tau_ddmcd = tau_ddmc; + [[maybe_unused]] const auto lam_extd = lam_ext; if constexpr (FT == FrequencyType::gray) { ss_l = mscatter.RosselandMeanTotalScatteringCoefficient(rho_l, temp_l); aa_l = mopac.AbsorptionCoefficient(rho_l, temp_l, gmode2d); @@ -822,8 +828,8 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { const Real tau_umin = dx_umin * (ss_u + aa_u); Real tau_l = dx_lz * (ss_l + aa_l); Real tau_u = dx_uz * (ss_u + aa_u); - tau_l = tau_lmin > tau_ddmc ? tau_l : 2.0 * lam_ext; - tau_u = tau_umin > tau_ddmc ? tau_u : 2.0 * lam_ext; + tau_l = tau_lmin > tau_ddmcd ? tau_l : 2.0 * lam_extd; + tau_u = tau_umin > tau_ddmcd ? tau_u : 2.0 * lam_extd; // set probability (face DDMC albedo); for grey mode these are copies vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), k, j, i) = @@ -841,7 +847,7 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { dlnud, hd, sbd, - tau_ddmc, + tau_ddmcd, dx_lmin, dx_umin, dx_lz, From 74d48a8addaca08ddea3f4bc0660b5364e8d3668 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Wed, 17 Jun 2026 18:09:41 -0600 Subject: [PATCH 06/15] + Update singularity-opac and parse power-law opacity in mcblock. + Separate mcblock absorption/scattering input nodes (Artemis-like). --- analysis/plot.py | 27 ++++++++++++++++- external/singularity-opac | 2 +- inputs/inf.in | 12 +++++--- inputs/inf_mg.in | 10 ++++-- inputs/inf_stiff.in | 12 +++++--- inputs/stepdiff.in | 10 ++++-- inputs/stepdiff_ddmc.in | 10 ++++-- inputs/stepdiff_smr.in | 10 ++++-- inputs/stepdiff_smr_ddmc.in | 10 ++++-- inputs/stepdiff_smr_hybrid.in | 10 ++++-- src/mcblock/mcblock.cpp | 57 ++++++++++++++++++++++------------- src/mcblock/opacity.hpp | 2 ++ 12 files changed, 123 insertions(+), 49 deletions(-) diff --git a/analysis/plot.py b/analysis/plot.py index 94748eb..d921ea0 100644 --- a/analysis/plot.py +++ b/analysis/plot.py @@ -23,7 +23,16 @@ # Plot a 1D profile of a variable def plot_1d( - fig, ax, filenames, variable_name, draw_meshblocks, vmin, vmax, coords, scale + fig, + ax, + filenames, + variable_name, + draw_meshblocks, + vmin, + vmax, + data_vbnd, + coords, + scale, ): for filename in filenames: @@ -39,6 +48,11 @@ def plot_1d( if scale == "log": variable = np.log10(variable) + if data_vbnd: + for b in range(dump.NumBlocks): + vmin = min(vmin, np.min(variable[b, idx_k, idx_j, :])) + vmax = max(vmax, np.max(variable[b, idx_k, idx_j, :])) + for b in range(dump.NumBlocks): ax.plot(dump.xc[b, idx_k, idx_j, :], variable[b, idx_k, idx_j, :]) @@ -54,6 +68,7 @@ def plot_2d( draw_meshblocks, vmin, vmax, + data_vbnd, coords, scale, ): @@ -70,6 +85,11 @@ def plot_2d( if scale == "log": variable = np.log10(variable) + if data_vbnd: + for b in range(dump.NumBlocks): + vmin = min(vmin, np.min(variable[b, idx_k, :, :])) + vmax = max(vmax, np.max(variable[b, idx_k, :, :])) + for b in range(dump.NumBlocks): ax.pcolormesh( dump.xc[b, idx_k, :, :], @@ -161,6 +181,9 @@ def plot_2d( parser.add_argument( "--vmax", type=float, default=0, help="Maximum value of colorbar" ) + parser.add_argument( + "--data_vbnd", action="store_true", help="Use min/max data for colorbar" + ) parser.add_argument( "--scale", type=str, @@ -187,6 +210,7 @@ def plot_2d( args.meshblocks, args.vmin, args.vmax, + args.data_vbnd, args.coords, args.scale, ) @@ -201,6 +225,7 @@ def plot_2d( args.meshblocks, args.vmin, args.vmax, + args.data_vbnd, args.coords, args.scale, ) diff --git a/external/singularity-opac b/external/singularity-opac index cdd365f..60b7411 160000 --- a/external/singularity-opac +++ b/external/singularity-opac @@ -1 +1 @@ -Subproject commit cdd365ffd56ccdab1ae4b4410d276d0a46a4a314 +Subproject commit 60b7411fae78175c26e71e6dfa804bd1b1cc0e1d diff --git a/inputs/inf.in b/inputs/inf.in index 970c5ce..9897be6 100644 --- a/inputs/inf.in +++ b/inputs/inf.in @@ -59,15 +59,19 @@ seed = 349857 frequency_type = gray -opacity_model = constant -opacity_constant_value = 1.0 # cm^2/g -scattering_model = constant -scattering_constant_value = 1.0e5 specific_heat = 1.0e8 # erg/K/g initial_density = 1.0 # g cm^-3 initial_temperature = 1.0e5 # K initial_radiation = thermal + +opacity_model = constant +constant_value = 1.0 # cm^2/g + + +opacity_model = constant +constant_value = 1.0e5 + file_type = hdf5 dt = 0.1 diff --git a/inputs/inf_mg.in b/inputs/inf_mg.in index 583d5d9..fc50ce0 100644 --- a/inputs/inf_mg.in +++ b/inputs/inf_mg.in @@ -54,14 +54,18 @@ seed = 349851 frequency_type = multigroup -opacity_model = ep_bremss -scattering_model = constant -scattering_constant_value = 1e-1 specific_heat = 1.0e8 # erg/K/g initial_density = 1.0e-5 # g cm^-3 initial_temperature = 1.0e5 # K initial_radiation = thermal + +opacity_model = ep_bremss + + +scattering_model = constant +scattering_constant_value = 1e-1 + file_type = hdf5 dt = 0.1 diff --git a/inputs/inf_stiff.in b/inputs/inf_stiff.in index e0adf6a..636576e 100644 --- a/inputs/inf_stiff.in +++ b/inputs/inf_stiff.in @@ -54,15 +54,19 @@ seed = 349856 frequency_type = gray -opacity_model = constant -opacity_constant_value = 1000.0 # cm^2/g -scattering_model = constant -scattering_constant_value = 0.0 #1.0e5 specific_heat = 1.0e7 # erg/K/g initial_density = 1.0 # g cm^-3 initial_temperature = 1.0e5 # K initial_radiation = thermal + +opacity_model = constant +constant_value = 1000.0 # cm^2/g + + +opacity_model = constant +constant_value = 0.0 #1.0e5 + file_type = hdf5 dt = 1e-11 diff --git a/inputs/stepdiff.in b/inputs/stepdiff.in index 83beb5d..f243756 100644 --- a/inputs/stepdiff.in +++ b/inputs/stepdiff.in @@ -66,14 +66,18 @@ seed = 349857 frequency_type = gray -opacity_model = none -scattering_model = constant -scattering_constant_value = 1.0e3 specific_heat = 1.0e8 # erg/K/g initial_density = 1.0 # g cm^-3 initial_temperature = 1.0e5 # K initial_radiation = thermal + +opacity_model = none + + +opacity_model = constant +constant_value = 1.0e3 + file_type = hdf5 dt = 3.335641e-10 diff --git a/inputs/stepdiff_ddmc.in b/inputs/stepdiff_ddmc.in index f86ce3c..508619b 100644 --- a/inputs/stepdiff_ddmc.in +++ b/inputs/stepdiff_ddmc.in @@ -67,14 +67,18 @@ seed = 349857 frequency_type = gray -opacity_model = none -scattering_model = constant -scattering_constant_value = 1.0e3 specific_heat = 1.0e8 # erg/K/g initial_density = 1.0 # g cm^-3 initial_temperature = 1.0e5 # K initial_radiation = thermal + +opacity_model = none + + +opacity_model = constant +constant_value = 1.0e3 + file_type = hdf5 dt = 3.335641e-10 diff --git a/inputs/stepdiff_smr.in b/inputs/stepdiff_smr.in index 9985f2e..fcbbcb0 100644 --- a/inputs/stepdiff_smr.in +++ b/inputs/stepdiff_smr.in @@ -76,14 +76,18 @@ seed = 349857 frequency_type = gray -opacity_model = none -scattering_model = constant -scattering_constant_value = 1.0e3 specific_heat = 1.0e8 # erg/K/g initial_density = 1.0 # g cm^-3 initial_temperature = 1.0e5 # K initial_radiation = thermal + +opacity_model = none + + +opacity_model = constant +constant_value = 1.0e3 + file_type = hdf5 dt = 3.335641e-10 diff --git a/inputs/stepdiff_smr_ddmc.in b/inputs/stepdiff_smr_ddmc.in index 57142ca..a1ce0c7 100644 --- a/inputs/stepdiff_smr_ddmc.in +++ b/inputs/stepdiff_smr_ddmc.in @@ -78,14 +78,18 @@ seed = 349857 frequency_type = gray -opacity_model = none -scattering_model = constant -scattering_constant_value = 1.0e3 specific_heat = 1.0e8 # erg/K/g initial_density = 1.0 # g cm^-3 initial_temperature = 1.0e5 # K initial_radiation = thermal + +opacity_model = none + + +opacity_model = constant +constant_value = 1.0e3 + file_type = hdf5 dt = 3.335641e-10 diff --git a/inputs/stepdiff_smr_hybrid.in b/inputs/stepdiff_smr_hybrid.in index f5ee961..b3d2dd1 100644 --- a/inputs/stepdiff_smr_hybrid.in +++ b/inputs/stepdiff_smr_hybrid.in @@ -78,14 +78,18 @@ diagnostic_level = 1 frequency_type = gray -opacity_model = none -scattering_model = constant -scattering_constant_value = 1.0e3 specific_heat = 1.0e8 # erg/K/g initial_density = 1.0 # g cm^-3 initial_temperature = 1.0e5 # K initial_radiation = thermal + +opacity_model = none + + +opacity_model = constant +constant_value = 1.0e3 + file_type = hdf5 dt = 3.335641e-10 diff --git a/src/mcblock/mcblock.cpp b/src/mcblock/mcblock.cpp index e87a108..2473760 100644 --- a/src/mcblock/mcblock.cpp +++ b/src/mcblock/mcblock.cpp @@ -110,23 +110,23 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Absorption opacity model Opacity opacity; MeanOpacity mopacity; - std::string opacity_model = pin->GetString("mcblock", "opacity_model"); + std::string abs_model = pin->GetString("mcblock/absorption", "opacity_model"); if (frequency_type == FrequencyType::gray) { - if (opacity_model == "none") { + if (abs_model == "none") { auto opac = singularity::photons::Gray(1.e-100); mopacity = singularity::photons::MeanNonCGSUnits( singularity::photons::MeanOpacityBase(opac, -1, 1, 2, -1, 1, 2), time_scale, mass_scale, length_scale, temperature_scale); - } else if (opacity_model == "constant") { - Real kappa = pin->GetReal("mcblock", "opacity_constant_value"); + } else if (abs_model == "constant") { + Real kappa = pin->GetReal("mcblock/absorption", "constant_value"); auto opac = singularity::photons::Gray(kappa); mopacity = singularity::photons::MeanNonCGSUnits( singularity::photons::MeanOpacityBase(opac, -1, 1, 2, -1, 1, 2), time_scale, mass_scale, length_scale, temperature_scale); - } else if (opacity_model == "table") { - std::string table_filename = pin->GetString("mcblock", "opacity_table"); + } else if (abs_model == "table") { + std::string table_filename = pin->GetString("mcblock/absorption", "opacity_table"); mopacity = singularity::photons::MeanNonCGSUnits( singularity::photons::MeanOpacityBase(table_filename), time_scale, @@ -135,16 +135,25 @@ std::shared_ptr Initialize(ParameterInput *pin) { PARTHENON_FAIL("Only none, constant, or table opacity models supported!"); } } else if (frequency_type == FrequencyType::multigroup) { - if (opacity_model == "none") { + if (abs_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"); + } else if (abs_model == "constant") { + Real kappa = pin->GetReal("mcblock/absorption", "constant_value"); opacity = singularity::photons::NonCGSUnits( singularity::photons::Gray(kappa), time_scale, mass_scale, length_scale, temperature_scale); - } else if (opacity_model == "ep_bremss") { + } else if (abs_model == "powerlaw") { + Real kappa0 = pin->GetReal("mcblock/absorption", "kappa0"); + Real rho_exp = pin->GetReal("mcblock/absorption", "rho_exp"); + Real temp_exp = pin->GetReal("mcblock/absorption", "temp_exp"); + Real nu_exp = pin->GetReal("mcblock/absorption", "nu_exp"); + Real nu_ref = pin->GetReal("mcblock/absorption", "nu_ref"); + opacity = singularity::photons::NonCGSUnits( + singularity::photons::PowerLaw(kappa0, rho_exp, temp_exp, nu_exp, nu_ref), + time_scale, mass_scale, length_scale, temperature_scale); + } else if (abs_model == "ep_bremss") { opacity = singularity::photons::NonCGSUnits( singularity::photons::EPBremss(), time_scale, mass_scale, length_scale, temperature_scale); @@ -154,7 +163,6 @@ std::shared_ptr Initialize(ParameterInput *pin) { } } pkg->AddParam<>("frequency_type", frequency_type); - pkg->AddParam<>("opacity_model", opacity_model); pkg->AddParam<>("opacity_h", opacity); pkg->AddParam<>("mopacity_h", mopacity); @@ -162,20 +170,19 @@ std::shared_ptr Initialize(ParameterInput *pin) { // 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); Scattering scattering; MeanScattering mscattering; - std::string scattering_model = - pin->GetOrAddString("mcblock", "scattering_model", "none"); + std::string sct_model = + pin->GetOrAddString("mcblock/scattering", "opacity_model", "none"); if (frequency_type == FrequencyType::gray) { - if (scattering_model == "none") { + if (sct_model == "none") { auto sopac = singularity::photons::GrayS(0.0, apm); mscattering = singularity::photons::MeanNonCGSUnitsS( 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"); + } else if (sct_model == "constant") { + Real kappa_s = pin->GetReal("mcblock/scattering", "constant_value"); auto sopac = singularity::photons::GrayS(kappa_s, apm); mscattering = singularity::photons::MeanNonCGSUnitsS( @@ -185,20 +192,28 @@ std::shared_ptr Initialize(ParameterInput *pin) { PARTHENON_FAIL("Only none or constant scattering models supported!"); } } else if (frequency_type == FrequencyType::multigroup) { - if (scattering_model == "none") { + if (sct_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"); + } else if (sct_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 if (sct_model == "powerlaw") { + Real kappa0 = pin->GetReal("mcblock/scattering", "kappa0"); + Real rho_exp = pin->GetReal("mcblock/scattering", "rho_exp"); + Real temp_exp = pin->GetReal("mcblock/scattering", "temp_exp"); + Real nu_exp = pin->GetReal("mcblock/scattering", "nu_exp"); + Real nu_ref = pin->GetReal("mcblock/scattering", "nu_ref"); + scattering = singularity::photons::NonCGSUnitsS( + singularity::photons::PowerLawS(kappa0, rho_exp, temp_exp, nu_exp, nu_ref), + 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); diff --git a/src/mcblock/opacity.hpp b/src/mcblock/opacity.hpp index 4396d09..4ddbb40 100644 --- a/src/mcblock/opacity.hpp +++ b/src/mcblock/opacity.hpp @@ -24,6 +24,7 @@ namespace mcblock { // Reduced absorption variant just for jaybenne using Opacity = singularity::photons::impl::Variant< singularity::photons::NonCGSUnits, + singularity::photons::NonCGSUnits, singularity::photons::NonCGSUnits>; using MeanOpacity = @@ -32,6 +33,7 @@ using MeanOpacity = // Reduced scattering variant just for jaybenne using Scattering = singularity::photons::impl::S_Variant< singularity::photons::NonCGSUnitsS, + singularity::photons::NonCGSUnitsS, singularity::photons::NonCGSUnitsS>; using MeanScattering = From 40ad19c98f1571c3179ccfbd60ca8f9338c51566 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Thu, 18 Jun 2026 15:02:26 -0600 Subject: [PATCH 07/15] + Add extended power-law input (references and offsets) to mcblock. + Build and store vector of log-space group frequency midpoints. + Demo use of frequency group/bin points in sourcing.cpp. --- src/jaybenne/jaybenne.cpp | 14 +++++++++ src/jaybenne/sourcing.cpp | 64 +++++++++++++++++++-------------------- src/mcblock/mcblock.cpp | 42 +++++++++++++++++-------- 3 files changed, 75 insertions(+), 45 deletions(-) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index dc99f2e..ccbcfff 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -380,6 +380,20 @@ std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacit int n_nubins = pin->GetInteger(block_name, "n_nubins"); pkg->AddParam<>("n_nubins", n_nubins); + // Construct and store frequency grid + // NOTE: these are group interior points, not edges + std::vector nu_grid(n_nubins, 0.0); + numin *= time; + numax *= time; + // 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; ++n) { + nu_grid[n] = numin * std::exp((n + 0.5) * dlnu); + } + // store the grid in the parameter input + pkg->AddParam<>("dlnu", dlnu); + pkg->AddParam<>("nu_grid", nu_grid); + pkg->AddParam<>("frequency_type", FrequencyType::multigroup); // Emission CDF diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 765492c..c3084b2 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -44,8 +44,9 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { const Real h = jb_pkg->template Param("planck_constant"); auto &eos = jb_pkg->template Param("eos_d"); int n_nubins = JaybenneNull(); - Real numin = JaybenneNull(); - Real numax = JaybenneNull(); + Real dlnu = JaybenneNull(); + std::vector nu_grid = JaybenneNull>(); + ParArray1D nu_bins; MeanOpacity mopacity; Opacity opacity; if constexpr (FT == FrequencyType::gray) { @@ -53,8 +54,14 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { } 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"); + dlnu = jb_pkg->template Param("dlnu"); + nu_grid = jb_pkg->template Param>("nu_grid"); + nu_bins = ParArray1D("nu_bins", n_nubins); + auto nu_bins_h = nu_bins.GetHostMirror(); + for (int n = 0; n < n_nubins; ++n) { + nu_bins_h(n) = nu_grid[n]; + } + nu_bins.DeepCopy(nu_bins_h); } auto &do_emission = jb_pkg->template Param("do_emission"); auto &source_strategy = jb_pkg->template Param("source_strategy"); @@ -146,14 +153,14 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { const Real &sie = vmesh(b, fjh::sie(), k, j, i); const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie); [[maybe_unused]] const auto gmoded = gmode; - [[maybe_unused]] const Real &sbd = sb; - [[maybe_unused]] const Real &vvd = vv; - [[maybe_unused]] const Real &dtd = dt; + [[maybe_unused]] const auto &sbd = sb; + [[maybe_unused]] const auto &vvd = vv; + [[maybe_unused]] const auto &dtd = dt; [[maybe_unused]] auto mopac = mopacity; [[maybe_unused]] auto opac = opacity; - [[maybe_unused]] const auto numind = numin; - [[maybe_unused]] const auto numaxd = numax; [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto dlnud = dlnu; + [[maybe_unused]] const auto &nu_binsd = nu_bins; Real erad = JaybenneNull(); if constexpr (ST == SourceType::thermal) { erad = (4.0 * sbd / vvd) * std::pow(temp, 4.0) * dv; @@ -163,26 +170,19 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { emis = mopac.Emissivity(rho, temp, gmode); } 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) - Real nu = numind * std::exp(0.5 * dlnu); - Real dnu = nu * dlnu; + // calculate bin width (assuming log bin width) + Real dnu = dlnud * nu_binsd(0); vmesh(b, fj::emission_cdf(0), k, j, i) = - opac.EmissivityPerNu(rho, temp, numind * std::exp(0.5 * dlnu)) * - dnu; + opac.EmissivityPerNu(rho, temp, nu_binsd(0)) * 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) - nu = numind * std::exp((n + 0.5) * dlnu); - dnu = dlnu * nu; + dnu = dlnud * nu_binsd(n); vmesh(b, fj::emission_cdf(n), k, j, i) = - opac.EmissivityPerNu(rho, temp, nu) * dnu + + opac.EmissivityPerNu(rho, temp, nu_binsd(n)) * dnu + vmesh(b, fj::emission_cdf(n - 1), k, j, i); } emis = 0.0; for (int n = 0; n < n_nubinsd; n++) { - const Real dnu = dlnu * numind * std::exp((n + 0.5) * dlnu); + const Real dnu = dlnud * nu_binsd(n); // Get total emissivity emis += vmesh(b, fj::emission_cdf(n), k, j, i); // Normalize emission CDF @@ -266,13 +266,13 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { const Real y_min = coords.template Xc(jb.s) - 0.5 * dx_j; const Real z_min = coords.template Xc(kb.s) - 0.5 * dx_k; const int cell_idx_1d = (k - kb.s) * (nx1 * nx2) + (j - jb.s) * nx1 + (i - ib.s); - [[maybe_unused]] const Real &sbd = sb; - [[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; + [[maybe_unused]] const auto &sbd = sb; + [[maybe_unused]] const auto &dtd = dt; + [[maybe_unused]] const auto &t_startd = t_start; + [[maybe_unused]] const auto hd = h; + [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto dlnud = dlnu; + [[maybe_unused]] const auto &nu_binsd = nu_bins; // Starting index and length of particles in this cell const int &pstart_idx = prefix_sum(b, cell_idx_1d); @@ -313,7 +313,7 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { if constexpr (FT == FrequencyType::gray) { ppack_r(b, ph::energy(), n) = sample_Planck_energy(rng_gen, sbd, temp); } else if constexpr (FT == FrequencyType::multigroup) { - // Sample energy from CDF + // Sample energy (particle frequency) from CDF const Real rand = rng_gen.drand(); int n; for (n = 0; n < n_nubinsd; n++) { @@ -321,9 +321,7 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { break; } } - const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd; - const Real nu = numind * std::exp((n + 0.5) * dlnu); - ppack_r(b, ph::energy(), n) = hd * nu; + ppack_r(b, ph::energy(), n) = hd * nu_binsd(n); } if constexpr (ST == SourceType::emission) { diff --git a/src/mcblock/mcblock.cpp b/src/mcblock/mcblock.cpp index 2473760..6a64704 100644 --- a/src/mcblock/mcblock.cpp +++ b/src/mcblock/mcblock.cpp @@ -145,13 +145,24 @@ std::shared_ptr Initialize(ParameterInput *pin) { singularity::photons::Gray(kappa), time_scale, mass_scale, length_scale, temperature_scale); } else if (abs_model == "powerlaw") { - Real kappa0 = pin->GetReal("mcblock/absorption", "kappa0"); - Real rho_exp = pin->GetReal("mcblock/absorption", "rho_exp"); - Real temp_exp = pin->GetReal("mcblock/absorption", "temp_exp"); - Real nu_exp = pin->GetReal("mcblock/absorption", "nu_exp"); - Real nu_ref = pin->GetReal("mcblock/absorption", "nu_ref"); + // NOTE: reference values (ref) and offsets (off) must always be in cgs units + const Real kappa0 = pin->GetReal("mcblock/absorption", "kappa0"); + const Real rho_exp = pin->GetReal("mcblock/absorption", "rho_exp"); + const Real temp_exp = pin->GetReal("mcblock/absorption", "temp_exp"); + const Real nu_exp = pin->GetReal("mcblock/absorption", "nu_exp"); + const Real nu_ref = pin->GetReal("mcblock/absorption", "nu_ref"); + const Real nu_off = pin->GetReal("mcblock/absorption", "nu_off"); + const Real rho_ref = pin->GetReal("mcblock/absorption", "rho_ref"); + const Real rho_off = pin->GetReal("mcblock/absorption", "rho_off"); + const Real temp_ref = pin->GetReal("mcblock/absorption", "temp_ref"); + const Real temp_off = pin->GetReal("mcblock/absorption", "temp_off"); + const bool do_stim_emit = + pin->GetOrAddBoolean("mcblock/absorption", "do_stim_emit", false); + // const bool do_stim_emit = pin->Get(); opacity = singularity::photons::NonCGSUnits( - singularity::photons::PowerLaw(kappa0, rho_exp, temp_exp, nu_exp, nu_ref), + singularity::photons::PowerLaw(kappa0, rho_exp, temp_exp, nu_exp, nu_ref, + nu_off, rho_ref, rho_off, temp_ref, temp_off, + do_stim_emit), time_scale, mass_scale, length_scale, temperature_scale); } else if (abs_model == "ep_bremss") { opacity = singularity::photons::NonCGSUnits( @@ -202,13 +213,20 @@ std::shared_ptr Initialize(ParameterInput *pin) { singularity::photons::GrayS(kappa_s, apm), time_scale, mass_scale, length_scale, temperature_scale); } else if (sct_model == "powerlaw") { - Real kappa0 = pin->GetReal("mcblock/scattering", "kappa0"); - Real rho_exp = pin->GetReal("mcblock/scattering", "rho_exp"); - Real temp_exp = pin->GetReal("mcblock/scattering", "temp_exp"); - Real nu_exp = pin->GetReal("mcblock/scattering", "nu_exp"); - Real nu_ref = pin->GetReal("mcblock/scattering", "nu_ref"); + // NOTE: reference values (ref) and offsets (off) must always be in cgs units + const Real kappa0 = pin->GetReal("mcblock/scattering", "kappa0"); + const Real rho_exp = pin->GetReal("mcblock/scattering", "rho_exp"); + const Real temp_exp = pin->GetReal("mcblock/scattering", "temp_exp"); + const Real nu_exp = pin->GetReal("mcblock/scattering", "nu_exp"); + const Real nu_ref = pin->GetReal("mcblock/scattering", "nu_ref"); + const Real nu_off = pin->GetReal("mcblock/scattering", "nu_off"); + const Real rho_ref = pin->GetReal("mcblock/scattering", "rho_ref"); + const Real rho_off = pin->GetReal("mcblock/scattering", "rho_off"); + const Real temp_ref = pin->GetReal("mcblock/scattering", "temp_ref"); + const Real temp_off = pin->GetReal("mcblock/scattering", "temp_off"); scattering = singularity::photons::NonCGSUnitsS( - singularity::photons::PowerLawS(kappa0, rho_exp, temp_exp, nu_exp, nu_ref), + singularity::photons::PowerLawS(kappa0, rho_exp, temp_exp, nu_exp, nu_ref, + nu_off, rho_ref, rho_off, temp_ref, temp_off), time_scale, mass_scale, length_scale, temperature_scale); } else { PARTHENON_FAIL("Only none or constant scattering models supported!"); From 09353a8741169553e07525c765d0ff6eb0254e2a Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Thu, 18 Jun 2026 16:12:55 -0600 Subject: [PATCH 08/15] Incorporate nu_bins into scattering, transport, and ddmc_mg_utils. --- src/jaybenne/ddmc_mg_utils.hpp | 32 ++++++++++++++---------------- src/jaybenne/jaybenne.cpp | 30 ++++++++++++++++++---------- src/jaybenne/sample_ddmc_bface.cpp | 1 - src/jaybenne/scattering.hpp | 7 ++----- src/jaybenne/transport.cpp | 18 +++++++++++++++-- src/jaybenne/transport_ddmc.cpp | 24 +++++++++++++++------- 6 files changed, 70 insertions(+), 42 deletions(-) diff --git a/src/jaybenne/ddmc_mg_utils.hpp b/src/jaybenne/ddmc_mg_utils.hpp index 1007bae..2153dd6 100644 --- a/src/jaybenne/ddmc_mg_utils.hpp +++ b/src/jaybenne/ddmc_mg_utils.hpp @@ -43,7 +43,7 @@ struct ddmc_mg_cell_args { template KOKKOS_FORCEINLINE_FUNCTION std::pair calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args &dmg, - const bool &use_lo) { + const ParArray1D &nu_bins, const bool &use_lo) { // define extrapolation distance (Habetler & Matkowski 1975) constexpr Real lam_ext = 0.7104; @@ -57,8 +57,7 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args for (int n = 0; n < dmg.n_nubins; ++n) { // this is the mid-point of group n in log-space: - // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real ee = dmg.hd * dmg.numind * std::exp((n + 0.5) * dmg.dlnu); + const Real ee = dmg.hd * nu_bins(n); const Real dee = ee * dmg.dlnu; // evaluate face opacities at nu @@ -106,8 +105,9 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args template KOKKOS_FORCEINLINE_FUNCTION Real calc_ddmc_mg_leakprob(const OP &abs, const SC &sct, const ddmc_mg_leak_args &dmg, + const ParArray1D &nu_bins, const bool &use_lo) { - const auto leak_sum_pair = calc_ddmc_mg_leak_numdenom(abs, sct, dmg, use_lo); + const auto leak_sum_pair = calc_ddmc_mg_leak_numdenom(abs, sct, dmg, nu_bins, use_lo); if (!(leak_sum_pair.second > 0.0)) return 0.0; // normalize sum (so that leakage probability is an average) return leak_sum_pair.first / leak_sum_pair.second; @@ -117,11 +117,12 @@ KOKKOS_FORCEINLINE_FUNCTION Real calc_ddmc_mg_leakprob(const OP &abs, const SC & template KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &sct, const ddmc_mg_leak_args &dmg, + const ParArray1D &nu_bins, const bool &use_lo, RngGen &rng_gen) { // first get CDF totals - const auto leak_sum_pair = calc_ddmc_mg_leak_numdenom(abs, sct, dmg, use_lo); + const auto leak_sum_pair = calc_ddmc_mg_leak_numdenom(abs, sct, dmg, nu_bins, use_lo); const Real &leak_tot_sum = leak_sum_pair.first; const Real &planck_tot_sum = leak_sum_pair.second; @@ -139,8 +140,7 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &s for (int n = 0; n < dmg.n_nubins; ++n) { // this is the mid-point of group n in log-space: - // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real ee = dmg.hd * dmg.numind * std::exp((n + 0.5) * dmg.dlnu); + const Real ee = dmg.hd * nu_bins(n); const Real dee = ee * dmg.dlnu; // evaluate face opacities at nu @@ -191,7 +191,8 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &s // calculate: absorption, outscatter probability template KOKKOS_FORCEINLINE_FUNCTION std::pair -calc_ddmc_mg_probs(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc) { +calc_ddmc_mg_probs(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc, + const ParArray1D &nu_bins) { // initialize unnnormalized Planck integral and probability Real planck_sum = 0.0; @@ -202,8 +203,7 @@ calc_ddmc_mg_probs(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc) for (int n = 0; n < dmgc.n_nubins; ++n) { // this is the mid-point of group n in log-space: - // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real ee = dmgc.hd * dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); + const Real ee = dmgc.hd * nu_bins(n); const Real dee = ee * dmgc.dlnu; // evaluate face opacities at nu @@ -240,9 +240,9 @@ calc_ddmc_mg_probs(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc) // sample out-scatter IMC group // NOTE(MGDDMC): this routine is assuming Kirchhoff's Law for emissivity (LTE) template -KOKKOS_FORCEINLINE_FUNCTION Real sample_ddmc2imc_outscatter(const OP &abs, const SC &sct, - const ddmc_mg_cell_args &dmgc, - RngGen &rng_gen) { +KOKKOS_FORCEINLINE_FUNCTION Real +sample_ddmc2imc_outscatter(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc, + const ParArray1D &nu_bins, RngGen &rng_gen) { Real scat_out_tot_sum = 0.0; @@ -250,8 +250,7 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_ddmc2imc_outscatter(const OP &abs, const for (int n = 0; n < dmgc.n_nubins; ++n) { // this is the mid-point of group n in log-space: - // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real ee = dmgc.hd * dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); + const Real ee = dmgc.hd * nu_bins(n); const Real dee = ee * dmgc.dlnu; // evaluate face opacities at nu @@ -276,8 +275,7 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_ddmc2imc_outscatter(const OP &abs, const for (int n = 0; n < dmgc.n_nubins; ++n) { // this is the mid-point of group n in log-space: - // nu=exp(log(numin)+(n+0.5)*dlnu) - const Real ee = dmgc.hd * dmgc.numind * std::exp((n + 0.5) * dmgc.dlnu); + const Real ee = dmgc.hd * nu_bins(n); const Real dee = ee * dmgc.dlnu; // evaluate face opacities at nu diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index ccbcfff..50d94da 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -446,10 +446,11 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { MeanScattering mscattering; int n_nubins = -1; Real numin = -1.0; - Real numax = -1.0; Real dlnu = -1.0; Real h = -1.0; Real sb = -1.0; + std::vector nu_grid = JaybenneNull>(); + ParArray1D nu_bins; // get opacity average indicators const auto &use_planck = jbn->template Param("use_planck"); @@ -467,11 +468,17 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { scattering = jbn->template Param("scattering_d"); n_nubins = jbn->template Param("n_nubins"); numin = jbn->template Param("numin"); - numax = jbn->template Param("numax"); h = jbn->template Param("planck_constant"); sb = jbn->template Param("boltzmann"); - // initialize (assumed) log spacing - dlnu = (std::log(numax) - std::log(numin)) / n_nubins; + // initialize (assumed) log-spaced frequency grid + dlnu = jbn->template Param("dlnu"); + nu_grid = jbn->template Param>("nu_grid"); + nu_bins = ParArray1D("nu_bins", n_nubins); + auto nu_bins_h = nu_bins.GetHostMirror(); + for (int n = 0; n < n_nubins; ++n) { + nu_bins_h(n) = nu_grid[n]; + } + nu_bins.DeepCopy(nu_bins_h); } const auto &ib = md->GetBoundsI(IndexDomain::interior); @@ -591,6 +598,7 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] const auto numind = numin; [[maybe_unused]] const auto dlnud = dlnu; [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto &nu_binsd = nu_bins; [[maybe_unused]] const auto hd = h; [[maybe_unused]] const auto sbd = sb; [[maybe_unused]] const auto tau_ddmcd = tau_ddmc; @@ -641,11 +649,11 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { // integrate lo-x leakage probability vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), k, j, i) = - calc_ddmc_mg_leakprob(opac, scatter, dmg, use_lo); + calc_ddmc_mg_leakprob(opac, scatter, dmg, nu_binsd, use_lo); // integrate hi-x leakage probability vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), k, j, i) = - calc_ddmc_mg_leakprob(opac, scatter, dmg, use_hi); + calc_ddmc_mg_leakprob(opac, scatter, dmg, nu_binsd, use_hi); } }); @@ -708,6 +716,7 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] const auto numind = numin; [[maybe_unused]] const auto dlnud = dlnu; [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto &nu_binsd = nu_bins; [[maybe_unused]] const auto hd = h; [[maybe_unused]] const auto sbd = sb; [[maybe_unused]] const auto tau_ddmcd = tau_ddmc; @@ -759,11 +768,11 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { // integrate lo-y leakage probability vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), k, j, i) = - calc_ddmc_mg_leakprob(opac, scatter, dmg, use_lo); + calc_ddmc_mg_leakprob(opac, scatter, dmg, nu_binsd, use_lo); // integrate hi-y leakage probability vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), k, j, i) = - calc_ddmc_mg_leakprob(opac, scatter, dmg, use_hi); + calc_ddmc_mg_leakprob(opac, scatter, dmg, nu_binsd, use_hi); } }); } @@ -827,6 +836,7 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] const auto numind = numin; [[maybe_unused]] const auto dlnud = dlnu; [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto &nu_binsd = nu_bins; [[maybe_unused]] const auto hd = h; [[maybe_unused]] const auto sbd = sb; [[maybe_unused]] const auto tau_ddmcd = tau_ddmc; @@ -878,11 +888,11 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { // integrate lo-z leakage probability vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), k, j, i) = - calc_ddmc_mg_leakprob(opac, scatter, dmg, use_lo); + calc_ddmc_mg_leakprob(opac, scatter, dmg, nu_binsd, use_lo); // integrate hi-z leakage probability vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), k, j, i) = - calc_ddmc_mg_leakprob(opac, scatter, dmg, use_hi); + calc_ddmc_mg_leakprob(opac, scatter, dmg, nu_bins, use_hi); } }); } diff --git a/src/jaybenne/sample_ddmc_bface.cpp b/src/jaybenne/sample_ddmc_bface.cpp index 289a040..06b3bea 100644 --- a/src/jaybenne/sample_ddmc_bface.cpp +++ b/src/jaybenne/sample_ddmc_bface.cpp @@ -96,7 +96,6 @@ TaskStatus SampleDDMCBlockFace(MeshData *md) { const Real vv = jb_pkg->template Param("speed_of_light"); // get dimension indicators (for face probabilities) - const bool multi_d = (pm->ndim > 1); const bool three_d = (pm->ndim > 2); // Create SparsePack diff --git a/src/jaybenne/scattering.hpp b/src/jaybenne/scattering.hpp index 7358755..9547847 100644 --- a/src/jaybenne/scattering.hpp +++ b/src/jaybenne/scattering.hpp @@ -27,8 +27,6 @@ struct cell_scat_args { const Real &ss; // scattering opacity (1/length) // frequency data const int &n_nubinsd; // number of frequency groups - const Real &numind; // mininmum frequency - const Real &numaxd; // maximum frequency const Real &hd; // Planck constant }; struct ptcl_scat_args { @@ -57,6 +55,7 @@ void sample_vol_iso_dir(ptcl_scat_args sa) { //! TODO(BRR): template on scattering kernel type? template KOKKOS_FORCEINLINE_FUNCTION void scatter_kernel(const T &vmesh, cell_scat_args csa, + const ParArray1D &nu_bins, ptcl_scat_args psa) { namespace fj = field::jaybenne; @@ -75,11 +74,9 @@ KOKKOS_FORCEINLINE_FUNCTION void scatter_kernel(const T &vmesh, cell_scat_args c break; } } - const Real dlnu = (std::log(csa.numaxd) - std::log(csa.numind)) / csa.n_nubinsd; - const Real nu = csa.numind * std::exp((n + 0.5) * dlnu); // reset particle frequency - psa.ee = csa.hd * nu; + psa.ee = csa.hd * nu_bins(n); } } diff --git a/src/jaybenne/transport.cpp b/src/jaybenne/transport.cpp index c268ee0..6e50cbb 100644 --- a/src/jaybenne/transport.cpp +++ b/src/jaybenne/transport.cpp @@ -43,12 +43,24 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d int n_nubins = JaybenneNull(); Real numin = JaybenneNull(); Real numax = JaybenneNull(); + Real dlnu = JaybenneNull(); + std::vector nu_grid = JaybenneNull>(); + ParArray1D nu_bins; 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"); + // initialize (assumed) log-spaced frequency bins + dlnu = jb_pkg->template Param("dlnu"); + nu_grid = jb_pkg->template Param>("nu_grid"); + nu_bins = ParArray1D("nu_bins", n_nubins); + auto nu_bins_h = nu_bins.GetHostMirror(); + for (int n = 0; n < n_nubins; ++n) { + nu_bins_h(n) = nu_grid[n]; + } + nu_bins.DeepCopy(nu_bins_h); } auto &rng_pool = jb_pkg->template Param("rng_pool"); const Real vv = jb_pkg->template Param("speed_of_light"); @@ -93,6 +105,8 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d [[maybe_unused]] const auto numind = numin; [[maybe_unused]] const auto numaxd = numax; [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto dlnud = dlnu; + [[maybe_unused]] const auto nu_binsd = nu_bins; auto &coords = vmesh.GetCoordinates(b); const Real &dx_i = coords.template Dxc(0, 0, 0); @@ -236,11 +250,11 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d // clang-format off cell_scat_args csa{b, ip, jp, kp, ff, aa, ss, - n_nubinsd, numind, numaxd, hd}; + n_nubinsd, hd}; // clang-format on // invoke frequency-dependent scattering kernel - scatter_kernel(vmesh, csa, psa); + scatter_kernel(vmesh, csa, nu_binsd, psa); } } diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index 423a90c..dea65a1 100644 --- a/src/jaybenne/transport_ddmc.cpp +++ b/src/jaybenne/transport_ddmc.cpp @@ -52,14 +52,23 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R Real numin = JaybenneNull(); Real numax = JaybenneNull(); Real dlnu = JaybenneNull(); + std::vector nu_grid = JaybenneNull>(); + ParArray1D nu_bins; Opacity opacity; Scattering scattering; if constexpr (FT == FrequencyType::multigroup) { n_nubins = jb_pkg->template Param("n_nubins"); numin = jb_pkg->template Param("numin"); numax = jb_pkg->template Param("numax"); - // initialize (assumed) log spacing - dlnu = (std::log(numax) - std::log(numin)) / n_nubins; + // initialize (assumed) log-spaced frequency bins + dlnu = jb_pkg->template Param("dlnu"); + nu_grid = jb_pkg->template Param>("nu_grid"); + nu_bins = ParArray1D("nu_bins", n_nubins); + auto nu_bins_h = nu_bins.GetHostMirror(); + for (int n = 0; n < n_nubins; ++n) { + nu_bins_h(n) = nu_grid[n]; + } + nu_bins.DeepCopy(nu_bins_h); // set opacity objects opacity = jb_pkg->template Param("opacity_d"); scattering = jb_pkg->template Param("scattering_d"); @@ -110,6 +119,7 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R [[maybe_unused]] const auto numaxd = numax; [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto dlnud = dlnu; + [[maybe_unused]] const auto nu_binsd = nu_bins; auto &coords = vmesh.GetCoordinates(b); const Real &dx_i = coords.template Dxc(0, 0, 0); @@ -235,7 +245,7 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R rho, temp}; // clang-format on - const auto aagm = calc_ddmc_mg_probs(opac, scatter, dmgc); + const auto aagm = calc_ddmc_mg_probs(opac, scatter, dmgc, nu_binsd); aa_g = aagm.first; gm_g = aagm.second; } @@ -316,7 +326,7 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // clang-format off // sample particle frequency (energy units) - ee = sample_leakage_group(opac, scatter, dmg, use_lo, rng_gen); + ee = sample_leakage_group(opac, scatter, dmg, nu_binsd, use_lo, rng_gen); } } @@ -410,7 +420,7 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R rho, temp}; // clang-format on - ee = sample_ddmc2imc_outscatter(opac, scatter, dmgc, rng_gen); + ee = sample_ddmc2imc_outscatter(opac, scatter, dmgc, nu_bins, rng_gen); // resample particle direction sample_vol_iso_dir(psa); @@ -421,11 +431,11 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // clang-format off cell_scat_args csa{b, ip, jp, kp, ff, aa, ss, - n_nubinsd, numind, numaxd, hd}; + n_nubinsd, hd}; // clang-format on // invoke frequency-dependent scattering kernel - scatter_kernel(vmesh, csa, psa); + scatter_kernel(vmesh, csa, nu_binsd, psa); } } } From 2c698d7a497d006b3e0cc8e1d3c7269577be9b83 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Thu, 18 Jun 2026 16:28:36 -0600 Subject: [PATCH 09/15] Remove numin and numax from storage in jaybenne parameters. --- src/jaybenne/ddmc_mg_utils.hpp | 2 -- src/jaybenne/jaybenne.cpp | 16 ++++------------ src/jaybenne/transport.cpp | 6 ------ src/jaybenne/transport_ddmc.cpp | 9 --------- 4 files changed, 4 insertions(+), 29 deletions(-) diff --git a/src/jaybenne/ddmc_mg_utils.hpp b/src/jaybenne/ddmc_mg_utils.hpp index 2153dd6..6434f06 100644 --- a/src/jaybenne/ddmc_mg_utils.hpp +++ b/src/jaybenne/ddmc_mg_utils.hpp @@ -11,7 +11,6 @@ namespace jaybenne { // helper struct to encapsulate data needed to integrate DDMC probabilities over groups struct ddmc_mg_leak_args { const int &n_nubins; // number of frequency groups (or bins) - const Real &numind; // minimum frequency const Real &dlnu; // log-spacing of frequency groups const Real &hd; // Planck constant const Real &sbd; // Stefan-Boltzmann constant @@ -29,7 +28,6 @@ struct ddmc_mg_leak_args { // helper struct for in-cell MG DDMC processes struct ddmc_mg_cell_args { const int &n_nubins; // number of frequency groups (or bins) - const Real &numind; // minimum frequency const Real &dlnu; // log-spacing of frequency groups const Real &hd; // Planck constant const Real &sbd; // Stefan-Boltzmann constant diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 50d94da..8c63d95 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -374,17 +374,17 @@ std::shared_ptr Initialize(ParameterInput *pin, Opacity &opacit // Frequency discretization auto time = units.time; Real numin = pin->GetReal(block_name, "numin"); // in Hz - pkg->AddParam<>("numin", numin * time); // in code units Real numax = pin->GetReal(block_name, "numax"); // in Hz - pkg->AddParam<>("numax", numax * time); // in code units int n_nubins = pin->GetInteger(block_name, "n_nubins"); pkg->AddParam<>("n_nubins", n_nubins); + // assume units.time = [s/code time] = [code freq/Hz] + numin *= time; // in code units + numax *= time; // in code units + // Construct and store frequency grid // NOTE: these are group interior points, not edges std::vector nu_grid(n_nubins, 0.0); - numin *= time; - numax *= time; // 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; ++n) { @@ -445,7 +445,6 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { MeanOpacity mopacity; MeanScattering mscattering; int n_nubins = -1; - Real numin = -1.0; Real dlnu = -1.0; Real h = -1.0; Real sb = -1.0; @@ -467,7 +466,6 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { opacity = jbn->template Param("opacity_d"); scattering = jbn->template Param("scattering_d"); n_nubins = jbn->template Param("n_nubins"); - numin = jbn->template Param("numin"); h = jbn->template Param("planck_constant"); sb = jbn->template Param("boltzmann"); // initialize (assumed) log-spaced frequency grid @@ -595,7 +593,6 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] auto mscatter = mscattering; [[maybe_unused]] auto opac = opacity; [[maybe_unused]] auto scatter = scattering; - [[maybe_unused]] const auto numind = numin; [[maybe_unused]] const auto dlnud = dlnu; [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto &nu_binsd = nu_bins; @@ -628,7 +625,6 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { // create DDMC MG leakage data helper argument (defined in ddmc_mg_utils.hpp) // clang-format off const ddmc_mg_leak_args dmg{n_nubinsd, - numind, dlnud, hd, sbd, @@ -713,7 +709,6 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] auto mscatter = mscattering; [[maybe_unused]] auto opac = opacity; [[maybe_unused]] auto scatter = scattering; - [[maybe_unused]] const auto numind = numin; [[maybe_unused]] const auto dlnud = dlnu; [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto &nu_binsd = nu_bins; @@ -747,7 +742,6 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { // ddmc_mg_utils.hpp) // clang-format off const ddmc_mg_leak_args dmg{n_nubins, - numind, dlnud, hd, sbd, @@ -833,7 +827,6 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] auto mscatter = mscattering; [[maybe_unused]] auto opac = opacity; [[maybe_unused]] auto scatter = scattering; - [[maybe_unused]] const auto numind = numin; [[maybe_unused]] const auto dlnud = dlnu; [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto &nu_binsd = nu_bins; @@ -867,7 +860,6 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { // ddmc_mg_utils.hpp) // clang-format off const ddmc_mg_leak_args dmg{n_nubins, - numind, dlnud, hd, sbd, diff --git a/src/jaybenne/transport.cpp b/src/jaybenne/transport.cpp index 6e50cbb..1da5127 100644 --- a/src/jaybenne/transport.cpp +++ b/src/jaybenne/transport.cpp @@ -41,8 +41,6 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d Opacity opacity; Scattering scattering; int n_nubins = JaybenneNull(); - Real numin = JaybenneNull(); - Real numax = JaybenneNull(); Real dlnu = JaybenneNull(); std::vector nu_grid = JaybenneNull>(); ParArray1D nu_bins; @@ -50,8 +48,6 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real 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"); // initialize (assumed) log-spaced frequency bins dlnu = jb_pkg->template Param("dlnu"); nu_grid = jb_pkg->template Param>("nu_grid"); @@ -102,8 +98,6 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d // frequency data, needed for multigroup [[maybe_unused]] const auto hd = h; - [[maybe_unused]] const auto numind = numin; - [[maybe_unused]] const auto numaxd = numax; [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto dlnud = dlnu; [[maybe_unused]] const auto nu_binsd = nu_bins; diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index dea65a1..1bb2f99 100644 --- a/src/jaybenne/transport_ddmc.cpp +++ b/src/jaybenne/transport_ddmc.cpp @@ -49,8 +49,6 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R const Real h = jb_pkg->template Param("planck_constant"); const Real sb = jb_pkg->template Param("boltzmann"); int n_nubins = JaybenneNull(); - Real numin = JaybenneNull(); - Real numax = JaybenneNull(); Real dlnu = JaybenneNull(); std::vector nu_grid = JaybenneNull>(); ParArray1D nu_bins; @@ -58,8 +56,6 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R Scattering scattering; if constexpr (FT == FrequencyType::multigroup) { n_nubins = jb_pkg->template Param("n_nubins"); - numin = jb_pkg->template Param("numin"); - numax = jb_pkg->template Param("numax"); // initialize (assumed) log-spaced frequency bins dlnu = jb_pkg->template Param("dlnu"); nu_grid = jb_pkg->template Param>("nu_grid"); @@ -115,8 +111,6 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // frequency data, needed for multigroup [[maybe_unused]] const auto hd = h; [[maybe_unused]] const auto sbd = sb; - [[maybe_unused]] const auto numind = numin; - [[maybe_unused]] const auto numaxd = numax; [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto dlnud = dlnu; [[maybe_unused]] const auto nu_binsd = nu_bins; @@ -236,7 +230,6 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R if constexpr (FT == FrequencyType::multigroup) { // clang-format off ddmc_mg_cell_args dmgc{n_nubinsd, - numind, dlnud, hd, sbd, @@ -310,7 +303,6 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // clang-format off const ddmc_mg_leak_args dmg{n_nubinsd, - numind, dlnud, hd, sbd, @@ -411,7 +403,6 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R if (is_ddmc_step) { // clang-format off ddmc_mg_cell_args dmgc{n_nubinsd, - numind, dlnud, hd, sbd, From f795b1ade843d7e21d3e78c871325ecd6db1379d Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Thu, 18 Jun 2026 19:08:42 -0600 Subject: [PATCH 10/15] + Fix group-particle index shadowing error in sourcing.cpp. + Make powerlaw offsets optional. --- src/jaybenne/planck.hpp | 4 ++-- src/jaybenne/sourcing.cpp | 14 ++++++++------ src/mcblock/mcblock.cpp | 12 ++++++------ 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/jaybenne/planck.hpp b/src/jaybenne/planck.hpp index 903f663..d52d4aa 100644 --- a/src/jaybenne/planck.hpp +++ b/src/jaybenne/planck.hpp @@ -42,7 +42,7 @@ Real midpoint_Planck(const Real &temp, const Real &nu, const Real &dnu) { //! rng_gen: RNG pool //! T: distribution temperature KOKKOS_FORCEINLINE_FUNCTION -Real sample_Planck_energy(RngGen &rng_gen, const Real &sb, const Real &temp) { +Real sample_Planck_energy(RngGen &rng_gen, const Real &kbolt, const Real &temp) { // Sampling method from Everett & Cashwell 1972 const Real xi0 = rng_gen.drand(); const Real rhs = xi0 * std::pow(M_PI, 4.0) / 90.0; @@ -64,7 +64,7 @@ Real sample_Planck_energy(RngGen &rng_gen, const Real &sb, const Real &temp) { const Real xi2 = rng_gen.drand(); const Real xi3 = rng_gen.drand(); const Real xi4 = rng_gen.drand(); - return -(1.0 / ll) * std::log(xi1 * xi2 * xi3 * xi4) * sb * temp; + return -(1.0 / ll) * std::log(xi1 * xi2 * xi3 * xi4) * kbolt * temp; } } // namespace jaybenne diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index c3084b2..5c33ae8 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -83,6 +83,7 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { const Real &emit_temp_th = jb_pkg->template Param("emit_temp_threshold"); const Real &vv = jb_pkg->template Param("speed_of_light"); const Real &sb = jb_pkg->template Param("stefan_boltzmann"); + const Real &kbolt = jb_pkg->template Param("boltzmann"); // Create pack static auto desc = @@ -164,6 +165,7 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { Real erad = JaybenneNull(); if constexpr (ST == SourceType::thermal) { erad = (4.0 * sbd / vvd) * std::pow(temp, 4.0) * dv; + // TODO: add multigroup thermal mode } else if constexpr (ST == SourceType::emission) { Real emis = JaybenneNull(); if constexpr (FT == FrequencyType::gray) { @@ -266,7 +268,7 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { const Real y_min = coords.template Xc(jb.s) - 0.5 * dx_j; const Real z_min = coords.template Xc(kb.s) - 0.5 * dx_k; const int cell_idx_1d = (k - kb.s) * (nx1 * nx2) + (j - jb.s) * nx1 + (i - ib.s); - [[maybe_unused]] const auto &sbd = sb; + [[maybe_unused]] const auto &kboltd = kbolt; [[maybe_unused]] const auto &dtd = dt; [[maybe_unused]] const auto &t_startd = t_start; [[maybe_unused]] const auto hd = h; @@ -311,17 +313,17 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { const Real &sie = vmesh(b, fjh::sie(), k, j, i); const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie); if constexpr (FT == FrequencyType::gray) { - ppack_r(b, ph::energy(), n) = sample_Planck_energy(rng_gen, sbd, temp); + ppack_r(b, ph::energy(), n) = sample_Planck_energy(rng_gen, kboltd, temp); } else if constexpr (FT == FrequencyType::multigroup) { // Sample energy (particle frequency) 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) { + int g; + for (g = 0; g < n_nubinsd; ++g) { + if (vmesh(b, fj::emission_cdf(g), k, j, i) >= rand) { break; } } - ppack_r(b, ph::energy(), n) = hd * nu_binsd(n); + ppack_r(b, ph::energy(), n) = hd * nu_binsd(g); } if constexpr (ST == SourceType::emission) { diff --git a/src/mcblock/mcblock.cpp b/src/mcblock/mcblock.cpp index 6a64704..8a2e340 100644 --- a/src/mcblock/mcblock.cpp +++ b/src/mcblock/mcblock.cpp @@ -151,11 +151,11 @@ std::shared_ptr Initialize(ParameterInput *pin) { const Real temp_exp = pin->GetReal("mcblock/absorption", "temp_exp"); const Real nu_exp = pin->GetReal("mcblock/absorption", "nu_exp"); const Real nu_ref = pin->GetReal("mcblock/absorption", "nu_ref"); - const Real nu_off = pin->GetReal("mcblock/absorption", "nu_off"); + const Real nu_off = pin->GetOrAddReal("mcblock/absorption", "nu_off", 0.0); const Real rho_ref = pin->GetReal("mcblock/absorption", "rho_ref"); - const Real rho_off = pin->GetReal("mcblock/absorption", "rho_off"); + const Real rho_off = pin->GetOrAddReal("mcblock/absorption", "rho_off", 0.0); const Real temp_ref = pin->GetReal("mcblock/absorption", "temp_ref"); - const Real temp_off = pin->GetReal("mcblock/absorption", "temp_off"); + const Real temp_off = pin->GetOrAddReal("mcblock/absorption", "temp_off", 0.0); const bool do_stim_emit = pin->GetOrAddBoolean("mcblock/absorption", "do_stim_emit", false); // const bool do_stim_emit = pin->Get(); @@ -219,11 +219,11 @@ std::shared_ptr Initialize(ParameterInput *pin) { const Real temp_exp = pin->GetReal("mcblock/scattering", "temp_exp"); const Real nu_exp = pin->GetReal("mcblock/scattering", "nu_exp"); const Real nu_ref = pin->GetReal("mcblock/scattering", "nu_ref"); - const Real nu_off = pin->GetReal("mcblock/scattering", "nu_off"); + const Real nu_off = pin->GetOrAddReal("mcblock/scattering", "nu_off", 0.0); const Real rho_ref = pin->GetReal("mcblock/scattering", "rho_ref"); - const Real rho_off = pin->GetReal("mcblock/scattering", "rho_off"); + const Real rho_off = pin->GetOrAddReal("mcblock/scattering", "rho_off", 0.0); const Real temp_ref = pin->GetReal("mcblock/scattering", "temp_ref"); - const Real temp_off = pin->GetReal("mcblock/scattering", "temp_off"); + const Real temp_off = pin->GetOrAddReal("mcblock/scattering", "temp_off", 0.0); scattering = singularity::photons::NonCGSUnitsS( singularity::photons::PowerLawS(kappa0, rho_exp, temp_exp, nu_exp, nu_ref, nu_off, rho_ref, rho_off, temp_ref, temp_off), From 41018207062b635e36047954a909de3709efa5b9 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Fri, 19 Jun 2026 13:38:49 -0600 Subject: [PATCH 11/15] Set emission_cdf as discrete Planck function for thermal sourcing. --- src/jaybenne/sourcing.cpp | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 5c33ae8..6e3f831 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -155,6 +155,8 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie); [[maybe_unused]] const auto gmoded = gmode; [[maybe_unused]] const auto &sbd = sb; + [[maybe_unused]] const auto &kboltd = kbolt; + [[maybe_unused]] const auto hd = h; [[maybe_unused]] const auto &vvd = vv; [[maybe_unused]] const auto &dtd = dt; [[maybe_unused]] auto mopac = mopacity; @@ -165,7 +167,26 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { Real erad = JaybenneNull(); if constexpr (ST == SourceType::thermal) { erad = (4.0 * sbd / vvd) * std::pow(temp, 4.0) * dv; - // TODO: add multigroup thermal mode + // leverage emission_cdf for initial Planck sampling + if constexpr (FT == FrequencyType::multigroup) { + // calculate bin width (assuming log bin width) + Real ee = hd * nu_binsd(0); + Real dee = dlnud * ee; + vmesh(b, fj::emission_cdf(0), k, j, i) = + jaybenne::midpoint_Planck(kboltd * temp, ee, dee); + for (int n = 1; n < n_nubinsd; n++) { + ee = hd * nu_binsd(n); + dee = dlnud * ee; + vmesh(b, fj::emission_cdf(n), k, j, i) = + jaybenne::midpoint_Planck(kboltd * temp, ee, dee) + + vmesh(b, fj::emission_cdf(n - 1), k, j, i); + } + for (int n = 0; n < n_nubinsd; n++) { + // Normalize emission CDF + vmesh(b, fj::emission_cdf(n), k, j, i) /= + vmesh(b, fj::emission_cdf(n_nubinsd - 1), k, j, i); + } + } } else if constexpr (ST == SourceType::emission) { Real emis = JaybenneNull(); if constexpr (FT == FrequencyType::gray) { @@ -184,7 +205,6 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { } emis = 0.0; for (int n = 0; n < n_nubinsd; n++) { - const Real dnu = dlnud * nu_binsd(n); // Get total emissivity emis += vmesh(b, fj::emission_cdf(n), k, j, i); // Normalize emission CDF From 1e32d4af4efbf79202d51506e85a143ad1676bbf Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Sat, 20 Jun 2026 19:58:42 -0600 Subject: [PATCH 12/15] Fix multifrequency units bug and add DDMC elastic scattering mode. + Use frequency in 1/time in singularity-opac API functions. + Check heuristic if elastic scattering is dominant (is_elastic). + Avoid frequency redistribution during cell leakage if is_elastic. Note: the elastic scattering-dominant mode is a provisional solution, which needs more work for better consistency with the original inelastic mode. --- src/jaybenne/ddmc_mg_utils.hpp | 72 ++++++++++++++++----------------- src/jaybenne/jaybenne.cpp | 10 ++++- src/jaybenne/transport.cpp | 6 ++- src/jaybenne/transport_ddmc.cpp | 52 +++++++++++++++++------- 4 files changed, 85 insertions(+), 55 deletions(-) diff --git a/src/jaybenne/ddmc_mg_utils.hpp b/src/jaybenne/ddmc_mg_utils.hpp index 6434f06..993b013 100644 --- a/src/jaybenne/ddmc_mg_utils.hpp +++ b/src/jaybenne/ddmc_mg_utils.hpp @@ -54,15 +54,11 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args // integrate for (int n = 0; n < dmg.n_nubins; ++n) { - // this is the mid-point of group n in log-space: - const Real ee = dmg.hd * nu_bins(n); - const Real dee = ee * dmg.dlnu; - - // evaluate face opacities at nu - const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, ee); - const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, ee); - const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, ee); - const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, ee); + // evaluate face opacities at nu_bins(n) + const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, nu_bins(n)); + const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, nu_bins(n)); + const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, nu_bins(n)); + const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, nu_bins(n)); // calculate optical thicknesses from lower and upper cell const Real tau_lmin = dmg.dx_lmin * (ss_l + aa_l); @@ -85,6 +81,10 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args // evaluate a face temperature (note this is typically a 4-norm average) const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); + // convert to energy units for Planck integral + const Real ee = dmg.hd * nu_bins(n); + const Real dee = ee * dmg.dlnu; + // get an unnormalized, non-dimensional Planck integral over group const Real bg = midpoint_Planck(dmg.sbd * temp_f, ee, dee); @@ -137,15 +137,11 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &s // integrate for (int n = 0; n < dmg.n_nubins; ++n) { - // this is the mid-point of group n in log-space: - const Real ee = dmg.hd * nu_bins(n); - const Real dee = ee * dmg.dlnu; - - // evaluate face opacities at nu - const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, ee); - const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, ee); - const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, ee); - const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, ee); + // evaluate face opacities at nu_bins(n) + const Real ss_l = sct.TotalScatteringCoefficient(dmg.rho_l, dmg.temp_l, nu_bins(n)); + const Real aa_l = abs.AbsorptionCoefficient(dmg.rho_l, dmg.temp_l, nu_bins(n)); + const Real ss_u = sct.TotalScatteringCoefficient(dmg.rho_u, dmg.temp_u, nu_bins(n)); + const Real aa_u = abs.AbsorptionCoefficient(dmg.rho_u, dmg.temp_u, nu_bins(n)); // calculate optical thicknesses from lower and upper cell const Real tau_lmin = dmg.dx_lmin * (ss_l + aa_l); @@ -166,6 +162,10 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &s // evaluate a face temperature (note this is typically a 4-norm average) const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); + // convert to energy units for Planck integral + const Real ee = dmg.hd * nu_bins(n); + const Real dee = ee * dmg.dlnu; + // get an unnormalized, non-dimensional Planck integral over group const Real bg = midpoint_Planck(dmg.sbd * temp_f, ee, dee); @@ -200,14 +200,14 @@ calc_ddmc_mg_probs(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc, // integrate for (int n = 0; n < dmgc.n_nubins; ++n) { - // this is the mid-point of group n in log-space: + // evaluate face opacities at nu_bins(n) + const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, nu_bins(n)); + const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu_bins(n)); + + // convert to energy units for Planck integral const Real ee = dmgc.hd * nu_bins(n); const Real dee = ee * dmgc.dlnu; - // evaluate face opacities at nu - const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, ee); - const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, ee); - // get an unnormalized, non-dimensional Planck integral over group const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, ee, dee); @@ -247,17 +247,15 @@ sample_ddmc2imc_outscatter(const OP &abs, const SC &sct, const ddmc_mg_cell_args // integrate total cdf value for (int n = 0; n < dmgc.n_nubins; ++n) { - // this is the mid-point of group n in log-space: - const Real ee = dmgc.hd * nu_bins(n); - const Real dee = ee * dmgc.dlnu; - - // evaluate face opacities at nu - const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, ee); - const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, ee); + // evaluate face opacities at nu_bins(n) + const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, nu_bins(n)); + const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu_bins(n)); // check group exclusion if (!(dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc)) { // get an unnormalized, non-dimensional Planck integral over group + const Real ee = dmgc.hd * nu_bins(n); + const Real dee = ee * dmgc.dlnu; const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, ee, dee); scat_out_tot_sum += bg * aa; } @@ -272,17 +270,17 @@ sample_ddmc2imc_outscatter(const OP &abs, const SC &sct, const ddmc_mg_cell_args // find bin for sample for (int n = 0; n < dmgc.n_nubins; ++n) { - // this is the mid-point of group n in log-space: - const Real ee = dmgc.hd * nu_bins(n); - const Real dee = ee * dmgc.dlnu; - - // evaluate face opacities at nu - const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, ee); - const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, ee); + // evaluate face opacities at nu_bins(n) + const Real ss = sct.TotalScatteringCoefficient(dmgc.rho, dmgc.temp, nu_bins(n)); + const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu_bins(n)); // check group exclusion if (!(dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc)) { + // convert to energy units for Planck integral + const Real ee = dmgc.hd * nu_bins(n); + const Real dee = ee * dmgc.dlnu; + // get an unnormalized, non-dimensional Planck integral over group const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, ee, dee); diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 8c63d95..9433d67 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -501,10 +501,18 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { [[maybe_unused]] const auto gmoded = gmode; [[maybe_unused]] auto mopac = mopacity; [[maybe_unused]] auto opac = opacity; + [[maybe_unused]] const auto n_nubinsd = n_nubins; + [[maybe_unused]] const auto &nu_binsd = nu_bins; + [[maybe_unused]] const auto dlnud = dlnu; if constexpr (FT == FrequencyType::gray) { emis = mopac.Emissivity(rho, temp, gmoded); } else if constexpr (FT == FrequencyType::multigroup) { - emis = opac.Emissivity(rho, temp); + // NOTE: 'emis = opac.Emissivity(rho, temp);' may not integrate + emis = 0.0; + for (int n = 0; n < n_nubinsd; n++) { + const Real dnu = dlnud * nu_binsd(n); + emis += opac.EmissivityPerNu(rho, temp, nu_binsd(n)) * dnu; + } } vmesh(b, fj::fleck_factor(), k, j, i) = 1.0 / (1.0 + (4.0 * emis / (rho * cv * temp)) * dt); diff --git a/src/jaybenne/transport.cpp b/src/jaybenne/transport.cpp index 1da5127..f0c5cd4 100644 --- a/src/jaybenne/transport.cpp +++ b/src/jaybenne/transport.cpp @@ -37,6 +37,7 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d auto &resolved_pkgs = pm->resolved_packages; auto &jb_pkg = pm->packages.Get("jaybenne"); const Real h = jb_pkg->template Param("planck_constant"); + const Real hinv = 1.0 / h; auto &eos = jb_pkg->template Param("eos_d"); Opacity opacity; Scattering scattering; @@ -98,6 +99,7 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d // frequency data, needed for multigroup [[maybe_unused]] const auto hd = h; + [[maybe_unused]] const auto hinvd = hinv; [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto dlnud = dlnu; [[maybe_unused]] const auto nu_binsd = nu_bins; @@ -164,8 +166,8 @@ TaskStatus TransportPhotons(MeshData *md, const Real t_start, const Real d const Real &rho = vmesh(b, fjh::density(), kp, jp, ip); const Real &sie = vmesh(b, fjh::sie(), kp, jp, ip); const Real temp = eost.TemperatureFromDensityInternalEnergy(rho, sie); - ss = scatter.TotalScatteringCoefficient(rho, temp, ee); - aa = opac.AbsorptionCoefficient(rho, temp, ee); + ss = scatter.TotalScatteringCoefficient(rho, temp, hinvd * ee); + aa = opac.AbsorptionCoefficient(rho, temp, hinvd * ee); } // reset collision indicators diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index 1bb2f99..ff06bd8 100644 --- a/src/jaybenne/transport_ddmc.cpp +++ b/src/jaybenne/transport_ddmc.cpp @@ -47,6 +47,7 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // data needed for multigroup frequency sampling const Real h = jb_pkg->template Param("planck_constant"); + const Real hinv = 1.0 / h; const Real sb = jb_pkg->template Param("boltzmann"); int n_nubins = JaybenneNull(); Real dlnu = JaybenneNull(); @@ -110,6 +111,7 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // frequency data, needed for multigroup [[maybe_unused]] const auto hd = h; + [[maybe_unused]] const auto hinvd = hinv; [[maybe_unused]] const auto sbd = sb; [[maybe_unused]] const auto n_nubinsd = n_nubins; [[maybe_unused]] const auto dlnud = dlnu; @@ -179,8 +181,8 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R rho = vmesh(b, fjh::density(), kp, jp, ip); const Real &sie = vmesh(b, fjh::sie(), kp, jp, ip); temp = eost.TemperatureFromDensityInternalEnergy(rho, sie); - ss = scatter.TotalScatteringCoefficient(rho, temp, ee); - aa = opac.AbsorptionCoefficient(rho, temp, ee); + ss = scatter.TotalScatteringCoefficient(rho, temp, hinvd * ee); + aa = opac.AbsorptionCoefficient(rho, temp, hinvd * ee); } // reset collision indicators @@ -193,6 +195,10 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R Real e_abs = 0.0; if (is_ddmc_step) { + + // sample if particle is undergoing an elastic event + const bool is_elastic = rng_gen.drand() > aa / (ss + aa); + // Update cell of particle swarm_d.Xtoijk(x, y, z, ip, jp, kp); @@ -204,20 +210,35 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R 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_hi_face_prob(), kp, jp, ip); - const Real &Px_u = - vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), kp, jp, ip + 1); - const Real &Py_l = - multi_d ? vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), kp, jp, ip) : 0.0; - const Real &Py_u = - multi_d ? vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), kp, jp + 1, ip) + // get face probabilities, if no absorption, use per-group values + const Real Px_l = + is_elastic ? 1.0 / (3.0 * ss * dx_i) + : vmesh(b, TE::F1, fj::ddmc_hi_face_prob(), kp, jp, ip); + const Real Px_u = + is_elastic ? 1.0 / (3.0 * ss * dx_i) + : vmesh(b, TE::F1, fj::ddmc_lo_face_prob(), kp, jp, ip + 1); + const Real Py_l = + multi_d ? (is_elastic + ? 1.0 / (3.0 * ss * dx_j) + : vmesh(b, TE::F2, fj::ddmc_hi_face_prob(), kp, jp, ip)) : 0.0; - const Real &Pz_l = - three_d ? vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), kp, jp, ip) : 0.0; - const Real &Pz_u = - three_d ? vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), kp + 1, jp, ip) + const Real Py_u = + multi_d + ? (is_elastic + ? 1.0 / (3.0 * ss * dx_j) + : vmesh(b, TE::F2, fj::ddmc_lo_face_prob(), kp, jp + 1, ip)) + : 0.0; + const Real Pz_l = + three_d ? (is_elastic + ? 1.0 / (3.0 * ss * dx_k) + : vmesh(b, TE::F3, fj::ddmc_hi_face_prob(), kp, jp, ip)) : 0.0; + const Real Pz_u = + three_d + ? (is_elastic + ? 1.0 / (3.0 * ss * dx_k) + : vmesh(b, TE::F3, fj::ddmc_lo_face_prob(), kp + 1, jp, ip)) + : 0.0; // store old cell indices (this is only needed for multigroup) const int ip_old = ip; @@ -269,7 +290,8 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // This is hopefully a minor error. // particle must have leaked if nothing else - if (!(is_absorbed || is_scattered || is_census || is_rejected)) { + if (!(is_absorbed || is_scattered || is_census || is_rejected || + is_elastic)) { // only one index should be +/-1 of the current index const int ip_u = (ip_old == ip + 1 ? ip_old : ip); From c7566c8e6f9efc369bb98554946504ef49a4ec8a Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Mon, 22 Jun 2026 18:48:35 -0600 Subject: [PATCH 13/15] + Fix recently-added emission energy error in sourcing.cpp. + Suppress (for now) population control in multigroup mode. + Update metric for elastic-scattering based DDMC leakage. Note: the code with the new metric is motivated in a comment. --- src/jaybenne/jaybenne.cpp | 7 +++++-- src/jaybenne/sourcing.cpp | 5 ++--- src/jaybenne/transport_ddmc.cpp | 7 ++++++- 3 files changed, 13 insertions(+), 6 deletions(-) diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 9433d67..9a26ba1 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -178,8 +178,11 @@ TaskCollection RadiationStep(Mesh *pmesh, const SimTime &tm, const Real dt) { auto update_fluid = tl.AddTask(eval_rad, jaybenne::UpdateFluid, base.get()); // Control particle population - auto control_pop = tl.AddTask(update_fluid, jaybenne::ControlPopulation, base.get(), - ncycle, ncycle_out); + // TODO: implement population control in multigroup mode + if (fd == FrequencyType::gray) { + auto control_pop = tl.AddTask(update_fluid, jaybenne::ControlPopulation, base.get(), + ncycle, ncycle_out); + } // TODO: Defrag particles? Verify parth swarm defrag mechanics before uncommenting // auto defrag_pop = tl.AddTask(control_pop, jaybenne::DefragParticles, base.get()); diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 6e3f831..62e331b 100644 --- a/src/jaybenne/sourcing.cpp +++ b/src/jaybenne/sourcing.cpp @@ -203,10 +203,9 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) { opac.EmissivityPerNu(rho, temp, nu_binsd(n)) * dnu + vmesh(b, fj::emission_cdf(n - 1), k, j, i); } - emis = 0.0; + // Get total emissivity (before normalizing the CDF) + emis = vmesh(b, fj::emission_cdf(n_nubinsd - 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_nubinsd - 1), k, j, i); diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index ff06bd8..7eb19a1 100644 --- a/src/jaybenne/transport_ddmc.cpp +++ b/src/jaybenne/transport_ddmc.cpp @@ -197,7 +197,12 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R if (is_ddmc_step) { // sample if particle is undergoing an elastic event - const bool is_elastic = rng_gen.drand() > aa / (ss + aa); + // NOTE: this uses the fact that the number of random walks in a cell scales + // like tau^2, and at each event the probability of an elastic scatter is ss + // / (ss + aa). + const Real tau_min = dx_push * (ss + aa); + const bool is_elastic = + (rng_gen.drand() < std::pow(ss / (ss + aa), tau_min * tau_min)); // Update cell of particle swarm_d.Xtoijk(x, y, z, ip, jp, kp); From 62a1289966234be0a63d05ef00eed343811408b0 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 23 Jun 2026 17:52:57 -0600 Subject: [PATCH 14/15] Remove outscatter probability from leakage sample parameter (bug). + Add is_leaked parameter to DDMC step arguments. + Move face temperature evaluation out of leakage integrals. + Remove erroneous MG if-check around population control. --- src/jaybenne/ddmc_mg_utils.hpp | 24 ++++++------ src/jaybenne/jaybenne.cpp | 7 +--- src/jaybenne/transport_ddmc.cpp | 63 ++++++++++++++++++++++---------- src/jaybenne/transport_utils.hpp | 8 ++-- 4 files changed, 64 insertions(+), 38 deletions(-) diff --git a/src/jaybenne/ddmc_mg_utils.hpp b/src/jaybenne/ddmc_mg_utils.hpp index 993b013..0a6c157 100644 --- a/src/jaybenne/ddmc_mg_utils.hpp +++ b/src/jaybenne/ddmc_mg_utils.hpp @@ -46,6 +46,9 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args // define extrapolation distance (Habetler & Matkowski 1975) constexpr Real lam_ext = 0.7104; + // evaluate a face temperature for Planck evaluation + const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); + // initialize unnnormalized Planck integral and probability Real planck_sum = 0.0; Real leak_sum = 0.0; @@ -78,9 +81,6 @@ calc_ddmc_mg_leak_numdenom(const OP &abs, const SC &sct, const ddmc_mg_leak_args // calculate per-nu-bin leakage const Real Pg = 2.0 / (3.0 * (mtau_l + mtau_u)); - // evaluate a face temperature (note this is typically a 4-norm average) - const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); - // convert to energy units for Planck integral const Real ee = dmg.hd * nu_bins(n); const Real dee = ee * dmg.dlnu; @@ -127,6 +127,9 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &s // define extrapolation distance (Habetler & Matkowski 1975) constexpr Real lam_ext = 0.7104; + // evaluate a face temperature for Planck evaluation + const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); + Real leak_sum = 0.0; Real planck_sum = 0.0; @@ -159,9 +162,6 @@ KOKKOS_FORCEINLINE_FUNCTION Real sample_leakage_group(const OP &abs, const SC &s // calculate per-nu-bin leakage const Real Pg = 2.0 / (3.0 * (mtau_l + mtau_u)); - // evaluate a face temperature (note this is typically a 4-norm average) - const Real temp_f = std::max(dmg.temp_l, dmg.temp_u); - // convert to energy units for Planck integral const Real ee = dmg.hd * nu_bins(n); const Real dee = ee * dmg.dlnu; @@ -238,9 +238,9 @@ calc_ddmc_mg_probs(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc, // sample out-scatter IMC group // NOTE(MGDDMC): this routine is assuming Kirchhoff's Law for emissivity (LTE) template -KOKKOS_FORCEINLINE_FUNCTION Real -sample_ddmc2imc_outscatter(const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc, - const ParArray1D &nu_bins, RngGen &rng_gen) { +KOKKOS_FORCEINLINE_FUNCTION Real sample_ddmc2imc_outscatter( + const OP &abs, const SC &sct, const ddmc_mg_cell_args &dmgc, + const ParArray1D &nu_bins, RngGen &rng_gen, const bool stay_in = false) { Real scat_out_tot_sum = 0.0; @@ -252,7 +252,8 @@ sample_ddmc2imc_outscatter(const OP &abs, const SC &sct, const ddmc_mg_cell_args const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu_bins(n)); // check group exclusion - if (!(dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc)) { + const bool is_ddmc_grp = dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc; + if (stay_in ? is_ddmc_grp : !is_ddmc_grp) { // get an unnormalized, non-dimensional Planck integral over group const Real ee = dmgc.hd * nu_bins(n); const Real dee = ee * dmgc.dlnu; @@ -275,7 +276,8 @@ sample_ddmc2imc_outscatter(const OP &abs, const SC &sct, const ddmc_mg_cell_args const Real aa = abs.AbsorptionCoefficient(dmgc.rho, dmgc.temp, nu_bins(n)); // check group exclusion - if (!(dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc)) { + const bool is_ddmc_grp = dmgc.dx_min * (ss + aa) > dmgc.tau_ddmc; + if (stay_in ? is_ddmc_grp : !is_ddmc_grp) { // convert to energy units for Planck integral const Real ee = dmgc.hd * nu_bins(n); diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 9a26ba1..9433d67 100644 --- a/src/jaybenne/jaybenne.cpp +++ b/src/jaybenne/jaybenne.cpp @@ -178,11 +178,8 @@ TaskCollection RadiationStep(Mesh *pmesh, const SimTime &tm, const Real dt) { auto update_fluid = tl.AddTask(eval_rad, jaybenne::UpdateFluid, base.get()); // Control particle population - // TODO: implement population control in multigroup mode - if (fd == FrequencyType::gray) { - auto control_pop = tl.AddTask(update_fluid, jaybenne::ControlPopulation, base.get(), - ncycle, ncycle_out); - } + auto control_pop = tl.AddTask(update_fluid, jaybenne::ControlPopulation, base.get(), + ncycle, ncycle_out); // TODO: Defrag particles? Verify parth swarm defrag mechanics before uncommenting // auto defrag_pop = tl.AddTask(control_pop, jaybenne::DefragParticles, base.get()); diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index 7eb19a1..19fb361 100644 --- a/src/jaybenne/transport_ddmc.cpp +++ b/src/jaybenne/transport_ddmc.cpp @@ -200,9 +200,12 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // NOTE: this uses the fact that the number of random walks in a cell scales // like tau^2, and at each event the probability of an elastic scatter is ss // / (ss + aa). - const Real tau_min = dx_push * (ss + aa); - const bool is_elastic = - (rng_gen.drand() < std::pow(ss / (ss + aa), tau_min * tau_min)); + bool is_elastic = false; + if constexpr (FT == FrequencyType::multigroup) { + const Real tau_min = dx_push * (ss + aa); + is_elastic = + rng_gen.drand() < std::pow(ss / (ss + aa), tau_min * tau_min); + } // Update cell of particle swarm_d.Xtoijk(x, y, z, ip, jp, kp); @@ -254,21 +257,27 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R Real aa_g = aa; Real gm_g = 0.0; if constexpr (FT == FrequencyType::multigroup) { - // clang-format off - ddmc_mg_cell_args dmgc{n_nubinsd, - dlnud, - hd, - sbd, - tau_ddmc, - dx_push, - rho, - temp}; - // clang-format on - const auto aagm = calc_ddmc_mg_probs(opac, scatter, dmgc, nu_binsd); - aa_g = aagm.first; - gm_g = aagm.second; + if (is_elastic) { + aa_g = 0.0; + } else { + // clang-format off + ddmc_mg_cell_args dmgc{n_nubinsd, + dlnud, + hd, + sbd, + tau_ddmc, + dx_push, + rho, + temp}; + + const auto aagm = calc_ddmc_mg_probs(opac, scatter, dmgc, nu_binsd); + aa_g = aagm.first; + gm_g = aagm.second; + } } + bool is_leaked = false; + // create DDMC step argument list // clang-format off ddmc_step_args dia{ // constants @@ -281,7 +290,7 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R // updated by push t, x, y, z, vx, vy, vz, ip, jp, kp, ww, fraction, e_abs, - is_absorbed, is_scattered, is_census}; + is_absorbed, is_scattered, is_census, is_leaked}; // clang-format on // check for IMC-DDMC albedo rejection if particle arrived from IMC region @@ -290,13 +299,29 @@ TaskStatus TransportPhotons_DDMC(MeshData *md, const Real t_start, const R if (!is_rejected) ptcl_ddmc_step(dia, cutoff); if constexpr (FT == FrequencyType::multigroup) { + + if (is_census && !is_elastic) { + // clang-format off + ddmc_mg_cell_args dmgc{n_nubinsd, + dlnud, + hd, + sbd, + tau_ddmc, + dx_push, + rho, + temp}; + // clang-format on + // the final argument tells it to stay in DDMC groups + ee = sample_ddmc2imc_outscatter(opac, scatter, dmgc, nu_bins, rng_gen, + true); + } + // NOTE: these ddmc_mg_leak_args do not use adjacent cell, so // the face CDF does not sum here to ddmc_(hi|lo)_face_prob. // This is hopefully a minor error. // particle must have leaked if nothing else - if (!(is_absorbed || is_scattered || is_census || is_rejected || - is_elastic)) { + if (is_leaked && !is_elastic) { // only one index should be +/-1 of the current index const int ip_u = (ip_old == ip + 1 ? ip_old : ip); diff --git a/src/jaybenne/transport_utils.hpp b/src/jaybenne/transport_utils.hpp index ddd50cf..5d34b10 100644 --- a/src/jaybenne/transport_utils.hpp +++ b/src/jaybenne/transport_utils.hpp @@ -118,6 +118,7 @@ struct ddmc_step_args { bool &is_absorbed; // indicator for absorption in the step bool &is_scattered; // indicator for scattering in the step bool &is_census; // indicator for end of census + bool &is_leaked; // indicator that leakage occurred between cells }; KOKKOS_FORCEINLINE_FUNCTION @@ -295,9 +296,10 @@ void ptcl_ddmc_step(ddmc_step_args dia, const double cutoff) { } else if (xi < an_abs + sct_out + leak_tot) { // TODO(RTW): only sample direction if adjacent cell is below tau_ddmc + dia.is_leaked = true; - // particle will leak to an adjacent cell - const Real xim = xi - an_abs; + // particle will leak to an adjacent cell, reduce sample to leakage prob. + const Real xim = xi - an_abs - sct_out; if (xim < leakx_l) { // leak in negative x/X1 direction dia.ip -= 1; @@ -343,7 +345,7 @@ void ptcl_ddmc_step(ddmc_step_args dia, const double cutoff) { dia.y = dia.yl + 0.5 * dy; // sample direction sample_face_iso_dir(-dia.vv, dia.rng_gen, dia.vz, dia.vx, dia.vy); - } else if (xim <= leak_tot) { + } else if (xim <= leak_tot + rmin) { // leak in positive z/X3 direction dia.kp += dia.three_d; // update position From bebe963e757015994053140e7e9c84ea8fed02b8 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Wed, 24 Jun 2026 13:04:29 -0600 Subject: [PATCH 15/15] Add power-law form of stepdiff problem input with DDMC. + Use pure (elastic) scattering power-law. + Add analytic solution comparison to CI workflows and CI-launcher. --- .github/workflows/ci.yml | 2 + inputs/stepdiff_mg_plaw_ddmc.in | 96 +++++++++++++++++++++++++++++ tst/launch_ci_runner.py | 5 +- tst/stepdiff.py | 2 +- tst/stepdiff_mg_plaw.py | 106 ++++++++++++++++++++++++++++++++ 5 files changed, 209 insertions(+), 2 deletions(-) create mode 100644 inputs/stepdiff_mg_plaw_ddmc.in create mode 100755 tst/stepdiff_mg_plaw.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f3e3072..75ede63 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -141,3 +141,5 @@ jobs: ./stepdiff_smr.py --executable $JAYBENNE_EXECUTABLE \ --input ../inputs/stepdiff_smr_hybrid.in --use_mpiexec \ --mpi_nthreads 8 --mpi_oversubscribe + ./stepdiff_mg_plaw.py --executable $JAYBENNE_EXECUTABLE \ + --input ../inputs/stepdiff_mg_plaw_ddmc.in --use_mpiexec \ No newline at end of file diff --git a/inputs/stepdiff_mg_plaw_ddmc.in b/inputs/stepdiff_mg_plaw_ddmc.in new file mode 100644 index 0000000..13d6436 --- /dev/null +++ b/inputs/stepdiff_mg_plaw_ddmc.in @@ -0,0 +1,96 @@ +# ======================================================================================== +# (C) (or copyright) 2023-2024. 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. +# ======================================================================================== + + +problem_id = stepdiff + + +refinement = none + +nx1 = 100 +x1min = -0.5 +x1max = 0.5 +ix1_bc = outflow +ox1_bc = outflow + +nx2 = 1 +x2min = -0.5 +x2max = 0.5 +ix2_bc = periodic +ox2_bc = periodic + +nx3 = 1 +x3min = -0.5 +x3max = 0.5 +ix3_bc = periodic +ox3_bc = periodic + + +ix1_bc = jaybenne_reflecting +ox1_bc = jaybenne_reflecting +ix2_bc = periodic +ox2_bc = periodic +ix3_bc = periodic +ox3_bc = periodic + + +nx1 = 100 +nx2 = 1 +nx3 = 1 + + +tlim = 3.335641e-10 +integrator = rk1 + + +use_ddmc = true +num_particles = 100000 +dt = 3.335641e-11 +n_nubins = 16 +numin = 1.e12 +numax = 1.e17 +do_emission = false +do_feedback = false +transport_model = zone_size +tracking_algo = history # event +source_strategy = uniform +seed = 349857 + + +frequency_type = multigroup +specific_heat = 1.0e8 # erg/K/g +initial_density = 1.0 # g cm^-3 +initial_temperature = 1.0e5 # K +initial_radiation = thermal + + +opacity_model = none + + +opacity_model = powerlaw +kappa0 = 0.75e3 +rho_exp = 0.0 +temp_exp = 0.0 +nu_exp = -0.25 +nu_ref = 8.0584e15 +nu_off = 4.5316e14 +rho_ref = 1.0 +temp_ref = 1.0 + + +file_type = hdf5 +dt = 3.335641e-10 +variables = field.material.density, & + field.material.sie, & + field.material.internal_energy, & + field.jaybenne.energy_tally diff --git a/tst/launch_ci_runner.py b/tst/launch_ci_runner.py index b97f21b..52fdcb7 100755 --- a/tst/launch_ci_runner.py +++ b/tst/launch_ci_runner.py @@ -112,7 +112,10 @@ def run_tests_in_temp_dir(pr_number, head_repo, head_ref, output_dir): + " --restart ./stepdiff.out1.00001.rhdf --use_mpiexec" + " && ./stepdiff_smr.py --executable " + os.path.join(build_dir, "mcblock") - + " --input ../inputs/stepdiff_smr_hybrid.in --use_mpiexec --mpi_nthreads 8", + + " --input ../inputs/stepdiff_smr_hybrid.in --use_mpiexec --mpi_nthreads 8" + + " && ./stepdiff_mg_plaw.py --executable " + + os.path.join(build_dir, "mcblock") + + " --input ../inputs/stepdiff_mg_plaw_ddmc.in --use_mpiexec", ] ret = subprocess.run(test_command, check=True) diff --git a/tst/stepdiff.py b/tst/stepdiff.py index 24e145d..da6d4d9 100755 --- a/tst/stepdiff.py +++ b/tst/stepdiff.py @@ -30,7 +30,7 @@ modified_inputs["parthenon/meshblock/nx1"] = 128 # -- Analytic solution -tau = 1.000692e-7 +tau = 1.000692e-7 # = 3 * sigma_t / c ur0 = 7.5646e5 shift = 0.5 diff --git a/tst/stepdiff_mg_plaw.py b/tst/stepdiff_mg_plaw.py new file mode 100755 index 0000000..b1d7ca1 --- /dev/null +++ b/tst/stepdiff_mg_plaw.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python +# ======================================================================================== +# (C) (or copyright) 2023-2024. 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. +# ======================================================================================== + +import sys + +sys.dont_write_bytecode = True + +import os +import regression_test as rt +import numpy as np +from scipy.special import erf + +# -- constants +hp = 6.626e-27 #[erg-s] +kb = 1.381e-16 #[erg/K] +cl = 2.998e10 #[cm/s] + +# -- helper functions +def planck(T, nu): + x = hp * nu / (kb * T) + efac = np.exp(-x) + return x * x * x * efac / (1.0 - efac) + +class Plaw_Opac: + def __init__(self, kappa0_, rho_exp_, temp_exp_, nu_exp_, + nu_ref_, nu_off_, rho_ref_, rho_off_, + temp_ref_, temp_off_): + self.kappa0 = kappa0_ + self.rho_exp = rho_exp_ + self.temp_exp = temp_exp_ + self.nu_exp = nu_exp_ + self.nu_ref = nu_ref_ + self.nu_off = nu_off_ + self.rho_ref = rho_ref_ + self.rho_off = rho_off_ + self.temp_ref = temp_ref_ + self.temp_off = temp_off_ + + def plaw_opac(self, rho, T, nu): + xrho = (rho + self.rho_off) / self.rho_ref + xT = (T + self.temp_off) / self.temp_ref + xnu = (nu + self.nu_off) / self.nu_ref + return self.kappa0 * xrho**self.rho_exp * xT**self.temp_exp * xnu**self.nu_exp + +# -- parser +parser = rt.get_default_parser() +args = parser.parse_args() + +modified_inputs = {} +modified_inputs["parthenon/mesh/nx1"] = 128 +modified_inputs["parthenon/meshblock/nx1"] = 128 + +# -- frequency grid +n_nubins = 16 +numin = 1.e12 #[Hz] +numax = 1.e17 #[Hz] +dlnu = np.log(numax / numin) / n_nubins +nu_grid = np.array([numin * np.exp((ig + 0.5) * dlnu) for ig in range(n_nubins)]) +dnu = nu_grid * dlnu + +# -- Analytic solution +nur = nu_grid[12] +nuo = nu_grid[8] +pops = Plaw_Opac(0.75e3, 0.0, 0.0, -0.25, nur, nuo, 1.0, 0.0, 1.0, 0.0) +tau = 3.0 * pops.plaw_opac(1.0, 1.0, nu_grid) / cl # = 3 * sigma_t / c +T0 = 1e5 +plnk = planck(T0, nu_grid) * dnu +plnk /= np.sum(plnk) +ur0 = 7.5646e5 * plnk +shift = 0.5 + + +def urnu_solution(t, x, y, z, ig): + return ( + ur0[ig] + / 2.0 + * ( + erf(((x + shift) + 0.5) / (2.0 * np.sqrt(t / tau[ig]))) + - erf(((x + shift) - 0.5) / (2.0 * np.sqrt(t / tau[ig]))) + ) + ) + +def ur_solution(t, x, y, z): + ur = sum([urnu_solution(t, x, y, z, ig) for ig in range(n_nubins)]) + return ur + +code = rt.analytic_comparison( + args=args, + variables=["field.jaybenne.energy_tally"], + solutions=[ur_solution], + modified_inputs=modified_inputs, + tolerance=0.05, +) + +sys.exit(code)