Gray vs multigroup

src/CMakeLists.txt

@@ -31,6 +31,12 @@ if(JAYBENNE_STANDALONE_MODE)
   set(JAYBENNE_SCATTERING_OPACITY_TYPE
       "mcblock::Scattering"
       CACHE STRING "Namespace for scattering opacity Singularity type" FORCE)
+  set(JAYBENNE_ABSORPTION_MEAN_OPACITY_TYPE
+      "mcblock::MeanOpacity"
+      CACHE STRING "Namespace for absorption mean opacity Singularity type" FORCE)
+  set(JAYBENNE_SCATTERING_MEAN_OPACITY_TYPE
+      "mcblock::MeanScattering"
+      CACHE STRING "Namespace for scattering mean opacity Singularity type" FORCE)
   set(JAYBENNE_HOST_VARIABLE_HEADER
       "../mcblock/mcblock_variables.hpp"
       CACHE STRING "Path to host code variables header file" FORCE)

src/jaybenne/jaybenne.cpp

@@ -71,6 +71,7 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) {
   auto &jb_pkg = pmesh->packages.Get("jaybenne");
   const auto &max_transport_iterations = jb_pkg->Param<int>("max_transport_iterations");
   const bool &use_ddmc = jb_pkg->Param<bool>("use_ddmc");
+  const auto &fd = jb_pkg->Param<FrequencyType>("frequency_type");
 
   // MeshData subsets
   auto ddmc_field_names = std::vector<std::string>{fj::ddmc_face_prob::name()};
@@ -102,9 +103,20 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) {
 
     // prepare for iterative transport loop
     auto derived = tl.AddTask(none, UpdateDerivedTransportFields, base.get(), dt);
-    auto source = tl.AddTask(
-        derived, jaybenne::SourcePhotons<MeshData<Real>, jaybenne::SourceType::emission>,
-        base.get(), t_start, dt);
+    TaskID source = derived;
+    if (fd == FrequencyType::gray) {
+      auto source = tl.AddTask(
+          derived,
+          jaybenne::SourcePhotons<MeshData<Real>, jaybenne::SourceType::emission,
+                                  FrequencyType::gray>,
+          base.get(), t_start, dt);
+    } else if (fd == FrequencyType::multigroup) {
+      auto source = tl.AddTask(
+          derived,
+          jaybenne::SourcePhotons<MeshData<Real>, jaybenne::SourceType::emission,
+                                  FrequencyType::multigroup>,
+          base.get(), t_start, dt);
+    }
     auto bcs = use_ddmc ? parthenon::AddBoundaryExchangeTasks(source, tl, md_ddmc,
                                                               pmesh->multilevel)
                         : source;
@@ -151,12 +163,13 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) {
 }
 
 //----------------------------------------------------------------------------------------
-//! \fn  StateDescriptor Jaybenne::Initialize
+//! \fn  StateDescriptor Jaybenne::Initialize_impl
 //! \brief Initialize the Jaybenne physics package. This function defines and sets the
 //! parameters associated with Jaybenne, and enrolls the data variables associated with
-//! this physics package.
-std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin, Opacity &opacity,
-                                            Scattering &scattering, EOS &eos) {
+//! this physics package, for everything not related to frequency discretization.
+std::shared_ptr<StateDescriptor>
+Initialize_impl(ParameterInput *pin, EOS &eos,
+                singularity::RuntimePhysicalConstants units) {
   auto pkg = std::make_shared<StateDescriptor>("jaybenne");
 
   // Total number of particles
@@ -172,16 +185,11 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin, Opacity &opacit
                     "Minimum allowable swarm occupancy must be >= 0 and less than 1");
   pkg->AddParam<>("min_swarm_occupancy", min_swarm_occupancy);
 
-  // Frequency range
-  Real numin = pin->GetOrAddReal("jaybenne", "numin", std::numeric_limits<Real>::min());
-  pkg->AddParam<>("numin", numin);
-  Real numax = pin->GetOrAddReal("jaybenne", "numax", std::numeric_limits<Real>::max());
-  pkg->AddParam<>("numax", numax);
-
   // Physical constants
-  const auto units = opacity.GetRuntimePhysicalConstants();
+  // const auto units = opacity.GetRuntimePhysicalConstants();
   pkg->AddParam<>("speed_of_light", units.c);
   pkg->AddParam<>("stefan_boltzmann", units.sb);
+  pkg->AddParam<>("planck_constant", units.h);
 
   // RNG
   bool unique_rank_seeds = pin->GetOrAddBoolean("jaybenne", "unique_rank_seeds", true);
@@ -226,12 +234,6 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin, Opacity &opacit
   // Equation of state model
   pkg->AddParam<>("eos_d", eos.GetOnDevice());
 
-  // Opacity model
-  pkg->AddParam<>("opacity_d", opacity.GetOnDevice());
-
-  // Scattering model
-  pkg->AddParam<>("scattering_d", scattering.GetOnDevice());
-
   // Swarm and swarm variables
   Metadata swarm_metadata({Metadata::Provides, Metadata::None});
   pkg->AddSwarm(photons_swarm_name, swarm_metadata);
@@ -265,6 +267,63 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin, Opacity &opacit
   return pkg;
 }
 
