Power law opacity

CMakeLists.txt

@@ -123,3 +123,4 @@ feature_summary(WHAT ALL
                 DESCRIPTION "Enabled Features:"
                 VAR enabledFeaturesText)
 message(STATUS "${enabledFeaturesText}")
+message(STATUS "NOW 127 Compile Definitions for ${PROJECT_NAME}: ${defs}")

cmake/SetupDeps.cmake

@@ -22,6 +22,7 @@ target_link_libraries(singularity-opac::flags INTERFACE PortsofCall::PortsofCall
 # - provides Kokkos::kokkos
 #=======================================
 if (NOT TARGET Kokkos::kokkos)
+  add_subdirectory(utils/kokkos)
   find_package(Kokkos QUIET)
 else()
   message(status "Kokkos::kokkos provided by parent package")
@@ -90,3 +91,4 @@ cmake_dependent_option(SINGULARITY_USE_MPI
 if (NOT TARGET Catch2::Catch2)
   find_package(Catch2 QUIET)
 endif()
+

cmake/SetupFlags.cmake

@@ -72,6 +72,7 @@ INTERFACE
     $<${with_fmath}:
       SINGULARITY_USE_FMATH
     >
+    PORTABILITY_STRATEGY_NONE
   >
 )
 

cmake/SetupOptions.cmake

@@ -51,3 +51,9 @@ if(SINGULARITY_USE_CUDA)
     message(FATAL_ERROR "Cuda without kokkos is not currently supported")
   endif()
 endif()
+
+if(NOT SINGULARITY_USE_KOKKOS)
+  message(FATAL_ERROR
+    "For unknown reasons PORTABILITY_STRATEGY_KOKKOS will always be provided to ports-of-call. "
+    "For now you must use Kokkos with singularity-opac.")
+endif()

singularity-opac/photons/opac_photons.hpp

@@ -20,6 +20,7 @@
 #include <singularity-opac/photons/gray_opacity_photons.hpp>
 #include <singularity-opac/photons/non_cgs_photons.hpp>
 #include <singularity-opac/photons/photon_variant.hpp>
+#include <singularity-opac/photons/powerlaw_opacity_photons.hpp>
 #include <singularity-opac/photons/thermal_distributions_photons.hpp>
 
 #include <singularity-opac/photons/mean_opacity_photons.hpp>
@@ -29,10 +30,13 @@ namespace photons {
 
 using ScaleFree = GrayOpacity<PhysicalConstantsUnity>;
 using Gray = GrayOpacity<PhysicalConstantsCGS>;
+using PowerLawScaleFree = PowerLawOpacity<PhysicalConstantsUnity>;
+using PowerLaw = PowerLawOpacity<PhysicalConstantsCGS>;
 using EPBremss = EPBremsstrahlungOpacity<PhysicalConstantsCGS>;
 
-using Opacity = impl::Variant<ScaleFree, Gray, EPBremss, NonCGSUnits<Gray>,
-                              NonCGSUnits<EPBremss>>;
+using Opacity = impl::Variant<ScaleFree, Gray, PowerLawScaleFree, PowerLaw,
+                              EPBremss, NonCGSUnits<Gray>,
+                              NonCGSUnits<PowerLaw>, NonCGSUnits<EPBremss>>;
 
 } // namespace photons
 } // namespace singularity

singularity-opac/photons/powerlaw_opacity_photons.hpp

