From b93916beb4d961cd156465e7e20ed40cb471abb1 Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 19 May 2017 12:50:32 -0400 Subject: [PATCH 1/5] fixed several bugs in multipole python api module 1. read non-fissionable nuclide 2. evaluate xs correctly (same as the function in openmc) --- openmc/data/multipole.py | 42 +++++++++++++++++++++++++--------------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/openmc/data/multipole.py b/openmc/data/multipole.py index e4049c33c..4ffcce0bc 100644 --- a/openmc/data/multipole.py +++ b/openmc/data/multipole.py @@ -77,7 +77,7 @@ def _faddeeva(z): if np.angle(z) > 0: return wofz(z) else: - return -np.conj(wofz(z)) + return -np.conj(wofz(np.conj(z))) def _broaden_wmp_polynomials(E, dopp, n): @@ -112,7 +112,8 @@ def _broaden_wmp_polynomials(E, dopp, n): erf_beta = 1.0 exp_m_beta2 = 0.0 else: - erf_beta = np.erf(beta) + from scipy.special import erf + erf_beta = erf(beta) exp_m_beta2 = np.exp(-beta**2) # Assume that, for sure, we'll use a second order (1/E, 1/V, const) @@ -128,7 +129,7 @@ def _broaden_wmp_polynomials(E, dopp, n): # Perform recursive broadening of high order components. range(1, n-4) # replaces a do i = 1, n=3. All indices are reduced by one due to the # 1-based vs. 0-based indexing. - for i in range(1, n-4): + for i in range(1, n-2): if i != 1: factors[i+2] = (-factors[i-2] * (i - 1.0) * i * quarter_inv_dopp4 + factors[i] * (E + (1.0 + 2.0 * i) * half_inv_dopp2)) @@ -339,9 +340,9 @@ class WindowedMultipole(EqualityMixin): cv.check_type('data', data, np.ndarray) if len(data.shape) != 2: raise ValueError('Multipole data arrays must be 2D') - if data.shape[1] not in (4, 5): # 4 for RM, 5 for MLBW + if data.shape[1] not in (3, 4, 5): # 3 or 4 for RM, 4 or 5 for MLBW raise ValueError('The second dimension of multipole data arrays' - ' must have a length of either 4 or 5') + ' must have a length of 3, 4 or 5') if not np.issubdtype(data.dtype, complex): raise TypeError('Multipole data arrays must be complex dtype') self._data = data @@ -405,9 +406,9 @@ class WindowedMultipole(EqualityMixin): cv.check_type('curvefit', curvefit, np.ndarray) if len(curvefit.shape) != 3: raise ValueError('Multipole curvefit arrays must be 3D') - if curvefit.shape[2] != 3: # One each for sigT, sigA, sigF + if curvefit.shape[2] not in (2, 3): # sigT, sigA (and maybe sigF) raise ValueError('The third dimension of multipole curvefit' - ' arrays must have a length of 3') + ' arrays must have a length of 2 or 2 or 3') if not np.issubdtype(curvefit.dtype, float): raise TypeError('Multipole curvefit arrays must be float dtype') self._curvefit = curvefit @@ -434,7 +435,10 @@ class WindowedMultipole(EqualityMixin): group = group_or_filename else: h5file = h5py.File(group_or_filename, 'r') - version = h5file['version'].value[0].decode() + try: + version = h5file['version'].value.decode() + except: + version = h5file['version'].value[0].decode() if version != WMP_VERSION: raise ValueError('The given WMP data uses version ' + version + ' whereas your installation of the OpenMC ' @@ -520,7 +524,7 @@ class WindowedMultipole(EqualityMixin): """ if E < self.start_E: return (0, 0, 0) - if E >= self.end_E: return (0, 0, 0) + if E > self.end_E: return (0, 0, 0) # ====================================================================== # Bookkeeping @@ -578,14 +582,16 @@ class WindowedMultipole(EqualityMixin): * broadened_polynomials[i_poly]) sigA += (self.curvefit[i_window, i_poly, _FIT_A] * broadened_polynomials[i_poly]) - sigF += (self.curvefit[i_window, i_poly, _FIT_F] - * broadened_polynomials[i_poly]) + if self.fissionable: + sigF += (self.curvefit[i_window, i_poly, _FIT_F] + * broadened_polynomials[i_poly]) else: temp = invE for i_poly in range(self.fit_order+1): sigT += self.curvefit[i_window, i_poly, _FIT_T] * temp sigA += self.curvefit[i_window, i_poly, _FIT_A] * temp - sigF += self.curvefit[i_window, i_poly, _FIT_F] * temp + if self.fissionable: + sigF += self.curvefit[i_window, i_poly, _FIT_F] * temp temp *= sqrtE # ====================================================================== @@ -601,12 +607,14 @@ class WindowedMultipole(EqualityMixin): sigT_factor[self.l_value[i_pole]-1]).real + (self.data[i_pole, _MLBW_RX] * c_temp).real) sigA += (self.data[i_pole, _MLBW_RA] * c_temp).real - sigF += (self.data[i_pole, _MLBW_RF] * c_temp).real + if self.fissionable: + sigF += (self.data[i_pole, _MLBW_RF] * c_temp).real elif self.formalism == 'RM': sigT += (self.data[i_pole, _RM_RT] * c_temp * sigT_factor[self.l_value[i_pole]-1]).real sigA += (self.data[i_pole, _RM_RA] * c_temp).real - sigF += (self.data[i_pole, _RM_RF] * c_temp).real + if self.fissionable: + sigF += (self.data[i_pole, _RM_RF] * c_temp).real else: raise ValueError('Unrecognized/Unsupported R-matrix' ' formalism') @@ -621,12 +629,14 @@ class WindowedMultipole(EqualityMixin): sigT_factor[self.l_value[i_pole]-1] + self.data[i_pole, _MLBW_RX]) * w_val).real sigA += (self.data[i_pole, _MLBW_RA] * w_val).real - sigF += (self.data[i_pole, _MLBW_RF] * w_val).real + if self.fissionable: + sigF += (self.data[i_pole, _MLBW_RF] * w_val).real elif self.formalism == 'RM': sigT += (self.data[i_pole, _RM_RT] * w_val * sigT_factor[self.l_value[i_pole]-1]).real sigA += (self.data[i_pole, _RM_RA] * w_val).real - sigF += (self.data[i_pole, _RM_RF] * w_val).real + if self.fissionable: + sigF += (self.data[i_pole, _RM_RF] * w_val).real else: raise ValueError('Unrecognized/Unsupported R-matrix' ' formalism') From 31675e7759373a8dd13e30982c444b3dee6271c7 Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 19 May 2017 13:57:10 -0400 Subject: [PATCH 2/5] multipole evaluation xs for non-fissionable nuclide --- src/cross_section.F90 | 34 +++++++++++++++++++++++++--------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 19e44538c..ce7d97133 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -649,15 +649,19 @@ contains * broadened_polynomials(i_poly) sigA = sigA + multipole % curvefit(FIT_A, i_poly, i_window) & * broadened_polynomials(i_poly) - sigF = sigF + multipole % curvefit(FIT_F, i_poly, i_window) & - * broadened_polynomials(i_poly) + if (multipole % fissionable) then + sigF = sigF + multipole % curvefit(FIT_F, i_poly, i_window) & + * broadened_polynomials(i_poly) + end if end do else ! Evaluate as if it were a polynomial temp = invE do i_poly = 1, multipole % fit_order+1 sigT = sigT + multipole % curvefit(FIT_T, i_poly, i_window) * temp sigA = sigA + multipole % curvefit(FIT_A, i_poly, i_window) * temp - sigF = sigF + multipole % curvefit(FIT_F, i_poly, i_window) * temp + if (multipole % fissionable) then + sigF = sigF + multipole % curvefit(FIT_F, i_poly, i_window) * temp + end if temp = temp * sqrtE end do end if @@ -675,12 +679,16 @@ contains sigT_factor(multipole % l_value(i_pole))) & + real(multipole % data(MLBW_RX, i_pole) * c_temp) sigA = sigA + real(multipole % data(MLBW_RA, i_pole) * c_temp) - sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * c_temp) + if (multipole % fissionable) then + sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * c_temp) + end if else if (multipole % formalism == FORM_RM) then sigT = sigT + real(multipole % data(RM_RT, i_pole) * c_temp * & sigT_factor(multipole % l_value(i_pole))) sigA = sigA + real(multipole % data(RM_RA, i_pole) * c_temp) - sigF = sigF + real(multipole % data(RM_RF, i_pole) * c_temp) + if (multipole % fissionable) then + sigF = sigF + real(multipole % data(RM_RF, i_pole) * c_temp) + end if end if end do else @@ -694,12 +702,16 @@ contains sigT_factor(multipole % l_value(i_pole)) + & multipole % data(MLBW_RX, i_pole)) * w_val) sigA = sigA + real(multipole % data(MLBW_RA, i_pole) * w_val) - sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val) + if (multipole % fissionable) then + sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val) + end if else if (multipole % formalism == FORM_RM) then sigT = sigT + real(multipole % data(RM_RT, i_pole) * w_val * & sigT_factor(multipole % l_value(i_pole))) sigA = sigA + real(multipole % data(RM_RA, i_pole) * w_val) - sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val) + if (multipole % fissionable) then + sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val) + end if end if end do end if @@ -780,12 +792,16 @@ contains sigT_factor(multipole%l_value(i_pole)) + & multipole % data(MLBW_RX, i_pole)) * w_val) sigA = sigA + real(multipole % data(MLBW_RA, i_pole) * w_val) - sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val) + if (multipole % fissionable) then + sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val) + end if else if (multipole % formalism == FORM_RM) then sigT = sigT + real(multipole % data(RM_RT, i_pole) * w_val * & sigT_factor(multipole % l_value(i_pole))) sigA = sigA + real(multipole % data(RM_RA, i_pole) * w_val) - sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val) + if (multipole % fissionable) then + sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val) + end if end if end do sigT = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigT From 031ae7351612a9d74e3a5c6495c58ecdfcb7aa1d Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 19 May 2017 16:12:24 -0400 Subject: [PATCH 3/5] addressed @paul's comments --- openmc/data/multipole.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/openmc/data/multipole.py b/openmc/data/multipole.py index 4ffcce0bc..5ca72aaa6 100644 --- a/openmc/data/multipole.py +++ b/openmc/data/multipole.py @@ -77,7 +77,7 @@ def _faddeeva(z): if np.angle(z) > 0: return wofz(z) else: - return -np.conj(wofz(np.conj(z))) + return -np.conj(wofz(z.conjugate())) def _broaden_wmp_polynomials(E, dopp, n): @@ -126,16 +126,14 @@ def _broaden_wmp_polynomials(E, dopp, n): factors[2] = (factors[0] * (half_inv_dopp2 + E) + exp_m_beta2 / (beta * np.sqrt(np.pi))) - # Perform recursive broadening of high order components. range(1, n-4) - # replaces a do i = 1, n=3. All indices are reduced by one due to the + # Perform recursive broadening of high order components. range(1, n-2) + # replaces a do i = 1, n-3. All indices are reduced by one due to the # 1-based vs. 0-based indexing. for i in range(1, n-2): if i != 1: factors[i+2] = (-factors[i-2] * (i - 1.0) * i * quarter_inv_dopp4 + factors[i] * (E + (1.0 + 2.0 * i) * 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+2] = factors[i]*(E + (1.0 + 2.0 * i) * half_inv_dopp2) return factors @@ -408,7 +406,7 @@ class WindowedMultipole(EqualityMixin): raise ValueError('Multipole curvefit arrays must be 3D') if curvefit.shape[2] not in (2, 3): # sigT, sigA (and maybe sigF) raise ValueError('The third dimension of multipole curvefit' - ' arrays must have a length of 2 or 2 or 3') + ' arrays must have a length of 2 or 3') if not np.issubdtype(curvefit.dtype, float): raise TypeError('Multipole curvefit arrays must be float dtype') self._curvefit = curvefit From 379603a03f3890c0ecbe648616a589e026e26f7e Mon Sep 17 00:00:00 2001 From: jingang Date: Mon, 22 May 2017 16:12:36 -0400 Subject: [PATCH 4/5] imported exp, erf, pi funtion from math package --- openmc/data/multipole.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/data/multipole.py b/openmc/data/multipole.py index 5ca72aaa6..977c35b7e 100644 --- a/openmc/data/multipole.py +++ b/openmc/data/multipole.py @@ -2,6 +2,7 @@ from numbers import Integral, Real import h5py import numpy as np +from math import exp, erf, pi from six import string_types from . import WMP_VERSION @@ -112,9 +113,8 @@ def _broaden_wmp_polynomials(E, dopp, n): erf_beta = 1.0 exp_m_beta2 = 0.0 else: - from scipy.special import erf erf_beta = erf(beta) - exp_m_beta2 = np.exp(-beta**2) + exp_m_beta2 = exp(-beta**2) # Assume that, for sure, we'll use a second order (1/E, 1/V, const) # fit, and no less. @@ -124,7 +124,7 @@ def _broaden_wmp_polynomials(E, dopp, n): factors[0] = erf_beta / E factors[1] = 1.0 / sqrtE factors[2] = (factors[0] * (half_inv_dopp2 + E) - + exp_m_beta2 / (beta * np.sqrt(np.pi))) + + exp_m_beta2 / (beta * np.sqrt(pi))) # Perform recursive broadening of high order components. range(1, n-2) # replaces a do i = 1, n-3. All indices are reduced by one due to the @@ -435,7 +435,7 @@ class WindowedMultipole(EqualityMixin): h5file = h5py.File(group_or_filename, 'r') try: version = h5file['version'].value.decode() - except: + except AttributeError: version = h5file['version'].value[0].decode() if version != WMP_VERSION: raise ValueError('The given WMP data uses version ' @@ -621,7 +621,7 @@ class WindowedMultipole(EqualityMixin): # At temperature, use Faddeeva function-based form. for i_pole in range(startw, endw): Z = (sqrtE - self.data[i_pole, _MP_EA]) * dopp - w_val = _faddeeva(Z) * dopp * invE * np.sqrt(np.pi) + w_val = _faddeeva(Z) * dopp * invE * np.sqrt(pi) if self.formalism == 'MLBW': sigT += ((self.data[i_pole, _MLBW_RT] * sigT_factor[self.l_value[i_pole]-1] + From 8c4d71db9ef2789bbc28a09c2c88b4cdd90cb9d8 Mon Sep 17 00:00:00 2001 From: jingang Date: Mon, 22 May 2017 17:14:25 -0400 Subject: [PATCH 5/5] imported sqrt from math as well --- openmc/data/multipole.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/openmc/data/multipole.py b/openmc/data/multipole.py index 977c35b7e..8e7f7e18d 100644 --- a/openmc/data/multipole.py +++ b/openmc/data/multipole.py @@ -1,8 +1,8 @@ from numbers import Integral, Real +from math import exp, erf, pi, sqrt import h5py import numpy as np -from math import exp, erf, pi from six import string_types from . import WMP_VERSION @@ -102,7 +102,7 @@ def _broaden_wmp_polynomials(E, dopp, n): The value of each Doppler-broadened curvefit polynomial term. """ - sqrtE = np.sqrt(E) + sqrtE = sqrt(E) beta = sqrtE * dopp half_inv_dopp2 = 0.5 / dopp**2 quarter_inv_dopp4 = half_inv_dopp2**2 @@ -124,7 +124,7 @@ def _broaden_wmp_polynomials(E, dopp, n): factors[0] = erf_beta / E factors[1] = 1.0 / sqrtE factors[2] = (factors[0] * (half_inv_dopp2 + E) - + exp_m_beta2 / (beta * np.sqrt(pi))) + + exp_m_beta2 / (beta * sqrt(pi))) # Perform recursive broadening of high order components. range(1, n-2) # replaces a do i = 1, n-3. All indices are reduced by one due to the @@ -528,8 +528,8 @@ class WindowedMultipole(EqualityMixin): # Bookkeeping # Define some frequently used variables. - sqrtkT = np.sqrt(K_BOLTZMANN * T) - sqrtE = np.sqrt(E) + sqrtkT = sqrt(K_BOLTZMANN * T) + sqrtE = sqrt(E) invE = 1.0 / E dopp = self.sqrtAWR / sqrtkT @@ -537,7 +537,7 @@ class WindowedMultipole(EqualityMixin): # the 1-based vs. 0-based indexing. Similarly startw needs to be # decreased by 1. endw does not need to be decreased because # range(startw, endw) does not include endw. - i_window = int(np.floor((sqrtE - np.sqrt(self.start_E)) / self.spacing)) + i_window = int(np.floor((sqrtE - sqrt(self.start_E)) / self.spacing)) startw = self.w_start[i_window] - 1 endw = self.w_end[i_window] @@ -621,7 +621,7 @@ class WindowedMultipole(EqualityMixin): # At temperature, use Faddeeva function-based form. for i_pole in range(startw, endw): Z = (sqrtE - self.data[i_pole, _MP_EA]) * dopp - w_val = _faddeeva(Z) * dopp * invE * np.sqrt(pi) + w_val = _faddeeva(Z) * dopp * invE * sqrt(pi) if self.formalism == 'MLBW': sigT += ((self.data[i_pole, _MLBW_RT] * sigT_factor[self.l_value[i_pole]-1] +