Add scattering

singularity-opac/base/sp5.hpp

@@ -33,6 +33,11 @@ constexpr char PlanckMeanOpacity[] = "Planck mean opacity";
 constexpr char RosselandMeanOpacity[] = "Rosseland mean opacity";
 } // namespace MeanOpac
 
+namespace MeanSOpac {
+constexpr char PlanckMeanSOpacity[] = "Planck mean scattering opacity";
+constexpr char RosselandMeanSOpacity[] = "Rosseland mean scattering opacity";
+} // namespace MeanSOpac
+
 } // namespace SP5
 
 #endif // SINGULARITY_OPAC_BASE_SP5_

singularity-opac/neutrinos/gray_s_opacity_neutrinos.hpp

@@ -0,0 +1,75 @@
+// ======================================================================
+// © 2022. 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 SINGULARITY_OPAC_NEUTRINOS_GRAY_S_OPACITY_NEUTRINOS_
+#define SINGULARITY_OPAC_NEUTRINOS_GRAY_S_OPACITY_NEUTRINOS_
+
+#include <cassert>
+#include <cmath>
+#include <cstdio>
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-opac/base/opac_error.hpp>
+
+namespace singularity {
+namespace neutrinos {
+
+template <typename pc = PhysicalConstantsCGS>
+class GraySOpacity {
+ public:
+  GraySOpacity() = default;
+  GraySOpacity(const Real sigma, const Real avg_particle_mass)
+      : sigma_(sigma), apm_(avg_particle_mass) {}
+
+  GraySOpacity GetOnDevice() { return *this; }
+  PORTABLE_INLINE_FUNCTION
+  int nlambda() const noexcept { return 0; }
+  PORTABLE_INLINE_FUNCTION
+  void PrintParams() const noexcept {
+    printf("Gray scattering opacity. sigma = %g avg particle mass = %g\n",
+           sigma_, apm_);
+  }
+  inline void Finalize() noexcept {}
+
+  PORTABLE_INLINE_FUNCTION
+  Real TotalCrossSection(const Real rho, const Real temp, const Real Ye,
+                         const RadiationType type, const Real nu,
+                         Real *lambda = nullptr) const {
+    return sigma_;
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real DifferentialCrossSection(const Real rho, const Real temp, const Real Ye,
+                                const RadiationType type, const Real nu,
+                                const Real mu, Real *lambda = nullptr) const {
+    return sigma_ / (4. * M_PI);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real TotalScatteringCoefficient(const Real rho, const Real temp,
+                                  const Real Ye, const RadiationType type,
+                                  const Real nu, Real *lambda = nullptr) const {
+    return (rho / apm_) * sigma_;
+  }
+
+ private:
+  Real sigma_; // Scattering cross section. Units of cm^2
+  Real apm_;   // Mean molecular weight. Units of g
+};
+
+} // namespace neutrinos
+} // namespace singularity
+
+#endif // SINGULARITY_OPAC_NEUTRINOS_GRAY_S_OPACITY_NEUTRINOS_

singularity-opac/neutrinos/mean_neutrino_s_variant.hpp

@@ -0,0 +1,104 @@
+// ======================================================================
+// © 2022. 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 SINGULARITY_OPAC_NEUTRINOS_MEAN_NEUTRINO_S_VARIANT_
+#define SINGULARITY_OPAC_NEUTRINOS_MEAN_NEUTRINO_S_VARIANT_
+
+#include <utility>
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-opac/base/opac_error.hpp>
+#include <singularity-opac/base/radiation_types.hpp>
+#include <singularity-opac/neutrinos/neutrino_s_variant.hpp>
+#include <variant/include/mpark/variant.hpp>
+
+namespace singularity {
+namespace neutrinos {
+namespace impl {
+
+template <typename... SOpacs>
+class MeanSVariant {
+ private:
+  s_opac_variant<SOpacs...> s_opac_;
+
+ public:
+  template <
+      typename Choice,
+      typename std::enable_if<
+          !std::is_same<MeanSVariant, typename std::decay<Choice>::type>::value,
+          bool>::type = true>
+  PORTABLE_FUNCTION MeanSVariant(Choice &&choice)
+      : s_opac_(std::forward<Choice>(choice)) {}
+
+  PORTABLE_FUNCTION
+  MeanSVariant() noexcept = default;
+
+  template <
+      typename Choice,
+      typename std::enable_if<
+          !std::is_same<MeanSVariant, typename std::decay<Choice>::type>::value,
+          bool>::type = true>
+  PORTABLE_FUNCTION MeanSVariant &operator=(Choice &&s_opac) {
+    s_opac_ = std::forward<Choice>(s_opac);
+    return *this;
+  }
+
+  template <
+      typename Choice,
+      typename std::enable_if<
+          !std::is_same<MeanSVariant, typename std::decay<Choice>::type>::value,
+          bool>::type = true>
+  Choice get() {
+    return mpark::get<Choice>(s_opac_);
+  }
+
+  MeanSVariant GetOnDevice() {
+    return mpark::visit(
+        [](auto &s_opac) {
+          return s_opac_variant<SOpacs...>(s_opac.GetOnDevice());
+        },
+        s_opac_);
+  }
+
+  PORTABLE_INLINE_FUNCTION Real PlanckMeanTotalScatteringCoefficient(
+      const Real rho, const Real temp, const Real Ye,
+      const RadiationType type) const {
+    return mpark::visit(
+        [=](const auto &s_opac) {
+          return s_opac.PlanckMeanTotalScatteringCoefficient(rho, temp, Ye, type);
+        },
+        s_opac_);
+  }
+  PORTABLE_INLINE_FUNCTION Real RosselandMeanTotalScatteringCoefficient(
+      const Real rho, const Real temp, const Real Ye,
+      const RadiationType type) const {
+    return mpark::visit(
+        [=](const auto &s_opac) {
+          return s_opac.RosselandMeanTotalScatteringCoefficient(rho, temp, Ye, type);
+        },
+        s_opac_);
+  }
+
+  inline void Finalize() noexcept {
+    return mpark::visit([](auto &s_opac) { return s_opac.Finalize(); },
+                        s_opac_);
+  }
+};
+
+} // namespace impl
+} // namespace neutrinos
+} // namespace singularity
+
+#endif // SINGULARITY_OPAC_NEUTRINOS_MEAN_NEUTRINO_S_VARIANT_