From 500843bce431d9efaf9e4ef275501a6535024fc8 Mon Sep 17 00:00:00 2001 From: walshjon Date: Thu, 19 Mar 2015 03:19:04 -0700 Subject: [PATCH 01/24] 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 4c1f3c244..9b6259a57 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 a9311ec64..e4cc56ed0 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 fca735fd2..8dd271400 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 69bbcff34..624421b00 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 02/24] 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 624421b00..5691c9c0f 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 From 99dd8386449424a6c084cffe3aa56ebe07e4912e Mon Sep 17 00:00:00 2001 From: walshjon Date: Sun, 26 Apr 2015 23:35:36 -0400 Subject: [PATCH 03/24] change to input file specification --- src/input_xml.F90 | 11 ++++++----- src/physics.F90 | 6 +++--- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9b6259a57..0df77fe37 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1733,14 +1733,15 @@ contains 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 + if (check_for_node(node_nuc, "scattering")) then + call get_node_value(node_nuc, "scattering", temp_str) + if (trim(adjustl(to_lower(temp_str))) == "lab") then call list_iso_lab % append(.true.) - else if (trim(adjustl(to_lower(temp_str))) == "false") then + else if (trim(adjustl(to_lower(temp_str))) == "ace") then call list_iso_lab % append(.false.) else - call fatal_error("Isotropic lab scattering must be true or false") + call fatal_error("Scattering must be isotropic in lab or follow& + & the ACE file data") end if else call list_iso_lab % append(.false.) diff --git a/src/physics.F90 b/src/physics.F90 index 5691c9c0f..6cda2f9e4 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -337,8 +337,8 @@ contains 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") + 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, & @@ -351,7 +351,7 @@ contains ! Perform collision physics for elastic scattering call elastic_scatter(i_nuclide, rxn, & p % E, p % coord0 % uvw, p % mu, p % wgt, & - & materials(p % material) % p0(i_nuc_mat)) + materials(p % material) % p0(i_nuc_mat)) end if p % event_MT = ELASTIC From 4e808a4257bbb1d399d511307597168a4bb1b2f3 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Mon, 27 Apr 2015 00:39:02 -0400 Subject: [PATCH 04/24] Added paragraph on new scattering subelement for nuclides --- docs/source/usersguide/input.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index d98804c86..c1231e9a8 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1015,6 +1015,15 @@ Each ``material`` element can have the following attributes or sub-elements: .. note:: If one nuclide is specified in atom percent, all others must also be given in atom percent. The same applies for weight percentages. + An optional attribute/sub-element for each nuclide is ``scattering``. This + attribute may be set to "ace" to use the scattering laws specified in the + ACE files (default). Alternatively, when set to "lab", the ACE scattering + laws are used to sample the outgoing energy but an isotropic-in-lab + distribution is used to sample the outgoing angle at each scattering + interaction. The ``scattering`` attribute may be most useful when using + OpenMC to compute multi-group cross-sections for deterministic transport + codes and to quantify the effects of anisotropic scattering. + *Default*: None :element: From 25d9f8b2a3adcaae699813be8921b6b33445ff0f Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Mon, 27 Apr 2015 00:52:49 -0400 Subject: [PATCH 05/24] Added routines to Python API for nuclide-based scattering laws --- src/utils/openmc/material.py | 14 ++++++++++++++ src/utils/openmc/nuclide.py | 24 ++++++++++++++++++++++-- 2 files changed, 36 insertions(+), 2 deletions(-) diff --git a/src/utils/openmc/material.py b/src/utils/openmc/material.py index 87255083b..c7292416e 100644 --- a/src/utils/openmc/material.py +++ b/src/utils/openmc/material.py @@ -271,6 +271,12 @@ class Material(object): self._sab.append((name, xs)) + def make_isotropic_in_lab(self): + + for nuclide_name in self._nuclides: + self._nuclides[nuclide_name][0].make_isotropic_in_lab() + + def get_all_nuclides(self): nuclides = {} @@ -331,6 +337,9 @@ class Material(object): if not nuclide[0]._xs is None: xml_element.set("xs", nuclide[0]._xs) + if not nuclide[0]._scattering is None: + xml_element.set("scattering", nuclide[0]._scattering) + return xml_element @@ -496,6 +505,11 @@ class MaterialsFile(object): self._materials.remove(material) + def make_isotropic_in_lab(self): + + for material in self._materials: + materials.make_isotropic_in_lab() + def create_material_subelements(self): diff --git a/src/utils/openmc/nuclide.py b/src/utils/openmc/nuclide.py index c2287b23d..d9d49d039 100644 --- a/src/utils/openmc/nuclide.py +++ b/src/utils/openmc/nuclide.py @@ -9,6 +9,7 @@ class Nuclide(object): self._name = '' self._xs = None self._zaid = None + self._scattering = None # Set the Material class attributes self.name = name @@ -54,6 +55,11 @@ class Nuclide(object): return self._zaid + @property + def scattering(self): + return self._scattering + + @name.setter def name(self, name): @@ -87,10 +93,24 @@ class Nuclide(object): self._zaid = zaid + @scattering.setter + def scattering(self, scattering): + + if not scattering in ['ace', 'lab']: + msg = 'Unable to set scattering for Nuclide to {0} ' \ + 'which is not "ace" or "lab"'.format(scattering) + raise ValueError(msg) + + self._scattering = scattering + + def __repr__(self): string = 'Nuclide - {0}\n'.format(self._name) string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs) if self._zaid is not None: - string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self._zaid) - return string \ No newline at end of file + string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self._zaid) + if self._scattering is not None: + string += '{0: <16}{1}{2}\n'.format('\tscattering', '=\t', + self._scattering) + return string From e44e253993a89f44b826f28869b8b1d9b78b2da5 Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 27 Apr 2015 15:28:43 -0700 Subject: [PATCH 06/24] re-enable S(a,b) w/ forced iso lab scattering --- src/physics.F90 | 39 +++++++++++++++++++++++++++------------ 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 394277ff3..35cd00de5 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -332,13 +332,11 @@ 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) + p % E, p % coord0 % uvw, p % mu, & + materials(p % material) % p0(i_nuc_mat)) else ! get pointer to elastic scattering reaction @@ -492,13 +490,15 @@ contains ! according to a specified S(a,b) table. !=============================================================================== - subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu) + subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu_lab, iso_lab) integer, intent(in) :: i_nuclide ! index in micro_xs integer, intent(in) :: i_sab ! index in sab_tables real(8), intent(inout) :: E ! incoming/outgoing energy real(8), intent(inout) :: uvw(3) ! directional cosines - real(8), intent(out) :: mu ! scattering cosine + real(8) :: uvw_in(3) ! incoming direction + real(8), intent(out) :: mu_lab ! cosine of polar angle in lab system + logical, intent(in) :: iso_lab integer :: i ! incoming energy bin integer :: j ! outgoing energy bin @@ -522,10 +522,14 @@ contains real(8) :: c_j, c_j1 ! cumulative probability real(8) :: frac ! interpolation factor on outgoing energy real(8) :: r1 ! RNG for outgoing energy + real(8) :: phi ! azimuthal angle ! Get pointer to S(a,b) table sab => sab_tables(i_sab) + ! incoming direction + uvw_in(:) = uvw(:) + ! Determine whether inelastic or elastic scattering will occur if (prn() < micro_xs(i_nuclide) % elastic_sab / & micro_xs(i_nuclide) % elastic) then @@ -555,7 +559,7 @@ contains mu_i1jk = sab % elastic_mu(k,i+1) ! Cosine of angle between incoming and outgoing neutron - mu = (1 - f)*mu_ijk + f*mu_i1jk + mu_lab = (1 - f)*mu_ijk + f*mu_i1jk elseif (sab % elastic_mode == SAB_ELASTIC_EXACT) then ! This treatment is used for data derived in the coherent @@ -571,7 +575,7 @@ contains end if ! Characteristic scattering cosine for this Bragg edge - mu = ONE - TWO*sab % elastic_e_in(k) / E + mu_lab = ONE - TWO*sab % elastic_e_in(k) / E end if @@ -646,7 +650,7 @@ contains 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 + mu_lab = (1 - f)*mu_ijk + f*mu_i1jk else if (sab % secondary_mode == SAB_SECONDARY_CONT) then ! Continuous secondary energy - this is to be similar to @@ -723,7 +727,7 @@ contains ! Will use mu from the randomly chosen incoming and closest outgoing ! energy bins - mu = sab % inelastic_data(l) % mu(k, j) + mu_lab = sab % inelastic_data(l) % mu(k, j) else call fatal_error("Invalid secondary energy mode on S(a,b) table " & @@ -731,8 +735,19 @@ contains end if ! (inelastic secondary energy treatment) end if ! (elastic or inelastic) - ! 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_lab = dot_product(uvw_in, uvw) + else + ! change direction of particle + uvw = rotate_angle(uvw, mu_lab) + end if end subroutine sab_scatter From 77b4b727bdcb41625c176c6c94c0b13b6c62a5b1 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Wed, 29 Apr 2015 23:51:30 -0400 Subject: [PATCH 07/24] Removed ListLog and made iso_lab list ListInt --- src/input_xml.F90 | 18 ++-- src/list_header.F90 | 246 -------------------------------------------- 2 files changed, 11 insertions(+), 253 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a9c761842..01379a231 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, ListLog, ListReal + use list_header, only: ListChar, ListInt, ListReal use mesh_header, only: StructuredMesh use output, only: write_message use plot_header @@ -1581,7 +1581,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(ListInt) :: 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() @@ -1740,15 +1740,15 @@ contains if (check_for_node(node_nuc, "scattering")) then call get_node_value(node_nuc, "scattering", temp_str) if (trim(adjustl(to_lower(temp_str))) == "lab") then - call list_iso_lab % append(.true.) + call list_iso_lab % append(1) else if (trim(adjustl(to_lower(temp_str))) == "ace") then - call list_iso_lab % append(.false.) + call list_iso_lab % append(0) else call fatal_error("Scattering must be isotropic in lab or follow& & the ACE file data") end if else - call list_iso_lab % append(.false.) + call list_iso_lab % append(0) end if ! store full name @@ -1880,8 +1880,12 @@ contains 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) + ! Cast integer isotropic lab scattering flag to boolean + if (list_iso_lab % get_item(j) == 1) then + mat % p0(j) = .true. + else + mat % p0(j) = .false. + end if end do ALL_NUCLIDES diff --git a/src/list_header.F90 b/src/list_header.F90 index e4cc56ed0..9d6c13b2c 100644 --- a/src/list_header.F90 +++ b/src/list_header.F90 @@ -34,12 +34,6 @@ 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 @@ -114,28 +108,6 @@ 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 !=============================================================================== @@ -218,31 +190,6 @@ 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 !=============================================================================== @@ -325,32 +272,6 @@ 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). @@ -383,15 +304,6 @@ 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. @@ -496,39 +408,6 @@ 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. @@ -597,27 +476,6 @@ 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. @@ -789,62 +647,6 @@ 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. @@ -967,45 +769,6 @@ 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 !=============================================================================== @@ -1037,13 +800,4 @@ 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 From 08ad9d7425082a314036fda6cf82a378c09b94b3 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 26 Jul 2015 17:59:11 -0700 Subject: [PATCH 08/24] Removed remaining errant merge conflict metadata from nuclide.py --- openmc/nuclide.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 7de1e900c..25abe815b 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -92,8 +92,6 @@ class Nuclide(object): check_type('zaid', zaid, Integral) self._zaid = zaid -<<<<<<< HEAD - @scattering.setter def scattering(self, scattering): @@ -104,18 +102,12 @@ class Nuclide(object): self._scattering = scattering - -======= ->>>>>>> develop def __repr__(self): string = 'Nuclide - {0}\n'.format(self._name) string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs) if self._zaid is not None: string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self._zaid) -<<<<<<< HEAD if self._scattering is not None: string += '{0: <16}{1}{2}\n'.format('\tscattering', '=\t', self._scattering) -======= ->>>>>>> develop return string From 857a2c5734dfdc709885c2892bbc3038c5ab9eec Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Tue, 3 Nov 2015 10:19:42 -0500 Subject: [PATCH 09/24] All object names in Python API now initialized to an empty string --- openmc/material.py | 2 +- openmc/mesh.py | 2 +- openmc/surface.py | 2 +- openmc/tallies.py | 2 +- openmc/universe.py | 6 +++--- tests/test_track_output/results_test.dat | 9 +++++++++ 6 files changed, 16 insertions(+), 7 deletions(-) create mode 100644 tests/test_track_output/results_test.dat diff --git a/openmc/material.py b/openmc/material.py index 4100c26ec..9fc99c26a 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -193,7 +193,7 @@ class Material(object): name, basestring) self._name = name else: - self._name = None + self._name = '' def set_density(self, units, density=NO_DENSITY): """Set the density of the material diff --git a/openmc/mesh.py b/openmc/mesh.py index 822370c01..3b66076b7 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -155,7 +155,7 @@ class Mesh(object): name, basestring) self._name = name else: - self._name = None + self._name = '' @type.setter def type(self, meshtype): diff --git a/openmc/surface.py b/openmc/surface.py index ed446c374..279246b02 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -128,7 +128,7 @@ class Surface(object): check_type('surface name', name, basestring) self._name = name else: - self._name = None + self._name = '' @boundary_type.setter def boundary_type(self, boundary_type): diff --git a/openmc/tallies.py b/openmc/tallies.py index 575518734..0661ab67b 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -420,7 +420,7 @@ class Tally(object): cv.check_type('tally name', name, basestring) self._name = name else: - self._name = None + self._name = '' def add_filter(self, filter): """Add a filter to the tally diff --git a/openmc/universe.py b/openmc/universe.py index b38ff2c4b..b74dde656 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -153,7 +153,7 @@ class Cell(object): cv.check_type('cell name', name, basestring) self._name = name else: - self._name = None + self._name = '' @fill.setter def fill(self, fill): @@ -472,7 +472,7 @@ class Universe(object): cv.check_type('universe name', name, basestring) self._name = name else: - self._name = None + self._name = '' def add_cell(self, cell): """Add a cell to the universe. @@ -732,7 +732,7 @@ class Lattice(object): cv.check_type('lattice name', name, basestring) self._name = name else: - self._name = None + self._name = '' @outer.setter def outer(self, outer): diff --git a/tests/test_track_output/results_test.dat b/tests/test_track_output/results_test.dat new file mode 100644 index 000000000..6ded87a0e --- /dev/null +++ b/tests/test_track_output/results_test.dat @@ -0,0 +1,9 @@ + + + + + + + + + From 664a2a0e83b53fd53e77d576f1ef635c17b78cd2 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Tue, 3 Nov 2015 10:24:17 -0500 Subject: [PATCH 10/24] Now cast number of subdomains used in MGXS.get_xs(...) to an int --- openmc/mgxs/mgxs.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index d20ad2acf..599078519 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -673,7 +673,7 @@ class MGXS(object): num_groups = len(groups) # Reshape tally data array with separate axes for domain and energy - num_subdomains = xs.shape[0] / num_groups + num_subdomains = int(xs.shape[0] / num_groups) new_shape = (num_subdomains, num_groups) + xs.shape[1:] xs = np.reshape(xs, new_shape) @@ -1844,7 +1844,7 @@ class ScatterMatrixXS(MGXS): num_out_groups = len(out_groups) # Reshape tally data array with separate axes for domain and energy - num_subdomains = xs.shape[0] / (num_in_groups * num_out_groups) + num_subdomains = int(xs.shape[0] / (num_in_groups * num_out_groups)) new_shape = (num_subdomains, num_in_groups, num_out_groups) new_shape += xs.shape[1:] xs = np.reshape(xs, new_shape) @@ -2194,7 +2194,7 @@ class Chi(MGXS): num_groups = self.num_groups else: num_groups = len(groups) - num_subdomains = xs.shape[0] / num_groups + num_subdomains = int(xs.shape[0] / num_groups) new_shape = (num_subdomains, num_groups) + xs.shape[1:] xs = np.reshape(xs, new_shape) From c34f668a3f0a614ba981236fe66d0228cf3dd6bc Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Tue, 3 Nov 2015 10:30:30 -0500 Subject: [PATCH 11/24] Removed test_track_output/results_test.dat --- tests/test_track_output/results_test.dat | 9 --------- 1 file changed, 9 deletions(-) delete mode 100644 tests/test_track_output/results_test.dat diff --git a/tests/test_track_output/results_test.dat b/tests/test_track_output/results_test.dat deleted file mode 100644 index 6ded87a0e..000000000 --- a/tests/test_track_output/results_test.dat +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - From a11ccce8a3d8224413ce31fc2e2e0773bb00c50b Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 16:06:34 -0500 Subject: [PATCH 12/24] Changed isotropic in lab scattering tag to iso-in-lab --- docs/source/usersguide/input.rst | 12 ++++++------ openmc/nuclide.py | 4 ++-- src/input_xml.F90 | 2 +- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 7815d8810..a20ceda18 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1138,12 +1138,12 @@ Each ``material`` element can have the following attributes or sub-elements: An optional attribute/sub-element for each nuclide is ``scattering``. This attribute may be set to "ace" to use the scattering laws specified in the - ACE files (default). Alternatively, when set to "lab", the ACE scattering - laws are used to sample the outgoing energy but an isotropic-in-lab - distribution is used to sample the outgoing angle at each scattering - interaction. The ``scattering`` attribute may be most useful when using - OpenMC to compute multi-group cross-sections for deterministic transport - codes and to quantify the effects of anisotropic scattering. + ACE files (default). Alternatively, when set to "iso-in-lab", the ACE + scattering laws are used to sample the outgoing energy but an + isotropic-in-lab distribution is used to sample the outgoing angle at each + scattering interaction. The ``scattering`` attribute may be most useful + when using OpenMC to compute multi-group cross-sections for deterministic + transport codes and to quantify the effects of anisotropic scattering. *Default*: None diff --git a/openmc/nuclide.py b/openmc/nuclide.py index c2a4d8455..84ad112bc 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -102,9 +102,9 @@ class Nuclide(object): @scattering.setter def scattering(self, scattering): - if not scattering in ['ace', 'lab']: + if not scattering in ['ace', 'iso-in-lab']: msg = 'Unable to set scattering for Nuclide to {0} ' \ - 'which is not "ace" or "lab"'.format(scattering) + 'which is not "ace" or "iso-in-lab"'.format(scattering) raise ValueError(msg) self._scattering = scattering diff --git a/src/input_xml.F90 b/src/input_xml.F90 index b498b4ca7..70b9302df 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1938,7 +1938,7 @@ contains ! Check enforced isotropic lab scattering if (check_for_node(node_nuc, "scattering")) then call get_node_value(node_nuc, "scattering", temp_str) - if (trim(adjustl(to_lower(temp_str))) == "lab") then + if (trim(adjustl(to_lower(temp_str))) == "iso-in-lab") then call list_iso_lab % append(1) else if (trim(adjustl(to_lower(temp_str))) == "ace") then call list_iso_lab % append(0) From 0360983ecd335ef924f8f0df1aaefd4304a411d0 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 16:11:15 -0500 Subject: [PATCH 13/24] Removed unnecessary whitespace --- src/list_header.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/list_header.F90 b/src/list_header.F90 index 27721a7a6..8020b96b1 100644 --- a/src/list_header.F90 +++ b/src/list_header.F90 @@ -34,7 +34,6 @@ module list_header type(ListElemChar), pointer :: prev => null() end type ListElemChar - !=============================================================================== ! LIST* types contain the linked list with convenience methods. We originally ! considered using unlimited polymorphism to provide a single type, but compiler From 151fa4d47942a642ed8962d3904f326850233fe9 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 16:57:29 -0500 Subject: [PATCH 14/24] Refactored iso-in-lab scattering in physics.F90 to consolidate within scattering routine --- openmc/material.py | 4 +- src/physics.F90 | 135 ++++++++++++++++++--------------------------- 2 files changed, 57 insertions(+), 82 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 8ef4c8c1b..e100b7539 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -373,7 +373,7 @@ class Material(object): def make_isotropic_in_lab(self): for nuclide_name in self._nuclides: - self._nuclides[nuclide_name][0].make_isotropic_in_lab() + self._nuclides[nuclide_name][0].scattering = 'iso-in-lab' def get_all_nuclides(self): """Returns all nuclides in the material @@ -599,7 +599,7 @@ class MaterialsFile(object): def make_isotropic_in_lab(self): for material in self._materials: - materials.make_isotropic_in_lab() + material.make_isotropic_in_lab() def _create_material_subelements(self): subelement = ET.SubElement(self._materials_file, "default_xs") diff --git a/src/physics.F90 b/src/physics.F90 index 54528b5ba..03fc8a2de 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -314,6 +314,13 @@ contains real(8) :: cutoff type(Nuclide), pointer :: nuc type(Reaction), pointer :: rxn + real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering + real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering + real(8) :: mu_lab ! polar angle cosine for iso-in-lab scattering + real(8) :: phi ! azimuthal angle for iso-in-lab scattering + + ! copy incoming direction + uvw_old(:) = p % coord(1) % uvw ! Get pointer to nuclide and grid index/interpolation factor nuc => nuclides(i_nuclide) @@ -332,11 +339,9 @@ contains ! ELASTIC SCATTERING if (micro_xs(i_nuclide) % index_sab /= NONE) then - ! S(a,b) scattering call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, & - p % E, p % coord(1) % uvw, p % mu, & - materials(p % material) % p0(i_nuc_mat)) + p % E, p % coord(1) % uvw, p % mu) else ! get pointer to elastic scattering reaction @@ -344,8 +349,7 @@ contains ! Perform collision physics for elastic scattering call elastic_scatter(i_nuclide, rxn, & - p % E, p % coord(1) % uvw, p % mu, p % wgt, & - materials(p % material) % p0(i_nuc_mat)) + p % E, p % coord(1) % uvw, p % mu, p % wgt) end if p % event_MT = ELASTIC @@ -387,7 +391,7 @@ contains end do ! Perform collision physics for inelastic scattering - call inelastic_scatter(nuc, rxn, p, materials(p % material) % p0(i_nuc_mat)) + call inelastic_scatter(nuc, rxn, p) p % event_MT = rxn % MT end if @@ -395,6 +399,20 @@ contains ! Set event component p % event = EVENT_SCATTER + ! sample new outgoing angle for isotropic in lab scattering + if (materials(p % material) % p0(i_nuc_mat)) then + + ! sample isotropic-in-lab outgoing direction + uvw_new(1) = TWO * prn() - ONE + phi = TWO * PI * prn() + uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) + uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) + mu_lab = dot_product(uvw_old, uvw_new) + + ! change direction of particle + p % coord(1) % uvw = rotate_angle(uvw_new, mu_lab) + end if + end subroutine scatter !=============================================================================== @@ -402,24 +420,21 @@ contains ! target. !=============================================================================== - subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt, iso_lab) + subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt) 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 - 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) :: 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 :: nuc @@ -432,9 +447,6 @@ 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, & @@ -470,17 +482,11 @@ contains vel = sqrt(E) ! 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)) - else - uvw = v_n / vel - end if + ! neutron's pre- and post-collision angle + mu_lab = dot_product(uvw, v_n) / vel - mu_lab = dot_product(uvw_in, uvw) + ! Set energy and direction of particle in LAB frame + uvw = v_n / vel ! Because of floating-point roundoff, it may be possible for mu_lab to be ! outside of the range [-1,1). In these cases, we just set mu_lab to exactly @@ -495,15 +501,13 @@ contains ! according to a specified S(a,b) table. !=============================================================================== - subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu_lab, iso_lab) + subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu) integer, intent(in) :: i_nuclide ! index in micro_xs integer, intent(in) :: i_sab ! index in sab_tables real(8), intent(inout) :: E ! incoming/outgoing energy real(8), intent(inout) :: uvw(3) ! directional cosines - real(8) :: uvw_in(3) ! incoming direction - real(8), intent(out) :: mu_lab ! cosine of polar angle in lab system - logical, intent(in) :: iso_lab + real(8), intent(out) :: mu ! scattering cosine integer :: i ! incoming energy bin integer :: j ! outgoing energy bin @@ -527,14 +531,10 @@ contains real(8) :: c_j, c_j1 ! cumulative probability real(8) :: frac ! interpolation factor on outgoing energy real(8) :: r1 ! RNG for outgoing energy - real(8) :: phi ! azimuthal angle ! Get pointer to S(a,b) table sab => sab_tables(i_sab) - ! incoming direction - uvw_in(:) = uvw(:) - ! Determine whether inelastic or elastic scattering will occur if (prn() < micro_xs(i_nuclide) % elastic_sab / & micro_xs(i_nuclide) % elastic) then @@ -564,7 +564,7 @@ contains mu_i1jk = sab % elastic_mu(k,i+1) ! Cosine of angle between incoming and outgoing neutron - mu_lab = (1 - f)*mu_ijk + f*mu_i1jk + mu = (1 - f)*mu_ijk + f*mu_i1jk elseif (sab % elastic_mode == SAB_ELASTIC_EXACT) then ! This treatment is used for data derived in the coherent @@ -580,7 +580,7 @@ contains end if ! Characteristic scattering cosine for this Bragg edge - mu_lab = ONE - TWO*sab % elastic_e_in(k) / E + mu = ONE - TWO*sab % elastic_e_in(k) / E end if @@ -655,7 +655,7 @@ contains mu_i1jk = sab % inelastic_mu(k,j,i+1) ! Cosine of angle between incoming and outgoing neutron - mu_lab = (1 - f)*mu_ijk + f*mu_i1jk + mu = (1 - f)*mu_ijk + f*mu_i1jk else if (sab % secondary_mode == SAB_SECONDARY_CONT) then ! Continuous secondary energy - this is to be similar to @@ -732,7 +732,7 @@ contains ! Will use mu from the randomly chosen incoming and closest outgoing ! energy bins - mu_lab = sab % inelastic_data(l) % mu(k, j) + mu = sab % inelastic_data(l) % mu(k, j) else call fatal_error("Invalid secondary energy mode on S(a,b) table " & @@ -744,20 +744,10 @@ contains ! outside of the range [-1,1). In these cases, we just set mu to exactly ! -1 or 1 - if (abs(mu_lab) > ONE) mu_lab = sign(ONE,mu_lab) + if (abs(mu) > ONE) mu = sign(ONE,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_lab = dot_product(uvw_in, uvw) - else - ! change direction of particle - uvw = rotate_angle(uvw, mu_lab) - end if + ! change direction of particle + uvw = rotate_angle(uvw, mu) end subroutine sab_scatter @@ -1333,28 +1323,23 @@ contains ! than fission), i.e. level scattering, (n,np), (n,na), etc. !=============================================================================== - subroutine inelastic_scatter(nuc, rxn, p, iso_lab) - - type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn + subroutine inelastic_scatter(nuc, rxn, p) + type(Nuclide), pointer :: nuc + type(Reaction), pointer :: rxn type(Particle), intent(inout) :: p - logical, intent(in) :: iso_lab integer :: i ! loop index - integer :: law ! secondary energy distribution law - real(8) :: mu ! cosine of scattering angle in lab - real(8) :: A ! atomic weight ratio of nuclide + integer :: law ! secondary energy distribution law real(8) :: E ! energy in lab (incoming/outgoing) - real(8) :: E_in ! incoming energy - 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 + real(8) :: mu ! cosine of scattering angle in lab + real(8) :: A ! atomic weight ratio of nuclide + real(8) :: E_in ! incoming energy + real(8) :: E_cm ! outgoing energy in center-of-mass + real(8) :: Q ! Q-value of reaction + real(8) :: yield ! neutron yield - ! copy energy, direction of neutron + ! copy energy of neutron E_in = p % E - uvw_in(:) = p % coord(1) % uvw(:) ! determine A and Q A = nuc % awr @@ -1399,22 +1384,12 @@ contains if (abs(mu) > ONE) mu = sign(ONE,mu) end if - ! 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 - p % coord(1) % uvw(1) = TWO * prn() - ONE - phi = TWO * PI * prn() - p % coord(1) % uvw(2) = cos(phi) * sqrt(ONE - uvw_in(1) * uvw_in(1)) - p % coord(1) % uvw(3) = sin(phi) * sqrt(ONE - uvw_in(1) * uvw_in(1)) - mu = dot_product(uvw_in, p % coord(1) % uvw) - else - ! Set outgoing energy and scattering angle - p % E = E - p % mu = mu + ! Set outgoing energy and scattering angle + p % E = E + p % mu = mu - ! change direction of particle - p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) - end if + ! change direction of particle + p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) ! change weight of particle based on yield if (rxn % multiplicity_with_E) then From 6e827c25e7333893e100c7b5dd276ccfca9a6a10 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 17:02:55 -0500 Subject: [PATCH 15/24] Removed unnecessary trailing whitespace from physics.F90 --- src/physics.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics.F90 b/src/physics.F90 index 03fc8a2de..00db08950 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -316,7 +316,7 @@ contains type(Reaction), pointer :: rxn real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering - real(8) :: mu_lab ! polar angle cosine for iso-in-lab scattering + real(8) :: mu_lab ! polar angle cosine for iso-in-lab scattering real(8) :: phi ! azimuthal angle for iso-in-lab scattering ! copy incoming direction From 81116a20496689a53de9ec60546f17b76a885c65 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 17:34:04 -0500 Subject: [PATCH 16/24] Elements can now use iso-in-lab scattering in adddition to Nuclides --- openmc/element.py | 23 ++++++++++++++++++++++- openmc/material.py | 13 +++++++++---- openmc/nuclide.py | 2 ++ src/input_xml.F90 | 29 ++++++++++++++++++++++++++++- 4 files changed, 61 insertions(+), 6 deletions(-) diff --git a/openmc/element.py b/openmc/element.py index 56821b5d2..f39bcbfef 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -24,6 +24,8 @@ class Element(object): Chemical symbol of the element, e.g. Pu xs : str Cross section identifier, e.g. 71c + scattering : 'ace' or 'iso-in-lab' or None + The type of angular scattering distribution to use """ @@ -31,6 +33,7 @@ class Element(object): # Initialize class attributes self._name = '' self._xs = None + self._scattering = None # Set class attributes self.name = name @@ -60,6 +63,10 @@ class Element(object): def __repr__(self): string = 'Element - {0}\n'.format(self._name) string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs) + if self._scattering is not None: + string += '{0: <16}{1}{2}\n'.format('\tscattering', '=\t', + self._scattering) + return string @property @@ -70,6 +77,10 @@ class Element(object): def name(self): return self._name + @property + def scattering(self): + return self._scattering + @xs.setter def xs(self, xs): check_type('cross section identifier', xs, basestring) @@ -78,4 +89,14 @@ class Element(object): @name.setter def name(self, name): check_type('name', name, basestring) - self._name = name \ No newline at end of file + self._name = name + + @scattering.setter + def scattering(self, scattering): + + if not scattering in ['ace', 'iso-in-lab']: + msg = 'Unable to set scattering for Element to {0} ' \ + 'which is not "ace" or "iso-in-lab"'.format(scattering) + raise ValueError(msg) + + self._scattering = scattering diff --git a/openmc/material.py b/openmc/material.py index e100b7539..98d0d3f87 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -374,6 +374,8 @@ class Material(object): def make_isotropic_in_lab(self): for nuclide_name in self._nuclides: self._nuclides[nuclide_name][0].scattering = 'iso-in-lab' + for element_name in self._elements: + self._element[element_name][0].scattering = 'iso-in-lab' def get_all_nuclides(self): """Returns all nuclides in the material @@ -405,11 +407,11 @@ class Material(object): else: xml_element.set("wo", str(nuclide[1])) - if nuclide[0]._xs is not None: - xml_element.set("xs", nuclide[0]._xs) + if nuclide[0].xs is not None: + xml_element.set("xs", nuclide[0].xs) - if not nuclide[0]._scattering is None: - xml_element.set("scattering", nuclide[0]._scattering) + if not nuclide[0].scattering is None: + xml_element.set("scattering", nuclide[0].scattering) return xml_element @@ -423,6 +425,9 @@ class Material(object): else: xml_element.set("wo", str(element[1])) + if not element[0].scattering is None: + xml_element.set("scattering", element[0].scattering) + return xml_element def _get_nuclides_xml(self, nuclides, distrib=False): diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 84ad112bc..1c926dc3d 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -26,6 +26,8 @@ class Nuclide(object): zaid : int 1000*(atomic number) + mass number. As an example, the zaid of U-235 would be 92235. + scattering : 'ace' or 'iso-in-lab' or None + The type of angular scattering distribution to use """ diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 70b9302df..1e1492efc 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1759,8 +1759,10 @@ contains integer :: i ! loop index for materials integer :: j ! loop index for nuclides + integer :: k ! llop index for elements integer :: n ! number of nuclides integer :: n_sab ! number of sab tables for a material + integer :: n_nuc_ele ! number of nuclides in an element integer :: index_list ! index in xs_listings array integer :: index_nuclide ! index in nuclides integer :: index_sab ! index in sab_tables @@ -2021,6 +2023,9 @@ contains &element: " // trim(name)) end if + ! Get current number of nuclides + n_nuc_ele = list_names % size() + ! Expand element into naturally-occurring isotopes if (check_for_node(node_ele, "ao")) then call get_node_value(node_ele, "ao", temp_dble) @@ -2030,6 +2035,29 @@ contains call fatal_error("The ability to expand a natural element based on & &weight percentage is not yet supported.") end if + + ! Compute number of new nuclides from the natural element expansion + n_nuc_ele = list_names % size() - n_nuc_ele + + ! Check enforced isotropic lab scattering + if (check_for_node(node_ele, "scattering")) then + call get_node_value(node_ele, "scattering", temp_str) + else + temp_str = "ace" + end if + + ! Set ace or iso-in-lab scattering for each nuclide in element + do k = 1, n_nuc_ele + if (trim(adjustl(to_lower(temp_str))) == "iso-in-lab") then + call list_iso_lab % append(1) + else if (trim(adjustl(to_lower(temp_str))) == "ace") then + call list_iso_lab % append(0) + else + call fatal_error("Scattering must be isotropic in lab or follow& + & the ACE file data") + end if + end do + end do NATURAL_ELEMENTS ! ======================================================================== @@ -4213,7 +4241,6 @@ contains call list_density % append(density * 0.999885_8) call list_names % append('1002.' // xs) call list_density % append(density * 0.000115_8) - case ('he') call list_names % append('2003.' // xs) call list_density % append(density * 0.00000134_8) From 96578a515a63f23d220d2839e968e56de7de1cdd Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 17:35:22 -0500 Subject: [PATCH 17/24] Updated docs with iso-in-lab scattering for elements --- docs/source/usersguide/input.rst | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index a20ceda18..ce1d07ca0 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1141,7 +1141,7 @@ Each ``material`` element can have the following attributes or sub-elements: ACE files (default). Alternatively, when set to "iso-in-lab", the ACE scattering laws are used to sample the outgoing energy but an isotropic-in-lab distribution is used to sample the outgoing angle at each - scattering interaction. The ``scattering`` attribute may be most useful + scattering interaction. The ``scattering`` attribute may be most useful when using OpenMC to compute multi-group cross-sections for deterministic transport codes and to quantify the effects of anisotropic scattering. @@ -1171,6 +1171,16 @@ Each ``material`` element can have the following attributes or sub-elements: *Default*: None + An optional attribute/sub-element for each element is ``scattering``. This + attribute may be set to "ace" to use the scattering laws specified in the + ACE files (default). Alternatively, when set to "iso-in-lab", the ACE + scattering laws are used to sample the outgoing energy but an + isotropic-in-lab distribution is used to sample the outgoing angle at each + scattering interaction. The ``scattering`` attribute may be most useful + when using OpenMC to compute multi-group cross-sections for deterministic + transport codes and to quantify the effects of anisotropic scattering. + + *Default*: None :sab: Associates an S(a,b) table with the material. This element has From 3e34c5882c79ada0552a2492937fda149e39c6af Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 17:36:59 -0500 Subject: [PATCH 18/24] Now using property getters in Nuclide and Element __repr__ routines --- openmc/element.py | 4 ++-- openmc/nuclide.py | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/openmc/element.py b/openmc/element.py index f39bcbfef..3a15a8f63 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -63,9 +63,9 @@ class Element(object): def __repr__(self): string = 'Element - {0}\n'.format(self._name) string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs) - if self._scattering is not None: + if self.scattering is not None: string += '{0: <16}{1}{2}\n'.format('\tscattering', '=\t', - self._scattering) + self.scattering) return string diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 1c926dc3d..c5bc4e078 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -113,10 +113,10 @@ class Nuclide(object): def __repr__(self): string = 'Nuclide - {0}\n'.format(self._name) - string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs) - if self._zaid is not None: - string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self._zaid) - if self._scattering is not None: + string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self.xs) + if self.zaid is not None: + string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self.zaid) + if self.scattering is not None: string += '{0: <16}{1}{2}\n'.format('\tscattering', '=\t', - self._scattering) + self.scattering) return string From 3bbefb7169ac07c9e5129600f271e999be68add7 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 17:39:13 -0500 Subject: [PATCH 19/24] Fixed mis-spelled comment in physics.F90 --- src/input_xml.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1e1492efc..1ee2a060e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1759,7 +1759,7 @@ contains integer :: i ! loop index for materials integer :: j ! loop index for nuclides - integer :: k ! llop index for elements + integer :: k ! loop index for elements integer :: n ! number of nuclides integer :: n_sab ! number of sab tables for a material integer :: n_nuc_ele ! number of nuclides in an element From 4b70c4abc16714ecb1fe1cfff70912a1c0b93f16 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 17:56:42 -0500 Subject: [PATCH 20/24] Removed trailing whitespace in input_xml.F90 --- src/input_xml.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1ee2a060e..df10d0f70 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2047,7 +2047,7 @@ contains end if ! Set ace or iso-in-lab scattering for each nuclide in element - do k = 1, n_nuc_ele + do k = 1, n_nuc_ele if (trim(adjustl(to_lower(temp_str))) == "iso-in-lab") then call list_iso_lab % append(1) else if (trim(adjustl(to_lower(temp_str))) == "ace") then From 4db72c90af2dd62d52a199f4e86ef26676abab53 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Sun, 8 Nov 2015 20:53:18 -0500 Subject: [PATCH 21/24] Fixed issues in iso-in-lab scattering in physics.F90 per comments by @walshjon --- src/physics.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 00db08950..9f8fb535b 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -316,7 +316,6 @@ contains type(Reaction), pointer :: rxn real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering - real(8) :: mu_lab ! polar angle cosine for iso-in-lab scattering real(8) :: phi ! azimuthal angle for iso-in-lab scattering ! copy incoming direction @@ -407,10 +406,10 @@ contains phi = TWO * PI * prn() uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) - mu_lab = dot_product(uvw_old, uvw_new) + p % mu = dot_product(uvw_old, uvw_new) ! change direction of particle - p % coord(1) % uvw = rotate_angle(uvw_new, mu_lab) + p % coord(1) % uvw = uvw_new end if end subroutine scatter From d03272867fbb4a908e74290fad04943517a26014 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Mon, 9 Nov 2015 00:11:50 -0500 Subject: [PATCH 22/24] Made MGXS subdomain averaging use scalar flux weighting --- openmc/mgxs/mgxs.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 599078519..1296f195c 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -824,9 +824,9 @@ class MGXS(object): mean = tally.get_reshaped_data(value='mean') std_dev = tally.get_reshaped_data(value='std_dev') - # Get the mean of the mean, std. dev. across requested subdomains - mean = np.mean(mean[subdomains, ...], axis=0) - std_dev = np.mean(std_dev[subdomains, ...]**2, axis=0) + # Get the mean, std. dev. across requested subdomains + mean = np.sum(mean[subdomains, ...], axis=0) + std_dev = np.sum(std_dev[subdomains, ...]**2, axis=0) std_dev = np.sqrt(std_dev) # If domain is distribcell, make subdomain-averaged a 'cell' domain From baa0b4c10cba06e174bec57fc41f73682642a5f2 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Mon, 9 Nov 2015 09:26:48 -0500 Subject: [PATCH 23/24] Now using a "data" descriptor rather than "ace" for the nuclide/element scattering XML tag --- docs/source/usersguide/input.rst | 12 ++++++------ openmc/element.py | 6 +++--- openmc/material.py | 4 ++-- openmc/nuclide.py | 6 +++--- src/input_xml.F90 | 10 +++++----- 5 files changed, 19 insertions(+), 19 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index ce1d07ca0..04a6406dc 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1137,9 +1137,9 @@ Each ``material`` element can have the following attributes or sub-elements: be given in atom percent. The same applies for weight percentages. An optional attribute/sub-element for each nuclide is ``scattering``. This - attribute may be set to "ace" to use the scattering laws specified in the - ACE files (default). Alternatively, when set to "iso-in-lab", the ACE - scattering laws are used to sample the outgoing energy but an + attribute may be set to "data" to use the scattering laws specified by the + cross section library (default). Alternatively, when set to "iso-in-lab", + the scattering laws are used to sample the outgoing energy but an isotropic-in-lab distribution is used to sample the outgoing angle at each scattering interaction. The ``scattering`` attribute may be most useful when using OpenMC to compute multi-group cross-sections for deterministic @@ -1172,9 +1172,9 @@ Each ``material`` element can have the following attributes or sub-elements: *Default*: None An optional attribute/sub-element for each element is ``scattering``. This - attribute may be set to "ace" to use the scattering laws specified in the - ACE files (default). Alternatively, when set to "iso-in-lab", the ACE - scattering laws are used to sample the outgoing energy but an + attribute may be set to "data" to use the scattering laws specified by the + cross section library (default). Alternatively, when set to "iso-in-lab", + the scattering laws are used to sample the outgoing energy but an isotropic-in-lab distribution is used to sample the outgoing angle at each scattering interaction. The ``scattering`` attribute may be most useful when using OpenMC to compute multi-group cross-sections for deterministic diff --git a/openmc/element.py b/openmc/element.py index 3a15a8f63..d395b434f 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -24,7 +24,7 @@ class Element(object): Chemical symbol of the element, e.g. Pu xs : str Cross section identifier, e.g. 71c - scattering : 'ace' or 'iso-in-lab' or None + scattering : 'data' or 'iso-in-lab' or None The type of angular scattering distribution to use """ @@ -94,9 +94,9 @@ class Element(object): @scattering.setter def scattering(self, scattering): - if not scattering in ['ace', 'iso-in-lab']: + if not scattering in ['data', 'iso-in-lab']: msg = 'Unable to set scattering for Element to {0} ' \ - 'which is not "ace" or "iso-in-lab"'.format(scattering) + 'which is not "data" or "iso-in-lab"'.format(scattering) raise ValueError(msg) self._scattering = scattering diff --git a/openmc/material.py b/openmc/material.py index 98d0d3f87..292fc82ca 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -33,8 +33,8 @@ NO_DENSITY = 99999. class Material(object): - """A material composed of a collection of nuclides/elements that can be assigned - to a region of space. + """A material composed of a collection of nuclides/elements that can be + assigned to a region of space. Parameters ---------- diff --git a/openmc/nuclide.py b/openmc/nuclide.py index c5bc4e078..2dd8eb153 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -26,7 +26,7 @@ class Nuclide(object): zaid : int 1000*(atomic number) + mass number. As an example, the zaid of U-235 would be 92235. - scattering : 'ace' or 'iso-in-lab' or None + scattering : 'data' or 'iso-in-lab' or None The type of angular scattering distribution to use """ @@ -104,9 +104,9 @@ class Nuclide(object): @scattering.setter def scattering(self, scattering): - if not scattering in ['ace', 'iso-in-lab']: + if not scattering in ['data', 'iso-in-lab']: msg = 'Unable to set scattering for Nuclide to {0} ' \ - 'which is not "ace" or "iso-in-lab"'.format(scattering) + 'which is not "data" or "iso-in-lab"'.format(scattering) raise ValueError(msg) self._scattering = scattering diff --git a/src/input_xml.F90 b/src/input_xml.F90 index df10d0f70..ae01b0b5c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1940,9 +1940,9 @@ contains ! Check enforced isotropic lab scattering if (check_for_node(node_nuc, "scattering")) then call get_node_value(node_nuc, "scattering", temp_str) - if (trim(adjustl(to_lower(temp_str))) == "iso-in-lab") then + if (adjustl(to_lower(temp_str)) == "iso-in-lab") then call list_iso_lab % append(1) - else if (trim(adjustl(to_lower(temp_str))) == "ace") then + else if (adjustl(to_lower(temp_str)) == "data") then call list_iso_lab % append(0) else call fatal_error("Scattering must be isotropic in lab or follow& @@ -2043,14 +2043,14 @@ contains if (check_for_node(node_ele, "scattering")) then call get_node_value(node_ele, "scattering", temp_str) else - temp_str = "ace" + temp_str = "data" end if ! Set ace or iso-in-lab scattering for each nuclide in element do k = 1, n_nuc_ele - if (trim(adjustl(to_lower(temp_str))) == "iso-in-lab") then + if (adjustl(to_lower(temp_str)) == "iso-in-lab") then call list_iso_lab % append(1) - else if (trim(adjustl(to_lower(temp_str))) == "ace") then + else if (adjustl(to_lower(temp_str)) == "data") then call list_iso_lab % append(0) else call fatal_error("Scattering must be isotropic in lab or follow& From a21cb9d003ebf98374ea20cf7872dd6dbc75caea Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Mon, 9 Nov 2015 09:54:11 -0500 Subject: [PATCH 24/24] Updated relaxng schema for nuclide/element scattering tag --- src/relaxng/materials.rnc | 4 ++++ src/relaxng/materials.rng | 46 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+) diff --git a/src/relaxng/materials.rnc b/src/relaxng/materials.rnc index ae8565497..da68ebf9e 100644 --- a/src/relaxng/materials.rnc +++ b/src/relaxng/materials.rnc @@ -15,6 +15,8 @@ element materials { attribute name { xsd:string { maxLength = "7" } }) & (element xs { xsd:string { maxLength = "3" } } | attribute xs { xsd:string { maxLength = "3" } })? & + (element scattering { ( "data" | "iso-in-lab" ) } | + attribute scattering { ( "data" | "iso-in-lab" ) })? & ( (element ao { xsd:double } | attribute ao { xsd:double }) | (element wo { xsd:double } | attribute wo { xsd:double }) @@ -26,6 +28,8 @@ element materials { attribute name { xsd:string { maxLength = "2" } }) & (element xs { xsd:string { maxLength = "3" } } | attribute xs { xsd:string { maxLength = "3" } })? & + (element scattering { ( "data" | "iso-in-lab" ) } | + attribute scattering { ( "data" | "iso-in-lab" ) })? & ( (element ao { xsd:double } | attribute ao { xsd:double }) | (element wo { xsd:double } | attribute wo { xsd:double }) diff --git a/src/relaxng/materials.rng b/src/relaxng/materials.rng index 4c9496eae..7aba5f738 100644 --- a/src/relaxng/materials.rng +++ b/src/relaxng/materials.rng @@ -12,6 +12,20 @@ + + + + + 52 + + + + + 52 + + + + @@ -67,6 +81,22 @@ + + + + + data + iso-in-lab + + + + + data + iso-in-lab + + + + @@ -117,6 +147,22 @@ + + + + + data + iso-in-lab + + + + + data + iso-in-lab + + + +