Add diff tally for keff---density or nuc. density

This commit is contained in:
Sterling Harper 2015-11-08 17:24:16 -05:00
parent 39e246ed38
commit 199d4af9c5
22 changed files with 732 additions and 4 deletions

View file

@ -259,7 +259,7 @@ module constants
EVENT_ABSORB = 2
! Tally score type
integer, parameter :: N_SCORE_TYPES = 22
integer, parameter :: N_SCORE_TYPES = 23
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
@ -282,7 +282,8 @@ module constants
SCORE_NU_SCATTER_YN = -19, & ! angular flux-weighted nu-scattering moment (0:N)
SCORE_EVENTS = -20, & ! number of events
SCORE_DELAYED_NU_FISSION = -21, & ! delayed neutron production rate
SCORE_INVERSE_VELOCITY = -22 ! flux-weighted inverse velocity
SCORE_INVERSE_VELOCITY = -22, & ! flux-weighted inverse velocity
SCORE_KEFF = -23 ! multiplication factor
! Maximum scattering order supported
integer, parameter :: MAX_ANG_ORDER = 10
@ -348,6 +349,11 @@ module constants
K_TRACKLENGTH = 3, &
LEAKAGE = 4
! Differential tally dependent variables
integer, parameter :: &
DIFF_DENSITY = 1, &
DIFF_NUCLIDE_DENSITY = 2
! ============================================================================
! RANDOM NUMBER STREAM CONSTANTS

View file

@ -831,6 +831,7 @@ contains
! Change density in g/cm^3 to atom/b-cm. Since all values are now in atom
! percent, the sum needs to be re-evaluated as 1/sum(x*awr)
if (.not. density_in_atom) then
mat % density_gpcc = -mat % density
sum_percent = ZERO
do j = 1, mat%n_nuclides
index_list = xs_listing_dict%get_key(mat%names(j))

View file