@@ -0,0 +1,185 @@
+// ======================================================================
+// © 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.
+// ======================================================================
+
+#ifndef SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_
+#define SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_
+
+#include <cassert>
+#include <cmath>
+#include <cstdio>
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-opac/base/opac_error.hpp>
+#include <singularity-opac/photons/thermal_distributions_photons.hpp>
+
+namespace singularity {
+namespace photons {
+
+template <typename pc = PhysicalConstantsCGS>
+class PowerLawOpacity {
+ public:
+  PowerLawOpacity() = default;
+  PowerLawOpacity(const Real kappa0, const Real rho_exp, const Real temp_exp)
+      : kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp) {}
+  PowerLawOpacity(const PlanckDistribution<pc> &dist, const Real kappa0,
+                  const Real rho_exp, const Real temp_exp)
+      : dist_(dist), kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp) {}
+
+  PowerLawOpacity GetOnDevice() { return *this; }
+  PORTABLE_INLINE_FUNCTION
+  int nlambda() const noexcept { return 0; }
+  PORTABLE_INLINE_FUNCTION
+  void PrintParams() const noexcept {
+    printf("Power law opacity. kappa0 = %g rho_exp = %g temp_exp = %g\n",
+           kappa0_, rho_exp_, temp_exp_);
+  }
+  inline void Finalize() noexcept {}
+
+  PORTABLE_INLINE_FUNCTION
+  Real AbsorptionCoefficient(const Real rho, const Real temp, const Real nu,
+                             Real *lambda = nullptr) const {
+    return dist_.AbsorptionCoefficientFromKirkhoff(*this, rho, temp, nu,
+                                                   lambda);
+  }
+
+  template <typename FrequencyIndexer, typename DataIndexer>
+  PORTABLE_INLINE_FUNCTION void
+  AbsorptionCoefficient(const Real rho, const Real temp,
+                        FrequencyIndexer &nu_bins, DataIndexer &coeffs,
+                        const int nbins, Real *lambda = nullptr) const {
+    for (int i = 0; i < nbins; ++i) {
+      coeffs[i] = AbsorptionCoefficient(rho, temp, nu_bins[i], lambda);
+    }
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real AngleAveragedAbsorptionCoefficient(const Real rho, const Real temp,
+                                          const Real nu,
+                                          Real *lambda = nullptr) const {
+    return dist_.AngleAveragedAbsorptionCoefficientFromKirkhoff(
+        *this, rho, temp, nu, lambda);
+  }
+
+  template <typename FrequencyIndexer, typename DataIndexer>
+  PORTABLE_INLINE_FUNCTION void AngleAveragedAbsorptionCoefficient(
+      const Real rho, const Real temp, FrequencyIndexer &nu_bins,
+      DataIndexer &coeffs, const int nbins, Real *lambda = nullptr) const {
+    for (int i = 0; i < nbins; ++i) {
+      coeffs[i] =
+          AngleAveragedAbsorptionCoefficient(rho, temp, nu_bins[i], lambda);
+    }
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real EmissivityPerNuOmega(const Real rho, const Real temp, const Real nu,
+                            Real *lambda = nullptr) const {
+    Real Bnu = dist_.ThermalDistributionOfTNu(temp, nu, lambda);
+    return rho *
+           (kappa0_ * std::pow(rho, rho_exp_) * std::pow(temp, temp_exp_)) *
+           Bnu;
+  }
+
+  template <typename FrequencyIndexer, typename DataIndexer>
+  PORTABLE_INLINE_FUNCTION void
+  EmissivityPerNuOmega(const Real rho, const Real temp,
+                       FrequencyIndexer &nu_bins, DataIndexer &coeffs,
+                       const int nbins, Real *lambda = nullptr) const {
+    for (int i = 0; i < nbins; ++i) {
+      coeffs[i] = EmissivityPerNuOmega(rho, temp, nu_bins[i], lambda);
+    }
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real EmissivityPerNu(const Real rho, const Real temp, const Real nu,
+                       Real *lambda = nullptr) const {
+    return 4 * M_PI * EmissivityPerNuOmega(rho, temp, nu, lambda);
+  }
+
+  template <typename FrequencyIndexer, typename DataIndexer>
+  PORTABLE_INLINE_FUNCTION void
+  EmissivityPerNu(const Real rho, const Real temp, FrequencyIndexer &nu_bins,
+                  DataIndexer &coeffs, const int nbins,
+                  Real *lambda = nullptr) const {
+    for (int i = 0; i < nbins; ++i) {
+      coeffs[i] = EmissivityPerNu(rho, temp, nu_bins[i], lambda);
+    }
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real Emissivity(const Real rho, const Real temp,
+                  Real *lambda = nullptr) const {
+    Real B = dist_.ThermalDistributionOfT(temp, lambda);
+    return rho *
+           (kappa0_ * std::pow(rho, rho_exp_) * std::pow(temp, temp_exp_)) * B;
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real NumberEmissivity(const Real rho, const Real temp,
+                        Real *lambda = nullptr) const {
+    return (kappa0_ * std::pow(rho, rho_exp_) * std::pow(temp, temp_exp_)) *
+           dist_.ThermalNumberDistributionOfT(temp, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real ThermalDistributionOfTNu(const Real temp, const Real nu,
+                                Real *lambda = nullptr) const {
+    return dist_.ThermalDistributionOfTNu(temp, nu, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real DThermalDistributionOfTNuDT(const Real temp, const Real nu,
+                                   Real *lambda = nullptr) const {
+    return dist_.DThermalDistributionOfTNuDT(temp, nu, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real ThermalDistributionOfT(const Real temp, Real *lambda = nullptr) const {
+    return dist_.ThermalDistributionOfT(temp, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION Real
+  ThermalNumberDistributionOfT(const Real temp, Real *lambda = nullptr) const {
+    return dist_.ThermalNumberDistributionOfT(temp, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real EnergyDensityFromTemperature(const Real temp,
+                                    Real *lambda = nullptr) const {
+    return dist_.EnergyDensityFromTemperature(temp, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real TemperatureFromEnergyDensity(const Real er,
+                                    Real *lambda = nullptr) const {
+    return dist_.TemperatureFromEnergyDensity(er, lambda);
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real NumberDensityFromTemperature(const Real temp,
+                                    Real *lambda = nullptr) const {
+    return dist_.NumberDensityFromTemperature(temp, lambda);
+  }
+
+ private:
+  Real kappa0_;   // Opacity scale. Units of cm^2/g
+  Real rho_exp_;  // Power law index of density
+  Real temp_exp_; // Power law index of temperature
+  PlanckDistribution<pc> dist_;
+};
+
+} // namespace photons
+} // namespace singularity
+
+#endif // SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_

test/CMakeLists.txt

@@ -39,9 +39,10 @@ PRIVATE
   test_thermal_dist.cpp
   test_scalefree_opacities.cpp
   test_gray_opacities.cpp
-  test_epbremsstrahlung_opacities.cpp
   test_gray_s_opacities.cpp
+  test_epbremsstrahlung_opacities.cpp
   test_brt_opacities.cpp
+  test_powerlaw_opacities.cpp
   test_thomson_s_opacities.cpp
   test_chebyshev.cpp
   test_spiner_opac_neutrinos.cpp