diff --git a/docs/source/io_formats/data_wmp.rst b/docs/source/io_formats/data_wmp.rst index 8d7cf6ba0..e7261f3fb 100644 --- a/docs/source/io_formats/data_wmp.rst +++ b/docs/source/io_formats/data_wmp.rst @@ -17,13 +17,13 @@ Windowed Multipole Library Format If 1, Doppler broaden curve fit for window with corresponding index. If 0, do not. - **curvefit** (*double[][][]*) - Curve fit coefficients. Indexed by (reaction type, coefficient index, - window index). + Curve fit coefficients. Indexed by (window index, coefficient index, + reaction type). - **data** (*complex[][]*) Complex poles and residues. Each pole has a corresponding set of residues. For example, the :math:`i`-th pole and corresponding residues are stored as - + .. math:: \text{data}[:,i] = [\text{pole},~\text{residue}_1,~\text{residue}_2, ~\ldots] diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index e2d776dc2..79cd06325 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -203,19 +203,6 @@ extern "C" double normal_variate(double mean, double std_dev, uint64_t* seed); extern "C" double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed); -//============================================================================== -//! Doppler broadens the windowed multipole curvefit. -//! -//! The curvefit is a polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E)... -//! -//! \param E The energy to evaluate the broadening at -//! \param dopp sqrt(atomic weight ratio / kT) with kT given in eV -//! \param n The number of components to the polynomial -//! \param factors The output leading coefficient -//============================================================================== - -extern "C" void broaden_wmp_polynomials(double E, double dopp, int n, double factors[]); - //============================================================================== //! Constructs a natural cubic spline. //! diff --git a/include/openmc/wmp.h b/include/openmc/wmp.h index 2378f89da..014ee7fe7 100644 --- a/include/openmc/wmp.h +++ b/include/openmc/wmp.h @@ -35,6 +35,13 @@ constexpr std::array WMP_VERSION {1, 1}; class WindowedMultipole { public: + // Types + struct WindowInfo { + int index_start; // Index of starting pole + int index_end; // Index of ending pole + bool broaden_poly; // Whether to broaden polynomial curvefit + }; + // Constructors, destructors WindowedMultipole(hid_t group); @@ -46,7 +53,7 @@ public: //! \param E Incident neutron energy in [eV] //! \param sqrtkT Square root of temperature times Boltzmann constant //! \return Tuple of elastic scattering, absorption, and fission cross sections in [b] - std::tuple evaluate(double E, double sqrtkT); + std::tuple evaluate(double E, double sqrtkT) const; //! \brief Evaluates the windowed multipole equations for the derivative of //! cross sections in the resolved resonance regions with respect to @@ -56,20 +63,19 @@ public: //! \param sqrtkT Square root of temperature times Boltzmann constant //! \return Tuple of derivatives of elastic scattering, absorption, and //! fission cross sections in [b/K] - std::tuple evaluate_deriv(double E, double sqrtkT); + std::tuple evaluate_deriv(double E, double sqrtkT) const; // Data members std::string name_; //!< Name of nuclide - bool fissionable_; //!< Is the nuclide fissionable? - xt::xtensor, 2> data_; //!< Poles and residues - double sqrt_awr_; //!< Square root of atomic weight ratio double E_min_; //!< Minimum energy in [eV] double E_max_; //!< Maximum energy in [eV] - double spacing_; //!< Spacing in sqrt(E) space + double sqrt_awr_; //!< Square root of atomic weight ratio + double inv_spacing_; //!< 1 / spacing in sqrt(E) space int fit_order_; //!< Order of the fit - xt::xtensor windows_; //!< Indices of pole at start/end of window - xt::xtensor curvefit_; //!< Fitting function (reaction, coeff index, window index) - xt::xtensor broaden_poly_; //!< Whether to broaden curvefit + bool fissionable_; //!< Is the nuclide fissionable? + std::vector window_info_; // Information about a window + xt::xtensor curvefit_; // Curve fit coefficients (window, poly order, reaction) + xt::xtensor, 2> data_; //!< Poles and residues // Constant data static constexpr int MAX_POLY_COEFFICIENTS = @@ -91,6 +97,20 @@ void check_wmp_version(hid_t file); //! \param[in] i_nuclide Index in global nuclides array void read_multipole_data(int i_nuclide); + +//============================================================================== +//! Doppler broadens the windowed multipole curvefit. +//! +//! The curvefit is a polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E)... +//! +//! \param E The energy to evaluate the broadening at +//! \param dopp sqrt(atomic weight ratio / kT) with kT given in eV +//! \param n The number of components to the polynomial +//! \param factors The output leading coefficient +//============================================================================== + +extern "C" void broaden_wmp_polynomials(double E, double dopp, int n, double factors[]); + } // namespace openmc #endif // OPENMC_WMP_H diff --git a/src/cell.cpp b/src/cell.cpp index 6ca39749f..418258d8b 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -514,7 +514,7 @@ CSGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) cons // Check if this distance is the new minimum. if (d < min_dist) { - if (std::abs(d - min_dist) / min_dist >= FP_PRECISION) { + if (min_dist - d >= FP_PRECISION*min_dist) { min_dist = d; i_surf = -token; } diff --git a/src/math_functions.cpp b/src/math_functions.cpp index a3168fca1..ac60c6ce0 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -718,51 +718,6 @@ double watt_spectrum(double a, double b, uint64_t* seed) { } -void broaden_wmp_polynomials(double E, double dopp, int n, double factors[]) -{ - // Factors is already pre-allocated - double sqrtE = std::sqrt(E); - double beta = sqrtE * dopp; - double half_inv_dopp2 = 0.5 / (dopp * dopp); - double quarter_inv_dopp4 = half_inv_dopp2 * half_inv_dopp2; - - double erf_beta; // error function of beta - double exp_m_beta2; // exp(-beta**2) - if (beta > 6.0) { - // Save time, ERF(6) is 1 to machine precision. - // beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon. - erf_beta = 1.; - exp_m_beta2 = 0.; - } else { - erf_beta = std::erf(beta); - exp_m_beta2 = std::exp(-beta * beta); - } - - // Assume that, for sure, we'll use a second order (1/E, 1/V, const) - // fit, and no less. - - factors[0] = erf_beta / E; - factors[1] = 1. / sqrtE; - factors[2] = factors[0] * (half_inv_dopp2 + E) + exp_m_beta2 / - (beta * SQRT_PI); - - // Perform recursive broadening of high order components - for (int i = 0; i < n - 3; i++) { - double ip1_dbl = i + 1; - if (i != 0) { - factors[i + 3] = -factors[i - 1] * (ip1_dbl - 1.) * ip1_dbl * - quarter_inv_dopp4 + factors[i + 1] * - (E + (1. + 2. * ip1_dbl) * half_inv_dopp2); - } else { - // Although it's mathematically identical, factors[0] will contain - // nothing, and we don't want to have to worry about memory. - factors[i + 3] = factors[i + 1] * - (E + (1. + 2. * ip1_dbl) * half_inv_dopp2); - } - } -} - - void spline(int n, const double x[], const double y[], double z[]) { std::vector c_new(n-1); diff --git a/src/physics.cpp b/src/physics.cpp index 6e4960513..4d692abef 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -973,7 +973,7 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_ if (prn(seed) < alpha) { // With probability alpha, we sample the distribution p(y) = - // y*e^(-y). This can be done with sampling scheme C45 frmo the Monte + // y*e^(-y). This can be done with sampling scheme C45 from the Monte // Carlo sampler beta_vt_sq = -std::log(r1*r2); diff --git a/src/wmp.cpp b/src/wmp.cpp index 708cf923f..6a56ddd0f 100644 --- a/src/wmp.cpp +++ b/src/wmp.cpp @@ -24,7 +24,8 @@ WindowedMultipole::WindowedMultipole(hid_t group) name_ = object_name(group).substr(1); // Read scalar values. - read_dataset(group, "spacing", spacing_); + read_dataset(group, "spacing", inv_spacing_); + inv_spacing_ = 1.0 / inv_spacing_; read_dataset(group, "sqrtAWR", sqrt_awr_); read_dataset(group, "E_min", E_min_); read_dataset(group, "E_max", E_max_); @@ -39,19 +40,22 @@ WindowedMultipole::WindowedMultipole(hid_t group) // Read the "windows" array and use its shape to figure out the number of // windows. - read_dataset(group, "windows", windows_); - int n_windows = windows_.shape()[0]; + xt::xtensor windows; + read_dataset(group, "windows", windows); + int n_windows = windows.shape()[0]; + windows -= 1; // Adjust to 0-based indices // Read the "broaden_poly" arrays. - read_dataset(group, "broaden_poly", broaden_poly_); - if (n_windows != broaden_poly_.shape()[0]) { + xt::xtensor broaden_poly; + read_dataset(group, "broaden_poly", broaden_poly); + if (n_windows != broaden_poly.shape()[0]) { fatal_error("broaden_poly array shape is not consistent with the windows " "array shape in WMP library for " + name_ + "."); } // Read the "curvefit" array. read_dataset(group, "curvefit", curvefit_); - if (n_windows != broaden_poly_.shape()[0]) { + if (n_windows != curvefit_.shape()[0]) { fatal_error("curvefit array shape is not consistent with the windows " "array shape in WMP library for " + name_ + "."); } @@ -63,10 +67,18 @@ WindowedMultipole::WindowedMultipole(hid_t group) "Need to compile with WindowedMultipole::MAX_POLY_COEFFICIENTS = {}", fit_order_ + 1)); } + + // Move window information into a vector + window_info_.resize(n_windows); + for (int i = 0; i < n_windows; ++i) { + window_info_[i].index_start = windows(i, 0); + window_info_[i].index_end = windows(i, 1); + window_info_[i].broaden_poly = broaden_poly[i]; + } } std::tuple -WindowedMultipole::evaluate(double E, double sqrtkT) +WindowedMultipole::evaluate(double E, double sqrtkT) const { using namespace std::complex_literals; @@ -78,10 +90,9 @@ WindowedMultipole::evaluate(double E, double sqrtkT) double invE = 1.0 / E; // Locate window containing energy - int i_window = std::min(windows_.shape()[0] - 1, - static_cast((sqrtE - std::sqrt(E_min_)) / spacing_)); - int startw = windows_(i_window, 0) - 1; - int endw = windows_(i_window, 1) - 1; + int i_window = std::min(window_info_.size() - 1, + static_cast((sqrtE - std::sqrt(E_min_)) * inv_spacing_)); + const auto& window {window_info_[i_window]}; // Initialize the ouptut cross sections double sig_s = 0.0; @@ -91,7 +102,7 @@ WindowedMultipole::evaluate(double E, double sqrtkT) // ========================================================================== // Add the contribution from the curvefit polynomial. - if (sqrtkT > 0.0 && broaden_poly_(i_window)) { + if (sqrtkT > 0.0 && window.broaden_poly) { // Broaden the curvefit. double dopp = sqrt_awr_ / sqrtkT; std::array broadened_polynomials; @@ -121,9 +132,9 @@ WindowedMultipole::evaluate(double E, double sqrtkT) if (sqrtkT == 0.0) { // If at 0K, use asymptotic form. - for (int i_pole = startw; i_pole <= endw; ++i_pole) { + for (int i_pole = window.index_start; i_pole <= window.index_end; ++i_pole) { std::complex psi_chi = -1.0i / (data_(i_pole, MP_EA) - sqrtE); - std::complex c_temp = psi_chi / E; + std::complex c_temp = psi_chi * invE; sig_s += (data_(i_pole, MP_RS) * c_temp).real(); sig_a += (data_(i_pole, MP_RA) * c_temp).real(); if (fissionable_) { @@ -133,15 +144,13 @@ WindowedMultipole::evaluate(double E, double sqrtkT) } else { // At temperature, use Faddeeva function-based form. double dopp = sqrt_awr_ / sqrtkT; - if (endw >= startw) { - for (int i_pole = startw; i_pole <= endw; ++i_pole) { - std::complex z = (sqrtE - data_(i_pole, MP_EA)) * dopp; - std::complex w_val = faddeeva(z) * dopp * invE * SQRT_PI; - sig_s += (data_(i_pole, MP_RS) * w_val).real(); - sig_a += (data_(i_pole, MP_RA) * w_val).real(); - if (fissionable_) { - sig_f += (data_(i_pole, MP_RF) * w_val).real(); - } + for (int i_pole = window.index_start; i_pole <= window.index_end; ++i_pole) { + std::complex z = (sqrtE - data_(i_pole, MP_EA)) * dopp; + std::complex w_val = faddeeva(z) * dopp * invE * SQRT_PI; + sig_s += (data_(i_pole, MP_RS) * w_val).real(); + sig_a += (data_(i_pole, MP_RA) * w_val).real(); + if (fissionable_) { + sig_f += (data_(i_pole, MP_RF) * w_val).real(); } } } @@ -150,7 +159,7 @@ WindowedMultipole::evaluate(double E, double sqrtkT) } std::tuple -WindowedMultipole::evaluate_deriv(double E, double sqrtkT) +WindowedMultipole::evaluate_deriv(double E, double sqrtkT) const { // ========================================================================== // Bookkeeping @@ -166,9 +175,8 @@ WindowedMultipole::evaluate_deriv(double E, double sqrtkT) } // Locate us - int i_window = (sqrtE - std::sqrt(E_min_)) / spacing_; - int startw = windows_(i_window, 0) - 1; - int endw = windows_(i_window, 1) - 1; + int i_window = (sqrtE - std::sqrt(E_min_)) * inv_spacing_; + const auto& window {window_info_[i_window]}; // Initialize the ouptut cross sections. double sig_s = 0.0; @@ -184,20 +192,19 @@ WindowedMultipole::evaluate_deriv(double E, double sqrtkT) // Add the contribution from the poles in this window. double dopp = sqrt_awr_ / sqrtkT; - if (endw >= startw) { - for (int i_pole = startw; i_pole <= endw; ++i_pole) { - std::complex z = (sqrtE - data_(i_pole, MP_EA)) * dopp; - std::complex w_val = -invE * SQRT_PI * 0.5 * w_derivative(z, 2); - sig_s += (data_(i_pole, MP_RS) * w_val).real(); - sig_a += (data_(i_pole, MP_RA) * w_val).real(); - if (fissionable_) { - sig_f += (data_(i_pole, MP_RF) * w_val).real(); - } + for (int i_pole = window.index_start; i_pole <= window.index_end; ++i_pole) { + std::complex z = (sqrtE - data_(i_pole, MP_EA)) * dopp; + std::complex w_val = -invE * SQRT_PI * 0.5 * w_derivative(z, 2); + sig_s += (data_(i_pole, MP_RS) * w_val).real(); + sig_a += (data_(i_pole, MP_RA) * w_val).real(); + if (fissionable_) { + sig_f += (data_(i_pole, MP_RF) * w_val).real(); } - sig_s *= -0.5*sqrt_awr_ / std::sqrt(K_BOLTZMANN) * std::pow(T, -1.5); - sig_a *= -0.5*sqrt_awr_ / std::sqrt(K_BOLTZMANN) * std::pow(T, -1.5); - sig_f *= -0.5*sqrt_awr_ / std::sqrt(K_BOLTZMANN) * std::pow(T, -1.5); } + double norm = -0.5*sqrt_awr_ / std::sqrt(K_BOLTZMANN) * std::pow(T, -1.5); + sig_s *= norm; + sig_a *= norm; + sig_f *= norm; return std::make_tuple(sig_s, sig_a, sig_f); } @@ -250,4 +257,46 @@ void read_multipole_data(int i_nuclide) file_close(file); } +void broaden_wmp_polynomials(double E, double dopp, int n, double factors[]) +{ + // Broadening of polynomials follows procedure outlined in C. Josey, P. Ducru, + // B. Forget, and K. Smith, "Windowed multipole for cross section Doppler + // broadening," J. Comput. Phys., 307, 715-727 (2016). + // https://doi.org/10.1016/j.jcp.2015.08.013 + + // Factors is already pre-allocated + double sqrtE = std::sqrt(E); + double beta = sqrtE * dopp; + double half_inv_dopp2 = 0.5 / (dopp * dopp); + double quarter_inv_dopp4 = half_inv_dopp2 * half_inv_dopp2; + + double erf_beta; // error function of beta + double exp_m_beta2; // exp(-beta**2) + if (beta > 6.0) { + // Save time, ERF(6) is 1 to machine precision. + // beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon. + erf_beta = 1.; + exp_m_beta2 = 0.; + } else { + erf_beta = std::erf(beta); + exp_m_beta2 = std::exp(-beta * beta); + } + + // Assume that, for sure, we'll use a second order (1/E, 1/V, const) + // fit, and no less. + factors[0] = erf_beta / E; + factors[1] = 1. / sqrtE; + factors[2] = factors[0] * (half_inv_dopp2 + E) + exp_m_beta2 / + (beta * SQRT_PI); + if (n > 3) factors[3] = factors[1] * (E + 3.0 * half_inv_dopp2); + + // Perform recursive broadening of high order components (Eq. 16) + for (int i = 1; i < n - 3; i++) { + double ip1_dbl = i + 1; + factors[i + 3] = -factors[i - 1] * (ip1_dbl - 1.) * ip1_dbl * + quarter_inv_dopp4 + factors[i + 1] * + (E + (1. + 2. * ip1_dbl) * half_inv_dopp2); + } +} + } // namespace openmc