From 425de8943fe29e3ebbfa96fcd9288e08823f636c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 26 Apr 2019 11:04:59 -0500 Subject: [PATCH] Minor optimizations for WindowedMultipole --- docs/source/io_formats/data_wmp.rst | 6 +++--- include/openmc/wmp.h | 2 +- src/physics.cpp | 2 +- src/wmp.cpp | 20 +++++++++++--------- 4 files changed, 16 insertions(+), 14 deletions(-) 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/wmp.h b/include/openmc/wmp.h index 2378f89da..6bf14e51b 100644 --- a/include/openmc/wmp.h +++ b/include/openmc/wmp.h @@ -65,7 +65,7 @@ public: 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 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) diff --git a/src/physics.cpp b/src/physics.cpp index 2636207fb..201a893c6 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -972,7 +972,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..10921fca3 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_); @@ -41,6 +42,7 @@ WindowedMultipole::WindowedMultipole(hid_t group) // 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_); @@ -51,7 +53,7 @@ WindowedMultipole::WindowedMultipole(hid_t group) // 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_ + "."); } @@ -79,9 +81,9 @@ WindowedMultipole::evaluate(double E, double sqrtkT) // 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; + static_cast((sqrtE - std::sqrt(E_min_)) * inv_spacing_)); + int startw = windows_(i_window, 0); + int endw = windows_(i_window, 1); // Initialize the ouptut cross sections double sig_s = 0.0; @@ -123,7 +125,7 @@ WindowedMultipole::evaluate(double E, double sqrtkT) // If at 0K, use asymptotic form. for (int i_pole = startw; i_pole <= endw; ++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_) { @@ -166,9 +168,9 @@ 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_; + int startw = windows_(i_window, 0); + int endw = windows_(i_window, 1); // Initialize the ouptut cross sections. double sig_s = 0.0;