Merge pull request #331 from wbinventor/python-api

OpenMC Python API
This commit is contained in:
Paul Romano 2015-04-23 21:27:53 +07:00
commit f3cd40ec8c
151 changed files with 12179 additions and 1746 deletions

6
.gitignore vendored
View file

@ -11,6 +11,9 @@
# OpenMC executable
src/openmc
# Inputs generated from Python API
examples/python/**/*.xml
# emacs backups
*~
@ -41,3 +44,6 @@ results_error.dat
# Data downloaded from NNDC
data/nndc
.idea
.idea/*

View file

@ -784,7 +784,7 @@ The following quadratic surfaces can be modeled:
Each ``<cell>`` element can have the following attributes or sub-elements:
:id:
A unique integer that can be used to identify the surface.
A unique integer that can be used to identify the cell.
*Default*: None

View file

@ -0,0 +1,141 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 15
inactive = 5
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate some Nuclides
h1 = openmc.Nuclide('H-1')
o16 = openmc.Nuclide('O-16')
u235 = openmc.Nuclide('U-235')
# Instantiate some Materials and register the appropriate Nuclides
moderator = openmc.Material(material_id=41, name='moderator')
moderator.set_density('g/cc', 1.0)
moderator.add_nuclide(h1, 2.)
moderator.add_nuclide(o16, 1.)
moderator.add_s_alpha_beta('HH2O', '71t')
fuel = openmc.Material(material_id=40, name='fuel')
fuel.set_density('g/cc', 4.5)
fuel.add_nuclide(u235, 1.)
# Instantiate a MaterialsFile, register all Materials, and export to XML
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '71c'
materials_file.add_materials([moderator, fuel])
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml File
###############################################################################
# Instantiate ZCylinder surfaces
surf1 = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=7, name='surf 1')
surf2 = openmc.ZCylinder(surface_id=2, x0=0, y0=0, R=9, name='surf 2')
surf3 = openmc.ZCylinder(surface_id=3, x0=0, y0=0, R=11, name='surf 3')
surf3.boundary_type = 'vacuum'
# Instantiate Cells
cell1 = openmc.Cell(cell_id=1, name='cell 1')
cell2 = openmc.Cell(cell_id=100, name='cell 2')
cell3 = openmc.Cell(cell_id=101, name='cell 3')
cell4 = openmc.Cell(cell_id=2, name='cell 4')
# Register Surfaces with Cells
cell1.add_surface(surface=surf2, halfspace=-1)
cell2.add_surface(surface=surf1, halfspace=-1)
cell3.add_surface(surface=surf1, halfspace=+1)
cell4.add_surface(surface=surf2, halfspace=+1)
cell4.add_surface(surface=surf3, halfspace=-1)
# Register Materials with Cells
cell2.fill = fuel
cell3.fill = moderator
cell4.fill = moderator
# Instantiate Universes
universe1 = openmc.Universe(universe_id=37)
root = openmc.Universe(universe_id=0, name='root universe')
cell1.fill = universe1
# Register Cells with Universes
universe1.add_cells([cell2, cell3])
root.add_cells([cell1, cell4])
# Instantiate a Geometry and register the root Universe
geometry = openmc.Geometry()
geometry.root_universe = root
# Instantiate a GeometryFile, register Geometry, and export to XML
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
geometry_file.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a SettingsFile, set all runtime parameters, and export to XML
settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-4, -4, -4, 4, 4, 4])
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml File
###############################################################################
# Instantiate some tally Filters
cell_filter = openmc.Filter(type='cell', bins=100)
energy_filter = openmc.Filter(type='energy', bins=[0., 20.])
energyout_filter = openmc.Filter(type='energyout', bins=[0., 20.])
# Instantiate the first Tally
first_tally = openmc.Tally(tally_id=1, label='first tally')
first_tally.add_filter(cell_filter)
scores = ['total', 'scatter', 'nu-scatter', \
'absorption', 'fission', 'nu-fission']
for score in scores:
first_tally.add_score(score)
# Instantiate the second Tally
second_tally = openmc.Tally(tally_id=2, label='second tally')
second_tally.add_filter(cell_filter)
second_tally.add_filter(energy_filter)
scores = ['total', 'scatter', 'nu-scatter', \
'absorption', 'fission', 'nu-fission']
for score in scores:
second_tally.add_score(score)
# Instantiate the third Tally
third_tally = openmc.Tally(tally_id=3, label='third tally')
third_tally.add_filter(cell_filter)
third_tally.add_filter(energy_filter)
third_tally.add_filter(energyout_filter)
scores = ['scatter', 'nu-scatter', 'nu-fission']
for score in scores:
third_tally.add_score(score)
# Instantiate a TalliesFile, register all Tallies, and export to XML
tallies_file = openmc.TalliesFile()
tallies_file.add_tally(first_tally)
tallies_file.add_tally(second_tally)
tallies_file.add_tally(third_tally)
tallies_file.export_to_xml()

View file

@ -0,0 +1,158 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 20
inactive = 10
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate some Nuclides
h1 = openmc.Nuclide('H-1')
o16 = openmc.Nuclide('O-16')
u235 = openmc.Nuclide('U-235')
fe56 = openmc.Nuclide('Fe-56')
# Instantiate some Materials and register the appropriate Nuclides
fuel = openmc.Material(material_id=1, name='fuel')
fuel.set_density('g/cc', 4.5)
fuel.add_nuclide(u235, 1.)
moderator = openmc.Material(material_id=2, name='moderator')
moderator.set_density('g/cc', 1.0)
moderator.add_nuclide(h1, 2.)
moderator.add_nuclide(o16, 1.)
moderator.add_s_alpha_beta('HH2O', '71t')
iron = openmc.Material(material_id=3, name='iron')
iron.set_density('g/cc', 7.9)
iron.add_nuclide(fe56, 1.)
# Instantiate a MaterialsFile, register all Materials, and export to XML
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '71c'
materials_file.add_materials([moderator, fuel, iron])
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml File
###############################################################################
# Instantiate Surfaces
left = openmc.XPlane(surface_id=1, x0=-3, name='left')
right = openmc.XPlane(surface_id=2, x0=3, name='right')
bottom = openmc.YPlane(surface_id=3, y0=-4, name='bottom')
top = openmc.YPlane(surface_id=4, y0=4, name='top')
fuel_surf = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4)
left.boundary_type = 'vacuum'
right.boundary_type = 'vacuum'
top.boundary_type = 'vacuum'
bottom.boundary_type = 'vacuum'
# Instantiate Cells
cell1 = openmc.Cell(cell_id=1, name='Cell 1')
cell2 = openmc.Cell(cell_id=101, name='cell 2')
cell3 = openmc.Cell(cell_id=102, name='cell 3')
cell4 = openmc.Cell(cell_id=500, name='cell 4')
cell5 = openmc.Cell(cell_id=600, name='cell 5')
cell6 = openmc.Cell(cell_id=601, name='cell 6')
# Register Surfaces with Cells
cell1.add_surface(left, halfspace=+1)
cell1.add_surface(right, halfspace=-1)
cell1.add_surface(bottom, halfspace=+1)
cell1.add_surface(top, halfspace=-1)
cell2.add_surface(fuel_surf, halfspace=-1)
cell3.add_surface(fuel_surf, halfspace=+1)
cell5.add_surface(fuel_surf, halfspace=-1)
cell6.add_surface(fuel_surf, halfspace=+1)
# Register Materials with Cells
cell2.fill = fuel
cell3.fill = moderator
cell4.fill = moderator
cell5.fill = iron
cell6.fill = moderator
# Instantiate Universe
univ1 = openmc.Universe(universe_id=1)
univ2 = openmc.Universe(universe_id=3)
univ3 = openmc.Universe(universe_id=4)
root = openmc.Universe(universe_id=0, name='root universe')
# Register Cells with Universe
univ1.add_cells([cell2, cell3])
univ2.add_cells([cell4])
univ3.add_cells([cell5, cell6])
root.add_cell(cell1)
# Instantiate a Lattice
lattice = openmc.HexLattice(lattice_id=5)
lattice.center = [0., 0., 0.]
lattice.pitch = [1., 2.]
lattice.universes = \
[ [ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ],
[ [univ2] + [univ1]*11, [univ2] + [univ1]*5, [univ1] ],
[ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ] ]
lattice.outer = univ2
# Fill Cell with the Lattice
cell1.fill = lattice
# Instantiate a Geometry and register the root Universe
geometry = openmc.Geometry()
geometry.root_universe = root
# Instantiate a GeometryFile, register Geometry, and export to XML
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
geometry_file.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a SettingsFile, set all runtime parameters, and export to XML
settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC plots.xml File
###############################################################################
plot_xy = openmc.Plot(plot_id=1)
plot_xy.filename = 'plot_xy'
plot_xy.origin = [0, 0, 0]
plot_xy.width = [6, 6]
plot_xy.pixels = [400, 400]
plot_xy.color = 'mat'
plot_yz = openmc.Plot(plot_id=2)
plot_yz.filename = 'plot_yz'
plot_yz.basis = 'yz'
plot_yz.origin = [0, 0, 0]
plot_yz.width = [8, 8]
plot_yz.pixels = [400, 400]
plot_yz.color = 'mat'
# Instantiate a PlotsFile, add Plot, and export to XML
plot_file = openmc.PlotsFile()
plot_file.add_plot(plot_xy)
plot_file.add_plot(plot_yz)
plot_file.export_to_xml()

View file

@ -0,0 +1,189 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 20
inactive = 10
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate some Nuclides
h1 = openmc.Nuclide('H-1')
o16 = openmc.Nuclide('O-16')
u235 = openmc.Nuclide('U-235')
# Instantiate some Materials and register the appropriate Nuclides
fuel = openmc.Material(material_id=1, name='fuel')
fuel.set_density('g/cc', 4.5)
fuel.add_nuclide(u235, 1.)
moderator = openmc.Material(material_id=2, name='moderator')
moderator.set_density('g/cc', 1.0)
moderator.add_nuclide(h1, 2.)
moderator.add_nuclide(o16, 1.)
moderator.add_s_alpha_beta('HH2O', '71t')
# Instantiate a MaterialsFile, register all Materials, and export to XML
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '71c'
materials_file.add_materials([moderator, fuel])
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml File
###############################################################################
# Instantiate Surfaces
left = openmc.XPlane(surface_id=1, x0=-2, name='left')
right = openmc.XPlane(surface_id=2, x0=2, name='right')
bottom = openmc.YPlane(surface_id=3, y0=-2, name='bottom')
top = openmc.YPlane(surface_id=4, y0=2, name='top')
fuel1 = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4)
fuel2 = openmc.ZCylinder(surface_id=6, x0=0, y0=0, R=0.3)
fuel3 = openmc.ZCylinder(surface_id=7, x0=0, y0=0, R=0.2)
left.boundary_type = 'vacuum'
right.boundary_type = 'vacuum'
top.boundary_type = 'vacuum'
bottom.boundary_type = 'vacuum'
# Instantiate Cells
cell1 = openmc.Cell(cell_id=1, name='Cell 1')
cell2 = openmc.Cell(cell_id=2, name='Cell 2')
cell3 = openmc.Cell(cell_id=101, name='cell 3')
cell4 = openmc.Cell(cell_id=102, name='cell 4')
cell5 = openmc.Cell(cell_id=201, name='cell 5')
cell6 = openmc.Cell(cell_id=202, name='cell 6')
cell7 = openmc.Cell(cell_id=301, name='cell 7')
cell8 = openmc.Cell(cell_id=302, name='cell 8')
# Register Surfaces with Cells
cell1.add_surface(left, halfspace=+1)
cell1.add_surface(right, halfspace=-1)
cell1.add_surface(bottom, halfspace=+1)
cell1.add_surface(top, halfspace=-1)
cell2.add_surface(left, halfspace=+1)
cell2.add_surface(right, halfspace=-1)
cell2.add_surface(bottom, halfspace=+1)
cell2.add_surface(top, halfspace=-1)
cell3.add_surface(fuel1, halfspace=-1)
cell4.add_surface(fuel1, halfspace=+1)
cell5.add_surface(fuel2, halfspace=-1)
cell6.add_surface(fuel2, halfspace=+1)
cell7.add_surface(fuel3, halfspace=-1)
cell8.add_surface(fuel3, halfspace=+1)
# Register Materials with Cells
cell3.fill = fuel
cell4.fill = moderator
cell5.fill = fuel
cell6.fill = moderator
cell7.fill = fuel
cell8.fill = moderator
# Instantiate Universe
univ1 = openmc.Universe(universe_id=1)
univ2 = openmc.Universe(universe_id=2)
univ3 = openmc.Universe(universe_id=3)
univ4 = openmc.Universe(universe_id=5)
root = openmc.Universe(universe_id=0, name='root universe')
# Register Cells with Universe
univ1.add_cells([cell3, cell4])
univ2.add_cells([cell5, cell6])
univ3.add_cells([cell7, cell8])
root.add_cell(cell1)
univ4.add_cell(cell2)
# Instantiate nested Lattices
lattice1 = openmc.RectLattice(lattice_id=4, name='4x4 assembly')
lattice1.dimension = [2, 2]
lattice1.lower_left = [-1., -1.]
lattice1.pitch = [1., 1.]
lattice1.universes = [[univ1, univ2],
[univ2, univ3]]
lattice2 = openmc.RectLattice(lattice_id=6, name='4x4 core')
lattice2.dimension = [2, 2]
lattice2.lower_left = [-2., -2.]
lattice2.pitch = [2., 2.]
lattice2.universes = [[univ4, univ4],
[univ4, univ4]]
# Fill Cell with the Lattice
cell1.fill = lattice2
cell2.fill = lattice1
# Instantiate a Geometry and register the root Universe
geometry = openmc.Geometry()
geometry.root_universe = root
# Instantiate a GeometryFile, register Geometry, and export to XML
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
geometry_file.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a SettingsFile, set all runtime parameters, and export to XML
settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC plots.xml File
###############################################################################
plot = openmc.Plot(plot_id=1)
plot.origin = [0, 0, 0]
plot.width = [4, 4]
plot.pixels = [400, 400]
plot.color = 'mat'
# Instantiate a PlotsFile, add Plot, and export to XML
plot_file = openmc.PlotsFile()
plot_file.add_plot(plot)
plot_file.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml File
###############################################################################
# Instantiate a tally mesh
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'rectangular'
mesh.dimension = [4, 4]
mesh.lower_left = [-2, -2]
mesh.width = [1, 1]
# Instantiate tally Filter
mesh_filter = openmc.Filter()
mesh_filter.mesh = mesh
# Instantiate the Tally
tally = openmc.Tally(tally_id=1)
tally.add_filter(mesh_filter)
tally.add_score('total')
# Instantiate a TalliesFile, register Tally/Mesh, and export to XML
tallies_file = openmc.TalliesFile()
tallies_file.add_mesh(mesh)
tallies_file.add_tally(tally)
tallies_file.export_to_xml()

View file

@ -0,0 +1,176 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 20
inactive = 10
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate some Nuclides
h1 = openmc.Nuclide('H-1')
o16 = openmc.Nuclide('O-16')
u235 = openmc.Nuclide('U-235')
# Instantiate some Materials and register the appropriate Nuclides
fuel = openmc.Material(material_id=1, name='fuel')
fuel.set_density('g/cc', 4.5)
fuel.add_nuclide(u235, 1.)
moderator = openmc.Material(material_id=2, name='moderator')
moderator.set_density('g/cc', 1.0)
moderator.add_nuclide(h1, 2.)
moderator.add_nuclide(o16, 1.)
moderator.add_s_alpha_beta('HH2O', '71t')
# Instantiate a MaterialsFile, register all Materials, and export to XML
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '71c'
materials_file.add_materials([moderator, fuel])
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml File
###############################################################################
# Instantiate Surfaces
left = openmc.XPlane(surface_id=1, x0=-2, name='left')
right = openmc.XPlane(surface_id=2, x0=2, name='right')
bottom = openmc.YPlane(surface_id=3, y0=-2, name='bottom')
top = openmc.YPlane(surface_id=4, y0=2, name='top')
fuel1 = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4)
fuel2 = openmc.ZCylinder(surface_id=6, x0=0, y0=0, R=0.3)
fuel3 = openmc.ZCylinder(surface_id=7, x0=0, y0=0, R=0.2)
left.boundary_type = 'vacuum'
right.boundary_type = 'vacuum'
top.boundary_type = 'vacuum'
bottom.boundary_type = 'vacuum'
# Instantiate Cells
cell1 = openmc.Cell(cell_id=1, name='Cell 1')
cell2 = openmc.Cell(cell_id=101, name='cell 2')
cell3 = openmc.Cell(cell_id=102, name='cell 3')
cell4 = openmc.Cell(cell_id=201, name='cell 4')
cell5 = openmc.Cell(cell_id=202, name='cell 5')
cell6 = openmc.Cell(cell_id=301, name='cell 6')
cell7 = openmc.Cell(cell_id=302, name='cell 7')
# Register Surfaces with Cells
cell1.add_surface(left, halfspace=+1)
cell1.add_surface(right, halfspace=-1)
cell1.add_surface(bottom, halfspace=+1)
cell1.add_surface(top, halfspace=-1)
cell2.add_surface(fuel1, halfspace=-1)
cell3.add_surface(fuel1, halfspace=+1)
cell4.add_surface(fuel2, halfspace=-1)
cell5.add_surface(fuel2, halfspace=+1)
cell6.add_surface(fuel3, halfspace=-1)
cell7.add_surface(fuel3, halfspace=+1)
# Register Materials with Cells
cell2.fill = fuel
cell3.fill = moderator
cell4.fill = fuel
cell5.fill = moderator
cell6.fill = fuel
cell7.fill = moderator
# Instantiate Universe
univ1 = openmc.Universe(universe_id=1)
univ2 = openmc.Universe(universe_id=2)
univ3 = openmc.Universe(universe_id=3)
root = openmc.Universe(universe_id=0, name='root universe')
# Register Cells with Universe
univ1.add_cells([cell2, cell3])
univ2.add_cells([cell4, cell5])
univ3.add_cells([cell6, cell7])
root.add_cell(cell1)
# Instantiate a Lattice
lattice = openmc.RectLattice(lattice_id=5)
lattice.dimension = [4, 4]
lattice.lower_left = [-2., -2.]
lattice.pitch = [1., 1.]
lattice.universes = [[univ1, univ2, univ1, univ2],
[univ2, univ3, univ2, univ3],
[univ1, univ2, univ1, univ2],
[univ2, univ3, univ2, univ3]]
# Fill Cell with the Lattice
cell1.fill = lattice
# Instantiate a Geometry and register the root Universe
geometry = openmc.Geometry()
geometry.root_universe = root
# Instantiate a GeometryFile, register Geometry, and export to XML
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
geometry_file.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a SettingsFile, set all runtime parameters, and export to XML
settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC plots.xml File
###############################################################################
plot = openmc.Plot(plot_id=1)
plot.origin = [0, 0, 0]
plot.width = [4, 4]
plot.pixels = [400, 400]
plot.color = 'mat'
# Instantiate a PlotsFile, add Plot, and export to XML
plot_file = openmc.PlotsFile()
plot_file.add_plot(plot)
plot_file.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml File
###############################################################################
# Instantiate a tally mesh
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'rectangular'
mesh.dimension = [4, 4]
mesh.lower_left = [-2, -2]
mesh.width = [1, 1]
# Instantiate tally Filter
mesh_filter = openmc.Filter()
mesh_filter.mesh = mesh
# Instantiate the Tally
tally = openmc.Tally(tally_id=1)
tally.add_filter(mesh_filter)
tally.add_score('total')
# Instantiate a TalliesFile, register Tally/Mesh, and export to XML
tallies_file = openmc.TalliesFile()
tallies_file.add_mesh(mesh)
tallies_file.add_tally(tally)
tallies_file.export_to_xml()

View file

@ -0,0 +1,214 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 100
inactive = 10
particles = 1000
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate some Nuclides
h1 = openmc.Nuclide('H-1')
h2 = openmc.Nuclide('H-2')
he4 = openmc.Nuclide('He-4')
b10 = openmc.Nuclide('B-10')
b11 = openmc.Nuclide('B-11')
o16 = openmc.Nuclide('O-16')
o17 = openmc.Nuclide('O-17')
cr50 = openmc.Nuclide('Cr-50')
cr52 = openmc.Nuclide('Cr-52')
cr53 = openmc.Nuclide('Cr-53')
cr54 = openmc.Nuclide('Cr-54')
fe54 = openmc.Nuclide('Fe-54')
fe56 = openmc.Nuclide('Fe-56')
fe57 = openmc.Nuclide('Fe-57')
fe58 = openmc.Nuclide('Fe-58')
zr90 = openmc.Nuclide('Zr-90')
zr91 = openmc.Nuclide('Zr-91')
zr92 = openmc.Nuclide('Zr-92')
zr94 = openmc.Nuclide('Zr-94')
zr96 = openmc.Nuclide('Zr-96')
sn112 = openmc.Nuclide('Sn-112')
sn114 = openmc.Nuclide('Sn-114')
sn115 = openmc.Nuclide('Sn-115')
sn116 = openmc.Nuclide('Sn-116')
sn117 = openmc.Nuclide('Sn-117')
sn118 = openmc.Nuclide('Sn-118')
sn119 = openmc.Nuclide('Sn-119')
sn120 = openmc.Nuclide('Sn-120')
sn122 = openmc.Nuclide('Sn-122')
sn124 = openmc.Nuclide('Sn-124')
u234 = openmc.Nuclide('U-234')
u235 = openmc.Nuclide('U-235')
u238 = openmc.Nuclide('U-238')
# Instantiate some Materials and register the appropriate Nuclides
uo2 = openmc.Material(material_id=1, name='UO2 fuel at 2.4% wt enrichment')
uo2.set_density('g/cm3', 10.29769)
uo2.add_nuclide(u234, 4.4843e-6)
uo2.add_nuclide(u235, 5.5815e-4)
uo2.add_nuclide(u238, 2.2408e-2)
uo2.add_nuclide(o16, 4.5829e-2)
uo2.add_nuclide(o17, 1.1164e-4)
helium = openmc.Material(material_id=2, name='Helium for gap')
helium.set_density('g/cm3', 0.001598)
helium.add_nuclide(he4, 2.4044e-4)
zircaloy = openmc.Material(material_id=3, name='Zircaloy 4')
zircaloy.set_density('g/cm3', 6.55)
zircaloy.add_nuclide(o16, 3.0743e-4)
zircaloy.add_nuclide(o17, 7.4887e-7)
zircaloy.add_nuclide(cr50, 3.2962e-6)
zircaloy.add_nuclide(cr52, 6.3564e-5)
zircaloy.add_nuclide(cr53, 7.2076e-6)
zircaloy.add_nuclide(cr54, 1.7941e-6)
zircaloy.add_nuclide(fe54, 8.6699e-6)
zircaloy.add_nuclide(fe56, 1.3610e-4)
zircaloy.add_nuclide(fe57, 3.1431e-6)
zircaloy.add_nuclide(fe58, 4.1829e-7)
zircaloy.add_nuclide(zr90, 2.1827e-2)
zircaloy.add_nuclide(zr91, 4.7600e-3)
zircaloy.add_nuclide(zr92, 7.2758e-3)
zircaloy.add_nuclide(zr94, 7.3734e-3)
zircaloy.add_nuclide(zr96, 1.1879e-3)
zircaloy.add_nuclide(sn112, 4.6735e-6)
zircaloy.add_nuclide(sn114, 3.1799e-6)
zircaloy.add_nuclide(sn115, 1.6381e-6)
zircaloy.add_nuclide(sn116, 7.0055e-5)
zircaloy.add_nuclide(sn117, 3.7003e-5)
zircaloy.add_nuclide(sn118, 1.1669e-4)
zircaloy.add_nuclide(sn119, 4.1387e-5)
zircaloy.add_nuclide(sn120, 1.5697e-4)
zircaloy.add_nuclide(sn122, 2.2308e-5)
zircaloy.add_nuclide(sn124, 2.7897e-5)
borated_water = openmc.Material(material_id=4, name='Borated water at 975 ppm')
borated_water.set_density('g/cm3', 0.740582)
borated_water.add_nuclide(b10, 8.0042e-6)
borated_water.add_nuclide(b11, 3.2218e-5)
borated_water.add_nuclide(h1, 4.9457e-2)
borated_water.add_nuclide(h2, 7.4196e-6)
borated_water.add_nuclide(o16, 2.4672e-2)
borated_water.add_nuclide(o17, 6.0099e-5)
borated_water.add_s_alpha_beta('HH2O', '71t')
# Instantiate a MaterialsFile, register all Materials, and export to XML
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '71c'
materials_file.add_materials([uo2, helium, zircaloy, borated_water])
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml File
###############################################################################
# Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, R=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(surface_id=3, x0=0, y0=0, R=0.45720, name='Clad OR')
left = openmc.XPlane(surface_id=4, x0=-0.62992, name='left')
right = openmc.XPlane(surface_id=5, x0=0.62992, name='right')
bottom = openmc.YPlane(surface_id=6, y0=-0.62992, name='bottom')
top = openmc.YPlane(surface_id=7, y0=0.62992, name='top')
left.boundary_type = 'reflective'
right.boundary_type = 'reflective'
top.boundary_type = 'reflective'
bottom.boundary_type = 'reflective'
# Instantiate Cells
fuel = openmc.Cell(cell_id=1, name='cell 1')
gap = openmc.Cell(cell_id=2, name='cell 2')
clad = openmc.Cell(cell_id=3, name='cell 3')
water = openmc.Cell(cell_id=4, name='cell 4')
# Register Surfaces with Cells
fuel.add_surface(fuel_or, halfspace=-1)
gap.add_surface(fuel_or, halfspace=+1)
gap.add_surface(clad_ir, halfspace=-1)
clad.add_surface(clad_ir, halfspace=+1)
clad.add_surface(clad_or, halfspace=-1)
water.add_surface(clad_or, halfspace=+1)
water.add_surface(left, halfspace=+1)
water.add_surface(right, halfspace=-1)
water.add_surface(bottom, halfspace=+1)
water.add_surface(top, halfspace=-1)
# Register Materials with Cells
fuel.fill = uo2
gap.fill = helium
clad.fill = zircaloy
water.fill = borated_water
# Instantiate Universe
root = openmc.Universe(universe_id=0, name='root universe')
# Register Cells with Universe
root.add_cells([fuel, gap, clad, water])
# Instantiate a Geometry and register the root Universe
geometry = openmc.Geometry()
geometry.root_universe = root
# Instantiate a GeometryFile, register Geometry, and export to XML
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
geometry_file.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a SettingsFile, set all runtime parameters, and export to XML
settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-0.62992, -0.62992, -1, \
0.62992, 0.62992, 1])
settings_file.entropy_lower_left = [-0.39218, -0.39218, -1.e50]
settings_file.entropy_upper_right = [0.39218, 0.39218, 1.e50]
settings_file.entropy_dimension = [10, 10, 1]
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml File
###############################################################################
# Instantiate a tally mesh
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'rectangular'
mesh.dimension = [100, 100, 1]
mesh.lower_left = [-0.62992, -0.62992, -1.e50]
mesh.upper_right = [0.62992, 0.62992, 1.e50]
# Instantiate some tally Filters
energy_filter = openmc.Filter(type='energy', bins=[0., 4.e-6, 20.])
mesh_filter = openmc.Filter()
mesh_filter.mesh = mesh
# Instantiate the Tally
tally = openmc.Tally(tally_id=1)
tally.add_filter(energy_filter)
tally.add_filter(mesh_filter)
tally.add_score('flux')
tally.add_score('fission')
tally.add_score('nu-fission')
# Instantiate a TalliesFile, register all Tallies, and export to XML
tallies_file = openmc.TalliesFile()
tallies_file.add_mesh(mesh)
tallies_file.add_tally(tally)
tallies_file.export_to_xml()

View file

@ -0,0 +1,92 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 500
inactive = 10
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate a Nuclides
u235 = openmc.Nuclide('U-235')
# Instantiate a Material and register the Nuclide
fuel = openmc.Material(material_id=1, name='fuel')
fuel.set_density('g/cc', 4.5)
fuel.add_nuclide(u235, 1.)
# Instantiate a MaterialsFile, register Material, and export to XML
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '71c'
materials_file.add_material(fuel)
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml File
###############################################################################
# Instantiate Surfaces
surf1 = openmc.XPlane(surface_id=1, x0=-1, name='surf 1')
surf2 = openmc.XPlane(surface_id=2, x0=+1, name='surf 2')
surf3 = openmc.YPlane(surface_id=3, y0=-1, name='surf 3')
surf4 = openmc.YPlane(surface_id=4, y0=+1, name='surf 4')
surf5 = openmc.ZPlane(surface_id=5, z0=-1, name='surf 5')
surf6 = openmc.ZPlane(surface_id=6, z0=+1, name='surf 6')
surf1.boundary_type = 'vacuum'
surf2.boundary_type = 'vacuum'
surf3.boundary_type = 'reflective'
surf4.boundary_type = 'reflective'
surf5.boundary_type = 'reflective'
surf6.boundary_type = 'reflective'
# Instantiate Cell
cell = openmc.Cell(cell_id=1, name='cell 1')
# Register Surfaces with Cell
cell.add_surface(surface=surf1, halfspace=+1)
cell.add_surface(surface=surf2, halfspace=-1)
cell.add_surface(surface=surf3, halfspace=+1)
cell.add_surface(surface=surf4, halfspace=-1)
cell.add_surface(surface=surf5, halfspace=+1)
cell.add_surface(surface=surf6, halfspace=-1)
# Register Material with Cell
cell.fill = fuel
# Instantiate Universes
root = openmc.Universe(universe_id=0, name='root universe')
# Register Cell with Universe
root.add_cell(cell)
# Instantiate a Geometry and register the root Universe
geometry = openmc.Geometry()
geometry.root_universe = root
# Instantiate a GeometryFile, register Geometry, and export to XML
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
geometry_file.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a SettingsFile, set all runtime parameters, and export to XML
settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.export_to_xml()

File diff suppressed because one or more lines are too long

View file

@ -190,7 +190,6 @@ contains
! Deallocate temporary arrays for names of nuclides and S(a,b) tables
if (allocated(mat % names)) deallocate(mat % names)
if (allocated(mat % sab_names)) deallocate(mat % sab_names)
end do MATERIAL_LOOP2

View file

@ -207,9 +207,9 @@ contains
end if
! Apply rotation
if (allocated(c % rotation)) then
p % coord % xyz = matmul(c % rotation, p % coord % xyz)
p % coord % uvw = matmul(c % rotation, p % coord % uvw)
if (allocated(c % rotation_matrix)) then
p % coord % xyz = matmul(c % rotation_matrix, p % coord % xyz)
p % coord % uvw = matmul(c % rotation_matrix, p % coord % uvw)
p % coord % rotated = .true.
end if

View file

@ -141,8 +141,9 @@ module geometry_header
! here too
! Rotation matrix and translation vector
real(8), allocatable :: rotation(:,:)
real(8), allocatable :: translation(:)
real(8), allocatable :: rotation(:)
real(8), allocatable :: rotation_matrix(:,:)
end type Cell
! array index of universe 0

View file

@ -206,7 +206,7 @@ contains
logical :: status ! does the group exist
! Check if group exists
call h5ltpath_valid_f(hdf5_fh, trim(group), .true., status, hdf5_err)
call h5ltpath_valid_f(hdf5_fh, trim(group), .true., status, hdf5_err)
! Either create or open group
if (status) then
@ -257,7 +257,7 @@ contains
integer(HID_T), intent(in) :: group ! name of group
character(*), intent(in) :: name ! name of data
integer, intent(inout) :: buffer ! read data to here
integer, intent(inout) :: buffer ! read data to here
integer :: buffer_copy(1) ! need an array for read
@ -461,7 +461,7 @@ contains
integer(HID_T), intent(in) :: group ! name of group
character(*), intent(in) :: name ! name of data
real(8), intent(inout) :: buffer ! read data to here
real(8), intent(inout) :: buffer ! read data to here
real(8) :: buffer_copy(1) ! need an array for read
@ -726,7 +726,7 @@ contains
call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err)
! Set up dimesnions of string to write
dims2 = (/length, 1/) ! full array of strings to write
dims2 = (/length, 1/) ! full array of strings to write
dims1(1) = length ! length of string
! Copy over string buffer to a rank 1 array
@ -760,7 +760,7 @@ contains
! character(len=length, kind=c_char), pointer :: chr_ptr
! f_ptr = c_loc(buf_ptr(1))
! call h5dread_f(dset, H5T_STRING, f_ptr, hdf5_err, xfer_prp=plist)
! call c_f_pointer(buf_ptr(1), chr_ptr)
! call c_f_pointer(buf_ptr(1), chr_ptr)
! buffer = chr_ptr
! nullify(chr_ptr)
@ -838,7 +838,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -916,7 +916,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -996,7 +996,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1077,7 +1077,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1159,7 +1159,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1238,7 +1238,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1316,7 +1316,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1396,7 +1396,7 @@ contains
f_ptr = c_loc(buffer(1,1))
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1477,7 +1477,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1559,7 +1559,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1639,7 +1639,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, long_type, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1726,7 +1726,7 @@ contains
call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err)
! Set up dimesnions of string to write
dims2 = (/length, 1/) ! full array of strings to write
dims2 = (/length, 1/) ! full array of strings to write
dims1(1) = length ! length of string
! Copy over string buffer to a rank 1 array
@ -1762,7 +1762,7 @@ contains
! character(len=length, kind=c_char), pointer :: chr_ptr
! f_ptr = c_loc(buf_ptr(1))
! call h5dread_f(dset, H5T_STRING, f_ptr, hdf5_err, xfer_prp=plist)
! call c_f_pointer(buf_ptr(1), chr_ptr)
! call c_f_pointer(buf_ptr(1), chr_ptr)
! buffer = chr_ptr
! nullify(chr_ptr)

View file

@ -123,6 +123,10 @@ contains
CELL_LOOP: do i = 1, n_cells
c => cells(i)
! Write internal OpenMC index for this cell
call su % write_data(i, "index", &
group="geometry/cells/cell " // trim(to_str(c % id)))
! Write universe for this cell
call su % write_data(universes(c % universe) % id, "universe", &
group="geometry/cells/cell " // trim(to_str(c % id)))
@ -139,11 +143,33 @@ contains
call su % write_data(materials(c % material) % id, "material", &
group="geometry/cells/cell " // trim(to_str(c % id)))
end if
case (CELL_FILL)
call su % write_data("universe", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call su % write_data(universes(c % fill) % id, "material", &
call su % write_data(universes(c % fill) % id, "fill", &
group="geometry/cells/cell " // trim(to_str(c % id)))
if (allocated(c % translation)) then
call su % write_data(1, "translated", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call su % write_data(c % translation, "translation", length=3, &
group="geometry/cells/cell " // trim(to_str(c % id)))
else
call su % write_data(0, "translated", &
group="geometry/cells/cell " // trim(to_str(c % id)))
end if
if (allocated(c % rotation_matrix)) then
call su % write_data(1, "rotated", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call su % write_data(c % rotation, "rotation", length=3, &
group="geometry/cells/cell " // trim(to_str(c % id)))
else
call su % write_data(0, "rotated", &
group="geometry/cells/cell " // trim(to_str(c % id)))
end if
case (CELL_LATTICE)
call su % write_data("lattice", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
@ -170,6 +196,10 @@ contains
SURFACE_LOOP: do i = 1, n_surfaces
s => surfaces(i)
! Write internal OpenMC index for this surface
call su % write_data(i, "index", &
group="geometry/surfaces/surface " // trim(to_str(s % id)))
! Write surface type
select case (s % type)
case (SURF_PX)
@ -254,6 +284,10 @@ contains
UNIVERSE_LOOP: do i = 1, n_universes
u => universes(i)
! Write internal OpenMC index for this universe
call su % write_data(i, "index", &
group="geometry/universes/universe " // trim(to_str(u % id)))
! Write list of cells in this universe
if (u % n_cells > 0) then
call su % write_data(u % cells, "cells", length=u % n_cells, &
@ -273,17 +307,21 @@ contains
LATTICE_LOOP: do i = 1, n_lattices
lat => lattices(i) % obj
! Write internal OpenMC index for this lattice
call su % write_data(i, "index", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
! Write lattice type
select type (lat)
type is (RectLattice)
! Write lattice type.
call su % write_data("rectangular", "type", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
! Write number of lattice cells.
call su % write_data(lat % n_cells, "n_cells", length=3, &
! Write lattice dimensions, lower left corner, and pitch
call su % write_data(lat % n_cells, "dimension", length=3, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
! Write lattice lower-left.
if (lat % is_3d) then
call su % write_data(lat % lower_left, "lower_left", length=3, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
@ -292,7 +330,6 @@ contains
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
! Write lattice pitch.
if (lat % is_3d) then
call su % write_data(lat % pitch, "pitch", length=3, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
@ -301,6 +338,9 @@ contains
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
call su % write_data(lat % outer, "outer", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
! Write lattice universes.
allocate(lattice_universes(lat % n_cells(1), lat % n_cells(2), &
&lat % n_cells(3)))
@ -327,23 +367,25 @@ contains
call su % write_data(lat % n_rings, "n_axial", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
! Write lattice center.
! Write lattice center, pitch and outer universe.
if (lat % is_3d) then
call su % write_data(lat % center, "center", length=3, &
call su % write_data(lat % center, "center", length=3, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
else
call su % write_data(lat % center, "center", length=2, &
call su % write_data(lat % center, "center", length=2, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
if (lat % is_3d) then
call su % write_data(lat % pitch, "pitch", length=3, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
else
call su % write_data(lat % pitch, "pitch", length=2, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
! Write lattice pitch.
if (lat % is_3d) then
call su % write_data(lat % pitch, "pitch", length=2, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
else
call su % write_data(lat % pitch, "pitch", length=1, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
call su % write_data(lat % outer, "outer", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
! Write lattice universes.
allocate(lattice_universes(2*lat % n_rings - 1, 2*lat % n_rings - 1, &
@ -391,6 +433,10 @@ contains
do i = 1, n_materials
m => materials(i)
! Write internal OpenMC index for this material
call su % write_data(i, "index", &
group="materials/material " // trim(to_str(m % id)))
! Write atom density with units
call su % write_data(m % density, "atom_density", &
group="materials/material " // trim(to_str(m % id)))
@ -416,6 +462,9 @@ contains
group="materials/material " // trim(to_str(m % id)))
! Write S(a,b) information if present
call su % write_data(m % n_sab, "n_sab", &
group="materials/material " // trim(to_str(m % id)))
if (m % n_sab > 0) then
call su % write_data(m % i_sab_nuclides, "i_sab_nuclides", &
length=m % n_sab, &
@ -423,6 +472,11 @@ contains
call su % write_data(m % i_sab_tables, "i_sab_tables", &
length=m % n_sab, &
group="materials/material " // trim(to_str(m % id)))
do j = 1, m % n_sab
call su % write_data(m % sab_names(j), to_str(j), &
group="materials/material " // &
trim(to_str(m % id)) // "/sab_tables")
end do
end if
end do
@ -601,6 +655,10 @@ contains
NUCLIDE_LOOP: do i = 1, n_nuclides_total
nuc => nuclides(i)
! Write internal OpenMC index for this nuclide
call su % write_data(i, "index", &
group="nuclides/" // trim(nuc % name))
! Determine size of cross-sections
size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8
size_total = size_xs
@ -608,6 +666,8 @@ contains
! Write some basic attributes
call su % write_data(nuc % zaid, "zaid", &
group="nuclides/" // trim(nuc % name))
call su % write_data(xs_listings(nuc % listing) % alias, "alias", &
group="nuclides/" // trim(nuc % name))
call su % write_data(nuc % awr, "awr", &
group="nuclides/" // trim(nuc % name))
call su % write_data(nuc % kT, "kT", &

View file

@ -1057,13 +1057,14 @@ contains
! Copy rotation angles in x,y,z directions
call get_node_array(node_cell, "rotation", temp_double_array3)
c % rotation = temp_double_array3
phi = -temp_double_array3(1) * PI/180.0_8
theta = -temp_double_array3(2) * PI/180.0_8
psi = -temp_double_array3(3) * PI/180.0_8
! Calculate rotation matrix based on angles given
allocate(c % rotation(3,3))
c % rotation = reshape((/ &
allocate(c % rotation_matrix(3,3))
c % rotation_matrix = reshape((/ &
cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), &
-cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), &
cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), &

View file

@ -21,6 +21,7 @@ module mpiio_interface
module procedure mpi_write_integer_4Darray
module procedure mpi_write_long
module procedure mpi_write_string
!module procedure mpi_write_string_1Darray
end interface mpi_write_data
! Generic HDF5 read procedure interface
@ -37,6 +38,7 @@ module mpiio_interface
module procedure mpi_read_integer_4Darray
module procedure mpi_read_long
module procedure mpi_read_string
!module procedure mpi_read_string_1Darray
end interface mpi_read_data
contains

View file

@ -1761,6 +1761,91 @@ contains
#endif
end subroutine read_string
!===============================================================================
! WRITE_STRING_1DARRAY writes 1-D string data
!===============================================================================
subroutine write_string_1Darray(self, buffer, name, group, length, collect)
integer, intent(in) :: length ! length of array to write
character(*), intent(in) :: buffer(:) ! data to write
character(*), intent(in) :: name ! name of data
character(*), intent(in), optional :: group ! HDF5 group name
logical, intent(in), optional :: collect ! collective I/O
class(BinaryOutput) :: self
character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name
character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name
integer :: n
logical :: collect_
! Get string length
n = len_trim(buffer(1))
! Set name
name_ = trim(name)
! Set group
if (present(group)) then
group_ = trim(group)
end if
! Set up collective vs. independent I/O
if (present(collect)) then
collect_ = collect
else
collect_ = .true.
end if
#ifdef HDF5
! Check if HDF5 group should be created/opened
if (present(group)) then
call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp)
else
self % hdf5_grp = self % hdf5_fh
endif
if (present(group)) call hdf5_close_group(self % hdf5_grp)
#endif
end subroutine write_string_1Darray
!===============================================================================
! READ_STRING_1DARRAY reads 1-D string data
!===============================================================================
subroutine read_string_1Darray(self, buffer, name, group, length, collect)
integer, intent(in) :: length ! length of array to write
character(*), intent(inout) :: buffer(:) ! data to write
character(*), intent(in) :: name ! name of data
character(*), intent(in), optional :: group ! HDF5 group name
logical, intent(in), optional :: collect ! collective I/O
class(BinaryOutput) :: self
character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name
character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name
integer :: n
logical :: collect_
! Get string length
n = len(buffer(1))
! Set name
name_ = trim(name)
! Set group
if (present(group)) then
group_ = trim(group)
end if
! Set up collective vs. independent I/O
if (present(collect)) then
collect_ = collect
else
collect_ = .true.
end if
end subroutine read_string_1Darray
!===============================================================================
! WRITE_ATTRIBUTE_STRING

View file

@ -19,6 +19,8 @@ module state_point
use string, only: to_str, zero_padded, count_digits
use output_interface
use tally_header, only: TallyObject
use mesh_header, only: StructuredMesh
use dict_header, only: ElemKeyValueII, ElemKeyValueCI
#ifdef MPI
use mpi
@ -36,15 +38,17 @@ contains
subroutine write_state_point()
character(MAX_FILE_LEN) :: filename
integer :: i
integer :: j
integer :: k
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
character(8) :: moment_name ! name of moment (e.g, P3, Y-1,1)
integer, allocatable :: temp_array(:)
type(TallyObject), pointer :: t => null()
character(MAX_FILE_LEN) :: filename
integer :: i, j, k
integer, allocatable :: id_array(:)
integer, allocatable :: key_array(:)
type(StructuredMesh), pointer :: mesh
type(TallyObject), pointer :: tally
type(ElemKeyValueII), pointer :: current
type(ElemKeyValueII), pointer :: next
character(8) :: moment_name ! name of moment (e.g, P3)
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
! Set filename for state point
filename = trim(path_output) // 'statepoint.' // &
@ -87,17 +91,26 @@ contains
! Write run information
call sp % write_data(run_mode, "run_mode")
call sp % write_data(n_particles, "n_particles")
call sp % write_data(n_batches, "n_batches")
! Write out current batch number
call sp % write_data(current_batch, "current_batch")
! Indicate whether source bank is stored in statepoint
if (source_separate) then
call sp % write_data(0, "source_present")
else
call sp % write_data(1, "source_present")
end if
! Write out information for eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call sp % write_data(n_inactive, "n_inactive")
call sp % write_data(gen_per_batch, "gen_per_batch")
call sp % write_data(k_generation, "k_generation", &
length=current_batch*gen_per_batch)
call sp % write_data(entropy, "entropy", length=current_batch*gen_per_batch)
call sp % write_data(entropy, "entropy", &
length=current_batch*gen_per_batch)
call sp % write_data(k_col_abs, "k_col_abs")
call sp % write_data(k_col_tra, "k_col_tra")
call sp % write_data(k_abs_tra, "k_abs_tra")
@ -105,6 +118,10 @@ contains
! Write out CMFD info
if (cmfd_on) then
#ifdef HDF5
call sp % open_group("cmfd")
call sp % close_group()
#endif
call sp % write_data(1, "cmfd_on")
call sp % write_data(cmfd % indices, "indices", length=4, group="cmfd")
call sp % write_data(cmfd % k_cmfd, "k_cmfd", length=current_batch, &
@ -114,7 +131,7 @@ contains
cmfd % indices(2), cmfd % indices(3)/), &
group="cmfd")
call sp % write_data(cmfd % entropy, "cmfd_entropy", &
length=current_batch, group="cmfd")
length=current_batch, group="cmfd")
call sp % write_data(cmfd % balance, "cmfd_balance", &
length=current_batch, group="cmfd")
call sp % write_data(cmfd % dom, "cmfd_dominance", &
@ -126,151 +143,202 @@ contains
end if
end if
! Write number of meshes
call sp % write_data(n_meshes, "n_meshes", group="tallies")
#ifdef HDF5
call sp % open_group("tallies")
call sp % close_group()
#endif
! Write information for meshes
MESH_LOOP: do i = 1, n_meshes
call sp % write_data(meshes(i) % id, "id", &
group="tallies/mesh" // to_str(i))
call sp % write_data(meshes(i) % type, "type", &
group="tallies/mesh" // to_str(i))
call sp % write_data(meshes(i) % n_dimension, "n_dimension", &
group="tallies/mesh" // to_str(i))
call sp % write_data(meshes(i) % dimension, "dimension", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
call sp % write_data(meshes(i) % lower_left, "lower_left", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
call sp % write_data(meshes(i) % upper_right, "upper_right", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
call sp % write_data(meshes(i) % width, "width", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
end do MESH_LOOP
! Write number of meshes
call sp % write_data(n_meshes, "n_meshes", group="tallies/meshes")
if (n_meshes > 0) then
! Print list of mesh IDs
current => mesh_dict % keys()
allocate(id_array(n_meshes))
allocate(key_array(n_meshes))
i = 1
do while (associated(current))
key_array(i) = current % key
id_array(i) = current % value
! Move to next mesh
next => current % next
deallocate(current)
current => next
i = i + 1
end do
call sp % write_data(id_array, "ids", &
group="tallies/meshes", length=n_meshes)
call sp % write_data(key_array, "keys", &
group="tallies/meshes", length=n_meshes)
deallocate(key_array)
! Write information for meshes
MESH_LOOP: do i = 1, n_meshes
mesh => meshes(id_array(i))
call sp % write_data(mesh % id, "id", &
group="tallies/meshes/mesh " // trim(to_str(mesh % id)))
call sp % write_data(mesh % type, "type", &
group="tallies/meshes/mesh " // trim(to_str(mesh % id)))
call sp % write_data(mesh % n_dimension, "n_dimension", &
group="tallies/meshes/mesh " // trim(to_str(mesh % id)))
call sp % write_data(mesh % dimension, "dimension", &
group="tallies/meshes/mesh " // trim(to_str(mesh % id)), &
length=mesh % n_dimension)
call sp % write_data(mesh % lower_left, "lower_left", &
group="tallies/meshes/mesh " // trim(to_str(mesh % id)), &
length=mesh % n_dimension)
call sp % write_data(mesh % upper_right, "upper_right", &
group="tallies/meshes/mesh " // trim(to_str(mesh % id)), &
length=mesh % n_dimension)
call sp % write_data(mesh % width, "width", &
group="tallies/meshes/mesh " // trim(to_str(mesh % id)), &
length=mesh % n_dimension)
end do MESH_LOOP
deallocate(id_array)
end if
! Write number of tallies
call sp % write_data(n_tallies, "n_tallies", group="tallies")
! Write all tally information except results
TALLY_METADATA: do i = 1, n_tallies
!Get pointer to tally
t => tallies(i)
if (n_tallies > 0) then
! Write id
call sp % write_data(t % id, "id", group="tallies/tally" // to_str(i))
! Print list of tally IDs
allocate(id_array(n_tallies))
allocate(key_array(n_tallies))
! Write number of realizations
call sp % write_data(t % n_realizations, "n_realizations", &
group="tallies/tally" // to_str(i))
! Write all tally information except results
do i = 1, n_tallies
tally => tallies(i)
key_array(i) = tally % id
id_array(i) = i
end do
! Write size of each tally
call sp % write_data(t % total_score_bins, "total_score_bins", &
group="tallies/tally" // to_str(i))
call sp % write_data(t % total_filter_bins, "total_filter_bins", &
group="tallies/tally" // to_str(i))
call sp % write_data(id_array, "ids", &
group="tallies", length=n_tallies)
call sp % write_data(key_array, "keys", &
group="tallies", length=n_tallies)
! Write number of filters
call sp % write_data(t % n_filters, "n_filters", &
group="tallies/tally" // to_str(i))
deallocate(key_array)
! Write filter information
FILTER_LOOP: do j = 1, t % n_filters
! Write all tally information except results
TALLY_METADATA: do i = 1, n_tallies
! Write type of filter
call sp % write_data(t % filters(j) % type, "type", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j))
! Get pointer to tally
tally => tallies(i)
! Write number of bins for this filter
call sp % write_data(t % filters(j) % n_bins, "n_bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j))
call sp % write_data(len(tally % label), "label_size", &
group="tallies/tally " // trim(to_str(tally % id)))
if (len(tally % label) > 0) then
call sp % write_data(tally % label, "label", &
group="tallies/tally " // trim(to_str(tally % id)))
endif
! Write bins
if (t % filters(j) % type == FILTER_ENERGYIN .or. &
t % filters(j) % type == FILTER_ENERGYOUT) then
call sp % write_data(t % filters(j) % real_bins, "bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
length=size(t % filters(j) % real_bins))
else
call sp % write_data(t % filters(j) % int_bins, "bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
length=size(t % filters(j) % int_bins))
end if
call sp % write_data(tally % estimator, "estimator", &
group="tallies/tally " // trim(to_str(tally % id)))
call sp % write_data(tally % n_realizations, "n_realizations", &
group="tallies/tally " // trim(to_str(tally % id)))
call sp % write_data(tally % n_filters, "n_filters", &
group="tallies/tally " // trim(to_str(tally % id)))
end do FILTER_LOOP
! Write filter information
FILTER_LOOP: do j = 1, tally % n_filters
! Write number of nuclide bins
call sp % write_data(t % n_nuclide_bins, "n_nuclide_bins", &
group="tallies/tally" // to_str(i))
call sp % write_data(tally % filters(j) % type, "type", &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j))
call sp % write_data(tally % filters(j) % n_bins, "n_bins", &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j))
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT) then
call sp % write_data(tally % filters(j) % real_bins, "bins", &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j), &
length=size(tally % filters(j) % real_bins))
else
call sp % write_data(tally % filters(j) % int_bins, "bins", &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j), &
length=size(tally % filters(j) % int_bins))
end if
! Set up nuclide bin array and then write
allocate(temp_array(t % n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins
if (t % nuclide_bins(j) > 0) then
temp_array(j) = nuclides(t % nuclide_bins(j)) % zaid
else
temp_array(j) = t % nuclide_bins(j)
end if
end do NUCLIDE_LOOP
call sp % write_data(temp_array, "nuclide_bins", &
group="tallies/tally" // to_str(i), length=t % n_nuclide_bins)
deallocate(temp_array)
end do FILTER_LOOP
! Write number of score bins, score bins, user score bins
call sp % write_data(t % n_score_bins, "n_score_bins", &
group="tallies/tally" // to_str(i))
call sp % write_data(t % score_bins, "score_bins", &
group="tallies/tally" // to_str(i), length=t % n_score_bins)
call sp % write_data(t % n_user_score_bins, "n_user_score_bins", &
group="tallies/tally" // to_str(i))
call sp % write_data(tally % n_nuclide_bins, "n_nuclides", &
group="tallies/tally " // trim(to_str(tally % id)))
! Write explicit moment order strings for each score bin
k = 1
MOMENT_LOOP: do j = 1, t % n_user_score_bins
select case(t % score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
moment_name = 'P' // to_str(t % moment_order(k))
call sp % write_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally" // trim(to_str(i)) // "/moments")
k = k + 1
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
do n_order = 0, t % moment_order(k)
moment_name = 'P' // trim(to_str(n_order))
! Set up nuclide bin array and then write
allocate(key_array(tally % n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
key_array(j) = tally % nuclide_bins(j)
end do NUCLIDE_LOOP
call sp % write_data(key_array, "nuclides", &
group="tallies/tally " // trim(to_str(tally % id)), &
length=tally % n_nuclide_bins)
deallocate(key_array)
call sp % write_data(tally % n_score_bins, "n_score_bins", &
group="tallies/tally " // trim(to_str(tally % id)))
call sp % write_data(tally % score_bins, "score_bins", &
group="tallies/tally " // trim(to_str(tally % id)), &
length=tally % n_score_bins)
call sp % write_data(tally % n_user_score_bins, "n_user_score_bins", &
group="tallies/tally " // to_str(tally % id))
! Write explicit moment order strings for each score bin
k = 1
MOMENT_LOOP: do j = 1, tally % n_user_score_bins
select case(tally % score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
moment_name = 'P' // to_str(tally % moment_order(k))
call sp % write_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally" // trim(to_str(i)) // "/moments")
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
k = k + 1
end do
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
SCORE_TOTAL_YN)
do n_order = 0, t % moment_order(k)
do nm_order = -n_order, n_order
moment_name = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order))
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
do n_order = 0, tally % moment_order(k)
moment_name = 'P' // trim(to_str(n_order))
call sp % write_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally" // trim(to_str(i)) // "/moments")
k = k + 1
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
k = k + 1
end do
end do
case default
moment_name = ''
call sp % write_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally" // trim(to_str(i)) // "/moments")
k = k + 1
end select
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
SCORE_TOTAL_YN)
do n_order = 0, tally % moment_order(k)
do nm_order = -n_order, n_order
moment_name = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order))
call sp % write_data(moment_name, "order" // &
trim(to_str(k)), &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
k = k + 1
end do
end do
case default
moment_name = ''
call sp % write_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
k = k + 1
end select
end do MOMENT_LOOP
end do MOMENT_LOOP
end do TALLY_METADATA
end do TALLY_METADATA
! Indicate where source bank is stored in statepoint
if (source_separate) then
call sp % write_data(0, "source_present")
else
call sp % write_data(1, "source_present")
end if
end if
! Check for the no-tally-reduction method
@ -300,12 +368,12 @@ contains
TALLY_RESULTS: do i = 1, n_tallies
! Set point to current tally
t => tallies(i)
tally => tallies(i)
! Write sum and sum_sq for each bin
call sp % write_tally_result(t % results, "results", &
group="tallies/tally" // to_str(i), &
n1=size(t % results, 1), n2=size(t % results, 2))
call sp % write_tally_result(tally % results, "results", &
group="tallies/tally " // trim(to_str(tally % id)), &
n1=size(tally % results, 1), n2=size(tally % results, 2))
end do TALLY_RESULTS
@ -321,6 +389,10 @@ contains
end if
if (master .and. n_tallies > 0) then
deallocate(id_array)
end if
end subroutine write_state_point
!===============================================================================
@ -427,7 +499,10 @@ contains
#ifdef MPI
real(8) :: dummy ! temporary receive buffer for non-root reduces
#endif
type(TallyObject), pointer :: t => null()
integer, allocatable :: id_array(:)
type(ElemKeyValueII), pointer :: current
type(ElemKeyValueII), pointer :: next
type(TallyObject), pointer :: tally
type(TallyResult), allocatable :: tallyresult_temp(:,:)
! ==========================================================================
@ -484,22 +559,38 @@ contains
! Indicate that tallies are on
if (master) then
call sp % write_data(1, "tallies_present", group="tallies")
! Build list of tally IDs
current => tally_dict % keys()
allocate(id_array(n_tallies))
i = 1
do while (associated(current))
id_array(i) = current % value
! Move to next tally
next => current % next
deallocate(current)
current => next
i = i + 1
end do
end if
! Write all tally results
TALLY_RESULTS: do i = 1, n_tallies
t => tallies(i)
tally => tallies(i)
! Determine size of tally results array
m = size(t % results, 1)
n = size(t % results, 2)
m = size(tally % results, 1)
n = size(tally % results, 2)
n_bins = m*n*2
! Allocate array for storing sums and sums of squares, but
! contiguously in memory for each
allocate(tally_temp(2,m,n))
tally_temp(1,:,:) = t % results(:,:) % sum
tally_temp(2,:,:) = t % results(:,:) % sum_sq
tally_temp(1,:,:) = tally % results(:,:) % sum
tally_temp(2,:,:) = tally % results(:,:) % sum_sq
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into
@ -508,11 +599,12 @@ contains
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, &
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
#endif
! At the end of the simulation, store the results back in the
! regular TallyResults array
if (current_batch == n_batches) then
t % results(:,:) % sum = tally_temp(1,:,:)
t % results(:,:) % sum_sq = tally_temp(2,:,:)
tally % results(:,:) % sum = tally_temp(1,:,:)
tally % results(:,:) % sum_sq = tally_temp(2,:,:)
end if
! Put in temporary tally result
@ -521,8 +613,8 @@ contains
tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:)
! Write reduced tally results to file
call sp % write_tally_result(t % results, "results", &
group="tallies/tally" // to_str(i), n1=m, n2=n)
call sp % write_tally_result(tally % results, "results", &
group="tallies/tally " // trim(to_str(tally % id)), n1=m, n2=n)
! Deallocate temporary tally result
deallocate(tallyresult_temp)
@ -537,6 +629,9 @@ contains
! Deallocate temporary copy of tally results
deallocate(tally_temp)
end do TALLY_RESULTS
deallocate(id_array)
else
if (master) then
! Indicate that tallies are off
@ -552,20 +647,23 @@ contains
subroutine load_state_point()
character(MAX_FILE_LEN) :: path_temp
character(19) :: current_time
integer :: i
integer :: j
integer :: k
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
integer :: length(4)
integer :: int_array(3)
integer, allocatable :: temp_array(:)
logical :: source_present
real(8) :: real_array(3)
character(8) :: moment_name ! name of moment (e.g, P3, Y-1,1)
type(TallyObject), pointer :: t => null()
character(MAX_FILE_LEN) :: path_temp
character(19) :: current_time
character(52) :: label
integer :: i, j, k
integer :: length(4)
integer :: int_array(3)
integer, allocatable :: id_array(:)
integer, allocatable :: key_array(:)
integer :: curr_key
integer, allocatable :: temp_array(:)
logical :: source_present
real(8) :: real_array(3)
type(StructuredMesh), pointer :: mesh
type(TallyObject), pointer :: tally
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
character(8) :: moment_name ! name of moment (e.g, P3, Y-1,1)
! Write message
call write_message("Loading state point " // trim(path_state_point) &
@ -607,10 +705,22 @@ contains
! Read and overwrite run information except number of batches
call sp % read_data(run_mode, "run_mode")
call sp % read_data(n_particles, "n_particles")
call sp % read_data(int_array(1), "n_batches")
! Take maximum of statepoint n_batches and input n_batches
n_batches = max(n_batches, int_array(1))
! Read batch number to restart at
call sp % read_data(restart_batch, "current_batch")
! Check for source in statepoint if needed
call sp % read_data(int_array(1), "source_present")
if (int_array(1) == 1) then
source_present = .true.
else
source_present = .false.
end if
if (restart_batch > n_batches) then
call fatal_error("The number batches specified in settings.xml is fewer &
& than the number of batches in the given statepoint file.")
@ -634,7 +744,7 @@ contains
! Read in to see if CMFD was on
call sp % read_data(int_array(1), "cmfd_on")
! Write out CMFD info
! Read in CMFD info
if (int_array(1) == 1) then
call sp % read_data(cmfd % indices, "indices", length=4, group="cmfd")
call sp % read_data(cmfd % k_cmfd, "k_cmfd", length=restart_batch, &
@ -643,7 +753,7 @@ contains
call sp % read_data(cmfd % cmfd_src, "cmfd_src", &
length=length, group="cmfd")
call sp % read_data(cmfd % entropy, "cmfd_entropy", &
length=restart_batch, group="cmfd")
length=restart_batch, group="cmfd")
call sp % read_data(cmfd % balance, "cmfd_balance", &
length=restart_batch, group="cmfd")
call sp % read_data(cmfd % dom, "cmfd_dominance", &
@ -654,139 +764,158 @@ contains
end if
! Read number of meshes
call sp % read_data(n_meshes, "n_meshes", group="tallies")
call sp % read_data(n_meshes, "n_meshes", group="tallies/meshes")
! Read and overwrite mesh information
MESH_LOOP: do i = 1, n_meshes
call sp % read_data(meshes(i) % id, "id", &
group="tallies/mesh" // to_str(i))
call sp % read_data(meshes(i) % type, "type", &
group="tallies/mesh" // to_str(i))
call sp % read_data(meshes(i) % n_dimension, "n_dimension", &
group="tallies/mesh" // to_str(i))
call sp % read_data(meshes(i) % dimension, "dimension", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
call sp % read_data(meshes(i) % lower_left, "lower_left", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
call sp % read_data(meshes(i) % upper_right, "upper_right", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
call sp % read_data(meshes(i) % width, "width", &
group="tallies/mesh" // to_str(i), &
length=meshes(i) % n_dimension)
end do MESH_LOOP
if (n_meshes > 0) then
! Read list of mesh keys-> IDs
allocate(id_array(n_meshes))
allocate(key_array(n_meshes))
call sp % read_data(id_array, "ids", &
group="tallies/meshes", length=n_meshes)
call sp % read_data(key_array, "keys", &
group="tallies/meshes", length=n_meshes)
! Read and overwrite mesh information
MESH_LOOP: do i = 1, n_meshes
mesh => meshes(id_array(i))
curr_key = key_array(id_array(i))
call sp % read_data(mesh % id, "id", &
group="tallies/meshes/mesh " // trim(to_str(curr_key)))
call sp % read_data(mesh % type, "type", &
group="tallies/meshes/mesh " // trim(to_str(curr_key)))
call sp % read_data(mesh % n_dimension, "n_dimension", &
group="tallies/meshes/mesh " // trim(to_str(meshes(i) % id)))
call sp % read_data(mesh % dimension, "dimension", &
group="tallies/meshes/mesh " // trim(to_str(curr_key)), &
length=mesh % n_dimension)
call sp % read_data(mesh % lower_left, "lower_left", &
group="tallies/meshes/mesh " // trim(to_str(curr_key)), &
length=mesh % n_dimension)
call sp % read_data(mesh % upper_right, "upper_right", &
group="tallies/meshes/mesh " // trim(to_str(curr_key)), &
length=mesh % n_dimension)
call sp % read_data(mesh % width, "width", &
group="tallies/meshes/mesh " // trim(to_str(curr_key)), &
length=meshes(i) % n_dimension)
end do MESH_LOOP
deallocate(id_array)
deallocate(key_array)
end if
! Read and overwrite number of tallies
call sp % read_data(n_tallies, "n_tallies", group="tallies")
! Read list of tally keys-> IDs
allocate(id_array(n_tallies))
allocate(key_array(n_tallies))
call sp % read_data(id_array, "ids", group="tallies", length=n_tallies)
call sp % read_data(key_array, "keys", group="tallies", length=n_tallies)
! Read in tally metadata
TALLY_METADATA: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
tally => tallies(i)
curr_key = key_array(id_array(i))
! Read tally id
call sp % read_data(t % id, "id", group="tallies/tally" // to_str(i))
! Read number of realizations
call sp % read_data(t % n_realizations, "n_realizations", &
group="tallies/tally" // to_str(i))
! Read size of tally results
call sp % read_data(int_array(1), "total_score_bins", &
group="tallies/tally" // to_str(i))
call sp % read_data(int_array(2), "total_filter_bins", &
group="tallies/tally" // to_str(i))
! Check size of tally results array
if (int_array(1) /= t % total_score_bins .and. &
int_array(2) /= t % total_filter_bins) then
call fatal_error("Input file tally structure is different from &
&restart.")
call sp % read_data(j, "label_size", group="tallies/tally " // &
trim(to_str(curr_key)))
if (j > 0) then
call sp % read_data(label, "label", group="tallies/tally " // &
trim(to_str(curr_key)))
end if
! Read number of filters
call sp % read_data(t % n_filters, "n_filters", &
group="tallies/tally" // to_str(i))
call sp % read_data(tally % estimator, "estimator", &
group="tallies/tally " // trim(to_str(curr_key)))
call sp % read_data(tally % n_realizations, "n_realizations", &
group="tallies/tally " // trim(to_str(curr_key)))
call sp % read_data(tally % n_filters, "n_filters", &
group="tallies/tally " // trim(to_str(curr_key)))
! Read filter information
FILTER_LOOP: do j = 1, t % n_filters
! Read type of filter
call sp % read_data(t % filters(j) % type, "type", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j))
! Read number of bins for this filter
call sp % read_data(t % filters(j) % n_bins, "n_bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j))
! Read bins
if (t % filters(j) % type == FILTER_ENERGYIN .or. &
t % filters(j) % type == FILTER_ENERGYOUT) then
call sp % read_data(t % filters(j) % real_bins, "bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
length=size(t % filters(j) % real_bins))
FILTER_LOOP: do j = 1, tally % n_filters
call sp % read_data(tally % filters(j) % type, "type", &
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j))
call sp % read_data(tally % filters(j) % n_bins, "n_bins", &
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j))
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT) then
call sp % read_data(tally % filters(j) % real_bins, "bins", &
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j), &
length=size(tally % filters(j) % real_bins))
else
call sp % read_data(t % filters(j) % int_bins, "bins", &
group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), &
length=size(t % filters(j) % int_bins))
call sp % read_data(tally % filters(j) % int_bins, "bins", &
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j), &
length=size(tally % filters(j) % int_bins))
end if
end do FILTER_LOOP
! Read number of nuclide bins
call sp % read_data(t % n_nuclide_bins, "n_nuclide_bins", &
group="tallies/tally" // to_str(i))
call sp % read_data(tally % n_nuclide_bins, "n_nuclides", &
group="tallies/tally " // trim(to_str(curr_key)))
! Set up nuclide bin array and then write
allocate(temp_array(t % n_nuclide_bins))
call sp % read_data(temp_array, "nuclide_bins", &
group="tallies/tally" // to_str(i), length=t % n_nuclide_bins)
NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins
! Set up nuclide bin array and then read
allocate(temp_array(tally % n_nuclide_bins))
call sp % read_data(temp_array, "nuclides", &
group="tallies/tally " // trim(to_str(curr_key)), &
length=tally % n_nuclide_bins)
NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
if (temp_array(j) > 0) then
nuclides(t % nuclide_bins(j)) % zaid = temp_array(j)
tally % nuclide_bins(j) = temp_array(j)
else
t % nuclide_bins(j) = temp_array(j)
tally % nuclide_bins(j) = temp_array(j)
end if
end do NUCLIDE_LOOP
deallocate(temp_array)
! Write number of score bins, score bins, user score bins
call sp % read_data(t % n_score_bins, "n_score_bins", &
group="tallies/tally" // to_str(i))
call sp % read_data(t % score_bins, "score_bins", &
group="tallies/tally" // to_str(i), length=t % n_score_bins)
call sp % read_data(t % n_user_score_bins, "n_user_score_bins", &
group="tallies/tally" // to_str(i))
call sp % read_data(tally % n_score_bins, "n_score_bins", &
group="tallies/tally " // trim(to_str(curr_key)))
call sp % read_data(tally % score_bins, "score_bins", &
group="tallies/tally " // trim(to_str(curr_key)), &
length=tally % n_score_bins)
call sp % read_data(tally % n_user_score_bins, "n_user_score_bins", &
group="tallies/tally " // trim(to_str(curr_key)))
! Read explicit moment order strings for each score bin
k = 1
MOMENT_LOOP: do j = 1, t % n_user_score_bins
select case(t % score_bins(k))
MOMENT_LOOP: do j = 1, tally % n_user_score_bins
select case(tally % score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
call sp % read_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally" // trim(to_str(i)) // "/moments")
group="tallies/tally " // trim(to_str(curr_key)) // "/moments")
k = k + 1
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
do n_order = 0, t % moment_order(k)
do n_order = 0, tally % moment_order(k)
call sp % read_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally" // trim(to_str(i)) // "/moments")
group="tallies/tally " // trim(to_str(curr_key)) // "/moments")
k = k + 1
end do
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
SCORE_TOTAL_YN)
do n_order = 0, t % moment_order(k)
do n_order = 0, tally % moment_order(k)
do nm_order = -n_order, n_order
call sp % read_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally" // trim(to_str(i)) // "/moments")
group="tallies/tally " // trim(to_str(curr_key)) // &
"/moments")
k = k + 1
end do
end do
case default
call sp % read_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally" // trim(to_str(i)) // "/moments")
group="tallies/tally " // trim(to_str(curr_key)) // "/moments")
k = k + 1
end select
@ -794,14 +923,6 @@ contains
end do TALLY_METADATA
! Check for source in statepoint if needed
call sp % read_data(int_array(1), "source_present")
if (int_array(1) == 1) then
source_present = .true.
else
source_present = .false.
end if
! Check to make sure source bank is present
if (path_source_point == path_state_point .and. .not. source_present) then
call fatal_error("Source bank must be contained in statepoint restart &
@ -826,24 +947,30 @@ contains
n1=N_GLOBAL_TALLIES, n2=1)
! Check if tally results are present
call sp % read_data(int_array(1), "tallies_present", group="tallies", collect=.false.)
call sp % read_data(int_array(1), "tallies_present", &
group="tallies", collect=.false.)
! Read in sum and sum squared
if (int_array(1) == 1) then
TALLY_RESULTS: do i = 1, n_tallies
! Set pointer to tally
t => tallies(i)
tally => tallies(i)
curr_key = key_array(id_array(i))
! Read sum and sum_sq for each bin
call sp % read_tally_result(t % results, "results", &
group="tallies/tally" // to_str(i), &
n1=size(t % results, 1), n2=size(t % results, 2))
call sp % read_tally_result(tally % results, "results", &
group="tallies/tally " // trim(to_str(curr_key)), &
n1=size(tally % results, 1), n2=size(tally % results, 2))
end do TALLY_RESULTS
end if
end if
deallocate(id_array)
deallocate(key_array)
! Read source if in eigenvalue mode
if (run_mode == MODE_EIGENVALUE) then

View file

@ -182,7 +182,7 @@ contains
if (coord % next % rotated) then
! If next level is rotated, apply rotation matrix
coord % next % uvw = matmul(cells(coord % cell) % &
rotation, coord % uvw)
rotation_matrix, coord % uvw)
else
! Otherwise, copy this level's direction
coord % next % uvw = coord % uvw

View file

@ -1,48 +0,0 @@
#!/usr/bin/env python2
from xml.dom.minidom import parse
class Geometry(object):
def __init__(self, filename):
dom = parse(filename)
rootElement = dom.firstChild
cells = rootElement.getElementsByTagName('cell')
surfaces = rootElement.getElementsByTagName('surfaces')
lattices = rootElement.getElementsByTagName('lattices')
self.cells = [Cell(elem) for elem in cells]
self.surfaces = [Surface(elem) for elem in surfaces]
self.lattices = [Lattice(elem) for elem in lattices]
class Cell(object):
def __init__(self, domElement):
self.parse(domElement)
def parse(self, element):
for attribute in ['uid', 'universe', 'fill', 'material', 'surfaces']:
if element.hasAttribute(attribute):
setattr(self, attribute, element.getAttribute(attribute))
# Split strings into lists where necessary
if hasattr(self, 'surfaces'):
self.surfaces = self.surfaces.split()
class Surface(object):
def __init__(self, domElement):
self.parse(domElement)
def parse(self, element):
for attribute in ['uid', 'type', 'coeffs', 'boundary']:
if element.hasAttribute(attribute):
setattr(self, attribute, element.getAttribute(attribute))
# Split strings into lists where necessary
if hasattr(self, 'coeffs'):
self.surfaces = self.surfaces.split()

View file

@ -0,0 +1,17 @@
from openmc.element import *
from openmc.geometry import *
from openmc.nuclide import *
from openmc.material import *
from openmc.plots import *
from openmc.settings import *
from openmc.surface import *
from openmc.universe import *
from openmc.tallies import *
from openmc.cmfd import *
from openmc.executor import *
#from statepoint import *
try:
from openmc.opencg_compatible import *
except ImportError:
pass

View file

@ -0,0 +1,13 @@
import numpy as np
def is_integer(val):
return isinstance(val, (int, np.int32, np.int64))
def is_float(val):
return isinstance(val, (float, np.float32, np.float64))
def is_string(val):
return isinstance(val, (str, np.str))

View file

@ -0,0 +1,93 @@
def sort_xml_elements(tree):
# Retrieve all children of the root XML node in the tree
elements = tree.getchildren()
# Initialize empty lists for the sorted and comment elements
sorted_elements = []
# Initialize an empty set of tags (e.g., Surface, Cell, and Lattice)
tags = set()
# Find the unique tags in the tree
for element in elements:
tags.add(element.tag)
# Initialize an empty list for the comment elements
comment_elements = []
# Find the comment elements and record their ordering within the
# tree using a precedence with respect to the subsequent nodes
for index, element in enumerate(elements):
next_element = None
if 'Comment' in str(element.tag):
if index < len(elements)-1:
next_element = elements[index+1]
comment_elements.append((element, next_element))
# Now iterate over all tags and order the elements within each tag
for tag in tags:
# Retrieve all of the elements for this tag
try:
tagged_elements = tree.findall(tag)
except:
continue
# Initialize an empty list of tuples to sort (id, element)
tagged_data = []
# Retrieve the IDs for each of the elements
for element in tagged_elements:
key = element.get('id')
# If this element has an "ID" tag, append it to the list to sort
if not key is None:
tagged_data.append((key, element))
# Sort the elements according to the IDs for this tag
tagged_data.sort()
sorted_elements.extend(list(item[-1] for item in tagged_data))
# Add the comment elements while preserving the original precedence
for element, next_element in comment_elements:
index = sorted_elements.index(next_element)
sorted_elements.insert(index, element)
# Remove all of the sorted elements from the tree
for element in sorted_elements:
tree.remove(element)
# Add the sorted elements back to the tree in the proper order
tree.extend(sorted_elements)
def clean_xml_indentation(element, level=0):
"""
copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint
it basically walks your tree and adds spaces and newlines so the tree is
printed in a nice way
"""
i = "\n" + level*" "
if len(element):
if not element.text or not element.text.strip():
element.text = i + " "
if not element.tail or not element.tail.strip():
element.tail = i
for element in element:
clean_xml_indentation(element, level+1)
if not element.tail or not element.tail.strip():
element.tail = i
else:
if level and (not element.tail or not element.tail.strip()):
element.tail = i

717
src/utils/openmc/cmfd.py Normal file
View file

@ -0,0 +1,717 @@
from xml.etree import ElementTree as ET
import numpy as np
from openmc.checkvalue import *
from openmc.clean_xml import *
class CMFDMesh(object):
def __init__(self):
self._lower_left = None
self._upper_right = None
self._dimension = None
self._width = None
self._energy = None
self._albedo = None
self._map = None
@property
def lower_left(self):
return self._lower_left
@property
def upper_right(self):
return self._upper_right
@property
def dimension(self):
return self._dimension
@property
def width(self):
return self._width
@property
def energy(self):
return self._energy
@property
def albedo(self):
return self._albedo
@property
def map(self):
return self._mape
@lower_left.setter
def lower_left(self, lower_left):
if not isinstance(lower_left, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh with lower_left {0} which is ' \
'not a Python list, tuple or NumPy array'.format(lower_left)
raise ValueError(msg)
elif len(lower_left) != 2 and len(lower_left) != 3:
msg = 'Unable to set CMFD Mesh with lower_left {0} since it ' \
'must include 2 or 3 dimensions'.format(lower_left)
raise ValueError(msg)
for coord in lower_left:
if not is_integer(coord) and not is_float(coord):
msg = 'Unable to set CMFD Mesh with lower_left {0} which is ' \
'not an integer or a floating point value'.format(coord)
raise ValueError(msg)
self._lower_left = lower_left
@upper_right.setter
def upper_right(self, upper_right):
if not isinstance(upper_right, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh with upper_right {0} which is ' \
'not a Python list, tuple or NumPy array'.format(upper_right)
raise ValueError(msg)
if len(upper_right) != 2 and len(upper_right) != 3:
msg = 'Unable to set CMFD Mesh with upper_right {0} since it ' \
'must include 2 or 3 dimensions'.format(upper_right)
raise ValueError(msg)
for coord in upper_right:
if not is_integer(coord) and not is_float(coord):
msg = 'Unable to set CMFD Mesh with upper_right {0} which ' \
'is not an integer or floating point value'.format(coord)
raise ValueError(msg)
self._upper_right = upper_right
@dimension.setter
def dimension(self, dimension):
if not isinstance(dimension, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh with dimension {0} which is ' \
'not a Python list, tuple or NumPy array'.format(dimension)
raise ValueError(msg)
elif len(dimension) != 2 and len(dimension) != 3:
msg = 'Unable to set CMFD Mesh with dimension {0} since it ' \
'must include 2 or 3 dimensions'.format(dimension)
raise ValueError(msg)
for dim in dimension:
if not is_integer(dim):
msg = 'Unable to set CMFD Mesh with dimension {0} which ' \
'is a non-integer'.format(dim)
raise ValueError(msg)
self._dimension = dimension
@width.setter
def width(self, width):
if not width is None:
if not isinstance(width, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh with width {0} which ' \
'is not a Python list, tuple or NumPy array'.format(width)
raise ValueError(msg)
if len(width) != 2 and len(width) != 3:
msg = 'Unable to set CMFD Mesh with width {0} since it must ' \
'include 2 or 3 dimensions'.format(width)
raise ValueError(msg)
for dim in width:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to set CMFD Mesh with width {0} which is ' \
'not an integer or floating point value'.format(width)
raise ValueError(msg)
self._width = width
@energy.setter
def energy(self, energy):
if not isinstance(energy, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh energy to {0} which is not ' \
'a Python tuple/list or NumPy array'.format(energy)
raise ValueError(msg)
for e in energy:
if not is_integer(e) and not is_float(e):
msg = 'Unable to set CMFD Mesh energy to {0} which is not ' \
'an integer or floating point value'.format(e)
raise ValueError(msg)
elif e < 0:
msg = 'Unable to set CMFD Mesh energy to {0} which is ' \
'is a negative integer'.format(e)
raise ValueError(msg)
self._energy = energy
@albedo.setter
def albedo(self, albedo):
if not isinstance(albedo, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh albedo to {0} which is not ' \
'a Python tuple/list or NumPy array'.format(albedo)
raise ValueError(msg)
if not len(albedo) == 6:
msg = 'Unable to set CMFD Mesh albedo to {0} which is not ' \
'length 6 for +/-x,y,z'.format(albedo)
raise ValueError(msg)
for a in albedo:
if not is_integer(a) and not is_float(a):
msg = 'Unable to set CMFD Mesh albedo to {0} which is not ' \
'an integer or floating point value'.format(a)
raise ValueError(msg)
elif a < 0 or a > 1:
msg = 'Unable to set CMFD Mesh albedo to {0} which is ' \
'is not in [0,1]'.format(a)
raise ValueError(msg)
self._albedo = albedo
@map.setter
def map(self, map):
if not isinstance(map, (tuple, list, np.ndarray)):
msg = 'Unable to set CMFD Mesh map to {0} which is not ' \
'a Python tuple/list or NumPy array'.format(map)
raise ValueError(msg)
for m in map:
if m != 1 and m != 2:
msg = 'Unable to set CMFD Mesh map to {0} which is ' \
'is not 1 or 2'.format(m)
raise ValueError(msg)
self._map = map
def get_mesh_xml(self):
element = ET.Element("mesh")
if len(self._lower_left) == 2:
subelement = ET.SubElement(element, "lower_left")
subelement.text = '{0} {1}'.format(self._lower_left[0],
self._lower_left[1])
else:
subelement = ET.SubElement(element, "lower_left")
subelement.text = '{0} {1} {2}'.format(self._lower_left[0],
self._lower_left[1],
self._lower_left[2])
if not self._upper_right is None:
if len(self._upper_right) == 2:
subelement = ET.SubElement(element, "upper_right")
subelement.text = '{0} {1}'.format(self._upper_right[0],
self._upper_right[1])
else:
subelement = ET.SubElement(element, "upper_right")
subelement.text = '{0} {1} {2}'.format(self._upper_right[0],
self._upper_right[1],
self._upper_right[2])
if len(self._dimension) == 2:
subelement = ET.SubElement(element, "dimension")
subelement.text = '{0} {1}'.format(self._dimension[0],
self._dimension[1])
else:
subelement = ET.SubElement(element, "dimension")
subelement.text = '{0} {1} {2}'.format(self._dimension[0],
self._dimension[1],
self._dimension[2])
if not self._width is None:
if len(self._width) == 2:
subelement = ET.SubElement(element, "width")
subelement.text = '{0} {1}'.format(self._width[0],
self._width[1])
else:
subelement = ET.SubElement(element, "width")
subelement.text = '{0} {1} {2}'.format(self._width[0],
self._width[1],
self._width[2])
if not self._energy is None:
subelement = ET.SubElement(element, "energy")
energy = ''
for e in self._energy:
energy += '{0} '.format(e)
subelement.set("energy", energy.rstrip(' '))
if not self._albedo is None:
subelement = ET.SubElement(element, "albedo")
albedo = ''
for a in self._albedo:
albedo += '{0} '.format(a)
subelement.set("albedo", albedo.rstrip(' '))
if not self._map is None:
subelement = ET.SubElement(element, "map")
map = ''
for m in self._map:
map += '{0} '.format(m)
subelement.set("map", map.rstrip(' '))
return element
class CMFDFile(object):
def __init__(self):
self._active_flush = None
self._begin = None
self._display = None
self._feedback = None
self._inactive = None
self._inactive_flush = None
self._ksp_monitor = None
self._cmfd_mesh = None
self._norm = None
self._num_flushes = None
self._power_monitor = None
self._run_adjoint = None
self._snes_monitor = None
self._solver = None
self._write_matrices = None
self._cmfd_file = ET.Element("cmfd")
self._cmfd_mesh_element = None
@property
def active_flush(self):
return self._active_flush
@property
def begin(self):
return self._begin
@property
def display(self):
return self._display
@property
def feedback(self):
return self._feedback
@property
def inactive(self):
return self._inactive
@property
def inactive_flush(self):
return self._inactive_flush
@property
def ksp_monitor(self):
return self._ksp_monitor
@property
def cmfd_mesh(self):
return self._cmfd_mesh
@property
def norm(self):
return self._norm
@property
def num_flushes(self):
return self._num_flushes
@property
def power_monitor(self):
return self._power_monitor
@property
def run_adjoint(self):
return self._run_adjoint
@property
def snes_monitor(self):
return self._snes_monitor
@property
def solver(self):
return self._solver
@property
def write_matrices(self):
return self._write_matrices
@active_flush.setter
def active_flush(self, active_flush):
if not is_integer(active_flush):
msg = 'Unable to set CMFD active flush batch to a non-integer ' \
'value {0}'.format(active_flush)
raise ValueError(msg)
if active_flush < 0:
msg = 'Unable to set CMFD active flush batch to a negative ' \
'value {0}'.format(active_flush)
raise ValueError(msg)
self._active_flush = active_flush
@begin.setter
def begin(self, begin):
if not is_integer(begin):
msg = 'Unable to set CMFD begin batch to a non-integer ' \
'value {0}'.format(begin)
raise ValueError(msg)
if begin <= 0:
msg = 'Unable to set CMFD begin batch batch to a negative ' \
'value {0}'.format(begin)
raise ValueError(msg)
self._begin = begin
@display.setter
def display(self, display):
if not is_string(display):
msg = 'Unable to set CMFD display to a non-string ' \
'value'.format(display)
raise ValueError(msg)
if display not in ['balance', 'dominance', 'entropy', 'source']:
msg = 'Unable to set CMFD display to {0} which is ' \
'not an accepted value'.format(display)
raise ValueError(msg)
self._display = display
@feedback.setter
def feedback(self, feedback):
if not isinstance(feedback, bool):
msg = 'Unable to set CMFD feedback to {0} which is ' \
'a non-boolean value'.format(feedback)
raise ValueError(msg)
self._feedback = feedback
@inactive.setter
def inactive(self, inactive):
if not isinstance(inactive, bool):
msg = 'Unable to set CMFD inactive batch to {0} which is ' \
' a non-boolean value'.format(inactive)
raise ValueError(msg)
self._inactive = inactive
@inactive_flush.setter
def inactive_flush(self, inactive_flush):
if not is_integer(inactive_flush):
msg = 'Unable to set CMFD inactive flush batch to {0} which is ' \
'a non-integer value'.format(inactive_flush)
raise ValueError(msg)
if inactive_flush <= 0:
msg = 'Unable to set CMFD inactive flush batch to {0} which is ' \
'a negative value {0}'.format(inactive_flush)
raise ValueError(msg)
self._inactive_flush = inactive_flush
@ksp_monitor.setter
def ksp_monitor(self, ksp_monitor):
if not isinstance(ksp_monitor, bool):
msg = 'Unable to set CMFD ksp monitor to {0} which is a ' \
'non-boolean value'.format(ksp_monitor)
raise ValueError(msg)
self._ksp_monitor = ksp_monitor
@cmfd_mesh.setter
def cmfd_mesh(self, mesh):
if not isinstance(mesh, CMFDMesh):
msg = 'Unable to set CMFD mesh to {0} which is not a ' \
'CMFDMesh object'.format(mesh)
raise ValueError(msg)
self._mesh = mesh
@norm.setter
def norm(self, norm):
if not is_integer(norm) and not is_float(norm):
msg = 'Unable to set the CMFD norm to {0} which is not ' \
'an integer or floating point value'.format(norm)
raise ValueError(msg)
self._norm = norm
@snes_monitor.setter
def snum_flushes(self, num_flushes):
if not is_integer(num_flushes):
msg = 'Unable to set the CMFD number of flushes to {0} ' \
'which is not an integer value'.format(num_flushes)
raise ValueError(msg)
if num_flushes < 0:
msg = 'Unable to set CMFD number of flushes to a negative ' \
'value {0}'.format(num_flushes)
raise ValueError(msg)
self._num_flushes = num_flushes
@power_monitor.setter
def power_monitor(self, power_monitor):
if not isinstance(power_monitor, bool):
msg = 'Unable to set CMFD power monitor to {0} which is a ' \
'non-boolean value'.format(power_monitor)
raise ValueError(msg)
self._power_monitor = power_monitor
@run_adjoint.setter
def run_adjoint(self, run_adjoint):
if not isinstance(run_adjoint, bool):
msg = 'Unable to set CMFD run adjoint to {0} which is a ' \
'non-boolean value'.format(run_adjoint)
raise ValueError(msg)
self._run_adjoint = run_adjoint
@snes_monitor.setter
def snes_monitor(self, snes_monitor):
if not isinstance(snes_monitor, bool):
msg = 'Unable to set CMFD snes monitor to {0} which is a ' \
'non-boolean value'.format(snes_monitor)
raise ValueError(msg)
self._snes_monitor = snes_monitor
@solver.setter
def solver(self, solver):
if not solver in ['power', 'jfnk']:
msg = 'Unable to set CMFD solver to {0} which is not ' \
'"power" or "jfnk"'.format(solver)
raise ValueError(msg)
self._solver = solver
@write_matrices.setter
def write_matrices(self, write_matrices):
if not isinstance(write_matrices, bool):
msg = 'Unable to set CMFD write matrices to {0} which is a ' \
'non-boolean value'.format(write_matrices)
raise ValueError(msg)
self._write_matrices = write_matrices
def get_active_flush_subelement(self):
if not self._active_flush is None:
element = ET.SubElement(self._cmfd_file, "active_flush")
element.text = '{0}'.format(str(self._active_flush))
def get_begin_subelement(self):
if not self._begin is None:
element = ET.SubElement(self._cmfd_file, "begin")
element.text = '{0}'.format(str(self._begin))
def create_display_subelement(self):
if not self._display is None:
element = ET.SubElement(self._cmfd_file, "display")
element.text = '{0}'.format(str(self._display))
def create_feedback_subelement(self):
if not self._feedback is None:
element = ET.SubElement(self._cmfd_file, "feeback")
element.text = '{0}'.format(str(self._feedback).lower())
def create_inactive_subelement(self):
if not self._inactive is None:
element = ET.SubElement(self._cmfd_file, "inactive")
element.text = '{0}'.format(str(self._inactive).lower())
def create_inactive_flush_subelement(self):
if not self._inactive_flush is None:
element = ET.SubElement(self._cmfd_file, "inactive_flush")
element.text = '{0}'.format(str(self._inactive_flush))
def create_ksp_monitor_subelement(self):
if not self._ksp_monitor is None:
element = ET.SubElement(self._cmfd_file, "ksp_monitor")
element.text = '{0}'.format(str(self._ksp_monitor).lower())
def create_mesh_subelement(self):
if not self._mesh is None:
xml_element = self._mesh.get_mesh_xml()
self._cmfd_file.append(xml_element)
def create_norm_subelement(self):
if not self._num_flushes is None:
element = ET.SubElement(self._cmfd_file, "norm")
element.text = '{0}'.format(str(self._norm))
def create_num_flushes_subelement(self):
if not self._num_flushes is None:
element = ET.SubElement(self._cmfd_file, "num_flushes")
element.text = '{0}'.format(str(self._num_flushes))
def create_power_monitor_subelement(self):
if not self._power_monitor is None:
element = ET.SubElement(self._cmfd_file, "power_monitor")
element.text = '{0}'.format(str(self._power_monitor).lower())
def create_run_adjoint_subelement(self):
if not self._run_adjoint is None:
element = ET.SubElement(self._cmfd_file, "run_adjoint")
element.text = '{0}'.format(str(self._run_adjoint).lower())
def create_snes_monitor_subelement(self):
if not self._snes_monitor is None:
element = ET.SubElement(self._cmfd_file, "snes_monitor")
element.text = '{0}'.format(str(self._snes_monitor).lower())
def create_solver_subelement(self):
if not self._solver is None:
element = ET.SubElement(self._cmfd_file, "solver")
element.text = '{0}'.format(str(self._solver))
def create_write_matrices_subelement(self):
if not self._write_matrices is None:
element = ET.SubElement(self._cmfd_file, "write_matrices")
element.text = '{0}'.format(str(self._write_matrices).lower())
def export_to_xml(self):
self.create_active_flush_subelement()
self.create_begin_subelement()
self.create_display_subelement()
self.create_feedback_subelement()
self.create_inactive_subelement()
self.create_inactive_flush_subelement()
self.create_ksp_monitor_subelement()
self.create_mesh_subelement()
self.create_norm_subelement()
self.create_num_flushes_subelement()
self.create_power_monitor_subelement()
self.create_run_adjoint_subelement()
self.create_snes_monitor_subelement()
self.create_solver_subelement()
self.create_write_matrices_subelement()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._cmfd_file)
# Write the XML Tree to the cmfd.xml file
tree = ET.ElementTree(self._cmfd_file)
tree.write("cmfd.xml", xml_declaration=True,
encoding='utf-8', method="xml")

View file

@ -0,0 +1,126 @@
"""Dictionaries of integer-to-string mappings from openmc/src/constants.F90"""
SURFACE_TYPES = {1: 'x-plane',
2: 'y-plane',
3: 'z-plane',
4: 'plane',
5: 'x-cylinder',
6: 'y-cylinder',
7: 'z-cylinder',
8: 'sphere',
9: 'x-cone',
10: 'y-cone',
11: 'z-cone'}
BC_TYPES = {0: 'transmission',
1: 'vacuum',
2: 'reflective',
3: 'periodic'}
FILL_TYPES = {1: 'normal',
2: 'fill',
3: 'lattice'}
LATTICE_TYPES = {1: 'rectangular',
2: 'hexagonal'}
ESTIMATOR_TYPES = {1: 'analog',
2: 'tracklength'}
FILTER_TYPES = {1: 'universe',
2: 'material',
3: 'cell',
4: 'cellborn',
5: 'surface',
6: 'mesh',
7: 'energy',
8: 'energyout',
9: 'distribcell'}
SCORE_TYPES = {-1: 'flux',
-2: 'total',
-3: 'scatter',
-4: 'nu-scatter',
-5: 'scatter-n',
-6: 'scatter-pn',
-7: 'nu-scatter-n',
-8: 'nu-scatter-pn',
-9: 'transport',
-10: 'n1n',
-11: 'n2n',
-12: 'n3n',
-13: 'n4n',
-14: 'absorption',
-15: 'fission',
-16: 'nu-fission',
-17: 'kappa-fission',
-18: 'current',
-19: 'flux-yn',
-20: 'total-yn',
-21: 'scatter-yn',
-22: 'nu-scatter-yn',
-23: 'events',
1: '(n,total)',
2: '(n,elastic)',
4: '(n,level)',
11: '(n,2nd)',
16: '(n,2n)',
17: '(n,3n)',
18: '(n,fission)',
19: '(n,f)',
20: '(n,nf)',
21: '(n,2nf)',
22: '(n,na)',
23: '(n,n3a)',
24: '(n,2na)',
25: '(n,3na)',
28: '(n,np)',
29: '(n,n2a)',
30: '(n,2n2a)',
32: '(n,nd)',
33: '(n,nt)',
34: '(n,nHe-3)',
35: '(n,nd2a)',
36: '(n,nt2a)',
37: '(n,4n)',
38: '(n,3nf)',
41: '(n,2np)',
42: '(n,3np)',
44: '(n,n2p)',
45: '(n,npa)',
91: '(n,nc)',
101: '(n,disappear)',
102: '(n,gamma)',
103: '(n,p)',
104: '(n,d)',
105: '(n,t)',
106: '(n,3He)',
107: '(n,a)',
108: '(n,2a)',
109: '(n,3a)',
111: '(n,2p)',
112: '(n,pa)',
113: '(n,t2a)',
114: '(n,d2a)',
115: '(n,pd)',
116: '(n,pt)',
117: '(n,da)',
201: '(n,Xn)',
202: '(n,Xgamma)',
203: '(n,Xp)',
204: '(n,Xd)',
205: '(n,Xt)',
206: '(n,X3He)',
207: '(n,Xa)',
444: '(damage)',
649: '(n,pc)',
699: '(n,dc)',
749: '(n,tc)',
799: '(n,3Hec)',
849: '(n,tc)'}
SCORE_TYPES.update({MT: '(n,n' + str(MT-50) + ')' for MT in range(51,91)})
SCORE_TYPES.update({MT: '(n,p' + str(MT-600) + ')' for MT in range(600,649)})
SCORE_TYPES.update({MT: '(n,d' + str(MT-650) + ')' for MT in range(650,699)})
SCORE_TYPES.update({MT: '(n,t' + str(MT-700) + ')' for MT in range(700,749)})
SCORE_TYPES.update({MT: '(n,3He' + str(MT-750) + ')' for MT in range(750,649)})
SCORE_TYPES.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)})

View file

@ -0,0 +1,79 @@
from openmc.checkvalue import *
class Element(object):
def __init__(self, name='', xs=None):
# Initialize class attributes
self._name = ''
self._xs = None
# Set the Material class attributes
self.name = name
if not xs is None:
self.xs = xs
def __eq__(self, element2):
# Check type
if not isinstance(element2, Element):
return False
# Check name
if self._name != element2._name:
return False
# Check xs
elif self._xs != element2._xs:
return False
else:
return True
def __hash__(self):
hashable = []
hashable.append(self._name)
hashable.append(self._xs)
return hash(tuple(hashable))
@property
def xs(self):
return self._xs
@property
def name(self):
return self._name
@xs.setter
def xs(self, xs):
if not is_string(xs):
msg = 'Unable to set cross-section identifier xs for Element ' \
'with a non-string value {0}'.format(xs)
raise ValueError(msg)
self._xs = xs
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Element with a non-string ' \
'value {0}'.format(name)
raise ValueError(msg)
self._name = name
def __repr__(self):
string = 'Element - {0}\n'.format(self._name)
string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs)
return string

View file

@ -0,0 +1,79 @@
import subprocess
import os
from openmc.checkvalue import *
class Executor(object):
def __init__(self):
self._working_directory = '.'
@property
def working_directory(self):
return self._working_directory
@working_directory.setter
def working_directory(self, working_directory):
if not is_string(working_directory):
msg = 'Unable to set Executor\'s working directory to {0} ' \
'since it is not a string'.format(working_directory)
raise ValueError(msg)
elif not os.path.isdir(working_directory):
msg = 'Unable to set Executor\'s working directory to {0} ' \
'which does not exist'.format(working_directory)
raise ValueError(msg)
self._working_directory = working_directory
def plot_geometry(self, output=True):
if output:
subprocess.check_call('openmc -p', shell=True,
cwd=self._working_directory)
else:
subprocess.check_call('openmc -p', shell=True,
stdout=open(os.devnull, 'w'),
cwd=self._working_directory)
def run_simulation(self, particles=None, threads=None,
geometry_debug=False, restart_file=None,
tracks=False, mpi_procs=1, output=True):
post_args = ' '
pre_args = ''
if is_integer(particles) and particles > 0:
post_args += '-n {0} '.format(particles)
if is_integer(threads) and threads > 0:
post_args += '-s {0} '.format(threads)
if geometry_debug:
post_args += '-g '
if is_string(restart_file):
post_args += '-r {0} '.format(restart_file)
if tracks:
post_args += '-t'
if is_integer(mpi_procs) and mpi_procs > 1:
pre_args += 'mpirun -n {0} '.format(mpi_procs)
command = pre_args + 'openmc ' + post_args
if output:
subprocess.check_call(command, shell=True,
cwd=self._working_directory)
else:
subprocess.check_call(command, shell=True,
stdout=open(os.devnull, 'w'),
cwd=self._working_directory)

View file

@ -0,0 +1,174 @@
from xml.etree import ElementTree as ET
import openmc
from openmc.clean_xml import *
def reset_auto_ids():
openmc.reset_auto_material_id()
openmc.reset_auto_surface_id()
openmc.reset_auto_cell_id()
openmc.reset_auto_universe_id()
class Geometry(object):
def __init__(self):
# Initialize Geometry class attributes
self._root_universe = None
self._offsets = {}
@property
def root_universe(self):
return self._root_universe
@root_universe.setter
def root_universe(self, root_universe):
if not isinstance(root_universe, openmc.Universe):
msg = 'Unable to add root Universe {0} to Geometry since ' \
'it is not a Universe'.format(root_universe)
raise ValueError(msg)
elif root_universe._id != 0:
msg = 'Unable to add root Universe {0} to Geometry since ' \
'it has ID={1} instead of ' \
'ID=0'.format(root_universe, root_universe._id)
raise ValueError(msg)
self._root_universe = root_universe
def get_offset(self, path, filter_offset):
"""
Returns the corresponding location in the results array for a given
path and filter number.
Parameters
----------
path : list
A list of IDs that form the path to the target. It should begin
with 0 for the base universe, and should cover every universe,
cell, and lattice passed through. For the case of the lattice,
a tuple should be provided to indicate which coordinates in the
lattice should be entered. This should be in the
form: (lat_id, i_x, i_y, i_z)
filter_offset : int
An integer that specifies which offset map the filter is using
"""
# Return memoize'd offset if possible
if (path, filter_offset) in self._offsets:
offset = self._offsets[(path, filter_offset)]
# Begin recursive call to compute offset starting with the base Universe
else:
offset = self._root_universe.get_offset(path, filter_offset)
self._offsets[(path, filter_offset)] = offset
# Return the final offset
return offset
def get_all_cells(self):
return self._root_universe.get_all_cells()
def get_all_universes(self):
return self._root_universe.get_all_universes()
def get_all_nuclides(self):
nuclides = {}
materials = self.get_all_materials()
for material in materials:
nuclides.update(material.get_all_nuclides())
return nuclides
def get_all_materials(self):
material_cells = self.get_all_material_cells()
materials = set()
for cell in material_cells:
materials.add(cell._fill)
return list(materials)
def get_all_material_cells(self):
all_cells = self.get_all_cells()
material_cells = set()
for cell_id, cell in all_cells.items():
if cell._type == 'normal':
material_cells.add(cell)
return list(material_cells)
def get_all_material_universes(self):
all_universes = self.get_all_universes()
material_universes = set()
for universe_id, universe in all_universes.items():
cells = universe._cells
for cell_id, cell in cells.items():
if cell._type == 'normal':
material_universes.add(universe)
return list(material_universes)
class GeometryFile(object):
def __init__(self):
# Initialize GeometryFile class attributes
self._geometry = None
self._geometry_file = ET.Element("geometry")
@property
def geometry(self):
return self._geometry
@geometry.setter
def geometry(self, geometry):
if not isinstance(geometry, Geometry):
msg = 'Unable to set the Geometry to {0} for the GeometryFile ' \
'since it is not a Geometry object'.format(geometry)
raise ValueError(msg)
self._geometry = geometry
def export_to_xml(self):
root_universe = self._geometry._root_universe
root_universe.create_xml_subelement(self._geometry_file)
# Clean the indentation in the file to be user-readable
sort_xml_elements(self._geometry_file)
clean_xml_indentation(self._geometry_file)
# Write the XML Tree to the materials.xml file
tree = ET.ElementTree(self._geometry_file)
tree.write("geometry.xml", xml_declaration=True,
encoding='utf-8', method="xml")

View file

@ -0,0 +1,523 @@
from collections import MappingView
from copy import deepcopy
import warnings
from xml.etree import ElementTree as ET
import openmc
from openmc.checkvalue import *
from openmc.clean_xml import *
# A list of all IDs for all Materials created
MATERIAL_IDS = []
# A static variable for auto-generated Material IDs
AUTO_MATERIAL_ID = 10000
def reset_auto_material_id():
global AUTO_MATERIAL_ID, MATERIAL_IDS
AUTO_MATERIAL_ID = 10000
MATERIAL_IDS = []
# Units for density supported by OpenMC
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'at/b-cm', 'at/cm3', 'sum']
# ENDF temperatures
ENDF_TEMPS = np.array([300, 600, 700, 900, 1200, 1500])
# ENDF ZAIDs
ENDF_ZAIDS = np.array(['70c', '71c', '72c', '73c', '74c'])
# Constant for density when not needed
NO_DENSITY = 99999.
class Material(object):
def __init__(self, material_id=None, name=''):
# Initialize class attributes
self.id = material_id
self.name = name
self._density = None
self._density_units = ''
# A dictionary of Nuclides
# Keys - Nuclide names
# Values - tuple (nuclide, percent, percent type)
self._nuclides = {}
# A dictionary of Elements
# Keys - Element names
# Values - tuple (element, percent, percent type)
self._elements = {}
# If specified, a list of tuples of (table name, xs identifier)
self._sab = []
# If true, the material will be initialized as distributed
self._convert_to_distrib_comps = False
# If specified, this file will be used instead of composition values
self._distrib_otf_file = None
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@property
def density(self):
return self._density
@property
def density_units(self):
return self._density_units
@property
def convert_to_distrib_comps(self):
return self._convert_to_distrib_comps
@property
def distrib_otf_file(self):
return self._distrib_otf_file
@id.setter
def id(self, material_id):
global AUTO_MATERIAL_ID, MATERIAL_IDS
# If the Material already has an ID, remove it from global list
if hasattr(self, '_id') and not self._id is None:
MATERIAL_IDS.remove(self._id)
if material_id is None:
self._id = AUTO_MATERIAL_ID
MATERIAL_IDS.append(AUTO_MATERIAL_ID)
AUTO_MATERIAL_ID += 1
# Check that the ID is an integer and wasn't already used
elif not is_integer(material_id):
msg = 'Unable to set a non-integer Material ' \
'ID {0}'.format(material_id)
raise ValueError(msg)
elif material_id in MATERIAL_IDS:
msg = 'Unable to set Material ID to {0} since a Material with ' \
'this ID was already initialized'.format(material_id)
raise ValueError(msg)
elif material_id < 0:
msg = 'Unable to set Material ID to {0} since it must be a ' \
'non-negative integer'.format(material_id)
raise ValueError(msg)
else:
self._id = material_id
MATERIAL_IDS.append(material_id)
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Material ID={0} with a non-string ' \
'value {1}'.format(self._id, name)
raise ValueError(msg)
else:
self._name = name
def set_density(self, units, density=NO_DENSITY):
if not is_float(density):
msg = 'Unable to set the density for Material ID={0} to a ' \
'non-floating point value {1}'.format(self._id, density)
raise ValueError(msg)
elif not units in DENSITY_UNITS:
msg = 'Unable to set the density for Material ID={0} with ' \
'units {1}'.format(self._id, units)
raise ValueError(msg)
if density == NO_DENSITY and units is not 'sum':
msg = 'Unable to set the density Material ID={0} ' \
'because a density must be set when not using ' \
'sum unit'.format(self._id)
raise ValueError(msg)
self._density = density
self._density_units = units
@distrib_otf_file.setter
def distrib_otf_file(self, filename):
# TODO: remove this when distributed materials are merged
warnings.warn('This feature is not yet implemented in a release ' \
'version of openmc')
if not is_string(filename) and not filename is None:
msg = 'Unable to add OTF material file to Material ID={0} with a ' \
'non-string name {1}'.format(self._id, filename)
raise ValueError(msg)
self._distrib_otf_file = filename
@convert_to_distrib_comps.setter
def convert_to_distrib_comps(self):
# TODO: remove this when distributed materials are merged
warnings.warn('This feature is not yet implemented in a release ' \
'version of openmc')
self._convert_to_distrib_comps = True
def add_nuclide(self, nuclide, percent, percent_type='ao'):
if not isinstance(nuclide, openmc.Nuclide):
msg = 'Unable to add an Nuclide to Material ID={0} with a ' \
'non-Nuclide value {1}'.format(self._id, nuclide)
raise ValueError(msg)
elif not is_float(percent):
msg = 'Unable to add an Nuclide to Material ID={0} with a ' \
'non-floating point value {1}'.format(self._id, percent)
raise ValueError(msg)
elif not percent_type in ['ao', 'wo', 'at/g-cm']:
msg = 'Unable to add an Nuclide to Material ID={0} with a ' \
'percent type {1}'.format(self._id, percent_type)
raise ValueError(msg)
# Copy this Nuclide to separate it from the Nuclide in other Materials
nuclide = deepcopy(nuclide)
self._nuclides[nuclide._name] = (nuclide, percent, percent_type)
def remove_nuclide(self, nuclide):
if not isinstance(nuclide, openmc.Nuclide):
msg = 'Unable to remove a Nuclide {0} in Material ID={1} ' \
'since it is not a Nuclide'.format(self._id, nuclide)
raise ValueError(msg)
# If the Material contains the Nuclide, delete it
if nuclide._name in self._nuclides:
del self._nuclides[nuclide._name]
def add_element(self, element, percent, percent_type='ao'):
if not isinstance(element, openmc.Element):
msg = 'Unable to add an Element to Material ID={0} with a ' \
'non-Element value {1}'.format(self._id, element)
raise ValueError(msg)
if not is_float(percent):
msg = 'Unable to add an Element to Material ID={0} with a ' \
'non-floating point value {1}'.format(self._id, percent)
raise ValueError(msg)
if not percent_type in ['ao', 'wo']:
msg = 'Unable to add an Element to Material ID={0} with a ' \
'percent type {1}'.format(self._id, percent_type)
raise ValueError(msg)
# Copy this Element to separate it from same Element in other Materials
element = deepcopy(element)
self._elements[element._name] = (element, percent, percent_type)
def remove_element(self, element):
# If the Material contains the Element, delete it
if element._name in self._elements:
del self._elements[element._name]
def add_s_alpha_beta(self, name, xs):
if not is_string(name):
msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
'non-string table name {1}'.format(self._id, name)
raise ValueError(msg)
if not is_string(xs):
msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
'non-string cross-section identifier {1}'.format(self._id, xs)
raise ValueError(msg)
self._sab.append((name, xs))
def get_all_nuclides(self):
nuclides = {}
for nuclide_name, nuclide_tuple in self._nuclides.items():
nuclide = nuclide_tuple[0]
density = nuclide_tuple[1]
nuclides[nuclide._name] = (nuclide, density)
return nuclides
def _repr__(self):
string = 'Material\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
string += '{0: <16}{1}{2}'.format('\tDensity', '=\t', self._density)
string += ' [{0}]\n'.format(self._density_units)
string += '{0: <16}\n'.format('\tS(a,b) Tables')
for sab in self._sab:
string += '{0: <16}{1}[{2}{3}]\n'.format('\tS(a,b)', '=\t',
sab[0], sab[1])
string += '{0: <16}\n'.format('\tNuclides')
for nuclide in self._nuclides:
percent = self._nuclides[nuclide][1]
percent_type = self._nuclides[nuclide][2]
string += '{0: <16}'.format('\t{0}'.format(nuclide))
string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type)
string += '{0: <16}\n'.format('\tElements')
for element in self._elements:
percent = self._nuclides[element][1]
percent_type = self._nuclides[element][2]
string += '{0: >16}'.format('\t{0}'.format(element))
string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type)
return string
def get_nuclide_xml(self, nuclide, distrib=False):
xml_element = ET.Element("nuclide")
xml_element.set("name", nuclide[0]._name)
if not distrib:
if nuclide[2] is 'ao':
xml_element.set("ao", str(nuclide[1]))
else:
xml_element.set("wo", str(nuclide[1]))
if not nuclide[0]._xs is None:
xml_element.set("xs", nuclide[0]._xs)
return xml_element
def get_element_xml(self, element, distrib=False):
xml_element = ET.Element("element")
xml_element.set("name", str(element[0]._name))
if not distrib:
if element[2] is 'ao':
xml_element.set("ao", str(element[1]))
else:
xml_element.set("wo", str(element[1]))
return xml_element
def get_nuclides_xml(self, nuclides, distrib=False):
xml_elements = []
for nuclide in nuclides.values():
xml_elements.append(self.get_nuclide_xml(nuclide, distrib))
return xml_elements
def get_elements_xml(self, elements, distrib=False):
xml_elements = []
for element in elements.values():
xml_elements.append(self.get_element_xml(element, distrib))
return xml_elements
def get_material_xml(self):
# Create Material XML element
element = ET.Element("material")
element.set("id", str(self._id))
# Create density XML subelement
subelement = ET.SubElement(element, "density")
if self._density_units is not 'sum':
subelement.set("value", str(self._density))
subelement.set("units", self._density_units)
if not self._convert_to_distrib_comps:
# Create nuclide XML subelements
subelements = self.get_nuclides_xml(self._nuclides)
for subelement in subelements:
element.append(subelement)
# Create element XML subelements
subelements = self.get_elements_xml(self._elements)
for subelement in subelements:
element.append(subelement)
else:
subelement = ET.SubElement(element, "compositions")
comps = []
allnucs = self._nuclides.values() + self._elements.values()
dist_per_type = allnucs[0][2]
for nuc, per, typ in allnucs:
if not typ == dist_per_type:
msg = 'All nuclides and elements in a distributed ' \
'material must have the same type, either ao or wo'
raise ValueError(msg)
comps.append(per)
if self._distrib_otf_file is None:
# Create values and units subelements
subsubelement = ET.SubElement(subelement, "values")
subsubelement.text = ' '.join([str(c) for c in comps])
subsubelement = ET.SubElement(subelement, "units")
subsubelement.text = dist_per_type
else:
# Specify the materials file
subsubelement = ET.SubElement(subelement, "otf_file_path")
subsubelement.text = self._distrib_otf_file
# Create nuclide XML subelements
subelements = self.get_nuclides_xml(self._nuclides, distrib=True)
for subelement_nuc in subelements:
subelement.append(subelement_nuc)
# Create element XML subelements
subelements = self.get_elements_xml(self._elements, distrib=True)
for subelement_ele in subelements:
subelement.append(subelement_ele)
if len(self._sab) > 0:
for sab in self._sab:
subelement = ET.SubElement(element, "sab")
subelement.set("name", sab[0])
subelement.set("xs", sab[1])
return element
class MaterialsFile(object):
def __init__(self):
# Initialize MaterialsFile class attributes
self._materials = []
self._default_xs = None
self._materials_file = ET.Element("materials")
@property
def default_xs(self):
return self._default_xs
@default_xs.setter
def default_xs(self, xs):
if not is_string(xs):
msg = 'Unable to set default xs to a non-string value'.format(xs)
raise ValueError(msg)
self._default_xs = xs
def add_material(self, material):
if not isinstance(material, Material):
msg = 'Unable to add a non-Material {0} to the ' \
'MaterialsFile'.format(material)
raise ValueError(msg)
self._materials.append(material)
def add_materials(self, materials):
if not isinstance(materials, (tuple, list, MappingView)):
msg = 'Unable to create OpenMC materials.xml file from {0} which ' \
'is not a Python tuple/list'.format(materials)
raise ValueError(msg)
for material in materials:
self.add_material(material)
def remove_materials(self, material):
if not isinstance(material, Material):
msg = 'Unable to remove a non-Material {0} from the ' \
'MaterialsFile'.format(material)
raise ValueError(msg)
self._materials.remove(material)
def create_material_subelements(self):
subelement = ET.SubElement(self._materials_file, "default_xs")
if not self._default_xs is None:
subelement.text = self._default_xs
for material in self._materials:
xml_element = material.get_material_xml()
if len(material._name) > 0:
self._materials_file.append(ET.Comment(material._name))
self._materials_file.append(xml_element)
def export_to_xml(self):
self.create_material_subelements()
# Clean the indentation in the file to be user-readable
sort_xml_elements(self._materials_file)
clean_xml_indentation(self._materials_file)
# Write the XML Tree to the materials.xml file
tree = ET.ElementTree(self._materials_file)
tree.write("materials.xml", xml_declaration=True,
encoding='utf-8', method="xml")

View file

@ -0,0 +1,96 @@
from openmc.checkvalue import *
class Nuclide(object):
def __init__(self, name='', xs=None):
# Initialize class attributes
self._name = ''
self._xs = None
self._zaid = None
# Set the Material class attributes
self.name = name
if not xs is None:
self.xs = xs
def __eq__(self, nuclide2):
# Check type
if not isinstance(nuclide2, Nuclide):
return False
# Check name
elif self._name != nuclide2._name:
return False
# Check xs
elif self._xs != nuclide2._xs:
return False
else:
return True
def __hash__(self):
return hash((self._name, self._xs))
@property
def name(self):
return self._name
@property
def xs(self):
return self._xs
@property
def zaid(self):
return self._zaid
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Nuclide with a non-string ' \
'value {0}'.format(name)
raise ValueError(msg)
self._name = name
@xs.setter
def xs(self, xs):
if not is_string(xs):
msg = 'Unable to set cross-section identifier xs for Nuclide ' \
'with a non-string value {0}'.format(xs)
raise ValueError(msg)
self._xs = xs
@zaid.setter
def zaid(self, zaid):
if not is_integer(zaid):
msg = 'Unable to set zaid for Nuclide ' \
'with a non-integer {0}'.format(zaid)
raise ValueError(msg)
self._zaid = zaid
def __repr__(self):
string = 'Nuclide - {0}\n'.format(self._name)
string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs)
if self._zaid is not None:
string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self._zaid)
return string

View file

@ -0,0 +1,858 @@
import copy
import numpy as np
import opencg
import openmc
# A dictionary of all OpenMC Materials created
# Keys - Material IDs
# Values - Materials
OPENMC_MATERIALS = {}
# A dictionary of all OpenCG Materials created
# Keys - Material IDs
# Values - Materials
OPENCG_MATERIALS = {}
# A dictionary of all OpenMC Surfaces created
# Keys - Surface IDs
# Values - Surfaces
OPENMC_SURFACES = {}
# A dictionary of all OpenCG Surfaces created
# Keys - Surface IDs
# Values - Surfaces
OPENCG_SURFACES = {}
# A dictionary of all OpenMC Cells created
# Keys - Cell IDs
# Values - Cells
OPENMC_CELLS = {}
# A dictionary of all OpenCG Cells created
# Keys - Cell IDs
# Values - Cells
OPENCG_CELLS = {}
# A dictionary of all OpenMC Universes created
# Keys - Universes IDs
# Values - Universes
OPENMC_UNIVERSES = {}
# A dictionary of all OpenCG Universes created
# Keys - Universes IDs
# Values - Universes
OPENCG_UNIVERSES = {}
# A dictionary of all OpenMC Lattices created
# Keys - Lattice IDs
# Values - Lattices
OPENMC_LATTICES = {}
# A dictionary of all OpenCG Lattices created
# Keys - Lattice IDs
# Values - Lattices
OPENCG_LATTICES = {}
def get_opencg_material(openmc_material):
if not isinstance(openmc_material, openmc.Material):
msg = 'Unable to create an OpenCG Material from {0} ' \
'which is not an OpenMC Material'.format(openmc_material)
raise ValueError(msg)
global OPENCG_MATERIALS
material_id = openmc_material._id
# If this Material was already created, use it
if material_id in OPENCG_MATERIALS:
return OPENCG_MATERIALS[material_id]
# Create an OpenCG Material to represent this OpenMC Material
name = openmc_material._name
opencg_material = opencg.Material(material_id=material_id, name=name)
# Add the OpenMC Material to the global collection of all OpenMC Materials
OPENMC_MATERIALS[material_id] = openmc_material
# Add the OpenCG Material to the global collection of all OpenCG Materials
OPENCG_MATERIALS[material_id] = opencg_material
return opencg_material
def get_openmc_material(opencg_material):
if not isinstance(opencg_material, opencg.Material):
msg = 'Unable to create an OpenMC Material from {0} ' \
'which is not an OpenCG Material'.format(opencg_material)
raise ValueError(msg)
global OPENMC_MATERIALS
material_id = opencg_material._id
# If this Material was already created, use it
if material_id in OPENMC_MATERIALS:
return OPENMC_MATERIALS[material_id]
# Create an OpenMC Material to represent this OpenCG Material
name = opencg_material._name
openmc_material = openmc.Material(material_id=material_id, name=name)
# Add the OpenMC Material to the global collection of all OpenMC Materials
OPENMC_MATERIALS[material_id] = openmc_material
# Add the OpenCG Material to the global collection of all OpenCG Materials
OPENCG_MATERIALS[material_id] = opencg_material
return openmc_material
def is_opencg_surface_compatible(opencg_surface):
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to check if OpenCG Surface is compatible' \
'since {0} is not a Surface'.format(opencg_surface)
raise ValueError(msg)
if opencg_surface._type in ['x-squareprism',
'y-squareprism', 'z-squareprism']:
return False
else:
return True
def get_opencg_surface(openmc_surface):
if not isinstance(openmc_surface, openmc.Surface):
msg = 'Unable to create an OpenCG Surface from {0} ' \
'which is not an OpenMC Surface'.format(openmc_surface)
raise ValueError(msg)
global OPENCG_SURFACES
surface_id = openmc_surface._id
# If this Material was already created, use it
if surface_id in OPENCG_SURFACES:
return OPENCG_SURFACES[surface_id]
# Create an OpenCG Surface to represent this OpenMC Surface
name = openmc_surface._name
# Correct for OpenMC's syntax for Surfaces dividing Cells
boundary = openmc_surface._boundary_type
if boundary == 'transmission':
boundary = 'interface'
opencg_surface = None
if openmc_surface._type == 'plane':
A = openmc_surface._coeffs['A']
B = openmc_surface._coeffs['B']
C = openmc_surface._coeffs['C']
D = openmc_surface._coeffs['D']
opencg_surface = opencg.Plane(surface_id, name, boundary, A, B, C, D)
elif openmc_surface._type == 'x-plane':
x0 = openmc_surface._coeffs['x0']
opencg_surface = opencg.XPlane(surface_id, name, boundary, x0)
elif openmc_surface._type == 'y-plane':
y0 = openmc_surface._coeffs['y0']
opencg_surface = opencg.YPlane(surface_id, name, boundary, y0)
elif openmc_surface._type == 'z-plane':
z0 = openmc_surface._coeffs['z0']
opencg_surface = opencg.ZPlane(surface_id, name, boundary, z0)
elif openmc_surface._type == 'x-cylinder':
y0 = openmc_surface._coeffs['y0']
z0 = openmc_surface._coeffs['z0']
R = openmc_surface._coeffs['R']
opencg_surface = opencg.XCylinder(surface_id, name,
boundary, y0, z0, R)
elif openmc_surface._type == 'y-cylinder':
x0 = openmc_surface._coeffs['x0']
z0 = openmc_surface._coeffs['z0']
R = openmc_surface._coeffs['R']
opencg_surface = opencg.YCylinder(surface_id, name,
boundary, x0, z0, R)
elif openmc_surface._type == 'z-cylinder':
x0 = openmc_surface._coeffs['x0']
y0 = openmc_surface._coeffs['y0']
R = openmc_surface._coeffs['R']
opencg_surface = opencg.ZCylinder(surface_id, name,
boundary, x0, y0, R)
# Add the OpenMC Surface to the global collection of all OpenMC Surfaces
OPENMC_SURFACES[surface_id] = openmc_surface
# Add the OpenCG Surface to the global collection of all OpenCG Surfaces
OPENCG_SURFACES[surface_id] = opencg_surface
return opencg_surface
def get_openmc_surface(opencg_surface):
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from {0} which ' \
'is not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
global openmc_surface
surface_id = opencg_surface._id
# If this Surface was already created, use it
if surface_id in OPENMC_SURFACES:
return OPENMC_SURFACES[surface_id]
# Create an OpenMC Surface to represent this OpenCG Surface
name = opencg_surface._name
# Correct for OpenMC's syntax for Surfaces dividing Cells
boundary = opencg_surface._boundary_type
if boundary == 'interface':
boundary = 'transmission'
if opencg_surface._type == 'plane':
A = opencg_surface._coeffs['A']
B = opencg_surface._coeffs['B']
C = opencg_surface._coeffs['C']
D = opencg_surface._coeffs['D']
openmc_surface = openmc.Plane(surface_id, boundary, A, B, C, D, name)
elif opencg_surface._type == 'x-plane':
x0 = opencg_surface._coeffs['x0']
openmc_surface = openmc.XPlane(surface_id, boundary, x0, name)
elif opencg_surface._type == 'y-plane':
y0 = opencg_surface._coeffs['y0']
openmc_surface = openmc.YPlane(surface_id, boundary, y0, name)
elif opencg_surface._type == 'z-plane':
z0 = opencg_surface._coeffs['z0']
openmc_surface = openmc.ZPlane(surface_id, boundary, z0, name)
elif opencg_surface._type == 'x-cylinder':
y0 = opencg_surface._coeffs['y0']
z0 = opencg_surface._coeffs['z0']
R = opencg_surface._coeffs['R']
openmc_surface = openmc.XCylinder(surface_id, boundary, y0, z0, R, name)
elif opencg_surface._type == 'y-cylinder':
x0 = opencg_surface._coeffs['x0']
z0 = opencg_surface._coeffs['z0']
R = opencg_surface._coeffs['R']
openmc_surface = openmc.YCylinder(surface_id, boundary, x0, z0, R, name)
elif opencg_surface._type == 'z-cylinder':
x0 = opencg_surface._coeffs['x0']
y0 = opencg_surface._coeffs['y0']
R = opencg_surface._coeffs['R']
openmc_surface = openmc.ZCylinder(surface_id, boundary, x0, y0, R, name)
else:
msg = 'Unable to create an OpenMC Surface from an OpenCG ' \
'Surface of type {0} since it is not a compatible ' \
'Surface type in OpenMC'.format(opencg_surface._type)
raise ValueError(msg)
# Add the OpenMC Surface to the global collection of all OpenMC Surfaces
OPENMC_SURFACES[surface_id] = openmc_surface
# Add the OpenCG Surface to the global collection of all OpenCG Surfaces
OPENCG_SURFACES[surface_id] = opencg_surface
return openmc_surface
def get_compatible_opencg_surfaces(opencg_surface):
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from {0} which ' \
'is not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
global OPENMC_SURFACES
surface_id = opencg_surface._id
# If this Surface was already created, use it
if surface_id in OPENMC_SURFACES:
return OPENMC_SURFACES[surface_id]
# Create an OpenMC Surface to represent this OpenCG Surface
name = opencg_surface._name
boundary = opencg_surface._boundary_type
if opencg_surface._type == 'x-squareprism':
y0 = opencg_surface._coeffs['y0']
z0 = opencg_surface._coeffs['z0']
R = opencg_surface._coeffs['R']
# Create a list of the four planes we need
left = opencg.YPlane(name=name, boundary=boundary, y0=y0-R)
right = opencg.YPlane(name=name, boundary=boundary, y0=y0+R)
bottom = opencg.ZPlane(name=name, boundary=boundary, z0=z0-R)
top = opencg.ZPlane(name=name, boundary=boundary, z0=z0+R)
surfaces = [left, right, bottom, top]
elif opencg_surface._type == 'y-squareprism':
x0 = opencg_surface._coeffs['x0']
z0 = opencg_surface._coeffs['z0']
R = opencg_surface._coeffs['R']
# Create a list of the four planes we need
left = opencg.XPlane(name=name, boundary=boundary, x0=x0-R)
right = opencg.XPlane(name=name, boundary=boundary, x0=x0+R)
bottom = opencg.ZPlane(name=name, boundary=boundary, z0=z0-R)
top = opencg.ZPlane(name=name, boundary=boundary, z0=z0+R)
surfaces = [left, right, bottom, top]
elif opencg_surface._type == 'z-squareprism':
x0 = opencg_surface._coeffs['x0']
y0 = opencg_surface._coeffs['y0']
R = opencg_surface._coeffs['R']
# Create a list of the four planes we need
left = opencg.XPlane(name=name, boundary=boundary, x0=x0-R)
right = opencg.XPlane(name=name, boundary=boundary, x0=x0+R)
bottom = opencg.YPlane(name=name, boundary=boundary, y0=y0-R)
top = opencg.YPlane(name=name, boundary=boundary, y0=y0+R)
surfaces = [left, right, bottom, top]
else:
msg = 'Unable to create a compatible OpenMC Surface an OpenCG ' \
'Surface of type {0} since it already a compatible ' \
'Surface type in OpenMC'.format(opencg_surface._type)
raise ValueError(msg)
# Add the OpenMC Surface(s) to the global collection of all OpenMC Surfaces
OPENMC_SURFACES[surface_id] = surfaces
# Add the OpenCG Surface to the global collection of all OpenCG Surfaces
OPENCG_SURFACES[surface_id] = opencg_surface
return surfaces
def get_opencg_cell(openmc_cell):
if not isinstance(openmc_cell, openmc.Cell):
msg = 'Unable to create an OpenCG Cell from {0} which ' \
'is not an OpenMC Cell'.format(openmc_cell)
raise ValueError(msg)
global OPENCG_CELLS
cell_id = openmc_cell._id
# If this Cell was already created, use it
if cell_id in OPENCG_CELLS:
return OPENCG_CELLS[cell_id]
# Create an OpenCG Cell to represent this OpenMC Cell
name = openmc_cell._name
opencg_cell = opencg.Cell(cell_id, name)
fill = openmc_cell._fill
if (openmc_cell._type == 'normal'):
opencg_cell.setFill(get_opencg_material(fill))
elif (openmc_cell._type == 'fill'):
opencg_cell.setFill(get_opencg_universe(fill))
else:
opencg_cell.setFill(get_opencg_lattice(fill))
if openmc_cell._rotation is not None:
opencg_cell.setRotation(openmc_cell._rotation)
if openmc_cell._translation is not None:
opencg_cell.setTranslation(openmc_cell._translation)
surfaces = openmc_cell._surfaces
for surface_id in surfaces:
surface = surfaces[surface_id][0]
halfspace = surfaces[surface_id][1]
opencg_cell.addSurface(get_opencg_surface(surface), halfspace)
# Add the OpenMC Cell to the global collection of all OpenMC Cells
OPENMC_CELLS[cell_id] = openmc_cell
# Add the OpenCG Cell to the global collection of all OpenCG Cells
OPENCG_CELLS[cell_id] = opencg_cell
return opencg_cell
def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create compatible OpenMC Cell from {0} which ' \
'is not an OpenCG Cell'.format(opencg_cell)
raise ValueError(msg)
elif not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create compatible OpenMC Cell since {0} is ' \
'not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
elif not halfspace in [-1, +1]:
msg = 'Unable to create compatible Cell since {0}' \
'is not a +/-1 halfspace'.format(halfspace)
raise ValueError(msg)
# Initialize an empty list for the new compatible cells
compatible_cells = []
# SquarePrism Surfaces
if opencg_surface._type in ['x-squareprism',
'y-squareprism', 'z-squareprism']:
# Get the compatible Surfaces (XPlanes and YPlanes)
compatible_surfaces = get_compatible_opencg_surfaces(opencg_surface)
opencg_cell.removeSurface(opencg_surface)
# If Cell is inside SquarePrism, add "inside" of Surface halfspaces
if halfspace == -1:
opencg_cell.addSurface(compatible_surfaces[0], +1)
opencg_cell.addSurface(compatible_surfaces[1], -1)
opencg_cell.addSurface(compatible_surfaces[2], +1)
opencg_cell.addSurface(compatible_surfaces[3], -1)
compatible_cells.append(opencg_cell)
# If Cell is outside SquarePrism, add "outside" of Surface halfspaces
else:
# Create 8 Cell clones to represent each of the disjoint planar
# Surface halfspace intersections
num_clones = 8
for clone_id in range(num_clones):
# Create a cloned OpenCG Cell with Surfaces compatible with OpenMC
clone = opencg_cell.clone()
compatible_cells.append(clone)
# Top left subcell - add left XPlane, top YPlane
if clone_id == 0:
clone.add_surface(compatible_surfaces[0], -1)
clone.add_surface(compatible_surfaces[3], +1)
# Top center subcell - add top YPlane, left/right XPlanes
elif clone_id == 1:
clone.add_surface(compatible_surfaces[0], +1)
clone.add_surface(compatible_surfaces[1], -1)
clone.add_surface(compatible_surfaces[3], +1)
# Top right subcell - add top YPlane, right XPlane
elif clone_id == 2:
clone.add_surface(compatible_surfaces[1], +1)
clone.add_surface(compatible_surfaces[3], +1)
# Right center subcell - add right XPlane, top/bottom YPlanes
elif clone_id == 3:
clone.add_surface(compatible_surfaces[1], +1)
clone.add_surface(compatible_surfaces[3], -1)
clone.add_surface(compatible_surfaces[2], +1)
# Bottom right subcell - add right XPlane, bottom YPlane
elif clone_id == 4:
clone.add_surface(compatible_surfaces[1], +1)
clone.add_surface(compatible_surfaces[2], -1)
# Bottom center subcell - add bottom YPlane, left/right XPlanes
elif clone_id == 5:
clone.add_surface(compatible_surfaces[0], +1)
clone.add_surface(compatible_surfaces[1], -1)
clone.add_surface(compatible_surfaces[2], -1)
# Bottom left subcell - add bottom YPlane, left XPlane
elif clone_id == 6:
clone.add_surface(compatible_surfaces[0], -1)
clone.add_surface(compatible_surfaces[2], -1)
# Left center subcell - add left XPlane, top/bottom YPlanes
elif clone_id == 7:
clone.add_surface(compatible_surfaces[0], -1)
clone.add_surface(compatible_surfaces[3], -1)
clone.add_surface(compatible_surfaces[2], +1)
# Remove redundant Surfaces from the Cells
for cell in compatible_cells:
cell.removeRedundantSurfaces()
# Return the list of OpenMC compatible OpenCG Cells
return compatible_cells
def make_opencg_cells_compatible(opencg_universe):
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to make compatible OpenCG Cells for {0} which ' \
'is not an OpenCG Universe'.format(opencg_universe)
raise ValueError(msg)
# Check all OpenCG Cells in this Universe for compatibility with OpenMC
opencg_cells = opencg_universe._cells
for cell_id, opencg_cell in opencg_cells.items():
# Check each of the OpenCG Surfaces for OpenMC compatibility
surfaces = opencg_cell._surfaces
for surface_id in surfaces:
surface = surfaces[surface_id][0]
halfspace = surfaces[surface_id][1]
# If this Surface is not compatible with OpenMC, create compatible
# OpenCG cells with a compatible version of this OpenCG Surface
if not is_opencg_surface_compatible(surface):
# Get one or more OpenCG Cells that are compatible with OpenMC
# NOTE: This does not necessarily make OpenCG fully compatible.
# It only removes the incompatible Surface and replaces it with
# compatible OpenCG Surface(s). The recursive call at the end
# of this block is necessary in the event that there are more
# incompatible Surfaces in this Cell that are not accounted for.
cells = get_compatible_opencg_cells(opencg_cell,
surface, halfspace)
# Remove the non-compatible OpenCG Cell from the Universe
opencg_universe.removeCell(opencg_cell)
# Add the compatible OpenCG Cells to the Universe
opencg_universe.addCells(cells)
# Make recursive call to look at the updated state of the
# OpenCG Universe and return
return make_opencg_cells_compatible(opencg_universe)
# If all OpenCG Cells in the OpenCG Universe are compatible, return
return
def get_openmc_cell(opencg_cell):
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create an OpenMC Cell from {0} which ' \
'is not an OpenCG Cell'.format(opencg_cell)
raise ValueError(msg)
global OPENMC_CELLS
cell_id = opencg_cell._id
# If this Cell was already created, use it
if cell_id in OPENMC_CELLS:
return OPENMC_CELLS[cell_id]
# Create an OpenCG Cell to represent this OpenMC Cell
name = opencg_cell._name
openmc_cell = openmc.Cell(cell_id, name)
fill = opencg_cell._fill
if (opencg_cell._type == 'universe'):
openmc_cell.fill = get_openmc_universe(fill)
elif (opencg_cell._type == 'lattice'):
openmc_cell.fill = get_openmc_lattice(fill)
else:
openmc_cell.fill = get_openmc_material(fill)
if opencg_cell._rotation:
rotation = np.asarray(opencg_cell._rotation, dtype=np.int)
openmc_cell.rotation = rotation
if opencg_cell._translation:
translation = np.asarray(opencg_cell._translation, dtype=np.float64)
openmc_cell.setTranslation(translation)
surfaces = opencg_cell._surfaces
for surface_id in surfaces:
surface = surfaces[surface_id][0]
halfspace = surfaces[surface_id][1]
openmc_cell.add_surface(get_openmc_surface(surface), halfspace)
# Add the OpenMC Cell to the global collection of all OpenMC Cells
OPENMC_CELLS[cell_id] = openmc_cell
# Add the OpenCG Cell to the global collection of all OpenCG Cells
OPENCG_CELLS[cell_id] = opencg_cell
return openmc_cell
def get_opencg_universe(openmc_universe):
if not isinstance(openmc_universe, openmc.Universe):
msg = 'Unable to create an OpenCG Universe from {0} which ' \
'is not an OpenMC Universe'.format(openmc_universe)
raise ValueError(msg)
global OPENCG_UNIVERSES
universe_id = openmc_universe._id
# If this Universe was already created, use it
if universe_id in OPENCG_UNIVERSES:
return OPENCG_UNIVERSES[universe_id]
# Create an OpenCG Universe to represent this OpenMC Universe
name = openmc_universe._name
opencg_universe = opencg.Universe(universe_id, name)
# Convert all OpenMC Cells in this Universe to OpenCG Cells
openmc_cells = openmc_universe._cells
for cell_id, openmc_cell in openmc_cells.items():
opencg_cell = get_opencg_cell(openmc_cell)
opencg_universe.addCell(opencg_cell)
# Add the OpenMC Universe to the global collection of all OpenMC Universes
OPENMC_UNIVERSES[universe_id] = openmc_universe
# Add the OpenCG Universe to the global collection of all OpenCG Universes
OPENCG_UNIVERSES[universe_id] = opencg_universe
return opencg_universe
def get_openmc_universe(opencg_universe):
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to create an OpenMC Universe from {0} which ' \
'is not an OpenCG Universe'.format(opencg_universe)
raise ValueError(msg)
global OPENMC_UNIVERSES
universe_id = opencg_universe._id
# If this Universe was already created, use it
if universe_id in OPENMC_UNIVERSES:
return OPENMC_UNIVERSES[universe_id]
# Make all OpenCG Cells and Surfaces in this Universe compatible with OpenMC
make_opencg_cells_compatible(opencg_universe)
# Create an OpenMC Universe to represent this OpenCSg Universe
name = opencg_universe._name
openmc_universe = openmc.Universe(universe_id, name)
# Convert all OpenCG Cells in this Universe to OpenMC Cells
opencg_cells = opencg_universe._cells
for cell_id, opencg_cell in opencg_cells.items():
openmc_cell = get_openmc_cell(opencg_cell)
openmc_universe.add_cell(openmc_cell)
# Add the OpenMC Universe to the global collection of all OpenMC Universes
OPENMC_UNIVERSES[universe_id] = openmc_universe
# Add the OpenCG Universe to the global collection of all OpenCG Universes
OPENCG_UNIVERSES[universe_id] = opencg_universe
return openmc_universe
def get_opencg_lattice(openmc_lattice):
if not isinstance(openmc_lattice, openmc.Lattice):
msg = 'Unable to create an OpenCG Lattice from {0} which ' \
'is not an OpenMC Lattice'.format(openmc_lattice)
raise ValueError(msg)
global OPENCG_LATTICES
lattice_id = openmc_lattice._id
# If this Lattice was already created, use it
if lattice_id in OPENCG_LATTICES:
return OPENCG_LATTICES[lattice_id]
# Create an OpenCG Lattice to represent this OpenMC Lattice
name = openmc_lattice.name
dimension = openmc_lattice.dimension
pitch = openmc_lattice.pitch
lower_left = openmc_lattice.lower_left
universes = openmc_lattice.universes
if len(pitch) == 2:
new_pitch = np.ones(3, dtype=np.float64)
new_pitch[:2] = pitch
pitch = new_pitch
if len(lower_left) == 2:
new_lower_left = np.ones(3, dtype=np.float64)
new_lower_left[:2] = lower_left
lower_left = new_lower_left
# Initialize an empty array for the OpenCG nested Universes in this Lattice
universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), \
dtype=opencg.Universe)
# Create OpenCG Universes for each unique nested Universe in this Lattice
unique_universes = openmc_lattice.get_unique_universes()
for universe_id, universe in unique_universes.items():
unique_universes[universe_id] = get_opencg_universe(universe)
# Build the nested Universe array
for z in range(dimension[2]):
for y in range(dimension[1]):
for x in range(dimension[0]):
universe_id = universes[x][dimension[1]-y-1][z]._id
universe_array[z][y][x] = unique_universes[universe_id]
opencg_lattice = opencg.Lattice(lattice_id, name)
opencg_lattice.setDimension(dimension)
opencg_lattice.setWidth(pitch)
opencg_lattice.setUniverses(universe_array)
offset = np.array(lower_left, dtype=np.float64) - \
((np.array(pitch, dtype=np.float64) * \
np.array(dimension, dtype=np.float64))) / -2.0
opencg_lattice.setOffset(offset)
# Add the OpenMC Lattice to the global collection of all OpenMC Lattices
OPENMC_LATTICES[lattice_id] = openmc_lattice
# Add the OpenCG Lattice to the global collection of all OpenCG Lattices
OPENCG_LATTICES[lattice_id] = opencg_lattice
return opencg_lattice
def get_openmc_lattice(opencg_lattice):
if not isinstance(opencg_lattice, opencg.Lattice):
msg = 'Unable to create an OpenMC Lattice from {0} which ' \
'is not an OpenCG Lattice'.format(opencg_lattice)
raise ValueError(msg)
global OPENMC_LATTICES
lattice_id = opencg_lattice._id
# If this Lattice was already created, use it
if lattice_id in OPENMC_LATTICES:
return OPENMC_LATTICES[lattice_id]
dimension = opencg_lattice._dimension
width = opencg_lattice._width
offset = opencg_lattice._offset
universes = opencg_lattice._universes
# Initialize an empty array for the OpenMC nested Universes in this Lattice
universe_array = np.ndarray(tuple(np.array(dimension)), \
dtype=openmc.Universe)
# Create OpenMC Universes for each unique nested Universe in this Lattice
unique_universes = opencg_lattice.getUniqueUniverses()
for universe_id, universe in unique_universes.items():
unique_universes[universe_id] = get_openmc_universe(universe)
# Build the nested Universe array
for z in range(dimension[2]):
for y in range(dimension[1]):
for x in range(dimension[0]):
universe_id = universes[z][y][x]._id
universe_array[x][y][z] = unique_universes[universe_id]
# Reverse y-dimension in array to match ordering in OpenCG
universe_array = universe_array[:,::-1,:]
lower_left = np.array(offset, dtype=np.float64) + \
((np.array(width, dtype=np.float64) * \
np.array(dimension, dtype=np.float64))) / -2.0
openmc_lattice = openmc.RectLattice(lattice_id=lattice_id)
openmc_lattice.dimension = dimension
openmc_lattice.pitch = width
openmc_lattice.universes = universe_array
openmc_lattice.lower_left = lower_left
# Add the OpenMC Lattice to the global collection of all OpenMC Lattices
OPENMC_LATTICES[lattice_id] = openmc_lattice
# Add the OpenCG Lattice to the global collection of all OpenCG Lattices
OPENCG_LATTICES[lattice_id] = opencg_lattice
return openmc_lattice
def get_opencg_geometry(openmc_geometry):
if not isinstance(openmc_geometry, openmc.Geometry):
msg = 'Unable to get OpenCG geometry from {0} which is ' \
'not an OpenMC Geometry object'.format(openmc_geometry)
raise ValueError(msg)
# Clear dictionaries and auto-generated IDs
OPENMC_SURFACES.clear()
OPENCG_SURFACES.clear()
OPENMC_CELLS.clear()
OPENCG_CELLS.clear()
OPENMC_UNIVERSES.clear()
OPENCG_UNIVERSES.clear()
OPENMC_LATTICES.clear()
OPENCG_LATTICES.clear()
openmc_root_universe = openmc_geometry._root_universe
opencg_root_universe = get_opencg_universe(openmc_root_universe)
opencg_geometry = opencg.Geometry()
opencg_geometry.setRootUniverse(opencg_root_universe)
opencg_geometry.initializeCellOffsets()
return opencg_geometry
def get_openmc_geometry(opencg_geometry):
if not isinstance(opencg_geometry, opencg.Geometry):
msg = 'Unable to get OpenMC geometry from {0} which is ' \
'not an OpenCG Geometry object'.format(opencg_geometry)
raise ValueError(msg)
# Deep copy the goemetry since it may be modified to make all Surfaces
# compatible with OpenMC's specifications
opencg_geometry.assignAutoIds()
opencg_geometry = copy.deepcopy(opencg_geometry)
# Update Cell bounding boxes in Geometry
opencg_geometry.updateBoundingBoxes()
# Clear dictionaries and auto-generated ID
OPENMC_SURFACES.clear()
OPENCG_SURFACES.clear()
OPENMC_CELLS.clear()
OPENCG_CELLS.clear()
OPENMC_UNIVERSES.clear()
OPENCG_UNIVERSES.clear()
OPENMC_LATTICES.clear()
OPENCG_LATTICES.clear()
# Make the entire geometry "compatible" before assigning auto IDs
universes = opencg_geometry.getAllUniverses()
for universe_id, universe in universes.items():
if not isinstance(universe, opencg.Lattice):
make_opencg_cells_compatible(universe)
opencg_geometry.assignAutoIds()
opencg_root_universe = opencg_geometry._root_universe
openmc_root_universe = get_openmc_universe(opencg_root_universe)
openmc_geometry = openmc.Geometry()
openmc_geometry.root_universe = openmc_root_universe
return openmc_geometry

517
src/utils/openmc/plots.py Normal file
View file

@ -0,0 +1,517 @@
from xml.etree import ElementTree as ET
import numpy as np
from openmc.checkvalue import *
from openmc.clean_xml import *
# A static variable for auto-generated Plot IDs
AUTO_PLOT_ID = 10000
def reset_auto_plot_id():
global AUTO_PLOT_ID
AUTO_PLOT_ID = 10000
BASES = ['xy', 'xz', 'yz']
class Plot(object):
def __init__(self, plot_id=None, name=''):
# Initialize Plot class attributes
self.id = plot_id
self.name = name
self._width = [4.0, 4.0]
self._pixels = [1000, 1000]
self._origin = [0., 0., 0.]
self._filename = 'plot'
self._color = 'cell'
self._type = 'slice'
self._basis = 'xy'
self._background = [0, 0, 0]
self._mask_components = None
self._mask_background = None
self._col_spec = None
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@property
def width(self):
return self._width
@property
def pixels(self):
return self._pixels
@property
def origin(self):
return self._origin
@property
def filename(self):
return self._filename
@property
def color(self):
return self._color
@property
def type(self):
return self._type
@property
def basis(self):
return self._basis
@property
def background(self):
return self._background
@property
def mask_componenets(self):
return self._mask_components
@property
def mask_background(self):
return self._mask_background
@property
def col_spec(self):
return self._col_spec
@id.setter
def id(self, plot_id):
if plot_id is None:
global AUTO_PLOT_ID
self._id = AUTO_PLOT_ID
AUTO_PLOT_ID += 1
# Check that the ID is an integer and wasn't already used
elif not is_integer(plot_id):
msg = 'Unable to set a non-integer Plot ID {0}'.format(plot_id)
raise ValueError(msg)
elif plot_id < 0:
msg = 'Unable to set Plot ID to {0} since it must be a ' \
'non-negative integer'.format(plot_id)
raise ValueError(msg)
else:
self._id = plot_id
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Plot ID={0} with a non-string ' \
'value {1}'.format(self._id, name)
raise ValueError(msg)
else:
self._name = name
@width.setter
def width(self, width):
if not isinstance(width, (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with width {1} which is not ' \
'a Python tuple/list or NumPy array'.format(self._id, width)
raise ValueError(msg)
elif len(width) != 2 and len(width) != 3:
msg = 'Unable to create Plot ID={0} with width {1} since only 2D ' \
'and 3D plots are supported'.format(self._id, width)
raise ValueError(msg)
for dim in width:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to create Plot ID={0} with width {1} since ' \
'each element must be a floating point value or ' \
'integer'.format(self._id, width)
raise ValueError(msg)
self._width = width
@origin.setter
def origin(self, origin):
if not isinstance(origin, (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with origin {1} which is not ' \
'a Python tuple/list or NumPy array'.format(self._id, origin)
raise ValueError(msg)
elif len(origin) != 3:
msg = 'Unable to create Plot ID={0} with origin {1} since only ' \
'a 3D coordinate must be input'.format(self._id, origin)
raise ValueError(msg)
for dim in origin:
if not is_integer(dim) and not is_float(dim):
msg = 'Unable to create Plot ID={0} with origin {1} since ' \
'each element must be a floating point value or ' \
'integer'.format(self._id, origin)
raise ValueError(msg)
self._origin = origin
@pixels.setter
def pixels(self, pixels):
if not isinstance(pixels, (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with pixels {1} which is not ' \
'a Python tuple/list or NumPy array'.format(self._id, pixels)
raise ValueError(msg)
elif len(pixels) != 2 and len(pixels) != 3:
msg = 'Unable to create Plot ID={0} with pixels {1} since ' \
'only 2D and 3D plots are supported'.format(self._id, pixels)
raise ValueError(msg)
for dim in pixels:
if not is_integer(dim):
msg = 'Unable to create Plot ID={0} with pixel value {1} ' \
'which is not an integer'.format(self._id, dim)
raise ValueError(msg)
elif dim < 0:
msg = 'Unable to create Plot ID={0} with pixel value {1} ' \
'which is less than 0'.format(self._id, dim)
raise ValueError(msg)
self._pixels = pixels
@filename.setter
def filename(self, filename):
if not is_string(filename):
msg = 'Unable to create Plot ID={0} with filename {1} which is ' \
'not a string'.format(self._id, filename)
raise ValueError(msg)
self._filename = filename
@color.setter
def color(self, color):
if not is_string(color):
msg = 'Unable to create Plot ID={0} with color {1} which is not ' \
'a string'.format(self._id, color)
raise ValueError(msg)
elif not color in ['cell', 'mat']:
msg = 'Unable to create Plot ID={0} with color {1} which is not ' \
'a cell or mat'.format(self._id, color)
raise ValueError(msg)
self._color = color
@type.setter
def type(self, type):
if not is_string(type):
msg = 'Unable to create Plot ID={0} with type {1} which is not ' \
'a string'.format(self._id, type)
raise ValueError(msg)
elif not type in ['slice', 'voxel']:
msg = 'Unable to create Plot ID={0} with type {1} which is not ' \
'slice or voxel'.format(self._id, type)
raise ValueError(msg)
self._type = type
@basis.setter
def basis(self, basis):
if not is_string(basis):
msg = 'Unable to create Plot ID={0} with basis {1} which is not ' \
'a string'.format(self._id, basis)
raise ValueError(msg)
elif not basis in ['xy', 'xz', 'yz']:
msg = 'Unable to create Plot ID={0} with basis {1} which is not ' \
'xy, xz, or yz'.format(self._id, basis)
raise ValueError(msg)
self._basis = basis
@background.setter
def background(self, background):
if not isinstance(background, (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with background {1} ' \
'which is not a Python tuple/list or NumPy ' \
'array'.format(self._id, background)
raise ValueError(msg)
elif len(background) != 3:
msg = 'Unable to create Plot ID={0} with background {1} ' \
'which is not 3 integer RGB ' \
'values'.format(self._id, background)
raise ValueError(msg)
for rgb in background:
if not is_integer(rgb):
msg = 'Unable to create Plot ID={0} with background RGB ' \
'value {1} which is not an integer'.format(self._id, rgb)
raise ValueError(msg)
elif rgb < 0 or rgb > 255:
msg = 'Unable to create Plot ID={0} with background RGB value ' \
'{1} which is not between 0 and 255'.format(self._id, rgb)
raise ValueError(msg)
self._background = background
@col_spec.setter
def col_spec(self, col_spec):
if not isinstance(col_spec, dict):
msg= 'Unable to create Plot ID={0} with col_spec parameter {1} ' \
'which is not a Python dictionary of IDs to ' \
'pixels'.format(self._id, col_spec)
raise ValueError(msg)
for key in col_spec:
if not is_integer(key):
msg = 'Unable to create Plot ID={0} with col_spec ID {1} ' \
'which is not an integer'.format(self._id, key)
raise ValueError(msg)
elif key < 0:
msg = 'Unable to create Plot ID={0} with col_spec ID {1} ' \
'which is less than 0'.format(self._id, key)
raise ValueError(msg)
elif not isinstance(col_spec[key], (tuple, list, np.ndarray)):
msg = 'Unable to create Plot ID={0} with col_spec RGB ' \
'values {1} which is not a Python tuple/list or NumPy ' \
'array'.format(self._id, col_spec[key])
raise ValueError(msg)
elif len(col_spec[key]) != 3:
msg = 'Unable to create Plot ID={0} with col_spec RGB ' \
'values of length {1} since 3 values must be ' \
'input'.format(self._id, len(col_spec[key]))
raise ValueError(msg)
self._col_spec = col_spec
@mask_componenets.setter
def mask_components(self, mask_components):
if not isinstance(mask_components, (list, tuple, np.ndarray)):
msg = 'Unable to create Plot ID={0} with mask components {1} ' \
'which is not a Python tuple/list or NumPy ' \
'array'.format(self._id, mask_components)
raise ValueError(msg)
for component in mask_components:
if not is_integer(component):
msg = 'Unable to create Plot ID={0} with mask component {1} ' \
'which is not an integer'.format(self._id, component)
raise ValueError(msg)
elif component < 0:
msg = 'Unable to create Plot ID={0} with mask component {1} ' \
'which is less than 0'.format(self._id, component)
raise ValueError(msg)
self._mask_components = mask_components
@mask_background.setter
def mask_background(self, mask_background):
if not isinstance(mask_background, (list, tuple, np.ndarray)):
msg = 'Unable to create Plot ID={0} with mask background {1} ' \
'which is not a Python tuple/list or NumPy ' \
'array'.format(self._id, mask_background)
raise ValueError(msg)
elif len(mask_background) != 3 and len(mask_background) != 0:
msg = 'Unable to create Plot ID={0} with mask background ' \
'{1} since 3 RGB values must be ' \
'input'.format(self._id, mask_background)
raise ValueError(msg)
for rgb in mask_background:
if not is_integer(rgb):
msg = 'Unable to create Plot ID={0} with mask background RGB ' \
'value {1} which is not an integer'.format(self._id, rgb)
raise ValueError(msg)
elif rgb < 0 or rgb > 255:
msg = 'Unable to create Plot ID={0} with mask bacground ' \
'RGB value {1} which is not between 0 and ' \
'255'.format(self._id, rgb)
raise ValueError(msg)
self._mask_background = mask_background
def __repr__(self):
string = 'Plot\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
string += '{0: <16}{1}{2}\n'.format('\tFilename', '=\t', self._filename)
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type)
string += '{0: <16}{1}{2}\n'.format('\tBasis', '=\t', self._basis)
string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self._width)
string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._origin)
string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._origin)
string += '{0: <16}{1}{2}\n'.format('\tColor', '=\t', self._color)
string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t',
self._mask_components)
string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t',
self._mask_background)
string += '{0: <16}{1}{2}\n'.format('\tCol Spec', '=\t', self._col_spec)
return string
def get_plot_xml(self):
element = ET.Element("plot")
element.set("id", str(self._id))
element.set("filename", self._filename)
element.set("color", self._color)
element.set("type", self._type)
if self._type is 'slice':
element.set("basis", self._basis)
subelement = ET.SubElement(element, "origin")
text = ''
for coord in self._origin:
text += str(coord) + ' '
subelement.text = text.rstrip(' ')
subelement = ET.SubElement(element, "width")
text = ''
for dim in self._width:
text += str(dim) + ' '
subelement.text = text.rstrip(' ')
subelement = ET.SubElement(element, "pixels")
text = ''
for dim in self._pixels:
text += str(dim) + ' '
subelement.text = text.rstrip(' ')
if not self._mask_background is None:
subelement = ET.SubElement(element, "background")
text = ''
for rgb in self._background:
text += str(rgb) + ' '
subelement.text = text.rstrip(' ')
if not self._col_spec is None:
for key in self._col_spec:
subelement = ET.SubElement(element, "col_spec")
subelement.set("id", '{0}'.format(key))
value = self._col_spec[key]
subelement.set("rgb",'{0} {1} {2}'.format(value[0],
value[1], value[2]))
if not self._mask_components is None:
subelement = ET.SubElement(element, "mask")
text = ''
for id in self._mask_components:
text += str(id) + ' '
subelement.set("components", text.rstrip(' '))
rgb = self._mask_background
subelement.set("background", '{0} {1} {2}'.format(rgb[0],
rgb[1], rgb[2]))
return element
class PlotsFile(object):
def __init__(self):
# Initialize PlotsFile class attributes
self._plots = []
self._plots_file = ET.Element("plots")
def add_plot(self, plot):
if not isinstance(plot, Plot):
msg = 'Unable to add a non-Plot {0} to the PlotsFile'.format(plot)
raise ValueError(msg)
self._plots.append(plot)
def remove_plot(self, plot):
self._plots.remove(plot)
def create_plot_subelements(self):
for plot in self._plots:
xml_element = plot.get_plot_xml()
if len(plot._name) > 0:
self._plots_file.append(ET.Comment(plot._name))
self._plots_file.append(xml_element)
def export_to_xml(self):
self.create_plot_subelements()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._plots_file)
# Write the XML Tree to the plots.xml file
tree = ET.ElementTree(self._plots_file)
tree.write("plots.xml", xml_declaration=True,
encoding='utf-8', method="xml")

1514
src/utils/openmc/settings.py Normal file

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,753 @@
import copy
import struct
import numpy as np
import scipy.stats
import openmc
from openmc.constants import *
class SourceSite(object):
def __init__(self):
self._weight = None
self._xyz = None
self._uvw = None
self._E = None
def __repr__(self):
string = 'SourceSite\n'
string += '{0: <16}{1}{2}\n'.format('\tweight', '=\t', self._weight)
string += '{0: <16}{1}{2}\n'.format('\tE', '=\t', self._E)
string += '{0: <16}{1}{2}\n'.format('\t(x,y,z)', '=\t', self._xyz)
string += '{0: <16}{1}{2}\n'.format('\t(u,v,w)', '=\t', self._uvw)
return string
@property
def weight(self):
return self._weight
@property
def xyz(self):
return self._xyz
@property
def uvw(self):
return self._uvw
@property
def E(self):
return self._E
class StatePoint(object):
def __init__(self, filename):
if filename.endswith('.h5'):
import h5py
self._f = h5py.File(filename, 'r')
self._hdf5 = True
else:
self._f = open(filename, 'rb')
self._hdf5 = False
# Set flags for what data has been read
self._results = False
self._source = False
self._with_summary = False
# Read all metadata
self._read_metadata()
# Read information about tally meshes
self._read_meshes()
# Read tally metadata
self._read_tallies()
def close(self):
self._f.close()
@property
def k_combined(self):
return self._k_combined
@property
def n_particles(self):
return self._n_particles
@property
def n_batches(self):
return self._n_batches
@property
def current_batch(self):
return self._current_batch
@property
def results(self):
return self._results
@property
def source(self):
return self._source
@property
def with_summary(self):
return self._with_summary
@property
def tallies(self):
return self._tallies
def _read_metadata(self):
# Read filetype
self._filetype = self._get_int(path='filetype')[0]
# Read statepoint revision
self._revision = self._get_int(path='revision')[0]
if self._revision != 13:
raise Exception('Statepoint Revision is not consistent.')
# Read OpenMC version
if self._hdf5:
self._version = [self._get_int(path='version_major')[0],
self._get_int(path='version_minor')[0],
self._get_int(path='version_release')[0]]
else:
self._version = self._get_int(3)
# Read date and time
self._date_and_time = self._get_string(19, path='date_and_time')
# Read path
self._path = self._get_string(255, path='path').strip()
# Read random number seed
self._seed = self._get_long(path='seed')[0]
# Read run information
self._run_mode = self._get_int(path='run_mode')[0]
self._n_particles = self._get_long(path='n_particles')[0]
self._n_batches = self._get_int(path='n_batches')[0]
# Read current batch
self._current_batch = self._get_int(path='current_batch')[0]
# Read whether or not the source site distribution is present
self._source_present = self._get_int(path='source_present')[0]
# Read criticality information
if self._run_mode == 2:
self._read_criticality()
def _read_criticality(self):
# Read criticality information
if self._run_mode == 2:
self._n_inactive = self._get_int(path='n_inactive')[0]
self._gen_per_batch = self._get_int(path='gen_per_batch')[0]
self._k_batch = self._get_double(
self._current_batch*self._gen_per_batch,
path='k_generation')
self._entropy = self._get_double(
self._current_batch*self._gen_per_batch, path='entropy')
self._k_col_abs = self._get_double(path='k_col_abs')[0]
self._k_col_tra = self._get_double(path='k_col_tra')[0]
self._k_abs_tra = self._get_double(path='k_abs_tra')[0]
self._k_combined = self._get_double(2, path='k_combined')
# Read CMFD information (if used)
self._read_cmfd()
def _read_cmfd(self):
base = 'cmfd'
# Read CMFD information
self._cmfd_on = self._get_int(path='cmfd_on')[0]
if self._cmfd_on == 1:
self._cmfd_indices = self._get_int(4, path='{0}/indices'.format(base))
self._k_cmfd = self._get_double(self._current_batch,
path='{0}/k_cmfd'.format(base))
self._cmfd_src = self._get_double_array(np.product(self._cmfd_indices),
path='{0}/cmfd_src'.format(base))
self._cmfd_src = np.reshape(self._cmfd_src, tuple(self._cmfd_indices),
order='F')
self._cmfd_entropy = self._get_double(self._current_batch,
path='{0}/cmfd_entropy'.format(base))
self._cmfd_balance = self._get_double(self._current_batch,
path='{0}/cmfd_balance'.format(base))
self._cmfd_dominance = self._get_double(self._current_batch,
path='{0}/cmfd_dominance'.format(base))
self._cmfd_srccmp = self._get_double(self._current_batch,
path='{0}/cmfd_srccmp'.format(base))
def _read_meshes(self):
# Initialize dictionaries for the Meshes
# Keys - Mesh IDs
# Values - Mesh objects
self._meshes = {}
# Read the number of Meshes
self._n_meshes = self._get_int(path='tallies/meshes/n_meshes')[0]
# Read a list of the IDs for each Mesh
if self._n_meshes > 0:
# OpenMC Mesh IDs (redefined internally from user definitions)
self._mesh_ids = self._get_int(self._n_meshes,
path='tallies/meshes/ids')
# User-defined Mesh IDs
self._mesh_keys = self._get_int(self._n_meshes,
path='tallies/meshes/keys')
else:
self._mesh_keys = []
self._mesh_ids = []
# Build dictionary of Meshes
base = 'tallies/meshes/mesh '
# Iterate over all Meshes
for mesh_key in self._mesh_keys:
# Read the user-specified Mesh ID and type
mesh_id = self._get_int(path='{0}{1}/id'.format(base, mesh_key))[0]
mesh_type = self._get_int(path='{0}{1}/type'.format(base, mesh_key))[0]
# Get the Mesh dimension
n_dimension = self._get_int(
path='{0}{1}/n_dimension'.format(base, mesh_key))[0]
# Read the mesh dimensions, lower-left coordinates,
# upper-right coordinates, and width of each mesh cell
dimension = self._get_int(
n_dimension, path='{0}{1}/dimension'.format(base, mesh_key))
lower_left = self._get_double(
n_dimension, path='{0}{1}/lower_left'.format(base, mesh_key))
upper_right = self._get_double(
n_dimension, path='{0}{1}/upper_right'.format(base, mesh_key))
width = self._get_double(
n_dimension, path='{0}{1}/width'.format(base, mesh_key))
# Create the Mesh and assign properties to it
mesh = openmc.Mesh(mesh_id)
mesh.dimension = dimension
mesh.width = width
mesh.lower_left = lower_left
mesh.upper_right = upper_right
#FIXME: Set the mesh type to 'rectangular' by default
mesh.type = 'rectangular'
# Add mesh to the global dictionary of all Meshes
self._meshes[mesh_id] = mesh
def _read_tallies(self):
# Initialize dictionaries for the Tallies
# Keys - Tally IDs
# Values - Tally objects
self._tallies = {}
# Read the number of tallies
self._n_tallies = self._get_int(path='/tallies/n_tallies')[0]
# Read a list of the IDs for each Tally
if self._n_tallies > 0:
# OpenMC Tally IDs (redefined internally from user definitions)
self._tally_ids = self._get_int(
self._n_tallies, path='tallies/ids')
# User-defined Tally IDs
self._tally_keys = self._get_int(
self._n_tallies, path='tallies/keys')
else:
self._tally_keys = []
self._tally_ids = []
base = 'tallies/tally '
# Iterate over all Tallies
for tally_key in self._tally_keys:
# Read user-specified Tally label (if specified)
label_size = self._get_int(
path='{0}{1}/label_size'.format(base, tally_key))[0]
if label_size > 0:
label = self._get_string(
label_size, path='{0}{1}/label'.format(base, tally_key))
# Remove leading and trailing characters from string label
label = label.lstrip('[\'')
label = label.rstrip('\']')
# Read integer Tally estimator type code (analog or tracklength)
estimator_type = self._get_int(
path='{0}{1}/estimator'.format(base, tally_key))[0]
# Read the Tally size specifications
n_realizations = self._get_int(
path='{0}{1}/n_realizations'.format(base, tally_key))[0]
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_key, label)
tally.estimator = ESTIMATOR_TYPES[estimator_type]
tally.num_realizations = n_realizations
# Read the number of Filters
n_filters = self._get_int(
path='{0}{1}/n_filters'.format(base, tally_key))[0]
subbase = '{0}{1}/filter '.format(base, tally_key)
# Initialize the stride
stride = 1
# Initialize all Filters
for j in range(1, n_filters+1):
# Read the integer Filter type code
filter_type = self._get_int(
path='{0}{1}/type'.format(subbase, j))[0]
n_bins = self._get_int(
path='{0}{1}/n_bins'.format(subbase, j))[0]
if n_bins <= 0:
msg = 'Unable to create Filter {0} for Tally ID={2} ' \
'since no bins were specified'.format(j, tally_key)
raise ValueError(msg)
# Read the bin values
if FILTER_TYPES[filter_type] in ['energy', 'energyout']:
bins = self._get_double(
n_bins+1, path='{0}{1}/bins'.format(subbase, j))
# FIXME
elif FILTER_TYPES[filter_type] in ['mesh', 'distribcell']:
bins = self._get_int(
path='{0}{1}/bins'.format(subbase, j))[0]
else:
bins = self._get_int(
n_bins, path='{0}{1}/bins'.format(subbase, j))
# Create Filter object
filter = openmc.Filter(FILTER_TYPES[filter_type], bins)
filter.stride = stride
filter.num_bins = n_bins
if FILTER_TYPES[filter_type] == 'mesh':
filter.mesh = self._meshes[bins]
# Add Filter to the Tally
tally.add_filter(filter)
# Update the stride for the next Filter
stride *= n_bins
# Read Nuclide bins
n_nuclides = self._get_int(
path='{0}{1}/n_nuclides'.format(base, tally_key))[0]
nuclide_zaids = self._get_int(
n_nuclides, path='{0}{1}/nuclides'.format(base, tally_key))
# Add all Nuclides to the Tally
for nuclide_zaid in nuclide_zaids:
tally.add_nuclide(nuclide_zaid)
# Read score bins
n_score_bins = self._get_int(
path='{0}{1}/n_score_bins'.format(base, tally_key))[0]
tally.num_score_bins = n_score_bins
scores = [SCORE_TYPES[j] for j in self._get_int(
n_score_bins, path='{0}{1}/score_bins'.format(base, tally_key))]
n_user_scores = self._get_int(
path='{0}{1}/n_user_score_bins'.format(base, tally_key))[0]
# Read scattering moment order strings (e.g., P3, Y-1,2, etc.)
moments = []
subbase = '{0}{1}/moments/'.format(base, tally_key)
# Extract the moment order string for each score
for k in range(len(scores)):
moment = self._get_string(8,
path='{0}order{1}'.format(subbase, k+1))
moment = moment.lstrip('[\'')
moment = moment.rstrip('\']')
# Remove extra whitespace
moment.replace(" ", "")
moments.append(moment)
# Add the scores to the Tally
for j, score in enumerate(scores):
# If this is a scattering moment, insert the scattering order
if '-n' in score:
score = score.replace('-n', '-' + str(moments[j]))
elif '-pn' in score:
score = score.replace('-pn', '-' + str(moments[j]))
elif '-yn' in score:
score = score.replace('-yn', '-' + str(moments[j]))
tally.add_score(score)
# Add Tally to the global dictionary of all Tallies
self._tallies[tally_key] = tally
def read_results(self):
# Number of realizations for global Tallies
self._n_realizations = self._get_int(path='n_realizations')[0]
# Read global Tallies
n_global_tallies = self._get_int(path='n_global_tallies')[0]
if self._hdf5:
data = self._f['global_tallies'].value
self._global_tallies = np.column_stack((data['sum'], data['sum_sq']))
else:
self._global_tallies = np.array(self._get_double(2*n_global_tallies))
self._global_tallies.shape = (n_global_tallies, 2)
# Flag indicating if Tallies are present
self._tallies_present = self._get_int(path='tallies/tallies_present')[0]
base = 'tallies/tally '
# Read Tally results
if self._tallies_present:
# Iterate over and extract the results for all Tallies
for tally_key in self._tally_keys:
# Get this Tally
tally = self._tallies[tally_key]
# Compute the total number of bins for this Tally
num_tot_bins = tally.num_bins
# Extract Tally data from the file
if self._hdf5:
data = self._f['{0}{1}/results'.format(base, tally_key)].value
sum = data['sum']
sum_sq = data['sum_sq']
else:
results = np.array(self._get_double(2*num_tot_bins))
sum = results[0::2]
sum_sq = results[1::2]
# Define a routine to convert 0 to 1
nonzero = lambda val: 1 if not val else val
# Reshape the results arrays
new_shape = (nonzero(tally.num_filter_bins),
nonzero(tally.num_nuclides),
nonzero(tally.num_score_bins))
sum = np.reshape(sum, new_shape)
sum_sq = np.reshape(sum_sq, new_shape)
# Set the data for this Tally
tally.set_results(sum=sum, sum_sq=sum_sq)
# Indicate that Tally results have been read
self._results = True
def read_source(self):
# Check whether Tally results have been read
if not self._results:
self.read_results()
# Check if source bank is in statepoint
if not self._source_present:
print('Unable to read source since it is not in statepoint file')
return
# Initialize a NumPy array for the source sites
self._source = np.empty(self._n_particles, dtype=SourceSite)
# For HDF5 state points, copy entire bank
if self._hdf5:
source_sites = self._f['source_bank'].value
# Initialize SourceSite object for each particle
for i in range(self._n_particles):
# Initialize new source site
site = SourceSite()
# Read position, angle, and energy
if self._hdf5:
site._weight, site._xyz, site._uvw, site._E = source_sites[i]
else:
site._weight = self._get_double()[0]
site._xyz = self._get_double(3)
site._uvw = self._get_double(3)
site._E = self._get_double()[0]
# Store the source site in the NumPy array
self._source[i] = site
def compute_ci(self, confidence=0.95):
"""Computes confidence intervals for each Tally bin."""
# Determine significance level and percentile for two-sided CI
alpha = 1 - confidence
percentile = 1 - alpha/2
# Calculate t-value
t_value = scipy.stats.t.ppf(percentile, self._n_realizations - 1)
self.compute_stdev(t_value)
def compute_stdev(self, t_value=1.0):
"""
Computes the sample mean and the standard deviation of the mean
for each Tally bin.
"""
# Determine number of realizations
n = self._n_realizations
# Calculate the standard deviation for each global tally
for i in range(len(self._global_tallies)):
# Get sum and sum of squares
s, s2 = self._global_tallies[i]
# Calculate sample mean and replace value
s /= n
self._global_tallies[i, 0] = s
# Calculate standard deviation
if s != 0.0:
self._global_tallies[i, 1] = t_value * np.sqrt((s2 / n - s**2) / (n-1))
# Calculate sample mean and standard deviation for user-defined Tallies
for tally_id, tally in self._tallies.items():
tally.compute_std_dev(t_value)
def get_tally(self, score, filters, nuclides,
label='', estimator='tracklength'):
"""Finds and returns a Tally object with certain properties.
Parameters
----------
score : str
The score string
filters : list
A list of Filter objects
nuclides : list
A list of Nuclide objects
label : str
The label specified for the Tally (default is '')
estimator: str
The type of estimator ('tracklength' (default) or 'analog')
"""
# Loop over the domain-to-tallies mapping to find the Tally
tally = None
# Iterate over all tallies to find the appropriate one
for tally_id, test_tally in self._tallies.items():
# Determine if the queried Tally label is the same as this Tally
if not label == test_tally._label:
continue
# Determine if the queried Tally estimator is the same as this Tally
if not estimator == test_tally._estimator:
continue
# Determine if the queried Tally scores are the same as this Tally
if not score in test_tally._scores:
continue
# Determine if queried Tally filters is the same length as this Tally
if len(filters) != len(test_tally._filters):
continue
# Determine if the queried Tally filters are the same as this Tally
contains_filters = True
# Iterate over the filters requested by the user
for filter in filters:
if not filter in test_tally._filters:
contains_filters = False
break
# Determine if the queried Nuclide is in this Tally
contains_nuclides = True
# Iterate over the Nuclides requested by the user
for nuclide in nuclides:
if not nuclide in test_tally._nuclides:
contains_nuclides = False
break
# If the Tally contained all Filters and Nuclides, return the Tally
if contains_filters and contains_nuclides:
tally = test_tally
break
# If we did not find the Tally, return an error message
if tally is None:
raise LookupError('Unable to get Tally')
return tally
def get_tally_id(self, score, filters, label='', estimator='tracklength'):
"""Retrieve the Tally ID for a given list of filters and score(s).
Parameters
----------
score : str
The score string
filters : list
A list of Filter objects
label : str
The label specified for the Tally (default is '')
estimator: str
The type of estimator ('tracklength' (default) or 'analog')
"""
tally = self.get_tally(score, filters, label, estimator)
return tally._id
def link_with_summary(self, summary):
if not isinstance(summary, openmc.summary.Summary):
msg = 'Unable to link statepoint with {0} which ' \
'is not a Summary object'.format(summary)
raise ValueError(msg)
for tally_id, tally in self._tallies.items():
nuclide_zaids = copy.deepcopy(tally._nuclides)
for nuclide_zaid in nuclide_zaids:
tally.remove_nuclide(nuclide_zaid)
if nuclide_zaid == -1:
tally.add_nuclide(openmc.Nuclide('total'))
else:
tally.add_nuclide(summary.nuclides[nuclide_zaid])
for filter in tally._filters:
if filter._type == 'surface':
surface_ids = []
for bin in filter._bins:
surface_ids.append(summary.surfaces[bin]._id)
filter.bins = surface_ids
if filter._type in ['cell', 'distribcell']:
distribcell_ids = []
for bin in filter._bins:
distribcell_ids.append(summary.cells[bin]._id)
filter.bins = distribcell_ids
if filter._type == 'universe':
universe_ids = []
for bin in filter._bins:
universe_ids.append(summary.universes[bin]._id)
filter.bins = universe_ids
if filter._type == 'material':
material_ids = []
for bin in filter._bins:
material_ids.append(summary.materials[bin]._id)
filter.bins = material_ids
self._with_summary = True
def _get_data(self, n, typeCode, size):
return list(struct.unpack('={0}{1}'.format(n,typeCode),
self._f.read(n*size)))
def _get_int(self, n=1, path=None):
if self._hdf5:
return [int(v) for v in self._f[path].value]
else:
return [int(v) for v in self._get_data(n, 'i', 4)]
def _get_long(self, n=1, path=None):
if self._hdf5:
return [long(v) for v in self._f[path].value]
else:
return [long(v) for v in self._get_data(n, 'q', 8)]
def _get_float(self, n=1, path=None):
if self._hdf5:
return [float(v) for v in self._f[path].value]
else:
return [float(v) for v in self._get_data(n, 'f', 4)]
def _get_double(self, n=1, path=None):
if self._hdf5:
return [float(v) for v in self._f[path].value]
else:
return [float(v) for v in self._get_data(n, 'd', 8)]
def _get_double_array(self, n=1, path=None):
if self._hdf5:
return self._f[path].value
else:
return self._get_data(n, 'd', 8)
def _get_string(self, n=1, path=None):
if self._hdf5:
return str(self._f[path].value)
else:
return str(self._get_data(n, 's', 1)[0])

518
src/utils/openmc/summary.py Normal file
View file

@ -0,0 +1,518 @@
import numpy as np
import openmc
from openmc.opencg_compatible import get_opencg_geometry
try:
import h5py
except ImportError:
msg = 'Unable to import h5py which is necessary for openmc.summary'
raise ImportError(msg)
class Summary(object):
def __init__(self, filename):
openmc.reset_auto_ids()
if not filename.endswith(('.h5', '.hdf5')):
msg = 'Unable to open "{0}" which is not an HDF5 summary file'
raise ValueError(msg)
self._f = h5py.File(filename, 'r')
self.openmc_geometry = None
self.opencg_geometry = None
self._read_metadata()
self._read_geometry()
def _read_metadata(self):
# Read OpenMC version
self.version = [self._f['version_major'][0],
self._f['version_minor'][0],
self._f['version_release'][0]]
# Read date and time
self.date_and_time = self._f['date_and_time'][...]
self.n_batches = self._f['n_batches'][0]
self.n_particles = self._f['n_particles'][0]
self.n_active = self._f['n_active'][0]
self.n_inactive = self._f['n_inactive'][0]
self.gen_per_batch = self._f['gen_per_batch'][0]
self.n_procs = self._f['n_procs'][0]
def _read_geometry(self):
# Read in and initialize the Materials and Geometry
self._read_nuclides()
self._read_materials()
self._read_surfaces()
self._read_cells()
self._read_universes()
self._read_lattices()
self._finalize_geometry()
def _read_nuclides(self):
self.n_nuclides = self._f['nuclides/n_nuclides'][0]
# Initialize dictionary for each Nuclide
# Keys - Nuclide ZAIDs
# Values - Nuclide objects
self.nuclides = {}
for key in self._f['nuclides'].keys():
if key == 'n_nuclides':
continue
index = self._f['nuclides'][key]['index'][0]
alias = self._f['nuclides'][key]['alias'][0]
zaid = self._f['nuclides'][key]['zaid'][0]
# Read the Nuclide's name (e.g., 'H-1' or 'U-235')
name = alias.split('.')[0]
# Read the Nuclide's cross-section identifier (e.g., '70c')
xs = alias.split('.')[1]
# Initialize this Nuclide and add to global dictionary of Nuclides
if 'nat' in name:
self.nuclides[zaid] = openmc.Element(name=name, xs=xs)
else:
self.nuclides[zaid] = openmc.Nuclide(name=name, xs=xs)
self.nuclides[zaid].zaid = zaid
def _read_materials(self):
self.n_materials = self._f['materials/n_materials'][0]
# Initialize dictionary for each Material
# Keys - Material keys
# Values - Material objects
self.materials = {}
for key in self._f['materials'].keys():
if key == 'n_materials':
continue
material_id = int(key.lstrip('material '))
index = self._f['materials'][key]['index'][0]
density = self._f['materials'][key]['atom_density'][0]
nuc_densities = self._f['materials'][key]['nuclide_densities'][...]
nuclides = self._f['materials'][key]['nuclides'][...]
n_sab = self._f['materials'][key]['n_sab'][0]
sab_names = []
sab_xs = []
# Read the names of the S(a,b) tables for this Material
for i in range(1, n_sab+1):
sab_table = self._f['materials'][key]['sab_tables'][str(i)][0]
# Read the cross-section identifiers for each S(a,b) table
sab_names.append(sab_table.split('.')[0])
sab_xs.append(sab_table.split('.')[1])
# Create the Material
material = openmc.Material(material_id=material_id)
# Set the Material's density to g/cm3 - this is what is used in OpenMC
material.set_density(density=density, units='g/cm3')
# Add all Nuclides to the Material
for i, zaid in enumerate(nuclides):
nuclide = self.get_nuclide_by_zaid(zaid)
density = nuc_densities[i]
if isinstance(nuclide, openmc.Nuclide):
material.add_nuclide(nuclide, percent=density, percent_type='ao')
elif isinstance(nuclide, openmc.Element):
material.add_element(nuclide, percent=density, percent_type='ao')
# Add S(a,b) table(s?) to the Material
for i in range(n_sab):
name = sab_names[i]
xs = sab_xs[i]
material.add_s_alpha_beta(name, xs)
# Add the Material to the global dictionary of all Materials
self.materials[index] = material
def _read_surfaces(self):
self.n_surfaces = self._f['geometry/n_surfaces'][0]
# Initialize dictionary for each Surface
# Keys - Surface keys
# Values - Surfacee objects
self.surfaces = {}
for key in self._f['geometry/surfaces'].keys():
if key == 'n_surfaces':
continue
surface_id = int(key.lstrip('surface '))
index = self._f['geometry/surfaces'][key]['index'][0]
surf_type = self._f['geometry/surfaces'][key]['type'][...][0]
bc = self._f['geometry/surfaces'][key]['boundary_condition'][...][0]
coeffs = self._f['geometry/surfaces'][key]['coefficients'][...]
# Create the Surface based on its type
if surf_type == 'X Plane':
x0 = coeffs[0]
surface = openmc.XPlane(surface_id, bc, x0)
elif surf_type == 'Y Plane':
y0 = coeffs[0]
surface = openmc.YPlane(surface_id, bc, y0)
elif surf_type == 'Z Plane':
z0 = coeffs[0]
surface = openmc.ZPlane(surface_id, bc, z0)
elif surf_type == 'Plane':
A = coeffs[0]
B = coeffs[1]
C = coeffs[2]
D = coeffs[3]
surface = openmc.Plane(surface_id, bc, A, B, C, D)
elif surf_type == 'X Cylinder':
y0 = coeffs[0]
z0 = coeffs[1]
R = coeffs[2]
surface = openmc.XCylinder(surface_id, bc, y0, z0, R)
elif surf_type == 'Y Cylinder':
x0 = coeffs[0]
z0 = coeffs[1]
R = coeffs[2]
surface = openmc.YCylinder(surface_id, bc, x0, z0, R)
elif surf_type == 'Z Cylinder':
x0 = coeffs[0]
y0 = coeffs[1]
R = coeffs[2]
surface = openmc.ZCylinder(surface_id, bc, x0, y0, R)
elif surf_type == 'Sphere':
x0 = coeffs[0]
y0 = coeffs[1]
z0 = coeffs[2]
R = coeffs[3]
surface = openmc.Sphere(surface_id, bc, x0, y0, z0, R)
elif surf_type in ['X Cone', 'Y Cone', 'Z Cone']:
x0 = coeffs[0]
y0 = coeffs[1]
z0 = coeffs[2]
R2 = coeffs[3]
if surf_type == 'X Cone':
surface = openmc.XCone(surface_id, bc, x0, y0, z0, R2)
if surf_type == 'Y Cone':
surface = openmc.YCone(surface_id, bc, x0, y0, z0, R2)
if surf_type == 'Z Cone':
surface = openmc.ZCone(surface_id, bc, x0, y0, z0, R2)
# Add Surface to global dictionary of all Surfaces
self.surfaces[index] = surface
def _read_cells(self):
self.n_cells = self._f['geometry/n_cells'][0]
# Initialize dictionary for each Cell
# Keys - Cell keys
# Values - Cell objects
self.cells = {}
# Initialize dictionary for each Cell's fill
# (e.g., Material, Universe or Lattice ID)
# This dictionary is used later to link the fills with
# the corresponding objects
# Keys - Cell keys
# Values - Filling Material, Universe or Lattice ID
self._cell_fills = {}
for key in self._f['geometry/cells'].keys():
if key == 'n_cells':
continue
cell_id = int(key.lstrip('cell '))
index = self._f['geometry/cells'][key]['index'][0]
fill_type = self._f['geometry/cells'][key]['fill_type'][...][0]
if fill_type == 'normal':
fill = self._f['geometry/cells'][key]['material'][0]
elif fill_type == 'universe':
fill = self._f['geometry/cells'][key]['fill'][0]
else:
fill = self._f['geometry/cells'][key]['lattice'][0]
if 'surfaces' in self._f['geometry/cells'][key].keys():
surfaces = self._f['geometry/cells'][key]['surfaces'][...]
else:
surfaces = []
# Create this Cell
cell = openmc.Cell(cell_id=cell_id)
if fill_type == 'universe':
translated = self._f['geometry/cells'][key]['translated'][0]
if translated:
translation = \
self._f['geometry/cells'][key]['translation'][...]
translation = np.asarray(translation, dtype=np.float64)
cell.translation = translation
rotated = self._f['geometry/cells'][key]['rotated'][0]
if rotated:
rotation = \
self._f['geometry/cells'][key]['rotation'][...]
rotation = np.asarray(rotation, dtype=np.int)
cell.rotation = rotation
# Store Cell fill information for after Universe/Lattice creation
self._cell_fills[index] = (fill_type, fill)
# Iterate over all Surfaces and add them to the Cell
for surface_halfspace in surfaces:
halfspace = np.sign(surface_halfspace)
surface_id = np.abs(surface_halfspace)
surface = self.surfaces[surface_id]
cell.add_surface(surface, halfspace)
# Add the Cell to the global dictionary of all Cells
self.cells[index] = cell
def _read_universes(self):
self.n_universes = self._f['geometry/n_universes'][0]
# Initialize dictionary for each Universe
# Keys - Universe keys
# Values - Universe objects
self.universes = {}
for key in self._f['geometry/universes'].keys():
if key == 'n_universes':
continue
universe_id = int(key.lstrip('universe '))
index = self._f['geometry/universes'][key]['index'][0]
cells = self._f['geometry/universes'][key]['cells'][...]
# Create this Universe
universe = openmc.Universe(universe_id=universe_id)
# Add each Cell to the Universe
for cell_id in cells:
cell = self.cells[cell_id]
universe.add_cell(cell)
# Add the Universe to the global list of Universes
self.universes[index] = universe
def _read_lattices(self):
self.n_lattices = self._f['geometry/n_lattices'][0]
# Initialize lattices for each Lattice
# Keys - Lattice keys
# Values - Lattice objects
self.lattices = {}
for key in self._f['geometry/lattices'].keys():
if key == 'n_lattices':
continue
lattice_id = int(key.lstrip('lattice '))
index = self._f['geometry/lattices'][key]['index'][0]
lattice_type = self._f['geometry/lattices'][key]['type'][...][0]
if lattice_type == 'rectangular':
dimension = self._f['geometry/lattices'][key]['dimension'][...]
lower_left = \
self._f['geometry/lattices'][key]['lower_left'][...]
pitch = self._f['geometry/lattices'][key]['pitch'][...]
outer = self._f['geometry/lattices'][key]['outer'][0]
universe_ids = \
self._f['geometry/lattices'][key]['universes'][...]
universe_ids = np.swapaxes(universe_ids, 0, 1)
universe_ids = np.swapaxes(universe_ids, 1, 2)
# Create the Lattice
lattice = openmc.RectLattice(lattice_id=lattice_id)
lattice.dimension = tuple(dimension)
lattice.lower_left = lower_left
lattice.pitch = pitch
# If the Universe specified outer the Lattice is not void (-22)
if outer != -22:
lattice.outer = self.universes[outer]
# Build array of Universe pointers for the Lattice
universes = \
np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe)
for x in range(universe_ids.shape[0]):
for y in range(universe_ids.shape[1]):
for z in range(universe_ids.shape[2]):
universes[x,y,z] = \
self.get_universe_by_id(universe_ids[x,y,z])
# Transpose, reverse y-dimension for appropriate ordering
shape = universes.shape
universes = np.transpose(universes, (1,0,2))
universes.shape = shape
universes = universes[:,::-1,:]
lattice.universes = universes
# Add the Lattice to the global dictionary of all Lattices
self.lattices[index] = lattice
if lattice_type == 'hexagonal':
n_rings = self._f['geometry/latties'][key]['n_rings'][0]
n_axial = self._f['geometry/latties'][key]['n_axial'][0]
center = self._f['geometry/lattices'][key]['center'][...]
pitch = self._f['geometry/lattices'][key]['pitch'][...]
outer = self._f['geometry/lattices'][key]['outer'][0]
universe_ids = \
self._f['geometry/lattices'][key]['universes'][...]
# Create the Lattice
lattice = openmc.HexLattice(lattice_id=lattice_id)
lattice.num_rings(n_rings)
lattice.num_axial = n_axial
lattice.center = center
lattice.pitch = pitch
# If the Universe specified outer the Lattice is not void (-22)
if outer != -22:
lattice.outer = self.universes[outer]
# Build array of Universe pointers for the Lattice
universes = \
np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe)
for i in range(universe_ids.shape[0]):
for j in range(universe_ids.shape[1]):
for k in range(universe_ids.shape[2]):
if universe_ids[i,j,k] != -1:
universes[i,j,k] = \
self.get_universe_by_id(universe_ids[i,j,k])
# Add the Lattice to the global dictionary of all Lattices
self.lattices[index] = lattice
def _finalize_geometry(self):
# Initialize Geometry object
self.openmc_geometry = openmc.Geometry()
# Iterate over all Cells and add fill Materials, Universes and Lattices
for cell_key in self._cell_fills.keys():
# Determine fill type ('normal', 'universe', or 'lattice') and ID
fill_type = self._cell_fills[cell_key][0]
fill_id = self._cell_fills[cell_key][1]
# Retrieve the object corresponding to the fill type and ID
if fill_type == 'normal':
if fill_id > 0:
fill = self.get_material_by_id(fill_id)
else:
fill = 'void'
elif fill_type == 'universe':
fill = self.get_universe_by_id(fill_id)
else:
fill = self.get_lattice_by_id(fill_id)
# Set the fill for the Cell
self.cells[cell_key].fill = fill
# Set the root universe for the Geometry
root_universe = self.get_universe_by_id(0)
self.openmc_geometry.root_universe = root_universe
def make_opencg_geometry(self):
if self.opencg_geometry is None:
self.opencg_geometry = get_opencg_geometry(self.openmc_geometry)
def get_nuclide_by_zaid(self, zaid):
for index, nuclide in self.nuclides.items():
if nuclide._zaid == zaid:
return nuclide
return None
def get_material_by_id(self, material_id):
for index, material in self.materials.items():
if material._id == material_id:
return material
return None
def get_surface_by_id(self, surface_id):
for index, surface in self.surfaces.items():
if surface._id == surface_id:
return surface
return None
def get_cell_by_id(self, cell_id):
for index, cell in self.cells.items():
if cell._id == cell_id:
return cell
return None
def get_universe_by_id(self, universe_id):
for index, universe in self.universes.items():
if universe._id == universe_id:
return universe
return None
def get_lattice_by_id(self, lattice_id):
for index, lattice in self.lattices.items():
if lattice._id == lattice_id:
return lattice
return None

731
src/utils/openmc/surface.py Normal file
View file

@ -0,0 +1,731 @@
from xml.etree import ElementTree as ET
from openmc.checkvalue import *
from openmc.constants import BC_TYPES
# A static variable for auto-generated Surface IDs
AUTO_SURFACE_ID = 10000
def reset_auto_surface_id():
global AUTO_SURFACE_ID
AUTO_SURFACE_ID = 10000
class Surface(object):
def __init__(self, surface_id=None, boundary_type='transmission', name=''):
# Initialize class attributes
self.id = surface_id
self.name = name
self._type = ''
self.boundary_type = boundary_type
# A dictionary of the quadratic surface coefficients
# Key - coefficeint name
# Value - coefficient value
self._coeffs = {}
# An ordered list of the coefficient names to export to XML in the
# proper order
self._coeff_keys = []
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@property
def type(self):
return self._type
@property
def boundary_type(self):
return self._boundary_type
@property
def coeffs(self):
return self._coeffs
@id.setter
def id(self, surface_id):
if surface_id is None:
global AUTO_SURFACE_ID
self._id = AUTO_SURFACE_ID
AUTO_SURFACE_ID += 1
# Check that the ID is an integer and wasn't already used
elif not is_integer(surface_id):
msg = 'Unable to set a non-integer Surface ' \
'ID {0}'.format(surface_id)
raise ValueError(msg)
elif surface_id < 0:
msg = 'Unable to set Surface ID to {0} since it must be a ' \
'non-negative integer'.format(surface_id)
raise ValueError(msg)
else:
self._id = surface_id
@name.setter
def name(self, name):
if not is_string(name):
msg = 'Unable to set name for Surface ID={0} with a non-string ' \
'value {1}'.format(self._id, name)
raise ValueError(msg)
else:
self._name = name
@boundary_type.setter
def boundary_type(self, boundary_type):
if not is_string(boundary_type):
msg = 'Unable to set boundary type for Surface ID={0} with a ' \
'non-string value {1}'.format(self._id, boundary_type)
raise ValueError(msg)
elif not boundary_type in BC_TYPES.values():
msg = 'Unable to set boundary type for Surface ID={0} to ' \
'{1} which is not trasmission, vacuum or ' \
'reflective'.format(boundary_type)
raise ValueError(msg)
else:
self._boundary_type = boundary_type
def __repr__(self):
string = 'Surface\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type)
string += '{0: <16}{1}{2}\n'.format('\tBoundary', '=\t', self._boundary_type)
coeffs = '{0: <16}'.format('\tCoefficients') + '\n'
for coeff in self._coeffs:
coeffs += '{0: <16}{1}{2}\n'.format(coeff, '=\t', self._coeffs[coeff])
string += coeffs
return string
def create_xml_subelement(self):
element = ET.Element("surface")
element.set("id", str(self._id))
element.set("type", self._type)
element.set("boundary", self._boundary_type)
coeffs = ''
for coeff in self._coeff_keys:
coeffs += '{0} '.format(self._coeffs[coeff])
element.set("coeffs", coeffs.rstrip(' '))
return element
class Plane(Surface):
def __init__(self, surface_id=None, boundary_type='transmission',
A=None, B=None, C=None, D=None, name='',):
# Initialize Plane class attributes
super(Plane, self).__init__(surface_id, boundary_type, name=name)
self._type = 'plane'
self._coeff_keys = ['A', 'B', 'C', 'D']
if not A is None:
self.a = A
if not B is None:
self.b = B
if not C is None:
self.c = C
if not D is None:
self.d = D
@property
def a(self):
return self.coeffs['A']
@property
def b(self):
return self.coeffs['B']
@property
def c(self):
return self.coeffs['C']
@property
def d(self):
return self.coeffs['D']
@a.setter
def a(self, A):
if not is_integer(A) and not is_float(A):
msg = 'Unable to set A coefficient for Plane ID={0} to a ' \
'non-integer value {1}'.format(self._id, A)
raise ValueError(msg)
self._coeffs['A'] = A
@b.setter
def b(self, B):
if not is_integer(B) and not is_float(B):
msg = 'Unable to set B coefficient for Plane ID={0} to a ' \
'non-integer value {1}'.format(self._id, B)
raise ValueError(msg)
self._coeffs['B'] = B
@c.setter
def c(self, C):
if not is_integer(C) and not is_float(C):
msg = 'Unable to set C coefficient for Plane ID={0} to a ' \
'non-integer value {1}'.format(self._id, C)
raise ValueError(msg)
self._coeffs['C'] = C
@d.setter
def d(self, D):
if not is_integer(D) and not is_float(D):
msg = 'Unable to set D coefficient for Plane ID={0} to a ' \
'non-integer value {1}'.format(self._id, D)
raise ValueError(msg)
self._coeffs['D'] = D
class XPlane(Plane):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=None, name=''):
# Initialize XPlane class attributes
super(XPlane, self).__init__(surface_id, boundary_type, name=name)
self._type = 'x-plane'
self._coeff_keys = ['x0']
if not x0 is None:
self.x0 = x0
@property
def x0(self):
return self.coeff['x0']
@x0.setter
def x0(self, x0):
if not is_integer(x0) and not is_float(x0):
msg = 'Unable to set x0 coefficient for XPlane ID={0} to a ' \
'non-integer value {1}'.format(self._id, x0)
raise ValueError(msg)
self._coeffs['x0'] = x0
class YPlane(Plane):
def __init__(self, surface_id=None, boundary_type='transmission',
y0=None, name=''):
# Initialize YPlane class attributes
super(YPlane, self).__init__(surface_id, boundary_type, name=name)
self._type = 'y-plane'
self._coeff_keys = ['y0']
if not y0 is None:
self.y0 = y0
@property
def y0(self):
return self.coeffs['y0']
@y0.setter
def y0(self, y0):
if not is_integer(y0) and not is_float(y0):
msg = 'Unable to set y0 coefficient for XPlane ID={0} to a ' \
'non-integer value {1}'.format(self._id, y0)
raise ValueError(msg)
self._coeffs['y0'] = y0
class ZPlane(Plane):
def __init__(self, surface_id=None, boundary_type='transmission',
z0=None, name=''):
# Initialize ZPlane class attributes
super(ZPlane, self).__init__(surface_id, boundary_type, name=name)
self._type = 'z-plane'
self._coeff_keys = ['z0']
if not z0 is None:
self.z0 = z0
@property
def z0(self):
return self.coeffs['z0']
@z0.setter
def z0(self, z0):
if not is_integer(z0) and not is_float(z0):
msg = 'Unable to set z0 coefficient for ZPlane ID={0} to a ' \
'non-integer value {1}'.format(self._id, z0)
raise ValueError(msg)
self._coeffs['z0'] = z0
class Cylinder(Surface):
def __init__(self, surface_id=None, boundary_type='transmission',
R=None, name=''):
# Initialize Cylinder class attributes
super(Cylinder, self).__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['R']
if not R is None:
self.r = R
@property
def r(self):
return self.coeffs['R']
@r.setter
def r(self, R):
if not is_integer(R) and not is_float(R):
msg = 'Unable to set R coefficient for Cylinder ID={0} to a ' \
'non-integer value {1}'.format(self._id, R)
raise ValueError(msg)
self._coeffs['R'] = R
class XCylinder(Cylinder):
def __init__(self, surface_id=None, boundary_type='transmission',
y0=None, z0=None, R=None, name=''):
# Initialize XCylinder class attributes
super(XCylinder, self).__init__(surface_id, boundary_type, R, name=name)
self._type = 'x-cylinder'
self._coeff_keys = ['y0', 'z0', 'R']
if not y0 is None:
self.y0 = y0
if not z0 is None:
self.z0 = z0
@property
def y0(self):
return self.coeffs['y0']
@property
def z0(self):
return self.coeffs['z0']
@y0.setter
def y0(self, y0):
if not is_integer(y0) and not is_float(y0):
msg = 'Unable to set y0 coefficient for XCylinder ID={0} to a ' \
'non-integer value {1}'.format(self._id, y0)
raise ValueError(msg)
self._coeffs['y0'] = y0
@z0.setter
def z0(self, z0):
if not is_integer(z0) and not is_float(z0):
msg = 'Unable to set z0 coefficient for XCylinder ID={0} to a ' \
'non-integer value {1}'.format(self._id, z0)
raise ValueError(msg)
self._coeffs['z0'] = z0
class YCylinder(Cylinder):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=None, z0=None, R=None, name=''):
# Initialize YCylinder class attributes
super(YCylinder, self).__init__(surface_id, boundary_type, R, name=name)
self._type = 'y-cylinder'
self._coeff_keys = ['x0', 'z0', 'R']
if not x0 is None:
self.x0 = x0
if not z0 is None:
self.z0 = z0
@property
def x0(self):
return self.coeffs['x0']
@property
def z0(self):
return self.coeffs['z0']
@x0.setter
def x0(self, x0):
if not is_integer(x0) and not is_float(x0):
msg = 'Unable to set x0 coefficient for YCylinder ID={0} to a ' \
'non-integer value {1}'.format(self._id, x0)
raise ValueError(msg)
self._coeffs['x0'] = x0
@z0.setter
def z0(self, z0):
if not is_integer(z0) and not is_float(z0):
msg = 'Unable to set z0 coefficient for YCylinder ID={0} to a ' \
'non-integer value {1}'.format(self._id, z0)
raise ValueError(msg)
self._coeffs['z0'] = z0
class ZCylinder(Cylinder):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=None, y0=None, R=None, name=''):
# Initialize ZCylinder class attributes
# Initialize YPlane class attributes
super(ZCylinder, self).__init__(surface_id, boundary_type, R, name=name)
self._type = 'z-cylinder'
self._coeff_keys = ['x0', 'y0', 'R']
if not x0 is None:
self.x0 = x0
if not y0 is None:
self.y0 = y0
@property
def x0(self):
return self.coeffs['x0']
@property
def y0(self):
return self.coeffs['y0']
@x0.setter
def x0(self, x0):
if not is_integer(x0) and not is_float(x0):
msg = 'Unable to set x0 coefficient for ZCylinder ID={0} to a ' \
'non-integer value {1}'.format(self._id, x0)
raise ValueError(msg)
self._coeffs['x0'] = x0
@y0.setter
def y0(self, y0):
if not is_integer(y0) and not is_float(y0):
msg = 'Unable to set y0 coefficient for ZCylinder ID={0} to a ' \
'non-integer value {1}'.format(self._id, y0)
raise ValueError(msg)
self._coeffs['y0'] = y0
class Sphere(Surface):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=None, y0=None, z0=None, R=None, name=''):
# Initialize Sphere class attributes
super(Sphere, self).__init__(surface_id, boundary_type, name=name)
self._type = 'sphere'
self._coeff_keys = ['x0', 'y0', 'z0', 'R']
if not x0 is None:
self.x0 = x0
if not y0 is None:
self.y0 = y0
if not z0 is None:
self.z0 = z0
if not R is None:
self.r = R
@property
def x0(self):
return self.coeffs['x0']
@property
def y0(self):
return self.coeffs['y0']
@property
def z0(self):
return self.coeffs['z0']
@property
def r(self):
return self.coeffs['R']
@x0.setter
def x0(self, x0):
if not is_integer(x0) and not is_float(x0):
msg = 'Unable to set x0 coefficient for Sphere ID={0} to a ' \
'non-integer value {1}'.format(self._id, x0)
raise ValueError(msg)
self._coeffs['x0'] = x0
@y0.setter
def y0(self, y0):
if not is_integer(y0) and not is_float(y0):
msg = 'Unable to set y0 coefficient for Sphere ID={0} to a ' \
'non-integer value {1}'.format(self._id, y0)
raise ValueError(msg)
self._coeffs['y0'] = y0
@z0.setter
def z0(self, z0):
if not is_integer(z0) and not is_float(z0):
msg = 'Unable to set z0 coefficient for Sphere ID={0} to a ' \
'non-integer value {1}'.format(self._id, z0)
raise ValueError(msg)
self._coeffs['z0'] = z0
@r.setter
def r(self, R):
if not is_integer(R) and not is_float(R):
msg = 'Unable to set R coefficient for Sphere ID={0} to a ' \
'non-integer value {1}'.format(self._id, R)
raise ValueError(msg)
self._coeffs['R'] = R
class Cone(Surface):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=None, y0=None, z0=None, R2=None, name=''):
# Initialize Cone class attributes
super(Cone, self).__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['x0', 'y0', 'z0', 'R2']
if not x0 is None:
self.x0 = x0
if not y0 is None:
self.y0 = y0
if not z0 is None:
self.z0 = z0
if not R2 is None:
self.r2 = R2
@property
def x0(self):
return self.coeffs['x0']
@property
def y0(self):
return self.coeffs['y0']
@property
def z0(self):
return self.coeffs['z0']
@property
def r2(self):
return self.coeffs['r2']
@x0.setter
def x0(self, x0):
if not is_integer(x0) and not is_float(x0):
msg = 'Unable to set x0 coefficient for Cone ID={0} to a ' \
'non-integer value {1}'.format(self._id, x0)
raise ValueError(msg)
self._coeffs['x0'] = x0
@y0.setter
def y0(self, y0):
if not is_integer(y0) and not is_float(y0):
msg = 'Unable to set y0 coefficient for Cone ID={0} to a ' \
'non-integer value {1}'.format(self._id, y0)
raise ValueError(msg)
self._coeffs['y0'] = y0
@z0.setter
def z0(self, z0):
if not is_integer(z0) and not is_float(z0):
msg = 'Unable to set z0 coefficient for Cone ID={0} to a ' \
'non-integer value {1}'.format(self._id, z0)
raise ValueError(msg)
self._coeffs['z0'] = z0
@r2.setter
def r2(self, R2):
if not is_integer(R2) and not is_float(R2):
msg = 'Unable to set R^2 coefficient for Cone ID={0} to a ' \
'non-integer value {1}'.format(self._id, R2)
raise ValueError(msg)
self._coeffs['R2'] = R2
class XCone(Cone):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=None, y0=None, z0=None, R2=None, name=''):
# Initialize XCone class attributes
super(XCone, self).__init__(surface_id, boundary_type, x0, y0,
z0, R2, name=name)
self._type = 'x-cone'
class YCone(Cone):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=None, y0=None, z0=None, R2=None, name=''):
# Initialize YCone class attributes
super(YCone, self).__init__(surface_id, boundary_type, x0, y0, z0,
R2, name=name)
self._type = 'y-cone'
class ZCone(Cone):
def __init__(self, surface_id=None, boundary_type='transmission',
x0=None, y0=None, z0=None, R2=None, name=''):
# Initialize ZCone class attributes
super(ZCone, self).__init__(surface_id, boundary_type, x0, y0, z0,
R2, name=name)
self._type = 'z-cone'

1365
src/utils/openmc/tallies.py Normal file

File diff suppressed because it is too large Load diff

1499
src/utils/openmc/universe.py Normal file

File diff suppressed because it is too large Load diff

View file

@ -2,10 +2,10 @@
from distutils.core import setup
setup(name='statepoint',
version='0.6.1',
description='OpenMC StatePoint',
author='Paul Romano',
author_email='paul.k.romano@gmail.com',
setup(name='openmc',
version='0.6.2',
description='OpenMC Python API',
author='Will Boyd',
author_email='wbinventor@gmail.com',
url='https://github.com/mit-crpg/openmc',
py_modules=['statepoint'])
packages=['openmc'])

View file

@ -1,645 +0,0 @@
#!/usr/bin/env python2
import struct
from collections import OrderedDict
import numpy as np
import scipy.stats
REVISION_STATEPOINT = 13
filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn',
5: 'surface', 6: 'mesh', 7: 'energyin', 8: 'energyout'}
score_types = {-1: 'flux',
-2: 'total',
-3: 'scatter',
-4: 'nu-scatter',
-5: 'scatter-n',
-6: 'scatter-pn',
-7: 'nu-scatter-n',
-8: 'nu-scatter-pn',
-9: 'transport',
-10: 'n1n',
-11: 'absorption',
-12: 'fission',
-13: 'nu-fission',
-14: 'kappa-fission',
-15: 'current',
-16: 'flux-yn',
-17: 'total-yn',
-18: 'scatter-yn',
-19: 'nu-scatter-yn',
-20: 'events',
1: '(n,total)',
2: '(n,elastic)',
4: '(n,level)',
11: '(n,2nd)',
16: '(n,2n)',
17: '(n,3n)',
18: '(n,fission)',
19: '(n,f)',
20: '(n,nf)',
21: '(n,2nf)',
22: '(n,na)',
23: '(n,n3a)',
24: '(n,2na)',
25: '(n,3na)',
28: '(n,np)',
29: '(n,n2a)',
30: '(n,2n2a)',
32: '(n,nd)',
33: '(n,nt)',
34: '(n,nHe-3)',
35: '(n,nd2a)',
36: '(n,nt2a)',
37: '(n,4n)',
38: '(n,3nf)',
41: '(n,2np)',
42: '(n,3np)',
44: '(n,n2p)',
45: '(n,npa)',
91: '(n,nc)',
101: '(n,disappear)',
102: '(n,gamma)',
103: '(n,p)',
104: '(n,d)',
105: '(n,t)',
106: '(n,3He)',
107: '(n,a)',
108: '(n,2a)',
109: '(n,3a)',
111: '(n,2p)',
112: '(n,pa)',
113: '(n,t2a)',
114: '(n,d2a)',
115: '(n,pd)',
116: '(n,pt)',
117: '(n,da)',
201: '(n,Xn)',
202: '(n,Xgamma)',
203: '(n,Xp)',
204: '(n,Xd)',
205: '(n,Xt)',
206: '(n,X3He)',
207: '(n,Xa)',
444: '(damage)',
649: '(n,pc)',
699: '(n,dc)',
749: '(n,tc)',
799: '(n,3Hec)',
849: '(n,tc)'}
score_types.update({MT: '(n,n' + str(MT-50) + ')' for MT in range(51,91)})
score_types.update({MT: '(n,p' + str(MT-600) + ')' for MT in range(600,649)})
score_types.update({MT: '(n,d' + str(MT-650) + ')' for MT in range(650,699)})
score_types.update({MT: '(n,t' + str(MT-700) + ')' for MT in range(700,749)})
score_types.update({MT: '(n,3He' + str(MT-750) + ')' for MT in range(750,649)})
score_types.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)})
class Mesh(object):
def __init__(self):
pass
def __repr__(self):
if hasattr(self, "dimension"):
return "<Mesh: {0}>".format(tuple(self.dimension))
else:
return "<Mesh>"
class Filter(object):
def __init__(self):
self.type = 0
self.bins = []
def __repr__(self):
return "<Filter: {0}>".format(self.type)
class Tally(object):
def __init__(self):
self.filters = OrderedDict()
class SourceSite(object):
def __init__(self):
self.weight = None
self.xyz = None
self.uvw = None
self.E = None
def __repr__(self):
return "<SourceSite: xyz={0} at E={1}>".format(self.xyz, self.E)
class StatePoint(object):
def __init__(self, filename):
if filename.endswith('.h5'):
import h5py
self._f = h5py.File(filename, 'r')
self._hdf5 = True
else:
self._f = open(filename, 'rb')
self._hdf5 = False
# Set flags for what data was read
self._metadata = False
self._results = False
self._source = False
# Initialize arrays for meshes and tallies
self.meshes = []
self.tallies = []
self.source = []
# Read all metadata
self._read_metadata()
def _read_metadata(self):
# Read filetype
self.filetype = self._get_int(path='filetype')[0]
# Read statepoint revision
self.revision = self._get_int(path='revision')[0]
if self.revision != REVISION_STATEPOINT:
raise Exception('Statepoint Revision is not consistent.')
# Read OpenMC version
if self._hdf5:
self.version = [self._get_int(path='version_major')[0],
self._get_int(path='version_minor')[0],
self._get_int(path='version_release')[0]]
else:
self.version = self._get_int(3)
# Read date and time
self.date_and_time = self._get_string(19, path='date_and_time')
# Read path
self.path = self._get_string(255, path='path').strip()
# Read random number seed
self.seed = self._get_long(path='seed')[0]
# Read run information
self.run_mode = self._get_int(path='run_mode')[0]
self.n_particles = self._get_long(path='n_particles')[0]
# Read current batch
self.current_batch = self._get_int(path='current_batch')[0]
# Read criticality information
if self.run_mode == 2:
self.n_inactive = self._get_int(path='n_inactive')[0]
self.gen_per_batch = self._get_int(path='gen_per_batch')[0]
self.k_batch = self._get_double(
self.current_batch*self.gen_per_batch, path='k_generation')
self.entropy = self._get_double(
self.current_batch*self.gen_per_batch, path='entropy')
self.k_col_abs = self._get_double(path='k_col_abs')[0]
self.k_col_tra = self._get_double(path='k_col_tra')[0]
self.k_abs_tra = self._get_double(path='k_abs_tra')[0]
self.k_combined = self._get_double(2, path='k_combined')
# Read CMFD information
cmfd_present = self._get_int(path='cmfd_on')[0]
if cmfd_present == 1:
self.cmfd_indices = self._get_int(4, path='cmfd/indices')
self.k_cmfd = self._get_double(self.current_batch,
path='cmfd/k_cmfd')
self.cmfd_src = self._get_double_array(np.product(self.cmfd_indices),
path='cmfd/cmfd_src')
self.cmfd_src = np.reshape(self.cmfd_src,
tuple(self.cmfd_indices), order='F')
self.cmfd_entropy = self._get_double(self.current_batch,
path='cmfd/cmfd_entropy')
self.cmfd_balance = self._get_double(self.current_batch,
path='cmfd/cmfd_balance')
self.cmfd_dominance = self._get_double(self.current_batch,
path='cmfd/cmfd_dominance')
self.cmfd_srccmp = self._get_double(self.current_batch,
path='cmfd/cmfd_srccmp')
# Read number of meshes
n_meshes = self._get_int(path='tallies/n_meshes')[0]
# Read meshes
for i in range(n_meshes):
m = Mesh()
self.meshes.append(m)
base = 'tallies/mesh' + str(i+1) + '/'
# Read id, mesh type, and number of dimensions
m.id = self._get_int(path=base+'id')[0]
m.type = self._get_int(path=base+'type')[0]
n = self._get_int(path=base+'n_dimension')[0]
# Read mesh size, lower-left coordinates, upper-right coordinates,
# and width of each mesh cell
m.dimension = self._get_int(n, path=base+'dimension')
m.lower_left = self._get_double(n, path=base+'lower_left')
m.upper_right = self._get_double(n, path=base+'upper_right')
m.width = self._get_double(n, path=base+'width')
# Read number of tallies
n_tallies = self._get_int(path='tallies/n_tallies')[0]
for i in range(n_tallies):
# Create Tally object and add to list of tallies
t = Tally()
self.tallies.append(t)
base = 'tallies/tally' + str(i+1) + '/'
# Read id and number of realizations
t.id = self._get_int(path=base+'id')[0]
t.n_realizations = self._get_int(path=base+'n_realizations')[0]
# Read sizes of tallies
t.total_score_bins = self._get_int(path=base+'total_score_bins')[0]
t.total_filter_bins = self._get_int(path=base+'total_filter_bins')[0]
# Read number of filters
n_filters = self._get_int(path=base+'n_filters')[0]
for j in range(n_filters):
# Create Filter object
f = Filter()
base = 'tallies/tally{0}/filter{1}/'.format(i+1, j+1)
# Get type of filter
f.type = filter_types[self._get_int(path=base+'type')[0]]
# Add to filter dictionary
t.filters[f.type] = f
# Determine how many bins are in this filter
f.length = self._get_int(path=base+'n_bins')[0]
assert f.length > 0
if f.type == 'energyin' or f.type == 'energyout':
f.bins = self._get_double(f.length + 1, path=base+'bins')
elif f.type == 'mesh':
f.bins = self._get_int(path=base+'bins')
else:
f.bins = self._get_int(f.length, path=base+'bins')
base = 'tallies/tally' + str(i+1) + '/'
# Read nuclide bins
n_nuclides = self._get_int(path=base+'n_nuclide_bins')[0]
t.n_nuclides = n_nuclides
t.nuclides = self._get_int(n_nuclides, path=base+'nuclide_bins')
# Read score bins and user score bins
t.n_scores = self._get_int(path=base+'n_score_bins')[0]
t.scores = [score_types[j] for j in self._get_int(
t.n_scores, path=base+'score_bins')]
t.n_user_scores = self._get_int(path=base+'n_user_score_bins')[0]
# Read scattering moment order strings (e.g., P3, Y-1,2, etc.)
t.moments = list()
base += 'moments/'
# Extract the moment order string for each score
for k in range(t.n_scores):
moment = self._get_string(8, path=base+'order{0}'.format(k+1))
moment = moment.lstrip('[\'')
moment = moment.rstrip('\']')
# Remove extra whitespace
moment.replace(" ", "")
# Add the moment to the tally's list
t.moments.append(moment)
# Set up stride
stride = 1
for f in list(t.filters.values())[::-1]:
f.stride = stride
stride *= f.length
# Source bank present
source_present = self._get_int(path='source_present')[0]
if source_present == 1:
self.source_present = True
else:
self.source_present = False
# Set flag indicating metadata has already been read
self._metadata = True
def read_results(self):
# Check whether metadata has been read
if not self._metadata:
self._read_metadata()
# Number of realizations for global tallies
self.n_realizations = self._get_int(path='n_realizations')[0]
# Read global tallies
n_global_tallies = self._get_int(path='n_global_tallies')[0]
if self._hdf5:
data = self._f['global_tallies'].value
self.global_tallies = np.column_stack((data['sum'], data['sum_sq']))
else:
self.global_tallies = np.array(self._get_double(2*n_global_tallies))
self.global_tallies.shape = (n_global_tallies, 2)
# Flag indicating if tallies are present
tallies_present = self._get_int(path='tallies/tallies_present')[0]
# Read tally results
if tallies_present:
for i, t in enumerate(self.tallies):
n = t.total_score_bins * t.total_filter_bins
if self._hdf5:
path = 'tallies/tally{0}/results'.format(i+1)
data = self._f[path].value
t.results = np.column_stack((data['sum'], data['sum_sq']))
t.results.shape = (t.total_filter_bins, t.total_score_bins, 2)
else:
t.results = np.array(self._get_double(2*n))
t.results.shape = (t.total_filter_bins, t.total_score_bins, 2)
# Indicate that tally results have been read
self._results = True
def read_source(self):
# Check whether tally results have been read
if not self._results:
self.read_results()
# Check if source bank is in statepoint
if not self.source_present:
print('Source not in statepoint file.')
return
# For HDF5 state points, copy entire bank
if self._hdf5:
source_sites = self._f['source_bank'].value
for i in range(self.n_particles):
s = SourceSite()
self.source.append(s)
# Read position, angle, and energy
if self._hdf5:
s.weight, s.xyz, s.uvw, s.E = source_sites[i]
else:
s.weight = self._get_double()[0]
s.xyz = self._get_double(3)
s.uvw = self._get_double(3)
s.E = self._get_double()[0]
def generate_ci(self, confidence=0.95):
"""Calculates confidence intervals for each tally bin."""
# Determine number of realizations
n = self.n_realizations
# Determine significance level and percentile for two-sided CI
alpha = 1 - confidence
percentile = 1 - alpha/2
# Calculate t-value
t_value = scipy.stats.t.ppf(percentile, n - 1)
self.generate_stdev(t_value)
def generate_stdev(self, t_value=1.0):
"""
Calculates the sample mean and standard deviation of the mean for each
tally bin.
"""
# Determine number of realizations
n = self.n_realizations
# Global tallies
for i in range(len(self.global_tallies)):
# Get sum and sum of squares
s, s2 = self.global_tallies[i]
# Calculate sample mean and replace value
s /= n
self.global_tallies[i,0] = s
# Calculate standard deviation
if s != 0.0:
self.global_tallies[i,1] = t_value*np.sqrt((s2/n - s*s)/(n-1))
# Regular tallies
for t in self.tallies:
for i in range(t.results.shape[0]):
for j in range(t.results.shape[1]):
# Get sum and sum of squares
s, s2 = t.results[i,j]
# Calculate sample mean and replace value
s /= n
t.results[i,j,0] = s
# Calculate standard deviation
if s != 0.0:
t.results[i,j,1] = t_value*np.sqrt((s2/n - s*s)/(n-1))
def get_value(self, tally_index, spec_list, score_index):
"""Returns a tally score given a list of filters to satisfy.
Parameters
----------
tally_index : int
Index for tally in StatePoint.tallies list
spec_list : list
A list of tuples where the first value in each tuple is the filter
type, e.g. 'cell', and the second value is the desired index. If the
first value in the tuple is 'mesh', the second value should be a
tuple with three integers specifying the mesh indices.
Example: [('cell', 1), ('mesh', (14,17,20)), ('energyin', 2)]
score_index : int
Index corresponding to score for tally, i.e. the second index in
Tally.results[:,:,:].
"""
# Get Tally object given the index
t = self.tallies[tally_index]
# Initialize index for filter in Tally.results[:,:,:]
filter_index = 0
# Loop over specified filters in spec_list
for f_type, f_index in spec_list:
# Treat mesh filter separately
if f_type == 'mesh':
# Get index in StatePoint.meshes
mesh_index = t.filters['mesh'].bins[0] - 1
# Get dimensions of corresponding mesh
nx, ny, nz = self.meshes[mesh_index].dimension
# Convert (x,y,z) to a single bin -- this is similar to
# subroutine mesh_indices_to_bin in openmc/src/mesh.F90.
value = ((f_index[0] - 1)*ny*nz +
(f_index[1] - 1)*nz +
(f_index[2] - 1))
filter_index += value*t.filters[f_type].stride
else:
filter_index += f_index*t.filters[f_type].stride
# Return the desired result from Tally.results. This could be the sum and
# sum of squares, or it could be mean and stdev if self.generate_stdev()
# has been called already.
return t.results[filter_index, score_index]
def extract_results(self, tally_id, score_str):
"""Returns a tally results dictionary given a tally_id and score string.
Parameters
----------
tally_id : int
Index for the tally in StatePoint.tallies list
score_str : string
Corresponds to the string entered for a score in tallies.xml.
For a flux score extraction it would be 'score'
"""
# get tally
try:
tally = self.tallies[tally_id-1]
except:
print('Tally does not exist')
return
# get the score index if it is present
try:
idx = tally.scores.index(score_str)
except ValueError:
print('Score does not exist')
print(tally.scores)
return
# create numpy array for mean and 95% CI
n_bins = len(tally.results)
n_filters = len(tally.filters)
n_scores = len(tally.scores)
meanv = np.zeros(n_bins)
unctv = np.zeros(n_bins)
filters = np.zeros((n_bins,n_filters))
filtmax = np.zeros(n_filters+1)
meshmax = np.zeros(4)
filtmax[0] = 1
meshmax[0] = 1
# get number of realizations
n = tally.n_realizations
# get t-value
t_value = scipy.stats.t.ppf(0.975, n - 1)
# calculate mean
meanv = tally.results[:,idx,0]
meanv = meanv / n
# calculate 95% two-sided CI
unctv = tally.results[:,idx,1]
unctv = t_value*np.sqrt((unctv/n - meanv*meanv)/(n-1))/meanv
# create output dictionary
data = {'mean':meanv,'CI95':unctv}
# get bounds of filter bins
for akey in tally.filters.keys():
idx = list(tally.filters.keys()).index(akey)
filtmax[n_filters - idx] = tally.filters[akey].length
# compute bin info
for i in range(n_filters):
# compute indices for filter combination
filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) %
np.prod(filtmax[0:i+2]))/(np.prod(filtmax[0:i+1]))) + 1
# append in dictionary bin with filter
data.update({list(tally.filters.keys())[n_filters - i - 1]:
filters[:,n_filters - i - 1]})
# check for mesh
if list(tally.filters.keys())[n_filters - i - 1] == 'mesh':
dims = list(self.meshes[tally.filters['mesh'].bins[0] - 1].dimension)
dims.reverse()
dims = np.asarray(dims)
if score_str == 'current':
dims += 1
meshmax[1:4] = dims
mesh_bins = np.zeros((n_bins,3))
mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) %
np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1
mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) %
np.prod(meshmax[0:3]))/(np.prod(meshmax[0:2]))) + 1
mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) %
np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1
data.update({'mesh': list(zip(mesh_bins[:,0], mesh_bins[:,1],
mesh_bins[:,2]))})
i += 1
# add in maximum bin filters and order
b = list(tally.filters.keys())
b.reverse()
filtmax = list(filtmax[1:])
try:
idx = b.index('mesh')
filtmax[idx] = np.max(mesh_bins[:,2])
filtmax.insert(idx,np.max(mesh_bins[:,1]))
filtmax.insert(idx,np.max(mesh_bins[:,0]))
except ValueError:
pass
data.update({'bin_order':b,'bin_max':filtmax})
return data
def _get_data(self, n, typeCode, size):
return list(struct.unpack('={0}{1}'.format(n,typeCode),
self._f.read(n*size)))
def _get_int(self, n=1, path=None):
if self._hdf5:
return [int(v) for v in self._f[path].value]
else:
return [int(v) for v in self._get_data(n, 'i', 4)]
def _get_long(self, n=1, path=None):
if self._hdf5:
return [int(v) for v in self._f[path].value]
else:
return [int(v) for v in self._get_data(n, 'q', 8)]
def _get_float(self, n=1, path=None):
if self._hdf5:
return [float(v) for v in self._f[path].value]
else:
return [float(v) for v in self._get_data(n, 'f', 4)]
def _get_double(self, n=1, path=None):
if self._hdf5:
return [float(v) for v in self._f[path].value]
else:
return [float(v) for v in self._get_data(n, 'd', 8)]
def _get_double_array(self, n=1, path=None):
if self._hdf5:
return self._f[path].value
else:
return self._get_data(n, 'd', 8)
def _get_string(self, n=1, path=None):
if self._hdf5:
return str(self._f[path].value)
else:
return str(self._get_data(n, 's', 1)[0])

View file

@ -7,7 +7,7 @@ import itertools
import re
import warnings
from statepoint import StatePoint
from openmc.statepoint import StatePoint
alphanum = re.compile(r"[\W_]+")

View file

@ -4,7 +4,7 @@ import sys
from numpy.testing import assert_allclose, assert_equal
from statepoint import StatePoint
from openmc.statepoint import StatePoint
if len(sys.argv) > 2:
path1 = sys.argv[1]

View file

@ -8,7 +8,7 @@ import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
from statepoint import StatePoint
from openmc.statepoint import StatePoint
# Get filename
filename = argv[1]

View file

@ -9,7 +9,7 @@ import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
from statepoint import StatePoint
from openmc.statepoint import StatePoint
# Get filename
filename = argv[1]

View file

@ -23,7 +23,7 @@ import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
from statepoint import StatePoint
from openmc.statepoint import StatePoint
##################################### USER OPTIONS

View file

@ -2,15 +2,16 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string

View file

@ -3,41 +3,51 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.20.binary')
sp = StatePoint('statepoint.20.binary')
sp.read_results()
# extract tally results and convert to vector
results1 = sp.tallies[0].results
shape1 = results1.shape
size1 = (np.product(shape1))
results1 = np.reshape(results1, size1)
results2 = sp.tallies[1].results
shape2 = results2.shape
size2 = (np.product(shape2))
results2 = np.reshape(results2, size2)
results3 = sp.tallies[2].results
shape3 = results3.shape
size3 = (np.product(shape3))
results3 = np.reshape(results3, size3)
results4 = sp.tallies[3].results
shape4 = results4.shape
size4 = (np.product(shape4))
results4 = np.reshape(results4, size4)
tally1 = sp._tallies[1]
results1 = np.zeros((tally1._sum.size + tally1._sum.size, ))
results1[0::2] = tally1._sum.ravel()
results1[1::2] = tally1._sum_sq.ravel()
for tally_id in sp._tallies:
if 'CMFD flux, total, scatter-1' in sp._tallies[tally_id]._label:
tally2 = sp._tallies[tally_id]
elif 'CMFD neutron production' in sp._tallies[tally_id]._label:
tally3 = sp._tallies[tally_id]
elif 'CMFD surface currents' in sp._tallies[tally_id]._label:
tally4 = sp._tallies[tally_id]
results2 = np.zeros((tally2._sum.size + tally2._sum.size, ))
results2[0::2] = tally2._sum.ravel()
results2[1::2] = tally2._sum_sq.ravel()
results3 = np.zeros((tally3._sum.size + tally3._sum.size, ))
results3[0::2] = tally3._sum.ravel()
results3[1::2] = tally3._sum_sq.ravel()
results4 = np.zeros((tally4._sum.size + tally4._sum.size, ))
results4[0::2] = tally4._sum.ravel()
results4[1::2] = tally4._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tally 1:\n'
@ -55,25 +65,25 @@ for item in results4:
# write out cmfd answers
outstr += 'cmfd indices\n'
for item in sp.cmfd_indices:
for item in sp._cmfd_indices:
outstr += "{0:12.6E}\n".format(item)
outstr += 'k cmfd\n'
for item in sp.k_cmfd:
for item in sp._k_cmfd:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd entropy\n'
for item in sp.cmfd_entropy:
for item in sp._cmfd_entropy:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd balance\n'
for item in sp.cmfd_balance:
for item in sp._cmfd_balance:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd dominance ratio\n'
for item in sp.cmfd_dominance:
for item in sp._cmfd_dominance:
outstr += "{0:10.3E}\n".format(item)
outstr += 'cmfd openmc source comparison\n'
for item in sp.cmfd_srccmp:
for item in sp._cmfd_srccmp:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd source\n'
cmfdsrc = np.reshape(sp.cmfd_src, np.product(sp.cmfd_indices), order='F')
cmfdsrc = np.reshape(sp._cmfd_src, np.product(sp._cmfd_indices), order='F')
for item in cmfdsrc:
outstr += "{0:12.6E}\n".format(item)

View file

@ -3,41 +3,51 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.20.binary')
sp = StatePoint('statepoint.20.binary')
sp.read_results()
# extract tally results and convert to vector
results1 = sp.tallies[0].results
shape1 = results1.shape
size1 = (np.product(shape1))
results1 = np.reshape(results1, size1)
results2 = sp.tallies[1].results
shape2 = results2.shape
size2 = (np.product(shape2))
results2 = np.reshape(results2, size2)
results3 = sp.tallies[2].results
shape3 = results3.shape
size3 = (np.product(shape3))
results3 = np.reshape(results3, size3)
results4 = sp.tallies[3].results
shape4 = results4.shape
size4 = (np.product(shape4))
results4 = np.reshape(results4, size4)
tally1 = sp._tallies[1]
results1 = np.zeros((tally1._sum.size + tally1._sum.size, ))
results1[0::2] = tally1._sum.ravel()
results1[1::2] = tally1._sum_sq.ravel()
for tally_id in sp._tallies:
if 'CMFD flux, total, scatter-1' in sp._tallies[tally_id]._label:
tally2 = sp._tallies[tally_id]
elif 'CMFD neutron production' in sp._tallies[tally_id]._label:
tally3 = sp._tallies[tally_id]
elif 'CMFD surface currents' in sp._tallies[tally_id]._label:
tally4 = sp._tallies[tally_id]
results2 = np.zeros((tally2._sum.size + tally2._sum.size, ))
results2[0::2] = tally2._sum.ravel()
results2[1::2] = tally2._sum_sq.ravel()
results3 = np.zeros((tally3._sum.size + tally3._sum.size, ))
results3[0::2] = tally3._sum.ravel()
results3[1::2] = tally3._sum_sq.ravel()
results4 = np.zeros((tally4._sum.size + tally4._sum.size, ))
results4[0::2] = tally4._sum.ravel()
results4[1::2] = tally4._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tally 1:\n'
@ -55,25 +65,25 @@ for item in results4:
# write out cmfd answers
outstr += 'cmfd indices\n'
for item in sp.cmfd_indices:
for item in sp._cmfd_indices:
outstr += "{0:12.6E}\n".format(item)
outstr += 'k cmfd\n'
for item in sp.k_cmfd:
for item in sp._k_cmfd:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd entropy\n'
for item in sp.cmfd_entropy:
for item in sp._cmfd_entropy:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd balance\n'
for item in sp.cmfd_balance:
for item in sp._cmfd_balance:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd dominance ratio\n'
for item in sp.cmfd_dominance:
for item in sp._cmfd_dominance:
outstr += "{0:10.3E}\n".format(item)
outstr += 'cmfd openmc source comparison\n'
for item in sp.cmfd_srccmp:
for item in sp._cmfd_srccmp:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd source\n'
cmfdsrc = np.reshape(sp.cmfd_src, np.product(sp.cmfd_indices), order='F')
cmfdsrc = np.reshape(sp._cmfd_src, np.product(sp._cmfd_indices), order='F')
for item in cmfdsrc:
outstr += "{0:12.6E}\n".format(item)

View file

@ -3,41 +3,51 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.20.binary')
sp = StatePoint('statepoint.20.binary')
sp.read_results()
# extract tally results and convert to vector
results1 = sp.tallies[0].results
shape1 = results1.shape
size1 = (np.product(shape1))
results1 = np.reshape(results1, size1)
results2 = sp.tallies[1].results
shape2 = results2.shape
size2 = (np.product(shape2))
results2 = np.reshape(results2, size2)
results3 = sp.tallies[2].results
shape3 = results3.shape
size3 = (np.product(shape3))
results3 = np.reshape(results3, size3)
results4 = sp.tallies[3].results
shape4 = results4.shape
size4 = (np.product(shape4))
results4 = np.reshape(results4, size4)
tally1 = sp._tallies[1]
results1 = np.zeros((tally1._sum.size + tally1._sum.size, ))
results1[0::2] = tally1._sum.ravel()
results1[1::2] = tally1._sum_sq.ravel()
for tally_id in sp._tallies:
if 'CMFD flux, total, scatter-1' in sp._tallies[tally_id]._label:
tally2 = sp._tallies[tally_id]
elif 'CMFD neutron production' in sp._tallies[tally_id]._label:
tally3 = sp._tallies[tally_id]
elif 'CMFD surface currents' in sp._tallies[tally_id]._label:
tally4 = sp._tallies[tally_id]
results2 = np.zeros((tally2._sum.size + tally2._sum.size, ))
results2[0::2] = tally2._sum.ravel()
results2[1::2] = tally2._sum_sq.ravel()
results3 = np.zeros((tally3._sum.size + tally3._sum.size, ))
results3[0::2] = tally3._sum.ravel()
results3[1::2] = tally3._sum_sq.ravel()
results4 = np.zeros((tally4._sum.size + tally4._sum.size, ))
results4[0::2] = tally4._sum.ravel()
results4[1::2] = tally4._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tally 1:\n'
@ -55,25 +65,25 @@ for item in results4:
# write out cmfd answers
outstr += 'cmfd indices\n'
for item in sp.cmfd_indices:
for item in sp._cmfd_indices:
outstr += "{0:12.6E}\n".format(item)
outstr += 'k cmfd\n'
for item in sp.k_cmfd:
for item in sp._k_cmfd:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd entropy\n'
for item in sp.cmfd_entropy:
for item in sp._cmfd_entropy:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd balance\n'
for item in sp.cmfd_balance:
for item in sp._cmfd_balance:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd dominance ratio\n'
for item in sp.cmfd_dominance:
for item in sp._cmfd_dominance:
outstr += "{0:10.3E}\n".format(item)
outstr += 'cmfd openmc source comparison\n'
for item in sp.cmfd_srccmp:
for item in sp._cmfd_srccmp:
outstr += "{0:12.6E}\n".format(item)
outstr += 'cmfd source\n'
cmfdsrc = np.reshape(sp.cmfd_src, np.product(sp.cmfd_indices), order='F')
cmfdsrc = np.reshape(sp._cmfd_src, np.product(sp._cmfd_indices), order='F')
for item in cmfdsrc:
outstr += "{0:12.6E}\n".format(item)

View file

@ -3,29 +3,31 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.zeros((tally._sum.size + tally._sum.size, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.7.binary')
sp = StatePoint('statepoint.7.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -3,27 +3,29 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'entropy:\n'
for item in sp.entropy:
for item in sp._entropy:
outstr += "{0:12.6E}\n".format(item)
# write results to file

View file

@ -3,29 +3,31 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.empty((tally._sum.size + tally._sum.size, ), dtype=np.float64)
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -3,29 +3,31 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.empty((tally._sum.size + tally._sum.size, ), dtype=np.float64)
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -3,29 +3,31 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.zeros((tally._sum.size + tally._sum.size, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -3,29 +3,30 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.zeros((tally._sum.size + tally._sum.size, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -3,29 +3,30 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.zeros((tally._sum.size + tally._sum.size, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -3,29 +3,31 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.zeros((tally._sum.size + tally._sum.size, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -3,29 +3,31 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.zeros((tally._sum.size + tally._sum.size, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -3,29 +3,31 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.zeros((tally._sum.size + tally._sum.size, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -3,29 +3,31 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.zeros((tally._sum.size + tally._sum.size, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write out tally results
outstr += 'tallies:\n'

View file

@ -3,26 +3,28 @@
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
results = sp.tallies[0].results
shape = results.shape
size = (np.product(shape))
results = np.reshape(results, size)
tally = sp._tallies[1]
results = np.zeros((tally._sum.size + tally._sum.size, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
# set up output string
outstr = ''
# write out tally results
outstr += 'tallies:\n'
for item in results:

View file

@ -2,23 +2,24 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -1,24 +1,27 @@
#!/usr/bin/env python
import sys
import numpy as np
sys.path.insert(0, '../../src/utils')
# import statepoint
sys.path.append('../../src/utils')
import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -1,16 +1,19 @@
#!/usr/bin/env python
import sys
import numpy as np
sys.path.insert(0, '../../src/utils')
# import statepoint
sys.path.append('../../src/utils')
import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
@ -18,7 +21,7 @@ outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,8 +2,9 @@
import sys
# import particle restart
sys.path.insert(0, '../../src/utils')
# import particle restart
import particle_restart as pr
# read in particle restart file

View file

@ -2,8 +2,9 @@
import sys
# import particle restart
sys.path.insert(0, '../../src/utils')
# import particle restart
import particle_restart as pr
# read in particle restart file

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

View file

@ -2,23 +2,25 @@
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# import statepoint
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
sp = StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
outstr += "{0:12.6E} {1:12.6E}\n".format(sp._k_combined[0], sp._k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:

Some files were not shown because too many files have changed in this diff Show more