+//----------------------------------------------------------------------------------------
+//! \fn  StateDescriptor Jaybenne::Initialize
+//! \brief Initialize the Jaybenne physics package. This function defines and sets the
+//! parameters associated with Jaybenne, and enrolls the data variables associated with
+//! this physics package, specific to the gray frequency discretization.
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin, MeanOpacity &mopacity,
+                                            MeanScattering &mscattering, EOS &eos) {
+  const auto units = mopacity.GetRuntimePhysicalConstants();
+
+  auto pkg = Initialize_impl(pin, eos, units);
+
+  pkg->AddParam<>("frequency_type", FrequencyType::multigroup);
+
+  // Opacity model
+  pkg->AddParam<>("mopacity_d", mopacity.GetOnDevice());
+
+  // Scattering model
+  pkg->AddParam<>("mscattering_d", mscattering.GetOnDevice());
+
+  return pkg;
+}
+
+//----------------------------------------------------------------------------------------
+//! \fn  StateDescriptor Jaybenne::Initialize
+//! \brief Initialize the Jaybenne physics package. This function defines and sets the
+//! parameters associated with Jaybenne, and enrolls the data variables associated with
+//! this physics package, specific to the multigroup frequency discretization.
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin, Opacity &opacity,
+                                            Scattering &scattering, EOS &eos) {
+  const auto units = opacity.GetRuntimePhysicalConstants();
+
+  auto pkg = Initialize_impl(pin, eos, units);
+
+  // Frequency discretization
+  auto time = units.time;
+  Real numin = pin->GetReal("jaybenne", "numin"); // in Hz
+  pkg->AddParam<>("numin", numin * time);         // in code units
+  Real numax = pin->GetReal("jaybenne", "numax"); // in Hz
+  pkg->AddParam<>("numax", numax * time);         // in code units
+  int n_nubins = pin->GetInteger("jaybenne", "n_nubins");
+  pkg->AddParam<>("n_nubins", n_nubins);
+
+  pkg->AddParam<>("frequency_type", FrequencyType::multigroup);
+
+  // Emission CDF
+  Metadata m_onecopy({Metadata::Cell, Metadata::OneCopy}, std::vector<int>({n_nubins}));
+  pkg->AddField(field::jaybenne::emission_cdf::name(), m_onecopy);
+
+  // Opacity model
+  pkg->AddParam<>("opacity_d", opacity.GetOnDevice());
+
+  // Scattering model
+  pkg->AddParam<>("scattering_d", scattering.GetOnDevice());
+
+  return pkg;
+}
+
 //----------------------------------------------------------------------------------------
 //! \fn  Real Jaybenne::EstimateTimestepMesh
 //! \brief Compute radiation timestep
