mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
update for new wmp library v1.0 - openmc code
This commit is contained in:
jingang
parent
3b8d1ffb4b
commit
60a88b9cab
5 changed files with 140 additions and 259 deletions
|
|
@ -25,7 +25,7 @@ module constants
|
|||
integer, parameter :: VERSION_VOLUME(2) = [1, 0]
|
||||
integer, parameter :: VERSION_VOXEL(2) = [1, 0]
|
||||
integer, parameter :: VERSION_MGXS_LIBRARY(2) = [1, 0]
|
||||
character(10), parameter :: VERSION_MULTIPOLE = "v0.2"
|
||||
character(10), parameter :: VERSION_MULTIPOLE = "v1.0"
|
||||
|
||||
! ============================================================================
|
||||
! ADJUSTABLE PARAMETERS
|
||||
|
|
|
|||
|
|
@ -10,27 +10,16 @@ module multipole_header
|
|||
!========================================================================
|
||||
! Multipole related constants
|
||||
|
||||
! Formalisms
|
||||
integer, parameter :: FORM_MLBW = 2, &
|
||||
FORM_RM = 3, &
|
||||
FORM_RML = 7
|
||||
|
||||
! Constants that determine which value to access
|
||||
integer, parameter :: MP_EA = 1 ! Pole
|
||||
|
||||
! Reich-Moore indices
|
||||
integer, parameter :: RM_RT = 2, & ! Residue total
|
||||
RM_RA = 3, & ! Residue absorption
|
||||
RM_RF = 4 ! Residue fission
|
||||
|
||||
! Multi-level Breit Wigner indices
|
||||
integer, parameter :: MLBW_RT = 2, & ! Residue total
|
||||
MLBW_RX = 3, & ! Residue compettitive
|
||||
MLBW_RA = 4, & ! Residue absorption
|
||||
MLBW_RF = 5 ! Residue fission
|
||||
! Residue indices
|
||||
integer, parameter :: MP_RS = 2, & ! Residue scattering
|
||||
MP_RA = 3, & ! Residue absorption
|
||||
MP_RF = 4 ! Residue fission
|
||||
|
||||
! Polynomial fit indices
|
||||
integer, parameter :: FIT_T = 1, & ! Total
|
||||
integer, parameter :: FIT_S = 1, & ! Scattering
|
||||
FIT_A = 2, & ! Absorption
|
||||
FIT_F = 3 ! Fission
|
||||
|
||||
|
|
@ -46,33 +35,22 @@ module multipole_header
|
|||
! Isotope Properties
|
||||
|
||||
logical :: fissionable ! Is this isotope fissionable?
|
||||
integer, allocatable :: l_value(:) ! The l index of the pole
|
||||
integer :: num_l ! Number of unique l values
|
||||
real(8), allocatable :: pseudo_k0RS(:) ! The value (sqrt(2*mass neutron
|
||||
! /reduced planck constant)
|
||||
! * AWR/(AWR + 1)
|
||||
! * scattering radius for
|
||||
! each l
|
||||
complex(8), allocatable :: data(:,:) ! Poles and residues
|
||||
real(8) :: sqrtAWR ! Square root of the atomic
|
||||
! weight ratio
|
||||
integer :: formalism ! R-matrix formalism
|
||||
|
||||
!=========================================================================
|
||||
! Windows
|
||||
|
||||
integer :: fit_order ! Order of the fit. 1 linear,
|
||||
! 2 quadratic, etc.
|
||||
real(8) :: start_E ! Start energy for the windows
|
||||
real(8) :: end_E ! End energy for the windows
|
||||
real(8) :: E_min ! Start energy for the windows
|
||||
real(8) :: E_max ! End energy for the windows
|
||||
real(8) :: spacing ! The actual spacing in sqrt(E)
|
||||
! space.
|
||||
! spacing = sqrt(multipole_w % endE - multipole_w % startE)
|
||||
! / multipole_w % windows
|
||||
integer, allocatable :: w_start(:) ! Contains the index of the pole at
|
||||
! the start of the window
|
||||
integer, allocatable :: w_end(:) ! Contains the index of the pole at
|
||||
! the end of the window
|
||||
integer, allocatable :: windows(:, :) ! Contains the indexes of the poles
|
||||
! at the start and end of the
|
||||
! window
|
||||
real(8), allocatable :: curvefit(:,:,:) ! Contains the fitting function.
|
||||
! (reaction type, coeff index,
|
||||
! window index)
|
||||
|
|
@ -96,7 +74,7 @@ contains
|
|||
character(len=*), intent(in) :: filename
|
||||
|
||||
character(len=10) :: version
|
||||
integer :: i, n_poles, n_residue_types, n_windows
|
||||
integer :: i, n_poles, n_residues, n_windows
|
||||
integer(HSIZE_T) :: dims_1d(1), dims_2d(2), dims_3d(3)
|
||||
integer(HID_T) :: file_id
|
||||
integer(HID_T) :: group_id
|
||||
|
|
@ -114,65 +92,49 @@ contains
|
|||
// trim(filename) // " uses version " // trim(version) // ".")
|
||||
|
||||
! Read scalar values.
|
||||
call read_dataset(this % formalism, group_id, "formalism")
|
||||
call read_dataset(this % spacing, group_id, "spacing")
|
||||
call read_dataset(this % sqrtAWR, group_id, "sqrtAWR")
|
||||
call read_dataset(this % start_E, group_id, "start_E")
|
||||
call read_dataset(this % end_E, group_id, "end_E")
|
||||
call read_dataset(this % E_min, group_id, "E_min")
|
||||
call read_dataset(this % E_max, group_id, "E_max")
|
||||
|
||||
! Read the "data" array. Use its shape to figure out the number of poles
|
||||
! and residue types in this data.
|
||||
dset = open_dataset(group_id, "data")
|
||||
call get_shape(dset, dims_2d)
|
||||
n_residue_types = int(dims_2d(1), 4) - 1
|
||||
n_residues = int(dims_2d(1), 4) - 1
|
||||
n_poles = int(dims_2d(2), 4)
|
||||
allocate(this % data(n_residue_types+1, n_poles))
|
||||
call read_dataset(this % data, dset)
|
||||
allocate(this % data(n_residues+1, n_poles))
|
||||
if (n_poles > 0) call read_dataset(this % data, dset)
|
||||
call close_dataset(dset)
|
||||
|
||||
! Check to see if this data includes fission residues.
|
||||
if (this % formalism == FORM_RM) then
|
||||
this % fissionable = (n_residue_types == 3)
|
||||
else
|
||||
! Assume FORM_MLBW.
|
||||
this % fissionable = (n_residue_types == 4)
|
||||
end if
|
||||
this % fissionable = (n_residues == 3)
|
||||
|
||||
! Read the "l_value" array.
|
||||
allocate(this % l_value(n_poles))
|
||||
call read_dataset(this % l_value, group_id, "l_value")
|
||||
|
||||
! Figure out the number of unique l values in the l_value array.
|
||||
do i = 1, n_poles
|
||||
if (.not. l_val_dict % has(this % l_value(i))) then
|
||||
call l_val_dict % set(this % l_value(i), 0)
|
||||
end if
|
||||
end do
|
||||
this % num_l = l_val_dict % size()
|
||||
call l_val_dict % clear()
|
||||
|
||||
! Read the "pseudo_K0RS" array.
|
||||
allocate(this % pseudo_k0RS(this % num_l))
|
||||
call read_dataset(this % pseudo_k0RS, group_id, "pseudo_K0RS")
|
||||
|
||||
! Read the "w_start" array and use its shape to figure out the number of
|
||||
! Read the "windows" array and use its shape to figure out the number of
|
||||
! windows.
|
||||
dset = open_dataset(group_id, "w_start")
|
||||
call get_shape(dset, dims_1d)
|
||||
n_windows = int(dims_1d(1), 4)
|
||||
allocate(this % w_start(n_windows))
|
||||
call read_dataset(this % w_start, dset)
|
||||
dset = open_dataset(group_id, "windows")
|
||||
call get_shape(dset, dims_2d)
|
||||
n_windows = int(dims_2d(2), 4)
|
||||
allocate(this % windows(dims_2d(1), n_windows))
|
||||
call read_dataset(this % windows, dset)
|
||||
call close_dataset(dset)
|
||||
|
||||
! Read the "w_end" and "broaden_poly" arrays.
|
||||
allocate(this % w_end(n_windows))
|
||||
call read_dataset(this % w_end, group_id, "w_end")
|
||||
! Read the "broaden_poly" arrays.
|
||||
dset = open_dataset(group_id, "broaden_poly")
|
||||
call get_shape(dset, dims_1d)
|
||||
if (dims_1d(1) /= n_windows) call fatal_error("broaden_poly array shape is&
|
||||
¬ consistent with the windows array shape in multipole library"&
|
||||
// trim(filename) // ".")
|
||||
allocate(this % broaden_poly(n_windows))
|
||||
call read_dataset(this % broaden_poly, group_id, "broaden_poly")
|
||||
call read_dataset(this % broaden_poly, dset)
|
||||
call close_dataset(dset)
|
||||
|
||||
! Read the "curvefit" array.
|
||||
dset = open_dataset(group_id, "curvefit")
|
||||
call get_shape(dset, dims_3d)
|
||||
if (dims_3d(3) /= n_windows) call fatal_error("curvefit array shape is not&
|
||||
&consistent with the windows array shape in multipole library"&
|
||||
// trim(filename) // ".")
|
||||
allocate(this % curvefit(dims_3d(1), dims_3d(2), dims_3d(3)))
|
||||
call read_dataset(this % curvefit, dset)
|
||||
call close_dataset(dset)
|
||||
|
|
|
|||
|
|
@ -12,11 +12,9 @@ module nuclide_header
|
|||
use hdf5_interface
|
||||
use math, only: faddeeva, w_derivative, &
|
||||
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
|
||||
use multipole_header, only: MP_EA, MP_RS, MP_RA, MP_RF, &
|
||||
FIT_S, FIT_A, FIT_F, MultipoleArray
|
||||
use message_passing
|
||||
use multipole_header, only: MultipoleArray
|
||||
use random_lcg, only: prn, future_prn, prn_set_stream
|
||||
use reaction_header, only: Reaction
|
||||
use sab_header, only: SAlphaBeta, sab_tables
|
||||
|
|
@ -849,7 +847,7 @@ contains
|
|||
integer :: threshold ! threshold energy index
|
||||
real(8) :: f ! interp factor on nuclide energy grid
|
||||
real(8) :: kT ! temperature in eV
|
||||
real(8) :: sig_t, sig_a, sig_f ! Intermediate multipole variables
|
||||
real(8) :: sig_s, sig_a, sig_f ! Intermediate multipole variables
|
||||
|
||||
! Initialize cached cross sections to zero
|
||||
micro_xs % elastic = CACHE_INVALID
|
||||
|
|
@ -859,8 +857,8 @@ contains
|
|||
! Check to see if there is multipole data present at this energy
|
||||
use_mp = .false.
|
||||
if (this % mp_present) then
|
||||
if (E >= this % multipole % start_E .and. &
|
||||
E <= this % multipole % end_E) then
|
||||
if (E >= this % multipole % E_min .and. &
|
||||
E <= this % multipole % E_max) then
|
||||
use_mp = .true.
|
||||
end if
|
||||
end if
|
||||
|
|
@ -868,9 +866,10 @@ contains
|
|||
! Evaluate multipole or interpolate
|
||||
if (use_mp) then
|
||||
! Call multipole kernel
|
||||
call multipole_eval(this % multipole, E, sqrtkT, sig_t, sig_a, sig_f)
|
||||
call multipole_eval(this % multipole, E, sqrtkT, sig_s, sig_a, sig_f)
|
||||
|
||||
micro_xs % total = sig_t
|
||||
micro_xs % total = sig_s + sig_a
|
||||
micro_xs % elastic = sig_s
|
||||
micro_xs % absorption = sig_a
|
||||
micro_xs % fission = sig_f
|
||||
|
||||
|
|
@ -1076,9 +1075,6 @@ contains
|
|||
micro_xs % elastic = (ONE - f) * rx % xs(i_temp, i_grid) + &
|
||||
f * rx % xs(i_temp, i_grid + 1)
|
||||
end associate
|
||||
else
|
||||
! For multipole, elastic is total - absorption
|
||||
micro_xs % elastic = micro_xs % total - micro_xs % absorption
|
||||
end if
|
||||
end subroutine nuclide_calculate_elastic_xs
|
||||
|
||||
|
|
@ -1130,7 +1126,7 @@ contains
|
|||
! sections in the resolved resonance regions
|
||||
!===============================================================================
|
||||
|
||||
subroutine multipole_eval(multipole, E, sqrtkT, sig_t, sig_a, sig_f)
|
||||
subroutine multipole_eval(multipole, E, sqrtkT, sig_s, 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
|
||||
|
|
@ -1138,7 +1134,7 @@ contains
|
|||
real(8), intent(in) :: sqrtkT ! The temperature in the form
|
||||
! sqrt(kT), at which
|
||||
! to evaluate the XS.
|
||||
real(8), intent(out) :: sig_t ! Total cross section
|
||||
real(8), intent(out) :: sig_s ! Scattering 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
|
||||
|
|
@ -1146,7 +1142,6 @@ contains
|
|||
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) :: 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
|
||||
|
|
@ -1166,18 +1161,13 @@ contains
|
|||
invE = ONE / E
|
||||
|
||||
! Locate us.
|
||||
i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing &
|
||||
i_window = floor((sqrtE - sqrt(multipole % E_min)) / multipole % spacing &
|
||||
+ ONE)
|
||||
startw = multipole % w_start(i_window)
|
||||
endw = multipole % w_end(i_window)
|
||||
|
||||
! Fill in factors.
|
||||
if (startw <= endw) then
|
||||
call compute_sig_t_factor(multipole, sqrtE, sig_t_factor)
|
||||
end if
|
||||
startw = multipole % windows(1, i_window)
|
||||
endw = multipole % windows(2, i_window)
|
||||
|
||||
! Initialize the ouptut cross sections.
|
||||
sig_t = ZERO
|
||||
sig_s = ZERO
|
||||
sig_a = ZERO
|
||||
sig_f = ZERO
|
||||
|
||||
|
|
@ -1190,7 +1180,7 @@ contains
|
|||
call broaden_wmp_polynomials(E, dopp, multipole % fit_order + 1, &
|
||||
broadened_polynomials)
|
||||
do i_poly = 1, multipole % fit_order+1
|
||||
sig_t = sig_t + multipole % curvefit(FIT_T, i_poly, i_window) &
|
||||
sig_s = sig_s + multipole % curvefit(FIT_S, i_poly, i_window) &
|
||||
* broadened_polynomials(i_poly)
|
||||
sig_a = sig_a + multipole % curvefit(FIT_A, i_poly, i_window) &
|
||||
* broadened_polynomials(i_poly)
|
||||
|
|
@ -1202,7 +1192,7 @@ contains
|
|||
else ! Evaluate as if it were a polynomial
|
||||
temp = invE
|
||||
do i_poly = 1, multipole % fit_order+1
|
||||
sig_t = sig_t + multipole % curvefit(FIT_T, i_poly, i_window) * temp
|
||||
sig_s = sig_s + multipole % curvefit(FIT_S, i_poly, i_window) * temp
|
||||
sig_a = sig_a + multipole % curvefit(FIT_A, i_poly, i_window) * temp
|
||||
if (multipole % fissionable) then
|
||||
sig_f = sig_f + multipole % curvefit(FIT_F, i_poly, i_window) * temp
|
||||
|
|
@ -1219,21 +1209,10 @@ contains
|
|||
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
|
||||
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
|
||||
sig_f = sig_f + real(multipole % data(MLBW_RF, i_pole) * c_temp)
|
||||
end if
|
||||
else if (multipole % formalism == FORM_RM) then
|
||||
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
|
||||
sig_f = sig_f + real(multipole % data(RM_RF, i_pole) * c_temp)
|
||||
end if
|
||||
sig_s = sig_s + real(multipole % data(MP_RS, i_pole) * c_temp)
|
||||
sig_a = sig_a + real(multipole % data(MP_RA, i_pole) * c_temp)
|
||||
if (multipole % fissionable) then
|
||||
sig_f = sig_f + real(multipole % data(MP_RF, i_pole) * c_temp)
|
||||
end if
|
||||
end do
|
||||
else
|
||||
|
|
@ -1243,21 +1222,10 @@ contains
|
|||
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
|
||||
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
|
||||
sig_f = sig_f + real(multipole % data(MLBW_RF, i_pole) * w_val)
|
||||
end if
|
||||
else if (multipole % formalism == FORM_RM) then
|
||||
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
|
||||
sig_f = sig_f + real(multipole % data(RM_RF, i_pole) * w_val)
|
||||
end if
|
||||
sig_s = sig_s + real(multipole % data(MP_RS, i_pole) * w_val)
|
||||
sig_a = sig_a + real(multipole % data(MP_RA, i_pole) * w_val)
|
||||
if (multipole % fissionable) then
|
||||
sig_f = sig_f + real(multipole % data(MP_RF, i_pole) * w_val)
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
|
|
@ -1270,7 +1238,7 @@ contains
|
|||
! temperature.
|
||||
!===============================================================================
|
||||
|
||||
subroutine multipole_deriv_eval(multipole, E, sqrtkT, sig_t, sig_a, sig_f)
|
||||
subroutine multipole_deriv_eval(multipole, E, sqrtkT, sig_s, 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
|
||||
|
|
@ -1278,12 +1246,11 @@ contains
|
|||
real(8), intent(in) :: sqrtkT ! The temperature in the form
|
||||
! sqrt(kT), at which to
|
||||
! evaluate the XS.
|
||||
real(8), intent(out) :: sig_t ! Total cross section
|
||||
real(8), intent(out) :: sig_s ! Scattering 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) :: 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)
|
||||
|
|
@ -1305,18 +1272,13 @@ contains
|
|||
&derivatives are not implemented for 0 Kelvin cross sections.")
|
||||
|
||||
! Locate us
|
||||
i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing &
|
||||
i_window = floor((sqrtE - sqrt(multipole % E_min)) / multipole % spacing &
|
||||
+ ONE)
|
||||
startw = multipole % w_start(i_window)
|
||||
endw = multipole % w_end(i_window)
|
||||
|
||||
! Fill in factors.
|
||||
if (startw <= endw) then
|
||||
call compute_sig_t_factor(multipole, sqrtE, sig_t_factor)
|
||||
end if
|
||||
startw = multipole % windows(1, i_window)
|
||||
endw = multipole % windows(2, i_window)
|
||||
|
||||
! Initialize the ouptut cross sections.
|
||||
sig_t = ZERO
|
||||
sig_s = ZERO
|
||||
sig_a = ZERO
|
||||
sig_f = ZERO
|
||||
|
||||
|
|
@ -1333,61 +1295,18 @@ contains
|
|||
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
|
||||
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
|
||||
sig_f = sig_f + real(multipole % data(MLBW_RF, i_pole) * w_val)
|
||||
end if
|
||||
else if (multipole % formalism == FORM_RM) then
|
||||
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
|
||||
sig_f = sig_f + real(multipole % data(RM_RF, i_pole) * w_val)
|
||||
end if
|
||||
sig_s = sig_s + real(multipole % data(MP_RS, i_pole) * w_val)
|
||||
sig_a = sig_a + real(multipole % data(MP_RA, i_pole) * w_val)
|
||||
if (multipole % fissionable) then
|
||||
sig_f = sig_f + real(multipole % data(MP_RF, i_pole) * w_val)
|
||||
end if
|
||||
end do
|
||||
sig_t = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_t
|
||||
sig_s = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sig_s
|
||||
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_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_sig_t_factor(multipole, sqrtE, sig_t_factor)
|
||||
type(MultipoleArray), intent(in) :: multipole
|
||||
real(8), intent(in) :: sqrtE
|
||||
complex(8), intent(out) :: sig_t_factor(multipole % num_l)
|
||||
|
||||
integer :: iL
|
||||
real(8) :: twophi(multipole % num_l)
|
||||
real(8) :: arg
|
||||
|
||||
do iL = 1, multipole % num_l
|
||||
twophi(iL) = multipole % pseudo_k0RS(iL) * sqrtE
|
||||
if (iL == 2) then
|
||||
twophi(iL) = twophi(iL) - atan(twophi(iL))
|
||||
else if (iL == 3) then
|
||||
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)
|
||||
twophi(iL) = twophi(iL) - atan(arg)
|
||||
end if
|
||||
end do
|
||||
|
||||
twophi = 2.0_8 * twophi
|
||||
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
|
||||
! for a given nuclide at the trial relative energy used in resonance scattering
|
||||
|
|
|
|||
|
|
@ -526,8 +526,8 @@ contains
|
|||
! Check to see if we are in a windowed multipole range. WMP only supports
|
||||
! the first fission reaction.
|
||||
if (nuc % mp_present) then
|
||||
if (E >= nuc % multipole % start_E .and. &
|
||||
E <= nuc % multipole % end_E) then
|
||||
if (E >= nuc % multipole % E_min .and. &
|
||||
E <= nuc % multipole % E_max) then
|
||||
i_reaction = nuc % index_fission(1)
|
||||
return
|
||||
end if
|
||||
|
|
|
|||
|
|
@ -3009,7 +3009,7 @@ contains
|
|||
integer :: l
|
||||
logical :: scoring_diff_nuclide
|
||||
real(8) :: flux_deriv
|
||||
real(8) :: dsig_t, dsig_a, dsig_f, cum_dsig
|
||||
real(8) :: dsig_s, dsig_a, dsig_f, cum_dsig
|
||||
|
||||
if (score == ZERO) return
|
||||
|
||||
|
|
@ -3250,17 +3250,18 @@ contains
|
|||
if (mat % nuclide(l) == p % event_nuclide) exit
|
||||
end do
|
||||
|
||||
dsig_t = ZERO
|
||||
dsig_s = 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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsig_t * mat % atom_density(l) / material_xs % total)
|
||||
+ (dsig_s + dsig_a) * mat % atom_density(l) &
|
||||
/ material_xs % total)
|
||||
end associate
|
||||
else
|
||||
score = score * flux_deriv
|
||||
|
|
@ -3276,18 +3277,17 @@ contains
|
|||
if (mat % nuclide(l) == p % event_nuclide) exit
|
||||
end do
|
||||
|
||||
dsig_t = ZERO
|
||||
dsig_s = 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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv + (dsig_t - dsig_a) &
|
||||
* mat % atom_density(l) / &
|
||||
score = score * (flux_deriv + dsig_s * mat % atom_density(l) / &
|
||||
(material_xs % total - material_xs % absorption))
|
||||
end associate
|
||||
else
|
||||
|
|
@ -3306,10 +3306,10 @@ contains
|
|||
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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv + dsig_a * mat % atom_density(l) &
|
||||
|
|
@ -3331,10 +3331,10 @@ contains
|
|||
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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
|
|
@ -3356,10 +3356,10 @@ contains
|
|||
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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
|
|
@ -3392,12 +3392,13 @@ contains
|
|||
do l = 1, mat % n_nuclides
|
||||
associate (nuc => nuclides(mat % nuclide(l)))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E .and. &
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max .and. &
|
||||
micro_xs(mat % nuclide(l)) % total > ZERO) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig + dsig_t * mat % atom_density(l)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig + (dsig_s + dsig_a) &
|
||||
* mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
|
|
@ -3406,17 +3407,17 @@ contains
|
|||
+ cum_dsig / material_xs % total)
|
||||
else if (materials(p % material) % id == deriv % diff_material &
|
||||
.and. material_xs % total > ZERO) then
|
||||
dsig_t = ZERO
|
||||
dsig_s = 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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
+ dsig_t / micro_xs(i_nuclide) % total)
|
||||
+ (dsig_s + dsig_a) / micro_xs(i_nuclide) % total)
|
||||
else
|
||||
score = score * flux_deriv
|
||||
end if
|
||||
|
|
@ -3430,14 +3431,13 @@ contains
|
|||
do l = 1, mat % n_nuclides
|
||||
associate (nuc => nuclides(mat % nuclide(l)))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E .and. &
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max .and. &
|
||||
(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, dsig_t, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig &
|
||||
+ (dsig_t - dsig_a) * mat % atom_density(l)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig + dsig_s * mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
|
|
@ -3447,17 +3447,17 @@ contains
|
|||
else if ( materials(p % material) % id == deriv % diff_material &
|
||||
.and. (material_xs % total - material_xs % absorption) > ZERO)&
|
||||
then
|
||||
dsig_t = ZERO
|
||||
dsig_s = 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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv + (dsig_t - dsig_a) &
|
||||
score = score * (flux_deriv + dsig_s &
|
||||
/ (micro_xs(i_nuclide) % total &
|
||||
- micro_xs(i_nuclide) % absorption))
|
||||
else
|
||||
|
|
@ -3473,11 +3473,11 @@ contains
|
|||
do l = 1, mat % n_nuclides
|
||||
associate (nuc => nuclides(mat % nuclide(l)))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E .and. &
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max .and. &
|
||||
micro_xs(mat % nuclide(l)) % absorption > ZERO) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig + dsig_a * mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
|
|
@ -3490,10 +3490,10 @@ contains
|
|||
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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
|
|
@ -3511,11 +3511,11 @@ contains
|
|||
do l = 1, mat % n_nuclides
|
||||
associate (nuc => nuclides(mat % nuclide(l)))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E .and. &
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max .and. &
|
||||
micro_xs(mat % nuclide(l)) % fission > ZERO) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
cum_dsig = cum_dsig + dsig_f * mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
|
|
@ -3528,10 +3528,10 @@ contains
|
|||
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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
|
|
@ -3549,11 +3549,11 @@ contains
|
|||
do l = 1, mat % n_nuclides
|
||||
associate (nuc => nuclides(mat % nuclide(l)))
|
||||
if (nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E .and. &
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max .and. &
|
||||
micro_xs(mat % nuclide(l)) % nu_fission > ZERO) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, 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
|
||||
|
|
@ -3568,10 +3568,10 @@ contains
|
|||
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
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
call multipole_deriv_eval(nuc % multipole, p % last_E, &
|
||||
p % sqrtkT, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
end if
|
||||
end associate
|
||||
score = score * (flux_deriv &
|
||||
|
|
@ -3603,7 +3603,7 @@ contains
|
|||
real(8), intent(in) :: distance ! Neutron flight distance
|
||||
|
||||
integer :: i, l
|
||||
real(8) :: dsig_t, dsig_a, dsig_f
|
||||
real(8) :: dsig_s, dsig_a, dsig_f
|
||||
|
||||
! A void material cannot be perturbed so it will not affect flux derivatives
|
||||
if (p % material == MATERIAL_VOID) return
|
||||
|
|
@ -3640,15 +3640,15 @@ contains
|
|||
do l=1, mat % n_nuclides
|
||||
associate (nuc => nuclides(mat % nuclide(l)))
|
||||
if (nuc % mp_present .and. &
|
||||
p % E >= nuc % multipole % start_E .and. &
|
||||
p % E <= nuc % multipole % end_E) then
|
||||
p % E >= nuc % multipole % E_min .and. &
|
||||
p % E <= nuc % multipole % E_max) then
|
||||
! phi is proportional to e^(-Sigma_tot * dist)
|
||||
! (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, dsig_t, dsig_a, dsig_f)
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
deriv % flux_deriv = deriv % flux_deriv &
|
||||
- distance * dsig_t * mat % atom_density(l)
|
||||
- distance * (dsig_s + dsig_a) * mat % atom_density(l)
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
|
|
@ -3679,7 +3679,7 @@ contains
|
|||
type(Particle), intent(in) :: p
|
||||
|
||||
integer :: i, j, l
|
||||
real(8) :: dsig_t, dsig_a, dsig_f
|
||||
real(8) :: dsig_s, dsig_a, dsig_f
|
||||
|
||||
! A void material cannot be perturbed so it will not affect flux derivatives
|
||||
if (p % material == MATERIAL_VOID) return
|
||||
|
|
@ -3727,14 +3727,14 @@ contains
|
|||
associate (nuc => nuclides(mat % nuclide(l)))
|
||||
if (mat % nuclide(l) == p % event_nuclide .and. &
|
||||
nuc % mp_present .and. &
|
||||
p % last_E >= nuc % multipole % start_E .and. &
|
||||
p % last_E <= nuc % multipole % end_E) then
|
||||
p % last_E >= nuc % multipole % E_min .and. &
|
||||
p % last_E <= nuc % multipole % E_max) then
|
||||
! phi is proportional to Sigma_s
|
||||
! (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, dsig_t, dsig_a, dsig_f)
|
||||
deriv % flux_deriv = deriv % flux_deriv + (dsig_t - dsig_a)&
|
||||
p % sqrtkT, dsig_s, dsig_a, dsig_f)
|
||||
deriv % flux_deriv = deriv % flux_deriv + (dsig_s + 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