mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Improve multipole readability
This commit is contained in:
Sterling Harper
parent
83773b0321
commit
afccdba6d5
3 changed files with 88 additions and 88 deletions
|
|
@ -6,7 +6,7 @@ module cross_section
|
|||
use global
|
||||
use list_header, only: ListElemInt
|
||||
use material_header, only: Material
|
||||
use math, only: w, broaden_n_polynomials
|
||||
use math, only: faddeeva, broaden_wmp_polynomials
|
||||
use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, RM_RF, &
|
||||
MLBW_RT, MLBW_RX, MLBW_RA, MLBW_RF, FIT_T, FIT_A,&
|
||||
FIT_F, MultipoleArray, max_poly, max_L, max_poles
|
||||
|
|
@ -588,123 +588,119 @@ contains
|
|||
real(8), intent(out) :: sigT ! Total cross section
|
||||
real(8), intent(out) :: sigA ! Absorption cross section
|
||||
real(8), intent(out) :: sigF ! Fission cross section
|
||||
complex(8) :: psi_ki ! The value of the psi-ki function for the asymptotic
|
||||
! form
|
||||
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)
|
||||
real(8) :: sqrtE ! sqrt(E), eV
|
||||
real(8) :: invE ! 1/E, eV
|
||||
real(8) :: dopp ! sqrt(atomic weight ratio / kT)
|
||||
real(8) :: dopp_ecoef ! sqrt(atomic weight ratio * pi / kT) / E
|
||||
real(8) :: dopp ! sqrt(atomic weight ratio / kT) = 1 / (2 sqrt(xi))
|
||||
real(8) :: temp ! real temporary value
|
||||
real(8) :: E ! energy, eV
|
||||
real(8) :: sqrtkT ! sqrt(kT (in eV))
|
||||
integer :: iP ! index of pole
|
||||
integer :: iC ! index of curvefit
|
||||
integer :: iW ! index of window
|
||||
integer :: i_pole ! index of pole
|
||||
integer :: i_poly ! index of curvefit
|
||||
integer :: i_window ! index of window
|
||||
integer :: startw ! window start pointer (for poles)
|
||||
integer :: startw_1 ! window start pointer - 1
|
||||
integer :: startw_endw ! window start pointer - window end pointer
|
||||
integer :: endw ! window end pointer
|
||||
|
||||
! Convert to eV
|
||||
! ==========================================================================
|
||||
! Bookkeeping
|
||||
|
||||
! Convert to eV.
|
||||
E = Emev * 1.0e6_8
|
||||
sqrtkT = sqrtkT_ * 1.0e3_8
|
||||
|
||||
! Define some frequently used variables.
|
||||
sqrtE = sqrt(E)
|
||||
invE = ONE/E
|
||||
invE = ONE / E
|
||||
dopp = multipole % sqrtAWR / sqrtkT
|
||||
|
||||
if(.not. mp_already_alloc) then
|
||||
if (.not. mp_already_alloc) then
|
||||
call multipole_eval_allocate()
|
||||
end if
|
||||
|
||||
! Locate us
|
||||
iW = floor((sqrtE - sqrt(multipole % start_E))/multipole % spacing + ONE)
|
||||
! Locate us.
|
||||
i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing &
|
||||
+ ONE)
|
||||
startw = multipole % w_start(i_window)
|
||||
endw = multipole % w_end(i_window)
|
||||
|
||||
startw = multipole % w_start(iW)
|
||||
startw_1 = startw - 1 ! This is an index shift parameter.
|
||||
endw = multipole % w_end(iW)
|
||||
startw_endw = endw - startw + 1
|
||||
|
||||
! Fill in factors
|
||||
! Fill in factors.
|
||||
if (startw <= endw) then
|
||||
call fill_factors(multipole, sqrtE, sigT_factor, twophi, multipole % num_l)
|
||||
call fill_factors(multipole, sqrtE, sigT_factor, twophi, &
|
||||
multipole % num_l)
|
||||
end if
|
||||
|
||||
! Generate some doppler broadening parameters
|
||||
|
||||
! dopp_ecoef is inverse of dopp, divided by E, multiplied by sqrt(pi).
|
||||
dopp = multipole % sqrtAWR / sqrtKT
|
||||
dopp_ecoef = dopp * invE * SQRT_PI
|
||||
|
||||
! Initialize the ouptut cross sections.
|
||||
sigT = ZERO
|
||||
sigA = ZERO
|
||||
sigF = ZERO
|
||||
! Evaluate linefit first
|
||||
|
||||
if(sqrtkT /= 0 .and. multipole % broaden_poly(iW) == 1) then ! Broaden the curvefit.
|
||||
call broaden_n_polynomials(E, dopp, multipole % fit_order + 1, broadened_polynomials)
|
||||
! ==========================================================================
|
||||
! Add the contribution from the curvefit polynomial.
|
||||
|
||||
do iC = 1, multipole % fit_order+1
|
||||
sigT = sigT + multipole % curvefit(FIT_T, iC, iW)*broadened_polynomials(iC)
|
||||
sigA = sigA + multipole % curvefit(FIT_A, iC, iW)*broadened_polynomials(iC)
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + multipole % curvefit(FIT_F, iC, iW)*broadened_polynomials(iC)
|
||||
end if
|
||||
if (sqrtkT /= ZERO .and. multipole % broaden_poly(i_window) == 1) then
|
||||
! Broaden the curvefit.
|
||||
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) &
|
||||
* 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)
|
||||
end do
|
||||
else ! Evaluate as if it were a polynomial
|
||||
temp = invE
|
||||
do iC = 1, multipole % fit_order+1
|
||||
|
||||
sigT = sigT + multipole % curvefit(FIT_T, iC, iW)*temp
|
||||
sigA = sigA + multipole % curvefit(FIT_A, iC, iW)*temp
|
||||
if (multipole % fissionable) then
|
||||
sigF = sigF + multipole % curvefit(FIT_F, iC, iW)*temp
|
||||
end if
|
||||
|
||||
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
|
||||
temp = temp * sqrtE
|
||||
end do
|
||||
end if
|
||||
|
||||
! Then get the poles we want and broaden them.
|
||||
! ==========================================================================
|
||||
! Add the contribution from the poles in this window.
|
||||
|
||||
if (sqrtkT == ZERO) then
|
||||
! If at 0K, use asymptotic form.
|
||||
do iP = startw, endw
|
||||
psi_ki = -ONEI/(multipole % data(MP_EA, iP) - sqrtE)
|
||||
c_temp = psi_ki/E
|
||||
do i_pole = startw, endw
|
||||
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, iP) * c_temp * &
|
||||
sigT_factor(multipole % l_value(iP))) &
|
||||
+ real(multipole % data(MLBW_RX, iP) * c_temp)
|
||||
sigA = sigA + real(multipole % data(MLBW_RA, iP) * c_temp)
|
||||
sigF = sigF + real(multipole % data(MLBW_RF, iP) * c_temp)
|
||||
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)
|
||||
sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * c_temp)
|
||||
else if (multipole % formalism == FORM_RM) then
|
||||
sigT = sigT + real(multipole % data(RM_RT, iP) * c_temp* &
|
||||
sigT_factor(multipole % l_value(iP)))
|
||||
sigA = sigA + real(multipole % data(RM_RA, iP) * c_temp)
|
||||
sigF = sigF + real(multipole % data(RM_RF, iP) * c_temp)
|
||||
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)
|
||||
end if
|
||||
end do
|
||||
else
|
||||
! At temperature, use Faddeeva function-based form.
|
||||
if(endw >= startw) then
|
||||
do iP = startw, endw
|
||||
Z = (sqrtE - multipole % data(MP_EA, iP)) * dopp
|
||||
w_val = w(Z) * dopp_ecoef
|
||||
|
||||
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, iP) * &
|
||||
sigT_factor(multipole%l_value(iP)) + &
|
||||
multipole % data(MLBW_RX, iP)) * w_val)
|
||||
sigA = sigA + real(multipole % data(MLBW_RA, iP) * w_val)
|
||||
sigF = sigF + real(multipole % data(MLBW_RF, iP) * w_val)
|
||||
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)
|
||||
sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val)
|
||||
else if (multipole % formalism == FORM_RM) then
|
||||
sigT = sigT + real(multipole % data(RM_RT, iP) * w_val * &
|
||||
sigT_factor(multipole % l_value(iP)))
|
||||
sigA = sigA + real(multipole % data(RM_RA, iP) * w_val)
|
||||
sigF = sigF + real(multipole % data(RM_RF, iP) * w_val)
|
||||
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)
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
|
|
@ -735,13 +731,14 @@ contains
|
|||
arg = 3.0_8 * twophi(iL) / (3.0_8 - twophi(iL)**2)
|
||||
twophi(iL) = twophi(iL) - atan(arg)
|
||||
else if (iL == 4) then
|
||||
arg = twophi(iL) * (15.0_8 - twophi(iL)**2) / (15.0_8 - 6.0_8 * twophi(iL)**2)
|
||||
arg = twophi(iL) * (15.0_8 - twophi(iL)**2) &
|
||||
/ (15.0_8 - 6.0_8 * twophi(iL)**2)
|
||||
twophi(iL) = twophi(iL) - atan(arg)
|
||||
end if
|
||||
end do
|
||||
|
||||
twophi = 2.0_8 * twophi
|
||||
sigT_factor = cmplx(cos(twophi),-sin(twophi), KIND=8)
|
||||
sigT_factor = cmplx(cos(twophi), -sin(twophi), KIND=8)
|
||||
end subroutine
|
||||
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -1178,7 +1178,7 @@ contains
|
|||
do k = 2, mat % n_nuclides
|
||||
! Warn the user if the nuclides don't have identical temperatues.
|
||||
if (nuclides(mat % nuclide(k)) % kT /= min_temp &
|
||||
.and. .not. warning_given) then
|
||||
.and. .not. warning_given .and. multipole_active) then
|
||||
call warning("OpenMC cannot &
|
||||
&identify the temperature of at least one cell. For the &
|
||||
&purposes of multipole cross section evaluations, all cells &
|
||||
|
|
|
|||
29
src/math.F90
29
src/math.F90
|
|
@ -719,13 +719,14 @@ contains
|
|||
end function watt_spectrum
|
||||
|
||||
!===============================================================================
|
||||
! W acts as a front end to the MIT Faddeeva function, Faddeeva_w.
|
||||
! FADDEEVA the Faddeeva function, using Stephen Johnson's implementation
|
||||
!===============================================================================
|
||||
|
||||
function w(z) result(wv)
|
||||
function faddeeva(z) result(wv)
|
||||
complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at
|
||||
complex(8) :: wv ! The resulting w(z) value
|
||||
real(C_DOUBLE) :: relerr ! Target relative error in inner loop of MIT Faddeeva
|
||||
real(C_DOUBLE) :: relerr ! Target relative error in inner loop of MIT
|
||||
! Faddeeva
|
||||
|
||||
! Technically, the value we want is given by the equation:
|
||||
! w(z) = I/Pi * Integrate[Exp[-t^2]/(z-t), {t, -Infinity, Infinity}]
|
||||
|
|
@ -747,23 +748,24 @@ contains
|
|||
wv = -conjg(faddeeva_w(conjg(z), relerr))
|
||||
end if
|
||||
|
||||
end function w
|
||||
end function faddeeva
|
||||
|
||||
!===============================================================================
|
||||
! BROADEN_N_POLYNOMIALS doppler broadens polynomials of the form
|
||||
! BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit. The
|
||||
! curvefit is a polynomial of the form
|
||||
! a/En + b/sqrt(En) + c + d sqrt(En) ...
|
||||
! exactly and quickly.
|
||||
!===============================================================================
|
||||
|
||||
subroutine broaden_n_polynomials(En, dopp, n, factors)
|
||||
subroutine broaden_wmp_polynomials(En, dopp, n, factors)
|
||||
real(8), intent(in) :: En ! Energy to evaluate at
|
||||
real(8), intent(in) :: dopp ! sqrt(atomic weight ratio / kT), kT given in eV.
|
||||
real(8), intent(in) :: dopp ! sqrt(atomic weight ratio / kT),
|
||||
! kT given in eV.
|
||||
integer, intent(in) :: n ! number of components to polynomial
|
||||
real(8), intent(out):: factors(n) ! output leading coefficient
|
||||
|
||||
integer :: i
|
||||
|
||||
real(8) :: sqrtE ! Sqrt(energy)
|
||||
real(8) :: sqrtE ! sqrt(energy)
|
||||
real(8) :: beta ! sqrt(atomic weight ratio * E / kT)
|
||||
real(8) :: half_inv_dopp2 ! 0.5 / dopp**2
|
||||
real(8) :: quarter_inv_dopp4 ! 0.25 / dopp**4
|
||||
|
|
@ -790,20 +792,21 @@ contains
|
|||
|
||||
factors(1) = erfbeta / En
|
||||
factors(2) = ONE / sqrtE
|
||||
factors(3) = factors(1) * (half_inv_dopp2 + En) + exp_m_beta2 / (beta * SQRT_PI)
|
||||
factors(3) = factors(1) * (half_inv_dopp2 + En) &
|
||||
+ exp_m_beta2 / (beta * SQRT_PI)
|
||||
|
||||
! Perform recursive broadening of high order components
|
||||
do i = 1, n-3
|
||||
if (i /= 1) then
|
||||
factors(i+3) = -factors(i-1) * (i - ONE) * i * quarter_inv_dopp4 + &
|
||||
factors(i+1) * (En + (ONE + TWO * i) * half_inv_dopp2)
|
||||
factors(i+3) = -factors(i-1) * (i - ONE) * i * quarter_inv_dopp4 &
|
||||
+ factors(i+1) * (En + (ONE + TWO * 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+3) = factors(i+1)*(En + (ONE + TWO * i) * half_inv_dopp2)
|
||||
end if
|
||||
end do
|
||||
end subroutine broaden_n_polynomials
|
||||
end subroutine broaden_wmp_polynomials
|
||||
|
||||
!===============================================================================
|
||||
! find_angle finds the closest angle on the data grid and returns that index
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue