diff --git a/src/constants.F90 b/src/constants.F90 index ab85b277a6..a3118f84da 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 diff --git a/src/multipole_header.F90 b/src/multipole_header.F90 index f2d0c91062..05eb62e1aa 100644 --- a/src/multipole_header.F90 +++ b/src/multipole_header.F90 @@ -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) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 8cd225bc98..a2e47fa434 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 diff --git a/src/physics.F90 b/src/physics.F90 index d5f6302a82..fa75fa0ed9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 2c660ce87b..2365ec46f2 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -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