diff --git a/plan_histories/70.md b/plan_histories/70.md new file mode 100644 index 0000000..e760790 --- /dev/null +++ b/plan_histories/70.md @@ -0,0 +1,25 @@ +## Goal + +Extend photon `PowerLaw` opacities to support an optional frequency power law +while preserving the legacy gray behavior when no frequency exponent is +provided. + +## Functional Plan + +- Add optional `nu_exp` and `nu_ref` parameters to photon `PowerLaw`. +- Apply the frequency power law only to the monochromatic opacity and + monochromatic emissivity APIs. +- Preserve the existing legacy constructor calls by keeping the default + frequency behavior equivalent to `nu_exp = 0`. +- Keep the `NonCGSUnits` wrapper path working for the frequency-dependent + `PowerLaw`. +- Add focused regression coverage that checks: + - the monochromatic absorption coefficient matches the expected + `(\nu / \nu_ref)^{\nu_exp}` scaling + - the monochromatic emissivity remains `j_\nu = \alpha_\nu B_\nu` + - the same frequency dependence survives the `NonCGSUnits` wrapper + +## Scope Boundaries + +- Do not add integrated `Emissivity` or `NumberEmissivity` support for + `nu_exp != 0`. diff --git a/singularity-opac/photons/powerlaw_opacity_photons.hpp b/singularity-opac/photons/powerlaw_opacity_photons.hpp index 4a55177..69f5860 100644 --- a/singularity-opac/photons/powerlaw_opacity_photons.hpp +++ b/singularity-opac/photons/powerlaw_opacity_photons.hpp @@ -1,5 +1,5 @@ // ====================================================================== -// © 2024. Triad National Security, LLC. All rights reserved. This +// © 2024-2026. 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. @@ -16,6 +16,8 @@ #ifndef SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_ #define SINGULARITY_OPAC_PHOTONS_POWERLAW_OPACITY_PHOTONS_ +// This file was made in part with generative AI. + #include #include #include @@ -33,26 +35,61 @@ class PowerLawOpacity { using PC = pc; 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 Real kappa0, const Real rho_exp, const Real temp_exp, + const Real nu_exp = 0., const Real nu_ref = 1., const Real nu_off = 0., + const Real rho_ref = 1., const Real rho_off = 0., + const Real temp_ref = 1., const Real temp_off = 0., + const bool do_stim_emit = false) + : PowerLawOpacity(PlanckDistribution{}, kappa0, rho_exp, temp_exp, + nu_exp, nu_ref, nu_off, rho_ref, rho_off, temp_ref, + temp_off, do_stim_emit) {} 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) {} + const Real rho_exp, const Real temp_exp, + const Real nu_exp = 0., const Real nu_ref = 1., const Real nu_off = 0., + const Real rho_ref = 1., const Real rho_off = 0., + const Real temp_ref = 1., const Real temp_off = 0., + const bool do_stim_emit = false) + : dist_(dist), kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp), + rho_ref_(rho_ref), rho_off_(rho_off), temp_ref_(temp_ref), temp_off_(temp_off), + nu_exp_(nu_exp), nu_ref_(nu_ref), nu_off_(nu_off), do_stim_emit_(do_stim_emit) { + if (!(rho_ref_ > 0.)) { + OPAC_ERROR("PowerLawOpacity: rho_ref must be positive"); + } + if (!(temp_ref_ > 0.)) { + OPAC_ERROR("PowerLawOpacity: temp_ref must be positive"); + } + if (!(nu_ref_ > 0.)) { + OPAC_ERROR("PowerLawOpacity: nu_ref must be positive"); + } + if (rho_off_ < 0.) { + OPAC_ERROR("PowerLawOpacity: rho_off must be nonnegative"); + } + if (temp_off_ < 0.) { + OPAC_ERROR("PowerLawOpacity: temp_off must be nonnegative"); + } + if (nu_off_ < 0.) { + OPAC_ERROR("PowerLawOpacity: nu_off must be nonnegative"); + } + } 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_); + printf("Power law opacity. kappa0 = %g rho_exp = %g temp_exp = %g " + "nu_exp = %g nu_ref = %g nu_off = %g rho_ref = %g rho_off = %g" + "temp_ref = %g temp_off = %g do_stim_emit = %d\n", + kappa0_, rho_exp_, temp_exp_, nu_exp_, nu_ref_, + nu_off_, rho_ref_, rho_off_, temp_ref_, temp_off_, + do_stim_emit_); } inline void Finalize() noexcept {} PORTABLE_INLINE_FUNCTION Real AbsorptionCoefficient(const Real rho, const Real temp, const Real nu, Real *lambda = nullptr) const { - return rho * (kappa0_ * std::pow(rho, rho_exp_) * std::pow(temp, temp_exp_)); + return rho * OpacityScale_(rho, temp, nu); } template @@ -86,9 +123,7 @@ class PowerLawOpacity { 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; + return rho * OpacityScale_(rho, temp, nu) * Bnu; } template @@ -120,15 +155,34 @@ class PowerLawOpacity { 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; + // Once the opacity depends on frequency, the total emissivity is no longer + // the gray factorization alpha(T, rho) * ThermalDistributionOfT(T). This + // class intentionally supports only the monochromatic frequency-resolved + // emissivity APIs in that case. Supporting the fully integrated form would + // require carrying factorial/Gamma-style moments and Riemann zeta + // functions for the frequency-dependent power law. + if (nu_exp_ != 0.) { + OPAC_ERROR("PowerLawOpacity: total emissivity is only supported for " + "nu_exp = 0. Use EmissivityPerNuOmega or EmissivityPerNu " + "for frequency-dependent power laws."); + } + return rho * OpacityPrefactor_(rho, temp) * + dist_.ThermalDistributionOfT(temp, lambda); } 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_)) * + // Same limitation as Emissivity(): for frequency-dependent power laws, the + // integrated number emissivity is not treated as a gray closed-form API, + // and full support would likewise require factorial/Gamma-style moments + // and Riemann zeta functions. + if (nu_exp_ != 0.) { + OPAC_ERROR("PowerLawOpacity: total number emissivity is only supported " + "for nu_exp = 0. Use EmissivityPerNuOmega or EmissivityPerNu " + "for frequency-dependent power laws."); + } + return OpacityPrefactor_(rho, temp) * dist_.ThermalNumberDistributionOfT(temp, lambda); } @@ -178,9 +232,31 @@ class PowerLawOpacity { } 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 + PORTABLE_INLINE_FUNCTION + Real OpacityScale_(const Real rho, const Real temp, const Real nu) const { + const Real freq_plaw = std::pow((nu + nu_off_) / nu_ref_, nu_exp_); + const Real stim_fact = do_stim_emit_ ? -std::expm1(-(pc::h * nu / (pc::kb * temp))) : 1.0; + return OpacityPrefactor_(rho, temp) * freq_plaw * stim_fact; + } + + PORTABLE_INLINE_FUNCTION + Real OpacityPrefactor_(const Real rho, const Real temp) const { + const Real rhom = (rho + rho_off_) / rho_ref_; + const Real tempm = (temp + temp_off_) / temp_ref_; + return kappa0_ * std::pow(rhom, rho_exp_) * std::pow(tempm, temp_exp_); + } + + Real kappa0_; // Opacity scale. Units depend on nu_exp and nu_ref. + Real rho_exp_; // Power law index of density + Real temp_exp_; // Power law index of temperature + Real rho_ref_; // Density normalization for rho_exp. Units of g/cm^3 + Real rho_off_; // Density offset (same units as rho_ref). Units of g/cm^3 + Real temp_ref_; // Temperature normalization for temp_exp. Units of K + Real temp_off_; // Temperature offset (same units as temp_ref). Units of K + Real nu_exp_; // Power law index of frequency + Real nu_ref_; // Frequency normalization for nu_exp. Units of 1/s + Real nu_off_; // Frequency offset (same units as nu_ref). Units of 1/s + bool do_stim_emit_; // indicator to use stimulated (LTE) emission factor PlanckDistribution dist_; }; diff --git a/singularity-opac/photons/powerlaw_s_opacity_photons.hpp b/singularity-opac/photons/powerlaw_s_opacity_photons.hpp new file mode 100644 index 0000000..6185ea4 --- /dev/null +++ b/singularity-opac/photons/powerlaw_s_opacity_photons.hpp @@ -0,0 +1,125 @@ +// ====================================================================== +// © 2026. 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_S_OPACITY_PHOTONS_ +#define SINGULARITY_OPAC_PHOTONS_POWERLAW_S_OPACITY_PHOTONS_ + +// This file was partly copied from a file made in part with generative AI. + +#include +#include +#include + +#include +#include + +namespace singularity { +namespace photons { + +template +class PowerLawSOpacity { + public: + PowerLawSOpacity() = default; + PowerLawSOpacity(const Real kappa0, const Real rho_exp, const Real temp_exp, + const Real nu_exp = 0., const Real nu_ref = 1., const Real nu_off = 0., + const Real rho_ref = 1., const Real rho_off = 0., + const Real temp_ref = 1., const Real temp_off = 0., + const Real avg_particle_mass = 1.) + : kappa0_(kappa0), rho_exp_(rho_exp), temp_exp_(temp_exp), + rho_ref_(rho_ref), rho_off_(rho_off), temp_ref_(temp_ref), temp_off_(temp_off), + nu_exp_(nu_exp), nu_ref_(nu_ref), nu_off_(nu_off), apm_(avg_particle_mass) { + if (!(nu_ref_ > 0.)) { + OPAC_ERROR("PowerLawSOpacity: nu_ref must be positive"); + } + if (!(temp_ref_ > 0.)) { + OPAC_ERROR("PowerLawSOpacity: temp_ref must be positive"); + } + if (!(nu_ref_ > 0.)) { + OPAC_ERROR("PowerLawSOpacity: nu_ref must be positive"); + } + if (rho_off_ < 0.) { + OPAC_ERROR("PowerLawSOpacity: rho_off must be nonnegative"); + } + if (temp_off_ < 0.) { + OPAC_ERROR("PowerLawSOpacity: temp_off must be nonnegative"); + } + if (nu_off_ < 0.) { + OPAC_ERROR("PowerLawSOpacity: nu_off must be nonnegative"); + } + } + + PowerLawSOpacity GetOnDevice() { return *this; } + PORTABLE_INLINE_FUNCTION + int nlambda() const noexcept { return 0; } + PORTABLE_INLINE_FUNCTION + void PrintParams() const noexcept { + printf("Power law scattering opacity. kappa0 = %g rho_exp = %g temp_exp = %g " + "nu_exp = %g nu_ref = %g nu_off = %g rho_ref = %g rho_off = %g" + "temp_ref = %g temp_off = %g avg particle mass = %g\n", + kappa0_, rho_exp_, temp_exp_, nu_exp_, nu_ref_, + nu_off_, rho_ref_, rho_off_, temp_ref_, temp_off_, + apm_); + } + inline void Finalize() noexcept {} + + PORTABLE_INLINE_FUNCTION + Real TotalCrossSection(const Real rho, const Real temp, const Real nu, + Real *lambda = nullptr) const { + return OpacityScale_(rho, temp, nu); + } + + PORTABLE_INLINE_FUNCTION + Real DifferentialCrossSection(const Real rho, const Real temp, const Real nu, + const Real mu, Real *lambda = nullptr) const { + // assumed isotropic, elastic + return OpacityScale_(rho, temp, nu) / (4. * M_PI); + } + + PORTABLE_INLINE_FUNCTION + Real TotalScatteringCoefficient(const Real rho, const Real temp, + const Real nu, Real *lambda = nullptr) const { + return (rho / apm_) * OpacityScale_(rho, temp, nu); + } + + private: + PORTABLE_INLINE_FUNCTION + Real OpacityScale_(const Real rho, const Real temp, const Real nu) const { + return OpacityPrefactor_(rho, temp) * std::pow((nu + nu_off_) / nu_ref_, nu_exp_); + } + + PORTABLE_INLINE_FUNCTION + Real OpacityPrefactor_(const Real rho, const Real temp) const { + const Real rhom = (rho + rho_off_) / rho_ref_; + const Real tempm = (temp + temp_off_) / temp_ref_; + return kappa0_ * std::pow(rhom, rho_exp_) * std::pow(tempm, temp_exp_); + } + + Real kappa0_; // Opacity scale. Units depend on nu_exp and nu_ref. + Real rho_exp_; // Power law index of density + Real temp_exp_; // Power law index of temperature + Real rho_ref_; // Density normalization for rho_exp. Units of g/cm^3 + Real rho_off_; // Density offset (same units as rho_ref). Units of g/cm^3 + Real temp_ref_; // Temperature normalization for temp_exp. Units of K + Real temp_off_; // Temperature offset (same units as temp_ref). Units of K + Real nu_exp_; // Power law index of frequency + Real nu_ref_; // Frequency normalization for nu_exp. Units of 1/s + Real nu_off_; // Frequency offset (same units as nu_ref). Units of 1/s + Real apm_; // Mean molecular weight. Units of g +}; + +} // namespace photons +} // namespace singularity + +#endif // SINGULARITY_OPAC_PHOTONS_POWERLAW_S_OPACITY_PHOTONS_ diff --git a/singularity-opac/photons/s_opac_photons.hpp b/singularity-opac/photons/s_opac_photons.hpp index 472ffa8..257837e 100644 --- a/singularity-opac/photons/s_opac_photons.hpp +++ b/singularity-opac/photons/s_opac_photons.hpp @@ -21,6 +21,7 @@ #include #include #include +#include #include namespace singularity { @@ -29,9 +30,11 @@ namespace photons { using ScaleFreeS = GraySOpacity; using GrayS = GraySOpacity<>; using ThomsonS = ThomsonSOpacity<>; +using PowerLawS = PowerLawSOpacity<>; -using SOpacity = impl::S_Variant, NonCGSUnitsS>; +using SOpacity = impl::S_Variant, NonCGSUnitsS, + NonCGSUnitsS>; } // namespace photons } // namespace singularity diff --git a/test/test_powerlaw_opacities.cpp b/test/test_powerlaw_opacities.cpp index 73e39b1..a7c770a 100644 --- a/test/test_powerlaw_opacities.cpp +++ b/test/test_powerlaw_opacities.cpp @@ -1,5 +1,5 @@ // ====================================================================== -// © 2024. Triad National Security, LLC. All rights reserved. This +// © 2024-2026. 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. @@ -13,6 +13,8 @@ // publicly, and to permit others to do so. // ====================================================================== +// This file was made in part with generative AI. + #include #include @@ -23,6 +25,7 @@ #include #include +#include using namespace singularity; @@ -45,9 +48,11 @@ PORTABLE_INLINE_FUNCTION Real CalcFrequency(const int n, const Real nu_min, } constexpr Real EPS_TEST = 1e-3; +constexpr Real EPS_EXACT = 1.e-12; constexpr Real rho_exp = 2.; constexpr Real temp_exp = 2.5; +constexpr Real nu_exp = 2.; TEST_CASE("Scale free power law photon opacities", "[PowerLawScaleFreePhotonOpacities]") { @@ -196,3 +201,236 @@ TEST_CASE("CGS power law photon opacities", "[PowerLawCGSPhotonOpacities]") { } } +TEST_CASE("Frequency-dependent 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 / Hz^nu_exp + constexpr Real nu_ref = 1.e12; // Hz + + WHEN("We initialize a frequency-dependent CGS power law photon opacity") { + photons::PowerLaw opac_host(kappa0, rho_exp, temp_exp, nu_exp, nu_ref); + photons::Opacity opac = opac_host.GetOnDevice(); + + THEN("The monochromatic opacity and emissivity follow the frequency power " + "law") { + 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 frequency-dependent emissivities", 0, n_nu, + PORTABLE_LAMBDA(const int &i) { + const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu); + const Real alpha = opac.AbsorptionCoefficient(rho, temp, 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) * + std::pow(nu / nu_ref, nu_exp); + const Real alpha_expected = rho * kappa; + const Real jnu_expected = + alpha_expected * opac.ThermalDistributionOfTNu(temp, nu); + if (FractionalDifference(alpha, alpha_expected) > EPS_EXACT) { + n_wrong_d() += 1; + } + if (FractionalDifference(jnu, jnu_expected) > EPS_EXACT) { + n_wrong_d() += 1; + } + if (FractionalDifference(Jnu, 4 * M_PI * jnu) > EPS_EXACT) { + 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 frequency-dependent 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, nu_exp, nu_ref), time_unit, + mass_unit, length_unit, temp_unit); + photons::Opacity opac = opac_host.GetOnDevice(); + photons::PowerLaw opac_cgs_host(kappa0, rho_exp, temp_exp, nu_exp, nu_ref); + photons::Opacity opac_cgs = opac_cgs_host.GetOnDevice(); + + THEN("The frequency dependence is preserved through the non-CGS wrapper") { + 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 frequency-dependent non-cgs emissivities", 0, n_nu, + PORTABLE_LAMBDA(const int &i) { + const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu); + const Real alpha = opac.AbsorptionCoefficient( + rho / rho_unit, temp / temp_unit, nu * time_unit); + 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 alpha_cgs = + opac_cgs.AbsorptionCoefficient(rho, temp, nu); + const Real jnu_cgs = opac_cgs.EmissivityPerNuOmega(rho, temp, nu); + const Real Jnu_cgs = opac_cgs.EmissivityPerNu(rho, temp, nu); + if (FractionalDifference(alpha / length_unit, alpha_cgs) > + EPS_EXACT) { + n_wrong_d() += 1; + } + if (FractionalDifference(jnu * j_unit, jnu_cgs) > EPS_EXACT) { + n_wrong_d() += 1; + } + if (FractionalDifference(Jnu * j_unit, Jnu_cgs) > EPS_EXACT) { + 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("Frequency-dependent CGS power law photon scattering opacities", + "[PowerLawCGSPhotonSOpacities]") { + 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 / Hz^nu_exp + constexpr Real nu_ref = 1.e12; // Hz + + WHEN("We initialize a frequency-dependent CGS power law photon " + "scattering opacity") { + photons::PowerLawS opac_host(kappa0, rho_exp, temp_exp, nu_exp, nu_ref); + photons::SOpacity opac = opac_host.GetOnDevice(); + + THEN("The monochromatic opacity follows the frequency power law") { + 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 frequency-dependent power-law scattering", 0, n_nu, + PORTABLE_LAMBDA(const int &i) { + const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu); + const Real alpha = opac.TotalScatteringCoefficient(rho, temp, nu); + const Real kappa = opac.TotalCrossSection(rho, temp, nu); + const Real dcrss = opac.DifferentialCrossSection(rho, temp, nu, 1.0); + const Real kappa_expected = kappa0 * std::pow(rho, rho_exp) * + std::pow(temp, temp_exp) * + std::pow(nu / nu_ref, nu_exp); + const Real alpha_expected = rho * kappa_expected; + if (FractionalDifference(kappa, kappa_expected) > EPS_EXACT) { + n_wrong_d() += 1; + } + if (FractionalDifference(alpha, alpha_expected) > EPS_EXACT) { + n_wrong_d() += 1; + } + if (FractionalDifference(kappa, 4 * M_PI * dcrss) > EPS_EXACT) { + 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 frequency-dependent CGS power law photon scattering " + "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 kappa_unit = length_unit * length_unit; + + photons::NonCGSUnitsS opac_host( + photons::PowerLawS(kappa0, rho_exp, temp_exp, nu_exp, nu_ref), time_unit, + mass_unit, length_unit, temp_unit); + photons::SOpacity opac = opac_host.GetOnDevice(); + photons::PowerLawS opac_cgs_host(kappa0, rho_exp, temp_exp, nu_exp, nu_ref); + photons::SOpacity opac_cgs = opac_cgs_host.GetOnDevice(); + + THEN("The frequency dependence is preserved through the non-CGS wrapper") { + 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 frequency-dependent non-cgs power-law scattering", 0, n_nu, + PORTABLE_LAMBDA(const int &i) { + const Real nu = CalcFrequency(i, nu_min, nu_max, n_nu); + const Real alpha = opac.TotalScatteringCoefficient( + rho / rho_unit, temp / temp_unit, nu * time_unit); + const Real kappa = opac.TotalCrossSection( + rho / rho_unit, temp / temp_unit, nu * time_unit); + const Real dcrss = opac.DifferentialCrossSection( + rho / rho_unit, temp / temp_unit, nu * time_unit, 1.0); + const Real alpha_cgs = opac_cgs.TotalScatteringCoefficient( + rho, temp, nu); + const Real kappa_cgs = opac_cgs.TotalCrossSection( + rho, temp, nu); + const Real dcrss_cgs = opac_cgs.DifferentialCrossSection( + rho, temp, nu, 1.0); + if (FractionalDifference(kappa * kappa_unit, kappa_cgs) > EPS_EXACT) { + n_wrong_d() += 1; + } + if (FractionalDifference(alpha / length_unit, alpha_cgs) > EPS_EXACT) { + n_wrong_d() += 1; + } + if (FractionalDifference(dcrss * kappa_unit, dcrss_cgs) > EPS_EXACT) { + 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(); + opac_cgs.Finalize(); + } +}