Skip to content

Commit

Permalink
done
Browse files Browse the repository at this point in the history
  • Loading branch information
@brryan
brryan committed Sep 25, 2024
1 parent d218895 commit fdb05d9
Show file tree
Hide file tree
Showing 2 changed files with 32 additions and 26 deletions.
26 changes: 15 additions & 11 deletions singularity-opac/photons/powerlaw_opacity_photons.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,18 +31,19 @@ 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 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 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 {}

Expand Down Expand Up @@ -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 <typename FrequencyIndexer, typename DataIndexer>
Expand Down Expand Up @@ -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);
}

Expand Down Expand Up @@ -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<pc> dist_;
};

Expand Down
32 changes: 17 additions & 15 deletions test/test_powerlaw_opacities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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") {
Expand All @@ -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) {
Expand All @@ -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") {
Expand All @@ -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) {
Expand Down Expand Up @@ -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<photons::PowerLaw> 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") {
Expand All @@ -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) {
Expand Down

0 comments on commit fdb05d9

Please sign in to comment.