diff --git a/openmc/data/multipole.py b/openmc/data/multipole.py index e4049c33c..8e7f7e18d 100644 --- a/openmc/data/multipole.py +++ b/openmc/data/multipole.py @@ -1,4 +1,5 @@ from numbers import Integral, Real +from math import exp, erf, pi, sqrt import h5py import numpy as np @@ -77,7 +78,7 @@ def _faddeeva(z): if np.angle(z) > 0: return wofz(z) else: - return -np.conj(wofz(z)) + return -np.conj(wofz(z.conjugate())) def _broaden_wmp_polynomials(E, dopp, n): @@ -101,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 @@ -112,8 +113,8 @@ def _broaden_wmp_polynomials(E, dopp, n): erf_beta = 1.0 exp_m_beta2 = 0.0 else: - erf_beta = np.erf(beta) - exp_m_beta2 = np.exp(-beta**2) + erf_beta = erf(beta) + 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. @@ -123,18 +124,16 @@ 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 * sqrt(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-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)) 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 @@ -339,9 +338,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 +404,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 3') if not np.issubdtype(curvefit.dtype, float): raise TypeError('Multipole curvefit arrays must be float dtype') self._curvefit = curvefit @@ -434,7 +433,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 AttributeError: + 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,14 +522,14 @@ 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 # 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 @@ -535,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] @@ -578,14 +580,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 +605,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') @@ -615,18 +621,20 @@ 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 * sqrt(pi) if self.formalism == 'MLBW': sigT += ((self.data[i_pole, _MLBW_RT] * 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') 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