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);
+    auto 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;
@@ -115,9 +127,21 @@ TaskCollection RadiationStep(Mesh *pmesh, const Real t_start, const Real dt) {
       Kokkos::Profiling::pushRegion("Jaybenne::TransportLoop");
       return TaskStatus::complete;
     });
-    auto transport =
-        use_ddmc ? itl.AddTask(time_start, TransportPhotons_DDMC, base.get(), t_start, dt)
-                 : itl.AddTask(time_start, TransportPhotons, base.get(), t_start, dt);
+    auto transport = time_start;
+    if (fd == FrequencyType::gray) {
+      transport =
+          use_ddmc ? itl.AddTask(time_start, TransportPhotons_DDMC<FrequencyType::gray>,
+                                 base.get(), t_start, dt)
+                   : itl.AddTask(time_start, TransportPhotons<FrequencyType::gray>,
+                                 base.get(), t_start, dt);
+    } else if (fd == FrequencyType::multigroup) {
+      transport =
+          use_ddmc
+              ? itl.AddTask(time_start, TransportPhotons_DDMC<FrequencyType::multigroup>,
+                            base.get(), t_start, dt)
+              : itl.AddTask(time_start, TransportPhotons<FrequencyType::multigroup>,
+                            base.get(), t_start, dt);
+    }
     auto reset_comms = itl.AddTask(transport, MeshResetCommunication, base.get());
     auto send = itl.AddTask(reset_comms, MeshSend, base.get());
     auto receive = itl.AddTask(transport | send, MeshReceive, base.get());
@@ -151,12 +175,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 +197,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 +246,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 +279,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 +627,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,24 @@ 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
+template <FrequencyType FT>
 TaskStatus TransportPhotons(MeshData<Real> *md, const Real t_start, const Real dt);
+template <FrequencyType FT>
 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;