From 07feb281d1c89327f84dd999ebaa24fd17a0e432 Mon Sep 17 00:00:00 2001 From: "Alex R. Long" Date: Mon, 16 Aug 2021 09:37:49 -0600 Subject: [PATCH 1/3] Make basic CDI functions constexpr + Move functions that don't require any class member data outside CDI class + Mark basic integration functions as constexpr + Remove files that formerly held low-level planck integration functions --- src/cdi/CDI.cc | 8 +- src/cdi/CDI.hh | 317 +++++++++++------- ..._Integrate_Rosseland_Planckian_Spectrum.cc | 75 ----- src/cdi/CDI_integratePlanckSpectrum.cc | 54 --- src/cdi/CDI_integrateRosselandSpectrum.cc | 49 --- src/cdi/test/tCDI.cc | 52 +-- .../Pseudo_Line_Analytic_MultigroupOpacity.cc | 4 +- src/compton_tools/cskrw.cc | 4 +- 8 files changed, 230 insertions(+), 333 deletions(-) delete mode 100644 src/cdi/CDI_Integrate_Rosseland_Planckian_Spectrum.cc delete mode 100644 src/cdi/CDI_integratePlanckSpectrum.cc delete mode 100644 src/cdi/CDI_integrateRosselandSpectrum.cc diff --git a/src/cdi/CDI.cc b/src/cdi/CDI.cc index 4e294b496d..446045fa43 100644 --- a/src/cdi/CDI.cc +++ b/src/cdi/CDI.cc @@ -112,7 +112,7 @@ double CDI::integratePlanckSpectrum(size_t const groupIndex, double const T) { const double upper_bound = frequencyGroupBoundaries[groupIndex]; Check(upper_bound > lower_bound); - double integral = integratePlanckSpectrum(lower_bound, upper_bound, T); + double integral = rtt_cdi::integratePlanckSpectrum(lower_bound, upper_bound, T); Ensure(integral >= 0.0); Ensure(integral <= 1.0); @@ -136,7 +136,7 @@ double CDI::integratePlanckSpectrum(const double T) { Check(upper_bound > lower_bound); // calculate the integral - double integral = integratePlanckSpectrum(lower_bound, upper_bound, T); + double integral = rtt_cdi::integratePlanckSpectrum(lower_bound, upper_bound, T); Ensure(integral >= 0.0 && integral <= 1.0); @@ -165,7 +165,7 @@ double CDI::integrateRosselandSpectrum(size_t const groupIndex, double const T) const double highFreq = frequencyGroupBoundaries[groupIndex]; Check(highFreq > lowFreq); - double rosseland = integrateRosselandSpectrum(lowFreq, highFreq, T); + double rosseland = rtt_cdi::integrateRosselandSpectrum(lowFreq, highFreq, T); return rosseland; } @@ -196,7 +196,7 @@ void CDI::integrate_Rosseland_Planckian_Spectrum(const size_t groupIndex, const Check(highFreq > lowFreq); // Call the general frequency version - integrate_Rosseland_Planckian_Spectrum(lowFreq, highFreq, T, planck, rosseland); + rtt_cdi::integrate_Rosseland_Planckian_Spectrum(lowFreq, highFreq, T, planck, rosseland); return; } diff --git a/src/cdi/CDI.hh b/src/cdi/CDI.hh index 79e7e687f4..91c5d99e59 100644 --- a/src/cdi/CDI.hh +++ b/src/cdi/CDI.hh @@ -47,7 +47,9 @@ double constexpr coeff_21 = 43.867 / 107290978560589824.0; double constexpr coeff = 0.1539897338202651; // 15/pi^4 double constexpr NORM_FACTOR = 0.25 * coeff; // 15/(4*pi^4); +} // namespace +namespace rtt_cdi { //------------------------------------------------------------------------------------------------// /*! * \fn inline double taylor_series_planck(double x) @@ -76,7 +78,7 @@ double constexpr NORM_FACTOR = 0.25 * coeff; // 15/(4*pi^4); * \param[in] x The point at which the Planck integral is evaluated. * \return The integral value. */ -static inline double taylor_series_planck(double x) { +constexpr double taylor_series_planck(double x) { Require(x >= 0.0); double const xsqrd = x * x; @@ -103,7 +105,7 @@ static inline double taylor_series_planck(double x) { * \brief return the 10-term Polylogarithmic expansion (minus one) for the Planck integral given * \f$ x \f$ and \f$ e^{-x} \f$ (for efficiency) */ -static double polylog_series_minus_one_planck(double const x, double const eix) { +constexpr double polylog_series_minus_one_planck(double const x, double const eix) { Require(x >= 0.0); Require(x < 1.0e154); // value will be squared, make sure it's less than sqrt of max double Require(rtt_dsxx::soft_equiv(std::exp(-x), eix)); @@ -176,6 +178,43 @@ static double polylog_series_minus_one_planck(double const x, double const eix) return poly; } +//------------------------------------------------------------------------------------------------// +/** + * \brief Integrate the normalized Planckian spectrum from 0 to \f$ x (\frac{h\nu}{kT}) \f$. + * + * \param scaled_freq upper integration limit, scaled by the temperature. + * \param exp_scaled_freq upper integration limit, scaled by an exponential function. + * \return integrated normalized Planckian from 0 to x \f$(\frac{h\nu}{kT})\f$ + * + * There are 3 cases to consider: + * 1. nu/T is very large + * If nu/T is large enough, then the integral will be 1.0. + * 2. nu/T is small + * Represent the integral via Taylor series expansion (this will be more efficient than case 3). + * 3. All other cases. Use the polylog algorithm. + */ +constexpr double integrate_planck(double const scaled_freq) { + Require(scaled_freq >= 0); + + double const exp_scaled_freq = std::exp(-scaled_freq); + // Case 1: nu/T very large -> integral == 1.0 + if (scaled_freq > 1.0e100) + return 1.0; + + // Case 2: nu/T is sufficiently small + // FWIW the break is at about scaled_freq < 2.06192398071289 + double const taylor = taylor_series_planck(std::min(scaled_freq, 1.0e15)); + // Case 3: all other situations + double const poly = polylog_series_minus_one_planck(scaled_freq, exp_scaled_freq) + 1.0; + + // Choose between 2&3: For large enough nu/T, the next line will always select the polylog value. + double const integral = std::min(taylor, poly); + + Ensure(integral >= 0.0); + Ensure(integral <= 1.0); + return integral; +} + //------------------------------------------------------------------------------------------------// /*! * \brief Compute the difference between an integrated Planck and Rosseland curves over \f$ (0,\nu) @@ -198,7 +237,7 @@ static double polylog_series_minus_one_planck(double const x, double const eix) * \param[in] exp_freq exp(-freq) * \return The difference between the integrated Planck and Rosseland curves over \f$ (0,\nu) \f$. */ -static double Planck2Rosseland(double const freq, double const exp_freq) { +constexpr double Planck2Rosseland(double const freq, double const exp_freq) { Require(freq >= 0.0); Require(rtt_dsxx::soft_equiv(exp_freq, std::exp(-freq))); @@ -217,9 +256,155 @@ static double Planck2Rosseland(double const freq, double const exp_freq) { return NORM_FACTOR * exp_freq * freq_3 * freq / -std::expm1(-freq); } -} // end of unnamed namespace +//------------------------------------------------------------------------------------------------// +/*! \brief Integrate the normalized Planckian and Rosseland spectra from 0 to \f$ x + * (\frac{h\nu}{kT}) \f$. + * + * \param scaled_freq frequency upper integration limit scaled by temperature + * \param exp_scaled_freq upper integration limit, scaled by an exponential function. + * \param planck Variable to return the Planck integral + * \param rosseland Variable to return the Rosseland integral + */ +constexpr void integrate_planck_rosseland(double const scaled_freq, double const exp_scaled_freq, + double &planck, double &rosseland) { + Require(scaled_freq >= 0.0); + Require(rtt_dsxx::soft_equiv(exp_scaled_freq, std::exp(-scaled_freq))); -namespace rtt_cdi { + // Calculate the Planckian integral + planck = integrate_planck(scaled_freq); + + Require(planck >= 0.0); + Require(planck <= 1.0); + + double const diff_rosseland = Planck2Rosseland(scaled_freq, exp_scaled_freq); + + rosseland = planck - diff_rosseland; + + Ensure(rosseland >= 0.0); + Ensure(rosseland <= 1.0); + + return; +} + +//------------------------------------------------------------------------------------------------// +/*! + * \brief Integrate the Planckian spectrum over a frequency range. + * + * The arguments to this function must all be in consistent units. For example, if low and high are + * expressed in keV, then the temperature must also be expressed in keV. If low and high are in Hz + * and temperature is in K, then low and high must first be multiplied by Planck's constant and + * temperature by Boltzmann's constant before they are passed to this function. + * + * \param[in] low lower frequency bound. + * \param[in] high higher frequency bound. + * \param[in] T the temperature (must be greater than 0.0) + * + * \return integrated normalized Plankian from low to high + */ +constexpr double integratePlanckSpectrum(double low, double high, const double T) { + Require(low >= 0.0); + Require(high >= low); + Require(T >= 0.0); + + // high/T must be < numeric_limits::max(). So, if T ~< high*min, then return early with + // zero values (assuming max() ~ 1/min()). + if (T <= high * std::numeric_limits::min()) + return 0.0; + + // Sale the frequencies by temperature + low /= T; + high /= T; + + double integral = integrate_planck(high) - integrate_planck(low); + + Ensure(integral >= 0.0); + Ensure(integral <= 1.0); + + return integral; +} + +//------------------------------------------------------------------------------------------------// +/*! + * \brief Integrate the Planckian and Rosseland spectrum over a frequency range. + * + * The arguments to this function must all be in consistent units. For example, if low and high are + * expressed in keV, then the temperature must also be expressed in keV. If low and high are in Hz + * and temperature is in K, then low and high must first be multiplied by Planck's constant and + * temperature by Boltzmann's constant before they are passed to this function. + * + * \param low Lower limit of frequency range. + * \param high Higher limit of frequency range. + * \param[in] T Temperature (must be greater than 0.0). + * \param[out] planck On return, contains the integrated normalized Planckian from low to high. + * \param[out] rosseland On return, contains the integrated normalized Rosseland from low to high + * + * \f[ + * planck(T) = \int_{\nu_1}^{\nu_2}{B(\nu,T)d\nu} + * rosseland(T) = \int_{\nu_1}^{\nu_2}{\frac{\partial B(\nu,T)}{\partial T}d\nu} + * \f] + */ +constexpr void integrate_Rosseland_Planckian_Spectrum(double low, double high, double const T, + double &planck, double &rosseland) { + Require(low >= 0.0); + Require(high >= low); + Require(T >= 0.0); + + // high/T must be < numeric_limits::max(). So, if T ~< high*min, then return early with + // zero values (assuming max() ~ 1/min()). + if (T <= high * std::numeric_limits::min()) { + planck = 0.0; + rosseland = 0.0; + return; + } + + double planck_low = 0.0; + double planck_high = 0.0; + double rosseland_low = 0.0; + double rosseland_high = 0.0; + + // Sale the frequencies by temperature + low /= T; + high /= T; + + double const exp_low = std::exp(-low); + double const exp_high = std::exp(-high); + + integrate_planck_rosseland(low, exp_low, planck_low, rosseland_low); + integrate_planck_rosseland(high, exp_high, planck_high, rosseland_high); + + planck = planck_high - planck_low; + rosseland = rosseland_high - rosseland_low; + + return; +} + +//------------------------------------------------------------------------------------------------// +/*! + * \brief Integrate the Rosseland spectrum over a frequency range. + * + * The arguments to this function must all be in consistent units. For example, if low and high are + * expressed in keV, then the temperature must also be expressed in keV. If low and high are in Hz + * and temperature is in K, then low and high must first be multiplied by Planck's constant and + * temperature by Boltzmann's constant before they are passed to this function. + * + * \param low lower frequency bound. + * \param high higher frequency bound. + * \param T the temperature (must be greater than 0.0) + * + * \return integrated normalized Rosseland from low to high + */ +constexpr double integrateRosselandSpectrum(const double low, const double high, const double T) { + Require(low >= 0.0); + Require(high >= low); + Require(T >= 0.0); + + double planck = 0.0; + double rosseland = 0.0; + + integrate_Rosseland_Planckian_Spectrum(low, high, T, planck, rosseland); + + return rosseland; +} //================================================================================================// /*! @@ -262,11 +447,6 @@ namespace rtt_cdi { * * \arg integratePlanckSpectrum(double const freq, double const T) * - * \arg integrateRosselandSpectrum(double const lowf, double const hif, double T); - * - * \arg integrate_Rosseland_Planckian_Spectrum(double const lowf, const double hif, double const T, - * double& PL, double& ROSL); - * * \arg integratePlanckSpectrum(int const groupIndex, double const T); * * \arg integratePlanckSpectrum(double const T); @@ -498,18 +678,6 @@ class CDI { // IMPLELEMENTATION // ================ - //! Integrate the normalized Planckian from 0 to x (hnu/kT). - inline static double integrate_planck(double const scaled_frequency); - - //! Integrate the normalized Planckian from 0 to x (hnu/kT). - inline static double integrate_planck(double const scaled_frequency, - double const exp_scaled_freqeuency); - - //! Integrate the normalized Planckian and Rosseland from 0 to x (hnu/kT) - inline static void integrate_planck_rosseland(double const sclaed_frequency, - double const exp_scaled_frequency, double &planck, - double &rosseland); - public: // CONSTRUCTORS // ------------ @@ -599,26 +767,18 @@ public: // INTEGRATORS: // =========== - // Over a frequency range: - // ----------------------- - - //! Integrate the normalized Planckian over a frequency range. - static double integratePlanckSpectrum(double const low, double const high, double const T); - - //! Integrate the normalized Rosseland over a frequency range. - static double integrateRosselandSpectrum(double const low, double const high, double const T); - - //! Integrate the Planckian and Rosseland over a frequency range. - static void integrate_Rosseland_Planckian_Spectrum(double const low, double const high, - double const T, double &planck, - double &rosseland); - // Over a specific group: // --------------------- //! Integrate the normalized Planckian spectrum over a frequency group. static double integratePlanckSpectrum(size_t const groupIndex, double const T); + // Over the entire group spectrum: + // ------------------------------ + + //! Integrate the normalized Planckian spectrum over all frequency groups. + static double integratePlanckSpectrum(double const T); + //! Integrate the normalized Rosseland spectrum over a frequency group static double integrateRosselandSpectrum(size_t const groupIndex, double const T); @@ -626,12 +786,6 @@ public: static void integrate_Rosseland_Planckian_Spectrum(size_t const groupIndex, double const T, double &planck, double &rosseland); - // Over the entire group spectrum: - // ------------------------------ - - //! Integrate the normalized Planckian spectrum over all frequency groups. - static double integratePlanckSpectrum(double const T); - // Over a provided vector of frequency bounds at once: // -------------------------------------------------- @@ -657,85 +811,6 @@ public: // INLINE FUNCTIONS //------------------------------------------------------------------------------------------------// -//------------------------------------------------------------------------------------------------// -/** - * \brief Integrate the normalized Planckian spectrum from 0 to \f$ x (\frac{h\nu}{kT}) \f$. - * - * \param scaled_freq upper integration limit, scaled by the temperature. - * - * \return integrated normalized Planckian from 0 to x \f$(\frac{h\nu}{kT})\f$ - */ -double CDI::integrate_planck(double const scaled_freq) { - double const exp_scaled_freq = std::exp(-scaled_freq); - return CDI::integrate_planck(scaled_freq, exp_scaled_freq); -} - -//------------------------------------------------------------------------------------------------// -/** - * \brief Integrate the normalized Planckian spectrum from 0 to \f$ x (\frac{h\nu}{kT}) \f$. - * - * \param scaled_freq upper integration limit, scaled by the temperature. - * \param exp_scaled_freq upper integration limit, scaled by an exponential function. - * \return integrated normalized Planckian from 0 to x \f$(\frac{h\nu}{kT})\f$ - * - * There are 3 cases to consider: - * 1. nu/T is very large - * If nu/T is large enough, then the integral will be 1.0. - * 2. nu/T is small - * Represent the integral via Taylor series expansion (this will be more efficient than case 3). - * 3. All other cases. Use the polylog algorithm. - */ -double CDI::integrate_planck(double const scaled_freq, double const exp_scaled_freq) { - Require(scaled_freq >= 0); - - // Case 1: nu/T very large -> integral == 1.0 - if (scaled_freq > 1.0e100) - return 1.0; - - // Case 2: nu/T is sufficiently small - // FWIW the break is at about scaled_freq < 2.06192398071289 - double const taylor = taylor_series_planck(std::min(scaled_freq, 1.0e15)); - // Case 3: all other situations - double const poly = polylog_series_minus_one_planck(scaled_freq, exp_scaled_freq) + 1.0; - - // Choose between 2&3: For large enough nu/T, the next line will always select the polylog value. - double const integral = std::min(taylor, poly); - - Ensure(integral >= 0.0); - Ensure(integral <= 1.0); - return integral; -} - -//------------------------------------------------------------------------------------------------// -/*! \brief Integrate the normalized Planckian and Rosseland spectra from 0 to \f$ x - * (\frac{h\nu}{kT}) \f$. - * - * \param scaled_freq frequency upper integration limit scaled by temperature - * \param exp_scaled_freq upper integration limit, scaled by an exponential function. - * \param planck Variable to return the Planck integral - * \param rosseland Variable to return the Rosseland integral - */ -void CDI::integrate_planck_rosseland(double const scaled_freq, double const exp_scaled_freq, - double &planck, double &rosseland) { - Require(scaled_freq >= 0.0); - Require(rtt_dsxx::soft_equiv(exp_scaled_freq, std::exp(-scaled_freq))); - - // Calculate the Planckian integral - planck = integrate_planck(scaled_freq, exp_scaled_freq); - - Require(planck >= 0.0); - Require(planck <= 1.0); - - double const diff_rosseland = Planck2Rosseland(scaled_freq, exp_scaled_freq); - - rosseland = planck - diff_rosseland; - - Ensure(rosseland >= 0.0); - Ensure(rosseland <= 1.0); - - return; -} - } // end namespace rtt_cdi #endif // rtt_cdi_CDI_hh diff --git a/src/cdi/CDI_Integrate_Rosseland_Planckian_Spectrum.cc b/src/cdi/CDI_Integrate_Rosseland_Planckian_Spectrum.cc deleted file mode 100644 index 651a776ddf..0000000000 --- a/src/cdi/CDI_Integrate_Rosseland_Planckian_Spectrum.cc +++ /dev/null @@ -1,75 +0,0 @@ -//--------------------------------------------*-C++-*---------------------------------------------// -/*! - * \file cdi/CDI_Integrate_Rosseland_Planckian_Spectrum.cc - * \author Kelly Thompson - * \date Thu Jun 22 16:22:07 2000 - * \brief CDI class implementation file. - * \note Copyright (C) 2016-2020 Triad National Security, LLC., All rights reserved. */ -//------------------------------------------------------------------------------------------------// - -#include "CDI.hh" - -namespace rtt_cdi { - -//------------------------------------------------------------------------------------------------// -// Rosseland Spectrum Integrators -//------------------------------------------------------------------------------------------------// - -//------------------------------------------------------------------------------------------------// -/*! - * \brief Integrate the Planckian and Rosseland spectrum over a frequency range. - * - * The arguments to this function must all be in consistent units. For example, if low and high are - * expressed in keV, then the temperature must also be expressed in keV. If low and high are in Hz - * and temperature is in K, then low and high must first be multiplied by Planck's constant and - * temperature by Boltzmann's constant before they are passed to this function. - * - * \param low Lower limit of frequency range. - * \param high Higher limit of frequency range. - * \param[in] T Temperature (must be greater than 0.0). - * \param[out] planck On return, contains the integrated normalized Planckian from low to high. - * \param[out] rosseland On return, contains the integrated normalized Rosseland from low to high - * - * \f[ - * planck(T) = \int_{\nu_1}^{\nu_2}{B(\nu,T)d\nu} - * rosseland(T) = \int_{\nu_1}^{\nu_2}{\frac{\partial B(\nu,T)}{\partial T}d\nu} - * \f] - */ -void CDI::integrate_Rosseland_Planckian_Spectrum(double low, double high, double const T, - double &planck, double &rosseland) { - Require(low >= 0.0); - Require(high >= low); - Require(T >= 0.0); - - // high/T must be < numeric_limits::max(). So, if T ~< high*min, then return early with - // zero values (assuming max() ~ 1/min()). - if (T <= high * std::numeric_limits::min()) { - planck = 0.0; - rosseland = 0.0; - return; - } - - double planck_high, rosseland_high; - double planck_low, rosseland_low; - - // Sale the frequencies by temperature - low /= T; - high /= T; - - double const exp_low = std::exp(-low); - double const exp_high = std::exp(-high); - - integrate_planck_rosseland(low, exp_low, planck_low, rosseland_low); - integrate_planck_rosseland(high, exp_high, planck_high, rosseland_high); - - planck = planck_high - planck_low; - rosseland = rosseland_high - rosseland_low; - - return; -} - -} // end namespace rtt_cdi - -//------------------------------------------------------------------------------------------------// -// end of CDI_integrate_Rosseland_Planckian_Spectrum.cc -//------------------------------------------------------------------------------------------------// diff --git a/src/cdi/CDI_integratePlanckSpectrum.cc b/src/cdi/CDI_integratePlanckSpectrum.cc deleted file mode 100644 index b01e33caef..0000000000 --- a/src/cdi/CDI_integratePlanckSpectrum.cc +++ /dev/null @@ -1,54 +0,0 @@ -//--------------------------------------------*-C++-*---------------------------------------------// -/*! - * \file cdi/CDI_integratePlanckSpectrum.cc - * \author Kelly Thompson - * \date Thu Jun 22 16:22:07 2000 - * \brief CDI class implementation file. - * \note Copyright (C) 2016-2020 Triad National Security, LLC., All rights reserved. */ -//------------------------------------------------------------------------------------------------// - -#include "CDI.hh" - -namespace rtt_cdi { -//------------------------------------------------------------------------------------------------// -/*! - * \brief Integrate the Planckian spectrum over a frequency range. - * - * The arguments to this function must all be in consistent units. For example, if low and high are - * expressed in keV, then the temperature must also be expressed in keV. If low and high are in Hz - * and temperature is in K, then low and high must first be multiplied by Planck's constant and - * temperature by Boltzmann's constant before they are passed to this function. - * - * \param[in] low lower frequency bound. - * \param[in] high higher frequency bound. - * \param[in] T the temperature (must be greater than 0.0) - * - * \return integrated normalized Plankian from low to high - */ -double CDI::integratePlanckSpectrum(double low, double high, const double T) { - Require(low >= 0.0); - Require(high >= low); - Require(T >= 0.0); - - // high/T must be < numeric_limits::max(). So, if T ~< high*min, then return early with - // zero values (assuming max() ~ 1/min()). - if (T <= high * std::numeric_limits::min()) - return 0.0; - - // Sale the frequencies by temperature - low /= T; - high /= T; - - double integral = integrate_planck(high) - integrate_planck(low); - - Ensure(integral >= 0.0); - Ensure(integral <= 1.0); - - return integral; -} - -} // end namespace rtt_cdi - -//------------------------------------------------------------------------------------------------// -// end of CDI_integratePlanckSpectrum.cc -//------------------------------------------------------------------------------------------------// diff --git a/src/cdi/CDI_integrateRosselandSpectrum.cc b/src/cdi/CDI_integrateRosselandSpectrum.cc deleted file mode 100644 index 7d6d6d1856..0000000000 --- a/src/cdi/CDI_integrateRosselandSpectrum.cc +++ /dev/null @@ -1,49 +0,0 @@ -//--------------------------------------------*-C++-*---------------------------------------------// -/*! - * \file cdi/CDI_integrateRosselandSpectrum.cc - * \author Kelly Thompson - * \date Thu Jun 22 16:22:07 2000 - * \brief CDI class implementation file. - * \note Copyright (C) 2016-2020 Triad National Security, LLC., All rights reserved. */ -//------------------------------------------------------------------------------------------------// - -#include "CDI.hh" - -namespace rtt_cdi { - -//------------------------------------------------------------------------------------------------// -// Rosseland Spectrum Integrators -//------------------------------------------------------------------------------------------------// - -//------------------------------------------------------------------------------------------------// -/*! - * \brief Integrate the Rosseland spectrum over a frequency range. - * - * The arguments to this function must all be in consistent units. For example, if low and high are - * expressed in keV, then the temperature must also be expressed in keV. If low and high are in Hz - * and temperature is in K, then low and high must first be multiplied by Planck's constant and - * temperature by Boltzmann's constant before they are passed to this function. - * - * \param low lower frequency bound. - * \param high higher frequency bound. - * \param T the temperature (must be greater than 0.0) - * - * \return integrated normalized Rosseland from low to high - */ -double CDI::integrateRosselandSpectrum(const double low, const double high, const double T) { - Require(low >= 0.0); - Require(high >= low); - Require(T >= 0.0); - - double planck, rosseland; - - integrate_Rosseland_Planckian_Spectrum(low, high, T, planck, rosseland); - - return rosseland; -} - -} // end namespace rtt_cdi - -//------------------------------------------------------------------------------------------------// -// end of CDI_integrateRosselandSpectrum.cc -//------------------------------------------------------------------------------------------------// diff --git a/src/cdi/test/tCDI.cc b/src/cdi/test/tCDI.cc index 68bde2148a..9a4b14417f 100644 --- a/src/cdi/test/tCDI.cc +++ b/src/cdi/test/tCDI.cc @@ -163,7 +163,7 @@ void check_CDI(rtt_dsxx::UnitTest &ut, CDI const &cdi) { // Test Rosseland integration with MG opacities. { // integrate on the interval [0.5, 5.0] keV for T=1.0 keV. - double int_total1 = cdi.integrateRosselandSpectrum(0.5, 5.0, 1.0); + double int_total1 = rtt_cdi::integrateRosselandSpectrum(0.5, 5.0, 1.0); // group 1 has the same energy range. double int_total2 = cdi.integrateRosselandSpectrum(2, 1.0); if (soft_equiv(int_total1, int_total2, 1.0e-7)) { @@ -397,7 +397,7 @@ void test_planck_integration(rtt_dsxx::UnitTest &ut) { // integration is requested over a non-existent group bool caught = false; try { - CDI::integratePlanckSpectrum(1, 1.0); + rtt_cdi::CDI::integratePlanckSpectrum(1, 1.0); } catch (const rtt_dsxx::assertion &error) { ostringstream message; message << "Caught illegal Planck calculation exception: \n" @@ -410,7 +410,7 @@ void test_planck_integration(rtt_dsxx::UnitTest &ut) { } // check some planck integrals - double int_total = CDI::integratePlanckSpectrum(0.0, 100.0, 1.0); + double int_total = rtt_cdi::integratePlanckSpectrum(0.0, 100.0, 1.0); if (soft_equiv(int_total, 1.0, 3.5e-10)) { ostringstream message; message.precision(10); @@ -425,7 +425,7 @@ void test_planck_integration(rtt_dsxx::UnitTest &ut) { FAILMSG(message.str()); } - double int_1 = CDI::integratePlanckSpectrum(0.0, 5.0, 10.0); + double int_1 = rtt_cdi::integratePlanckSpectrum(0.0, 5.0, 10.0); if (soft_equiv(int_1, .00529316, 1.0e-6)) { ostringstream message; message.precision(10); @@ -440,7 +440,7 @@ void test_planck_integration(rtt_dsxx::UnitTest &ut) { FAILMSG(message.str()); } - double int_2 = CDI::integratePlanckSpectrum(0.0, .50, 10.0); + double int_2 = rtt_cdi::integratePlanckSpectrum(0.0, .50, 10.0); if (soft_equiv(int_2, 6.29674e-6, 1.0e-6)) { ostringstream message; message.precision(10); @@ -455,11 +455,11 @@ void test_planck_integration(rtt_dsxx::UnitTest &ut) { FAILMSG(message.str()); } - double int_0 = CDI::integratePlanckSpectrum(0.0, 0.0, 10.0); + double int_0 = rtt_cdi::integratePlanckSpectrum(0.0, 0.0, 10.0); if (!soft_equiv(int_0, 0.0, 1.e-6)) ITFAILS; - double int_range = CDI::integratePlanckSpectrum(0.1, 1.0, 1.0); + double int_range = rtt_cdi::integratePlanckSpectrum(0.1, 1.0, 1.0); if (!soft_equiv(int_range, 0.0345683, 3.e-5)) ITFAILS; @@ -470,13 +470,13 @@ void test_planck_integration(rtt_dsxx::UnitTest &ut) { // Extreme case 1. This should do the normal computation, but the result is zero. { - double const integrate_range_cold = CDI::integratePlanckSpectrum(0.001, 1.0, 1.0e-30); + double const integrate_range_cold = rtt_cdi::integratePlanckSpectrum(0.001, 1.0, 1.0e-30); if (!soft_equiv(integrate_range_cold, 0.0)) ITFAILS; } // Extreme case 2. T < numeric_limits::min() --> follow special logic and return zero. { - double const integrate_range_cold = CDI::integratePlanckSpectrum(0.1, 1.0, 1.0e-308); + double const integrate_range_cold = rtt_cdi::integratePlanckSpectrum(0.1, 1.0, 1.0e-308); if (!soft_equiv(integrate_range_cold, 0.0, tol)) ITFAILS; } @@ -541,10 +541,10 @@ void test_planck_integration(rtt_dsxx::UnitTest &ut) { if (CDI::getNumberFrequencyGroups() != 3) ITFAILS; - double g1_integral = CDI::integratePlanckSpectrum(1, 1.0); - double g2_integral = CDI::integratePlanckSpectrum(2, 1.0); - double g3_integral = CDI::integratePlanckSpectrum(3, 1.0); - double g_total = CDI::integratePlanckSpectrum(1.0); + double g1_integral = rtt_cdi::CDI::integratePlanckSpectrum(1, 1.0); + double g2_integral = rtt_cdi::CDI::integratePlanckSpectrum(2, 1.0); + double g3_integral = rtt_cdi::CDI::integratePlanckSpectrum(3, 1.0); + double g_total = rtt_cdi::CDI::integratePlanckSpectrum(1.0); if (soft_equiv(g1_integral, 0.00528686, 1.e-6)) PASSMSG("Group 1 integral within tolerance."); @@ -567,27 +567,27 @@ void test_planck_integration(rtt_dsxx::UnitTest &ut) { FAILMSG("Total integral over groups fails tolerance."); // Test that a zero temperature returns a zero. - if (soft_equiv(CDI::integratePlanckSpectrum(0.0, 100.0, 0.0), 0.0)) + if (soft_equiv(rtt_cdi::integratePlanckSpectrum(0.0, 100.0, 0.0), 0.0)) PASSMSG("Planck integral from hnu=0 to 100 at T=0 is zero: good!"); else FAILMSG("Planck integral from hnu=0 to 100 at T=0 is not zero: BAD!"); - if (soft_equiv(CDI::integratePlanckSpectrum(1, 0.0), 0.0)) + if (soft_equiv(rtt_cdi::CDI::integratePlanckSpectrum(1, 0.0), 0.0)) PASSMSG("Planck integral of group 1 at T=0 is zero: good!"); else FAILMSG("Planck integral of group 1 at T=0 is not zero: BAD!"); - if (soft_equiv(CDI::integratePlanckSpectrum(2, 0.0), 0.0)) + if (soft_equiv(rtt_cdi::CDI::integratePlanckSpectrum(2, 0.0), 0.0)) PASSMSG("Planck integral of group 2 at T=0 is zero: good!"); else FAILMSG("Planck integral of group 2 at T=0 is not zero: BAD!"); - if (soft_equiv(CDI::integratePlanckSpectrum(3, 0.0), 0.0)) + if (soft_equiv(rtt_cdi::CDI::integratePlanckSpectrum(3, 0.0), 0.0)) PASSMSG("Planck integral of group 3 at T=0 is zero: good!"); else FAILMSG("Planck integral of group 3 at T=0 is not zero: BAD!"); - if (soft_equiv(CDI::integratePlanckSpectrum(0.0), 0.0)) + if (soft_equiv(rtt_cdi::CDI::integratePlanckSpectrum(0.0), 0.0)) PASSMSG("Planck integral over all groups at T=0 is zero: good!"); else FAILMSG("Planck integral over all groups at T=0 is not zero: BAD!"); @@ -611,7 +611,7 @@ void test_planck_integration(rtt_dsxx::UnitTest &ut) { for (int group_index = 1; group_index <= 3; ++group_index) { - double planck_g = CDI::integratePlanckSpectrum(group_index, 1.0); + double planck_g = rtt_cdi::CDI::integratePlanckSpectrum(group_index, 1.0); if (!soft_equiv(planck[group_index - 1], planck_g)) ITFAILS; @@ -658,7 +658,7 @@ void test_rosseland_integration(rtt_dsxx::UnitTest &ut) { } // check some planck integrals - double int_total = CDI::integrateRosselandSpectrum(0, 100.0, 1.0); + double int_total = rtt_cdi::integrateRosselandSpectrum(0, 100.0, 1.0); if (soft_equiv(int_total, 1.0, 1.0e-7)) { ostringstream message; message.precision(10); @@ -673,7 +673,7 @@ void test_rosseland_integration(rtt_dsxx::UnitTest &ut) { FAILMSG(message.str()); } - double int_1 = CDI::integratePlanckSpectrum(0.1, 1.0, 1.0); + double int_1 = rtt_cdi::integratePlanckSpectrum(0.1, 1.0, 1.0); if (soft_equiv(int_1, .0345683, 1.0e-5)) { ostringstream message; message.precision(10); @@ -688,7 +688,7 @@ void test_rosseland_integration(rtt_dsxx::UnitTest &ut) { FAILMSG(message.str()); } - double int_2 = CDI::integrateRosselandSpectrum(0.1, 1.0, 1.0); + double int_2 = rtt_cdi::integrateRosselandSpectrum(0.1, 1.0, 1.0); if (soft_equiv(int_2, 0.01220025, 1.0e-5)) { ostringstream message; message.precision(10); @@ -704,7 +704,7 @@ void test_rosseland_integration(rtt_dsxx::UnitTest &ut) { } double PL, ROS; - CDI::integrate_Rosseland_Planckian_Spectrum(0.1, 1.0, 1.0, PL, ROS); + rtt_cdi::integrate_Rosseland_Planckian_Spectrum(0.1, 1.0, 1.0, PL, ROS); if (soft_equiv(PL, .0345683, 1.e-5)) { ostringstream message; message.precision(10); @@ -733,7 +733,7 @@ void test_rosseland_integration(rtt_dsxx::UnitTest &ut) { FAILMSG(message.str()); } - CDI::integrate_Rosseland_Planckian_Spectrum(0.1, 1.0, 1.0, PL, ROS); + rtt_cdi::integrate_Rosseland_Planckian_Spectrum(0.1, 1.0, 1.0, PL, ROS); if (soft_equiv(PL, .0345683, 1.e-5)) { ostringstream message; message.precision(10); @@ -767,11 +767,11 @@ void test_rosseland_integration(rtt_dsxx::UnitTest &ut) { } // Test that a zero temperature returns a zero - if (soft_equiv(CDI::integrateRosselandSpectrum(0.0, 100.0, 0.0), 0.0)) + if (soft_equiv(rtt_cdi::integrateRosselandSpectrum(0.0, 100.0, 0.0), 0.0)) PASSMSG("Rosseland integral from hnu=0 to 100 at T=0 is zero: good!"); else FAILMSG("Rosseland integral from hnu=0 to 100 at T=0 is not zero: BAD!"); - CDI::integrate_Rosseland_Planckian_Spectrum(0.0, 100.0, 0.0, PL, ROS); + rtt_cdi::integrate_Rosseland_Planckian_Spectrum(0.0, 100.0, 0.0, PL, ROS); if (soft_equiv(PL, 0.0)) PASSMSG("Rosseland call for Planck integral at T=0 is zero: good!"); else diff --git a/src/cdi_analytic/Pseudo_Line_Analytic_MultigroupOpacity.cc b/src/cdi_analytic/Pseudo_Line_Analytic_MultigroupOpacity.cc index 7a6858268c..d4bf083cf7 100644 --- a/src/cdi_analytic/Pseudo_Line_Analytic_MultigroupOpacity.cc +++ b/src/cdi_analytic/Pseudo_Line_Analytic_MultigroupOpacity.cc @@ -158,7 +158,7 @@ sf_double Pseudo_Line_Analytic_MultigroupOpacity::getOpacity(double T, double /* for (unsigned ig = 0; ig < qpoints_; ++ig) { double const x = (ig + 0.5) * (g1 - g0) / qpoints_ + g0; - double w = CDI::integrateRosselandSpectrum(g0, g1, T); + double w = rtt_cdi::integrateRosselandSpectrum(g0, g1, T); t += w / monoOpacity(x, T); b += w; @@ -194,7 +194,7 @@ sf_double Pseudo_Line_Analytic_MultigroupOpacity::getOpacity(double T, double /* for (unsigned ig = 0; ig < qpoints_; ++ig) { double const x = (ig + 0.5) * (g1 - g0) / qpoints_ + g0; - double w = CDI::integratePlanckSpectrum(g0, g1, T); + double w = rtt_cdi::integratePlanckSpectrum(g0, g1, T); t += w * monoOpacity(x, T); b += w; diff --git a/src/compton_tools/cskrw.cc b/src/compton_tools/cskrw.cc index 48f3cd0d2b..29d3982961 100644 --- a/src/compton_tools/cskrw.cc +++ b/src/compton_tools/cskrw.cc @@ -273,7 +273,7 @@ void Dense_Compton_Data::compute_nonlinear_difference() { for (UINT64 g = 0; g < numGroups; ++g) { const FP Elow = groupBdrs[g]; const FP Ehigh = groupBdrs[g + 1]; - bg[g] = rtt_cdi::CDI::integratePlanckSpectrum(Elow, Ehigh, T); + bg[g] = rtt_cdi::integratePlanckSpectrum(Elow, Ehigh, T); bgsum += bg[g]; } bgsum = bgsum > 0.0 ? bgsum : 1.0; @@ -917,7 +917,7 @@ void read_csk_files(std::string const &basename, int verbosity) { for (UINT64 g = 0; g < dat.numGroups; ++g) { const FP Elow = dat.groupBdrs[g]; const FP Ehigh = dat.groupBdrs[g + 1]; - bg[g] = rtt_cdi::CDI::integratePlanckSpectrum(Elow, Ehigh, T); + bg[g] = rtt_cdi::integratePlanckSpectrum(Elow, Ehigh, T); bgsum += bg[g]; if (verbosity > 1) From 000b26c91f13018ee81bd7609b38055fd533eb58 Mon Sep 17 00:00:00 2001 From: Alex Long Date: Tue, 17 Aug 2021 14:47:38 -0600 Subject: [PATCH 2/3] Don't use .data() yet, requires C++17 --- src/cdi/CDI.hh | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/cdi/CDI.hh b/src/cdi/CDI.hh index 91c5d99e59..29632c20fc 100644 --- a/src/cdi/CDI.hh +++ b/src/cdi/CDI.hh @@ -123,7 +123,10 @@ constexpr double polylog_series_minus_one_planck(double const x, double const ei 0.1111111111111111, // 1/9 0.1000000000000000 // 1/10 }; - double const *curr_inv = i_plus_two_inv.data(); + static_assert(i_plus_two_inv.size() > 0, "i_plus_two_inv must be non-zero size"); + // Note: in C++ 17 the array.data() method is constexpr and should be used here and the static + // assert can be removed + double const *curr_inv = &i_plus_two_inv[0]; // initialize to what would have been the calculation of the i=1 term. This saves a number of "mul // by one" ops. From 60200f6c55ec11d85add0e86b1a0511dc1041bb8 Mon Sep 17 00:00:00 2001 From: Alex Long Date: Tue, 17 Aug 2021 15:03:14 -0600 Subject: [PATCH 3/3] Make CDI functions available on device + Hard code std::pow of max double as there's no constexpr pow in msvc + Add device/config to Cmake for CDI --- src/cdi/CDI.hh | 47 ++++++++++++++++++++++-------------------- src/cdi/CMakeLists.txt | 5 +++++ 2 files changed, 30 insertions(+), 22 deletions(-) diff --git a/src/cdi/CDI.hh b/src/cdi/CDI.hh index 29632c20fc..da4a01090e 100644 --- a/src/cdi/CDI.hh +++ b/src/cdi/CDI.hh @@ -15,6 +15,7 @@ #include "EoS.hh" #include "GrayOpacity.hh" #include "MultigroupOpacity.hh" +#include "device/config.h" #include "ds++/Constexpr_Functions.hh" #include "ds++/FMA.hh" #include "ds++/Soft_Equivalence.hh" @@ -78,7 +79,7 @@ namespace rtt_cdi { * \param[in] x The point at which the Planck integral is evaluated. * \return The integral value. */ -constexpr double taylor_series_planck(double x) { +GPU_HOST_DEVICE inline double taylor_series_planck(double x) { Require(x >= 0.0); double const xsqrd = x * x; @@ -105,7 +106,7 @@ constexpr double taylor_series_planck(double x) { * \brief return the 10-term Polylogarithmic expansion (minus one) for the Planck integral given * \f$ x \f$ and \f$ e^{-x} \f$ (for efficiency) */ -constexpr double polylog_series_minus_one_planck(double const x, double const eix) { +GPU_HOST_DEVICE inline double polylog_series_minus_one_planck(double const x, double const eix) { Require(x >= 0.0); Require(x < 1.0e154); // value will be squared, make sure it's less than sqrt of max double Require(rtt_dsxx::soft_equiv(std::exp(-x), eix)); @@ -185,8 +186,7 @@ constexpr double polylog_series_minus_one_planck(double const x, double const ei /** * \brief Integrate the normalized Planckian spectrum from 0 to \f$ x (\frac{h\nu}{kT}) \f$. * - * \param scaled_freq upper integration limit, scaled by the temperature. - * \param exp_scaled_freq upper integration limit, scaled by an exponential function. + * \param[in] scaled_freq upper integration limit, scaled by the temperature. * \return integrated normalized Planckian from 0 to x \f$(\frac{h\nu}{kT})\f$ * * There are 3 cases to consider: @@ -196,7 +196,7 @@ constexpr double polylog_series_minus_one_planck(double const x, double const ei * Represent the integral via Taylor series expansion (this will be more efficient than case 3). * 3. All other cases. Use the polylog algorithm. */ -constexpr double integrate_planck(double const scaled_freq) { +GPU_HOST_DEVICE inline double integrate_planck(double const scaled_freq) { Require(scaled_freq >= 0); double const exp_scaled_freq = std::exp(-scaled_freq); @@ -240,13 +240,13 @@ constexpr double integrate_planck(double const scaled_freq) { * \param[in] exp_freq exp(-freq) * \return The difference between the integrated Planck and Rosseland curves over \f$ (0,\nu) \f$. */ -constexpr double Planck2Rosseland(double const freq, double const exp_freq) { +GPU_HOST_DEVICE inline double Planck2Rosseland(double const freq, double const exp_freq) { Require(freq >= 0.0); Require(rtt_dsxx::soft_equiv(exp_freq, std::exp(-freq))); // Case 1: if nu/T is sufficiently large, then the evaluation is 0.0. // this evaluation also prevents overflow when evaluating (nu/T)^4. - if (freq > std::pow(std::numeric_limits::max(), 1.0 / 4.0)) + if (freq > 1.15792038e77) // hard-coded pow(numeric_limits::max(), 1/4) return 0.0; double const freq_3 = freq * freq * freq; @@ -263,13 +263,14 @@ constexpr double Planck2Rosseland(double const freq, double const exp_freq) { /*! \brief Integrate the normalized Planckian and Rosseland spectra from 0 to \f$ x * (\frac{h\nu}{kT}) \f$. * - * \param scaled_freq frequency upper integration limit scaled by temperature - * \param exp_scaled_freq upper integration limit, scaled by an exponential function. - * \param planck Variable to return the Planck integral - * \param rosseland Variable to return the Rosseland integral + * \param[in] scaled_freq frequency upper integration limit scaled by temperature + * \param[in] exp_scaled_freq upper integration limit, scaled by an exponential function. + * \param[in] planck Variable to return the Planck integral + * \param[in] rosseland Variable to return the Rosseland integral */ -constexpr void integrate_planck_rosseland(double const scaled_freq, double const exp_scaled_freq, - double &planck, double &rosseland) { +GPU_HOST_DEVICE inline void integrate_planck_rosseland(double const scaled_freq, + double const exp_scaled_freq, double &planck, + double &rosseland) { Require(scaled_freq >= 0.0); Require(rtt_dsxx::soft_equiv(exp_scaled_freq, std::exp(-scaled_freq))); @@ -304,7 +305,7 @@ constexpr void integrate_planck_rosseland(double const scaled_freq, double const * * \return integrated normalized Plankian from low to high */ -constexpr double integratePlanckSpectrum(double low, double high, const double T) { +GPU_HOST_DEVICE inline double integratePlanckSpectrum(double low, double high, const double T) { Require(low >= 0.0); Require(high >= low); Require(T >= 0.0); @@ -335,8 +336,8 @@ constexpr double integratePlanckSpectrum(double low, double high, const double T * and temperature is in K, then low and high must first be multiplied by Planck's constant and * temperature by Boltzmann's constant before they are passed to this function. * - * \param low Lower limit of frequency range. - * \param high Higher limit of frequency range. + * \param[in] low Lower limit of frequency range. + * \param[in] high Higher limit of frequency range. * \param[in] T Temperature (must be greater than 0.0). * \param[out] planck On return, contains the integrated normalized Planckian from low to high. * \param[out] rosseland On return, contains the integrated normalized Rosseland from low to high @@ -346,8 +347,9 @@ constexpr double integratePlanckSpectrum(double low, double high, const double T * rosseland(T) = \int_{\nu_1}^{\nu_2}{\frac{\partial B(\nu,T)}{\partial T}d\nu} * \f] */ -constexpr void integrate_Rosseland_Planckian_Spectrum(double low, double high, double const T, - double &planck, double &rosseland) { +GPU_HOST_DEVICE inline void integrate_Rosseland_Planckian_Spectrum(double low, double high, + double const T, double &planck, + double &rosseland) { Require(low >= 0.0); Require(high >= low); Require(T >= 0.0); @@ -390,13 +392,14 @@ constexpr void integrate_Rosseland_Planckian_Spectrum(double low, double high, d * and temperature is in K, then low and high must first be multiplied by Planck's constant and * temperature by Boltzmann's constant before they are passed to this function. * - * \param low lower frequency bound. - * \param high higher frequency bound. - * \param T the temperature (must be greater than 0.0) + * \param[in] low lower frequency bound. + * \param[in] high higher frequency bound. + * \param[in] T the temperature (must be greater than 0.0) * * \return integrated normalized Rosseland from low to high */ -constexpr double integrateRosselandSpectrum(const double low, const double high, const double T) { +GPU_HOST_DEVICE inline double integrateRosselandSpectrum(const double low, const double high, + const double T) { Require(low >= 0.0); Require(high >= low); Require(T >= 0.0); diff --git a/src/cdi/CMakeLists.txt b/src/cdi/CMakeLists.txt index 2d90427242..7967ce3750 100755 --- a/src/cdi/CMakeLists.txt +++ b/src/cdi/CMakeLists.txt @@ -15,12 +15,17 @@ file(GLOB sources *.cc) file(GLOB explicit_instantiations *_pt.cc) file(GLOB headers *.hh) +set( include_dirs PUBLIC + $ +) + # ------------------------------------------------------------------------------------------------ # # Build package library # ------------------------------------------------------------------------------------------------ # add_component_library( TARGET Lib_cdi TARGET_DEPS Lib_dsxx + INCLUDE_DIRS "${include_dirs}" LIBRARY_NAME cdi HEADERS "${headers}" SOURCES "${sources}")