crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merged iso-in-lab

Browse source
This commit is contained in:
Adam Nelson 2015-11-10 05:32:52 -05:00
commit f7ffa71e4f
14 changed files with 246 additions and 35 deletions
Download patch file Download diff file Expand all files Collapse all files

19
docs/source/usersguide/input.rst
View file

@ -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

23
openmc/element.py
View file

@ -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
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

26
openmc/material.py
View file

@ -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")

2
openmc/mesh.py
View file

@ -155,7 +155,7 @@ class Mesh(object):
name, basestring)
self._name = name
else:
self._name = None
self._name = ''
@type.setter
def type(self, meshtype):

12
openmc/mgxs/mgxs.py
View file

@ -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)

29
openmc/nuclide.py
View file

@ -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
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

2
openmc/surface.py
View file

@ -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):

2
openmc/tallies.py
View file

@ -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

6
openmc/universe.py
View file

@ -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):

57
src/input_xml.F90
View file

@ -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 <sab> 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)

3
src/material_header.F90
View file

@ -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

50
src/physics.F90
View file

@ -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
!===============================================================================

4
src/relaxng/materials.rnc
View file

@ -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 })

46
src/relaxng/materials.rng
View file

@ -12,6 +12,20 @@
<data type="int"/>
</attribute>
</choice>
<optional>
<choice>
<element name="name">
<data type="string">
<param name="maxLength">52</param>
</data>
</element>
<attribute name="name">
<data type="string">
<param name="maxLength">52</param>
</data>
</attribute>
</choice>
</optional>
<element name="density">
<interleave>
<optional>
@ -67,6 +81,22 @@
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
</optional>
<choice>
<choice>
<element name="ao">
@ -117,6 +147,22 @@
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
</optional>
<choice>
<choice>
<element name="ao">