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/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_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/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/jaybenne/ddmc_mg_utils.hpp b/src/jaybenne/ddmc_mg_utils.hpp new file mode 100644 index 0000000..0a6c157 --- /dev/null +++ b/src/jaybenne/ddmc_mg_utils.hpp @@ -0,0 +1,304 @@ +#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 &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_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 + 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 &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 ParArray1D &nu_bins, const bool &use_lo) { + + // 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; + int nnubins_used = 0; + + // integrate + for (int n = 0; n < dmg.n_nubins; ++n) { + + // 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); + 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_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; + + // calculate per-nu-bin leakage + const Real Pg = 2.0 / (3.0 * (mtau_l + mtau_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); + + // aggregate values + planck_sum += bg; + leak_sum += bg * Pg; + } + } + + PARTHENON_REQUIRE(nnubins_used ? planck_sum > 0.0 : true, "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 ParArray1D &nu_bins, + const bool &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; +} + +//---------------------------------------------------------------------------------------- +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, nu_bins, 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; + + // 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; + + // sample + const Real rand1 = leak_tot_sum * rng_gen.drand(); + Real ee_sampled = -1.0; // poisoned initialization + + // integrate + for (int n = 0; n < dmg.n_nubins; ++n) { + + // 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); + 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_lmin > dmg.tau_ddmc : tau_umin > dmg.tau_ddmc); + + if (use_grp) { + 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)); + + // 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); + + // aggregate values + planck_sum += bg; + leak_sum += bg * Pg; + if (leak_sum > rand1) { + ee_sampled = ee; + break; + } + } + } + + PARTHENON_DEBUG_REQUIRE(ee_sampled > 0.0, "ee_sampled <= 0.0"); + + // return sampled energy value (in units of energy) + return ee_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, + const ParArray1D &nu_bins) { + + // 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) { + + // 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; + + // get an unnormalized, non-dimensional Planck integral over group + const Real bg = midpoint_Planck(dmgc.sbd * dmgc.temp, ee, dee); + + // 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; + } + } + + PARTHENON_DEBUG_REQUIRE(planck_sum > 0.0, "planck_sum = 0: no DDMC groups in DDMC."); + + // 1st entry = DDMC absorption, 2nd = out-scatter probability + 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}; + } +} + +//---------------------------------------------------------------------------------------- +// 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, const bool stay_in = false) { + + Real scat_out_tot_sum = 0.0; + + // integrate total cdf value + for (int n = 0; n < dmgc.n_nubins; ++n) { + + // 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 + 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; + 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 ee_sampled = -1.0; // poisoned initialization + + Real abs_sum = 0.0; + + // find bin for sample + for (int n = 0; n < dmgc.n_nubins; ++n) { + + // 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 + 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); + 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); + + // aggregate values + abs_sum += bg * aa; + + if (abs_sum > rand1) { + ee_sampled = ee; + break; + } + } + } + + return ee_sampled; +} + +} // namespace jaybenne + +#endif // JAYBENNE_DDMC_MG_UTILS_HPP_ diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp index 7ac2186..9433d67 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); @@ -232,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); @@ -324,7 +327,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; @@ -367,18 +371,29 @@ 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 - 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); + // 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 @@ -429,6 +444,12 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData *md, const Real dt) { Scattering scattering; MeanOpacity mopacity; MeanScattering mscattering; + int n_nubins = -1; + 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"); @@ -444,6 +465,18 @@ 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"); + h = jbn->template Param("planck_constant"); + sb = jbn->template Param("boltzmann"); + // 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); @@ -451,8 +484,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( @@ -468,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); @@ -497,8 +538,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; @@ -517,6 +556,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)); @@ -527,9 +568,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) @@ -551,27 +601,64 @@ 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 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; + [[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); 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 + 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_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) = + 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, + dlnud, + hd, + sbd, + tau_ddmcd, + dx_lmin, + dx_umin, + 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, 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, nu_binsd, use_hi); + } }); // set face probabilities in X2 direction @@ -584,7 +671,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)); @@ -595,9 +684,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) @@ -618,28 +717,65 @@ 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 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; + [[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); 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 + 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_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) - 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, + dlnud, + hd, + sbd, + tau_ddmcd, + dx_lmin, + dx_umin, + 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, 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, nu_binsd, use_hi); + } }); } @@ -653,6 +789,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 @@ -664,9 +802,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) @@ -687,28 +835,65 @@ 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 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; + [[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); 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 + 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_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) - 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, + dlnud, + hd, + sbd, + tau_ddmcd, + dx_lmin, + dx_umin, + 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, 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, nu_bins, 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..d52d4aa 100644 --- a/src/jaybenne/planck.hpp +++ b/src/jaybenne/planck.hpp @@ -18,13 +18,31 @@ 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; + 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); + } +} + //---------------------------------------------------------------------------------------- //! \fn Real sample_Planck_energy //! \brief Efficiently samples the Planck distribution for particle energy //! 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; @@ -46,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/sample_ddmc_bface.cpp b/src/jaybenne/sample_ddmc_bface.cpp index ffadb12..06b3bea 100644 --- a/src/jaybenne/sample_ddmc_bface.cpp +++ b/src/jaybenne/sample_ddmc_bface.cpp @@ -96,11 +96,11 @@ 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 - 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 +184,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 +217,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 +319,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 +366,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 +413,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..9547847 100644 --- a/src/jaybenne/scattering.hpp +++ b/src/jaybenne/scattering.hpp @@ -16,17 +16,68 @@ 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 &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, + const ParArray1D &nu_bins, + 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; + } + } + + // reset particle frequency + psa.ee = csa.hd * nu_bins(n); + } } } // namespace jaybenne diff --git a/src/jaybenne/sourcing.cpp b/src/jaybenne/sourcing.cpp index 765492c..62e331b 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"); @@ -76,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 = @@ -146,45 +154,58 @@ 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 &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; [[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; + // 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) { 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; + // 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++) { - const Real dnu = dlnu * numind * std::exp((n + 0.5) * dlnu); - // Get total emissivity - emis += vmesh(b, fj::emission_cdf(n), k, j, i); // Normalize emission CDF vmesh(b, fj::emission_cdf(n), k, j, i) /= vmesh(b, fj::emission_cdf(n_nubinsd - 1), k, j, i); @@ -266,13 +287,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 &kboltd = kbolt; + [[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); @@ -311,19 +332,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 from CDF + // 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; } } - 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(g); } if constexpr (ST == SourceType::emission) { diff --git a/src/jaybenne/transport.cpp b/src/jaybenne/transport.cpp index 1925e85..f0c5cd4 100644 --- a/src/jaybenne/transport.cpp +++ b/src/jaybenne/transport.cpp @@ -37,18 +37,27 @@ 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; 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"); @@ -90,9 +99,10 @@ 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 hinvd = hinv; [[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); @@ -150,15 +160,14 @@ 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) { 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 @@ -221,28 +230,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, hd}; + // clang-format on + + // invoke frequency-dependent scattering kernel + scatter_kernel(vmesh, csa, nu_binsd, psa); } } diff --git a/src/jaybenne/transport_ddmc.cpp b/src/jaybenne/transport_ddmc.cpp index be28bb6..19fb361 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,46 @@ 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 hinv = 1.0 / h; + const Real sb = jb_pkg->template Param("boltzmann"); + int n_nubins = 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"); + // 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"); + } + // Create SparsePack static auto desc = - MakePackDescriptor( resolved_pkgs.get()); auto vmesh = desc.GetPack(md); @@ -81,6 +108,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 hinvd = hinv; + [[maybe_unused]] const auto sbd = sb; + [[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); const Real &dx_j = coords.template Dxc(0, 0, 0); @@ -131,8 +167,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, hinvd * ee); + aa = opac.AbsorptionCoefficient(rho, temp, hinvd * ee); + } // reset collision indicators bool is_absorbed = false; @@ -144,6 +195,18 @@ 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 + // 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). + 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); @@ -155,31 +218,79 @@ 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_face_prob(), kp, jp, ip); - const Real &Px_u = vmesh(b, TE::F1, fj::ddmc_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; - const Real &Py_u = - multi_d ? vmesh(b, TE::F2, fj::ddmc_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; - const Real &Pz_u = - three_d ? vmesh(b, TE::F3, fj::ddmc_face_prob(), kp + 1, jp, ip) : 0.0; + // 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 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; + 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) { + 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 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, // 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 @@ -187,6 +298,82 @@ 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_leaked && !is_elastic) { + + // 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, + dlnud, + hd, + sbd, + tau_ddmc, + dx_f, + dx_f, + 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, nu_binsd, use_lo, rng_gen); + } + } + } else { // push particle @@ -253,10 +440,47 @@ 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, + dlnud, + hd, + sbd, + tau_ddmc, + dx_push, + rho, + temp}; + // clang-format on + ee = sample_ddmc2imc_outscatter(opac, scatter, dmgc, nu_bins, rng_gen); + + // resample particle direction + sample_vol_iso_dir(psa); + + } else { + + // form cell scattering argument struct + // clang-format off + cell_scat_args csa{b, ip, jp, kp, + ff, aa, ss, + n_nubinsd, hd}; + // clang-format on + + // invoke frequency-dependent scattering kernel + scatter_kernel(vmesh, csa, nu_binsd, psa); + } + } } if (is_census) { diff --git a/src/jaybenne/transport_utils.hpp b/src/jaybenne/transport_utils.hpp index cbc07ec..5d34b10 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 @@ -117,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 @@ -255,8 +257,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,12 +288,18 @@ 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 + 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; @@ -334,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 diff --git a/src/mcblock/mcblock.cpp b/src/mcblock/mcblock.cpp index fc110c4..8a2e340 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); @@ -111,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, @@ -136,16 +135,36 @@ 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") { + // 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->GetOrAddReal("mcblock/absorption", "nu_off", 0.0); + const Real rho_ref = pin->GetReal("mcblock/absorption", "rho_ref"); + 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->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(); + opacity = singularity::photons::NonCGSUnits( + 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( singularity::photons::EPBremss(), time_scale, mass_scale, length_scale, temperature_scale); @@ -155,7 +174,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); @@ -163,20 +181,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( @@ -186,20 +203,35 @@ 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") { + // 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->GetOrAddReal("mcblock/scattering", "nu_off", 0.0); + const Real rho_ref = pin->GetReal("mcblock/scattering", "rho_ref"); + 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->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), + 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 = 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)