From fdb05d9e31545c704cc7f8cd49763fa4ad9ce6fe Mon Sep 17 00:00:00 2001 From: Benjamin Ransom Ryan Date: Wed, 25 Sep 2024 12:31:37 -0600 Subject: [PATCH] done --- .../photons/powerlaw_opacity_photons.hpp | 26 ++++++++------- test/test_powerlaw_opacities.cpp | 32 ++++++++++--------- 2 files changed, 32 insertions(+), 26 deletions(-) diff --git a/singularity-opac/photons/powerlaw_opacity_photons.hpp b/singularity-opac/photons/powerlaw_opacity_photons.hpp index 736f639..7efe5e6 100644 --- a/singularity-opac/photons/powerlaw_opacity_photons.hpp +++ b/singularity-opac/photons/powerlaw_opacity_photons.hpp @@ -31,18 +31,19 @@ template class PowerLawOpacity { public: PowerLawOpacity() = default; - PowerLawOpacity(const Real kappa0, const Real A, const Real B) - : kappa0_(kappa0), A_(A), B_(B) {} + 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 A, const Real B) - : dist_(dist), kappa0_(kappa0), A_(A), B_(B) {} + 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 A = %g B = %g\n", kappa0_, A_, B_); + printf("Power law opacity. kappa0 = %g rho_exp = %g temp_exp = %g\n", + kappa0_, rho_exp_, temp_exp_); } inline void Finalize() noexcept {} @@ -85,7 +86,9 @@ 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, A_) * std::pow(temp, B_)) * Bnu; + return rho * + (kappa0_ * std::pow(rho, rho_exp_) * std::pow(temp, temp_exp_)) * + Bnu; } template @@ -118,13 +121,14 @@ class PowerLawOpacity { 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; + 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, A_) * std::pow(temp, B_)) * + return (kappa0_ * std::pow(rho, rho_exp_) * std::pow(temp, temp_exp_)) * dist_.ThermalNumberDistributionOfT(temp, lambda); } @@ -169,9 +173,9 @@ class PowerLawOpacity { } private: - Real kappa0_; // Opacity scale. Units of cm^2/g - Real A_; // Power law index of density - Real B_; // Power law index of temperature + 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_; }; diff --git a/test/test_powerlaw_opacities.cpp b/test/test_powerlaw_opacities.cpp index de83ba4..73e39b1 100644 --- a/test/test_powerlaw_opacities.cpp +++ b/test/test_powerlaw_opacities.cpp @@ -46,19 +46,20 @@ PORTABLE_INLINE_FUNCTION Real CalcFrequency(const int n, const Real nu_min, 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.e0; - constexpr Real temp = 1.e0; + 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; - constexpr Real A = 2.; - constexpr Real B = 2.5; - photons::PowerLawScaleFree opac_host(kappa0, A, B); + 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") { @@ -74,7 +75,8 @@ TEST_CASE("Scale free power law photon opacities", 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, A) * std::pow(temp, B); + 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) { @@ -96,17 +98,15 @@ TEST_CASE("Scale free power law photon opacities", } TEST_CASE("CGS power law photon opacities", "[PowerLawCGSPhotonOpacities]") { - constexpr Real rho = 1.e0; // g/cc + 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 - constexpr Real A = 2.; - constexpr Real B = 2.5; WHEN("We initialize a CGS power law photon opacity") { - photons::PowerLaw opac_host(kappa0, A, B); + 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") { @@ -122,7 +122,8 @@ TEST_CASE("CGS power law photon opacities", "[PowerLawCGSPhotonOpacities]") { 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, A) * std::pow(temp, B); + 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) { @@ -152,10 +153,10 @@ TEST_CASE("CGS power law photon opacities", "[PowerLawCGSPhotonOpacities]") { constexpr Real j_unit = mass_unit / (length_unit * time_unit * time_unit); photons::NonCGSUnits opac_host( - photons::PowerLaw(kappa0, A, B), time_unit, mass_unit, length_unit, - temp_unit); + 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, A, B); + 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") { @@ -173,7 +174,8 @@ TEST_CASE("CGS power law photon opacities", "[PowerLawCGSPhotonOpacities]") { 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, A) * std::pow(temp, B); + 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) {