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

Merge branch 'develop' of github.com:openmc-dev/openmc into develop

Browse source
This commit is contained in:
rockfool 2020-04-06 18:34:07 -04:00
commit e5052d1d0c
170 changed files with 32475 additions and 2725 deletions
Download patch file Download diff file Expand all files Collapse all files

7
CMakeLists.txt
View file

@ -74,6 +74,13 @@ if(HDF5_IS_PARALLEL)
message(STATUS "Using parallel HDF5")
endif()
# Version 1.12 of HDF5 deprecates the H5Oget_info_by_idx() interface.
# Thus, we give these flags to allow usage of the old interface in newer
# versions of HDF5.
if(NOT (${HDF5_VERSION} VERSION_LESS 1.12.0))
list(APPEND cxxflags -DH5Oget_info_by_idx_vers=1 -DH5O_info_t_vers=1)
endif()
#===============================================================================
# Set compile/link flags based on which compiler is being used
#===============================================================================

12
docs/source/examples/index.rst
View file

@ -23,7 +23,6 @@ General Usage
search
nuclear-data
nuclear-data-resonance-covariance
cad-geom
pincell-depletion
--------
@ -36,6 +35,7 @@ Geometry
hexagonal
triso
candu
cad-geom
------------------------------------
Multi-Group Cross Section Generation
@ -60,3 +60,13 @@ Multi-Group Mode
mg-mode-part-i
mg-mode-part-ii
mg-mode-part-iii
-----------------
Unstructured Mesh
-----------------
.. toctree::
:maxdepth: 1
unstructured-mesh-part-i
unstructured-mesh-part-ii

13
docs/source/examples/unstructured-mesh-part-i.rst Normal file
View file

@ -0,0 +1,13 @@
.. _notebook_unstructured_mesh_part_i:
===============================
Unstructured Mesh: Introduction
===============================
.. only:: html
.. notebook:: ../../../examples/jupyter/unstructured-mesh-part-i.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

13
docs/source/examples/unstructured-mesh-part-ii.rst Normal file
View file

@ -0,0 +1,13 @@
.. _notebook_unstructured_mesh_part_ii:
===================================================================================
Unstructured Mesh: Unstructured Mesh Tallies with CAD and Point Cloud Visualization
===================================================================================
.. only:: html
.. notebook:: ../../../examples/jupyter/unstructured-mesh-part-ii.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

4
docs/source/io_formats/statepoint.rst
View file

@ -77,6 +77,10 @@ The current version of the statepoint file format is 17.0.
of mesh.
- **width** (*double[]*) -- Width of each mesh cell in each
dimension.
- **Unstructured Mesh Only:**
- **volumes** (*double[]*) -- Volume of each mesh cell.
- **centroids** (*double[]*) -- Location of the mesh cell
centroids.
**/tallies/filters/**

33
docs/source/io_formats/tallies.rst
View file

@ -123,8 +123,8 @@ to the scored values. The ``filter`` element has the following
attributes/sub-elements:
:type:
The type of the filter. Accepted options are "cell", "cellfrom",
"cellborn", "surface", "material", "universe", "energy", "energyout", "mu",
The type of the filter. Accepted options are "cell", "cellfrom",
"cellborn", "surface", "material", "universe", "energy", "energyout", "mu",
"polar", "azimuthal", "mesh", "distribcell", "delayedgroup",
"energyfunction", and "particle".
@ -154,27 +154,27 @@ For each filter type, the following table describes what the ``bins`` attribute
should be set to:
:cell:
A list of unique IDs for cells in which the tally should be
A list of unique IDs for cells in which the tally should be
accumulated.
:surface:
This filter allows the tally to be scored when crossing a surface. A list of
surface IDs should be given. By default, net currents are tallied, and to
tally a partial current from one cell to another, this should be used in
This filter allows the tally to be scored when crossing a surface. A list of
surface IDs should be given. By default, net currents are tallied, and to
tally a partial current from one cell to another, this should be used in
combination with a cell or cell_from filter that defines the other cell.
This filter should not be used in combination with a meshfilter.
:cellfrom:
This filter allows the tally to be scored when crossing a surface and the
particle came from a specified cell. A list of cell IDs should be
This filter allows the tally to be scored when crossing a surface and the
particle came from a specified cell. A list of cell IDs should be
given.
To tally a partial current from a cell to another, this filter should be
To tally a partial current from a cell to another, this filter should be
used in combination with a cell filter, to define the other cell.
This filter should not be used in combination with a meshfilter.
:cellborn:
This filter allows the tally to be scored to only when particles were
originally born in a specified cell. A list of cell IDs should be
originally born in a specified cell. A list of cell IDs should be
given.
:material:
@ -276,7 +276,7 @@ should be set to:
<filter type="azimuthal" bins="2" />
:mesh:
The unique ID of a structured mesh to be tallied over.
The unique ID of a mesh to be tallied over.
:distribcell:
The single cell which should be tallied uniquely for all instances.
@ -307,12 +307,13 @@ should be set to:
``<mesh>`` Element
------------------
If a structured mesh is desired as a filter for a tally, it must be specified in
a separate element with the tag name ``<mesh>``. This element has the following
If a mesh is desired as a filter for a tally, it must be specified in a separate
element with the tag name ``<mesh>``. This element has the following
attributes/sub-elements:
:type:
The type of structured mesh. This can be either "regular" or "rectilinear".
The type of mesh. This can be either "regular", "rectilinear", or
"unstructured".
:dimension:
The number of mesh cells in each direction. (For regular mesh only.)
@ -337,6 +338,10 @@ attributes/sub-elements:
:z_grid:
The mesh divisions along the z-axis. (For rectilinear mesh only.)
:mesh_file:
The name of the mesh file to be loaded at runtime. (For unstructured mesh
only.)
.. note::
One of ``<upper_right>`` or ``<width>`` must be specified, but not both
(even if they are consistent with one another).

1
docs/source/pythonapi/base.rst
View file

@ -127,6 +127,7 @@ Constructing Tallies
openmc.ParticleFilter
openmc.RegularMesh
openmc.RectilinearMesh
openmc.UnstructuredMesh
openmc.Trigger
openmc.TallyDerivative
openmc.Tally

2
docs/source/usersguide/beginners.rst
View file

@ -59,7 +59,7 @@ deterministic methods:
Now let's look at the pros and cons of Monte Carlo methods:
- **Pro**: No mesh generation is required to build geometry. By using
`constructive solid geometry`_, it's possible to build arbitrarily complex
`constructive solid geometry`_, it's possible to build complex
models with curved surfaces.
- **Pro**: Monte Carlo methods can be used with either continuous-energy or

35
docs/source/usersguide/geometry.rst
View file

@ -185,7 +185,13 @@ the :class:`openmc.Cell` class::
In this example, an instance of :class:`openmc.Material` is assigned to the
:attr:`Cell.fill` attribute. One can also fill a cell with a :ref:`universe
<usersguide_universes>` or :ref:`lattice <usersguide_lattices>`.
<usersguide_universes>` or :ref:`lattice <usersguide_lattices>`. If you provide
no fill to a cell or assign a value of `None`, it will be treated as a "void"
cell with no material within. Particles are allowed to stream through the cell but
will undergo no collisions::
# This cell will be filled with void on export to XML
gap = openmc.Cell(region=pellet_gap)
The classes :class:`Halfspace`, :class:`Intersection`, :class:`Union`, and
:class:`Complement` and all instances of :class:`openmc.Region` and can be
@ -434,3 +440,30 @@ named ``dagmc.h5m``) when initializing a simulation. If a `geometry.xml
<https://svalinn.github.io/DAGMC/usersguide/tools.html#make-watertight>`_. Future
implementations of DAGMC geometry will support small volume overlaps and
un-merged surfaces.
-------------------------
Calculating Atoms Content
-------------------------
If the total volume occupied by all instances of a cell in the geometry is known
by the user, it is possible to assign this volume to a cell without performing a
:ref:`stochastic volume <usersguide_volume>` calculation::
from uncertainties import ufloat
# Set known total volume in [cc]
cell = openmc.Cell()
cell.volume = 17.0
# Set volume if it is known with some uncertainty
cell.volume = ufloat(17.0, 0.1)
Once a volume is set, and a cell is filled with a material or distributed
materials, it is possible to use the :func:`~openmc.Cell.atoms` method to obtain
a dictionary of nuclides and their total number of atoms in all instances
of a cell (e.g. ``{'H1': 1.0e22, 'O16': 0.5e22, ...}``)::
cell = openmc.Cell(fill = u02)
cell.volume = 17.0
O16_atoms = cell.atoms['O16']

6
docs/source/usersguide/processing.rst
View file

@ -41,9 +41,9 @@ Plotting in 2D
--------------
The :ref:`IPython notebook example <notebook_post_processing>` also demonstrates
how to plot a mesh tally in two dimensions using the Python API. One can also
use the :ref:`scripts_plot` script which provides an interactive GUI to explore
and plot mesh tallies for any scores and filter bins.
how to plot a structured mesh tally in two dimensions using the Python API. One
can also use the :ref:`scripts_plot` script which provides an interactive GUI to
explore and plot structured mesh tallies for any scores and filter bins.
.. image:: ../_images/plotmeshtally.png
:width: 400px

17
docs/source/usersguide/settings.rst
View file

@ -171,16 +171,17 @@ Shannon Entropy
To assess convergence of the source distribution, the scalar Shannon entropy
metric is often used in Monte Carlo codes. OpenMC also allows you to calculate
Shannon entropy at each generation over a specified mesh, created using the
:class:`openmc.Mesh` class. After instantiating a :class:`Mesh`, you need to
specify the lower-left coordinates of the mesh (:attr:`Mesh.lower_left`), the
number of mesh cells in each direction (:attr:`Mesh.dimension`) and either the
upper-right coordinates of the mesh (:attr:`Mesh.upper_right`) or the width of
each mesh cell (:attr:`Mesh.width`). Once you have a mesh, simply assign it to
the :attr:`Settings.entropy_mesh` attribute.
:class:`openmc.RegularMesh` class. After instantiating a :class:`RegularMesh`,
you need to specify the lower-left coordinates of the mesh
(:attr:`RegularMesh.lower_left`), the number of mesh cells in each direction
(:attr:`RegularMesh.dimension`) and either the upper-right coordinates of the
mesh (:attr:`RegularMesh.upper_right`) or the width of each mesh cell
(:attr:`RegularMesh.width`). Once you have a mesh, simply assign it to the
:attr:`Settings.entropy_mesh` attribute.
::
entropy_mesh = openmc.Mesh()
entropy_mesh = openmc.RegularMesh()
entropy_mesh.lower_left = (-50, -50, -25)
entropy_mesh.upper_right = (50, 50, 25)
entropy_mesh.dimension = (8, 8, 8)
@ -193,7 +194,7 @@ property::
geom = openmc.Geometry()
...
m = openmc.Mesh()
m = openmc.RegularMesh()
m.lower_left, m.upper_right = geom.bounding_box
m.dimension = (8, 8, 8)

140
examples/assembly/assembly.py Normal file
View file

@ -0,0 +1,140 @@
"""
This script builds a single PWR assembly and is a slightly more advanced
demonstration of model building using Python. The creation of two universes for
fuel pins and guide tube pins has been separated into functions, and then the
overall model is built by an `assembly` function. This script also demonstrates
the use of the `Model` class, which provides some extra convenience over using
`Geometry`, `Materials`, and `Settings` classes directly. Finally, the script
takes two command-line flags that indicate whether to build and/or run the
model.
"""
import argparse
from math import log10
import numpy as np
import openmc
# Define surfaces
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
# Define materials
fuel = openmc.Material(name='Fuel')
fuel.set_density('g/cm3', 10.29769)
fuel.add_nuclide('U234', 4.4843e-6)
fuel.add_nuclide('U235', 5.5815e-4)
fuel.add_nuclide('U238', 2.2408e-2)
fuel.add_nuclide('O16', 4.5829e-2)
clad = openmc.Material(name='Cladding')
clad.set_density('g/cm3', 6.55)
clad.add_nuclide('Zr90', 2.1827e-2)
clad.add_nuclide('Zr91', 4.7600e-3)
clad.add_nuclide('Zr92', 7.2758e-3)
clad.add_nuclide('Zr94', 7.3734e-3)
clad.add_nuclide('Zr96', 1.1879e-3)
hot_water = openmc.Material(name='Hot borated water')
hot_water.set_density('g/cm3', 0.740582)
hot_water.add_nuclide('H1', 4.9457e-2)
hot_water.add_nuclide('O16', 2.4672e-2)
hot_water.add_nuclide('B10', 8.0042e-6)
hot_water.add_nuclide('B11', 3.2218e-5)
hot_water.add_s_alpha_beta('c_H_in_H2O')
def fuel_pin():
"""Returns a fuel pin universe."""
fuel_cell = openmc.Cell(fill=fuel, region=-fuel_or)
clad_cell = openmc.Cell(fill=clad, region=+fuel_or & -clad_or)
hot_water_cell = openmc.Cell(fill=hot_water, region=+clad_or)
univ = openmc.Universe(name='Fuel Pin')
univ.add_cells([fuel_cell, clad_cell, hot_water_cell])
return univ
def guide_tube_pin():
"""Returns a control rod guide tube universe."""
gt_inner_cell = openmc.Cell(fill=hot_water, region=-fuel_or)
gt_clad_cell = openmc.Cell(fill=clad, region=+fuel_or & -clad_or)
gt_outer_cell = openmc.Cell(fill=hot_water, region=+clad_or)
univ = openmc.Universe(name='Guide Tube')
univ.add_cells([gt_inner_cell, gt_clad_cell, gt_outer_cell])
return univ
def assembly_model():
"""Returns a single PWR fuel assembly."""
model = openmc.model.Model()
# Create fuel assembly Lattice
pitch = 21.42
assembly = openmc.RectLattice(name='Fuel Assembly')
assembly.pitch = (pitch/17, pitch/17)
assembly.lower_left = (-pitch/2, -pitch/2)
# Create array indices for guide tube locations in lattice
gt_pos = np.array([
[2, 5], [2, 8], [2, 11],
[3, 3], [3, 13],
[5, 2], [5, 5], [5, 8], [5, 11], [5, 14],
[8, 2], [8, 5], [8, 8], [8, 11], [8, 14],
[11, 2], [11, 5], [11, 8], [11, 11], [11, 14],
[13, 3], [13, 13],
[14, 5], [14, 8], [14, 11]
])
# Create 17x17 array of universes. First we create a 17x17 array all filled
# with the fuel pin universe. Then, we replace the guide tube positions with
# the guide tube pin universe (note the use of numpy fancy indexing to
# achieve this).
assembly.universes = np.full((17, 17), fuel_pin())
assembly.universes[gt_pos[:, 0], gt_pos[:, 1]] = guide_tube_pin()
# Create outer boundary of the geometry to surround the lattice
outer_boundary = openmc.model.rectangular_prism(
pitch, pitch, boundary_type='reflective')
# Create a cell filled with the lattice
main_cell = openmc.Cell(fill=assembly, region=outer_boundary)
# Finally, create geometry by providing a list of cells that fill the root
# universe
model.geometry = openmc.Geometry([main_cell])
model.settings.batches = 150
model.settings.inactive = 50
model.settings.particles = 1000
model.settings.source = openmc.Source(space=openmc.stats.Box(
(-pitch/2, -pitch/2, -1),
(pitch/2, pitch/2, 1),
only_fissionable=True
))
# NOTE: We never actually created a Materials object. When you export/run
# using the Model object, if no materials were assigned it will look through
# the Geometry object and automatically export any materials that are
# necessary to build the model.
return model
if __name__ == '__main__':
# Set up command-line arguments for generating/running the model
parser = argparse.ArgumentParser()
parser.add_argument('--generate', action='store_true')
parser.add_argument('--run', action='store_true')
args = parser.parse_args()
if args.generate or args.run:
model = assembly_model()
if args.generate:
model.export_to_xml()
if args.run:
model.run()

0
examples/xml/custom_source/CMakeLists.txt → examples/custom_source/CMakeLists.txt
View file

19
examples/custom_source/README.md Normal file
View file

@ -0,0 +1,19 @@
# Building a Custom Source
To run this example, you first need to compile the custom source library, which
requires headers from OpenMC. A CMakeLists.txt file has been set up for you that
will search for OpenMC and build the custom library. To build the source
library, you can run:
mkdir build && cd build
OPENMC_ROOT=<path_to_openmc_install> cmake ..
make
After this, you can build the model by running `python build_xml.py`. In the XML
files that are created, you should see a reference to build/libsource.so, the
custom source library that was built by CMake. The model is also set up with a
mesh tally of the flux, so once you run `openmc`, you will get a statepoint file
with the tally results in it. Running `python show_flux.py` will pull in the
results from the statepoint file and display them. If all worked well, you
should see a ring "imprint" as well as a higher flux to the right side (since
the custom source has all particles moving in the positive x direction).

36
examples/custom_source/build_xml.py Normal file
View file

@ -0,0 +1,36 @@
import openmc
# Create a single material
iron = openmc.Material()
iron.set_density('g/cm3', 5.0)
iron.add_element('Fe', 1.0)
mats = openmc.Materials([iron])
mats.export_to_xml()
# Create a 5 cm x 5 cm box filled with iron
box = openmc.model.rectangular_prism(10.0, 10.0, boundary_type='vacuum')
cell = openmc.Cell(fill=iron, region=box)
geometry = openmc.Geometry([cell])
geometry.export_to_xml()
# Tell OpenMC we're going to use our custom source
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.batches = 10
settings.particles = 1000
source = openmc.Source()
source.library = 'build/libsource.so'
settings.source = source
settings.export_to_xml()
# Finally, define a mesh tally so that we can see the resulting flux
mesh = openmc.RegularMesh()
mesh.lower_left = (-5.0, -5.0)
mesh.upper_right = (5.0, 5.0)
mesh.dimension = (50, 50)
tally = openmc.Tally()
tally.filters = [openmc.MeshFilter(mesh)]
tally.scores = ['flux']
tallies = openmc.Tallies([tally])
tallies.export_to_xml()

14
examples/custom_source/show_flux.py Normal file
View file

@ -0,0 +1,14 @@
import matplotlib.pyplot as plt
import openmc
# Get the flux from the statepoint
with openmc.StatePoint('statepoint.10.h5') as sp:
flux = sp.tallies[1].mean
flux.shape = (50, 50)
# Plot the flux
fig, ax = plt.subplots()
ax.imshow(flux, origin='lower', extent=(-5.0, 5.0, -5.0, 5.0))
ax.set_xlabel('x [cm]')
ax.set_ylabel('y [cm]')
plt.show()

0
examples/xml/custom_source/source_ring.cpp → examples/custom_source/source_ring.cpp
View file

33
examples/jezebel/jezebel.py Normal file
View file

@ -0,0 +1,33 @@
import openmc
# Create plutonium metal material
pu = openmc.Material()
pu.set_density('sum')
pu.add_nuclide('Pu239', 3.7047e-02)
pu.add_nuclide('Pu240', 1.7512e-03)
pu.add_nuclide('Pu241', 1.1674e-04)
pu.add_element('Ga', 1.3752e-03)
mats = openmc.Materials([pu])
mats.export_to_xml()
# Create a single cell filled with the Pu metal
sphere = openmc.Sphere(r=6.3849, boundary_type='vacuum')
cell = openmc.Cell(fill=pu, region=-sphere)
geom = openmc.Geometry([cell])
geom.export_to_xml()
# Finally, define some run settings
settings = openmc.Settings()
settings.batches = 200
settings.inactive = 10
settings.particles = 10000
settings.export_to_xml()
# Run the simulation
openmc.run()
# Get the resulting k-effective value
n = settings.batches
with openmc.StatePoint(f'statepoint.{n}.h5') as sp:
keff = sp.k_combined
print(f'Final k-effective = {keff}')

BIN
examples/jupyter/images/manifold-cad.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 97 KiB

BIN
examples/jupyter/images/manifold_flux.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 60 KiB

BIN
examples/jupyter/images/manifold_pnt_cld.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 118 KiB

BIN
examples/jupyter/images/pin_mesh.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 129 KiB

BIN
examples/jupyter/images/umesh_flux.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 194 KiB

BIN
examples/jupyter/images/umesh_heating.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 129 KiB

BIN
examples/jupyter/images/umesh_w_assembly.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 110 KiB

704
examples/jupyter/unstructured-mesh-part-i.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

682
examples/jupyter/unstructured-mesh-part-ii.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

0
examples/python/lattice/hexagonal/build-xml.py → examples/lattice/hexagonal/build_xml.py
View file

0
examples/python/lattice/nested/build-xml.py → examples/lattice/nested/build_xml.py
View file

0
examples/python/lattice/simple/build-xml.py → examples/lattice/simple/build_xml.py
View file

112
examples/pincell/build_xml.py Normal file
View file

@ -0,0 +1,112 @@
from math import log10
import numpy as np
import openmc
###############################################################################
# Create materials for the problem
uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment')
uo2.set_density('g/cm3', 10.29769)
uo2.add_element('U', 1., enrichment=2.4)
uo2.add_element('O', 2.)
helium = openmc.Material(name='Helium for gap')
helium.set_density('g/cm3', 0.001598)
helium.add_element('He', 2.4044e-4)
zircaloy = openmc.Material(name='Zircaloy 4')
zircaloy.set_density('g/cm3', 6.55)
zircaloy.add_element('Sn', 0.014 , 'wo')
zircaloy.add_element('Fe', 0.00165, 'wo')
zircaloy.add_element('Cr', 0.001 , 'wo')
zircaloy.add_element('Zr', 0.98335, 'wo')
borated_water = openmc.Material(name='Borated water')
borated_water.set_density('g/cm3', 0.740582)
borated_water.add_element('B', 4.0e-5)
borated_water.add_element('H', 5.0e-2)
borated_water.add_element('O', 2.4e-2)
borated_water.add_s_alpha_beta('c_H_in_H2O')
# Collect the materials together and export to XML
materials = openmc.Materials([uo2, helium, zircaloy, borated_water])
materials.export_to_xml()
###############################################################################
# Define problem geometry
# Create cylindrical surfaces
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
# Create a region represented as the inside of a rectangular prism
pitch = 1.25984
box = openmc.rectangular_prism(pitch, pitch, boundary_type='reflective')
# Create cells, mapping materials to regions
fuel = openmc.Cell(fill=uo2, region=-fuel_or)
gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
water = openmc.Cell(fill=borated_water, region=+clad_or & box)
# Create a geometry and export to XML
geometry = openmc.Geometry([fuel, gap, clad, water])
geometry.export_to_xml()
###############################################################################
# Define problem settings
# Indicate how many particles to run
settings = openmc.Settings()
settings.batches = 100
settings.inactive = 10
settings.particles = 1000
# Create an initial uniform spatial source distribution over fissionable zones
lower_left = (-pitch/2, -pitch/2, -1)
upper_right = (pitch/2, pitch/2, 1)
uniform_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True)
settings.source = openmc.source.Source(space=uniform_dist)
# For source convergence checks, add a mesh that can be used to calculate the
# Shannon entropy
entropy_mesh = openmc.RegularMesh()
entropy_mesh.lower_left = (-fuel_or.r, -fuel_or.r)
entropy_mesh.upper_right = (fuel_or.r, fuel_or.r)
entropy_mesh.dimension = (10, 10)
settings.entropy_mesh = entropy_mesh
settings.export_to_xml()
###############################################################################
# Define tallies
# Create a mesh that will be used for tallying
mesh = openmc.RegularMesh()
mesh.dimension = (100, 100)
mesh.lower_left = (-pitch/2, -pitch/2)
mesh.upper_right = (pitch/2, pitch/2)
# Create a mesh filter that can be used in a tally
mesh_filter = openmc.MeshFilter(mesh)
# Now use the mesh filter in a tally and indicate what scores are desired
mesh_tally = openmc.Tally(name="Mesh tally")
mesh_tally.filters = [mesh_filter]
mesh_tally.scores = ['flux', 'fission', 'nu-fission']
# Let's also create a tally to get the flux energy spectrum. We start by
# creating an energy filter
e_min, e_max = 1e-5, 20.0e6
groups = 500
energies = np.logspace(log10(e_min), log10(e_max), groups + 1)
energy_filter = openmc.EnergyFilter(energies)
spectrum_tally = openmc.Tally(name="Flux spectrum")
spectrum_tally.filters = [energy_filter]
spectrum_tally.scores = ['flux']
# Instantiate a Tallies collection and export to XML
tallies = openmc.Tallies([mesh_tally, spectrum_tally])
tallies.export_to_xml()

23
examples/pincell/plot_spectrum.py Normal file
View file

@ -0,0 +1,23 @@
import matplotlib.pyplot as plt
import openmc
# Get results from statepoint
with openmc.StatePoint('statepoint.100.h5') as sp:
t = sp.get_tally(name="Flux spectrum")
# Get the energies from the energy filter
energy_filter = t.filters[0]
energies = energy_filter.bins[:, 0]
# Get the flux values
mean = t.get_values(value='mean').ravel()
uncertainty = t.get_values(value='std_dev').ravel()
# Plot flux spectrum
fix, ax = plt.subplots()
ax.loglog(energies, mean, drawstyle='steps-post')
ax.set_xlabel('Energy [eV]')
ax.set_ylabel('Flux')
ax.grid(True, which='both')
plt.show()

0
examples/python/pincell_depletion/chain_simple.xml → examples/pincell_depletion/chain_simple.xml
View file

0
examples/python/pincell_depletion/restart_depletion.py → examples/pincell_depletion/restart_depletion.py
View file

0
examples/python/pincell_depletion/run_depletion.py → examples/pincell_depletion/run_depletion.py
View file

131
examples/python/pincell_multigroup/build-xml.py → examples/pincell_multigroup/build_xml.py
View file

@ -1,20 +1,12 @@
from math import log10
import numpy as np
import openmc
import openmc.mgxs
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 100
inactive = 10
particles = 1000
###############################################################################
# Exporting to OpenMC mgxs.h5 file
###############################################################################
# Create multigroup data
# Instantiate the energy group data
groups = openmc.mgxs.EnergyGroups(group_edges=[
@ -69,107 +61,90 @@ mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata])
mg_cross_sections_file.export_to_hdf5()
###############################################################################
# Exporting to OpenMC materials.xml file
###############################################################################
# Create materials for the problem
# Instantiate some Macroscopic Data
uo2_data = openmc.Macroscopic('UO2')
h2o_data = openmc.Macroscopic('LWTR')
# Instantiate some Materials and register the appropriate Macroscopic objects
uo2 = openmc.Material(material_id=1, name='UO2 fuel')
uo2 = openmc.Material(name='UO2 fuel')
uo2.set_density('macro', 1.0)
uo2.add_macroscopic(uo2_data)
water = openmc.Material(material_id=2, name='Water')
water = openmc.Material(name='Water')
water.set_density('macro', 1.0)
water.add_macroscopic(h2o_data)
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([uo2, water])
materials_file.cross_sections = "./mgxs.h5"
materials_file.cross_sections = "mgxs.h5"
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml file
###############################################################################
# Define problem geometry
# Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.54, name='Fuel OR')
left = openmc.XPlane(surface_id=4, x0=-0.63, name='left')
right = openmc.XPlane(surface_id=5, x0=0.63, name='right')
bottom = openmc.YPlane(surface_id=6, y0=-0.63, name='bottom')
top = openmc.YPlane(surface_id=7, y0=0.63, name='top')
# Create a surface for the fuel outer radius
fuel_or = openmc.ZCylinder(r=0.54, name='Fuel OR')
left.boundary_type = 'reflective'
right.boundary_type = 'reflective'
top.boundary_type = 'reflective'
bottom.boundary_type = 'reflective'
# Create a region represented as the inside of a rectangular prism
pitch = 1.26
box = openmc.rectangular_prism(pitch, pitch, boundary_type='reflective')
# Instantiate Cells
fuel = openmc.Cell(cell_id=1, name='cell 1')
moderator = openmc.Cell(cell_id=2, name='cell 2')
fuel = openmc.Cell(fill=uo2, region=-fuel_or, name='fuel')
moderator = openmc.Cell(fill=water, region=+fuel_or & box, name='moderator')
# Use surface half-spaces to define regions
fuel.region = -fuel_or
moderator.region = +fuel_or & +left & -right & +bottom & -top
# Register Materials with Cells
fuel.fill = uo2
moderator.fill = water
# Instantiate Universe
root = openmc.Universe(universe_id=0, name='root universe')
# Register Cells with Universe
root.add_cells([fuel, moderator])
# Instantiate a Geometry, register the root Universe, and export to XML
geometry = openmc.Geometry(root)
# Create a geometry with the two cells and export to XML
geometry = openmc.Geometry([fuel, moderator])
geometry.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml file
###############################################################################
# Define problem settings
# Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.energy_mode = "multi-group"
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings = openmc.Settings()
settings.energy_mode = "multi-group"
settings.batches = 100
settings.inactive = 10
settings.particles = 1000
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-0.63, -0.63, -1, 0.63, 0.63, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()
lower_left = (-pitch/2, -pitch/2, -1)
upper_right = (pitch/2, pitch/2, 1)
uniform_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True)
settings.source = openmc.source.Source(space=uniform_dist)
settings.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml file
###############################################################################
# Define tallies
# Instantiate a tally mesh
mesh = openmc.RegularMesh(mesh_id=1)
mesh.dimension = [100, 100, 1]
mesh.lower_left = [-0.63, -0.63, -1.e50]
mesh.upper_right = [0.63, 0.63, 1.e50]
# Create a mesh that will be used for tallying
mesh = openmc.RegularMesh()
mesh.dimension = (100, 100)
mesh.lower_left = (-pitch/2, -pitch/2)
mesh.upper_right = (pitch/2, pitch/2)
# Instantiate some tally Filters
energy_filter = openmc.EnergyFilter([1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6,
1.0e6, 20.0e6])
# Create a mesh filter that can be used in a tally
mesh_filter = openmc.MeshFilter(mesh)
# Instantiate the Tally
tally = openmc.Tally(tally_id=1, name='tally 1')
tally.filters = [energy_filter, mesh_filter]
tally.scores = ['flux', 'fission', 'nu-fission']
# Now use the mesh filter in a tally and indicate what scores are desired
mesh_tally = openmc.Tally(name="Mesh tally")
mesh_tally.filters = [mesh_filter]
mesh_tally.scores = ['flux', 'fission', 'nu-fission']
# Instantiate a Tallies collection, register all Tallies, and export to XML
tallies_file = openmc.Tallies([tally])
tallies_file.export_to_xml()
# Let's also create a tally to get the flux energy spectrum. We start by
# creating an energy filter
e_min, e_max = 1e-5, 20.0e6
groups = 500
energies = np.logspace(log10(e_min), log10(e_max), groups + 1)
energy_filter = openmc.EnergyFilter(energies)
spectrum_tally = openmc.Tally(name="Flux spectrum")
spectrum_tally.filters = [energy_filter]
spectrum_tally.scores = ['flux']
# Instantiate a Tallies collection and export to XML
tallies = openmc.Tallies([mesh_tally, spectrum_tally])
tallies.export_to_xml()

121
examples/python/basic/build-xml.py
View file

@ -1,121 +0,0 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 15
inactive = 5
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml file
###############################################################################
# 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_element('H', 2.)
moderator.add_element('O', 1.)
moderator.add_s_alpha_beta('c_H_in_H2O')
fuel = openmc.Material(material_id=40, name='fuel')
fuel.set_density('g/cc', 4.5)
fuel.add_nuclide('U235', 1.)
# Instantiate a Materials collection and export to XML
materials_file = openmc.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')
# Use surface half-spaces to define regions
cell1.region = -surf2
cell2.region = -surf1
cell3.region = +surf1
cell4.region = +surf2 & -surf3
# 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, register the root Universe, and export to XML
geometry = openmc.Geometry(root)
geometry.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml file
###############################################################################
# Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-4., -4., -4., 4., 4., 4.]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml file
###############################################################################
# Instantiate some tally Filters
cell_filter = openmc.CellFilter(cell2)
energy_filter = openmc.EnergyFilter([0., 20.e6])
energyout_filter = openmc.EnergyoutFilter([0., 20.e6])
# Instantiate the first Tally
first_tally = openmc.Tally(tally_id=1, name='first tally')
first_tally.filters = [cell_filter]
scores = ['total', 'scatter', 'nu-scatter',
'absorption', 'fission', 'nu-fission']
first_tally.scores = scores
# Instantiate the second Tally
second_tally = openmc.Tally(tally_id=2, name='second tally')
second_tally.filters = [cell_filter, energy_filter]
second_tally.scores = scores
# Instantiate the third Tally
third_tally = openmc.Tally(tally_id=3, name='third tally')
third_tally.filters = [cell_filter, energy_filter, energyout_filter]
third_tally.scores = ['scatter', 'nu-scatter', 'nu-fission']
# Instantiate a Tallies collection and export to XML
tallies_file = openmc.Tallies((first_tally, second_tally, third_tally))
tallies_file.export_to_xml()

125
examples/python/boxes/build-xml.py
View file

@ -1,125 +0,0 @@
import numpy as np
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 15
inactive = 5
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate some Materials and register the appropriate Nuclides
fuel1 = openmc.Material(material_id=1, name='fuel')
fuel1.set_density('g/cc', 4.5)
fuel1.add_nuclide('U235', 1.)
fuel2 = openmc.Material(material_id=2, name='depleted fuel')
fuel2.set_density('g/cc', 4.5)
fuel2.add_nuclide('U238', 1.)
moderator = openmc.Material(material_id=3, name='moderator')
moderator.set_density('g/cc', 1.0)
moderator.add_element('H', 2.)
moderator.add_element('O', 1.)
moderator.add_s_alpha_beta('c_H_in_H2O')
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([fuel1, fuel2, moderator])
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml file
###############################################################################
# Instantiate planar surfaces
x1 = openmc.XPlane(surface_id=1, x0=-10)
x2 = openmc.XPlane(surface_id=2, x0=-7)
x3 = openmc.XPlane(surface_id=3, x0=-4)
x4 = openmc.XPlane(surface_id=4, x0=4)
x5 = openmc.XPlane(surface_id=5, x0=7)
x6 = openmc.XPlane(surface_id=6, x0=10)
y1 = openmc.YPlane(surface_id=11, y0=-10)
y2 = openmc.YPlane(surface_id=12, y0=-7)
y3 = openmc.YPlane(surface_id=13, y0=-4)
y4 = openmc.YPlane(surface_id=14, y0=4)
y5 = openmc.YPlane(surface_id=15, y0=7)
y6 = openmc.YPlane(surface_id=16, y0=10)
z1 = openmc.ZPlane(surface_id=21, z0=-10)
z2 = openmc.ZPlane(surface_id=22, z0=-7)
z3 = openmc.ZPlane(surface_id=23, z0=-4)
z4 = openmc.ZPlane(surface_id=24, z0=4)
z5 = openmc.ZPlane(surface_id=25, z0=7)
z6 = openmc.ZPlane(surface_id=26, z0=10)
# Set vacuum boundary conditions on outside
for surface in [x1, x6, y1, y6, z1, z6]:
surface.boundary_type = 'vacuum'
# Instantiate Cells
inner_box = openmc.Cell(cell_id=1, name='inner box')
middle_box = openmc.Cell(cell_id=2, name='middle box')
outer_box = openmc.Cell(cell_id=3, name='outer box')
# Use each set of six planes to create solid cube regions. We can then use these
# to create cubic shells.
inner_cube = +x3 & -x4 & +y3 & -y4 & +z3 & -z4
middle_cube = +x2 & -x5 & +y2 & -y5 & +z2 & -z5
outer_cube = +x1 & -x6 & +y1 & -y6 & +z1 & -z6
outside_inner_cube = -x3 | +x4 | -y3 | +y4 | -z3 | +z4
# Use surface half-spaces to define regions
inner_box.region = inner_cube
middle_box.region = middle_cube & outside_inner_cube
outer_box.region = outer_cube & ~middle_cube
# Register Materials with Cells
inner_box.fill = fuel1
middle_box.fill = fuel2
outer_box.fill = moderator
# Instantiate root universe
root = openmc.Universe(universe_id=0, name='root universe')
root.add_cells([inner_box, middle_box, outer_box])
# Instantiate a Geometry, register the root Universe, and export to XML
geometry = openmc.Geometry(root)
geometry.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
# Create an initial uniform spatial source distribution over fissionable zones
uniform_dist = openmc.stats.Box(*outer_cube.bounding_box, only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC plots.xml File
###############################################################################
plot = openmc.Plot(plot_id=1)
plot.origin = [0, 0, 0]
plot.width = [20, 20]
plot.pixels = [200, 200]
plot.color_by = 'cell'
# Instantiate a Plots collection and export to XML
plot_file = openmc.Plots([plot])
plot_file.export_to_xml()

138
examples/python/pincell/build-xml.py
View file

@ -1,138 +0,0 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 100
inactive = 10
particles = 1000
###############################################################################
# Exporting to OpenMC materials.xml file
###############################################################################
# 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_element('U', 1., enrichment=2.4)
uo2.add_element('O', 2.)
helium = openmc.Material(material_id=2, name='Helium for gap')
helium.set_density('g/cm3', 0.001598)
helium.add_element('He', 2.4044e-4)
zircaloy = openmc.Material(material_id=3, name='Zircaloy 4')
zircaloy.set_density('g/cm3', 6.55)
zircaloy.add_element('Sn', 0.014 , 'wo')
zircaloy.add_element('Fe', 0.00165, 'wo')
zircaloy.add_element('Cr', 0.001 , 'wo')
zircaloy.add_element('Zr', 0.98335, 'wo')
borated_water = openmc.Material(material_id=4, name='Borated water')
borated_water.set_density('g/cm3', 0.740582)
borated_water.add_element('B', 4.0e-5)
borated_water.add_element('H', 5.0e-2)
borated_water.add_element('O', 2.4e-2)
borated_water.add_s_alpha_beta('c_H_in_H2O')
# Instantiate a Materials collection and export to XML
materials_file = openmc.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')
# Use surface half-spaces to define regions
fuel.region = -fuel_or
gap.region = +fuel_or & -clad_ir
clad.region = +clad_ir & -clad_or
water.region = +clad_or & +left & -right & +bottom & -top
# 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, register the root Universe, and export to XML
geometry = openmc.Geometry(root)
geometry.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml file
###############################################################################
# Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
entropy_mesh = openmc.RegularMesh()
entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]
entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50]
entropy_mesh.dimension = [10, 10, 1]
settings_file.entropy_mesh = entropy_mesh
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml file
###############################################################################
# Instantiate a tally mesh
mesh = openmc.RegularMesh()
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.EnergyFilter([0., 4., 20.e6])
mesh_filter = openmc.MeshFilter(mesh)
# Instantiate the Tally
tally = openmc.Tally(tally_id=1, name='tally 1')
tally.filters = [energy_filter, mesh_filter]
tally.scores = ['flux', 'fission', 'nu-fission']
# Instantiate a Tallies collection and export to XML
tallies_file = openmc.Tallies([tally])
tallies_file.export_to_xml()

82
examples/python/reflective/build-xml.py
View file

@ -1,82 +0,0 @@
import numpy as np
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 500
inactive = 10
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml file
###############################################################################
# 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 Materials collection and export to XML
materials_file = openmc.Materials([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')
# Use surface half-spaces to define region
cell.region = +surf1 & -surf2 & +surf3 & -surf4 & +surf5 & -surf6
# 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, register the root Universe, and export to XML
geometry = openmc.Geometry(root)
geometry.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml file
###############################################################################
# Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
# Create an initial uniform spatial source distribution over fissionable zones
uniform_dist = openmc.stats.Box(*cell.region.bounding_box,
only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()

15
examples/xml/basic/geometry.xml
View file

@ -1,15 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Definition of Cells -->
<cell id="1" universe="0" fill="37" region="-2" />
<cell id="100" universe="37" material="40" region="-1" />
<cell id="101" universe="37" material="41" region="1" />
<cell id="2" universe="0" material="41" region="2 -3" />
<!-- Defition of Surfaces -->
<surface id="1" type="z-cylinder" coeffs="0 0 7" />
<surface id="2" type="z-cylinder" coeffs="0 0 9" />
<surface id="3" type="z-cylinder" coeffs="0 0 11" boundary="vacuum" />
</geometry>

16
examples/xml/basic/materials.xml
View file

@ -1,16 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
<material id="41">
<density value="1.0" units="g/cc" />
<nuclide name="H1" ao="2.0" />
<nuclide name="O16" ao="1.0" />
<sab name="c_H_in_H2O"/>
</material>
</materials>

16
examples/xml/basic/settings.xml
View file

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

31
examples/xml/basic/tallies.xml
View file

@ -1,31 +0,0 @@
<?xml version="1.0"?>
<tallies>
<filter id="1" type="cell">
<bins>100</bins>
</filter>
<filter id="2" type="energy">
<bins>0 20.0e6</bins>
</filter>
<filter id="3" type="energyout">
<bins>0 20.0e6</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>total scatter nu-scatter absorption fission nu-fission</scores>
</tally>
<tally id="2">
<filters>1 2</filters>
<scores>total scatter nu-scatter absorption fission nu-fission</scores>
</tally>
<tally id="3">
<filters>1 2 3</filters>
<scores>scatter nu-scatter nu-fission</scores>
</tally>
</tallies>

39
examples/xml/boxes/geometry.xml
View file

@ -1,39 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!--
This example consists of three nested boxes, and is meant to show how to
use Boolean operators to construct complex cell regions.
-->
<surface id="1" type="x-plane" coeffs="-10" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="-7" />
<surface id="3" type="x-plane" coeffs="-4" />
<surface id="4" type="x-plane" coeffs="4" />
<surface id="5" type="x-plane" coeffs="7" />
<surface id="6" type="x-plane" coeffs="10" boundary="vacuum" />
<surface id="11" type="y-plane" coeffs="-10" boundary="vacuum" />
<surface id="12" type="y-plane" coeffs="-7" />
<surface id="13" type="y-plane" coeffs="-4" />
<surface id="14" type="y-plane" coeffs="4" />
<surface id="15" type="y-plane" coeffs="7" />
<surface id="16" type="y-plane" coeffs="10" boundary="vacuum" />
<surface id="21" type="z-plane" coeffs="-10" boundary="vacuum" />
<surface id="22" type="z-plane" coeffs="-7" />
<surface id="23" type="z-plane" coeffs="-4" />
<surface id="24" type="z-plane" coeffs="4" />
<surface id="25" type="z-plane" coeffs="7" />
<surface id="26" type="z-plane" coeffs="10" boundary="vacuum" />
<!-- Innermost cube -->
<cell id="1" material="1" region="3 -4 13 -14 23 -24" />
<!-- Middle cubic shell -->
<cell id="2" material="2" region="2 -5 12 -15 22 -25 (-3 | 4 | -13 | 14 | -23 | 24)" />
<!-- Outermost cubic shell -->
<cell id="3" material="3" region="1 -6 11 -16 21 -26 ~(2 -5 12 -15 22 -25)" />
</geometry>

21
examples/xml/boxes/materials.xml
View file

@ -1,21 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
<material id="2">
<density value="4.5" units="g/cc" />
<nuclide name="U238" ao="1.0" />
</material>
<material id="3">
<density value="1.0" units="g/cc" />
<nuclide name="O16" ao="1.0" />
<nuclide name="H1" ao="2.0" />
<sab name="c_H_in_H2O" />
</material>
</materials>

9
examples/xml/boxes/plots.xml
View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<plots>
<plot id="1" type="slice">
<color_by>cell</color_by>
<origin>0. 0. 0.</origin>
<width>20. 20.</width>
<pixels>200 200</pixels>
</plot>
</plots>

15
examples/xml/boxes/settings.xml
View file

@ -1,15 +0,0 @@
<?xml version="1.0"?>
<settings>
<!-- Parameters for k-eigenvalue calculation -->
<run_mode>eigenvalue</run_mode>
<batches>15</batches>
<inactive>5</inactive>
<particles>10000</particles>
<!-- Starting source -->
<source>
<space type="box" parameters="-10. -10. -10. 10. 10. 10." />
</source>
</settings>

15
examples/xml/custom_source/geometry.xml
View file

@ -1,15 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Definition of Cells -->
<cell id="1" universe="0" fill="37" region="-2" />
<cell id="100" universe="37" material="40" region="-1" />
<cell id="101" universe="37" material="41" region="1" />
<cell id="2" universe="0" material="41" region="2 -3" />
<!-- Defition of Surfaces -->
<surface id="1" type="z-cylinder" coeffs="0 0 7" />
<surface id="2" type="z-cylinder" coeffs="0 0 9" />
<surface id="3" type="z-cylinder" coeffs="0 0 11" boundary="vacuum" />
</geometry>

16
examples/xml/custom_source/materials.xml
View file

@ -1,16 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
<material id="41">
<density value="1.0" units="g/cc" />
<nuclide name="H1" ao="2.0" />
<nuclide name="O16" ao="1.0" />
<sab name="c_H_in_H2O"/>
</material>
</materials>

14
examples/xml/custom_source/settings.xml
View file

@ -1,14 +0,0 @@
<?xml version="1.0"?>
<settings>
<run_mode>fixed source</run_mode>
<batches>10</batches>
<inactive>0</inactive>
<particles>100000</particles>
<!-- Starting source -->
<source>
<library>build/libsource.so</library>
</source>
</settings>

17
examples/xml/custom_source/tallies.xml
View file

@ -1,17 +0,0 @@
<?xml version="1.0"?>
<tallies>
<filter id="1" type="cell">
<bins>100</bins>
</filter>
<filter id="2" type="energy">
<bins>0 20.0e6</bins>
</filter>
<tally id="3">
<filters>1 2 </filters>
<scores>flux</scores>
</tally>
</tallies>

43
examples/xml/lattice/nested/geometry.xml
View file

@ -1,43 +0,0 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="6" region="1 -2 3 -4" />
<cell id="2" universe="5" fill="4" region="1 -2 3 -4" />
<cell id="101" universe="1" material="1" region="-5" />
<cell id="102" universe="1" material="2" region="5" />
<cell id="201" universe="2" material="1" region="-6" />
<cell id="202" universe="2" material="2" region="6" />
<cell id="301" universe="3" material="1" region="-7" />
<cell id="302" universe="3" material="2" region="7" />
<!-- 4 x 4 assembly -->
<lattice id="4">
<dimension>2 2</dimension>
<lower_left>-1.0 -1.0</lower_left>
<pitch>1.0 1.0</pitch>
<universes>
1 2
2 3
</universes>
</lattice>
<!-- 4 x 4 core -->
<lattice id="6">
<dimension>2 2</dimension>
<lower_left>-2.0 -2.0</lower_left>
<pitch>2.0 2.0</pitch>
<universes>
5 5
5 5
</universes>
</lattice>
<surface id="1" type="x-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="2.0" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="2.0" boundary="vacuum" />
<surface id="5" type="z-cylinder" coeffs="0.0 0.0 0.4" />
<surface id="6" type="z-cylinder" coeffs="0.0 0.0 0.3" />
<surface id="7" type="z-cylinder" coeffs="0.0 0.0 0.2" />
</geometry>

17
examples/xml/lattice/nested/materials.xml
View file

@ -1,17 +0,0 @@
<?xml version="1.0"?>
<materials>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
<material id="2">
<density value="1.0" units="g/cc" />
<nuclide name="H1" ao="2.0" />
<nuclide name="O16" ao="1.0" />
<sab name="c_H_in_H2O" />
</material>
</materials>

10
examples/xml/lattice/nested/plots.xml
View file

@ -1,10 +0,0 @@
<?xml version="1.0"?>
<plots>
<plot id="1" color_by="material">
<origin>0. 0. 0.</origin>
<width>4.0 4.0</width>
<pixels>400 400</pixels>
</plot>
</plots>

17
examples/xml/lattice/nested/settings.xml
View file

@ -1,17 +0,0 @@
<?xml version="1.0"?>
<settings>
<!-- Parameters for k-eigenvalue calculation -->
<run_mode>eigenvalue</run_mode>
<batches>20</batches>
<inactive>10</inactive>
<particles>10000</particles>
<!-- Starting source -->
<source>
<space type="box">
<parameters>-1 -1 -1 1 1 1</parameters>
</space>
</source>
</settings>

20
examples/xml/lattice/nested/tallies.xml
View file

@ -1,20 +0,0 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1">
<type>regular</type>
<dimension>4 4</dimension>
<lower_left>-2.0 -2.0</lower_left>
<width>1.0 1.0</width>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>total</scores>
</tally>
</tallies>

32
examples/xml/lattice/simple/geometry.xml
View file

@ -1,32 +0,0 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="5" region="1 -2 3 -4" />
<cell id="101" universe="1" material="1" region="-5" />
<cell id="102" universe="1" material="2" region="5" />
<cell id="201" universe="2" material="1" region="-6" />
<cell id="202" universe="2" material="2" region="6" />
<cell id="301" universe="3" material="1" region="-7" />
<cell id="302" universe="3" material="2" region="7" />
<lattice id="5">
<dimension>4 4</dimension>
<lower_left>-2.0 -2.0</lower_left>
<pitch>1.0 1.0</pitch>
<universes>
1 2 1 2
2 3 2 3
1 2 1 2
2 3 2 3
</universes>
</lattice>
<surface id="1" type="x-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="2.0" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="2.0" boundary="vacuum" />
<surface id="5" type="z-cylinder" coeffs="0.0 0.0 0.4" />
<surface id="6" type="z-cylinder" coeffs="0.0 0.0 0.3" />
<surface id="7" type="z-cylinder" coeffs="0.0 0.0 0.2" />
</geometry>

17
examples/xml/lattice/simple/materials.xml
View file

@ -1,17 +0,0 @@
<?xml version="1.0"?>
<materials>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
<material id="2">
<density value="1.0" units="g/cc" />
<nuclide name="H1" ao="2.0" />
<nuclide name="O16" ao="1.0" />
<sab name="c_H_in_H2O" />
</material>
</materials>

10
examples/xml/lattice/simple/plots.xml
View file

@ -1,10 +0,0 @@
<?xml version="1.0"?>
<plots>
<plot id="1" color_by="material">
<origin>0. 0. 0.</origin>
<width>4.0 4.0</width>
<pixels>400 400</pixels>
</plot>
</plots>

17
examples/xml/lattice/simple/settings.xml
View file

@ -1,17 +0,0 @@
<?xml version="1.0"?>
<settings>
<!-- Parameters for k-eigenvalue calculation -->
<run_mode>eigenvalue</run_mode>
<batches>20</batches>
<inactive>10</inactive>
<particles>10000</particles>
<!-- Starting source -->
<source>
<space type="box">
<parameters>-1 -1 -1 1 1 1</parameters>
</space>
</source>
</settings>

20
examples/xml/lattice/simple/tallies.xml
View file

@ -1,20 +0,0 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1">
<type>regular</type>
<dimension>4 4</dimension>
<lower_left>-2.0 -2.0</lower_left>
<width>1.0 1.0</width>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>total</scores>
</tally>
</tallies>

27
examples/xml/pincell/geometry.xml
View file

@ -1,27 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!--
This is a simple pin cell model based on dimensions from the MIT BEAVRS
(Benchmarking for Evaluation and Validation of Reactor Simulations)
benchmark.
-->
<!-- Surfaces for fuel, gap, cladding. Dimensions from Figure 2 in BEAVRS -->
<surface id="1" type="z-cylinder" coeffs="0. 0. 0.39218" /> <!-- Fuel OR -->
<surface id="2" type="z-cylinder" coeffs="0. 0. 0.40005" /> <!-- Clad IR -->
<surface id="3" type="z-cylinder" coeffs="0. 0. 0.45720" /> <!-- Clad OR -->
<!-- Reflective surfaces on outside of pin-cell. The lattice pitch is 1.25984
cm (taken from Table 2 in BEAVRS). -->
<surface id="4" type="x-plane" coeffs="-0.62992" boundary="reflective" />
<surface id="5" type="x-plane" coeffs=" 0.62992" boundary="reflective" />
<surface id="6" type="y-plane" coeffs="-0.62992" boundary="reflective" />
<surface id="7" type="y-plane" coeffs=" 0.62992" boundary="reflective" />
<cell id="1" material="1" region=" -1" /> <!-- UO2 Fuel -->
<cell id="2" material="2" region="1 -2" /> <!-- Helium gap -->
<cell id="3" material="3" region="2 -3" /> <!-- Zircaloy cladding -->
<cell id="4" material="4" region="3 4 -5 6 -7" /> <!-- Borated water -->
</geometry>

67
examples/xml/pincell/materials.xml
View file

@ -1,67 +0,0 @@
<?xml version="1.0"?>
<materials>
<!--
Since O-18 is not present in ENDF/B-VII, it was necessary to combine the
atom densities for O-17 and O-18 in any materials containing Oxygen.
-->
<!-- UO2 fuel at 2.4 wt% enrichment -->
<material id="1">
<density value="10.29769" units="g/cm3" />
<nuclide name="U234" ao="4.4843e-06" />
<nuclide name="U235" ao="5.5815e-04" />
<nuclide name="U238" ao="2.2408e-02" />
<nuclide name="O16" ao="4.5829e-02" />
<nuclide name="O17" ao="1.1164e-04" />
</material>
<!-- Helium for gap -->
<material id="2">
<density value="0.001598" units="g/cm3" />
<nuclide name="He4" ao="2.4044e-04" />
</material>
<!-- Zircaloy 4 -->
<material id="3">
<density value="6.55" units="g/cm3" />
<nuclide name="O16" ao="3.0743e-04" />
<nuclide name="O17" ao="7.4887e-07" />
<nuclide name="Cr50" ao="3.2962e-06" />
<nuclide name="Cr52" ao="6.3564e-05" />
<nuclide name="Cr53" ao="7.2076e-06" />
<nuclide name="Cr54" ao="1.7941e-06" />
<nuclide name="Fe54" ao="8.6699e-06" />
<nuclide name="Fe56" ao="1.3610e-04" />
<nuclide name="Fe57" ao="3.1431e-06" />
<nuclide name="Fe58" ao="4.1829e-07" />
<nuclide name="Zr90" ao="2.1827e-02" />
<nuclide name="Zr91" ao="4.7600e-03" />
<nuclide name="Zr92" ao="7.2758e-03" />
<nuclide name="Zr94" ao="7.3734e-03" />
<nuclide name="Zr96" ao="1.1879e-03" />
<nuclide name="Sn112" ao="4.6735e-06" />
<nuclide name="Sn114" ao="3.1799e-06" />
<nuclide name="Sn115" ao="1.6381e-06" />
<nuclide name="Sn116" ao="7.0055e-05" />
<nuclide name="Sn117" ao="3.7003e-05" />
<nuclide name="Sn118" ao="1.1669e-04" />
<nuclide name="Sn119" ao="4.1387e-05" />
<nuclide name="Sn120" ao="1.5697e-04" />
<nuclide name="Sn122" ao="2.2308e-05" />
<nuclide name="Sn124" ao="2.7897e-05" />
</material>
<!-- Borated water at 975 ppm -->
<material id="4">
<density value="0.740582" units="g/cm3" />
<nuclide name="B10" ao="8.0042e-06" />
<nuclide name="B11" ao="3.2218e-05" />
<nuclide name="H1" ao="4.9457e-02" />
<nuclide name="H2" ao="7.4196e-06" />
<nuclide name="O16" ao="2.4672e-02" />
<nuclide name="O17" ao="6.0099e-05" />
<sab name="c_H_in_H2O" />
</material>
</materials>

32
examples/xml/pincell/settings.xml
View file

@ -1,32 +0,0 @@
<?xml version="1.0"?>
<settings>
<!-- Define how many particles to run and for how many batches -->
<run_mode>eigenvalue</run_mode>
<batches>100</batches>
<inactive>10</inactive>
<particles>1000</particles>
<!-- The starting source is a uniform distribution over the entire pin
cell. Note that since this is effectively a 2D model, the z coordinates
are inconsequential -->
<source>
<space type="box">
<parameters>
-0.62992 -0.62992 -1.
0.62992 0.62992 1.
</parameters>
</space>
</source>
<!-- To assess convergence of the source distribution, we need to define the
bounds for a mesh over which the Shannon entropy should be
calculated. The extent in the z direction is made arbitrarily large. -->
<mesh id="1">
<lower_left>-0.39218 -0.39218 -1.e50</lower_left>
<upper_right>0.39218 0.39218 1.e50</upper_right>
<dimension>10 10 1</dimension>
</mesh>
<entropy_mesh>1</entropy_mesh>
</settings>

23
examples/xml/pincell/tallies.xml
View file

@ -1,23 +0,0 @@
<?xml version="1.0"?>
<tallies>
<mesh id="2" type="regular">
<dimension>100 100 1</dimension>
<lower_left>-0.62992 -0.62992 -1.e50</lower_left>
<upper_right>0.62992 0.62992 1.e50</upper_right>
</mesh>
<filter id="1" type="mesh">
<bins>2</bins>
</filter>
<filter id="2" type="energy">
<bins>0. 4. 20.0e6</bins>
</filter>
<tally id="1">
<filters>1 2</filters>
<scores>flux fission nu-fission</scores>
</tally>
</tallies>

10
examples/xml/pincell_multigroup/geometry.xml
View file

@ -1,10 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" name="cell 1" region="-1" universe="0" />
<cell id="2" material="2" name="cell 2" region="1 4 -5 6 -7" universe="0" />
<surface coeffs="0 0 0.54" id="1" name="Fuel OR" type="z-cylinder" />
<surface boundary="reflective" coeffs="-0.63" id="4" name="left" type="x-plane" />
<surface boundary="reflective" coeffs="0.63" id="5" name="right" type="x-plane" />
<surface boundary="reflective" coeffs="-0.63" id="6" name="bottom" type="y-plane" />
<surface boundary="reflective" coeffs="0.63" id="7" name="top" type="y-plane" />
</geometry>

12
examples/xml/pincell_multigroup/materials.xml
View file

@ -1,12 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>./mgxs.h5</cross_sections>
<material id="1" name="UO2 fuel">
<density units="macro" value="1.0" />
<macroscopic name="UO2" />
</material>
<material id="2" name="Water">
<density units="macro" value="1.0" />
<macroscopic name="LWTR" />
</material>
</materials>

BIN
examples/xml/pincell_multigroup/mgxs.h5
View file

Binary file not shown.

28
examples/xml/pincell_multigroup/plots.xml
View file

@ -1,28 +0,0 @@
<?xml version="1.0"?>
<plots>
<plot>
<id>1</id>
<filename>mat</filename>
<color_by>material</color_by>
<origin>0 0 0</origin>
<width>1.26 1.26</width>
<type>slice</type>
<pixels>1000 1000 </pixels>
<color id="1" rgb="255 0 0" />
<color id="2" rgb="0 0 0" />
<color id="3" rgb="0 255 0" />
<color id="4" rgb="0 0 255" />
</plot>
<plot>
<id>2</id>
<filename>cell</filename>
<color_by>cell</color_by>
<origin>0 0 0</origin>
<width>1.26 1.26</width>
<type>slice</type>
<pixels>1000 1000 </pixels>
</plot>
</plots>

13
examples/xml/pincell_multigroup/settings.xml
View file

@ -1,13 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>100</batches>
<inactive>10</inactive>
<source strength="1.0">
<space type="box">
<parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
</settings>

18
examples/xml/pincell_multigroup/tallies.xml
View file

@ -1,18 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" type="regular">
<dimension>100 100 1</dimension>
<lower_left>-0.63 -0.63 -1e+50</lower_left>
<upper_right>0.63 0.63 1e+50</upper_right>
</mesh>
<filter id="1" type="energy">
<bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
</filter>
<filter id="2" type="mesh">
<bins>1</bins>
</filter>
<tally id="1" name="tally 1">
<filters>1 2</filters>
<scores>flux fission nu-fission</scores>
</tally>
</tallies>

19
examples/xml/reflective/geometry.xml
View file

@ -1,19 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Definition of Cells -->
<cell id="1">
<universe>0</universe>
<material>1</material>
<region>1 -2 3 -4 5 -6</region>
</cell>
<!-- Defition of Surfaces -->
<surface id="1" type="x-plane" coeffs="-1" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="1" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-1" boundary="reflective" />
<surface id="4" type="y-plane" coeffs="1" boundary="reflective" />
<surface id="5" type="z-plane" coeffs="-1" boundary="reflective" />
<surface id="6" type="z-plane" coeffs="1" boundary="reflective" />
</geometry>

9
examples/xml/reflective/materials.xml
View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
</materials>

17
examples/xml/reflective/settings.xml
View file

@ -1,17 +0,0 @@
<?xml version="1.0"?>
<settings>
<!-- Parameters for k-eigenvalue calculation -->
<run_mode>eigenvalue</run_mode>
<batches>500</batches>
<inactive>10</inactive>
<particles>10000</particles>
<!-- Starting source -->
<source>
<space type="box">
<parameters>-1 -1 -1 1 1 1</parameters>
</space>
</source>
</settings>

5
include/openmc/distribution_multi.h
View file

@ -43,6 +43,11 @@ public:
//! \param seed Pseudorandom number seed pointer
//! \return Direction sampled
Direction sample(uint64_t* seed) const;
// Observing pointers
Distribution* mu() const { return mu_.get(); }
Distribution* phi() const { return phi_.get(); }
private:
UPtrDist mu_; //!< Distribution of polar angle
UPtrDist phi_; //!< Distribution of azimuthal angle

19
include/openmc/distribution_spatial.h
View file

@ -32,6 +32,11 @@ public:
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const;
// Observer pointers
Distribution* x() const { return x_.get(); }
Distribution* y() const { return x_.get(); }
Distribution* z() const { return x_.get(); }
private:
UPtrDist x_; //!< Distribution of x coordinates
UPtrDist y_; //!< Distribution of y coordinates
@ -50,6 +55,11 @@ public:
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const;
Distribution* r() const { return r_.get(); }
Distribution* phi() const { return phi_.get(); }
Distribution* z() const { return z_.get(); }
Position origin() const { return origin_; }
private:
UPtrDist r_; //!< Distribution of r coordinates
UPtrDist phi_; //!< Distribution of phi coordinates
@ -70,6 +80,11 @@ public:
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const;
Distribution* r() const { return r_.get(); }
Distribution* theta() const { return theta_.get(); }
Distribution* phi() const { return phi_.get(); }
Position origin () const { return origin_; }
private:
UPtrDist r_; //!< Distribution of r coordinates
UPtrDist theta_; //!< Distribution of theta coordinates
@ -92,6 +107,8 @@ public:
// Properties
bool only_fissionable() const { return only_fissionable_; }
Position lower_left() const { return lower_left_; }
Position upper_right() const { return upper_right_; }
private:
Position lower_left_; //!< Lower-left coordinates of box
Position upper_right_; //!< Upper-right coordinates of box
@ -112,6 +129,8 @@ public:
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const;
Position r() const { return r_; }
private:
Position r_; //!< Single position at which sites are generated
};

242
include/openmc/mesh.h
View file

@ -15,6 +15,13 @@
#include "openmc/particle.h"
#include "openmc/position.h"
#ifdef DAGMC
#include "moab/Core.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/Matrix3.hpp"
#include "moab/GeomUtil.hpp"
#endif
namespace openmc {
//==============================================================================
@ -30,6 +37,7 @@ extern std::unordered_map<int32_t, int32_t> mesh_map;
} // namespace model
class Mesh
{
public:
@ -61,6 +69,42 @@ public:
//! \return Mesh bin
virtual int get_bin(Position r) const = 0;
//! Get the number of mesh cells.
virtual int n_bins() const = 0;
//! Get the number of mesh cell surfaces.
virtual int n_surface_bins() const = 0;
//! Write mesh data to an HDF5 group
//
//! \param[in] group HDF5 group
virtual void to_hdf5(hid_t group) const = 0;
//! Find the mesh lines that intersect an axis-aligned slice plot
//
//! \param[in] plot_ll The lower-left coordinates of the slice plot.
//! \param[in] plot_ur The upper-right coordinates of the slice plot.
//! \return A pair of vectors indicating where the mesh lines lie along each
//! of the plot's axes. For example an xy-slice plot will get back a vector
//! of x-coordinates and another of y-coordinates. These vectors may be
//! empty for low-dimensional meshes.
virtual std::pair<std::vector<double>, std::vector<double>>
plot(Position plot_ll, Position plot_ur) const = 0;
//! Get a label for the mesh bin
virtual std::string bin_label(int bin) const = 0;
// Data members
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions
};
class StructuredMesh : public Mesh {
public:
StructuredMesh() = default;
StructuredMesh(pugi::xml_node node) : Mesh {node} {};
virtual ~StructuredMesh() = default;
//! Get bin given mesh indices
//
//! \param[in] Array of mesh indices
@ -80,32 +124,10 @@ public:
//! \param[out] ijk Mesh indices
virtual void get_indices_from_bin(int bin, int* ijk) const = 0;
//! Get the number of mesh cells.
virtual int n_bins() const = 0;
//! Get the number of mesh cell surfaces.
virtual int n_surface_bins() const = 0;
//! Find the mesh lines that intersect an axis-aligned slice plot
//
//! \param[in] plot_ll The lower-left coordinates of the slice plot.
//! \param[in] plot_ur The upper-right coordinates of the slice plot.
//! \return A pair of vectors indicating where the mesh lines lie along each
//! of the plot's axes. For example an xy-slice plot will get back a vector
//! of x-coordinates and another of y-coordinates. These vectors may be
//! empty for low-dimensional meshes.
virtual std::pair<std::vector<double>, std::vector<double>>
plot(Position plot_ll, Position plot_ur) const = 0;
//! Write mesh data to an HDF5 group
//
//! \param[in] group HDF5 group
virtual void to_hdf5(hid_t group) const = 0;
//! Get a label for the mesh bin
std::string bin_label(int bin) const override;
// Data members
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions
xt::xtensor<double, 1> lower_left_; //!< Lower-left coordinates of mesh
xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
};
@ -114,7 +136,7 @@ public:
//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes
//==============================================================================
class RegularMesh : public Mesh
class RegularMesh : public StructuredMesh
{
public:
// Constructors
@ -176,7 +198,8 @@ private:
bool intersects_3d(Position& r0, Position r1, int* ijk) const;
};
class RectilinearMesh : public Mesh
class RectilinearMesh : public StructuredMesh
{
public:
// Constructors
@ -218,13 +241,176 @@ public:
bool intersects(Position& r0, Position r1, int* ijk) const;
// Data members
xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
private:
std::vector<std::vector<double>> grid_;
};
#ifdef DAGMC
class UnstructuredMesh : public Mesh {
public:
UnstructuredMesh() = default;
UnstructuredMesh(pugi::xml_node);
~UnstructuredMesh() = default;
void bins_crossed(const Particle* p,
std::vector<int>& bins,
std::vector<double>& lengths) const override;
std::pair<std::vector<double>, std::vector<double>>
plot(Position plot_ll, Position plot_ur) const override;
//! Determine which surface bins were crossed by a particle.
//
//! \param[in] p Particle to check
//! \param[out] bins Surface bins that were crossed
void surface_bins_crossed(const Particle* p, std::vector<int>& bins) const;
//! Write mesh data to an HDF5 group.
//
//! \param[in] group HDF5 group
void to_hdf5(hid_t group) const;
//! Get bin at a given position.
//
//! \param[in] r Position to get bin for
//! \return Mesh bin
int get_bin(Position r) const;
int n_bins() const override;
int n_surface_bins() const override;
//! Retrieve a centroid for the mesh cell
//
// \param[in] tet MOAB EntityHandle of the tetrahedron
// \return The centroid of the element
Position centroid(moab::EntityHandle tet) const;
//! Return a string represntation of the mesh bin
//
//! \param[in] bin Mesh bin to generate a label for
std::string bin_label(int bin) const override;
//! Add a score to the mesh instance
void add_score(std::string score) const;
//! Set data for a score
void set_score_data(const std::string& score,
std::vector<double> values,
std::vector<double> std_dev) const;
//! Write the mesh with any current tally data
void write(std::string base_filename) const;
std::string filename_; //!< Path to unstructured mesh file
private:
//! Find all intersections with faces of the mesh.
//
//! \param[in] start Staring location
//! \param[in] dir Normalized particle direction
//! \param[in] track_len length of particle track
//! \param[out] Mesh intersections
void
intersect_track(const moab::CartVect& start,
const moab::CartVect& dir,
double track_len,
std::vector<double>& hits) const;
//! Calculate the volume for a given tetrahedron handle.
//
// \param[in] tet MOAB EntityHandle of the tetrahedron
double tet_volume(moab::EntityHandle tet) const;
//! Find the tetrahedron for the given location if
//! one exists
//
//! \param[in]
//! \return MOAB EntityHandle of tet
moab::EntityHandle get_tet(const Position& r) const;
//! Return the containing tet given a position
moab::EntityHandle get_tet(const moab::CartVect& r) const {
return get_tet(Position(r[0], r[1], r[2]));
};
//! Check for point containment within a tet; uses
//! pre-computed barycentric data.
//
//! \param[in] r Position to check
//! \param[in] MOAB terahedron to check
//! \return True if r is inside, False if r is outside
bool point_in_tet(const moab::CartVect& r,
moab::EntityHandle tet) const;
//! Compute barycentric coordinate data for all tetrahedra
//! in the mesh.
//
//! \param[in] tets MOAB Range of tetrahedral elements
void compute_barycentric_data(const moab::Range& tets);
//! Translate a MOAB EntityHandle to its corresponding bin.
//
//! \param[in] eh MOAB EntityHandle to translate
//! \return Mesh bin
int get_bin_from_ent_handle(moab::EntityHandle eh) const;
//! Translate a bin to its corresponding MOAB EntityHandle
//! for the tetrahedron representing that bin.
//
//! \param[in] bin Bin value to translate
//! \return MOAB EntityHandle of tet
moab::EntityHandle get_ent_handle_from_bin(int bin) const;
//! Get the bin for a given mesh cell index
//
//! \param[in] idx Index of the mesh cell.
//! \return Mesh bin
int get_bin_from_index(int idx) const;
//! Get the mesh cell index for a given position
//
//! \param[in] r Position to get index for
//! \param[in,out] in_mesh Whether position is in the mesh
int get_index(const Position& r, bool* in_mesh) const;
//! Get the mesh cell index from a bin
//
//! \param[in] bin Bin to get the index for
//! \return Index of the bin
int get_index_from_bin(int bin) const;
//! Build a KDTree for all tetrahedra in the mesh. All
//! triangles representing 2D faces of the mesh are
//! added to the tree as well.
//
//! \param[in] all_tets MOAB Range of tetrahedra for the tree
void build_kdtree(const moab::Range& all_tets);
//! Get the tags for a score from the mesh instance
//! or create them if they are not there
//
//! \param[in] score Name of the score
//! \return The MOAB value and error tag handles, respectively
std::pair<moab::Tag, moab::Tag>
get_score_tags(std::string score) const;
// data members
moab::Range ehs_; //!< Range of tetrahedra EntityHandle's in the mesh
moab::EntityHandle tetset_; //!< EntitySet containing all tetrahedra
moab::EntityHandle kdtree_root_; //!< Root of the MOAB KDTree
std::unique_ptr<moab::Interface> mbi_; //!< MOAB instance
std::unique_ptr<moab::AdaptiveKDTree> kdtree_; //!< MOAB KDTree instance
std::vector<moab::Matrix3> baryc_data_; //!< Barycentric data for tetrahedra
};
#endif
//==============================================================================
// Non-member functions
//==============================================================================
@ -239,6 +425,8 @@ void read_meshes(pugi::xml_node root);
//! \param[in] group HDF5 group
void meshes_to_hdf5(hid_t group);
RegularMesh* get_regular_mesh(int32_t index);
void free_memory_mesh();
} // namespace openmc

7
include/openmc/source.h
View file

@ -43,7 +43,14 @@ public:
Particle::Bank sample(uint64_t* seed) const;
// Properties
Particle::Type particle_type() const { return particle_; }
double strength() const { return strength_; }
// Make observing pointers available
SpatialDistribution* space() const { return space_.get(); }
UnitSphereDistribution* angle() const { return angle_.get(); }
Distribution* energy() const { return energy_.get(); }
private:
Particle::Type particle_ {Particle::Type::neutron}; //!< Type of particle emitted
double strength_ {1.0}; //!< Source strength

4
include/openmc/state_point.h
View file

@ -16,5 +16,9 @@ void read_source_bank(hid_t group_id);
void write_tally_results_nr(hid_t file_id);
void restart_set_keff();
#ifdef DAGMC
void write_unstructured_mesh_results();
#endif
} // namespace openmc
#endif // OPENMC_STATE_POINT_H

3
include/openmc/tallies/tally.h
View file

@ -72,6 +72,9 @@ public:
void accumulate();
//! A string representing the i-th score on this tally
std::string score_name(int score_idx) const;
//----------------------------------------------------------------------------
// Major public data members.

70
openmc/_utils.py
View file

@ -1,70 +0,0 @@
import hashlib
import os.path
from pathlib import Path
from urllib.parse import urlparse
from urllib.request import urlopen, Request
_BLOCK_SIZE = 16384
def download(url, checksum=None, as_browser=False, **kwargs):
"""Download file from a URL
Parameters
----------
url : str
URL from which to download
checksum : str or None
MD5 checksum to check against
as_browser : bool
Change User-Agent header to appear as a browser
kwargs : dict
Keyword arguments passed to :func:urllib.request.urlopen
Returns
-------
basename : str
Name of file written locally
"""
if as_browser:
page = Request(url, headers={'User-Agent': 'Mozilla/5.0'})
else:
page = url
with urlopen(page, **kwargs) as response:
# Get file size from header
file_size = response.length
# Check if file already downloaded
basename = Path(urlparse(url).path).name
if os.path.exists(basename):
if os.path.getsize(basename) == file_size:
print('Skipping {}, already downloaded'.format(basename))
return basename
# Copy file to disk in chunks
print('Downloading {}... '.format(basename), end='')
downloaded = 0
with open(basename, 'wb') as fh:
while True:
chunk = response.read(_BLOCK_SIZE)
if not chunk:
break
fh.write(chunk)
downloaded += len(chunk)
status = '{:10} [{:3.2f}%]'.format(
downloaded, downloaded * 100. / file_size)
print(status + '\b'*len(status), end='', flush=True)
print('')
if checksum is not None:
downloadsum = hashlib.md5(open(basename, 'rb').read()).hexdigest()
if downloadsum != checksum:
raise IOError("MD5 checksum for {} does not match. If this is "
"your first time receiving this message, please "
"re-run the script. Otherwise, please contact "
"OpenMC developers by emailing "
"openmc-users@googlegroups.com.".format(basename))
return basename

8
openmc/_xml.py
View file

@ -1,26 +1,20 @@
def clean_indentation(element, level=0, spaces_per_level=2):
"""
copy and paste from http://effbot.org/zone/element-lib.htm#prettyprint
copy and paste from https://effbot.org/zone/element-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*spaces_per_level*" "
if len(element):
if not element.text or not element.text.strip():
element.text = i + spaces_per_level*" "
if not element.tail or not element.tail.strip():
element.tail = i
for sub_element in element:
clean_indentation(sub_element, level+1, spaces_per_level)
if not sub_element.tail or not sub_element.tail.strip():
sub_element.tail = i
else:
if level and (not element.tail or not element.tail.strip()):
element.tail = i

167
openmc/arithmetic.py
View file

@ -17,7 +17,7 @@ _TALLY_ARITHMETIC_OPS = ['+', '-', '*', '/', '^']
_TALLY_AGGREGATE_OPS = ['sum', 'avg']
class CrossScore(object):
class CrossScore:
"""A special-purpose tally score used to encapsulate all combinations of two
tally's scores as an outer product for tally arithmetic.
@ -43,18 +43,10 @@ class CrossScore(object):
"""
def __init__(self, left_score=None, right_score=None, binary_op=None):
self._left_score = None
self._right_score = None
self._binary_op = None
if left_score is not None:
self.left_score = left_score
if right_score is not None:
self.right_score = right_score
if binary_op is not None:
self.binary_op = binary_op
def __init__(self, left_score, right_score, binary_op):
self.left_score = left_score
self.right_score = right_score
self.binary_op = binary_op
def __hash__(self):
return hash(repr(self))
@ -62,13 +54,9 @@ class CrossScore(object):
def __eq__(self, other):
return str(other) == str(self)
def __ne__(self, other):
return not self == other
def __repr__(self):
string = '({0} {1} {2})'.format(self.left_score,
self.binary_op, self.right_score)
return string
return '({} {} {})'.format(self.left_score, self.binary_op,
self.right_score)
@property
def left_score(self):
@ -101,7 +89,7 @@ class CrossScore(object):
self._binary_op = binary_op
class CrossNuclide(object):
class CrossNuclide:
"""A special-purpose nuclide used to encapsulate all combinations of two
tally's nuclides as an outer product for tally arithmetic.
@ -127,18 +115,10 @@ class CrossNuclide(object):
"""
def __init__(self, left_nuclide=None, right_nuclide=None, binary_op=None):
self._left_nuclide = None
self._right_nuclide = None
self._binary_op = None
if left_nuclide is not None:
self.left_nuclide = left_nuclide
if right_nuclide is not None:
self.right_nuclide = right_nuclide
if binary_op is not None:
self.binary_op = binary_op
def __init__(self, left_nuclide, right_nuclide, binary_op):
self.left_nuclide = left_nuclide
self.right_nuclide = right_nuclide
self.binary_op = binary_op
def __hash__(self):
return hash(repr(self))
@ -146,9 +126,6 @@ class CrossNuclide(object):
def __eq__(self, other):
return str(other) == str(self)
def __ne__(self, other):
return not self == other
def __repr__(self):
return self.name
@ -206,7 +183,7 @@ class CrossNuclide(object):
self._binary_op = binary_op
class CrossFilter(object):
class CrossFilter:
"""A special-purpose filter used to encapsulate all combinations of two
tally's filter bins as an outer product for tally arithmetic.
@ -239,26 +216,10 @@ class CrossFilter(object):
"""
def __init__(self, left_filter=None, right_filter=None, binary_op=None):
left_type = left_filter.type
right_type = right_filter.type
self._type = '({0} {1} {2})'.format(left_type, binary_op, right_type)
self._bins = {}
self._left_filter = None
self._right_filter = None
self._binary_op = None
if left_filter is not None:
self.left_filter = left_filter
self._bins['left'] = left_filter.bins
if right_filter is not None:
self.right_filter = right_filter
self._bins['right'] = right_filter.bins
if binary_op is not None:
self.binary_op = binary_op
def __init__(self, left_filter, right_filter, binary_op):
self.left_filter = left_filter
self.right_filter = right_filter
self.binary_op = binary_op
def __hash__(self):
return hash((self.left_filter, self.right_filter))
@ -266,21 +227,16 @@ class CrossFilter(object):
def __eq__(self, other):
return str(other) == str(self)
def __ne__(self, other):
return not self == other
def __repr__(self):
string = 'CrossFilter\n'
filter_type = '({0} {1} {2})'.format(self.left_filter.type,
self.binary_op,
self.right_filter.type)
filter_bins = '({0} {1} {2})'.format(self.left_filter.bins,
self.binary_op,
self.right_filter.bins)
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', filter_type)
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', filter_bins)
return string
filter_bins = '({} {} {})'.format(self.left_filter.bins,
self.binary_op,
self.right_filter.bins)
parts = [
'CrossFilter',
'{: <16}=\t{}'.format('\tType', self.type),
'{: <16}=\t{}'.format('\tBins', filter_bins)
]
return '\n'.join(parts)
@property
def left_filter(self):
@ -296,11 +252,13 @@ class CrossFilter(object):
@property
def type(self):
return self._type
left_type = self.left_filter.type
right_type = self.right_filter.type
return '({} {} {})'.format(left_type, self.binary_op, right_type)
@property
def bins(self):
return self._bins['left'], self._bins['right']
return self._left_filter.bins, self._right_filter.bins
@property
def num_bins(self):
@ -309,28 +267,17 @@ class CrossFilter(object):
else:
return 0
@type.setter
def type(self, filter_type):
if filter_type not in _FILTER_TYPES:
msg = 'Unable to set CrossFilter type to "{0}" since it ' \
'is not one of the supported types'.format(filter_type)
raise ValueError(msg)
self._type = filter_type
@left_filter.setter
def left_filter(self, left_filter):
cv.check_type('left_filter', left_filter,
(openmc.Filter, CrossFilter, AggregateFilter))
self._left_filter = left_filter
self._bins['left'] = left_filter.bins
@right_filter.setter
def right_filter(self, right_filter):
cv.check_type('right_filter', right_filter,
(openmc.Filter, CrossFilter, AggregateFilter))
self._right_filter = right_filter
self._bins['right'] = right_filter.bins
@binary_op.setter
def binary_op(self, binary_op):
@ -418,7 +365,8 @@ class CrossFilter(object):
return df
class AggregateScore(object):
class AggregateScore:
"""A special-purpose tally score used to encapsulate an aggregate of a
subset or all of tally's scores for tally aggregation.
@ -456,12 +404,9 @@ class AggregateScore(object):
def __eq__(self, other):
return str(other) == str(self)
def __ne__(self, other):
return not self == other
def __repr__(self):
string = ', '.join(map(str, self.scores))
string = '{0}({1})'.format(self.aggregate_op, string)
string = '{}({})'.format(self.aggregate_op, string)
return string
@property
@ -491,7 +436,7 @@ class AggregateScore(object):
self._aggregate_op = aggregate_op
class AggregateNuclide(object):
class AggregateNuclide:
"""A special-purpose tally nuclide used to encapsulate an aggregate of a
subset or all of tally's nuclides for tally aggregation.
@ -529,13 +474,10 @@ class AggregateNuclide(object):
def __eq__(self, other):
return str(other) == str(self)
def __ne__(self, other):
return not self == other
def __repr__(self):
# Append each nuclide in the aggregate to the string
string = '{0}('.format(self.aggregate_op)
string = '{}('.format(self.aggregate_op)
names = [nuclide.name if isinstance(nuclide, openmc.Nuclide)
else str(nuclide) for nuclide in self.nuclides]
string += ', '.join(map(str, names)) + ')'
@ -570,7 +512,7 @@ class AggregateNuclide(object):
self._aggregate_op = aggregate_op
class AggregateFilter(object):
class AggregateFilter:
"""A special-purpose tally filter used to encapsulate an aggregate of a
subset or all of a tally filter's bins for tally aggregation.
@ -600,17 +542,16 @@ class AggregateFilter(object):
"""
def __init__(self, aggregate_filter=None, bins=None, aggregate_op=None):
def __init__(self, aggregate_filter, bins=None, aggregate_op=None):
self._type = '{0}({1})'.format(aggregate_op,
aggregate_filter.short_name.lower())
self._type = '{}({})'.format(aggregate_op,
aggregate_filter.short_name.lower())
self._bins = None
self._aggregate_filter = None
self._aggregate_op = None
if aggregate_filter is not None:
self.aggregate_filter = aggregate_filter
self.aggregate_filter = aggregate_filter
if bins is not None:
self.bins = bins
if aggregate_op is not None:
@ -622,9 +563,6 @@ class AggregateFilter(object):
def __eq__(self, other):
return str(other) == str(self)
def __ne__(self, other):
return not self == other
def __gt__(self, other):
if self.type != other.type:
if self.aggregate_filter.type in _FILTER_TYPES and \
@ -641,10 +579,12 @@ class AggregateFilter(object):
return not self > other
def __repr__(self):
string = 'AggregateFilter\n'
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type)
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins)
return string
parts = [
'AggregateFilter',
'{: <16}=\t{}'.format('\tType', self.type),
'{: <16}=\t{}'.format('\tBins', self.bins)
]
return '\n'.join(parts)
@property
def aggregate_filter(self):
@ -669,7 +609,7 @@ class AggregateFilter(object):
@type.setter
def type(self, filter_type):
if filter_type not in _FILTER_TYPES:
msg = 'Unable to set AggregateFilter type to "{0}" since it ' \
msg = 'Unable to set AggregateFilter type to "{}" since it ' \
'is not one of the supported types'.format(filter_type)
raise ValueError(msg)
@ -722,7 +662,7 @@ class AggregateFilter(object):
if filter_bin not in self.bins:
msg = 'Unable to get the bin index for AggregateFilter since ' \
'"{0}" is not one of the bins'.format(filter_bin)
'"{}" is not one of the bins'.format(filter_bin)
raise ValueError(msg)
else:
return self.bins.index(filter_bin)
@ -800,12 +740,7 @@ class AggregateFilter(object):
return False
# None of the bins in this filter should match in the other filter
for bin in self.bins:
if bin in other.bins:
return False
# If all conditional checks passed then filters are mergeable
return True
return not any(b in other.bins for b in self.bins)
def merge(self, other):
"""Merge this aggregatefilter with another.
@ -823,8 +758,8 @@ class AggregateFilter(object):
"""
if not self.can_merge(other):
msg = 'Unable to merge "{0}" with "{1}" ' \
'filters'.format(self.type, other.type)
msg = 'Unable to merge "{}" with "{}" filters'.format(
self.type, other.type)
raise ValueError(msg)
# Create deep copy of filter to return as merged filter

74
openmc/cell.py
View file

@ -4,6 +4,7 @@ from copy import deepcopy
from math import cos, sin, pi
from numbers import Real
from xml.etree import ElementTree as ET
from uncertainties import UFloat
import sys
import warnings
@ -85,7 +86,12 @@ class Cell(IDManagerMixin):
volume : float
Volume of the cell in cm^3. This can either be set manually or
calculated in a stochastic volume calculation and added via the
:meth:`Cell.add_volume_information` method.
:meth:`Cell.add_volume_information` method. For 'distribmat' cells
it is the total volume of all instances.
atoms : collections.OrderedDict
Mapping of nuclides to the total number of atoms for each nuclide
present in the cell, or in all of its instances for a 'distribmat'
fill. For example, {'U235': 1.0e22, 'U238': 5.0e22, ...}.
"""
@ -134,6 +140,7 @@ class Cell(IDManagerMixin):
string += '\t{0: <15}=\t{1}\n'.format('Temperature',
self.temperature)
string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation)
string += '{: <16}=\t{}\n'.format('\tVolume', self.volume)
return string
@ -182,6 +189,55 @@ class Cell(IDManagerMixin):
def volume(self):
return self._volume
@property
def atoms(self):
if self._atoms is None:
if self._volume is None:
msg = ('Cannot calculate atom content becouse no volume '
'is set. Use Cell.volume to provide it or perform '
'a stochastic volume calculation.')
raise ValueError(msg)
elif self.fill_type == 'void':
msg = ('Cell is filled with void. It contains no atoms. '
'Material must be set to calculate atom content.')
raise ValueError(msg)
elif self.fill_type in ['lattice', 'universe']:
msg = ('Universe and Lattice cells can contain multiple '
'materials in diffrent proportions. Atom content must '
'be calculated with stochastic volume calculation.')
raise ValueError(msg)
elif self.fill_type == 'material':
# Get atomic densities
self._atoms = self._fill.get_nuclide_atom_densities()
# Convert to total number of atoms
for key, nuclide in self._atoms.items():
atom = nuclide[1] * self._volume * 1.0e+24
self._atoms[key] = atom
elif self.fill_type == 'distribmat':
# Assumes that volume is total volume of all instances
# Also assumes that all instances have the same volume
partial_volume = self.volume / len(self.fill)
self._atoms = OrderedDict()
for mat in self.fill:
for key, nuclide in mat.get_nuclide_atom_densities().items():
# To account for overlap of nuclides between distribmat
# we need to append new atoms to any existing value
# hence it is necessary to ask for default.
atom = self._atoms.setdefault(key, 0)
atom += nuclide[1] * partial_volume * 1.0e+24
self._atoms[key] = atom
else:
msg = 'Unrecognized fill_type: {}'.format(self.fill_type)
raise ValueError(msg)
return self._atoms
@property
def paths(self):
if self._paths is None:
@ -223,12 +279,16 @@ class Cell(IDManagerMixin):
elif not isinstance(fill, (openmc.Material, openmc.Lattice,
openmc.Universe)):
msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \
'Universe fill "{1}"'.format(self._id, fill)
msg = ('Unable to set Cell ID="{0}" to use a non-Material or '
'Universe fill "{1}"'.format(self._id, fill))
raise ValueError(msg)
self._fill = fill
# Info about atom content can now be invalid
# (since fill has just changed)
self._atoms = None
@rotation.setter
def rotation(self, rotation):
cv.check_length('cell rotation', rotation, 3)
@ -285,9 +345,15 @@ class Cell(IDManagerMixin):
@volume.setter
def volume(self, volume):
if volume is not None:
cv.check_type('cell volume', volume, Real)
cv.check_type('cell volume', volume, (Real, UFloat))
cv.check_greater_than('cell volume', volume, 0.0, equality=True)
self._volume = volume
# Info about atom content can now be invalid
# (sice volume has just changed)
self._atoms = None
def add_volume_information(self, volume_calc):
"""Add volume information to a cell.

64
openmc/checkvalue.py
View file

@ -4,7 +4,7 @@ from collections.abc import Iterable
import numpy as np
def check_type(name, value, expected_type, expected_iter_type=None):
def check_type(name, value, expected_type, expected_iter_type=None, *, none_ok=False):
"""Ensure that an object is of an expected type. Optionally, if the object is
iterable, check that each element is of a particular type.
@ -19,39 +19,44 @@ def check_type(name, value, expected_type, expected_iter_type=None):
expected_iter_type : type or Iterable of type or None, optional
Expected type of each element in value, assuming it is iterable. If
None, no check will be performed.
none_ok : bool, optional
Whether None is allowed as a value
"""
if none_ok and value is None:
return
if not isinstance(value, expected_type):
if isinstance(expected_type, Iterable):
msg = 'Unable to set "{0}" to "{1}" which is not one of the ' \
'following types: "{2}"'.format(name, value, ', '.join(
msg = 'Unable to set "{}" to "{}" which is not one of the ' \
'following types: "{}"'.format(name, value, ', '.join(
[t.__name__ for t in expected_type]))
else:
msg = 'Unable to set "{0}" to "{1}" which is not of type "{2}"'.format(
msg = 'Unable to set "{}" to "{}" which is not of type "{}"'.format(
name, value, expected_type.__name__)
raise TypeError(msg)
if expected_iter_type:
if isinstance(value, np.ndarray):
if not issubclass(value.dtype.type, expected_iter_type):
msg = 'Unable to set "{0}" to "{1}" since each item must be ' \
'of type "{2}"'.format(name, value,
expected_iter_type.__name__)
msg = 'Unable to set "{}" to "{}" since each item must be ' \
'of type "{}"'.format(name, value,
expected_iter_type.__name__)
raise TypeError(msg)
else:
return
for item in value:
if not isinstance(item, expected_iter_type):
if isinstance(expected_iter_type, Iterable):
msg = 'Unable to set "{0}" to "{1}" since each item must be ' \
'one of the following types: "{2}"'.format(
msg = 'Unable to set "{}" to "{}" since each item must be ' \
'one of the following types: "{}"'.format(
name, value, ', '.join([t.__name__ for t in
expected_iter_type]))
else:
msg = 'Unable to set "{0}" to "{1}" since each item must be ' \
'of type "{2}"'.format(name, value,
expected_iter_type.__name__)
msg = 'Unable to set "{}" to "{}" since each item must be ' \
'of type "{}"'.format(name, value,
expected_iter_type.__name__)
raise TypeError(msg)
@ -150,18 +155,18 @@ def check_length(name, value, length_min, length_max=None):
"""
if length_max is None:
if len(value) != length_min:
msg = 'Unable to set "{0}" to "{1}" since it must be of ' \
'length "{2}"'.format(name, value, length_min)
if len(value) < length_min:
msg = 'Unable to set "{}" to "{}" since it must be at least of ' \
'length "{}"'.format(name, value, length_min)
raise ValueError(msg)
elif not length_min <= len(value) <= length_max:
if length_min == length_max:
msg = 'Unable to set "{0}" to "{1}" since it must be of ' \
'length "{2}"'.format(name, value, length_min)
msg = 'Unable to set "{}" to "{}" since it must be of ' \
'length "{}"'.format(name, value, length_min)
else:
msg = 'Unable to set "{0}" to "{1}" since it must have length ' \
'between "{2}" and "{3}"'.format(name, value, length_min,
length_max)
msg = 'Unable to set "{}" to "{}" since it must have length ' \
'between "{}" and "{}"'.format(name, value, length_min,
length_max)
raise ValueError(msg)
@ -184,6 +189,7 @@ def check_value(name, value, accepted_values):
name, value, accepted_values)
raise ValueError(msg)
def check_less_than(name, value, maximum, equality=False):
"""Ensure that an object's value is less than a given value.
@ -211,6 +217,7 @@ def check_less_than(name, value, maximum, equality=False):
'or equal to "{2}"'.format(name, value, maximum)
raise ValueError(msg)
def check_greater_than(name, value, minimum, equality=False):
"""Ensure that an object's value is greater than a given value.
@ -265,13 +272,13 @@ def check_filetype_version(obj, expected_type, expected_version):
if this_version[0] != expected_version:
raise IOError('{} file has a version of {} which is not '
'consistent with the version expected by OpenMC, {}'
.format(this_filetype,
'.'.join(str(v) for v in this_version),
expected_version))
.format(this_filetype,
'.'.join(str(v) for v in this_version),
expected_version))
except AttributeError:
raise IOError('Could not read {} file. This most likely means the {} '
raise IOError('Could not read {} file. This most likely means the '
'file was produced by a different version of OpenMC than '
'the one you are using.'.format(expected_type))
'the one you are using.'.format(obj.filename))
class CheckedList(list):
@ -288,12 +295,13 @@ class CheckedList(list):
"""
def __init__(self, expected_type, name, items=[]):
def __init__(self, expected_type, name, items=None):
super().__init__()
self.expected_type = expected_type
self.name = name
for item in items:
self.append(item)
if items is not None:
for item in items:
self.append(item)
def __add__(self, other):
new_instance = copy.copy(self)

58
openmc/cmfd.py
View file

@ -55,7 +55,7 @@ _CURRENTS = {
}
class CMFDMesh(object):
class CMFDMesh:
"""A structured Cartesian mesh used for CMFD acceleration.
Attributes
@ -182,14 +182,14 @@ class CMFDMesh(object):
self._albedo = albedo
@map.setter
def map(self, meshmap):
check_type('CMFD mesh map', meshmap, Iterable, Integral)
for m in meshmap:
def map(self, mesh_map):
check_type('CMFD mesh map', mesh_map, Iterable, Integral)
for m in mesh_map:
check_value('CMFD mesh map', m, [0, 1])
self._map = meshmap
self._map = mesh_map
class CMFDRun(object):
class CMFDRun:
r"""Class for running CMFD acceleration through the C API.
Attributes
@ -535,7 +535,6 @@ class CMFDRun(object):
check_greater_than('CMFD feedback begin batch', begin, 0)
self._solver_begin = begin
@ref_d.setter
def ref_d(self, diff_params):
check_type('Reference diffusion params', diff_params,
@ -896,9 +895,8 @@ class CMFDRun(object):
cmfd_group.create_dataset('total_rate',
data=self._total_rate)
elif openmc.settings.verbosity >= 5:
print(' CMFD data not written to statepoint file'
'as it already exists in {}'.format(filename))
sys.stdout.flush()
print(' CMFD data not written to statepoint file as it '
'already exists in {}'.format(filename), flush=True)
def _initialize_linsolver(self):
# Determine number of rows in CMFD matrix
@ -1185,7 +1183,6 @@ class CMFDRun(object):
# Write CMFD output if CMFD on for current batch
self._write_cmfd_output()
def _cmfd_tally_reset(self):
"""Resets all CMFD tallies in memory"""
# Print message
@ -1753,7 +1750,7 @@ class CMFDRun(object):
# Compute new flux with C++ solver
innerits = openmc.lib._dll.openmc_run_linsolver(loss.data, s_o,
phi_n, toli)
phi_n, toli)
# Compute new source vector
s_n = prod.dot(phi_n)
@ -1785,7 +1782,7 @@ class CMFDRun(object):
# Update tolerance for inner iterations
toli = max(atoli, rtoli*norm_n)
def _check_convergence(self, s_n, s_o, k_n, k_o, iter, innerits):
def _check_convergence(self, s_n, s_o, k_n, k_o, iteration, innerits):
"""Checks the convergence of the CMFD problem
Parameters
@ -1798,9 +1795,9 @@ class CMFDRun(object):
K-effective from current iteration
k_o : float
K-effective from previous iteration
iter: int
iteration : int
Iteration number
innerits: int
innerits : int
Number of iterations required for convergence in inner GS loop
Returns
@ -1825,14 +1822,13 @@ class CMFDRun(object):
# Print out to user
if self._power_monitor and openmc.lib.master():
str1 = ' {:d}:'.format(iter)
str2 = 'k-eff: {:0.8f}'.format(k_n)
str3 = 'k-error: {:.5e}'.format(kerr)
str4 = 'src-error: {:.5e}'.format(serr)
str5 = ' {:d}'.format(innerits)
print('{:8s}{:20s}{:25s}{:s}{:s}'.format(str1, str2, str3, str4,
str5))
sys.stdout.flush()
print('{:8s}{:20s}{:25s}{:s}{:s}'.format(
' {:d}:'.format(iteration),
'k-eff: {:0.8f}'.format(k_n),
'k-error: {:.5e}'.format(kerr),
'src-error: {:.5e}'.format(serr),
' {:d}'.format(innerits)
), flush=True)
return iconv, serr
@ -1850,7 +1846,7 @@ class CMFDRun(object):
# Define coremap as cumulative sum over accelerated regions,
# otherwise set value to _CMFD_NOACCEL
self._coremap = np.where(self._coremap == 0, _CMFD_NOACCEL,
np.cumsum(self._coremap)-1)
np.cumsum(self._coremap) - 1)
# Reshape coremap to three dimensional array
# Indices of coremap in user input switched in x and z axes
@ -1929,19 +1925,19 @@ class CMFDRun(object):
') in group ' + str(ng-idx[-1])
raise OpenMCError(err_message)
# Get total rr from CMFD tally 0
# Get total reaction rate (rr) from CMFD tally 0
totalrr = tallies[tally_id].results[:,1,1]
# Reshape totalrr array to target shape. Swap x and z axes so that
# shape is now [nx, ny, nz, ng, 1]
# Reshape total reaction rate array to target shape. Swap x and z axes
# so that shape is now [nx, ny, nz, ng, 1]
reshape_totalrr = np.swapaxes(totalrr.reshape(target_tally_shape),
0, 2)
# Total rr is flipped in energy axis as tally results are given in
# reverse order of energy group
# Total reaction rate is flipped in energy axis as tally results are
# given in reverse order of energy group
reshape_totalrr = np.flip(reshape_totalrr, axis=3)
# Bank total rr to total_rate
# Bank total reaction rate to total_rate
self._total_rate = np.append(self._total_rate, reshape_totalrr,
axis=4)
@ -2047,7 +2043,7 @@ class CMFDRun(object):
axis=5)
# Compute current as aggregate of banked current_rate over tally window
self._current = np.where(is_accel[..., np.newaxis, np.newaxis],
self._current = np.where(is_accel[..., np.newaxis, np.newaxis],
np.sum(self._current_rate, axis=5), 0.0)
# Get p1 scatter rr from CMFD tally 3

4
openmc/data/endf.py
View file

@ -370,7 +370,7 @@ def get_evaluations(filename):
return evaluations
class Evaluation(object):
class Evaluation:
"""ENDF material evaluation with multiple files/sections
Parameters
@ -528,7 +528,7 @@ class Evaluation(object):
self.target['isomeric_state'])
class Tabulated2D(object):
class Tabulated2D:
"""Metadata for a two-dimensional function.
This is a dummy class that is not really used other than to store the

54
openmc/data/njoy.py
View file

@ -402,7 +402,6 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
# acer
if acer:
nacer_in = nlast
fname = '{}_{:.1f}'
for i, temperature in enumerate(temperatures):
# Extend input with an ACER run for each temperature
nace = nacer_in + 1 + 2*i
@ -411,8 +410,8 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
commands += _TEMPLATE_ACER.format(**locals())
# Indicate tapes to save for each ACER run
tapeout[nace] = fname.format("ace", temperature)
tapeout[ndir] = fname.format("xsdir", temperature)
tapeout[nace] = output_dir / "ace_{:.1f}".format(temperature)
tapeout[ndir] = output_dir / "xsdir_{:.1f}".format(temperature)
commands += 'stop\n'
run(commands, tapein, tapeout, **kwargs)
@ -422,7 +421,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file:
for temperature in temperatures:
# Get contents of ACE file
text = open(fname.format("ace", temperature), 'r').read()
text = (output_dir / "ace_{:.1f}".format(temperature)).read_text()
# If the target is metastable, make sure that ZAID in the ACE
# file reflects this by adding 400
@ -435,18 +434,18 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
ace_file.write(text)
# Concatenate into destination xsdir file
text = open(fname.format("xsdir", temperature), 'r').read()
xsdir_file.write(text)
xsdir_in = output_dir / "xsdir_{:.1f}".format(temperature)
xsdir_file.write(xsdir_in.read_text())
# Remove ACE/xsdir files for each temperature
for temperature in temperatures:
os.remove(fname.format("ace", temperature))
os.remove(fname.format("xsdir", temperature))
(output_dir / "ace_{:.1f}".format(temperature)).unlink()
(output_dir / "xsdir_{:.1f}".format(temperature)).unlink()
def make_ace_thermal(filename, filename_thermal, temperatures=None,
ace='ace', xsdir='xsdir', error=0.001, iwt=2,
evaluation=None, evaluation_thermal=None, **kwargs):
ace='ace', xsdir=None, output_dir=None, error=0.001,
iwt=2, evaluation=None, evaluation_thermal=None, **kwargs):
"""Generate thermal scattering ACE file from ENDF files
Parameters
@ -461,7 +460,12 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
ace : str, optional
Path of ACE file to write
xsdir : str, optional
Path of xsdir file to write
Path of xsdir file to write. Defaults to ``"xsdir"`` in the same
directory as ``ace``
output_dir : str, optional
Directory to write ace and xsdir files. If not provided, then write
output files to current directory. If given, must be a path to a
directory.
error : float, optional
Fractional error tolerance for NJOY processing
iwt : int
@ -481,6 +485,13 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
If the NJOY process returns with a non-zero status
"""
if output_dir is None:
output_dir = Path()
else:
output_dir = Path(output_dir)
if not output_dir.is_dir():
raise IOError("{} is not a directory".format(output_dir))
ev = evaluation if evaluation is not None else endf.Evaluation(filename)
mat = ev.material
zsymam = ev.target['zsymam']
@ -564,7 +575,6 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
# acer
nthermal_acer_in = nlast
fname = '{}_{:.1f}'
for i, temperature in enumerate(temperatures):
# Extend input with an ACER run for each temperature
nace = nthermal_acer_in + 1 + 2*i
@ -573,21 +583,23 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
commands += _THERMAL_TEMPLATE_ACER.format(**locals())
# Indicate tapes to save for each ACER run
tapeout[nace] = fname.format(ace, temperature)
tapeout[ndir] = fname.format(xsdir, temperature)
tapeout[nace] = output_dir / "ace_{:.1f}".format(temperature)
tapeout[ndir] = output_dir / "xsdir_{:.1f}".format(temperature)
commands += 'stop\n'
run(commands, tapein, tapeout, **kwargs)
with open(ace, 'w') as ace_file, open(xsdir, 'w') as xsdir_file:
ace = output_dir / ace
xsdir = (ace.parent / "xsdir") if xsdir is None else Path(xsdir)
with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file:
# Concatenate ACE and xsdir files together
for temperature in temperatures:
text = open(fname.format(ace, temperature), 'r').read()
ace_file.write(text)
ace_in = output_dir / "ace_{:.1f}".format(temperature)
ace_file.write(ace_in.read_text())
text = open(fname.format(xsdir, temperature), 'r').read()
xsdir_file.write(text)
xsdir_in = output_dir / "xsdir_{:.1f}".format(temperature)
xsdir_file.write(xsdir_in.read_text())
# Remove ACE/xsdir files for each temperature
for temperature in temperatures:
os.remove(fname.format(ace, temperature))
os.remove(fname.format(xsdir, temperature))
(output_dir / "ace_{:.1f}".format(temperature)).unlink()
(output_dir / "xsdir_{:.1f}".format(temperature)).unlink()

8
openmc/data/resonance.py
View file

@ -17,7 +17,7 @@ except ImportError:
import openmc.checkvalue as cv
class Resonances(object):
class Resonances:
"""Resolved and unresolved resonance data
Parameters
@ -119,7 +119,7 @@ class Resonances(object):
return cls(ranges)
class ResonanceRange(object):
class ResonanceRange:
"""Resolved resonance range
Parameters
@ -867,7 +867,7 @@ class RMatrixLimited(ResonanceRange):
return rml
class ParticlePair(object):
class ParticlePair:
def __init__(self, first, second, q_value, penetrability,
shift, mt):
self.first = first
@ -878,7 +878,7 @@ class ParticlePair(object):
self.mt = mt
class SpinGroup(object):
class SpinGroup:
"""Resonance spin group
Attributes

4
openmc/data/thermal.py
View file

@ -766,8 +766,8 @@ class ThermalScattering(EqualityMixin):
"""
with tempfile.TemporaryDirectory() as tmpdir:
# Run NJOY to create an ACE library
kwargs.setdefault('ace', os.path.join(tmpdir, 'ace'))
kwargs.setdefault('xsdir', os.path.join(tmpdir, 'xsdir'))
kwargs.setdefault('output_dir', tmpdir)
kwargs.setdefault('ace', os.path.join(kwargs['output_dir'], 'ace'))
kwargs['evaluation'] = evaluation
kwargs['evaluation_thermal'] = evaluation_thermal
make_ace_thermal(filename, filename_thermal, temperatures, **kwargs)

2
openmc/deplete/abc.py
View file

@ -190,7 +190,7 @@ class TransportOperator(ABC):
Returns
-------
volume : list of float
volume : dict of str to float
Volumes corresponding to materials in burn_list
nuc_list : list of str
A list of all nuclide names. Used for sorting the simulation.

2
openmc/deplete/atom_number.py
View file

@ -7,7 +7,7 @@ from collections import OrderedDict
import numpy as np
class AtomNumber(object):
class AtomNumber:
"""Stores local material compositions (atoms of each nuclide).
Parameters

192
openmc/deplete/chain.py
View file

@ -10,7 +10,7 @@ import math
import re
from collections import OrderedDict, defaultdict
from collections.abc import Mapping, Iterable
from numbers import Real
from numbers import Real, Integral
from warnings import warn
from openmc.checkvalue import check_type, check_greater_than
@ -141,7 +141,7 @@ _SECONDARY_PARTICLES = {
}
class Chain(object):
class Chain:
"""Full representation of a depletion chain.
A depletion chain can be created by using the :meth:`from_endf` method which
@ -500,7 +500,7 @@ class Chain(object):
# Gain
for _, target, branching_ratio in nuc.decay_modes:
# Allow for total annihilation for debug purposes
if target != 'Nothing':
if target is not None:
branch_val = branching_ratio * decay_constant
if branch_val != 0.0:
@ -525,17 +525,16 @@ class Chain(object):
matrix[i, i] -= path_rate
# Gain term; allow for total annihilation for debug purposes
if target != 'Nothing':
if r_type != 'fission':
if path_rate != 0.0:
k = self.nuclide_dict[target]
matrix[k, i] += path_rate * br
else:
for product, y in fission_yields[nuc.name].items():
yield_val = y * path_rate
if yield_val != 0.0:
k = self.nuclide_dict[product]
matrix[k, i] += yield_val
if r_type != 'fission':
if target is not None and path_rate != 0.0:
k = self.nuclide_dict[target]
matrix[k, i] += path_rate * br
else:
for product, y in fission_yields[nuc.name].items():
yield_val = y * path_rate
if yield_val != 0.0:
k = self.nuclide_dict[product]
matrix[k, i] += yield_val
# Clear set of reactions
reactions.clear()
@ -821,3 +820,168 @@ class Chain(object):
return stat
valid = valid and stat
return valid
def reduce(self, initial_isotopes, level=None):
"""Reduce the size of the chain by following transmutation paths
As an example, consider a simple chain with the following
isotopes and transmutation paths::
U235 (n,gamma) U236
(n,fission) (Xe135, I135, Cs135)
I135 (beta decay) Xe135 (beta decay) Cs135
Xe135 (n,gamma) Xe136
Calling ``chain.reduce(["I135"])`` will produce a depletion
chain that contains only isotopes that would originate from
I135: I135, Xe135, Cs135, and Xe136. U235 and U236 will not
be included, but multiple isotopes can be used to start
the search.
The ``level`` value controls the depth of the search.
``chain.reduce(["U235"], level=1)`` would return a chain
with all isotopes except Xe136, since it is two transmutations
removed from U235 in this case.
While targets will not be included in the new chain, the
total destruction rate and decay rate of included isotopes
will be preserved.
Parameters
----------
initial_isotopes : iterable of str
Start the search based on the contents of these isotopes
level : int, optional
Depth of transmuation path to follow. Must be greater than
or equal to zero. A value of zero returns a chain with
``initial_isotopes``. The default value of None implies
that all isotopes that appear in the transmutation paths
of the initial isotopes and their progeny should be
explored
Returns
-------
Chain
Depletion chain containing isotopes that would appear
after following up to ``level`` reactions and decay paths
"""
check_type("initial_isotopes", initial_isotopes, Iterable, str)
if level is None:
level = math.inf
else:
check_type("level", level, Integral)
check_greater_than("level", level, 0, equality=True)
all_isotopes = self._follow(set(initial_isotopes), level)
# Avoid re-sorting for fission yields
name_sort = sorted(all_isotopes)
nuclides = []
nuclide_dict = {}
reactions = set()
for idx, iso in enumerate(sorted(all_isotopes, key=openmc.data.zam)):
previous = self[iso]
new_nuclide = Nuclide(previous.name)
new_nuclide.half_life = previous.half_life
new_nuclide.decay_energy = new_nuclide.decay_energy
new_decay = []
for mode in previous.decay_modes:
if mode.target in all_isotopes:
new_decay.append(mode)
else:
new_decay.append(DecayTuple(
mode.type, None, mode.branching_ratio))
new_nuclide.decay_modes = new_decay
new_reactions = []
for rxn in previous.reactions:
if rxn.target in all_isotopes:
new_reactions.append(rxn)
reactions.add(rxn.type)
elif rxn.type == "fission":
new_yields = new_nuclide.yield_data = (
previous.yield_data.restrict_products(name_sort))
if new_yields is not None:
new_reactions.append(rxn)
reactions.add("fission")
# Maintain total destruction rates but set no target
else:
new_reactions.append(ReactionTuple(
rxn.type, None, rxn.Q, rxn.branching_ratio))
reactions.add(rxn.type)
new_nuclide.reactions = new_reactions
nuclides.append(new_nuclide)
nuclide_dict[iso] = idx
new_chain = type(self)()
new_chain.nuclides = nuclides
new_chain.nuclide_dict = nuclide_dict
# Doesn't appear that the ordering matters for the reactions,
# just the contents
new_chain.reactions = sorted(reactions)
return new_chain
def _follow(self, isotopes, level):
"""Return all isotopes present up to depth level"""
found = isotopes.copy()
remaining = set(self.nuclide_dict)
if not found.issubset(remaining):
raise IndexError(
"The following isotopes were not found in the chain: "
"{}".format(", ".join(found - remaining)))
if level == 0:
return found
remaining -= found
depth = 0
next_iso = set()
while depth < level and remaining:
# Exhaust all isotopes at this level
while isotopes:
iso = isotopes.pop()
found.add(iso)
nuclide = self[iso]
# Follow all transmutation paths for this nuclide
for rxn in nuclide.reactions + nuclide.decay_modes:
if rxn.type == "fission" or rxn.target is None:
continue
# Skip if we've already come across this isotope
elif (rxn.target in next_iso
or rxn.target in found or rxn.target in isotopes):
continue
next_iso.add(rxn.target)
if nuclide.yield_data is not None:
for product in nuclide.yield_data.products:
if (product in next_iso
or product in found or product in isotopes):
continue
next_iso.add(product)
if not next_iso:
# No additional isotopes to process, nor to update the
# current set of discovered isotopes
return found
# Prepare for next dig
depth += 1
isotopes |= next_iso
remaining -= next_iso
next_iso.clear()
# Process isotope that would have started next depth
found.update(isotopes)
return found

3
openmc/deplete/dummy_comm.py
View file

@ -1,6 +1,7 @@
import sys
class DummyCommunicator(object):
class DummyCommunicator:
rank = 0
size = 1

62
openmc/deplete/nuclide.py
View file

@ -13,7 +13,7 @@ try:
except ImportError:
import xml.etree.ElementTree as ET
from numpy import empty
from numpy import empty, searchsorted
from openmc.checkvalue import check_type
@ -30,8 +30,10 @@ Parameters
----------
type : str
Type of the decay mode, e.g., 'beta-'
target : str
Nuclide resulting from decay
target : str or None
Nuclide resulting from decay. A value of ``None`` implies the
target does not exist in the currently configured depletion
chain
branching_ratio : float
Branching ratio of the decay mode
@ -53,8 +55,11 @@ Parameters
----------
type : str
Type of the reaction, e.g., 'fission'
target : str
nuclide resulting from reaction
target : str or None
Nuclide resulting from reaction. A value of ``None``
implies either no single target, e.g. from fission,
or that the target nuclide is not considered
in the current depletion chain
Q : float
Q value of the reaction in [eV]
branching_ratio : float
@ -70,7 +75,7 @@ except AttributeError:
pass
class Nuclide(object):
class Nuclide:
"""Decay modes, reactions, and fission yields for a single nuclide.
Parameters
@ -179,6 +184,8 @@ class Nuclide(object):
for decay_elem in element.iter('decay'):
d_type = decay_elem.get('type')
target = decay_elem.get('target')
if target is not None and target.lower() == "nothing":
target = None
branching_ratio = float(decay_elem.get('branching_ratio'))
nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio))
@ -192,6 +199,8 @@ class Nuclide(object):
# just set null values
if r_type != 'fission':
target = reaction_elem.get('target')
if target is not None and target.lower() == "nothing":
target = None
else:
target = None
if fission_q is not None:
@ -226,7 +235,7 @@ class Nuclide(object):
for mode, daughter, br in self.decay_modes:
mode_elem = ET.SubElement(elem, 'decay')
mode_elem.set('type', mode)
mode_elem.set('target', daughter)
mode_elem.set('target', daughter or "Nothing")
mode_elem.set('branching_ratio', str(br))
elem.set('reactions', str(len(self.reactions)))
@ -234,7 +243,7 @@ class Nuclide(object):
rx_elem = ET.SubElement(elem, 'reaction')
rx_elem.set('type', rx)
rx_elem.set('Q', str(Q))
if rx != 'fission':
if rx != 'fission' or daughter is not None:
rx_elem.set('target', daughter)
if br != 1.0:
rx_elem.set('branching_ratio', str(br))
@ -397,9 +406,7 @@ class FissionYieldDistribution(Mapping):
for g_index, energy in enumerate(energies):
prod_map = fission_yields[energy]
for prod_ix, product in enumerate(ordered_prod):
yield_val = prod_map.get(product)
yield_matrix[g_index, prod_ix] = (
0.0 if yield_val is None else yield_val)
yield_matrix[g_index, prod_ix] = prod_map.get(product, 0.0)
self.energies = tuple(energies)
self.products = tuple(ordered_prod)
self.yield_matrix = yield_matrix
@ -435,7 +442,7 @@ class FissionYieldDistribution(Mapping):
FissionYieldDistribution
"""
all_yields = {}
for elem_index, yield_elem in enumerate(element.iter("fission_yields")):
for yield_elem in element.iter("fission_yields"):
energy = float(yield_elem.get("energy"))
products = yield_elem.find("products").text.split()
yields = map(float, yield_elem.find("data").text.split())
@ -460,6 +467,37 @@ class FissionYieldDistribution(Mapping):
data_elem = ET.SubElement(yield_element, "data")
data_elem.text = " ".join(map(str, yield_obj.yields))
def restrict_products(self, possible_products):
"""Return a new distribution with select products
Parameters
----------
possible_products : iterable of str
Candidate pool of fission products. Existing products
not contained here will not exist in the new instance
Returns
-------
FissionYieldDistribution or None
A value of None indicates no values in
``possible_products`` exist in :attr:`products`
"""
overlap = set(self.products).intersection(possible_products)
if not overlap:
return None
products = sorted(overlap)
indices = searchsorted(self.products, products)
# coerce back to dictionary to pass back to __init__
new_yields = {}
for ene, yields in zip(self.energies, self.yield_matrix.copy()):
new_yields[ene] = dict(zip(products, yields[indices]))
return type(self)(new_yields)
class FissionYield(Mapping):
"""Mapping for fission yields of a parent at a specific energy

23
openmc/deplete/operator.py
View file

@ -10,13 +10,10 @@ densities is all done in-memory instead of through the filesystem.
import sys
import copy
from collections import OrderedDict
from itertools import chain
import os
import time
import xml.etree.ElementTree as ET
from warnings import warn
import h5py
import numpy as np
from uncertainties import ufloat
@ -113,6 +110,13 @@ class Operator(TransportOperator):
``fission_yield_mode``. Will be passed directly on to the
helper. Passing a value of None will use the defaults for
the associated helper.
reduce_chain : bool, optional
If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
depletion chain up to ``reduce_chain_level``. Default is False.
reduce_chain_level : int, optional
Depth of the search when reducing the depletion chain. Only used
if ``reduce_chain`` evaluates to true. The default value of
``None`` implies no limit on the depth.
Attributes
----------
@ -158,7 +162,8 @@ class Operator(TransportOperator):
def __init__(self, geometry, settings, chain_file=None, prev_results=None,
diff_burnable_mats=False, energy_mode="fission-q",
fission_q=None, dilute_initial=1.0e3,
fission_yield_mode="constant", fission_yield_opts=None):
fission_yield_mode="constant", fission_yield_opts=None,
reduce_chain=False, reduce_chain_level=None):
if fission_yield_mode not in self._fission_helpers:
raise KeyError(
"fission_yield_mode must be one of {}, not {}".format(
@ -179,6 +184,16 @@ class Operator(TransportOperator):
self.geometry = geometry
self.diff_burnable_mats = diff_burnable_mats
# Reduce the chain before we create more materials
if reduce_chain:
all_isotopes = set()
for material in geometry.get_all_materials().values():
if not material.depletable:
continue
for name, _dens_percent, _dens_type in material.nuclides:
all_isotopes.add(name)
self.chain = self.chain.reduce(all_isotopes, reduce_chain_level)
# Differentiate burnable materials with multiple instances
if self.diff_burnable_mats:
self._differentiate_burnable_mats()

4
openmc/deplete/results.py
View file

@ -19,7 +19,7 @@ _VERSION_RESULTS = (1, 0)
__all__ = ["Results"]
class Results(object):
class Results:
"""Output of a depletion run
Attributes
@ -36,7 +36,7 @@ class Results(object):
Number of nuclides.
rates : list of ReactionRates
The reaction rates for each substep.
volume : OrderedDict of int to float
volume : OrderedDict of str to float
Dictionary mapping mat id to volume.
mat_to_ind : OrderedDict of str to int
A dictionary mapping mat ID as string to index.

32
openmc/deplete/results_list.py
View file

@ -2,7 +2,7 @@ import h5py
import numpy as np
from .results import Results, _VERSION_RESULTS
from openmc.checkvalue import check_filetype_version
from openmc.checkvalue import check_filetype_version, check_value
__all__ = ["ResultsList"]
@ -40,7 +40,7 @@ class ResultsList(list):
new.append(Results.from_hdf5(fh, i))
return new
def get_atoms(self, mat, nuc):
def get_atoms(self, mat, nuc, nuc_units="atoms", time_units="s"):
"""Get number of nuclides over time from a single material
.. note::
@ -57,15 +57,24 @@ class ResultsList(list):
Material name to evaluate
nuc : str
Nuclide name to evaluate
nuc_units : {"atoms", "atom/b-cm", "atom/cm3"}, optional
Units for the returned concentration. Default is ``"atoms"``
time_units : {"s", "min", "h", "d"}, optional
Units for the returned time array. Default is ``"s"`` to
return the value in seconds.
Returns
-------
time : numpy.ndarray
Array of times in [s]
Array of times in units of ``time_units``
concentration : numpy.ndarray
Total number of atoms for specified nuclide
Concentration of specified nuclide in units of ``nuc_units``
"""
check_value("time_units", time_units, {"s", "d", "min", "h"})
check_value("nuc_units", nuc_units,
{"atoms", "atom/b-cm", "atom/cm3"})
time = np.empty_like(self, dtype=float)
concentration = np.empty_like(self, dtype=float)
@ -74,6 +83,21 @@ class ResultsList(list):
time[i] = result.time[0]
concentration[i] = result[0, mat, nuc]
# Unit conversions
if time_units == "d":
time /= (60 * 60 * 24)
elif time_units == "h":
time /= (60 * 60)
elif time_units == "min":
time /= 60
if nuc_units != "atoms":
# Divide by volume to get density
concentration /= self[0].volume[mat]
if nuc_units == "atom/b-cm":
# 1 barn = 1e-24 cm^2
concentration *= 1e-24
return time, concentration
def get_reaction_rate(self, mat, nuc, rx):

5
openmc/element.py
View file

@ -40,7 +40,7 @@ class Element(str):
cross_sections=None):
"""Expand natural element into its naturally-occurring isotopes.
An optional cross_sections argument or the OPENMC_CROSS_SECTIONS
An optional cross_sections argument or the :envvar:`OPENMC_CROSS_SECTIONS`
environment variable is used to specify a cross_sections.xml file.
If the cross_sections.xml file is found, the element is expanded only
into the isotopes/nuclides present in cross_sections.xml. If no
@ -132,7 +132,6 @@ class Element(str):
# If a cross_sections library is present, check natural nuclides
# against the nuclides in the library
if cross_sections is not None:
library_nuclides = set()
tree = ET.parse(cross_sections)
root = tree.getroot()
@ -174,14 +173,12 @@ class Element(str):
# our knowledge of the common cross section libraries
# (ENDF, JEFF, and JENDL)
else:
# Add the mutual isotopes
for nuclide in mutual_nuclides:
abundances[nuclide] = NATURAL_ABUNDANCE[nuclide]
# Adjust the abundances for the absent nuclides
for nuclide in absent_nuclides:
if nuclide in ['O17', 'O18'] and 'O16' in mutual_nuclides:
abundances['O16'] += NATURAL_ABUNDANCE[nuclide]
elif nuclide == 'Ta180' and 'Ta181' in mutual_nuclides:

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