From 500843bce431d9efaf9e4ef275501a6535024fc8 Mon Sep 17 00:00:00 2001 From: walshjon Date: Thu, 19 Mar 2015 03:19:04 -0700 Subject: [PATCH 1/2] user option to enforce LAB isotropic elastic, disable S(a,b) --- src/input_xml.F90 | 23 +++- src/list_header.F90 | 247 ++++++++++++++++++++++++++++++++++++++++ src/material_header.F90 | 3 + src/physics.F90 | 52 ++++++--- 4 files changed, 306 insertions(+), 19 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4c1f3c244c..9b6259a57c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -7,7 +7,7 @@ module input_xml use error, only: fatal_error, warning use geometry_header, only: Cell, Surface, Lattice, RectLattice, HexLattice use global - use list_header, only: ListChar, ListReal + use list_header, only: ListChar, ListLog, ListReal use mesh_header, only: StructuredMesh use output, only: write_message use plot_header @@ -1577,6 +1577,7 @@ contains character(MAX_LINE_LEN) :: temp_str ! temporary string when reading type(ListChar) :: list_names ! temporary list of nuclide names type(ListReal) :: list_density ! temporary list of nuclide densities + type(ListLog) :: list_iso_lab ! temporary list of isotropic lab scatterers type(Material), pointer :: mat => null() type(Node), pointer :: doc => null() type(Node), pointer :: node_mat => null() @@ -1731,6 +1732,20 @@ contains end if end if + ! Check enforced isotropic lab scattering + if (check_for_node(node_nuc, "lab")) then + call get_node_value(node_nuc, "lab", temp_str) + if (trim(adjustl(to_lower(temp_str))) == "true") then + call list_iso_lab % append(.true.) + else if (trim(adjustl(to_lower(temp_str))) == "false") then + call list_iso_lab % append(.false.) + else + call fatal_error("Isotropic lab scattering must be true or false") + end if + else + call list_iso_lab % append(.false.) + end if + ! store full name call get_node_value(node_nuc, "name", temp_str) if (check_for_node(node_nuc, "xs")) & @@ -1822,6 +1837,7 @@ contains allocate(mat % names(n)) allocate(mat % nuclide(n)) allocate(mat % atom_density(n)) + allocate(mat % p0(n)) ALL_NUCLIDES: do j = 1, mat % n_nuclides ! Check that this nuclide is listed in the cross_sections.xml file @@ -1858,6 +1874,10 @@ contains ! Copy name and atom/weight percent mat % names(j) = name mat % atom_density(j) = list_density % get_item(j) + + ! Copy isotropic lab scattering flag + mat % p0(j) = list_iso_lab % get_item(j) + end do ALL_NUCLIDES ! Check to make sure either all atom percents or all weight percents are @@ -1874,6 +1894,7 @@ contains ! Clear lists call list_names % clear() call list_density % clear() + call list_iso_lab % clear() ! ======================================================================= ! READ AND PARSE TAG FOR S(a,b) DATA diff --git a/src/list_header.F90 b/src/list_header.F90 index a9311ec64a..e4cc56ed0a 100644 --- a/src/list_header.F90 +++ b/src/list_header.F90 @@ -34,6 +34,13 @@ module list_header type(ListElemChar), pointer :: prev => null() end type ListElemChar + type :: ListElemLog + logical :: data + type(ListElemLog), pointer :: next => null() + type(ListElemLog), pointer :: prev => null() + end type ListElemLog + + !=============================================================================== ! LIST* types contain the linked list with convenience methods. We originally ! considered using unlimited polymorphism to provide a single type, but compiler @@ -107,6 +114,28 @@ module list_header procedure :: size => list_size_char ! Size of list end type ListChar + type, public :: ListLog + private + integer :: count = 0 ! Number of elements in list + + ! Used in get_item for fast sequential lookups + integer :: last_index = huge(0) + type(ListElemLog), pointer :: last_elem => null() + + ! Pointers to beginning and end of list + type(ListElemLog), public, pointer :: head => null() + type(ListElemLog), public, pointer :: tail => null() + contains + procedure :: append => list_append_log ! Add item to end of list + procedure :: clear => list_clear_log ! Remove all items + procedure :: contains => list_contains_log ! Does list contain? + procedure :: get_item => list_get_item_log ! Get i-th item in list + procedure :: index => list_index_log ! Determine index of given item + procedure :: insert => list_insert_log ! Insert item in i-th position + procedure :: remove => list_remove_log ! Remove specified item + procedure :: size => list_size_log ! Size of list + end type ListLog + contains !=============================================================================== @@ -189,6 +218,31 @@ contains end subroutine list_append_char + subroutine list_append_log(this, data) + class(ListLog) :: this + logical :: data + + type(ListElemLog), pointer :: elem + + ! Create element and set dat + allocate(elem) + elem % data = data + + if (.not. associated(this % head)) then + ! If list is empty, set head and tail to new element + this % head => elem + this % tail => elem + else + ! Otherwise append element at end of list + this % tail % next => elem + elem % prev => this % tail + this % tail => this % tail % next + end if + + this % count = this % count + 1 + + end subroutine list_append_log + !=============================================================================== ! LIST_CLEAR removes all elements from the list !=============================================================================== @@ -271,6 +325,32 @@ contains end subroutine list_clear_char + subroutine list_clear_log(this) + class(ListLog) :: this + + type(ListElemLog), pointer :: current => null() + type(ListElemLog), pointer :: next => null() + + if (this % count > 0) then + current => this % head + do while (associated(current)) + ! Set pointer to next element + next => current % next + + ! Deallocate memory for current element + deallocate(current) + + ! Move to next element + current => next + end do + + nullify(this % head) + nullify(this % tail) + this % count = 0 + end if + + end subroutine list_clear_log + !=============================================================================== ! LIST_CONTAINS determines whether the list contains a specified item. Since it ! relies on the index method, it is O(n). @@ -303,6 +383,15 @@ contains end function list_contains_char + function list_contains_log(this, data) result(in_list) + class(ListLog) :: this + logical :: data + logical :: in_list + + in_list = (this % index(data) > 0) + + end function list_contains_log + !=============================================================================== ! LIST_GET_ITEM returns the item in the list at position 'i_list'. If the index ! is out of bounds, an error code is returned. @@ -407,6 +496,39 @@ contains end function list_get_item_char + function list_get_item_log(this, i_list) result(data) + class(ListLog) :: this + integer :: i_list + logical :: data + + integer :: last_index + + if (i_list < 1 .or. i_list > this % count) then + ! Check for index out of bounds + data = .false. + elseif (i_list == 1) then + data = this % head % data + this % last_index = 1 + this % last_elem => this % head + elseif (i_list == this % count) then + data = this % tail % data + this % last_index = this % count + this % last_elem => this % tail + else + if (i_list < this % last_index) then + this % last_index = 1 + this % last_elem => this % head + end if + + do last_index = this % last_index + 1, i_list + this % last_elem => this % last_elem % next + this % last_index = last_index + end do + data = this % last_elem % data + end if + + end function list_get_item_log + !=============================================================================== ! LIST_INDEX determines the first index in the list that contains 'data'. If ! 'data' is not present in the list, the return value is -1. @@ -475,6 +597,27 @@ contains end function list_index_char + function list_index_log(this, data) result(i_list) + + class(ListLog) :: this + logical :: data + integer :: i_list + + type(ListElemLog), pointer :: elem + + i_list = 0 + elem => this % head + do while (associated(elem)) + i_list = i_list + 1 + if (data .eqv. elem % data) exit + elem => elem % next + end do + + ! Check if we reached the end of the list + if (.not. associated(elem)) i_list = -1 + + end function list_index_log + !=============================================================================== ! LIST_INSERT inserts 'data' at index 'i_list' within the list. If 'i_list' ! exceeds the size of the list, the data is appends at the end of the list. @@ -646,6 +789,62 @@ contains end subroutine list_insert_char + subroutine list_insert_log(this, i_list, data) + + class(ListLog) :: this + integer :: i_list + logical :: data + + integer :: i + type(ListElemLog), pointer :: elem => null() + type(ListElemLog), pointer :: new_elem => null() + + if (i_list > this % count) then + ! Check whether specified index is greater than number of elements -- if + ! so, just append it to the end of the list + call this % append(data) + + else if (i_list == 1) then + ! Check for new head element + allocate(new_elem) + new_elem % data = data + new_elem % next => this % head + this % head => new_elem + this % count = this % count + 1 + + else + ! Default case with new element somewhere in middle of list + if (i_list >= this % last_index) then + i = this % last_index + elem => this % last_elem + else + i = 0 + elem => this % head + end if + do while (associated(elem)) + i = i + 1 + if (i == i_list - 1) then + ! Allocate new element + allocate(new_elem) + new_elem % data = data + + ! Put it before the i-th element + new_elem % prev => elem % prev + new_elem % next => elem + new_elem % prev % next => new_elem + new_elem % next % prev => new_elem + this % count = this % count + 1 + this % last_index = i_list + this % last_elem => new_elem + exit + end if + i = i + 1 + elem => elem % next + end do + end if + + end subroutine list_insert_log + !=============================================================================== ! LIST_REMOVE removes the first item in the list that contains 'data'. If 'data' ! is not in the list, no action is taken. @@ -768,6 +967,45 @@ contains end subroutine list_remove_char + subroutine list_remove_log(this, data) + + class(ListLog) :: this + logical :: data + + type(ListElemLog), pointer :: elem => null() + + elem => this % head + do while (associated(elem)) + ! Check for matching data + if (elem % data .eqv. data) then + + ! Determine whether the current element is the head, tail, or a middle + ! element + if (associated(elem, this % head)) then + this % head => elem % next + if (associated(elem, this % tail)) nullify(this % tail) + if (associated(this % head)) nullify(this % head % prev) + deallocate(elem) + else if (associated(elem, this % tail)) then + this % tail => elem % prev + deallocate(this % tail % next) + else + elem % prev % next => elem % next + elem % next % prev => elem % prev + deallocate(elem) + end if + + ! Decrease count and exit + this % count = this % count - 1 + exit + end if + + ! Advance pointers + elem => elem % next + end do + + end subroutine list_remove_log + !=============================================================================== ! LIST_SIZE returns the number of elements in the list !=============================================================================== @@ -799,4 +1037,13 @@ contains end function list_size_char + function list_size_log(this) result(size) + + class(ListLog) :: this + integer :: size + + size = this % count + + end function list_size_log + end module list_header diff --git a/src/material_header.F90 b/src/material_header.F90 index fca735fd2a..8dd2714006 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -25,6 +25,9 @@ module material_header ! Does this material contain fissionable nuclides? logical :: fissionable = .false. + ! enforce isotropic scattering in lab + logical, allocatable :: p0(:) + end type Material end module material_header diff --git a/src/physics.F90 b/src/physics.F90 index 69bbcff345..624421b00e 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -73,11 +73,12 @@ contains type(Particle), intent(inout) :: p integer :: i_nuclide ! index in nuclides array + integer :: i_nuc_mat ! index in material's nuclides array integer :: i_reaction ! index in nuc % reactions array type(Nuclide), pointer, save :: nuc => null() !$omp threadprivate(nuc) - i_nuclide = sample_nuclide(p, 'total ') + call sample_nuclide(p, 'total ', i_nuclide, i_nuc_mat) ! Get pointer to table nuc => nuclides(i_nuclide) @@ -107,7 +108,7 @@ contains ! Sample a scattering reaction and determine the secondary energy of the ! exiting neutron - call scatter(p, i_nuclide) + call scatter(p, i_nuclide, i_nuc_mat) ! Play russian roulette if survival biasing is turned on @@ -122,13 +123,13 @@ contains ! SAMPLE_NUCLIDE !=============================================================================== - function sample_nuclide(p, base) result(i_nuclide) + subroutine sample_nuclide(p, base, i_nuclide, i_nuc_mat) type(Particle), intent(in) :: p character(7), intent(in) :: base ! which reaction to sample based on - integer :: i_nuclide + integer, intent(out) :: i_nuclide + integer, intent(out) :: i_nuc_mat - integer :: i real(8) :: prob real(8) :: cutoff real(8) :: atom_density ! atom density of nuclide in atom/b-cm @@ -149,20 +150,20 @@ contains cutoff = prn() * material_xs % fission end select - i = 0 + i_nuc_mat = 0 prob = ZERO do while (prob < cutoff) - i = i + 1 + i_nuc_mat = i_nuc_mat + 1 ! Check to make sure that a nuclide was sampled - if (i > mat % n_nuclides) then + if (i_nuc_mat > mat % n_nuclides) then call write_particle_restart(p) call fatal_error("Did not sample any nuclide during collision.") end if ! Find atom density - i_nuclide = mat % nuclide(i) - atom_density = mat % atom_density(i) + i_nuclide = mat % nuclide(i_nuc_mat) + atom_density = mat % atom_density(i_nuc_mat) ! Determine microscopic cross section select case (base) @@ -179,7 +180,7 @@ contains prob = prob + sigma end do - end function sample_nuclide + end subroutine sample_nuclide !=============================================================================== ! SAMPLE_FISSION @@ -303,10 +304,11 @@ contains ! SCATTER !=============================================================================== - subroutine scatter(p, i_nuclide) + subroutine scatter(p, i_nuclide, i_nuc_mat) type(Particle), intent(inout) :: p integer, intent(in) :: i_nuclide + integer, intent(in) :: i_nuc_mat integer :: i integer :: i_grid @@ -334,6 +336,10 @@ contains ! ELASTIC SCATTERING if (micro_xs(i_nuclide) % index_sab /= NONE) then + if (materials(p % material) % p0(i_nuc_mat)) & + & call fatal_error("thermal scattering law data and isotropic lab& + & scattering specified for the same nuclide") + ! S(a,b) scattering call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, & p % E, p % coord0 % uvw, p % mu) @@ -344,7 +350,8 @@ contains ! Perform collision physics for elastic scattering call elastic_scatter(i_nuclide, rxn, & - p % E, p % coord0 % uvw, p % mu, p % wgt) + p % E, p % coord0 % uvw, p % mu, p % wgt, & + & materials(p % material) % p0(i_nuc_mat)) end if p % event_MT = ELASTIC @@ -402,17 +409,19 @@ contains ! target. !=============================================================================== - subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt) + subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt, iso_lab) integer, intent(in) :: i_nuclide type(Reaction), pointer :: rxn real(8), intent(inout) :: E real(8), intent(inout) :: uvw(3) - real(8), intent(out) :: mu_lab real(8), intent(inout) :: wgt + logical, intent(in) :: iso_lab real(8) :: awr ! atomic weight ratio of target real(8) :: mu_cm ! cosine of polar angle in center-of-mass + real(8), intent(out) :: mu_lab ! cosine of polar angle in lab system + real(8) :: phi ! azimuthal angle real(8) :: vel ! magnitude of velocity real(8) :: v_n(3) ! velocity of neutron real(8) :: v_cm(3) ! velocity of center-of-mass @@ -466,10 +475,17 @@ contains ! compute cosine of scattering angle in LAB frame by taking dot product of ! neutron's pre- and post-collision angle - mu_lab = dot_product(uvw, v_n) / vel + if (iso_lab) then + mu_lab = TWO * prn() - ONE + phi = TWO * PI * prn() + uvw = [mu_lab, cos(phi)*sqrt(ONE - mu_lab*mu_lab), & + & sin(phi)*sqrt(ONE - mu_lab*mu_lab)] + else + ! Set energy and direction of particle in LAB frame + uvw = v_n / vel + end if - ! Set energy and direction of particle in LAB frame - uvw = v_n / vel + mu_lab = dot_product(uvw, v_n) / vel end subroutine elastic_scatter From 77a995bbeae985d18a19b1eded7c2d5312f7dd9d Mon Sep 17 00:00:00 2001 From: walshjon Date: Sun, 22 Mar 2015 09:28:36 -0700 Subject: [PATCH 2/2] option to force all collisions to be lab isotropic --- src/physics.F90 | 40 ++++++++++++++++++++++++++++------------ 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 624421b00e..5691c9c0f1 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -394,7 +394,7 @@ contains ! Perform collision physics for inelastic scattering call inelastic_scatter(nuc, rxn, p % E, p % coord0 % uvw, & - p % mu, p % wgt) + p % mu, p % wgt, materials(p % material) % p0(i_nuc_mat)) p % event_MT = rxn % MT end if @@ -416,16 +416,17 @@ contains real(8), intent(inout) :: E real(8), intent(inout) :: uvw(3) real(8), intent(inout) :: wgt + real(8), intent(out) :: mu_lab ! cosine of polar angle in lab system logical, intent(in) :: iso_lab real(8) :: awr ! atomic weight ratio of target real(8) :: mu_cm ! cosine of polar angle in center-of-mass - real(8), intent(out) :: mu_lab ! cosine of polar angle in lab system real(8) :: phi ! azimuthal angle real(8) :: vel ! magnitude of velocity real(8) :: v_n(3) ! velocity of neutron real(8) :: v_cm(3) ! velocity of center-of-mass real(8) :: v_t(3) ! velocity of target nucleus + real(8) :: uvw_in(3) ! incoming direction real(8) :: uvw_cm(3) ! directional cosines in center-of-mass type(Nuclide), pointer, save :: nuc => null() !$omp threadprivate(nuc) @@ -439,6 +440,9 @@ contains ! Neutron velocity in LAB v_n = vel * uvw + ! incoming direction + uvw_in(:) = uvw(:) + ! Sample velocity of target nucleus if (.not. micro_xs(i_nuclide) % use_ptable) then call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, & @@ -474,18 +478,17 @@ contains vel = sqrt(E) ! compute cosine of scattering angle in LAB frame by taking dot product of - ! neutron's pre- and post-collision angle + ! neutron's pre- and post-collision unit vectors if (iso_lab) then - mu_lab = TWO * prn() - ONE + uvw(1) = TWO * prn() - ONE phi = TWO * PI * prn() - uvw = [mu_lab, cos(phi)*sqrt(ONE - mu_lab*mu_lab), & - & sin(phi)*sqrt(ONE - mu_lab*mu_lab)] + uvw(2) = cos(phi) * sqrt(ONE - uvw(1)*uvw(1)) + uvw(3) = sin(phi) * sqrt(ONE - uvw(1)*uvw(1)) else - ! Set energy and direction of particle in LAB frame uvw = v_n / vel end if - mu_lab = dot_product(uvw, v_n) / vel + mu_lab = dot_product(uvw_in, uvw) end subroutine elastic_scatter @@ -1293,7 +1296,7 @@ contains ! than fission), i.e. level scattering, (n,np), (n,na), etc. !=============================================================================== - subroutine inelastic_scatter(nuc, rxn, E, uvw, mu, wgt) + subroutine inelastic_scatter(nuc, rxn, E, uvw, mu, wgt, iso_lab) type(Nuclide), pointer :: nuc type(Reaction), pointer :: rxn @@ -1301,6 +1304,7 @@ contains real(8), intent(inout) :: uvw(3) ! directional cosines real(8), intent(out) :: mu ! cosine of scattering angle in lab real(8), intent(inout) :: wgt ! particle weight + logical, intent(in) :: iso_lab integer :: law ! secondary energy distribution law real(8) :: A ! atomic weight ratio of nuclide @@ -1308,9 +1312,12 @@ contains real(8) :: E_cm ! outgoing energy in center-of-mass real(8) :: Q ! Q-value of reaction real(8) :: yield ! neutron yield + real(8) :: uvw_in(3) ! incoming direction + real(8) :: phi ! azimuthal angle - ! copy energy of neutron + ! copy energy, direction of neutron E_in = E + uvw_in(:) = uvw(:) ! determine A and Q A = nuc % awr @@ -1344,8 +1351,17 @@ contains mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) end if - ! change direction of particle - uvw = rotate_angle(uvw, mu) + ! compute cosine of scattering angle in LAB frame by taking dot product of + ! neutron's pre- and post-collision unit vectors + if (iso_lab) then + uvw(1) = TWO * prn() - ONE + phi = TWO * PI * prn() + uvw(2) = cos(phi) * sqrt(ONE - uvw(1)*uvw(1)) + uvw(3) = sin(phi) * sqrt(ONE - uvw(1)*uvw(1)) + mu = dot_product(uvw_in, uvw) + else + uvw = rotate_angle(uvw_in, mu) + end if ! change weight of particle based on yield if (rxn % multiplicity_with_E) then