Update gray mode for compatibility with latest singularity-opac.

.gitmodules

@@ -1,9 +1,12 @@
 [submodule "external/parthenon"]
-	path = external/parthenon
-	url = https://github.com/parthenon-hpc-lab/parthenon.git
+  path = external/parthenon
+  url = https://github.com/parthenon-hpc-lab/parthenon.git
+  ignore = dirty
 [submodule "external/singularity-opac"]
-	path = external/singularity-opac
-	url = https://github.com/lanl/singularity-opac.git
+  path = external/singularity-opac
+  url = https://github.com/lanl/singularity-opac.git
+  ignore = dirty
 [submodule "external/singularity-eos"]
-	path = external/singularity-eos
-	url = https://github.com/lanl/singularity-eos.git
+  path = external/singularity-eos
+  url = https://github.com/lanl/singularity-eos.git
+  ignore = dirty

env/bash

@@ -108,11 +108,12 @@ elif [[ $PARTITION == "darwin-volta-x86" ]]; then
     module load cuda/12.3.1
     module load gcc/12.2.0
     module load hdf5/1.12.2
-    module load miniconda3
-    source /projects/opt/centos8/x86_64/miniconda3/py39_4.12.0/bin/activate /usr/projects/jovian/dependencies/python/volta-x86-py39
+    module load miniconda3/py311_23.11.0
+    source /usr/projects/jovian/dependencies/python/volta-x86-py311/bin/activate
     shorten_prompt
     module load cmake/3.26.3
     module list
+    export UCX_TLS=^gdr_copy # claude-opus suggestion for GPU buffer-pinning error (may just be node cn410)
     echo "...setup SUCCEEDED"
 elif [[ $PARTITION == "darwin-a100" ]]; then
     module purge

src/jaybenne/jaybenne.cpp

