crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #897 from smharper/partial_sab

Browse source 
Partial S(alpha, beta) tables
This commit is contained in:
Paul Romano 2017-07-13 15:02:11 -05:00 • committed by GitHub
commit b6a91cc25a
15 changed files with 326 additions and 179 deletions
Download patch file Download diff file Expand all files Collapse all files

6
docs/source/io_formats/materials.rst
View file

@ -107,9 +107,13 @@ Each ``material`` element can have the following attributes or sub-elements:
multi-group :ref:`energy_mode`.
:sab:
Associates an S(a,b) table with the material. This element has one
Associates an S(a,b) table with the material. This element has an
attribute/sub-element called ``name``. The ``name`` attribute
is the name of the S(a,b) table that should be associated with the material.
There is also an optional ``fraction`` element which indicates what fraction
of the relevant nuclides will be affected by the S(a,b) table (e.g. which
fraction of a material is crystalline versus amorphous). ``fraction``
defaults to unity.
*Default*: None

17
openmc/material.py
View file

@ -618,13 +618,18 @@ class Material(IDManagerMixin):
if element == elm[0]:
self._elements.remove(elm)
def add_s_alpha_beta(self, name):
def add_s_alpha_beta(self, name, fraction=1.0):
r"""Add an :math:`S(\alpha,\beta)` table to the material
Parameters
----------
name : str
Name of the :math:`S(\alpha,\beta)` table
fraction : float
The fraction of relevant nuclei that are affected by the
:math:`S(\alpha,\beta)` table. For example, if the material is a
block of carbon that is 60% graphite and 40% amorphous then add a
graphite :math:`S(\alpha,\beta)` table with fraction=0.6.
"""
@ -638,13 +643,17 @@ class Material(IDManagerMixin):
'non-string table name "{}"'.format(self._id, name)
raise ValueError(msg)
cv.check_type('S(a,b) fraction', fraction, Real)
cv.check_greater_than('S(a,b) fraction', fraction, 0.0, True)
cv.check_less_than('S(a,b) fraction', fraction, 1.0, True)
new_name = openmc.data.get_thermal_name(name)
if new_name != name:
msg = 'OpenMC S(a,b) tables follow the GND naming convention. ' \
'Table "{}" is being renamed as "{}".'.format(name, new_name)
warnings.warn(msg)
self._sab.append(new_name)
self._sab.append((new_name, fraction))
def make_isotropic_in_lab(self):
for nuclide, percent, percent_type in self._nuclides:
@ -972,7 +981,9 @@ class Material(IDManagerMixin):
if len(self._sab) > 0:
for sab in self._sab:
subelement = ET.SubElement(element, "sab")
subelement.set("name", sab)
subelement.set("name", sab[0])
if sab[1] != 1.0:
subelement.set("fraction", str(sab[1]))
return element

97
src/cross_section.F90
View file

