diff --git a/CMakeLists.txt b/CMakeLists.txt index e34cb89..b34b014 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -53,8 +53,8 @@ target_include_directories(singularity-opac::flags INTERFACE $) -include (SetupDeps) include (SetupOptions) +include (SetupDeps) include (SetupCompilers) include (SetupFlags) diff --git a/README.md b/README.md index 6a34156..459588f 100644 --- a/README.md +++ b/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 diff --git a/cmake/SetupDeps.cmake b/cmake/SetupDeps.cmake index 16455f5..99e1ba0 100644 --- a/cmake/SetupDeps.cmake +++ b/cmake/SetupDeps.cmake @@ -10,6 +10,25 @@ 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 @@ -17,16 +36,6 @@ endif() 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 # - cmake@3.20+ provides HDF5::HDF5, but diff --git a/cmake/SetupOptions.cmake b/cmake/SetupOptions.cmake index 21a857c..9c9e53d 100644 --- a/cmake/SetupOptions.cmake +++ b/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. diff --git a/singularity-opac/photons/opac_photons.hpp b/singularity-opac/photons/opac_photons.hpp index 2128a2c..50918a5 100644 --- a/singularity-opac/photons/opac_photons.hpp +++ b/singularity-opac/photons/opac_photons.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -29,10 +30,13 @@ namespace photons { using ScaleFree = GrayOpacity; using Gray = GrayOpacity; +using PowerLawScaleFree = PowerLawOpacity; +using PowerLaw = PowerLawOpacity; using EPBremss = EPBremsstrahlungOpacity; -using Opacity = impl::Variant, - NonCGSUnits>; +using Opacity = impl::Variant, + NonCGSUnits, NonCGSUnits>; } // namespace photons } // namespace singularity diff --git a/singularity-opac/photons/powerlaw_opacity_photons.hpp b/singularity-opac/photons/powerlaw_opacity_photons.hpp new file mode 100644 index 0000000..7efe5e6 --- /dev/null +++ b/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 +#include +#include + +#include +#include +#include + +namespace singularity { +namespace photons { + +template +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 &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 + 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 + 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 + 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 + 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 dist_; +}; + +} // namespace photons +} // namespace singularity + +#endif // SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 06f6a86..21be197 100644 --- a/test/CMakeLists.txt +++ b/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) diff --git a/test/test_powerlaw_opacities.cpp b/test/test_powerlaw_opacities.cpp new file mode 100644 index 0000000..73e39b1 --- /dev/null +++ b/test/test_powerlaw_opacities.cpp @@ -0,0 +1,198 @@ +// ====================================================================== +// © 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. +// ====================================================================== + +#include +#include + +#include + +#include +#include +#include + +#include + +using namespace singularity; + +using DataBox = Spiner::DataBox; + +#ifdef PORTABILITY_STRATEGY_KOKKOS +using atomic_view = Kokkos::MemoryTraits; +#endif + +template +PORTABLE_INLINE_FUNCTION T FractionalDifference(const T &a, const T &b) { + return 2 * std::abs(b - a) / (std::abs(a) + std::abs(b) + 1e-20); +} +PORTABLE_INLINE_FUNCTION Real CalcFrequency(const int n, const Real nu_min, + const Real nu_max, const int n_nu) { + const Real lnu_min = std::log(nu_min); + const Real lnu_max = std::log(nu_max); + const Real dlnu = (lnu_max - lnu_min) / static_cast(n_nu); + return std::exp(lnu_min + (n + 0.5) * dlnu); +} + +constexpr Real EPS_TEST = 1e-3; + +constexpr Real rho_exp = 2.; +constexpr Real temp_exp = 2.5; + +TEST_CASE("Scale free power law photon opacities", + "[PowerLawScaleFreePhotonOpacities]") { + WHEN("We initialize a scale-free power law photon opacity") { + constexpr Real rho = 1.5e0; + constexpr Real temp = 3.2e0; + constexpr Real nu_min = 1.e-1; + constexpr Real nu_max = 1.e1; + constexpr int n_nu = 100; + constexpr Real kappa0 = 1.5; + + photons::PowerLawScaleFree opac_host(kappa0, rho_exp, temp_exp); + photons::Opacity opac = opac_host.GetOnDevice(); + + THEN("The emissivity per nu omega is consistent with the emissity per nu") { + int n_wrong_h = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View n_wrong_d("wrong"); +#else + PortableMDArray n_wrong_d(&n_wrong_h, 1); +#endif + + portableFor( + "calc emissivities", 0, n_nu, PORTABLE_LAMBDA(const int &i) { + const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu); + const Real jnu = opac.EmissivityPerNuOmega(rho, temp, nu); + const Real Jnu = opac.EmissivityPerNu(rho, temp, nu); + const Real kappa = + kappa0 * std::pow(rho, rho_exp) * std::pow(temp, temp_exp); + if (FractionalDifference(jnu, rho * kappa * + opac.ThermalDistributionOfTNu( + temp, nu)) > EPS_TEST) { + n_wrong_d() += 1; + } + if (FractionalDifference(Jnu, 4 * M_PI * jnu) > EPS_TEST) { + n_wrong_d() += 1; + } + }); + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(n_wrong_h, n_wrong_d); +#endif + REQUIRE(n_wrong_h == 0); + } + + opac.Finalize(); + } +} + +TEST_CASE("CGS power law photon opacities", "[PowerLawCGSPhotonOpacities]") { + constexpr Real rho = 1.5e0; // g/cc + constexpr Real temp = 1.e3; // K + constexpr Real nu_min = 1.e10; // Hz + constexpr Real nu_max = 1.e14; // Hz + constexpr int n_nu = 100; + constexpr Real kappa0 = 1.5; // cm^2 / g + + WHEN("We initialize a CGS power law photon opacity") { + photons::PowerLaw opac_host(kappa0, rho_exp, temp_exp); + photons::Opacity opac = opac_host.GetOnDevice(); + + THEN("The emissivity per nu omega is consistent with the emissity per nu") { + int n_wrong_h = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View n_wrong_d("wrong"); +#else + PortableMDArray n_wrong_d(&n_wrong_h, 1); +#endif + + portableFor( + "calc emissivities", 0, n_nu, PORTABLE_LAMBDA(const int &i) { + const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu); + const Real jnu = opac.EmissivityPerNuOmega(rho, temp, nu); + const Real Jnu = opac.EmissivityPerNu(rho, temp, nu); + const Real kappa = + kappa0 * std::pow(rho, rho_exp) * std::pow(temp, temp_exp); + if (FractionalDifference(jnu, rho * kappa * + opac.ThermalDistributionOfTNu( + temp, nu)) > EPS_TEST) { + n_wrong_d() += 1; + } + if (FractionalDifference(Jnu, 4 * M_PI * jnu) > EPS_TEST) { + n_wrong_d() += 1; + } + }); + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(n_wrong_h, n_wrong_d); +#endif + REQUIRE(n_wrong_h == 0); + } + + opac.Finalize(); + } + + WHEN("We initialize a CGS power law photon opacity with non-CGS units") { + constexpr Real time_unit = 123.; + constexpr Real mass_unit = 456.; + constexpr Real length_unit = 789.; + constexpr Real temp_unit = 276.; + constexpr Real rho_unit = + mass_unit / (length_unit * length_unit * length_unit); + constexpr Real j_unit = mass_unit / (length_unit * time_unit * time_unit); + + photons::NonCGSUnits opac_host( + photons::PowerLaw(kappa0, rho_exp, temp_exp), time_unit, mass_unit, + length_unit, temp_unit); + photons::Opacity opac = opac_host.GetOnDevice(); + photons::PowerLaw opac_cgs_host(kappa0, rho_exp, temp_exp); + photons::Opacity opac_cgs = opac_cgs_host.GetOnDevice(); + + THEN("The emissivity per nu omega is consistent with the emissity per nu") { + int n_wrong_h = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View n_wrong_d("wrong"); +#else + PortableMDArray n_wrong_d(&n_wrong_h, 1); +#endif + + portableFor( + "calc emissivities", 0, n_nu, PORTABLE_LAMBDA(const int &i) { + const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu); + const Real jnu = opac.EmissivityPerNuOmega( + rho / rho_unit, temp / temp_unit, nu * time_unit); + const Real Jnu = opac.EmissivityPerNu( + rho / rho_unit, temp / temp_unit, nu * time_unit); + const Real kappa = + kappa0 * std::pow(rho, rho_exp) * std::pow(temp, temp_exp); + if (FractionalDifference( + jnu * j_unit, + opac_cgs.EmissivityPerNuOmega(rho, temp, nu)) > EPS_TEST) { + n_wrong_d() += 1; + } + if (FractionalDifference(Jnu, 4 * M_PI * jnu) > EPS_TEST) { + n_wrong_d() += 1; + } + }); + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(n_wrong_h, n_wrong_d); +#endif + REQUIRE(n_wrong_h == 0); + } + + opac.Finalize(); + } +} + diff --git a/utils/kokkos b/utils/kokkos index 1fb0c28..08ceff9 160000 --- a/utils/kokkos +++ b/utils/kokkos @@ -1 +1 @@ -Subproject commit 1fb0c284d458c75370094921d9f202c287502325 +Subproject commit 08ceff92bcf3a828844480bc1e6137eb74028517 diff --git a/utils/ports-of-call b/utils/ports-of-call index 993d820..58ce118 160000 --- a/utils/ports-of-call +++ b/utils/ports-of-call @@ -1 +1 @@ -Subproject commit 993d8209de98da186c5e0b445a4d6c359afa8c81 +Subproject commit 58ce1181b2d835bd32673ad70550c9130381f91b diff --git a/utils/spiner b/utils/spiner index 33862f8..bd57161 160000 --- a/utils/spiner +++ b/utils/spiner @@ -1 +1 @@ -Subproject commit 33862f831be65e2dd2284b42f05c477b6fe7367a +Subproject commit bd57161576a62a13341dada183f2b336a3e99b08 diff --git a/utils/variant b/utils/variant index 3c7fc82..23cb94f 160000 --- a/utils/variant +++ b/utils/variant @@ -1 +1 @@ -Subproject commit 3c7fc8266bb46046b42c2dc2663f9f505f0cec28 +Subproject commit 23cb94f027d4ef33bf48133acc2695c7e5c6f1e7