@@ -448,6 +448,7 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData<Real> *md, const Real dt) {
   Real dlnu = -1.0;
   Real h = -1.0;
   Real sb = -1.0;
+  Real ac = -1.0; // radiation constant times light speed (set below)
   std::vector<Real> nu_grid = JaybenneNull<std::vector<Real>>();
   ParArray1D<Real> nu_bins;
 
@@ -458,6 +459,7 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData<Real> *md, const Real dt) {
   const OpacityAveraging gmode = use_planck ? Planck : Rosseland;
 
   if constexpr (FT == FrequencyType::gray) {
+    ac = 4.0 * (jbn->template Param<Real>("stefan_boltzmann"));
     mopacity = jbn->template Param<MeanOpacity>("mopacity_d");
     mscattering = jbn->template Param<MeanScattering>("mscattering_d");
   } else if constexpr (FT == FrequencyType::multigroup) {
@@ -498,14 +500,16 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData<Real> *md, const Real dt) {
         const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie);
         const Real cv = eos.SpecificHeatFromDensityInternalEnergy(rho, sie);
         Real emis = JaybenneNull<Real>();
+        [[maybe_unused]] const auto acd = ac;
         [[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);
+          const Real T4 = SQR(SQR(temp));
+          emis = mopac.AbsorptionCoefficient(rho, temp, 0, gmoded) * acd * T4;
         } else if constexpr (FT == FrequencyType::multigroup) {
           // NOTE: 'emis = opac.Emissivity(rho, temp);' may not integrate
           emis = 0.0;
@@ -609,10 +613,10 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData<Real> *md, const Real dt) {
           [[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);
+            ss_l = mscatter.ScatteringCoefficient(rho_l, temp_l, 0, gmode2d);
+            aa_l = mopac.AbsorptionCoefficient(rho_l, temp_l, 0, gmode2d);
+            ss_u = mscatter.ScatteringCoefficient(rho_u, temp_u, 0, gmode2d);
+            aa_u = mopac.AbsorptionCoefficient(rho_u, temp_u, 0, gmode2d);
 
             // calculate optical thicknesses from lower and upper cell
             const Real tau_lmin = dx_lmin * (ss_l + aa_l);
@@ -725,10 +729,10 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData<Real> *md, const Real dt) {
             [[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);
+              ss_l = mscatter.ScatteringCoefficient(rho_l, temp_l, 0, gmode2d);
+              aa_l = mopac.AbsorptionCoefficient(rho_l, temp_l, 0, gmode2d);
+              ss_u = mscatter.ScatteringCoefficient(rho_u, temp_u, 0, gmode2d);
+              aa_u = mopac.AbsorptionCoefficient(rho_u, temp_u, 0, gmode2d);
 
               // calculate optical thicknesses from lower and upper cell
               const Real tau_lmin = dx_lmin * (ss_l + aa_l);
@@ -843,10 +847,10 @@ TaskStatus UpdateDerivedTransportFieldsImpl(MeshData<Real> *md, const Real dt) {
             [[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);
+              ss_l = mscatter.ScatteringCoefficient(rho_l, temp_l, 0, gmode2d);
+              aa_l = mopac.AbsorptionCoefficient(rho_l, temp_l, 0, gmode2d);
+              ss_u = mscatter.ScatteringCoefficient(rho_u, temp_u, 0, gmode2d);
+              aa_u = mopac.AbsorptionCoefficient(rho_u, temp_u, 0, gmode2d);
 
               // calculate optical thicknesses from lower and upper cell
               const Real tau_lmin = dx_lmin * (ss_l + aa_l);

src/jaybenne/sourcing.cpp

@@ -190,7 +190,9 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
               } else if constexpr (ST == SourceType::emission) {
                 Real emis = JaybenneNull<Real>();
                 if constexpr (FT == FrequencyType::gray) {
-                  emis = mopac.Emissivity(rho, temp, gmode);
+                  const Real T4 = SQR(SQR(temp));
+                  emis =
+                      mopac.AbsorptionCoefficient(rho, temp, 0, gmoded) * 4.0 * sbd * T4;
                 } else if constexpr (FT == FrequencyType::multigroup) {
                   // Construct emission CDF
                   // calculate bin width (assuming log bin width)

src/mcblock/mcblock.cpp

@@ -112,19 +112,25 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   MeanOpacity mopacity;
   std::string abs_model = pin->GetString("mcblock/absorption", "opacity_model");
   if (frequency_type == FrequencyType::gray) {
+
+    // hard-coded numbers in Hz
+    const std::array<Real, 2> gray_bounds = {time_scale * 1.e12, time_scale * 3.e20};
+
     if (abs_model == "none") {
       auto opac = singularity::photons::Gray(1.e-100);
       mopacity =
           singularity::photons::MeanNonCGSUnits<singularity::photons::MeanOpacityBase>(
-              singularity::photons::MeanOpacityBase(opac, -1, 1, 2, -1, 1, 2), time_scale,
-              mass_scale, length_scale, temperature_scale);
+              singularity::photons::MeanOpacityBase(opac, -1, 1, 2, -1, 1, 2, gray_bounds,
+                                                    1),
+              time_scale, mass_scale, length_scale, temperature_scale);
     } 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>(
-              singularity::photons::MeanOpacityBase(opac, -1, 1, 2, -1, 1, 2), time_scale,
-              mass_scale, length_scale, temperature_scale);
+              singularity::photons::MeanOpacityBase(opac, -1, 1, 2, -1, 1, 2, gray_bounds,
+                                                    1),
+              time_scale, mass_scale, length_scale, temperature_scale);
     } else if (abs_model == "table") {
       std::string table_filename = pin->GetString("mcblock/absorption", "opacity_table");
       mopacity =
@@ -186,18 +192,21 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   std::string sct_model =
       pin->GetOrAddString("mcblock/scattering", "opacity_model", "none");
   if (frequency_type == FrequencyType::gray) {
+    const std::array<Real, 2> gray_bounds = {1.e12, 3.e20};
     if (sct_model == "none") {
       auto sopac = singularity::photons::GrayS(0.0, apm);
       mscattering =
-          singularity::photons::MeanNonCGSUnitsS<singularity::photons::MeanSOpacityCGS>(
-              singularity::photons::MeanSOpacityCGS(sopac, -1., 1., 2, -1., 1., 2),
+          singularity::photons::MeanNonCGSUnitsS<singularity::photons::MeanSOpacityBase>(
+              singularity::photons::MeanSOpacityBase(sopac, -1., 1., 2, -1., 1., 2,
+                                                     gray_bounds, 1),
               time_scale, mass_scale, length_scale, 1.);
     } 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<singularity::photons::MeanSOpacityCGS>(
-              singularity::photons::MeanSOpacityCGS(sopac, -1., 1., 2, -1., 1., 2),
+          singularity::photons::MeanNonCGSUnitsS<singularity::photons::MeanSOpacityBase>(
+              singularity::photons::MeanSOpacityBase(sopac, -1., 1., 2, -1., 1., 2,
+                                                     gray_bounds, 1),
               time_scale, mass_scale, length_scale, 1.);
     } else {
       PARTHENON_FAIL("Only none or constant scattering models supported!");
@@ -321,8 +330,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
           const Real &rho = vmesh(b, fm::density(), k, j, i);
           const Real &sie = vmesh(b, fm::sie(), k, j, i);
           const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie);
-          const Real aa = mopacity.AbsorptionCoefficient(rho, temp, gmode);
-          const Real ss = mscattering.RosselandMeanTotalScatteringCoefficient(rho, temp);
+          const Real aa = mopacity.AbsorptionCoefficient(rho, temp, 0, gmode);
+          const Real ss = mscattering.ScatteringCoefficient(rho, temp, 0, gmode);
           vmesh(b, fm::absorption_opacity(), k, j, i) = aa;
           vmesh(b, fm::scattering_opacity(), k, j, i) = ss;
         });
@@ -338,6 +347,10 @@ void UpdateDerived(MeshData<Real> *md) {
   PARTHENON_INSTRUMENT
   namespace fm = field::material;
   using parthenon::MakePackDescriptor;
+  using singularity::photons::OpacityAveraging;
+  using singularity::photons::Planck;
+  using singularity::photons::Rosseland;
+
   auto pm = md->GetParentPointer();
   auto &resolved_pkgs = pm->resolved_packages;
   auto &jbn = pm->packages.Get("jaybenne");
@@ -376,15 +389,22 @@ void UpdateDerived(MeshData<Real> *md) {
 
   // update opacity (TODO: only gray for now)
   if (frequency_type == FrequencyType::gray) {
+
+    // get opacity average indicators
+    const auto &use_planck = jbn->template Param<bool>("use_planck");
+    const auto &use_rosseland = jbn->template Param<bool>("use_rosseland");
+    // set opacity mode for grey opacity
+    const OpacityAveraging gmode = (use_planck && !use_rosseland) ? Planck : Rosseland;
+
     parthenon::par_for(
         DEFAULT_LOOP_PATTERN, "Update opacity", parthenon::DevExecSpace(), 0,
         vmesh.GetNBlocks() - 1, kbe.s, kbe.e, jbe.s, jbe.e, ibe.s, ibe.e,
         KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) {
           const Real &rho = vmesh(b, fm::density(), k, j, i);
           const Real &sie = vmesh(b, fm::sie(), k, j, i);
           const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie);
-          const Real aa = mopacity.AbsorptionCoefficient(rho, temp);
-          const Real ss = mscattering.RosselandMeanTotalScatteringCoefficient(rho, temp);
+          const Real aa = mopacity.AbsorptionCoefficient(rho, temp, 0, gmode);
+          const Real ss = mscattering.ScatteringCoefficient(rho, temp, 0, gmode);
           vmesh(b, fm::absorption_opacity(), k, j, i) = aa;
           vmesh(b, fm::scattering_opacity(), k, j, i) = ss;
         });

src/mcblock/opacity.hpp

@@ -37,7 +37,7 @@ using Scattering = singularity::photons::impl::S_Variant<
     singularity::photons::NonCGSUnitsS<singularity::photons::ThomsonS>>;
 
 using MeanScattering =
-    singularity::photons::MeanNonCGSUnitsS<singularity::photons::MeanSOpacityCGS>;
+    singularity::photons::MeanNonCGSUnitsS<singularity::photons::MeanSOpacityBase>;
 
 } // namespace mcblock