@ -36,6 +36,7 @@ contains
integer :: i_grid ! index into logarithmic mapping array or material
! union grid
real(8) :: atom_density ! atom density of a nuclide
real(8) :: sab_frac ! fraction of atoms affected by S(a,b)
logical :: check_sab ! should we check for S(a,b) table?
! Set all material macroscopic cross sections to zero
@ -60,21 +61,25 @@ contains
! Add contribution from each nuclide in material
do i = 1, mat % n_nuclides
! ========================================================================
! ======================================================================
! CHECK FOR S(A,B) TABLE
i_sab = 0
sab_frac = ZERO
! Check if this nuclide matches one of the S(a,b) tables specified -- this
! relies on i_sab_nuclides being in sorted order
! Check if this nuclide matches one of the S(a,b) tables specified.
! This relies on i_sab_nuclides being in sorted order
if (check_sab) then
if (i == mat % i_sab_nuclides(j)) then
! Get index in sab_tables
i_sab = mat % i_sab_tables(j)
sab_frac = mat % sab_fracs(j)
! If particle energy is greater than the highest energy for the S(a,b)
! table, don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) i_sab = 0
! If particle energy is greater than the highest energy for the
! S(a,b) table, then don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then
i_sab = 0
end if
! Increment position in i_sab_nuclides
j = j + 1
@ -84,7 +89,7 @@ contains
end if
end if
! ========================================================================
! ======================================================================
! CALCULATE MICROSCOPIC CROSS SECTION
! Determine microscopic cross sections for this nuclide
@ -92,13 +97,14 @@ contains
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E &
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT)
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT &
.or. i_sab /= micro_xs(i_nuclide) % index_sab &
.or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, &
p % sqrtkT, sab_frac)
end if
! ========================================================================
! ======================================================================
! ADD TO MACROSCOPIC CROSS SECTION
! Copy atom density of nuclide in material
@ -133,13 +139,15 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT, &
sab_frac)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer, intent(in) :: i_log_union ! index into logarithmic mapping array or
! material union energy grid
real(8), intent(in) :: sqrtkT ! Square root of kT, material dependent
real(8), intent(in) :: sqrtkT ! square root of kT, material dependent
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
logical :: use_mp ! true if XS can be calculated with windowed multipole
integer :: i_temp ! index for temperature
@ -243,8 +251,10 @@ contains
micro_xs(i_nuclide) % interp_factor = f
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % thermal = ZERO
micro_xs(i_nuclide) % thermal_elastic = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) &
@ -271,20 +281,21 @@ contains
end if
! Initialize sab treatment to false
micro_xs(i_nuclide) % index_sab = NONE
micro_xs(i_nuclide) % elastic_sab = ZERO
micro_xs(i_nuclide) % index_sab = NONE
micro_xs(i_nuclide) % sab_frac = ZERO
! Initialize URR probability table treatment to false
micro_xs(i_nuclide) % use_ptable = .false.
micro_xs(i_nuclide) % use_ptable = .false.
! If there is S(a,b) data for this nuclide, we need to do a few
! things. Since the total cross section was based on non-S(a,b) data, we
! need to correct it by subtracting the non-S(a,b) elastic cross section and
! then add back in the calculated S(a,b) elastic+inelastic cross section.
! If there is S(a,b) data for this nuclide, we need to set the sab_scatter
! and sab_elastic cross sections and correct the total and elastic cross
! sections.
if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT)
if (i_sab > 0) then
call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT, sab_frac)
end if
! if the particle is in the unresolved resonance range and there are
! If the particle is in the unresolved resonance range and there are
! probability tables, we need to determine cross sections from the table
if (urr_ptables_on .and. nuc % urr_present .and. .not. use_mp) then
@ -295,7 +306,6 @@ contains
end if
micro_xs(i_nuclide) % last_E = E
micro_xs(i_nuclide) % last_index_sab = i_sab
micro_xs(i_nuclide) % last_sqrtkT = sqrtkT
end associate
@ -303,16 +313,16 @@ contains
!===============================================================================
! CALCULATE_SAB_XS determines the elastic and inelastic scattering
! cross-sections in the thermal energy range. These cross sections replace
! whatever data were taken from the normal Nuclide table.
! cross-sections in the thermal energy range. These cross sections replace a
! fraction of whatever data were taken from the normal Nuclide table.
!===============================================================================
subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT)
subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT, sab_frac)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
real(8), intent(in) :: sqrtkT ! temperature
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
integer :: i_grid ! index on S(a,b) energy grid
integer :: i_temp ! temperature index
@ -341,7 +351,8 @@ contains
! Randomly sample between temperature i and i+1
f = (kT - sab_tables(i_sab) % kTs(i_temp)) / &
(sab_tables(i_sab) % kTs(i_temp + 1) - sab_tables(i_sab) % kTs(i_temp))
(sab_tables(i_sab) % kTs(i_temp + 1) &
- sab_tables(i_sab) % kTs(i_temp))
if (f > prn()) i_temp = i_temp + 1
end if
@ -402,16 +413,20 @@ contains
end if
end associate
! Store the S(a,b) cross sections.
micro_xs(i_nuclide) % thermal = sab_frac * (elastic + inelastic)
micro_xs(i_nuclide) % thermal_elastic = sab_frac * elastic
! Correct total and elastic cross sections
micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % elastic + inelastic + elastic
micro_xs(i_nuclide) % elastic = inelastic + elastic
micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total &
+ micro_xs(i_nuclide) % thermal &
- sab_frac * micro_xs(i_nuclide) % elastic
micro_xs(i_nuclide) % elastic = micro_xs(i_nuclide) % thermal &
+ (ONE - sab_frac) * micro_xs(i_nuclide) % elastic
! Store S(a,b) elastic cross section for sampling later
micro_xs(i_nuclide) % elastic_sab = elastic
! Save temperature index
! Save temperature index and thermal fraction
micro_xs(i_nuclide) % index_temp_sab = i_temp
micro_xs(i_nuclide) % sab_frac = sab_frac
end subroutine calculate_sab_xs

