Power law opacity

singularity-opac/photons/powerlaw_opacity_photons.hpp

@@ -0,0 +1,181 @@
+// ======================================================================
+// © 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 A, const Real B)
+      : kappa0_(kappa0), A_(A), B_(B) {}
+  PowerLawOpacity(const PlanckDistribution<pc> &dist, const Real kappa0,
+                  const Real A, const Real B)
+      : dist_(dist), kappa0_(kappa0), A_(A), B_(B) {}
+
+  GrayOpacity 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 A = %g B = %g\n", kappa0_, A_, B_);
+  }
+  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, A_) * std::pow(temp, B_) * 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, A_) * std::pow(temp, B_)) * B;
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  Real NumberEmissivity(const Real rho, const Real temp,
+                        Real *lambda = nullptr) const {
+    return (kappa0_ * std::pow(rho, A_) * std::pow(temp, B_)) *
+           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 A_;      // Power law index of density
+  Real B_;      // 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