@@ -556,9 +615,17 @@ TaskStatus EvaluateRadiationEnergy(T *md) {
 //! \brief Initialize radiation based on material temperature and either thermal or
 //!        zero initial radiation.
 void InitializeRadiation(MeshBlockData<Real> *mbd, const bool is_thermal) {
+  auto &jb_pkg = mbd->GetBlockPointer()->packages.Get("jaybenne");
+  const auto &fd = jb_pkg->Param<FrequencyType>("frequency_type");
+
   if (is_thermal) {
-    jaybenne::SourcePhotons<MeshBlockData<Real>, jaybenne::SourceType::thermal>(mbd, 0.0,
-                                                                                0.0);
+    if (fd == FrequencyType::gray) {
+      jaybenne::SourcePhotons<MeshBlockData<Real>, jaybenne::SourceType::thermal,
+                              FrequencyType::gray>(mbd, 0.0, 0.0);
+    } else if (fd == FrequencyType::multigroup) {
+      jaybenne::SourcePhotons<MeshBlockData<Real>, jaybenne::SourceType::thermal,
+                              FrequencyType::multigroup>(mbd, 0.0, 0.0);
+    }
   }
 
   // Call this so radiation field variables are properly initialized for outputs

src/jaybenne/jaybenne.hpp

@@ -50,17 +50,22 @@ namespace jaybenne {
 std::shared_ptr<parthenon::StateDescriptor> Initialize(parthenon::ParameterInput *pin,
                                                        Opacity &opacity,
                                                        Scattering &scattering, EOS &eos);
+std::shared_ptr<parthenon::StateDescriptor> Initialize(parthenon::ParameterInput *pin,
+                                                       MeanOpacity &mopacity,
+                                                       MeanScattering &mscattering,
+                                                       EOS &eos);
 
 // Model enums
 enum class SourceStrategy { uniform, energy };
 enum class SourceType { thermal, emission };
+enum class FrequencyType { gray, multigroup };
 
 // Tasks
 TaskStatus TransportPhotons(MeshData<Real> *md, const Real t_start, const Real dt);
 TaskStatus TransportPhotons_DDMC(MeshData<Real> *md, const Real t_start, const Real dt);
 TaskStatus SampleDDMCBlockFace(MeshData<Real> *md);
 TaskStatus CheckCompletion(MeshData<Real> *md, const Real t_end);
-template <typename T, SourceType ST>
+template <typename T, SourceType ST, FrequencyType FT>
 TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt);
 TaskStatus DefragParticles(MeshBlock *pmb);
 TaskStatus UpdateDerivedTransportFields(MeshData<Real> *md, const Real dt);

src/jaybenne/jaybenne_config.hpp.in

@@ -22,6 +22,8 @@
 using @JAYBENNE_EOS_TYPE@;
 using @JAYBENNE_ABSORPTION_OPACITY_TYPE@;
 using @JAYBENNE_SCATTERING_OPACITY_TYPE@;
+using @JAYBENNE_ABSORPTION_MEAN_OPACITY_TYPE@;
+using @JAYBENNE_SCATTERING_MEAN_OPACITY_TYPE@;
 
 #include "@JAYBENNE_HOST_VARIABLE_HEADER@"
 

src/jaybenne/jaybenne_variables.hpp

@@ -37,6 +37,7 @@ JAYBENNE_FIELD_VARIABLE(field.jaybenne, fleck_factor);
 JAYBENNE_FIELD_VARIABLE(field.jaybenne, ddmc_face_prob);
 JAYBENNE_FIELD_VARIABLE(field.jaybenne, source_ew_per_cell);
 JAYBENNE_FIELD_VARIABLE(field.jaybenne, source_num_per_cell);
+JAYBENNE_FIELD_VARIABLE(field.jaybenne, emission_cdf);
 JAYBENNE_FIELD_VARIABLE(field.jaybenne, energy_delta);
 namespace host {
 typedef HOST_DENSITY density;

src/jaybenne/sourcing.cpp

@@ -22,7 +22,7 @@ namespace jaybenne {
 //! \fn  TaskStatus SourcePhotons
 //! \brief Create photons, either during initialization or during a timestep.
 //! TODO(BRR) modify interface so we don't need t_start, dt for initialization
-template <typename T, SourceType ST>
+template <typename T, SourceType ST, FrequencyType FT>
 TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
   namespace fj = field::jaybenne;
   namespace fjh = field::jaybenne::host;
@@ -31,8 +31,21 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
   auto pm = md->GetParentPointer();
   auto &resolved_pkgs = pm->resolved_packages;
   auto &jb_pkg = pm->packages.Get("jaybenne");
+  const Real h = jb_pkg->template Param<Real>("planck_constant");
   auto &eos = jb_pkg->template Param<EOS>("eos_d");
-  auto &opacity = jb_pkg->template Param<Opacity>("opacity_d");
+  const MeanOpacity *mopacity = nullptr;
+  const Opacity *opacity = nullptr;
+  int n_nubins = -1;
+  Real numin = -1.;
+  Real numax = -1.;
+  if constexpr (FT == FrequencyType::gray) {
+    mopacity = &(jb_pkg->template Param<MeanOpacity>("mopacity_d"));
+  } else if constexpr (FT == FrequencyType::multigroup) {
+    opacity = &(jb_pkg->template Param<Opacity>("opacity_d"));
+    n_nubins = jb_pkg->template Param<int>("n_nubins");
+    numin = jb_pkg->template Param<Real>("numin");
+    numax = jb_pkg->template Param<Real>("numax");
+  }
   auto &do_emission = jb_pkg->template Param<bool>("do_emission");
   auto &source_strategy = jb_pkg->template Param<SourceStrategy>("source_strategy");
   PARTHENON_REQUIRE(source_strategy != SourceStrategy::energy,
@@ -51,7 +64,8 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
   // Create pack
   static auto desc =
       MakePackDescriptor<fjh::density, fjh::sie, fj::fleck_factor, fj::source_ew_per_cell,
-                         fj::source_num_per_cell, fj::energy_delta>(resolved_pkgs.get());
+                         fj::source_num_per_cell, fj::emission_cdf, fj::energy_delta>(
+          resolved_pkgs.get());
   auto vmesh = desc.GetPack(md);
 
   // Indexing and dimensionality
@@ -87,12 +101,37 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
               [[maybe_unused]] const Real &sbd = sb;
               [[maybe_unused]] const Real &vvd = vv;
               [[maybe_unused]] const Real &dtd = dt;
-              [[maybe_unused]] auto &opac = opacity;
+              [[maybe_unused]] auto *opac = opacity;
+              [[maybe_unused]] auto *mopac = mopacity;
+              [[maybe_unused]] const auto numind = numin;
+              [[maybe_unused]] const auto numaxd = numax;
+              [[maybe_unused]] const auto n_nubinsd = n_nubins;
               Real erad = 0.0;
               if constexpr (ST == SourceType::thermal) {
                 erad = (4.0 * sbd / vvd) * std::pow(temp, 4.0) * dv;
               } else if constexpr (ST == SourceType::emission) {
-                const Real emis = opac.Emissivity(rho, temp);
+                Real emis = -1.;
+                if constexpr (FT == FrequencyType::gray) {
+                  emis = opac->Emissivity(rho, temp);
+                } else if constexpr (FT == FrequencyType::multigroup) {
+                  // Construct emission CDF
+                  const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd;
+                  vmesh(b, fj::emission_cdf(0), k, j, i) = opac->EmissivityPerNu(
+                      rho, temp, std::exp(std::log(numind) + 0.5 * dlnu));
+                  for (int n = 1; n < n_nubinsd; n++) {
+                    const Real nu = std::exp(std::log(numind) + (n + 0.5) * dlnu);
+                    vmesh(b, fj::emission_cdf(n), k, j, i) =
+                        opac->EmissivityPerNu(rho, temp, nu) +
+                        vmesh(b, fj::emission_cdf(n - 1), k, j, i);
+                  }
+                  for (int n = 0; n < n_nubinsd; n++) {
+                    // Get total emissivity
+                    emis += vmesh(b, fj::emission_cdf(n), k, j, i);
+                    // Normalize emission CDF
+                    vmesh(b, fj::emission_cdf(n), k, j, i) /=
+                        vmesh(b, fj::emission_cdf(n_nubins - 1), k, j, i);
+                  }
+                }
                 erad = vmesh(b, fj::fleck_factor(), k, j, i) * emis * dv * dtd;
               }
               // Set source_num_per_cell
@@ -156,6 +195,10 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
         auto rng_gen = rng_pool.get_state();
         [[maybe_unused]] const Real &dtd = dt;
         [[maybe_unused]] const Real &t_startd = t_start;
+        [[maybe_unused]] const Real hd = h;
+        [[maybe_unused]] const Real numaxd = numax;
+        [[maybe_unused]] const Real numind = numin;
+        [[maybe_unused]] const int n_nubinsd = n_nubins;
 
         // Starting index and length of particles in this cell
         const int &pstart_idx = prefix_sum(b, cell_idx_1d);
@@ -189,8 +232,23 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
           const Real &rho = vmesh(b, fjh::density(), k, j, i);
           const Real &sie = vmesh(b, fjh::sie(), k, j, i);
           const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie);
-          ppack_r(b, ph::energy(), n) = sample_Planck_energy(rng_gen, sb, temp);
-          ppack_r(b, ph::weight(), n) = vmesh(b, fj::source_ew_per_cell(), k, j, i);
+          if constexpr (FT == FrequencyType::gray) {
+            ppack_r(b, ph::energy(), n) = sample_Planck_energy(rng_gen, sb, temp);
+          } else if constexpr (FT == FrequencyType::multigroup) {
+            // Sample energy from CDF
+            const Real rand = rng_gen.drand();
+            int n;
+            for (n = 0; n < n_nubinsd; n++) {
+              if (vmesh(b, fj::emission_cdf(n), k, j, i) >= rand) {
+                break;
+              }
+            }
+            const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd;
+            const Real nu = std::exp(std::log(numind) + (n + 0.5) * dlnu);
+            ppack_r(b, ph::energy(), n) = hd * nu;
+          }
+          ppack_r(b, ph::weight(), n) =
+              vmesh(b, fj::source_ew_per_cell(), k, j, i) / ppack_r(b, ph::energy(), n);
 
           if constexpr (ST == SourceType::emission) {
             dejbn -= ppack_r(b, ph::weight(), n);
@@ -212,9 +270,22 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
 typedef MeshBlockData<Real> BD;
 typedef MeshData<Real> D;
 typedef SourceType ST;
-template TaskStatus SourcePhotons<BD, ST::thermal>(BD *md, const Real t0, const Real dt);
-template TaskStatus SourcePhotons<BD, ST::emission>(BD *md, const Real t0, const Real dt);
-template TaskStatus SourcePhotons<D, ST::thermal>(D *md, const Real t0, const Real dt);
-template TaskStatus SourcePhotons<D, ST::emission>(D *md, const Real t0, const Real dt);
+typedef FrequencyType FT;
+template TaskStatus SourcePhotons<BD, ST::thermal, FT::gray>(BD *md, const Real t0,
+                                                             const Real dt);
+template TaskStatus SourcePhotons<BD, ST::emission, FT::gray>(BD *md, const Real t0,
+                                                              const Real dt);
+template TaskStatus SourcePhotons<D, ST::thermal, FT::gray>(D *md, const Real t0,
+                                                            const Real dt);
+template TaskStatus SourcePhotons<D, ST::emission, FT::gray>(D *md, const Real t0,
+                                                             const Real dt);
+template TaskStatus SourcePhotons<BD, ST::thermal, FT::multigroup>(BD *md, const Real t0,
+                                                                   const Real dt);
+template TaskStatus SourcePhotons<BD, ST::emission, FT::multigroup>(BD *md, const Real t0,
+                                                                    const Real dt);
+template TaskStatus SourcePhotons<D, ST::thermal, FT::multigroup>(D *md, const Real t0,
+                                                                  const Real dt);
+template TaskStatus SourcePhotons<D, ST::emission, FT::multigroup>(D *md, const Real t0,
+                                                                   const Real dt);
 
 } // namespace jaybenne

src/mcblock/opacity.hpp

@@ -14,6 +14,8 @@
 #define MCBLOCK_OPACITY_HPP_
 
 // Singularity-opac includes
+#include <singularity-opac/photons/mean_opacity_photons.hpp>
+#include <singularity-opac/photons/mean_s_opacity_photons.hpp>
 #include <singularity-opac/photons/opac_photons.hpp>
 #include <singularity-opac/photons/s_opac_photons.hpp>
 
@@ -24,11 +26,15 @@ using Opacity = singularity::photons::impl::Variant<
     singularity::photons::NonCGSUnits<singularity::photons::Gray>,
     singularity::photons::NonCGSUnits<singularity::photons::EPBremss>>;
 
+using MeanOpacity = singularity::photons::MeanOpacity;
+
 // Reduced scattering variant just for jaybenne
 using Scattering = singularity::photons::impl::S_Variant<
     singularity::photons::NonCGSUnitsS<singularity::photons::GrayS>,
     singularity::photons::NonCGSUnitsS<singularity::photons::ThomsonS>>;
 
+using MeanScattering = singularity::photons::MeanSOpacity;
+
 } // namespace mcblock
 
 #endif // MCBLOCK_OPACITY_HPP_