35
src/input_xml.F90
View file

@ -2565,6 +2565,7 @@ contains
! table is indeed applied to multiple nuclides.
allocate(mat % sab_names(n_sab))
allocate(mat % i_sab_tables(n_sab))
allocate(mat % sab_fracs(n_sab))
do j = 1, n_sab
! Get pointer to S(a,b) table
@ -2578,6 +2579,13 @@ contains
name = trim(name)
mat % sab_names(j) = name
! Read the fraction of nuclei affected by this S(a,b) table
if (check_for_node(node_sab, "fraction")) then
call get_node_value(node_sab, "fraction", mat % sab_fracs(j))
else
mat % sab_fracs(j) = ONE
end if
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. library_dict % has_key(to_lower(name))) then
call fatal_error("Could not find S(a,b) table " // trim(name) &
@ -5159,9 +5167,11 @@ contains
integer :: m ! position for sorting
integer :: temp_nuclide ! temporary value for sorting
integer :: temp_table ! temporary value for sorting
real(8) :: temp_frac ! temporary value for sorting
logical :: found
type(VectorInt) :: i_sab_tables
type(VectorInt) :: i_sab_nuclides
type(VectorInt) :: i_sab_tables
type(VectorInt) :: i_sab_nuclides
type(VectorReal) :: sab_fracs
do i = 1, size(materials)
! Skip materials with no S(a,b) tables
@ -5179,6 +5189,7 @@ contains
if (any(sab % nuclides == nuclides(mat % nuclide(j)) % name)) then
call i_sab_tables % push_back(mat % i_sab_tables(k))
call i_sab_nuclides % push_back(j)
call sab_fracs % push_back(mat % sab_fracs(k))
found = .true.
end if
end do FIND_NUCLIDE
@ -5192,17 +5203,34 @@ contains
end if
end do ASSIGN_SAB
! Make sure each nuclide only appears in one table.
do j = 1, i_sab_nuclides % size()
do k = j+1, i_sab_nuclides % size()
if (i_sab_nuclides % data(j) == i_sab_nuclides % data(k)) then
call fatal_error(trim( &
nuclides(mat % nuclide(i_sab_nuclides % data(j))) % name) &
// " in material " // trim(to_str(mat % id)) // " was found &
&in multiple S(a,b) tables. Each nuclide can only appear in &
&one S(a,b) table per material.")
end if
end do
end do
! Update i_sab_tables and i_sab_nuclides
deallocate(mat % i_sab_tables)
deallocate(mat % sab_fracs)
m = i_sab_tables % size()
allocate(mat % i_sab_tables(m))
allocate(mat % i_sab_nuclides(m))
allocate(mat % sab_fracs(m))
mat % i_sab_tables(:) = i_sab_tables % data(1:m)
mat % i_sab_nuclides(:) = i_sab_nuclides % data(1:m)
mat % sab_fracs(:) = sab_fracs % data(1:m)
! Clear entries in vectors for next material
call i_sab_tables % clear()
call i_sab_nuclides % clear()
call sab_fracs % clear()
! If there are multiple S(a,b) tables, we need to make sure that the
! entries in i_sab_nuclides are sorted or else they won't be applied
@ -5215,6 +5243,7 @@ contains
m = k
temp_nuclide = mat % i_sab_nuclides(k)
temp_table = mat % i_sab_tables(k)
temp_frac = mat % i_sab_tables(k)
MOVE_OVER: do
! Check if insertion value is greater than (m-1)th value
@ -5223,6 +5252,7 @@ contains
! Move values over until hitting one that's not larger
mat % i_sab_nuclides(m) = mat % i_sab_nuclides(m-1)
mat % i_sab_tables(m) = mat % i_sab_tables(m-1)
mat % sab_fracs(m) = mat % sab_fracs(m-1)
m = m - 1
! Exit if we've reached the beginning of the list
@ -5232,6 +5262,7 @@ contains
! Put the original value into its new position
mat % i_sab_nuclides(m) = temp_nuclide
mat % i_sab_tables(m) = temp_table
mat % sab_fracs(m) = temp_frac
end do SORT_SAB
end if

3
src/material_header.F90
View file

@ -22,10 +22,11 @@ module material_header
! Unionized energy grid information
integer, allocatable :: nuclide_grid_index(:,:) ! nuclide e_grid pointers
! S(a,b) data references
! S(a,b) data
integer :: n_sab = 0 ! number of S(a,b) tables
integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide
integer, allocatable :: i_sab_tables(:) ! index in sab_tables
real(8), allocatable :: sab_fracs(:) ! how often to use S(a,b)
! Temporary names read during initialization
character(20), allocatable :: names(:) ! isotope names

44
src/nuclide_header.F90
View file

@ -102,33 +102,35 @@ module nuclide_header
end type Nuclide
!===============================================================================
! NUCLIDEMICROXS contains cached microscopic cross sections for a
! particular nuclide at the current energy
! NUCLIDEMICROXS contains cached microscopic cross sections for a particular
! nuclide at the current energy
!===============================================================================
type NuclideMicroXS
integer :: index_grid ! index on nuclide energy grid
integer :: index_temp ! temperature index for nuclide
real(8) :: last_E = ZERO ! last evaluated energy
real(8) :: interp_factor ! interpolation factor on nuc. energy grid
real(8) :: total ! microscropic total xs
real(8) :: elastic ! microscopic elastic scattering xs
real(8) :: absorption ! microscopic absorption xs
real(8) :: fission ! microscopic fission xs
real(8) :: nu_fission ! microscopic production xs
! Microscopic cross sections in barns
real(8) :: total
real(8) :: elastic ! If sab_frac is not 1 or 0, then this value is
! averaged over bound and non-bound nuclei
real(8) :: absorption
real(8) :: fission
real(8) :: nu_fission
real(8) :: thermal ! Bound thermal elastic & inelastic scattering
real(8) :: thermal_elastic ! Bound thermal elastic scattering
! Information for S(a,b) use
integer :: index_sab ! index in sab_tables (zero means no table)
integer :: last_index_sab = 0 ! index in sab_tables last used by this nuclide
integer :: index_temp_sab ! temperature index for sab_tables
real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table
! Indicies and factors needed to compute cross sections from the data tables
integer :: index_grid ! Index on nuclide energy grid
integer :: index_temp ! Temperature index for nuclide
real(8) :: interp_factor ! Interpolation factor on nuc. energy grid
integer :: index_sab = NONE ! Index in sab_tables
integer :: index_temp_sab ! Temperature index for sab_tables
real(8) :: sab_frac ! Fraction of atoms affected by S(a,b)
logical :: use_ptable ! In URR range with probability tables?
! Information for URR probability table use
logical :: use_ptable ! in URR range with probability tables?
! Information for Doppler broadening
! Energy and temperature last used to evaluate these cross sections. If
! these values have changed, then the cross sections must be re-evaluated.
real(8) :: last_E = ZERO ! Last evaluated energy
real(8) :: last_sqrtkT = ZERO ! Last temperature in sqrt(Boltzmann
! constant * temperature (eV))
! constant * temperature (eV))
end type NuclideMicroXS
!===============================================================================

62
src/physics.F90
View file

@ -324,6 +324,7 @@ contains
real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering
real(8) :: phi ! azimuthal angle for iso-in-lab scattering
real(8) :: kT ! temperature in eV
logical :: sampled ! whether or not a reaction type has been sampled
type(Nuclide), pointer :: nuc
! copy incoming direction
@ -337,36 +338,43 @@ contains
! For tallying purposes, this routine might be called directly. In that
! case, we need to sample a reaction via the cutoff variable
prob = ZERO
cutoff = prn() * (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption)
sampled = .false.
prob = prob + micro_xs(i_nuclide) % elastic
prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal
if (prob > cutoff) then
! =======================================================================
! 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)
! NON-S(A,B) ELASTIC SCATTERING
! Determine temperature
if (nuc % mp_present) then
kT = p % sqrtkT**2
else
! Determine temperature
if (nuc % mp_present) then
kT = p % sqrtkT**2
else
kT = nuc % kTs(micro_xs(i_nuclide) % index_temp)
end if
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, nuc % reactions(1), kT, &
p % E, p % coord(1) % uvw, p % mu, p % wgt)
kT = nuc % kTs(micro_xs(i_nuclide) % index_temp)
end if
p % event_MT = ELASTIC
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, nuc % reactions(1), kT, p % E, &
p % coord(1) % uvw, p % mu, p % wgt)
else
p % event_MT = ELASTIC
sampled = .true.
end if
prob = micro_xs(i_nuclide) % elastic
if (prob > cutoff .and. .not. sampled) then
! =======================================================================
! S(A,B) SCATTERING
call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, p % E, &
p % coord(1) % uvw, p % mu)
p % event_MT = ELASTIC
sampled = .true.
end if
if (.not. sampled) then
! =======================================================================
! INELASTIC SCATTERING
@ -388,7 +396,7 @@ contains
if (rx % MT == N_FISSION .or. rx % MT == N_F .or. rx % MT == N_NF &
.or. rx % MT == N_2NF .or. rx % MT == N_3NF) cycle
! some materials have gas production cross sections with MT > 200 that
! Some materials have gas production cross sections with MT > 200 that
! are duplicates. Also MT=4 is total level inelastic scattering which
! should be skipped
if (rx % MT >= 200 .or. rx % MT == N_LEVEL) cycle
@ -406,24 +414,24 @@ contains
! Perform collision physics for inelastic scattering
call inelastic_scatter(nuc, nuc%reactions(i), p)
p % event_MT = nuc%reactions(i)%MT
p % event_MT = nuc % reactions(i) % MT
end if
! Set event component
p % event = EVENT_SCATTER
! sample new outgoing angle for isotropic in lab scattering
! Sample new outgoing angle for isotropic in lab scattering
if (materials(p % material) % p0(i_nuc_mat)) then
! sample isotropic-in-lab outgoing direction
! 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
! Change direction of particle
p % coord(1) % uvw = uvw_new
end if
@ -559,8 +567,8 @@ contains
associate (sab => sab_tables(i_sab) % data(i_temp))
! Determine whether inelastic or elastic scattering will occur
if (prn() < micro_xs(i_nuclide) % elastic_sab / &
micro_xs(i_nuclide) % elastic) then
if (prn() < micro_xs(i_nuclide) % thermal_elastic / &
micro_xs(i_nuclide) % thermal) then
! elastic scattering
! Get index and interpolation factor for elastic grid

