mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge pull request #882 from liangjg/multipole_fix
Fix Multipole module
This commit is contained in:
Paul Romano
GitHub
commit
61ecb4757e
2 changed files with 60 additions and 36 deletions
|
|
@ -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')
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue