diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 3e7a6ec17..3da2382a8 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1175,6 +1175,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 "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 + transport codes and to quantify the effects of anisotropic scattering. + *Default*: None :element: @@ -1201,6 +1210,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 "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 + 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 diff --git a/openmc/element.py b/openmc/element.py index 56821b5d2..d395b434f 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 : 'data' 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 ['data', 'iso-in-lab']: + msg = 'Unable to set scattering for Element to {0} ' \ + '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 4100c26ec..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 ---------- @@ -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 @@ -371,6 +371,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].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 @@ -401,8 +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) return xml_element @@ -416,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): @@ -590,6 +602,10 @@ class MaterialsFile(object): self._materials.remove(material) + def make_isotropic_in_lab(self): + for material in self._materials: + material.make_isotropic_in_lab() + def _create_material_subelements(self): subelement = ET.SubElement(self._materials_file, "default_xs") 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/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index d20ad2acf..1296f195c 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) @@ -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 @@ -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) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index f24e62a48..2dd8eb153 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 : 'data' or 'iso-in-lab' or None + The type of angular scattering distribution to use """ @@ -34,6 +36,7 @@ class Nuclide(object): self._name = '' self._xs = None self._zaid = None + self._scattering = None # Set the Material class attributes self.name = name @@ -79,6 +82,10 @@ class Nuclide(object): def zaid(self): return self._zaid + @property + def scattering(self): + return self._scattering + @name.setter def name(self, name): check_type('name', name, basestring) @@ -92,4 +99,24 @@ class Nuclide(object): @zaid.setter def zaid(self, zaid): check_type('zaid', zaid, Integral) - self._zaid = zaid \ No newline at end of file + self._zaid = zaid + + @scattering.setter + def scattering(self, scattering): + + if not scattering in ['data', 'iso-in-lab']: + msg = 'Unable to set scattering for Nuclide to {0} ' \ + 'which is not "data" or "iso-in-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) + if self.scattering is not None: + string += '{0: <16}{1}{2}\n'.format('\tscattering', '=\t', + self.scattering) + return string 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/src/input_xml.F90 b/src/input_xml.F90 index 47370e39d..b95f9f483 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, Lattice, RectLattice, HexLattice use global - use list_header, only: ListChar, ListReal + use list_header, only: ListChar, ListInt, ListReal use mesh_header, only: RegularMesh use output, only: write_message use plot_header @@ -1802,8 +1802,10 @@ contains integer :: i ! loop index for materials integer :: j ! loop index for nuclides + 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 integer :: index_list ! index in xs_listings array integer :: index_nuclide ! index in nuclides integer :: index_sab ! index in sab_tables @@ -1818,6 +1820,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(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() @@ -1998,6 +2001,21 @@ contains end if end if + ! Check enforced isotropic lab scattering + if (check_for_node(node_nuc, "scattering")) then + call get_node_value(node_nuc, "scattering", temp_str) + if (adjustl(to_lower(temp_str)) == "iso-in-lab") then + call list_iso_lab % append(1) + 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& + & the ACE file data") + end if + else + call list_iso_lab % append(0) + end if + ! store full name call get_node_value(node_nuc, "name", temp_str) if (check_for_node(node_nuc, "xs")) & @@ -2119,6 +2137,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) @@ -2128,6 +2149,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 = "data" + end if + + ! Set ace or iso-in-lab scattering for each nuclide in element + do k = 1, n_nuc_ele + if (adjustl(to_lower(temp_str)) == "iso-in-lab") then + call list_iso_lab % append(1) + 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& + & the ACE file data") + end if + end do + end do NATURAL_ELEMENTS ! ======================================================================== @@ -2139,6 +2183,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 @@ -2177,6 +2222,14 @@ contains ! Copy name and atom/weight percent mat % names(j) = name mat % atom_density(j) = list_density % 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 ! Check to make sure either all atom percents or all weight percents are @@ -2193,6 +2246,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 @@ -4483,7 +4537,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) diff --git a/src/material_header.F90 b/src/material_header.F90 index c2f4b51bc..91c4cdfb8 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -33,6 +33,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 3ea28c228..672d26b38 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -75,10 +75,11 @@ 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_CE), pointer :: 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) @@ -108,7 +109,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 @@ -123,13 +124,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 + character(7), intent(in) :: base ! which reaction to sample based on + 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 @@ -281,10 +282,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_nuclide + integer, intent(in) :: i_nuc_mat integer :: i integer :: i_grid @@ -293,6 +295,12 @@ contains real(8) :: cutoff type(Nuclide_CE), 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) :: 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) @@ -371,6 +379,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)) + p % mu = dot_product(uvw_old, uvw_new) + + ! change direction of particle + p % coord(1) % uvw = uvw_new + end if + end subroutine scatter !=============================================================================== 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 + + + +