Power law opacity

CMakeLists.txt

@@ -53,8 +53,8 @@ target_include_directories(singularity-opac::flags
                            INTERFACE
                            $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}>)
 
-include (SetupDeps)
 include (SetupOptions)
+include (SetupDeps)
 include (SetupCompilers)
 include (SetupFlags)
 

README.md

@@ -57,6 +57,7 @@ A number of options are avaialable for compiling:
 | --------------------------------- | ------- | ------------------------------------------------------------------------------------ |
 | SINGULARITY_BUILD_TESTS           | OFF     | Build test infrastructure.                                                           |
 | SINGULARITY_USE_HDF5              | ON      | Enables HDF5. Required for Spiner opacities.                                         |
+| SINGULARITY_KOKKOS_IN_TREE        | OFF     | Force cmake to use Kokkos source included in tree.                                   |
 
 ## Copyright
 

cmake/SetupDeps.cmake

@@ -10,23 +10,32 @@ else()
   message(status "CUDA::toolkit provided by parent package")
 endif()
 
+#=======================================
+# Setup Kokkos
+# - provides Kokkos::kokkos
+#=======================================
+if (SINGULARITY_USE_KOKKOS)
+  if (NOT TARGET Kokkos::kokkos)
+    message(status "Kokkos::kokkos must be found")
+    if (SINGULARITY_KOKKOS_IN_TREE)
+      message(status "Using in-tree Kokkos")
+      add_subdirectory(utils/kokkos)
+    else()
+      message(status "Using system Kokkos if available")
+      find_package(Kokkos REQUIRED)
+    endif()
+  else()
+    message(status "Kokkos::kokkos provided by parent package")
+  endif()
+endif()
+
 #=======================================
 # Setup ports of call
 # - provides PortsofCall::PortsofCall
 #=======================================
 find_package(PortsofCall REQUIRED)
 target_link_libraries(singularity-opac::flags INTERFACE PortsofCall::PortsofCall)
 
-#=======================================
-# Setup Kokkos
-# - provides Kokkos::kokkos
-#=======================================
-if (NOT TARGET Kokkos::kokkos)
-  find_package(Kokkos QUIET)
-else()
-  message(status "Kokkos::kokkos provided by parent package")
-endif()
-
 #=======================================
 # Find HDF5
 # - [email protected]+ provides HDF5::HDF5, but

cmake/SetupOptions.cmake

@@ -17,11 +17,18 @@ option (SINGULARITY_USE_FMATH "Enable fast-math logarithms" ON)
 # Dependency options
 #=======================================
 # check for using in-tree third-party dependencies
-option (SINULARITY_KOKKOS_IN_TREE "Use in-tree dependencies" OFF)
-
 option (SINGULARITY_USE_KOKKOS "Use Kokkos for portability" OFF)
 option (SINGULARITY_USE_CUDA "Enable cuda support" OFF)
 option (SINGULARITY_USE_HDF5 "Pull in hdf5" OFF)
+cmake_dependent_option(SINULARITY_KOKKOS_IN_TREE
+  "Use in-tree dependencies" OFF
+  ${SINGULARITY_USE_KOKKOS} OFF)
+
+# Kill cmake's package registry because it can interfere
+if (SINGULARITY_KOKKOS_IN_TREE)
+  set(CMAKE_FIND_PACKAGE_NO_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE)
+  set(CMAKE_FIND_PACKAGE_NO_SYSTEM_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE)
+endif()
 
 # If the conditional is TRUE, it's the first default, else it's the
 # second.

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
@@ -58,7 +59,7 @@ PRIVATE
 # Ensure code works with C++11 and earlier
 # TODO(MM): Remove this later when it's not needed.
 set_target_properties(${PROJECT_NAME}_unit_tests
-  PROPERTIES CXX_STANDARD 14
+  PROPERTIES CXX_STANDARD 17
   CXX_STANDARD_REQUIRED YES
   CXX_EXTENSIONS NO)