@ -2197,6 +2197,7 @@ contains
type(Node), pointer :: node_tal => null()
type(Node), pointer :: node_filt => null()
type(Node), pointer :: node_trigger=>null()
type(Node), pointer :: node_deriv => null()
type(NodeList), pointer :: node_mesh_list => null()
type(NodeList), pointer :: node_tal_list => null()
type(NodeList), pointer :: node_filt_list => null()
@ -2407,7 +2408,7 @@ contains
t % estimator = ESTIMATOR_TRACKLENGTH
! Copy material id
! Copy tally id
if (check_for_node(node_tal, "id")) then
call get_node_value(node_tal, "id", t % id)
else
@ -2843,6 +2844,52 @@ contains
end do
end if
! =======================================================================
! READ DATA FOR DERIVATIVES
if (check_for_node(node_tal, "derivative")) then
call get_node_ptr(node_tal, "derivative", node_deriv)
allocate(t % deriv)
temp_str = ""
call get_node_value(node_deriv, "variable", temp_str)
temp_str = to_lower(temp_str)
select case(temp_str)
case("density")
t % deriv % dep_var = DIFF_DENSITY
call get_node_value(node_deriv, "material", t % deriv % diff_material)
case("nuclide_density")
t % deriv % dep_var = DIFF_NUCLIDE_DENSITY
call get_node_value(node_deriv, "material", t % deriv % diff_material)
call get_node_value(node_deriv, "nuclide", word)
word = trim(to_lower(word))
pair_list => nuclide_dict % keys()
do while (associated(pair_list))
if (starts_with(pair_list % key, word)) then
word = pair_list % key(1:150)
exit
end if
! Advance to next
pair_list => pair_list % next
end do
! Check if no nuclide was found
if (.not. associated(pair_list)) then
call fatal_error("Could not find the nuclide " &
&// trim(word) // " specified in tally " &
&// trim(to_str(t % id)) // " in any material.")
end if
deallocate(pair_list)
t % deriv % diff_nuclide = nuclide_dict % get_key(word)
end select
end if
! =======================================================================
! READ DATA FOR SCORES
@ -2991,6 +3038,11 @@ contains
call fatal_error("Cannot tally flux with an outgoing energy &
&filter.")
end if
if (allocated(t % deriv)) then
t % estimator = ESTIMATOR_COLLISION
end if
case ('flux-yn')
! Prohibit user from tallying flux for an individual nuclide
if (.not. (t % n_nuclide_bins == 1 .and. &
@ -3188,6 +3240,12 @@ contains
case ('events')
t % score_bins(j) = SCORE_EVENTS
case ('keff')
t % score_bins(j) = SCORE_KEFF
if (allocated(t % deriv)) then
t % estimator = ESTIMATOR_COLLISION
end if
case ('elastic', '(n,elastic)')
t % score_bins(j) = ELASTIC
case ('(n,2nd)')

View file

@ -13,6 +13,7 @@ module material_header
integer, allocatable :: nuclide(:) ! index in nuclides array
real(8) :: density ! total atom density in atom/b-cm
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
real(8) :: density_gpcc ! total density in g/cm^3
! Energy grid information
integer :: n_grid ! # of union material grid points

View file

@ -986,6 +986,7 @@ contains
score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment"
score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate"
score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity"
score_names(abs(SCORE_KEFF)) = "k_eff, multiplaction factor"
! Create filename for tally output
filename = trim(path_output) // "tallies.out"

View file

@ -654,6 +654,13 @@ contains
end if
end if
case (SCORE_KEFF)
if (t % estimator == ESTIMATOR_COLLISION) then
score = p % last_wgt * material_xs % nu_fission / material_xs % total
else if (t % estimator == ESTIMATOR_TRACKLENGTH) then
score = flux * material_xs % nu_fission
end if
case default
if (t % estimator == ESTIMATOR_ANALOG) then
! Any other score is assumed to be a MT number. Thus, we just need
@ -732,6 +739,35 @@ contains
end select
!#########################################################################
! Add derivative information on score for differential tallies.
if (allocated(t % deriv) .and. t % estimator == ESTIMATOR_COLLISION) then
select case (score_bin)
case (SCORE_KEFF)
select case (t % deriv % dep_var)
case (DIFF_DENSITY)
score = score * t % deriv % accumulator
case (DIFF_NUCLIDE_DENSITY)
mat => materials(p % material)
if (mat % id == t % deriv % diff_material .and. &
micro_xs(t % deriv % diff_nuclide) % nu_fission /= ZERO) then
do l = 1, mat % n_nuclides
if (mat % nuclide(l) == t % deriv % diff_nuclide) exit
end do
score = score * (t % deriv % accumulator + ONE &
/ mat % atom_density(l))
else
score = score * t % deriv % accumulator
end if
end select
end select
end if
!#########################################################################
! Expand score if necessary and add to tally results.
@ -2210,6 +2246,80 @@ contains
end subroutine score_surface_current
!===============================================================================
! SCORE_DIFF_ACCUMULATORS
!===============================================================================
subroutine score_diff_accumulators(p, distance, collision)
type(particle), intent(in) :: p
real(8), intent(in) :: distance
logical, intent(in) :: collision
integer :: i, j
type(TallyObject), pointer :: t
type(Material), pointer :: mat
do i = 1, n_user_tallies
t => user_tallies(i)
if (allocated(t % deriv)) then
select case (t % deriv % dep_var)
case (DIFF_DENSITY)
if (p % material == MATERIAL_VOID) cycle
mat => materials(p % material)
if (mat % id == t % deriv % diff_material) then
if (collision) then
t % deriv % accumulator = t % deriv % accumulator &
+ ONE / mat % density_gpcc &
- distance * material_xs % total / mat % density_gpcc
else
t % deriv % accumulator = t % deriv % accumulator &
- distance * material_xs % total / mat % density_gpcc
end if
end if
case (DIFF_NUCLIDE_DENSITY)
if (p % material == MATERIAL_VOID) cycle
mat => materials(p % material)
if (mat % id == t % deriv % diff_material) then
do j = 1, mat % n_nuclides
if (mat % nuclide(j) == t % deriv % diff_nuclide) exit
end do
if (mat % nuclide(j) /= t % deriv % diff_nuclide) then
call fatal_error("Couldn't find the right nuclide.")
end if
t % deriv % accumulator = t % deriv % accumulator &
- distance * micro_xs(t % deriv % diff_nuclide) % total
if (collision) then
if (p % event == EVENT_SCATTER) then
t % deriv % accumulator = t % deriv % accumulator &
+ micro_xs(t % deriv % diff_nuclide) % elastic &
/ material_xs % elastic
else if (p % event == EVENT_ABSORB) then
t % deriv % accumulator = t % deriv % accumulator &
+ micro_xs(t % deriv % diff_nuclide) % absorption &
/ material_xs % absorption
end if
end if
end if
end select
end if
end do
end subroutine
subroutine clear_diff_accumulators()
integer :: i
type(TallyObject), pointer :: t
do i = 1, n_user_tallies
t => user_tallies(i)
if (allocated(t % deriv)) then
t % deriv % accumulator = ZERO
end if
end do
end subroutine
!===============================================================================
! GET_NEXT_BIN determines the next scoring bin for a particular filter variable
!===============================================================================

View file

@ -64,6 +64,18 @@ module tally_header
procedure :: clear => tallyfilter_clear ! Deallocates TallyFilter
end type TallyFilter
!===============================================================================
! TALLYDERIVATIVE
!===============================================================================
type TallyDerivative
real(8) :: accumulator
integer :: dep_var
integer :: diff_material
integer :: diff_nuclide
end type TallyDerivative
!===============================================================================
! TALLYOBJECT describes a user-specified tally. The region of phase space to
! tally in is given by the TallyFilters and the results are stored in a
@ -130,6 +142,9 @@ module tally_header
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
! Derivative for differentially tallies
type(TallyDerivative), allocatable :: deriv
! Type-Bound procedures
contains
procedure :: clear => tallyobject_clear ! Deallocates TallyObject
@ -204,6 +219,9 @@ module tally_header
this % n_triggers = 0
if (allocated(this % deriv)) &
deallocate (this % deriv)
end subroutine tallyobject_clear
end module tally_header

View file

@ -13,7 +13,8 @@ module tracking
use random_lcg, only: prn
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_current
score_collision_tally, score_surface_current, &
score_diff_accumulators, clear_diff_accumulators
use track_output, only: initialize_particle_track, write_particle_track, &
add_particle_track, finalize_particle_track
@ -60,6 +61,9 @@ contains
call initialize_particle_track()
endif
! Every particle starts with no accumulated flux derivative.
call clear_diff_accumulators()
EVENT_LOOP: do
! If the cell hasn't been determined based on the particle's location,
! initiate a search for the current cell. This generally happens at the
@ -114,6 +118,9 @@ contains
distance * material_xs % nu_fission
end if
! Score flux derivative accumulators for differential tallies.
call score_diff_accumulators(p, distance, d_boundary > d_collision)
if (d_collision > d_boundary) then
! ====================================================================
! PARTICLE CROSSES SURFACE

View file

@ -0,0 +1,150 @@
import openmc
def make_mats(mod_density=0.7420582):
water = openmc.Material(name='water')
water.set_density('g/cm3', mod_density)
water.add_nuclide(openmc.Nuclide('H-1'), 2.0)
water.add_nuclide(openmc.Nuclide('O-16'), 1.0)
fuel = openmc.Material(name='fuel 2.4%')
fuel.set_density('g/cm3', 10.29769)
fuel.add_nuclide(openmc.Nuclide('U-234'), 5.7987e-06)
fuel.add_nuclide(openmc.Nuclide('U-235'), 7.2175e-04)
fuel.add_nuclide(openmc.Nuclide('U-238'), 2.2253e-02)
fuel.add_nuclide(openmc.Nuclide('O-16'), 4.5750e-02)
fuel.add_nuclide(openmc.Nuclide('O-17'), 9.4222e-05)
zirc4 = openmc.Material(name='zircaloy 4')
zirc4.set_density('g/cm3', 6.55)
zirc4.add_nuclide(openmc.Nuclide('O-16'), 3.0743e-04)
zirc4.add_nuclide(openmc.Nuclide('O-17'), 1.1711e-07)
# O-18 omitted
zirc4.add_nuclide(openmc.Nuclide('Cr-50'), 3.2962e-06)
zirc4.add_nuclide(openmc.Nuclide('Cr-52'), 6.3564e-05)
zirc4.add_nuclide(openmc.Nuclide('Cr-53'), 7.2076e-06)
zirc4.add_nuclide(openmc.Nuclide('Cr-54'), 1.7941e-06)
zirc4.add_nuclide(openmc.Nuclide('Fe-54'), 8.6699e-06)
zirc4.add_nuclide(openmc.Nuclide('Fe-56'), 1.3610e-04)
zirc4.add_nuclide(openmc.Nuclide('Fe-57'), 3.1431e-06)
zirc4.add_nuclide(openmc.Nuclide('Fe-58'), 4.1829e-07)
zirc4.add_nuclide(openmc.Nuclide('Zr-90'), 2.1827e-02)
zirc4.add_nuclide(openmc.Nuclide('Zr-91'), 4.7600e-03)
zirc4.add_nuclide(openmc.Nuclide('Zr-92'), 7.2758e-03)
zirc4.add_nuclide(openmc.Nuclide('Zr-94'), 7.3734e-03)
zirc4.add_nuclide(openmc.Nuclide('Zr-96'), 1.1879e-03)
zirc4.add_nuclide(openmc.Nuclide('Sn-112'), 4.6735e-06)
zirc4.add_nuclide(openmc.Nuclide('Sn-114'), 3.1799e-06)
zirc4.add_nuclide(openmc.Nuclide('Sn-115'), 1.6381e-06)
zirc4.add_nuclide(openmc.Nuclide('Sn-116'), 7.0055e-05)
zirc4.add_nuclide(openmc.Nuclide('Sn-117'), 3.7003e-05)
zirc4.add_nuclide(openmc.Nuclide('Sn-118'), 1.1669e-04)
zirc4.add_nuclide(openmc.Nuclide('Sn-119'), 4.1387e-05)
zirc4.add_nuclide(openmc.Nuclide('Sn-120'), 1.5697e-04)
zirc4.add_nuclide(openmc.Nuclide('Sn-122'), 2.2308e-05)
zirc4.add_nuclide(openmc.Nuclide('Sn-124'), 2.7897e-05)
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '71c'
materials_file.add_materials([water, fuel, zirc4])
return (materials_file,
{'mod':water,
'fuel':fuel,
'clad':zirc4})
def make_geom(mats):
# Instantiate surfaces.
rfo = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.39218, name='fuel outer')
rci = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.40005, name='clad inner')
rco = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.45720, name='clad outer')
x0 = openmc.XPlane(x0=-0.62992)
x1 = openmc.XPlane(x0=0.62992)
y0 = openmc.YPlane(y0=-0.62992)
y1 = openmc.YPlane(y0=0.62992)
z0 = openmc.ZPlane(z0=-1.0)
z1 = openmc.ZPlane(z0=1.0)
x0.boundary_type = 'reflective'
x1.boundary_type = 'reflective'
y0.boundary_type = 'reflective'
y1.boundary_type = 'reflective'
z0.boundary_type = 'reflective'
z1.boundary_type = 'reflective'
# Instantiate cells.
fuel_c = openmc.Cell(name='fuel')
fuel_c.add_surface(rfo, halfspace=-1)
fuel_c.add_surface(z0, halfspace=+1)
fuel_c.add_surface(z1, halfspace=-1)
gap = openmc.Cell(name='gap')
gap.add_surface(rci, halfspace=-1)
gap.add_surface(rfo, halfspace=+1)
gap.add_surface(z0, halfspace=+1)
gap.add_surface(z1, halfspace=-1)
clad_c = openmc.Cell(name='clad')
clad_c.add_surface(rco, halfspace=-1)
clad_c.add_surface(rci, halfspace=+1)
clad_c.add_surface(z0, halfspace=+1)
clad_c.add_surface(z1, halfspace=-1)
mod_c = openmc.Cell(name='moderator')
mod_c.add_surface(rco, halfspace=+1)
mod_c.add_surface(x0, halfspace=+1)
mod_c.add_surface(x1, halfspace=-1)
mod_c.add_surface(y0, halfspace=+1)
mod_c.add_surface(y1, halfspace=-1)
mod_c.add_surface(z0, halfspace=+1)
mod_c.add_surface(z1, halfspace=-1)
# Add materials to cells.
fuel_c.fill = mats['fuel']
gap.fill = 'void'
clad_c.fill = mats['clad']
mod_c.fill = mats['mod']
# Instantiate universes.
u0 = openmc.Universe(universe_id=0)
u0.add_cells([fuel_c, gap, clad_c, mod_c])
# Write the XML file.
geometry = openmc.Geometry()
geometry.root_universe = u0
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
return geometry_file
def make_settings(**kwargs):
if 'batches' in kwargs:
batches = kwargs['batches']
settings_file = openmc.SettingsFile()
settings_file.batches = kwargs.setdefault('batches', 100)
settings_file.inactive = kwargs.setdefault('inactive', 10)
settings_file.particles = kwargs.setdefault('particles', 1000)
settings_file.set_source_space('box', (-0.6, -0.6, -0.9,
0.6, 0.6, 0.9))
settings_file.entropy_lower_left = (-0.7, -0.7, -1.1)
settings_file.entropy_upper_right = (0.7, 0.7, 1.1)
settings_file.entropy_dimension = (10, 10, 10)
return settings_file
def make_inputs(mod_density=0.7420582, **kwargs):
(mat_file, mats) = make_mats(mod_density)
mat_file.export_to_xml()
geo_file = make_geom(mats)
geo_file.export_to_xml()
sets_file = make_settings(**kwargs)
sets_file.export_to_xml()
if __name__ == '__main__':
make_inputs(inactive=5, batches=10, particles=100)

View file

@ -0,0 +1,16 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="10000" material="10001" name="fuel" region="(-10000 10007 -10008)" universe="0" />
<cell id="10001" material="void" name="gap" region="(-10001 10000 10007 -10008)" universe="0" />
<cell id="10002" material="10002" name="clad" region="(-10002 10001 10007 -10008)" universe="0" />
<cell id="10003" material="10000" name="moderator" region="(10002 10003 -10004 10005 -10006 10007 -10008)" universe="0" />
<surface boundary="transmission" coeffs="0.0 0.0 0.39218" id="10000" name="fuel outer" type="z-cylinder" />
<surface boundary="transmission" coeffs="0.0 0.0 0.40005" id="10001" name="clad inner" type="z-cylinder" />
<surface boundary="transmission" coeffs="0.0 0.0 0.4572" id="10002" name="clad outer" type="z-cylinder" />
<surface boundary="reflective" coeffs="-0.62992" id="10003" type="x-plane" />
<surface boundary="reflective" coeffs="0.62992" id="10004" type="x-plane" />
<surface boundary="reflective" coeffs="-0.62992" id="10005" type="y-plane" />
<surface boundary="reflective" coeffs="0.62992" id="10006" type="y-plane" />
<surface boundary="reflective" coeffs="-1.0" id="10007" type="z-plane" />
<surface boundary="reflective" coeffs="1.0" id="10008" type="z-plane" />
</geometry>

View file

@ -0,0 +1,45 @@
<?xml version='1.0' encoding='utf-8'?>
<materials>
<default_xs>71c</default_xs>
<material id="10000" name="water">
<density units="g/cm3" value="0.7420582" />
<nuclide ao="2.0" name="H-1" />
<nuclide ao="1.0" name="O-16" />
</material>
<material id="10001" name="fuel 2.4%">
<density units="g/cm3" value="10.29769" />
<nuclide ao="5.7987e-06" name="U-234" />
<nuclide ao="0.00072175" name="U-235" />
<nuclide ao="0.022253" name="U-238" />
<nuclide ao="0.04575" name="O-16" />
<nuclide ao="9.4222e-05" name="O-17" />
</material>
<material id="10002" name="zircaloy 4">
<density units="g/cm3" value="6.55" />
<nuclide ao="0.00030743" name="O-16" />
<nuclide ao="1.1711e-07" name="O-17" />
<nuclide ao="3.2962e-06" name="Cr-50" />
<nuclide ao="6.3564e-05" name="Cr-52" />
<nuclide ao="7.2076e-06" name="Cr-53" />
<nuclide ao="1.7941e-06" name="Cr-54" />
<nuclide ao="8.6699e-06" name="Fe-54" />
<nuclide ao="0.0001361" name="Fe-56" />
<nuclide ao="3.1431e-06" name="Fe-57" />
<nuclide ao="4.1829e-07" name="Fe-58" />
<nuclide ao="0.021827" name="Zr-90" />
<nuclide ao="0.00476" name="Zr-91" />
<nuclide ao="0.0072758" name="Zr-92" />
<nuclide ao="0.0073734" name="Zr-94" />
<nuclide ao="0.0011879" name="Zr-96" />
<nuclide ao="4.6735e-06" name="Sn-112" />
<nuclide ao="3.1799e-06" name="Sn-114" />
<nuclide ao="1.6381e-06" name="Sn-115" />
<nuclide ao="7.0055e-05" name="Sn-116" />
<nuclide ao="3.7003e-05" name="Sn-117" />
<nuclide ao="0.00011669" name="Sn-118" />
<nuclide ao="4.1387e-05" name="Sn-119" />
<nuclide ao="0.00015697" name="Sn-120" />
<nuclide ao="2.2308e-05" name="Sn-122" />
<nuclide ao="2.7897e-05" name="Sn-124" />
</material>
</materials>

View file

@ -0,0 +1,5 @@
k-combined:
1.502605E+00 2.273838E-02
tally 1:
5.280104E+00
1.379258E+01

View file

@ -0,0 +1,18 @@
<?xml version='1.0' encoding='utf-8'?>
<settings>
<eigenvalue>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
</eigenvalue>
<source>
<space type="box">
<parameters>-0.6 -0.6 -0.9 0.6 0.6 0.9</parameters>
</space>
</source>
<entropy>
<dimension>10 10 10</dimension>
<lower_left>-0.7 -0.7 -1.1</lower_left>
<upper_right>0.7 0.7 1.1</upper_right>
</entropy>
</settings>

View file

@ -0,0 +1,7 @@
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<tally id="1">
<scores>keff</scores>
<derivative variable="density" material="10000"/>
</tally>
</tallies>

View file

@ -0,0 +1,11 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*', True)
harness.main()

View file

@ -0,0 +1,164 @@
import openmc
def make_mats(**kwargs):
water_dict = {
'B-10': 8.0042e-06,
'B-11': 3.2218e-05,
'H-1': 4.9457e-02,
'H-2': 7.4196e-06,
'O-16': 2.4672e-02,
'O-17': 9.3982e-06 + 5.0701e-05}
water = openmc.Material(name='water')
water.set_density('sum')
for nuclide in water_dict:
water.add_nuclide(nuclide, water_dict[nuclide])
U235_dens = kwargs.setdefault('U235_dens', 5.5814e-04)
fuel_dict = {
'B-10': 8.0042e-06,
'U-234': 4.4842e-06,
'U-235': U235_dens,
'U-238': 2.2407e-02,
'O-16': 4.5828e-02,
'O-17': 1.7457e-05 + 9.4176e-05}
fuel = openmc.Material(name='fuel 2.4%')
fuel.set_density('sum')
for nuclide in fuel_dict:
fuel.add_nuclide(nuclide, fuel_dict[nuclide])
zirc4_dict = {
'O-16': 3.0743e-04,
'O-17': 1.1711e-07 + 6.3176e-07,
'Cr-50': 3.2962e-06,
'Cr-52': 6.3564e-05,
'Cr-53': 7.2076e-06,
'Cr-54': 1.7941e-06,
'Fe-54': 8.6699e-06,
'Fe-56': 1.3610e-04,
'Fe-57': 3.1431e-06,
'Fe-58': 4.1829e-07,
'Zr-90': 2.1827e-02,
'Zr-91': 4.7600e-03,
'Zr-92': 7.2758e-03,
'Zr-94': 7.3734e-03,
'Zr-96': 1.1879e-03,
'Sn-112': 4.6735e-06,
'Sn-114': 3.1799e-06,
'Sn-115': 1.6381e-06,
'Sn-116': 7.0055e-05,
'Sn-117': 3.7003e-05,
'Sn-118': 1.1669e-04,
'Sn-119': 4.1387e-05,
'Sn-120': 1.5697e-04,
'Sn-122': 2.2308e-05,
'Sn-124': 2.7897e-05}
zirc4 = openmc.Material(name='zircaloy 4')
zirc4.set_density('sum')
for nuclide in zirc4_dict:
zirc4.add_nuclide(nuclide, zirc4_dict[nuclide])
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '71c'
materials_file.add_materials([water, fuel, zirc4])
return (materials_file,
{'mod':water,
'fuel':fuel,
'clad':zirc4})
def make_geom(mats):
# Instantiate surfaces.
rfo = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.39218, name='fuel outer')
rci = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.40005, name='clad inner')
rco = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.45720, name='clad outer')
x0 = openmc.XPlane(x0=-0.62992)
x1 = openmc.XPlane(x0=0.62992)
y0 = openmc.YPlane(y0=-0.62992)
y1 = openmc.YPlane(y0=0.62992)
z0 = openmc.ZPlane(z0=-1.0)
z1 = openmc.ZPlane(z0=1.0)
x0.boundary_type = 'reflective'
x1.boundary_type = 'reflective'
y0.boundary_type = 'reflective'
y1.boundary_type = 'reflective'
z0.boundary_type = 'reflective'
z1.boundary_type = 'reflective'
# Instantiate cells.
fuel_c = openmc.Cell(name='fuel')
fuel_c.add_surface(rfo, halfspace=-1)
fuel_c.add_surface(z0, halfspace=+1)
fuel_c.add_surface(z1, halfspace=-1)
gap = openmc.Cell(name='gap')
gap.add_surface(rci, halfspace=-1)
gap.add_surface(rfo, halfspace=+1)
gap.add_surface(z0, halfspace=+1)
gap.add_surface(z1, halfspace=-1)
clad_c = openmc.Cell(name='clad')
clad_c.add_surface(rco, halfspace=-1)
clad_c.add_surface(rci, halfspace=+1)
clad_c.add_surface(z0, halfspace=+1)
clad_c.add_surface(z1, halfspace=-1)
mod_c = openmc.Cell(name='moderator')
mod_c.add_surface(rco, halfspace=+1)
mod_c.add_surface(x0, halfspace=+1)
mod_c.add_surface(x1, halfspace=-1)
mod_c.add_surface(y0, halfspace=+1)
mod_c.add_surface(y1, halfspace=-1)
mod_c.add_surface(z0, halfspace=+1)
mod_c.add_surface(z1, halfspace=-1)
# Add materials to cells.
fuel_c.fill = mats['fuel']
gap.fill = 'void'
clad_c.fill = mats['clad']
mod_c.fill = mats['mod']
# Instantiate universes.
u0 = openmc.Universe(universe_id=0)
u0.add_cells([fuel_c, gap, clad_c, mod_c])
# Write the XML file.
geometry = openmc.Geometry()
geometry.root_universe = u0
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
return geometry_file
def make_settings(**kwargs):
if 'batches' in kwargs:
batches = kwargs['batches']
settings_file = openmc.SettingsFile()
settings_file.batches = kwargs.setdefault('batches', 1610)
settings_file.inactive = kwargs.setdefault('inactive', 10)
settings_file.particles = kwargs.setdefault('particles', 1000)
settings_file.set_source_space('box', (-0.6, -0.6, -0.9,
0.6, 0.6, 0.9))
settings_file.entropy_lower_left = (-0.7, -0.7, -1.1)
settings_file.entropy_upper_right = (0.7, 0.7, 1.1)
settings_file.entropy_dimension = (10, 10, 10)
return settings_file
def make_inputs(**kwargs):
(mat_file, mats) = make_mats(**kwargs)
mat_file.export_to_xml()
geo_file = make_geom(mats)
geo_file.export_to_xml()
sets_file = make_settings(**kwargs)
sets_file.export_to_xml()
if __name__ == '__main__':
make_inputs(inactive=5, batches=10, particles=100)

View file

@ -0,0 +1,16 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="10000" material="10001" name="fuel" region="(-10000 10007 -10008)" universe="0" />
<cell id="10001" material="void" name="gap" region="(-10001 10000 10007 -10008)" universe="0" />
<cell id="10002" material="10002" name="clad" region="(-10002 10001 10007 -10008)" universe="0" />
<cell id="10003" material="10000" name="moderator" region="(10002 10003 -10004 10005 -10006 10007 -10008)" universe="0" />
<surface boundary="transmission" coeffs="0.0 0.0 0.39218" id="10000" name="fuel outer" type="z-cylinder" />
<surface boundary="transmission" coeffs="0.0 0.0 0.40005" id="10001" name="clad inner" type="z-cylinder" />
<surface boundary="transmission" coeffs="0.0 0.0 0.4572" id="10002" name="clad outer" type="z-cylinder" />
<surface boundary="reflective" coeffs="-0.62992" id="10003" type="x-plane" />
<surface boundary="reflective" coeffs="0.62992" id="10004" type="x-plane" />
<surface boundary="reflective" coeffs="-0.62992" id="10005" type="y-plane" />
<surface boundary="reflective" coeffs="0.62992" id="10006" type="y-plane" />
<surface boundary="reflective" coeffs="-1.0" id="10007" type="z-plane" />
<surface boundary="reflective" coeffs="1.0" id="10008" type="z-plane" />
</geometry>

View file

@ -0,0 +1,50 @@
<?xml version='1.0' encoding='utf-8'?>
<materials>
<default_xs>71c</default_xs>
<material id="10000" name="water">
<density units="sum" />
<nuclide ao="6.00992e-05" name="O-17" />
<nuclide ao="0.024672" name="O-16" />
<nuclide ao="8.0042e-06" name="B-10" />
<nuclide ao="3.2218e-05" name="B-11" />
<nuclide ao="0.049457" name="H-1" />
<nuclide ao="7.4196e-06" name="H-2" />
</material>
<material id="10001" name="fuel 2.4%">
<density units="sum" />
<nuclide ao="0.022407" name="U-238" />
<nuclide ao="0.000111633" name="O-17" />
<nuclide ao="8.0042e-06" name="B-10" />
<nuclide ao="0.045828" name="O-16" />
<nuclide ao="4.4842e-06" name="U-234" />
<nuclide ao="0.00055814" name="U-235" />
</material>
<material id="10002" name="zircaloy 4">
<density units="sum" />
<nuclide ao="0.00011669" name="Sn-118" />
<nuclide ao="4.1387e-05" name="Sn-119" />
<nuclide ao="4.6735e-06" name="Sn-112" />
<nuclide ao="3.1799e-06" name="Sn-114" />
<nuclide ao="1.6381e-06" name="Sn-115" />
<nuclide ao="7.0055e-05" name="Sn-116" />
<nuclide ao="3.7003e-05" name="Sn-117" />
<nuclide ao="0.0001361" name="Fe-56" />
<nuclide ao="3.1431e-06" name="Fe-57" />
<nuclide ao="8.6699e-06" name="Fe-54" />
<nuclide ao="4.1829e-07" name="Fe-58" />
<nuclide ao="0.0072758" name="Zr-92" />
<nuclide ao="2.7897e-05" name="Sn-124" />
<nuclide ao="0.00015697" name="Sn-120" />
<nuclide ao="2.2308e-05" name="Sn-122" />
<nuclide ao="1.7941e-06" name="Cr-54" />
<nuclide ao="6.3564e-05" name="Cr-52" />
<nuclide ao="7.2076e-06" name="Cr-53" />
<nuclide ao="3.2962e-06" name="Cr-50" />
<nuclide ao="7.4887e-07" name="O-17" />
<nuclide ao="0.00030743" name="O-16" />
<nuclide ao="0.00476" name="Zr-91" />
<nuclide ao="0.021827" name="Zr-90" />
<nuclide ao="0.0011879" name="Zr-96" />
<nuclide ao="0.0073734" name="Zr-94" />
</material>
</materials>

View file

@ -0,0 +1,5 @@
k-combined:
1.055198E+00 5.785829E-02
tally 1:
3.081784E+03
1.905482E+06

View file

@ -0,0 +1,18 @@
<?xml version='1.0' encoding='utf-8'?>
<settings>
<eigenvalue>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
</eigenvalue>
<source>
<space type="box">
<parameters>-0.6 -0.6 -0.9 0.6 0.6 0.9</parameters>
</space>
</source>
<entropy>
<dimension>10 10 10</dimension>
<lower_left>-0.7 -0.7 -1.1</lower_left>
<upper_right>0.7 0.7 1.1</upper_right>
</entropy>
</settings>

View file

@ -0,0 +1,10 @@
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<tally id="1">
<scores>keff</scores>
<!--
<derivative variable="density" material="10000"/>
-->
<derivative variable="nuclide_density" material="10001" nuclide="U-235.71c"/>
</tally>
</tallies>

View file

@ -0,0 +1,11 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*', True)
harness.main()