From 85d8ff06ccffe93694b551671fcebb16d0ba966b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 28 Oct 2011 11:09:43 -0400 Subject: [PATCH] Incoherent inelastic S(a,b) complete. --- src/constants.f90 | 5 ++ src/cross_section.f90 | 3 + src/cross_section_header.f90 | 5 ++ src/physics.f90 | 147 ++++++++++++++++++++++++++++++++++- 4 files changed, 156 insertions(+), 4 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index 3ec09dc343..70017f7e9b 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -105,6 +105,11 @@ module constants & ANGLE_32_EQUI = 2, & ! 32 equiprobable bins & ANGLE_TABULAR = 3 ! Tabular angular distribution + ! Secondary energy mode for S(a,b) inelastic scattering + integer, parameter :: & + & SECONDARY_EQUAL = 0, & ! Equally-likely outgoing energy bins + & SECONDARY_SKEWED = 1 ! Skewed outgoing energy bins + ! Reaction types integer, parameter :: & TOTAL_XS = 1, & diff --git a/src/cross_section.f90 b/src/cross_section.f90 index 49906050b5..912d823e23 100644 --- a/src/cross_section.f90 +++ b/src/cross_section.f90 @@ -1205,6 +1205,9 @@ contains integer :: NMU ! number of outgoing angles integer :: JXS4 ! location of elastic energy table + ! read secondary energy mode for inelastic scattering + table % secondary_mode = NXS(7) + ! read number of inelastic energies and allocate arrays NE_in = XSS(JXS(1)) table % n_inelastic_e_in = NE_in diff --git a/src/cross_section_header.f90 b/src/cross_section_header.f90 index 4600d04a4a..5ec35d968f 100644 --- a/src/cross_section_header.f90 +++ b/src/cross_section_header.f90 @@ -132,6 +132,7 @@ module cross_section_header integer :: n_inelastic_e_in integer :: n_inelastic_e_out integer :: n_inelastic_mu + integer :: secondary_mode real(8), allocatable :: inelastic_e_in(:) real(8), allocatable :: inelastic_sigma(:) real(8), allocatable :: inelastic_e_out(:,:) @@ -178,6 +179,10 @@ module cross_section_header real(8) :: absorption real(8) :: fission real(8) :: nu_fission + + ! Information for S(a,b) use + logical :: use_sab + real(8) :: elastic_sab end type NuclideMicroXS !=============================================================================== diff --git a/src/physics.f90 b/src/physics.f90 index 09a73b5057..553cbc0e35 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -229,6 +229,10 @@ contains micro_xs(i) % index_grid = IE micro_xs(i) % interp_factor = f + ! Initialize sab treatment to false + micro_xs(i) % use_sab = .false. + micro_xs(i) % elastic_sab = ZERO + ! Initialize nuclide cross-sections to zero micro_xs(i) % fission = ZERO micro_xs(i) % nu_fission = ZERO @@ -261,6 +265,8 @@ contains ! then add back in the calculated S(a,b) elastic+inelastic cross section. if (index_sab > 0) then + micro_xs(i) % use_sab = .true. + ! Get pointer to S(a,b) table sab => sab_tables(index_sab) @@ -304,6 +310,9 @@ contains micro_xs(i) % total = micro_xs(i) % total - micro_xs(i) % elastic & + inelastic + elastic micro_xs(i) % elastic = inelastic + elastic + + ! Store ratio of elastic to elastic+inelastic for sampling later + micro_xs(i) % elastic_sab = elastic end if end subroutine calculate_nuclide_xs @@ -589,11 +598,17 @@ contains ! ======================================================================= ! ELASTIC SCATTERING - ! get pointer to elastic scattering reaction - rxn => nuc % reactions(1) + if (micro_xs(index_nuclide) % use_sab) then + ! S(a,b) scattering + call sab_scatter(p, index_nuclide, mat % sab_table) - ! Perform collision physics for elastic scattering - call elastic_scatter(p, nuc, rxn) + else + ! get pointer to elastic scattering reaction + rxn => nuc % reactions(1) + + ! Perform collision physics for elastic scattering + call elastic_scatter(p, nuc, rxn) + end if else ! ======================================================================= @@ -713,6 +728,130 @@ contains end subroutine elastic_scatter +!=============================================================================== +! SAB_SCATTER performs thermal scattering of a particle with a bound scatterer +! according to a specified S(a,b) table. Inelastic +! =============================================================================== + + subroutine sab_scatter(p, index_nuclide, index_sab) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index in micro_xs + integer, intent(in) :: index_sab ! index in sab_tables + + integer :: i ! incoming energy bin + integer :: j ! outgoing energy bin + integer :: k ! outgoing cosine bin + integer :: n_energy_out ! number of outgoing energy bins + real(8) :: f ! interpolation factor + real(8) :: r ! used for skewed sampling + real(8) :: E ! outgoing energy + real(8) :: E_ij ! outgoing energy j for E_in(i) + real(8) :: E_i1j ! outgoing energy j for E_in(i+1) + real(8) :: mu ! outgoing cosine + real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j) + real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j) + real(8) :: u, v, w ! directional cosines + character(MAX_LINE_LEN) :: msg + type(SAB_Table), pointer :: sab => null() + + ! Get pointer to S(a,b) table + sab => sab_tables(index_sab) + + ! Get index and interpolation factor for inelastic grid + if (p%E < sab % inelastic_e_in(1)) then + i = 1 + f = ZERO + else + i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) + f = (p%E - sab%inelastic_e_in(i)) / & + (sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i)) + end if + + ! Determine whether inelastic or elastic scattering will occur + if (rang() < micro_xs(index_nuclide) % elastic_sab / & + micro_xs(index_nuclide) % elastic) then + ! elastic scattering + + ! Outgoing energy is same as incoming energy + E = p % E + + msg = "Elastic scattering on S(a,b) table not yet supported" + call fatal_error(msg) + + else + ! Determine number of outgoing energy and angle bins + n_energy_out = sab % n_inelastic_e_out + + ! Now that we have an incoming energy bin, we need to determine the + ! outgoing energy bin. This will depend on the "secondary energy + ! mode". If the mode is 0, then the outgoing energy bin is chosen from a + ! set of equally-likely bins. However, if the mode is 1, then the first + ! two and last two bins are skewed to have lower probabilities than the + ! other bins (0.1 for the first and last bins and 0.4 for the second and + ! second to last bins, relative to a normal bin probability of 1) + + if (sab % secondary_mode == SECONDARY_EQUAL) then + ! All bins equally likely + j = 1 + rang() * n_energy_out + elseif (sab % secondary_mode == SECONDARY_SKEWED) then + r = rang() * (n_energy_out - 3) + if (r > ONE) then + ! equally likely N-4 middle bins + j = r + 2 + elseif (r > 0.6) then + ! second to last bin has relative probability of 0.4 + j = n_energy_out - 1 + elseif (r > 0.5) then + ! last bin has relative probability of 0.1 + j = n_energy_out + elseif (r > 0.1) then + ! second bin has relative probability of 0.4 + j = 2 + else + ! first bin has relative probability of 0.1 + j = 1 + end if + else + msg = "Invalid secondary energy mode on S(a,b) table " // & + trim(sab % name) + end if + + ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1) + E_ij = sab % inelastic_e_out(j,i) + E_i1j = sab % inelastic_e_out(j,i+1) + + ! Outgoing energy + E = (1 - f)*E_ij + f*E_i1j + + ! Sample outgoing cosine bin + k = 1 + rang() * sab % n_inelastic_mu + + ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) + mu_ijk = sab % inelastic_mu(k,j,i) + mu_i1jk = sab % inelastic_mu(k,j,i+1) + + ! Cosine of angle between incoming and outgoing neutron + mu = (1 - f)*mu_ijk + f*mu_i1jk + end if + + ! copy directional cosines + u = p % uvw(1) + v = p % uvw(2) + w = p % uvw(3) + + ! change direction of particle + call rotate_angle(u, v, w, mu) + p % uvw = (/ u, v, w /) + + ! change energy of particle + p % E = E + + ! Copy scattering cosine for tallies + p % mu = mu + + end subroutine sab_scatter + !=============================================================================== ! SAMPLE_TARGET_VELOCITY samples the target velocity based on the free gas ! scattering formulation used by most Monte Carlo codes. Excellent documentation