3
src/relaxng/materials.rnc
View file

@ -28,7 +28,8 @@ element materials {
}* &
element sab {
(element name { xsd:string } | attribute name { xsd:string })
(element name { xsd:string } | attribute name { xsd:string }) &
(element fraction { xsd:double } | attribute fraction { xsd:double })?
}*
}+ &

28
src/relaxng/materials.rng
View file

@ -119,14 +119,26 @@
</zeroOrMore>
<zeroOrMore>
<element name="sab">
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<interleave>
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="fraction">
<data type="double"/>
</element>
<attribute name="fraction">
<data type="double"/>
</attribute>
</choice>
</optional>
</interleave>
</element>
</zeroOrMore>
</interleave>

15
tests/test_salphabeta/geometry.xml
View file

@ -1,15 +0,0 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="x-plane" coeffs="-10" boundary="vacuum"/>
<surface id="2" type="x-plane" coeffs="-5" />
<surface id="3" type="x-plane" coeffs="0" />
<surface id="4" type="x-plane" coeffs="5" />
<surface id="5" type="x-plane" coeffs="10" boundary="vacuum"/>
<cell id="1" material="1" region="1 -2" />
<cell id="2" material="2" region="2 -3" />
<cell id="3" material="3" region="3 -4" />
<cell id="4" material="4" region="4 -5" />
</geometry>

60
tests/test_salphabeta/inputs_true.dat Normal file
View file

@ -0,0 +1,60 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2" universe="1" />
<cell id="2" material="2" region="2 -3" universe="1" />
<cell id="3" material="3" region="3 -4" universe="1" />
<cell id="4" material="4" region="4 -5" universe="1" />
<surface boundary="vacuum" coeffs="-10" id="1" type="x-plane" />
<surface coeffs="-5" id="2" type="x-plane" />
<surface coeffs="0" id="3" type="x-plane" />
<surface coeffs="5" id="4" type="x-plane" />
<surface boundary="vacuum" coeffs="10" id="5" type="x-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material id="1">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
<nuclide ao="1.0" name="H1" />
<sab fraction="0.5" name="c_H_in_H2O" />
</material>
<material id="2">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
<nuclide ao="1.0" name="C0" />
<sab name="c_Graphite" />
</material>
<material id="3">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
<nuclide ao="1.0" name="Be9" />
<nuclide ao="1.0" name="O16" />
<sab name="c_Be_in_BeO" />
<sab name="c_O_in_BeO" />
</material>
<material id="4">
<density units="g/cm3" value="5.90168" />
<nuclide ao="0.3" name="H1" />
<nuclide ao="0.15" name="Zr90" />
<nuclide ao="0.1" name="Zr91" />
<nuclide ao="0.1" name="Zr92" />
<nuclide ao="0.05" name="Zr94" />
<nuclide ao="0.05" name="Zr96" />
<nuclide ao="0.1" name="U235" />
<nuclide ao="0.15" name="U238" />
<sab name="c_Zr_in_ZrH" />
<sab name="c_H_in_ZrH" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>5</batches>
<inactive>0</inactive>
<source strength="1.0">
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

41
tests/test_salphabeta/materials.xml
View file

@ -1,41 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
<nuclide name="H1" ao="1.0" />
<sab name="c_H_in_H2O" />
</material>
<material id="2">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
<nuclide name="C0" ao="1.0" />
<sab name="c_Graphite" />
</material>
<material id="3">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
<nuclide name="Be9" ao="1.0" />
<nuclide name="O16" ao="1.0" />
<sab name="c_Be_in_BeO" />
<sab name="c_O_in_BeO" />
</material>
<material id="4">
<density value="5.90168" units="g/cm3" />
<nuclide name="H1" ao="0.3" />
<nuclide name="Zr90" ao="0.15" />
<nuclide name="Zr91" ao="0.1" />
<nuclide name="Zr92" ao="0.1" />
<nuclide name="Zr94" ao="0.05" />
<nuclide name="Zr96" ao="0.05" />
<nuclide name="U235" ao="0.1" />
<nuclide name="U238" ao="0.15" />
<sab name="c_Zr_in_ZrH" />
<sab name="c_H_in_ZrH" />
</material>
</materials>

2
tests/test_salphabeta/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
8.544160E-01 1.274133E-02
8.447580E-01 1.806149E-02

15
tests/test_salphabeta/settings.xml
View file

@ -1,15 +0,0 @@
<?xml version="1.0"?>
<settings>
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

77
tests/test_salphabeta/test_salphabeta.py
View file

@ -3,9 +3,82 @@
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
from testing_harness import PyAPITestHarness
import openmc
import openmc.model
def make_model():
model = openmc.model.Model()
# Materials
m1 = openmc.Material()
m1.set_density('g/cc', 4.5)
m1.add_nuclide('U235', 1.0)
m1.add_nuclide('H1', 1.0)
m1.add_s_alpha_beta('c_H_in_H2O', fraction=0.5)
m2 = openmc.Material()
m2.set_density('g/cc', 4.5)
m2.add_nuclide('U235', 1.0)
m2.add_nuclide('C0', 1.0)
m2.add_s_alpha_beta('c_Graphite')
m3 = openmc.Material()
m3.set_density('g/cc', 4.5)
m3.add_nuclide('U235', 1.0)
m3.add_nuclide('Be9', 1.0)
m3.add_nuclide('O16', 1.0)
m3.add_s_alpha_beta('c_Be_in_BeO')
m3.add_s_alpha_beta('c_O_in_BeO')
m4 = openmc.Material()
m4.set_density('g/cm3', 5.90168)
m4.add_nuclide('H1', 0.3)
m4.add_nuclide('Zr90', 0.15)
m4.add_nuclide('Zr91', 0.1)
m4.add_nuclide('Zr92', 0.1)
m4.add_nuclide('Zr94', 0.05)
m4.add_nuclide('Zr96', 0.05)
m4.add_nuclide('U235', 0.1)
m4.add_nuclide('U238', 0.15)
m4.add_s_alpha_beta('c_Zr_in_ZrH')
m4.add_s_alpha_beta('c_H_in_ZrH')
model.materials += [m1, m2, m3, m4]
# Geometry
x0 = openmc.XPlane(x0=-10, boundary_type='vacuum')
x1 = openmc.XPlane(x0=-5)
x2 = openmc.XPlane(x0=0)
x3 = openmc.XPlane(x0=5)
x4 = openmc.XPlane(x0=10, boundary_type='vacuum')
root_univ = openmc.Universe()
surfs = (x0, x1, x2, x3, x4)
mats = (m1, m2, m3, m4)
cells = []
for i in range(4):
cell = openmc.Cell()
cell.region = +surfs[i] & -surfs[i+1]
cell.fill = mats[i]
root_univ.add_cell(cell)
model.geometry.root_universe = root_univ
# Settings
model.settings.batches = 5
model.settings.inactive = 0
model.settings.particles = 1000
model.settings.source = openmc.Source(space=openmc.stats.Box(
[-4, -4, -4], [4, 4, 4]))
return model
if __name__ == '__main__':
harness = TestHarness('statepoint.10.h5')
model = make_model()
harness = PyAPITestHarness('statepoint.5.h5', model)
harness.main()