mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge pull request #954 from smharper/multipole_div_by_zero
Cleanup naming and a divide-by-zero for multipole
This commit is contained in:
Paul Romano
GitHub
commit
4be5397264
3 changed files with 158 additions and 158 deletions
|
|
@ -404,7 +404,7 @@ 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] not in (2, 3): # sigT, sigA (and maybe sigF)
|
||||
if curvefit.shape[2] not in (2, 3): # sig_t, sig_a (maybe sig_f)
|
||||
raise ValueError('The third dimension of multipole curvefit'
|
||||
' arrays must have a length of 2 or 3')
|
||||
if not np.issubdtype(curvefit.dtype, float):
|
||||
|
|
@ -531,8 +531,6 @@ class WindowedMultipole(EqualityMixin):
|
|||
sqrtkT = sqrt(K_BOLTZMANN * T)
|
||||
sqrtE = sqrt(E)
|
||||
invE = 1.0 / E
|
||||
if sqrtkT > 0.0:
|
||||
dopp = self.sqrtAWR / sqrtkT
|
||||
|
||||
# Locate us. The i_window calc omits a + 1 present in F90 because of
|
||||
# the 1-based vs. 0-based indexing. Similarly startw needs to be
|
||||
|
|
@ -547,7 +545,7 @@ class WindowedMultipole(EqualityMixin):
|
|||
# not appear in the absorption and fission equations.
|
||||
if startw <= endw:
|
||||
twophi = np.zeros(self.num_l, dtype=np.float)
|
||||
sigT_factor = np.zeros(self.num_l, dtype=np.cfloat)
|
||||
sig_t_factor = np.zeros(self.num_l, dtype=np.cfloat)
|
||||
|
||||
for iL in range(self.num_l):
|
||||
twophi[iL] = self.pseudo_k0RS[iL] * sqrtE
|
||||
|
|
@ -562,35 +560,36 @@ class WindowedMultipole(EqualityMixin):
|
|||
twophi[iL] = twophi[iL] - np.arctan(arg)
|
||||
|
||||
twophi = 2.0 * twophi
|
||||
sigT_factor = np.cos(twophi) - 1j*np.sin(twophi)
|
||||
sig_t_factor = np.cos(twophi) - 1j*np.sin(twophi)
|
||||
|
||||
# Initialize the ouptut cross sections.
|
||||
sigT = 0.0
|
||||
sigA = 0.0
|
||||
sigF = 0.0
|
||||
sig_t = 0.0
|
||||
sig_a = 0.0
|
||||
sig_f = 0.0
|
||||
|
||||
# ======================================================================
|
||||
# Add the contribution from the curvefit polynomial.
|
||||
|
||||
if sqrtkT != 0 and self.broaden_poly[i_window]:
|
||||
# Broaden the curvefit.
|
||||
dopp = self.sqrtAWR / sqrtkT
|
||||
broadened_polynomials = _broaden_wmp_polynomials(E, dopp,
|
||||
self.fit_order + 1)
|
||||
for i_poly in range(self.fit_order+1):
|
||||
sigT += (self.curvefit[i_window, i_poly, _FIT_T]
|
||||
* broadened_polynomials[i_poly])
|
||||
sigA += (self.curvefit[i_window, i_poly, _FIT_A]
|
||||
* broadened_polynomials[i_poly])
|
||||
sig_t += (self.curvefit[i_window, i_poly, _FIT_T]
|
||||
* broadened_polynomials[i_poly])
|
||||
sig_a += (self.curvefit[i_window, i_poly, _FIT_A]
|
||||
* broadened_polynomials[i_poly])
|
||||
if self.fissionable:
|
||||
sigF += (self.curvefit[i_window, i_poly, _FIT_F]
|
||||
* broadened_polynomials[i_poly])
|
||||
sig_f += (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
|
||||
sig_t += self.curvefit[i_window, i_poly, _FIT_T] * temp
|
||||
sig_a += self.curvefit[i_window, i_poly, _FIT_A] * temp
|
||||
if self.fissionable:
|
||||
sigF += self.curvefit[i_window, i_poly, _FIT_F] * temp
|
||||
sig_f += self.curvefit[i_window, i_poly, _FIT_F] * temp
|
||||
temp *= sqrtE
|
||||
|
||||
# ======================================================================
|
||||
|
|
@ -602,45 +601,46 @@ class WindowedMultipole(EqualityMixin):
|
|||
psi_chi = -1j / (self.data[i_pole, _MP_EA] - sqrtE)
|
||||
c_temp = psi_chi / E
|
||||
if self.formalism == 'MLBW':
|
||||
sigT += ((self.data[i_pole, _MLBW_RT] * c_temp *
|
||||
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
|
||||
sig_t += ((self.data[i_pole, _MLBW_RT] * c_temp *
|
||||
sig_t_factor[self.l_value[i_pole]-1]).real
|
||||
+ (self.data[i_pole, _MLBW_RX] * c_temp).real)
|
||||
sig_a += (self.data[i_pole, _MLBW_RA] * c_temp).real
|
||||
if self.fissionable:
|
||||
sigF += (self.data[i_pole, _MLBW_RF] * c_temp).real
|
||||
sig_f += (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
|
||||
sig_t += (self.data[i_pole, _RM_RT] * c_temp *
|
||||
sig_t_factor[self.l_value[i_pole]-1]).real
|
||||
sig_a += (self.data[i_pole, _RM_RA] * c_temp).real
|
||||
if self.fissionable:
|
||||
sigF += (self.data[i_pole, _RM_RF] * c_temp).real
|
||||
sig_f += (self.data[i_pole, _RM_RF] * c_temp).real
|
||||
else:
|
||||
raise ValueError('Unrecognized/Unsupported R-matrix'
|
||||
' formalism')
|
||||
|
||||
else:
|
||||
# At temperature, use Faddeeva function-based form.
|
||||
dopp = self.sqrtAWR / sqrtkT
|
||||
for i_pole in range(startw, endw):
|
||||
Z = (sqrtE - self.data[i_pole, _MP_EA]) * dopp
|
||||
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
|
||||
sig_t += ((self.data[i_pole, _MLBW_RT] *
|
||||
sig_t_factor[self.l_value[i_pole]-1] +
|
||||
self.data[i_pole, _MLBW_RX]) * w_val).real
|
||||
sig_a += (self.data[i_pole, _MLBW_RA] * w_val).real
|
||||
if self.fissionable:
|
||||
sigF += (self.data[i_pole, _MLBW_RF] * w_val).real
|
||||
sig_f += (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
|
||||
sig_t += (self.data[i_pole, _RM_RT] * w_val *
|
||||
sig_t_factor[self.l_value[i_pole]-1]).real
|
||||
sig_a += (self.data[i_pole, _RM_RA] * w_val).real
|
||||
if self.fissionable:
|
||||
sigF += (self.data[i_pole, _RM_RF] * w_val).real
|
||||
sig_f += (self.data[i_pole, _RM_RF] * w_val).real
|
||||
else:
|
||||
raise ValueError('Unrecognized/Unsupported R-matrix'
|
||||
' formalism')
|
||||
|
||||
return sigT, sigA, sigF
|
||||
return sig_t, sig_a, sig_f
|
||||
|
||||
def __call__(self, E, T):
|
||||
"""Compute total, absorption, and fission cross sections.
|
||||
|
|
|
|||
|
|
@ -157,10 +157,9 @@ contains
|
|||
integer :: i_high ! upper logarithmic mapping index
|
||||
integer :: i_rxn ! reaction index
|
||||
integer :: j ! index in DEPLETION_RX
|
||||
real(8) :: val ! temporary xs value
|
||||
real(8) :: f ! interp factor on nuclide energy grid
|
||||
real(8) :: kT ! temperature in eV
|
||||
real(8) :: sigT, sigA, sigF ! Intermediate multipole variables
|
||||
real(8) :: sig_t, sig_a, sig_f ! Intermediate multipole variables
|
||||
|
||||
! Initialize cached cross sections to zero
|
||||
micro_xs(i_nuclide) % thermal = ZERO
|
||||
|
|
@ -179,15 +178,15 @@ contains
|
|||
! Evaluate multipole or interpolate
|
||||
if (use_mp) then
|
||||
! Call multipole kernel
|
||||
call multipole_eval(nuc % multipole, E, sqrtkT, sigT, sigA, sigF)
|
||||
call multipole_eval(nuc % multipole, E, sqrtkT, sig_t, sig_a, sig_f)
|
||||
|
||||
micro_xs(i_nuclide) % total = sigT
|
||||
micro_xs(i_nuclide) % absorption = sigA
|
||||
micro_xs(i_nuclide) % elastic = sigT - sigA
|
||||
micro_xs(i_nuclide) % total = sig_t
|
||||
micro_xs(i_nuclide) % absorption = sig_a
|
||||
micro_xs(i_nuclide) % elastic = sig_t - sig_a
|
||||
|
||||
if (nuc % fissionable) then
|
||||
micro_xs(i_nuclide) % fission = sigF
|
||||
micro_xs(i_nuclide) % nu_fission = sigF * nuc % nu(E, EMISSION_TOTAL)
|
||||
micro_xs(i_nuclide) % fission = sig_f
|
||||
micro_xs(i_nuclide) % nu_fission = sig_f * nuc % nu(E, EMISSION_TOTAL)
|
||||
else
|
||||
micro_xs(i_nuclide) % fission = ZERO
|
||||
micro_xs(i_nuclide) % nu_fission = ZERO
|
||||
|
|
@ -198,7 +197,7 @@ contains
|
|||
micro_xs(i_nuclide) % reaction(:) = ZERO
|
||||
|
||||
! Only non-zero reaction is (n,gamma)
|
||||
micro_xs(i_nuclide) % reaction(4) = sigA - sigF
|
||||
micro_xs(i_nuclide) % reaction(4) = sig_a - sig_f
|
||||
end if
|
||||
|
||||
! Ensure these values are set
|
||||
|
|
@ -613,7 +612,7 @@ contains
|
|||
! sections in the resolved resonance regions
|
||||
!===============================================================================
|
||||
|
||||
subroutine multipole_eval(multipole, E, sqrtkT, sigT, sigA, sigF)
|
||||
subroutine multipole_eval(multipole, E, sqrtkT, sig_t, sig_a, sig_f)
|
||||
type(MultipoleArray), intent(in) :: multipole ! The windowed multipole
|
||||
! object to process.
|
||||
real(8), intent(in) :: E ! The energy at which to
|
||||
|
|
@ -621,15 +620,15 @@ contains
|
|||
real(8), intent(in) :: sqrtkT ! The temperature in the form
|
||||
! sqrt(kT), at which
|
||||
! to evaluate the XS.
|
||||
real(8), intent(out) :: sigT ! Total cross section
|
||||
real(8), intent(out) :: sigA ! Absorption cross section
|
||||
real(8), intent(out) :: sigF ! Fission cross section
|
||||
real(8), intent(out) :: sig_t ! Total cross section
|
||||
real(8), intent(out) :: sig_a ! Absorption cross section
|
||||
real(8), intent(out) :: sig_f ! Fission cross section
|
||||
complex(8) :: psi_chi ! The value of the psi-chi function for the
|
||||
! asymptotic form
|
||||
complex(8) :: c_temp ! complex temporary variable
|
||||
complex(8) :: w_val ! The faddeeva function evaluated at Z
|
||||
complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole)
|
||||
complex(8) :: sigT_factor(multipole % num_l)
|
||||
complex(8) :: sig_t_factor(multipole % num_l)
|
||||
real(8) :: broadened_polynomials(multipole % fit_order + 1)
|
||||
real(8) :: sqrtE ! sqrt(E), eV
|
||||
real(8) :: invE ! 1/E, eV
|
||||
|
|
@ -647,7 +646,6 @@ contains
|
|||
! Define some frequently used variables.
|
||||
sqrtE = sqrt(E)
|
||||
invE = ONE / E
|
||||
dopp = multipole % sqrtAWR / sqrtkT
|
||||
|
||||
! Locate us.
|
||||
i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing &
|
||||
|
|
@ -657,38 +655,39 @@ contains
|
|||
|
||||
! Fill in factors.
|
||||
if (startw <= endw) then
|
||||
call compute_sigT_factor(multipole, sqrtE, sigT_factor)
|
||||
call compute_sig_t_factor(multipole, sqrtE, sig_t_factor)
|
||||
end if
|
||||
|
||||
! Initialize the ouptut cross sections.
|
||||
sigT = ZERO
|
||||
sigA = ZERO
|
||||
sigF = ZERO
|
||||
sig_t = ZERO
|
||||
sig_a = ZERO
|
||||
sig_f = ZERO
|
||||
|
||||
! ==========================================================================
|
||||
! Add the contribution from the curvefit polynomial.
|
||||
|
||||
if (sqrtkT /= ZERO .and. multipole % broaden_poly(i_window) == 1) then
|
||||
! Broaden the curvefit.
|
||||
dopp = multipole % sqrtAWR / sqrtkT
|
||||
call broaden_wmp_polynomials(E, dopp, multipole % fit_order + 1, &
|
||||
broadened_polynomials)
|
||||
do i_poly = 1, multipole % fit_order+1
|
||||
sigT = sigT + multipole % curvefit(FIT_T, i_poly, i_window) &
|
||||
sig_t = sig_t + multipole % curvefit(FIT_T, i_poly, i_window) &
|
||||
* broadened_polynomials(i_poly)
|
||||
sigA = sigA + multipole % curvefit(FIT_A, i_poly, i_window) &
|
||||
sig_a = sig_a + multipole % curvefit(FIT_A, i_poly, i_window) &
|
||||
* broadened_polynomials(i_poly)
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + multipole % curvefit(FIT_F, i_poly, i_window) &
|
||||
sig_f = sig_f + 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
|
||||
sig_t = sig_t + multipole % curvefit(FIT_T, i_poly, i_window) * temp
|
||||
sig_a = sig_a + multipole % curvefit(FIT_A, i_poly, i_window) * temp
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + multipole % curvefit(FIT_F, i_poly, i_window) * temp
|
||||
sig_f = sig_f + multipole % curvefit(FIT_F, i_poly, i_window) * temp
|
||||
end if
|
||||
temp = temp * sqrtE
|
||||
end do
|
||||
|
|
@ -703,42 +702,43 @@ contains
|
|||
psi_chi = -ONEI / (multipole % data(MP_EA, i_pole) - sqrtE)
|
||||
c_temp = psi_chi / E
|
||||
if (multipole % formalism == FORM_MLBW) then
|
||||
sigT = sigT + real(multipole % data(MLBW_RT, i_pole) * c_temp * &
|
||||
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)
|
||||
sig_t = sig_t + real(multipole % data(MLBW_RT, i_pole) * c_temp * &
|
||||
sig_t_factor(multipole % l_value(i_pole))) &
|
||||
+ real(multipole % data(MLBW_RX, i_pole) * c_temp)
|
||||
sig_a = sig_a + real(multipole % data(MLBW_RA, i_pole) * c_temp)
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * c_temp)
|
||||
sig_f = sig_f + 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)
|
||||
sig_t = sig_t + real(multipole % data(RM_RT, i_pole) * c_temp * &
|
||||
sig_t_factor(multipole % l_value(i_pole)))
|
||||
sig_a = sig_a + real(multipole % data(RM_RA, i_pole) * c_temp)
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + real(multipole % data(RM_RF, i_pole) * c_temp)
|
||||
sig_f = sig_f + real(multipole % data(RM_RF, i_pole) * c_temp)
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
else
|
||||
! At temperature, use Faddeeva function-based form.
|
||||
dopp = multipole % sqrtAWR / sqrtkT
|
||||
if (endw >= startw) then
|
||||
do i_pole = startw, endw
|
||||
Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp
|
||||
w_val = faddeeva(Z) * dopp * invE * SQRT_PI
|
||||
if (multipole % formalism == FORM_MLBW) then
|
||||
sigT = sigT + real((multipole % data(MLBW_RT, i_pole) * &
|
||||
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)
|
||||
sig_t = sig_t + real((multipole % data(MLBW_RT, i_pole) * &
|
||||
sig_t_factor(multipole % l_value(i_pole)) + &
|
||||
multipole % data(MLBW_RX, i_pole)) * w_val)
|
||||
sig_a = sig_a + real(multipole % data(MLBW_RA, i_pole) * w_val)
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val)
|
||||
sig_f = sig_f + 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)
|
||||
sig_t = sig_t + real(multipole % data(RM_RT, i_pole) * w_val * &
|
||||
sig_t_factor(multipole % l_value(i_pole)))
|
||||
sig_a = sig_a + real(multipole % data(RM_RA, i_pole) * w_val)
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val)
|
||||
sig_f = sig_f + real(multipole % data(RM_RF, i_pole) * w_val)
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
|
@ -752,7 +752,7 @@ contains
|
|||
! temperature.
|
||||
!===============================================================================
|
||||
|
||||
subroutine multipole_deriv_eval(multipole, E, sqrtkT, sigT, sigA, sigF)
|
||||
subroutine multipole_deriv_eval(multipole, E, sqrtkT, sig_t, sig_a, sig_f)
|
||||
type(MultipoleArray), intent(in) :: multipole ! The windowed multipole
|
||||
! object to process.
|
||||
real(8), intent(in) :: E ! The energy at which to
|
||||
|
|
@ -760,12 +760,12 @@ contains
|
|||
real(8), intent(in) :: sqrtkT ! The temperature in the form
|
||||
! sqrt(kT), at which to
|
||||
! evaluate the XS.
|
||||
real(8), intent(out) :: sigT ! Total cross section
|
||||
real(8), intent(out) :: sigA ! Absorption cross section
|
||||
real(8), intent(out) :: sigF ! Fission cross section
|
||||
real(8), intent(out) :: sig_t ! Total cross section
|
||||
real(8), intent(out) :: sig_a ! Absorption cross section
|
||||
real(8), intent(out) :: sig_f ! Fission cross section
|
||||
complex(8) :: w_val ! The faddeeva function evaluated at Z
|
||||
complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole)
|
||||
complex(8) :: sigT_factor(multipole % num_l)
|
||||
complex(8) :: sig_t_factor(multipole % num_l)
|
||||
real(8) :: sqrtE ! sqrt(E), eV
|
||||
real(8) :: invE ! 1/E, eV
|
||||
real(8) :: dopp ! sqrt(atomic weight ratio / kT)
|
||||
|
|
@ -781,7 +781,6 @@ contains
|
|||
! Define some frequently used variables.
|
||||
sqrtE = sqrt(E)
|
||||
invE = ONE / E
|
||||
dopp = multipole % sqrtAWR / sqrtkT
|
||||
T = sqrtkT**2 / K_BOLTZMANN
|
||||
|
||||
if (sqrtkT == ZERO) call fatal_error("Windowed multipole temperature &
|
||||
|
|
@ -795,13 +794,13 @@ contains
|
|||
|
||||
! Fill in factors.
|
||||
if (startw <= endw) then
|
||||
call compute_sigT_factor(multipole, sqrtE, sigT_factor)
|
||||
call compute_sig_t_factor(multipole, sqrtE, sig_t_factor)
|
||||
end if
|
||||
|
||||
! Initialize the ouptut cross sections.
|
||||
sigT = ZERO
|
||||
sigA = ZERO
|
||||
sigF = ZERO
|
||||
sig_t = ZERO
|
||||
sig_a = ZERO
|
||||
sig_f = ZERO
|
||||
|
||||
! TODO Polynomials: Some of the curvefit polynomials Doppler broaden so
|
||||
! rigorously we should be computing the derivative of those. But in
|
||||
|
|
@ -811,42 +810,43 @@ contains
|
|||
! ==========================================================================
|
||||
! Add the contribution from the poles in this window.
|
||||
|
||||
dopp = multipole % sqrtAWR / sqrtkT
|
||||
if (endw >= startw) then
|
||||
do i_pole = startw, endw
|
||||
Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp
|
||||
w_val = -invE * SQRT_PI * HALF * w_derivative(Z, 2)
|
||||
if (multipole % formalism == FORM_MLBW) then
|
||||
sigT = sigT + real((multipole % data(MLBW_RT, i_pole) * &
|
||||
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)
|
||||
sig_t = sig_t + real((multipole % data(MLBW_RT, i_pole) * &
|
||||
sig_t_factor(multipole%l_value(i_pole)) + &
|
||||
multipole % data(MLBW_RX, i_pole)) * w_val)
|
||||
sig_a = sig_a + real(multipole % data(MLBW_RA, i_pole) * w_val)
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val)
|
||||
sig_f = sig_f + 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)
|
||||
sig_t = sig_t + real(multipole % data(RM_RT, i_pole) * w_val * &
|
||||
sig_t_factor(multipole % l_value(i_pole)))
|
||||
sig_a = sig_a + real(multipole % data(RM_RA, i_pole) * w_val)
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val)
|
||||
sig_f = sig_f + 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
|
||||
sigA = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigA
|
||||
sigF = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigF
|
||||
sig_t = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_t
|
||||
sig_a = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_a
|
||||
sig_f = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_f
|
||||
end if
|
||||
end subroutine multipole_deriv_eval
|
||||
|
||||
!===============================================================================
|
||||
! COMPUTE_SIGT_FACTOR calculates the sigT_factor, a factor inside of the sigT
|
||||
! equation not present in the sigA and sigF equations.
|
||||
! COMPUTE_SIG_T_FACTOR calculates the sig_t_factor, a factor inside of the sig_t
|
||||
! equation not present in the sig_a and sig_f equations.
|
||||
!===============================================================================
|
||||
|
||||
subroutine compute_sigT_factor(multipole, sqrtE, sigT_factor)
|
||||
subroutine compute_sig_t_factor(multipole, sqrtE, sig_t_factor)
|
||||
type(MultipoleArray), intent(in) :: multipole
|
||||
real(8), intent(in) :: sqrtE
|
||||
complex(8), intent(out) :: sigT_factor(multipole % num_l)
|
||||
complex(8), intent(out) :: sig_t_factor(multipole % num_l)
|
||||
|
||||
integer :: iL
|
||||
real(8) :: twophi(multipole % num_l)
|
||||
|
|
@ -867,8 +867,8 @@ contains
|
|||
end do
|
||||
|
||||
twophi = 2.0_8 * twophi
|
||||
sigT_factor = cmplx(cos(twophi), -sin(twophi), KIND=8)
|
||||
end subroutine compute_sigT_factor
|
||||
sig_t_factor = cmplx(cos(twophi), -sin(twophi), KIND=8)
|
||||
end subroutine compute_sig_t_factor
|
||||
|
||||
!===============================================================================
|
||||
! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
|
||||
|
|
|
|||
|
|
@ -2382,7 +2382,7 @@ contains
|
|||
|
||||
type(Particle), intent(in) :: p
|
||||
|
||||
integer :: i, j, m
|
||||
integer :: i, j
|
||||
integer :: i_tally
|
||||
integer :: i_filt
|
||||
integer :: i_bin
|
||||
|
|
@ -3472,7 +3472,7 @@ contains
|
|||
integer :: l
|
||||
logical :: scoring_diff_nuclide
|
||||
real(8) :: flux_deriv
|
||||
real(8) :: dsigT, dsigA, dsigF, cum_dsig
|
||||
real(8) :: dsig_t, dsig_a, dsig_f, cum_dsig
|
||||
|
||||
if (score == ZERO) return
|
||||
|
||||
|
|
@ -3713,17 +3713,17 @@ contains
|
|||
if (mat % nuclide(l) == p % event_nuclide) exit
|
||||
end do
|
||||
|
||||
dsigT = ZERO
|
||||
dsig_t = ZERO
|
||||
associate (nuc => nuclides(p % event_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsigT * mat % atom_density(l) / material_xs % total)
|
||||
+ dsig_t * mat % atom_density(l) / material_xs % total)
|
||||
end associate
|
||||
else
|
||||
score = score * flux_deriv
|
||||
|
|
@ -3739,17 +3739,17 @@ contains
|
|||
if (mat % nuclide(l) == p % event_nuclide) exit
|
||||
end do
|
||||
|
||||
dsigT = ZERO
|
||||
dsigA = ZERO
|
||||
dsig_t = ZERO
|
||||
dsig_a = ZERO
|
||||
associate (nuc => nuclides(p % event_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv + (dsigT - dsigA) &
|
||||
score = score * (flux_deriv + (dsig_t - dsig_a) &
|
||||
* mat % atom_density(l) / &
|
||||
(material_xs % total - material_xs % absorption))
|
||||
end associate
|
||||
|
|
@ -3766,17 +3766,17 @@ contains
|
|||
if (mat % nuclide(l) == p % event_nuclide) exit
|
||||
end do
|
||||
|
||||
dsigA = ZERO
|
||||
dsig_a = ZERO
|
||||
associate (nuc => nuclides(p % event_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsigA * mat % atom_density(l) / material_xs % absorption)
|
||||
score = score * (flux_deriv + dsig_a * mat % atom_density(l) &
|
||||
/ material_xs % absorption)
|
||||
end associate
|
||||
else
|
||||
score = score * flux_deriv
|
||||
|
|
@ -3791,17 +3791,17 @@ contains
|
|||
if (mat % nuclide(l) == p % event_nuclide) exit
|
||||
end do
|
||||
|
||||
dsigF = ZERO
|
||||
dsig_f = ZERO
|
||||
associate (nuc => nuclides(p % event_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsigF * mat % atom_density(l) / material_xs % fission)
|
||||
+ dsig_f * mat % atom_density(l) / material_xs % fission)
|
||||
end associate
|
||||
else
|
||||
score = score * flux_deriv
|
||||
|
|
@ -3816,17 +3816,17 @@ contains
|
|||
if (mat % nuclide(l) == p % event_nuclide) exit
|
||||
end do
|
||||
|
||||
dsigF = ZERO
|
||||
dsig_f = ZERO
|
||||
associate (nuc => nuclides(p % event_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsigF * mat % atom_density(l) / material_xs % nu_fission&
|
||||
+ dsig_f * mat % atom_density(l) / material_xs % nu_fission&
|
||||
* micro_xs(p % event_nuclide) % nu_fission &
|
||||
/ micro_xs(p % event_nuclide) % fission)
|
||||
end associate
|
||||
|
|
@ -3859,8 +3859,8 @@ contains
|
|||
p % last_E <= nuc % multipole % end_E .and. &
|
||||
micro_xs(mat % nuclide(l)) % total > ZERO) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
cum_dsig = cum_dsig + dsigT * mat % atom_density(l)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig + dsig_t * mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
|
|
@ -3869,17 +3869,17 @@ contains
|
|||
+ cum_dsig / material_xs % total)
|
||||
else if (materials(p % material) % id == deriv % diff_material &
|
||||
.and. material_xs % total > ZERO) then
|
||||
dsigT = ZERO
|
||||
dsig_t = ZERO
|
||||
associate (nuc => nuclides(i_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsigT / micro_xs(i_nuclide) % total)
|
||||
+ dsig_t / micro_xs(i_nuclide) % total)
|
||||
else
|
||||
score = score * flux_deriv
|
||||
end if
|
||||
|
|
@ -3898,9 +3898,9 @@ contains
|
|||
(micro_xs(mat % nuclide(l)) % total &
|
||||
- micro_xs(mat % nuclide(l)) % absorption) > ZERO) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig &
|
||||
+ (dsigT - dsigA) * mat % atom_density(l)
|
||||
+ (dsig_t - dsig_a) * mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
|
|
@ -3910,17 +3910,17 @@ contains
|
|||
else if ( materials(p % material) % id == deriv % diff_material &
|
||||
.and. (material_xs % total - material_xs % absorption) > ZERO)&
|
||||
then
|
||||
dsigT = ZERO
|
||||
dsigA = ZERO
|
||||
dsig_t = ZERO
|
||||
dsig_a = ZERO
|
||||
associate (nuc => nuclides(i_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv + (dsigT - dsigA) &
|
||||
score = score * (flux_deriv + (dsig_t - dsig_a) &
|
||||
/ (micro_xs(i_nuclide) % total &
|
||||
- micro_xs(i_nuclide) % absorption))
|
||||
else
|
||||
|
|
@ -3940,8 +3940,8 @@ contains
|
|||
p % last_E <= nuc % multipole % end_E .and. &
|
||||
micro_xs(mat % nuclide(l)) % absorption > ZERO) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
cum_dsig = cum_dsig + dsigA * mat % atom_density(l)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig + dsig_a * mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
|
|
@ -3950,17 +3950,17 @@ contains
|
|||
+ cum_dsig / material_xs % absorption)
|
||||
else if (materials(p % material) % id == deriv % diff_material &
|
||||
.and. material_xs % absorption > ZERO) then
|
||||
dsigA = ZERO
|
||||
dsig_a = ZERO
|
||||
associate (nuc => nuclides(i_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsigA / micro_xs(i_nuclide) % absorption)
|
||||
+ dsig_a / micro_xs(i_nuclide) % absorption)
|
||||
else
|
||||
score = score * flux_deriv
|
||||
end if
|
||||
|
|
@ -3978,8 +3978,8 @@ contains
|
|||
p % last_E <= nuc % multipole % end_E .and. &
|
||||
micro_xs(mat % nuclide(l)) % fission > ZERO) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
cum_dsig = cum_dsig + dsigF * mat % atom_density(l)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig + dsig_f * mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
|
|
@ -3988,17 +3988,17 @@ contains
|
|||
+ cum_dsig / material_xs % fission)
|
||||
else if (materials(p % material) % id == deriv % diff_material &
|
||||
.and. material_xs % fission > ZERO) then
|
||||
dsigF = ZERO
|
||||
dsig_f = ZERO
|
||||
associate (nuc => nuclides(i_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsigF / micro_xs(i_nuclide) % fission)
|
||||
+ dsig_f / micro_xs(i_nuclide) % fission)
|
||||
else
|
||||
score = score * flux_deriv
|
||||
end if
|
||||
|
|
@ -4016,8 +4016,8 @@ contains
|
|||
p % last_E <= nuc % multipole % end_E .and. &
|
||||
micro_xs(mat % nuclide(l)) % nu_fission > ZERO) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
cum_dsig = cum_dsig + dsigF * mat % atom_density(l) &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig + dsig_f * mat % atom_density(l) &
|
||||
* micro_xs(mat % nuclide(l)) % nu_fission &
|
||||
/ micro_xs(mat % nuclide(l)) % fission
|
||||
end if
|
||||
|
|
@ -4028,17 +4028,17 @@ contains
|
|||
+ cum_dsig / material_xs % nu_fission)
|
||||
else if (materials(p % material) % id == deriv % diff_material &
|
||||
.and. material_xs % nu_fission > ZERO) then
|
||||
dsigF = ZERO
|
||||
dsig_f = ZERO
|
||||
associate (nuc => nuclides(i_nuclide))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsigF / micro_xs(i_nuclide) % fission)
|
||||
+ dsig_f / micro_xs(i_nuclide) % fission)
|
||||
else
|
||||
score = score * flux_deriv
|
||||
end if
|
||||
|
|
@ -4066,7 +4066,7 @@ contains
|
|||
real(8), intent(in) :: distance ! Neutron flight distance
|
||||
|
||||
integer :: i, l
|
||||
real(8) :: dsigT, dsigA, dsigF
|
||||
real(8) :: dsig_t, dsig_a, dsig_f
|
||||
|
||||
! A void material cannot be perturbed so it will not affect flux derivatives
|
||||
if (p % material == MATERIAL_VOID) return
|
||||
|
|
@ -4109,9 +4109,9 @@ contains
|
|||
! (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist
|
||||
! (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist
|
||||
call multipole_deriv_eval(nuc % multipole, p % E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
deriv % flux_deriv = deriv % flux_deriv &
|
||||
- distance * dsigT * mat % atom_density(l)
|
||||
- distance * dsig_t * mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
|
|
@ -4142,7 +4142,7 @@ contains
|
|||
type(Particle), intent(in) :: p
|
||||
|
||||
integer :: i, j, l
|
||||
real(8) :: dsigT, dsigA, dsigF
|
||||
real(8) :: dsig_t, dsig_a, dsig_f
|
||||
|
||||
! A void material cannot be perturbed so it will not affect flux derivatives
|
||||
if (p % material == MATERIAL_VOID) return
|
||||
|
|
@ -4196,8 +4196,8 @@ contains
|
|||
! (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
|
||||
! (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsigT, dsigA, dsigF)
|
||||
deriv % flux_deriv = deriv % flux_deriv + (dsigT - dsigA)&
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
deriv % flux_deriv = deriv % flux_deriv + (dsig_t - dsig_a)&
|
||||
/ (micro_xs(mat % nuclide(l)) % total &
|
||||
- micro_xs(mat % nuclide(l)) % absorption)
|
||||
! Note that this is an approximation! The real scattering
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue