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

Merge remote-tracking branch 'upstream/develop' into delayed-tallies

Browse source
This commit is contained in:
Sam Shaner 2015-10-05 09:50:53 -04:00
commit 4c151d8a5c
367 changed files with 7183 additions and 16437 deletions
Download patch file Download diff file Expand all files Collapse all files

3
.gitignore vendored
View file

@ -36,6 +36,7 @@ src/xml-fortran/xmlreader
# Test results error file
results_error.dat
inputs_error.dat
# Test build files
tests/build/
@ -62,4 +63,4 @@ data/nndc
.idea/*
# IPython notebook checkpoints
.ipynb_checkpoints
.ipynb_checkpoints

57
docs/source/pythonapi/examples/pandas-dataframes.ipynb
View file

File diff suppressed because one or more lines are too long

68
docs/source/pythonapi/examples/post-processing.ipynb
View file

File diff suppressed because one or more lines are too long

45
docs/source/pythonapi/examples/tally-arithmetic.ipynb
View file

@ -182,20 +182,19 @@
"# Create fuel Cell\n",
"fuel_cell = openmc.Cell(name='1.6% Fuel')\n",
"fuel_cell.fill = fuel\n",
"fuel_cell.add_surface(fuel_outer_radius, halfspace=-1)\n",
"fuel_cell.region = -fuel_outer_radius\n",
"pin_cell_universe.add_cell(fuel_cell)\n",
"\n",
"# Create a clad Cell\n",
"clad_cell = openmc.Cell(name='1.6% Clad')\n",
"clad_cell.fill = zircaloy\n",
"clad_cell.add_surface(fuel_outer_radius, halfspace=+1)\n",
"clad_cell.add_surface(clad_outer_radius, halfspace=-1)\n",
"clad_cell.region = +fuel_outer_radius & -clad_outer_radius\n",
"pin_cell_universe.add_cell(clad_cell)\n",
"\n",
"# Create a moderator Cell\n",
"moderator_cell = openmc.Cell(name='1.6% Moderator')\n",
"moderator_cell.fill = water\n",
"moderator_cell.add_surface(clad_outer_radius, halfspace=+1)\n",
"moderator_cell.region = +clad_outer_radius\n",
"pin_cell_universe.add_cell(moderator_cell)"
]
},
@ -219,12 +218,7 @@
"root_cell.fill = pin_cell_universe\n",
"\n",
"# Add boundary planes\n",
"root_cell.add_surface(min_x, halfspace=+1)\n",
"root_cell.add_surface(max_x, halfspace=-1)\n",
"root_cell.add_surface(min_y, halfspace=+1)\n",
"root_cell.add_surface(max_y, halfspace=-1)\n",
"root_cell.add_surface(min_z, halfspace=+1)\n",
"root_cell.add_surface(max_z, halfspace=-1)\n",
"root_cell.region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z\n",
"\n",
"# Create root Universe\n",
"root_universe = openmc.Universe(universe_id=0, name='root universe')\n",
@ -358,7 +352,7 @@
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB98JFQMZGiFPL70AAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTUtMDktMjFUMTA6MDg6\nNTcrMDc6MDALr51VAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE1LTA5LTIxVDEwOjA4OjU3KzA3OjAw\nevIl6QAAAABJRU5ErkJggg==\n",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB98KAwI1HbUKyRQAAAPZSURBVGje7Zs7buMwEIZ9iey5\n0gyNjQpXKTYudIScgkdQYTfut1idwkdQkQNsYQO2Qj0sPiVK+mlQDmwgwIcgg8Cc4fCTSK5W4OeF\nkM8rHv+2I/rgxPZEPZgR7XtQxKdXYuUXJSUnBQ/9WCgo4vOSJ+WFUvF7E08mlia+rn7VcKXP8sRs\nzFX8b2MdX2y6v1Tw6MZUw4H4ojfIjD8mvn/qRL5p4+vvlMqvp2EhR8WBzfiz20hXORmP9fi/bM9E\neUFvV5H/0yRkeSbiGRfFJErxD9ENdz7Mbhig/h89fvtFdMiI/ePUIXV4lXju8K3DKv9NThOZ3q2K\nmUy6grxFES8rjeyic+FFQav+ncg3fXjH+Ts+/iibztFqOiZuZP/Z3OafPX40NGgST2r+uvQkXXp6\ncKvmr+r0e1Eef5um3+JHP3IFF1D/seNZJgaDmvY0Gav1s+2f1fqpIcublfKGt6apotG/NVx3SInW\ntLX+7Vg/Pv1YqOsnun6JSVdOXT/X7vk75f938QP+8OmSBs0fXtymMhJbf8qlPynYmpKCh7OB1fzN\nalOj1sl0ZAruHLiA+RM73pDe/VjMVP89+aTXwjyc/x5n+u991895/utrJTy8/06TXh0r/5JOa2Jm\nYmqi4r/vUm/H4wLmT+z4anhr05X+q6KUXhtzr/9qSff5L5uMT//V/NdU4YuBTPa/8P67l/6r44ds\n+hYuoP5jx9ciy6XTWlibBrmx8V/TdMfjkP+6pOsu/lvM9N90sf7r+f6m/65n+S8p/itN15v0UkW3\n/+48+PRfJX6S9Joo4g+G/1qYG9KroqP/WypcuvyXPf13wH89/hHef7MB6R3Cqn55U4rv4kfH3zaS\ngQuYP7HjVf89tXrbO+hfLdr+Ozv/SP1dgtQ/Ov8C+i/3+q/Zf2D/HWi6bjT6rym9I/v/03/b+LHS\n4cTg/utTsV7/net/Afzz4f0XGX84/2j9xZ4/sePR/of2X7D/o+vPo/sv6h9B/Bfxr9j1Hz2eN/hO\n8/wfff4A848+f/1A/530/I0+/8PvH9D3H9HnT+R49P0b+v4PfP/4E/wXfP8Mvf9G37/D/ovuP8Se\nP7Hj0f0vdP8tqP9O339cyv7p3P1fdP8Z3v9G999j13/seMax8x/o+ZN7+O+E8zdP/8XOf8Hnz9Dz\nb7HnT+x49PxlCp7/BM+fOv13wvnXBfivt2lMvD8TyH/Hnb+Gz3+j589jz5/Y8ej9h4D+W7qQmf57\nefqv239n3T+C7z+h969i13/seMax+3/o/cMcu/8Y2H9n3p+J6r98pv8m4fwXuH+M3n+OO3++AX9c\nlR+4PhbRAAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDE1LTEwLTAzVDA5OjUxOjU5KzA3OjAwJPCZIQAA\nACV0RVh0ZGF0ZTptb2RpZnkAMjAxNS0xMC0wM1QwOTo1MTo1OSswNzowMFWtIZ0AAAAASUVORK5C\nYII=\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
@ -569,8 +563,9 @@
" Copyright: 2011-2015 Massachusetts Institute of Technology\n",
" License: http://mit-crpg.github.io/openmc/license.html\n",
" Version: 0.7.0\n",
" Git SHA1: b167d70c877c516deca785801b9fa6f53fb0985b\n",
" Date/Time: 2015-09-21 10:25:26\n",
" Git SHA1: 71dfde8d12942170a9a8d91796ab40a6e9becaaf\n",
" Date/Time: 2015-10-03 09:53:29\n",
" OpenMP Threads: 4\n",
"\n",
" ===========================================================================\n",
" ========================> INITIALIZATION <=========================\n",
@ -625,20 +620,20 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 9.1800E-01 seconds\n",
" Reading cross sections = 6.5800E-01 seconds\n",
" Total time in simulation = 1.7037E+01 seconds\n",
" Time in transport only = 1.7024E+01 seconds\n",
" Time in inactive batches = 2.8600E+00 seconds\n",
" Time in active batches = 1.4177E+01 seconds\n",
" Time synchronizing fission bank = 4.0000E-03 seconds\n",
" Sampling source sites = 4.0000E-03 seconds\n",
" Total time for initialization = 3.7200E-01 seconds\n",
" Reading cross sections = 1.2400E-01 seconds\n",
" Total time in simulation = 4.7110E+00 seconds\n",
" Time in transport only = 4.6510E+00 seconds\n",
" Time in inactive batches = 6.1900E-01 seconds\n",
" Time in active batches = 4.0920E+00 seconds\n",
" Time synchronizing fission bank = 2.0000E-03 seconds\n",
" Sampling source sites = 2.0000E-03 seconds\n",
" SEND/RECV source sites = 0.0000E+00 seconds\n",
" Time accumulating tallies = 0.0000E+00 seconds\n",
" Total time for finalization = 1.0000E-03 seconds\n",
" Total time elapsed = 1.7971E+01 seconds\n",
" Calculation Rate (inactive) = 4370.63 neutrons/second\n",
" Calculation Rate (active) = 2645.13 neutrons/second\n",
" Total time for finalization = 2.0000E-03 seconds\n",
" Total time elapsed = 5.0950E+00 seconds\n",
" Calculation Rate (inactive) = 20193.9 neutrons/second\n",
" Calculation Rate (active) = 9164.22 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",

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

@ -720,9 +720,8 @@ Geometry Specification -- geometry.xml
The geometry in OpenMC is described using `constructive solid geometry`_ (CSG),
also sometimes referred to as combinatorial geometry. CSG allows a user to
create complex objects using Boolean operators on a set of simpler surfaces. In
the geometry model, each unique closed volume in defined by its bounding
surfaces. In OpenMC, most `quadratic surfaces`_ can be modeled and used as
bounding surfaces.
the geometry model, each unique volume is defined by its bounding surfaces. In
OpenMC, most `quadratic surfaces`_ can be modeled and used as bounding surfaces.
Every geometry.xml must have an XML declaration at the beginning of the file and
a root element named geometry. Within the root element the user can define any
@ -745,7 +744,7 @@ number of cells, surfaces, and lattices. Let us look at the following example:
<id>1</id>
<universe>0</universe>
<material>1</material>
<surfaces>-1</surfaces>
<region>-1</region>
</cell>
</geometry>
@ -763,7 +762,7 @@ could be written as:
<!-- This is a comment -->
<surface id="1" type="sphere" coeffs="0.0 0.0 0.0 5.0" boundary="vacuum" />
<cell id="1" universe="0" material="1" surfaces="-1" />
<cell id="1" universe="0" material="1" region="-1" />
</geometry>
@ -787,7 +786,8 @@ Each ``<surface>`` element can have the following attributes or sub-elements:
:type:
The type of the surfaces. This can be "x-plane", "y-plane", "z-plane",
"plane", "x-cylinder", "y-cylinder", "z-cylinder", or "sphere".
"plane", "x-cylinder", "y-cylinder", "z-cylinder", "sphere", "x-cone",
"y-cone", or "z-cone".
*Default*: None
@ -891,15 +891,29 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
*Default*: None
:surfaces:
A list of the ``ids`` for surfaces that bound this cell, e.g. if the cell
is on the negative side of surface 3 and the positive side of surface 5, the
bounding surfaces would be given as "-3 5".
:region:
A Boolean expression of half-spaces that defines the spatial region which
the cell occupies. Each half-space is identified by the unique ID of the
surface prefixed by `-` or `+` to indicate that it is the negative or
positive half-space, respectively. The `+` sign for a positive half-space
can be omitted. Valid Boolean operators are parentheses, union `|`,
complement `~`, and intersection. Intersection is implicit and indicated by
the presence of whitespace. The order of operator precedence is parentheses,
complement, intersection, and then union.
.. note:: The surface attribute/element can be omitted to make a cell fill
its entire universe.
As an example, the following code gives a cell that is the union of the
negative half-space of surface 3 and the complement of the intersection of
the positive half-space of surface 5 and the negative half-space of surface
2:
*Default*: No surfaces
.. code-block:: xml
<cell id="1" material="1" region="-3 | ~(5 -2)" />
.. note:: The ``region`` attribute/element can be omitted to make a cell
fill its entire universe.
*Default*: A region filling all space.
:rotation:
If the cell is filled with a universe, this element specifies the angles in

4
docs/source/usersguide/output/summary.rst
View file

@ -117,9 +117,9 @@ The current revision of the summary file format is 1.
Unique ID of the lattice which fills the cell. Only present if fill_type is
set to 'lattice'.
**/geometry/cells/cell <uid>/surfaces** (*int[]*)
**/geometry/cells/cell <uid>/region** (*char[]*)
Surface specification for the cell.
Region specification for the cell.
**/geometry/surfaces/surface <uid>/index** (*int*)

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

@ -1,6 +1,5 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
@ -54,12 +53,11 @@ cell2 = openmc.Cell(cell_id=100, name='cell 2')
cell3 = openmc.Cell(cell_id=101, name='cell 3')
cell4 = openmc.Cell(cell_id=2, name='cell 4')
# Register Surfaces with Cells
cell1.add_surface(surface=surf2, halfspace=-1)
cell2.add_surface(surface=surf1, halfspace=-1)
cell3.add_surface(surface=surf1, halfspace=+1)
cell4.add_surface(surface=surf2, halfspace=+1)
cell4.add_surface(surface=surf3, halfspace=-1)
# 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

134
examples/python/boxes/build-xml.py Normal file
View file

@ -0,0 +1,134 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 15
inactive = 5
particles = 10000
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate some Nuclides
h1 = openmc.Nuclide('H-1')
o16 = openmc.Nuclide('O-16')
u235 = openmc.Nuclide('U-235')
u238 = openmc.Nuclide('U-238')
# 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_nuclide(h1, 2.)
moderator.add_nuclide(o16, 1.)
moderator.add_s_alpha_beta('HH2O', '71t')
# Instantiate a MaterialsFile, register all Materials, and export to XML
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '71c'
materials_file.add_materials([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 and register the root Universe
geometry = openmc.Geometry()
geometry.root_universe = root
# Instantiate a GeometryFile, register Geometry, and export to XML
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
geometry_file.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a SettingsFile, set all runtime parameters, and export to XML
settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('point', [0., 0., 0.])
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 = 'cell'
# Instantiate a PlotsFile, add Plot, and export to XML
plot_file = openmc.PlotsFile()
plot_file.add_plot(plot)
plot_file.export_to_xml()

15
examples/python/lattice/hexagonal/build-xml.py
View file

@ -67,15 +67,12 @@ cell4 = openmc.Cell(cell_id=500, name='cell 4')
cell5 = openmc.Cell(cell_id=600, name='cell 5')
cell6 = openmc.Cell(cell_id=601, name='cell 6')
# Register Surfaces with Cells
cell1.add_surface(left, halfspace=+1)
cell1.add_surface(right, halfspace=-1)
cell1.add_surface(bottom, halfspace=+1)
cell1.add_surface(top, halfspace=-1)
cell2.add_surface(fuel_surf, halfspace=-1)
cell3.add_surface(fuel_surf, halfspace=+1)
cell5.add_surface(fuel_surf, halfspace=-1)
cell6.add_surface(fuel_surf, halfspace=+1)
# Use surface half-spaces to define regions
cell1.region = +left & -right & +bottom & -top
cell2.region = -fuel_surf
cell3.region = +fuel_surf
cell5.region = -fuel_surf
cell6.region = +fuel_surf
# Register Materials with Cells
cell2.fill = fuel

24
examples/python/lattice/nested/build-xml.py
View file

@ -66,21 +66,15 @@ cell6 = openmc.Cell(cell_id=202, name='cell 6')
cell7 = openmc.Cell(cell_id=301, name='cell 7')
cell8 = openmc.Cell(cell_id=302, name='cell 8')
# Register Surfaces with Cells
cell1.add_surface(left, halfspace=+1)
cell1.add_surface(right, halfspace=-1)
cell1.add_surface(bottom, halfspace=+1)
cell1.add_surface(top, halfspace=-1)
cell2.add_surface(left, halfspace=+1)
cell2.add_surface(right, halfspace=-1)
cell2.add_surface(bottom, halfspace=+1)
cell2.add_surface(top, halfspace=-1)
cell3.add_surface(fuel1, halfspace=-1)
cell4.add_surface(fuel1, halfspace=+1)
cell5.add_surface(fuel2, halfspace=-1)
cell6.add_surface(fuel2, halfspace=+1)
cell7.add_surface(fuel3, halfspace=-1)
cell8.add_surface(fuel3, halfspace=+1)
# Use surface half-space to define regions
cell1.region = +left & -right & +bottom & -top
cell2.region = +left & -right & +bottom & -top
cell3.region = -fuel1
cell4.region = +fuel1
cell5.region = -fuel2
cell6.region = +fuel2
cell7.region = -fuel3
cell8.region = +fuel3
# Register Materials with Cells
cell3.fill = fuel

20
examples/python/lattice/simple/build-xml.py
View file

@ -1,6 +1,5 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
@ -65,17 +64,14 @@ cell5 = openmc.Cell(cell_id=202, name='cell 5')
cell6 = openmc.Cell(cell_id=301, name='cell 6')
cell7 = openmc.Cell(cell_id=302, name='cell 7')
# Register Surfaces with Cells
cell1.add_surface(left, halfspace=+1)
cell1.add_surface(right, halfspace=-1)
cell1.add_surface(bottom, halfspace=+1)
cell1.add_surface(top, halfspace=-1)
cell2.add_surface(fuel1, halfspace=-1)
cell3.add_surface(fuel1, halfspace=+1)
cell4.add_surface(fuel2, halfspace=-1)
cell5.add_surface(fuel2, halfspace=+1)
cell6.add_surface(fuel3, halfspace=-1)
cell7.add_surface(fuel3, halfspace=+1)
# Use surface half-spaces to define regions
cell1.region = +left & -right & +bottom & -top
cell2.region = -fuel1
cell3.region = +fuel1
cell4.region = -fuel2
cell5.region = +fuel2
cell6.region = -fuel3
cell7.region = +fuel3
# Register Materials with Cells
cell2.fill = fuel

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

@ -1,6 +1,5 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
@ -132,17 +131,11 @@ gap = openmc.Cell(cell_id=2, name='cell 2')
clad = openmc.Cell(cell_id=3, name='cell 3')
water = openmc.Cell(cell_id=4, name='cell 4')
# Register Surfaces with Cells
fuel.add_surface(fuel_or, halfspace=-1)
gap.add_surface(fuel_or, halfspace=+1)
gap.add_surface(clad_ir, halfspace=-1)
clad.add_surface(clad_ir, halfspace=+1)
clad.add_surface(clad_or, halfspace=-1)
water.add_surface(clad_or, halfspace=+1)
water.add_surface(left, halfspace=+1)
water.add_surface(right, halfspace=-1)
water.add_surface(bottom, halfspace=+1)
water.add_surface(top, halfspace=-1)
# 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

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

@ -1,6 +1,5 @@
import openmc
###############################################################################
# Simulation Input File Parameters
###############################################################################
@ -52,13 +51,8 @@ surf6.boundary_type = 'reflective'
# Instantiate Cell
cell = openmc.Cell(cell_id=1, name='cell 1')
# Register Surfaces with Cell
cell.add_surface(surface=surf1, halfspace=+1)
cell.add_surface(surface=surf2, halfspace=-1)
cell.add_surface(surface=surf3, halfspace=+1)
cell.add_surface(surface=surf4, halfspace=-1)
cell.add_surface(surface=surf5, halfspace=+1)
cell.add_surface(surface=surf6, halfspace=-1)
# Use surface half-spaces to define region
cell.region = +surf1 & -surf2 & +surf3 & -surf4 & +surf5 & -surf6
# Register Material with Cell
cell.fill = fuel

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

@ -2,14 +2,14 @@
<geometry>
<!-- Definition of Cells -->
<cell id="1" universe="0" fill="37" surfaces="-2" />
<cell id="100" universe="37" material="40" surfaces="-1" />
<cell id="101" universe="37" material="41" surfaces="1" />
<cell id="2" universe="0" material="41" surfaces = "2 -3" />
<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>

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

@ -0,0 +1,39 @@
<?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>

23
examples/xml/boxes/materials.xml Normal file
View file

@ -0,0 +1,23 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" ao="1.0" />
</material>
<material id="2">
<density value="4.5" units="g/cc" />
<nuclide name="U-238" ao="1.0" />
</material>
<material id="3">
<density value="1.0" units="g/cc" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="H-1" ao="2.0" />
<sab name="HH2O" xs="71t" />
</material>
</materials>

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

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

16
examples/xml/boxes/settings.xml Normal file
View file

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

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

@ -1,14 +1,14 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="6" surfaces="1 -2 3 -4" />
<cell id="2" universe="5" fill="4" surfaces="1 -2 3 -4" />
<cell id="101" universe="1" material="1" surfaces="-5" />
<cell id="102" universe="1" material="2" surfaces="5" />
<cell id="201" universe="2" material="1" surfaces="-6" />
<cell id="202" universe="2" material="2" surfaces="6" />
<cell id="301" universe="3" material="1" surfaces="-7" />
<cell id="302" universe="3" material="2" surfaces="7" />
<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">

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

@ -1,13 +1,13 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="5" surfaces="1 -2 3 -4" />
<cell id="101" universe="1" material="1" surfaces="-5" />
<cell id="102" universe="1" material="2" surfaces="5" />
<cell id="201" universe="2" material="1" surfaces="-6" />
<cell id="202" universe="2" material="2" surfaces="6" />
<cell id="301" universe="3" material="1" surfaces="-7" />
<cell id="302" universe="3" material="2" surfaces="7" />
<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>

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

@ -19,9 +19,9 @@
<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" surfaces=" -1" /> <!-- UO2 Fuel -->
<cell id="2" material="2" surfaces="1 -2" /> <!-- Helium gap -->
<cell id="3" material="3" surfaces="2 -3" /> <!-- Zircaloy cladding -->
<cell id="4" material="4" surfaces="3 4 -5 6 -7" /> <!-- Borated water -->
<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>

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

@ -5,15 +5,15 @@
<cell id="1">
<universe>0</universe>
<material>1</material>
<surfaces>1 -2 3 -4 5 -6</surfaces>
<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="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="5" type="z-plane" coeffs="-1" boundary="reflective" />
<surface id="6" type="z-plane" coeffs="1" boundary="reflective" />
</geometry>

2
openmc/clean_xml.py
View file

@ -29,7 +29,7 @@ def sort_xml_elements(tree):
comment_elements.append((element, next_element))
# Now iterate over all tags and order the elements within each tag
for tag in tags:
for tag in sorted(list(tags)):
# Retrieve all of the elements for this tag
try:

10
openmc/material.py
View file

@ -1,4 +1,4 @@
from collections import Iterable
from collections import Iterable, OrderedDict
from copy import deepcopy
from numbers import Real, Integral
import warnings
@ -64,15 +64,15 @@ class Material(object):
self._density = None
self._density_units = ''
# A dictionary of Nuclides
# An ordered dictionary of Nuclides (order affects OpenMC results)
# Keys - Nuclide names
# Values - tuple (nuclide, percent, percent type)
self._nuclides = {}
self._nuclides = OrderedDict()
# A dictionary of Elements
# An ordered dictionary of Elements (order affects OpenMC results)
# Keys - Element names
# Values - tuple (element, percent, percent type)
self._elements = {}
self._elements = OrderedDict()
# If specified, a list of tuples of (table name, xs identifier)
self._sab = []

24
openmc/opencg_compatible.py
View file

@ -9,6 +9,8 @@ except ImportError:
raise ImportError(msg)
import openmc
from openmc.region import Intersection
from openmc.surface import Halfspace
# A dictionary of all OpenMC Materials created
@ -480,12 +482,24 @@ def get_opencg_cell(openmc_cell):
if openmc_cell.translation is not None:
opencg_cell.translation = openmc_cell.translation
surfaces = openmc_cell.surfaces
for surface_id in surfaces:
surface = surfaces[surface_id][0]
halfspace = surfaces[surface_id][1]
# Add surfaces to OpenCG cell from OpenMC cell region. Right now this only
# works if the region is a single half-space or an intersection of
# half-spaces, i.e., no complex cells.
region = openmc_cell.region
if isinstance(region, Halfspace):
surface = region.surface
halfspace = -1 if region.side == '-' else 1
opencg_cell.add_surface(get_opencg_surface(surface), halfspace)
elif isinstance(region, Intersection):
for node in region.nodes:
if not isinstance(node, Halfspace):
raise NotImplementedError("Complex cells not yet supported "
"in OpenCG.")
surface = node.surface
halfspace = -1 if node.side == '-' else 1
opencg_cell.add_surface(get_opencg_surface(surface), halfspace)
else:
raise NotImplementedError("Complex cells not yet supported in OpenCG.")
# Add the OpenMC Cell to the global collection of all OpenMC Cells
OPENMC_CELLS[cell_id] = openmc_cell

308
openmc/region.py Normal file
View file

@ -0,0 +1,308 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable
from openmc.checkvalue import check_type
class Region(object):
"""Region of space that can be assigned to a cell.
Region is an abstract base class that is inherited by Halfspace,
Intersection, Union, and Complement. Each of those respective classes are
typically not instantiated directly but rather are created through operators
of the Surface and Region classes.
"""
__metaclass__ = ABCMeta
def __and__(self, other):
return Intersection(self, other)
def __or__(self, other):
return Union(self, other)
def __invert__(self):
return Complement(self)
@abstractmethod
def __str__(self):
return ''
@staticmethod
def from_expression(expression, surfaces):
"""Generate a region given an infix expression.
Parameters
----------
expression : str
Boolean expression relating surface half-spaces. The possible
operators are union '|', intersection ' ', and complement '~'. For
example, '(1 -2) | 3 ~(4 -5)'.
surfaces : dict
Dictionary whose keys are suface IDs that appear in the Boolean
expression and whose values are Surface objects.
"""
# Strip leading and trailing whitespace
expression = expression.strip()
# Convert the string expression into a list of tokens, i.e., operators
# and surface half-spaces, representing the expression in infix
# notation.
i = 0
i_start = -1
tokens = []
while i < len(expression):
if expression[i] in '()|~ ':
# If special character appears immediately after a non-operator,
# create a token with the apporpriate half-space
if i_start >= 0:
j = int(expression[i_start:i])
if j < 0:
tokens.append(-surfaces[abs(j)])
else:
tokens.append(+surfaces[abs(j)])
if expression[i] in '()|~':
# For everything other than intersection, add the operator
# to the list of tokens
tokens.append(expression[i])
else:
# Find next non-space character
while expression[i+1] == ' ':
i += 1
# If previous token is a halfspace or right parenthesis and next token
# is not a left parenthese or union operator, that implies that the
# whitespace is to be interpreted as an intersection operator
if (i_start >= 0 or tokens[-1] == ')') and \
expression[i+1] not in ')|':
tokens.append(' ')
i_start = -1
else:
# Check for invalid characters
if expression[i] not in '-0123456789':
raise SyntaxError("Invalid character '{}' in expression"
.format(expression[i]))
# If we haven't yet reached the start of a word, start one
if i_start < 0:
i_start = i
i += 1
# If we've reached the end and we're still in a word, create a
# half-space token and add it to the list
if i_start >= 0:
j = int(expression[i_start:])
if j < 0:
tokens.append(-surfaces[abs(j)])
else:
tokens.append(+surfaces[abs(j)])
# The functions below are used to apply an operator to operands on the
# output queue during the shunting yard algorithm.
def can_be_combined(region):
return isinstance(region, Complement) or hasattr(region, 'surface')
def apply_operator(output, operator):
r2 = output.pop()
if operator == ' ':
r1 = output.pop()
if isinstance(r1, Intersection) and can_be_combined(r2):
r1.nodes.append(r2)
output.append(r1)
elif isinstance(r2, Intersection) and can_be_combined(r1):
r2.nodes.insert(0, r1)
output.append(r2)
elif isinstance(r1, Intersection) and isinstance(r2, Intersection):
r1.nodes += r2.nodes
output.append(r1)
else:
output.append(Intersection(r1, r2))
elif operator == '|':
r1 = output.pop()
if isinstance(r1, Union) and can_be_combined(r2):
r1.nodes.append(r2)
output.append(r1)
elif isinstance(r2, Union) and can_be_combined(r1):
r2.nodes.insert(0, r1)
output.append(r2)
elif isinstance(r1, Union) and isinstance(r2, Union):
r1.nodes += r2.nodes
output.append(r1)
else:
output.append(Union(r1, r2))
elif operator == '~':
output.append(Complement(r2))
# The following is an implementation of the shunting yard algorithm to
# generate an abstract syntax tree for the region expression.
output = []
stack = []
precedence = {'|': 1, ' ': 2, '~': 3}
associativity = {'|': 'left', ' ': 'left', '~': 'right'}
for token in tokens:
if token in (' ', '|', '~'):
# Normal operators
while stack:
op = stack[-1]
if (op not in ('(', ')') and
((associativity[token] == 'right' and
precedence[token] < precedence[op]) or
(associativity[token] == 'left' and
precedence[token] <= precedence[op]))):
apply_operator(output, stack.pop())
else:
break
stack.append(token)
elif token == '(':
# Left parentheses
stack.append(token)
elif token == ')':
# Right parentheses
while stack[-1] != '(':
apply_operator(output, stack.pop())
if len(stack) == 0:
raise SyntaxError('Mismatched parentheses in '
'region specification.')
stack.pop()
else:
# Surface halfspaces
output.append(token)
while stack:
if stack[-1] in '()':
raise SyntaxError('Mismatched parentheses in region '
'specification.')
apply_operator(output, stack.pop())
# Since we are generating an abstract syntax tree rather than a reverse
# Polish notation expression, the output queue should have a single item
# at the end
return output[0]
class Intersection(Region):
"""Intersection of two or more regions.
Instances of Intersection are generally created via the __and__ operator
applied to two instances of Region. This is illustrated in the following
example:
>>> equator = openmc.surface.ZPlane(z0=0.0)
>>> earth = openmc.surface.Sphere(R=637.1e6)
>>> northern_hemisphere = -earth & +equator
>>> southern_hemisphere = -earth & -equator
>>> type(northern_hemisphere)
<class 'openmc.region.Intersection'>
Parameters
----------
*nodes
Regions to take the intersection of
Attributes
----------
nodes : tuple of Region
Regions to take the intersection of
"""
def __init__(self, *nodes):
self.nodes = list(nodes)
@property
def nodes(self):
return self._nodes
@nodes.setter
def nodes(self, nodes):
check_type('nodes', nodes, Iterable, Region)
self._nodes = nodes
def __str__(self):
return '(' + ' '.join(map(str, self.nodes)) + ')'
class Union(Region):
"""Union of two or more regions.
Instances of Union are generally created via the __or__ operator applied to
two instances of Region. This is illustrated in the following example:
>>> s1 = openmc.surface.ZPlane(z0=0.0)
>>> s2 = openmc.surface.Sphere(R=637.1e6)
>>> type(-s2 | +s1)
<class 'openmc.region.Union'>
Parameters
----------
*nodes
Regions to take the union of
Attributes
----------
nodes : tuple of Region
Regions to take the union of
"""
def __init__(self, *nodes):
self.nodes = list(nodes)
@property
def nodes(self):
return self._nodes
@nodes.setter
def nodes(self, nodes):
check_type('nodes', nodes, Iterable, Region)
self._nodes = nodes
def __str__(self):
return '(' + ' | '.join(map(str, self.nodes)) + ')'
class Complement(Region):
"""Complement of a region.
The Complement of an existing Region can be created by using the __invert__
operator as the following example demonstrates:
>>> xl = openmc.surface.XPlane(x0=-10.0)
>>> xr = openmc.surface.XPlane(x0=10.0)
>>> yl = openmc.surface.YPlane(y0=-10.0)
>>> yr = openmc.surface.YPlane(y0=10.0)
>>> inside_box = +xl & -xr & +yl & -yl
>>> outside_box = ~inside_box
>>> type(outside_box)
<class 'openmc.region.Complement'>
Parameters
----------
node : Region
Region to take the complement of
Attributes
----------
node : Region
Regions to take the complement of
"""
def __init__(self, node):
self.node = node
@property
def node(self):
return self._node
@node.setter
def node(self, node):
check_type('node', node, Region)
self._node = node
def __str__(self):
return '~' + str(self.node)

18
openmc/summary.py
View file

@ -1,6 +1,7 @@
import numpy as np
import openmc
from openmc.region import Region
class Summary(object):
@ -213,10 +214,10 @@ class Summary(object):
else:
fill = self._f['geometry/cells'][key]['lattice'].value
if 'surfaces' in self._f['geometry/cells'][key].keys():
surfaces = self._f['geometry/cells'][key]['surfaces'][...]
if 'region' in self._f['geometry/cells'][key].keys():
region = self._f['geometry/cells'][key]['region'].value.decode()
else:
surfaces = []
region = []
# Create this Cell
cell = openmc.Cell(cell_id=cell_id, name=name)
@ -241,13 +242,10 @@ class Summary(object):
# Store Cell fill information for after Universe/Lattice creation
self._cell_fills[index] = (fill_type, fill)
# Iterate over all Surfaces and add them to the Cell
for surface_halfspace in surfaces:
halfspace = np.sign(surface_halfspace)
surface_id = abs(surface_halfspace)
surface = self.get_surface_by_id(surface_id)
cell.add_surface(surface, halfspace)
# Generate Region object given infix expression
if region:
cell.region = Region.from_expression(
region, {s.id: s for s in self.surfaces.values()})
# Add the Cell to the global dictionary of all Cells
self.cells[index] = cell

82
openmc/surface.py
View file

@ -4,6 +4,7 @@ from xml.etree import ElementTree as ET
import sys
from openmc.checkvalue import check_type, check_value, check_greater_than
from openmc.region import Region
if sys.version_info[0] >= 3:
basestring = str
@ -38,17 +39,17 @@ class Surface(object):
Attributes
----------
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
Boundary condition that defines the behavior for particles hitting the
surface.
coeffs : dict
Dictionary of surface coefficients
id : int
Unique identifier for the surface
name : str
Name of the surface
type : str
Type of the surface, e.g. 'x-plane'
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
Boundary condition that defines the behavior for particles hitting the
surface.
coeffs : dict
Dictionary of surface coefficients
"""
@ -68,6 +69,12 @@ class Surface(object):
# proper order
self._coeff_keys = []
def __neg__(self):
return Halfspace(self, '-')
def __pos__(self):
return Halfspace(self, '+')
@property
def id(self):
return self._id
@ -940,3 +947,68 @@ class ZCone(Cone):
R2, name=name)
self._type = 'z-cone'
class Halfspace(Region):
"""A positive or negative half-space region.
A half-space is either of the two parts into which a two-dimension surface
divides the three-dimensional Euclidean space. If the equation of the
surface is :math:`f(x,y,z) = 0`, the region for which :math:`f(x,y,z) < 0`
is referred to as the negative half-space and the region for which
:math:`f(x,y,z) > 0` is referred to as the positive half-space.
Instances of Halfspace are generally not instantiated directly. Rather, they
can be created from an existing Surface through the __neg__ and __pos__
operators, as the following example demonstrates:
>>> sphere = openmc.surface.Sphere(surface_id=1, R=10.0)
>>> inside_sphere = -sphere
>>> outside_sphere = +sphere
>>> type(inside_sphere)
<class 'openmc.surface.Halfspace'>
Parameters
----------
surface : Surface
Surface which divides Euclidean space.
side : {'+', '-'}
Indicates whether the positive or negative half-space is used.
Attributes
----------
surface : Surface
Surface which divides Euclidean space.
side : {'+', '-'}
Indicates whether the positive or negative half-space is used.
"""
def __init__(self, surface, side):
self.surface = surface
self.side = side
def __invert__(self):
return -self.surface if self.side == '+' else +self.surface
@property
def surface(self):
return self._surface
@surface.setter
def surface(self, surface):
check_type('surface', surface, Surface)
self._surface = surface
@property
def side(self):
return self._side
@side.setter
def side(self, side):
check_value('side', side, ('+', '-'))
self._side = side
def __str__(self):
return '-' + str(self.surface.id) if self.side == '-' \
else str(self.surface.id)

12
openmc/temp.py
View file

@ -1,12 +0,0 @@
from checkvalue import *
from checkvalue import _isinstance
import numpy as np
zs = np.zeros((2,))
print _isinstance(zs[0], Integral)
print _isinstance(zs[0], Real)
print _isinstance(zs[0], (Integral, Real))
print check_iterable_type('thing', zs, (Real, Integral))

107
openmc/universe.py
View file

@ -3,11 +3,14 @@ from collections import OrderedDict, Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
import warnings
import numpy as np
import openmc
import openmc.checkvalue as cv
from openmc.surface import Halfspace
from openmc.region import Region, Intersection, Complement
if sys.version_info[0] >= 3:
basestring = str
@ -45,10 +48,8 @@ class Cell(object):
Name of the cell
fill : Material or Universe or Lattice or 'void'
Indicates what the region of space is filled with
surfaces : dict
Dictionary whose keys are surface IDs and values are 2-tuples of a
Surface object and an integer identify whether the positive or negative
half-space is to be used
region : openmc.region.Region
Region of space that is assigned to the cell.
rotation : ndarray
If the cell is filled with a universe, this array specifies the angles
in degrees about the x, y, and z axes that the filled universe should be
@ -67,7 +68,7 @@ class Cell(object):
self.name = name
self._fill = None
self._type = None
self._surfaces = {}
self._region = None
self._rotation = None
self._translation = None
self._offsets = None
@ -96,8 +97,8 @@ class Cell(object):
return None
@property
def surfaces(self):
return self._surfaces
def region(self):
return self._region
@property
def rotation(self):
@ -173,6 +174,11 @@ class Cell(object):
cv.check_type('cell offsets', offsets, Iterable)
self._offsets = offsets
@region.setter
def region(self, region):
cv.check_type('cell region', region, Region)
self._region = region
def add_surface(self, surface, halfspace):
"""Add a half-space to the list of half-spaces whose intersection defines the
cell.
@ -186,6 +192,12 @@ class Cell(object):
"""
warnings.simplefilter('always', DeprecationWarning)
warnings.warn("Cell.add_surface(...) has been deprecated and may be "
"removed in a future version. The region for a Cell "
"should be defined using the region property directly.",
DeprecationWarning)
if not isinstance(surface, openmc.Surface):
msg = 'Unable to add Surface "{0}" to Cell ID="{1}" since it is ' \
'not a Surface object'.format(surface, self._id)
@ -196,28 +208,17 @@ class Cell(object):
'"{2}" since it is not +/-1'.format(surface, self._id, halfspace)
raise ValueError(msg)
# If the Cell does not already contain the Surface, add it
if surface._id not in self._surfaces:
self._surfaces[surface._id] = (surface, halfspace)
def remove_surface(self, surface):
"""Remove the half-space associated with a particular surface.
Parameters
----------
surface : openmc.surface.Surface
Surface to remove from definition
"""
if not isinstance(surface, openmc.Surface):
msg = 'Unable to remove Surface "{0}" from Cell ID="{1}" since it is ' \
'not a Surface object'.format(surface, self._id)
raise ValueError(msg)
# If the Cell contains the Surface, delete it
if surface._id in self._surfaces:
del self._surfaces[surface._id]
# If no region has been assigned, simply use the half-space. Otherwise,
# take the intersection of the current region and the half-space
# specified
region = +surface if halfspace == 1 else -surface
if self.region is None:
self.region = region
else:
if isinstance(self.region, Intersection):
self.region.nodes.append(region)
else:
self.region = Intersection(self.region, region)
def get_offset(self, path, filter_offset):
# Get the current element and remove it from the list
@ -311,13 +312,7 @@ class Cell(object):
else:
string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', self._fill)
string += '{0: <16}{1}\n'.format('\tSurfaces', '=\t')
for surface_id in self._surfaces:
halfspace = self._surfaces[surface_id][1]
string += '{0} '.format(halfspace * surface_id)
string = string.rstrip(' ') + '\n'
string += '{0: <16}{1}{2}\n'.format('\tRegion', '=\t', self._region)
string += '{0: <16}{1}{2}\n'.format('\tRotation', '=\t',
self._rotation)
@ -348,26 +343,30 @@ class Cell(object):
element.set("fill", str(self._fill))
self._fill.create_xml_subelement(xml_element)
if self._surfaces is not None:
surfaces = ''
if self.region is not None:
# Set the region attribute with the region specification
element.set("region", str(self.region))
for surface_id in self._surfaces:
# Determine if XML element already includes this Surface
path = './surface[@id=\'{0}\']'.format(surface_id)
test = xml_element.find(path)
# Only surfaces that appear in a region are added to the geometry
# file, so the appropriate check is performed here. First we create
# a function which is called recursively to navigate through the CSG
# tree. When it reaches a leaf (a Halfspace), it creates a <surface>
# element for the corresponding surface if none has been created
# thus far.
def create_surface_elements(node, element):
if isinstance(node, Halfspace):
path = './surface[@id=\'{0}\']'.format(node.surface.id)
if xml_element.find(path) is None:
surface_subelement = node.surface.create_xml_subelement()
xml_element.append(surface_subelement)
elif isinstance(node, Complement):
create_surface_elements(node.node, element)
else:
for subnode in node.nodes:
create_surface_elements(subnode, element)
# If the element does not contain the Surface subelement
if test is None:
# Create the XML subelement for this Surface
surface = self._surfaces[surface_id][0]
surface_subelement = surface.create_xml_subelement()
xml_element.append(surface_subelement)
# Append the halfspace and Surface ID
halfspace = self._surfaces[surface_id][1]
surfaces += '{0} '.format(halfspace * surface_id)
element.set("surfaces", surfaces.rstrip(' '))
# Call the recursive function from the top node
create_surface_elements(self.region, xml_element)
if self._translation is not None:
element.set("translation", ' '.join(map(str, self._translation)))

32
scripts/openmc-update-inputs
View file

@ -1,14 +1,15 @@
#!/usr/bin/env python
"""Update OpenMC's input XML files to the latest format.
Usage information can be obtained by running 'update_inputs.py --help':
Usage information can be obtained by running 'openmc-update-inputs --help':
usage: update_lattices.py [-h] IN [IN ...]
usage: openmc-update-inputs [-h] IN [IN ...]
Update lattices in geometry.xml files to the latest format. This will remove
'outside' attributes/elements and replace them with 'outer' attributes. Note
that this script will not delete the given files; it will append '.original'
to the given files and write new ones.
Update geometry.xml files to the latest format. This will remove 'outside'
attributes/elements from lattices and replace them with 'outer' attributes. For
'cell' elements, any 'surfaces' attributes/elements will be renamed
'region'. Note that this script will not delete the given files; it will append
'.original' to the given files and write new ones.
positional arguments:
IN Input geometry.xml file(s).
@ -35,9 +36,10 @@ they will be moved to a new file with '.original' appended to their name.
Formatting changes that will be made:
geometry.xml: Lattices containing 'outside' attributes/tags will be replaced
geometry.xml: Lattices containing 'outside' attributes/tags will be replaced
with lattices containing 'outer' attributes, and the appropriate
cells/universes will be added.
cells/universes will be added. Any 'surfaces' attributes/elements on a cell
will be renamed 'region'.
"""
@ -173,9 +175,6 @@ def update_geometry(geometry_root):
root = geometry_root
was_updated = False
# Ignore files that do not contain lattices.
if all([child.tag != 'lattice' for child in root]): return False
# Get a set of already-used universe and cell ids.
uids = get_universe_ids(root)
cids = get_cell_ids(root)
@ -233,6 +232,17 @@ def update_geometry(geometry_root):
del lat.attrib['width']
was_updated = True
# Change 'surfaces' to 'region' in cell definitions
for cell in root.iter('cell'):
elem = cell.find('surfaces')
if elem is not None:
elem.tag = 'region'
was_updated = True
if 'surfaces' in cell.attrib:
cell.attrib['region'] = cell.attrib['surfaces']
del cell.attrib['surfaces']
was_updated = True
return was_updated

9
src/constants.F90
View file

@ -87,10 +87,11 @@ module constants
! Logical operators for cell definitions
integer, parameter :: &
OP_LEFT_PAREN = huge(0), & ! Left parentheses
OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses
OP_UNION = huge(0) - 2, & ! Union operator
OP_DIFFERENCE = huge(0) - 3 ! Difference operator
OP_LEFT_PAREN = huge(0), & ! Left parentheses
OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses
OP_COMPLEMENT = huge(0) - 2, & ! Complement operator (~)
OP_INTERSECTION = huge(0) - 3, & ! Intersection operator
OP_UNION = huge(0) - 4 ! Union operator (^)
! Cell types
integer, parameter :: &

902
src/geometry.F90
View file

File diff suppressed because it is too large Load diff

27
src/geometry_header.F90
View file

@ -113,22 +113,6 @@ module geometry_header
class(Lattice), allocatable :: obj
end type LatticeContainer
!===============================================================================
! SURFACE type defines a first- or second-order surface that can be used to
! construct closed volumes (cells)
!===============================================================================
type Surface
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of surface
real(8), allocatable :: coeffs(:) ! Definition of surface
integer, allocatable :: &
neighbor_pos(:), & ! List of cells on positive side
neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
end type Surface
!===============================================================================
! CELL defines a closed volume by its bounding surfaces
!===============================================================================
@ -144,14 +128,13 @@ module geometry_header
! the geom
integer :: material ! Material within cell (0 for
! universe)
integer :: n_surfaces ! Number of surfaces within
integer, allocatable :: offset (:) ! Distribcell offset for tally
! counter
integer, allocatable :: &
& surfaces(:) ! List of surfaces bounding cell
! -- note that parentheses, union,
! etc operators will be listed here
! too
integer, allocatable :: region(:) ! Definition of spatial region as
! Boolean expression of half-spaces
integer, allocatable :: rpn(:) ! Reverse Polish notation for region
! expression
logical :: simple ! Is the region simple (intersections only)
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)

15
src/global.F90
View file

@ -6,11 +6,12 @@ module global
use cmfd_header
use constants
use dict_header, only: DictCharInt, DictIntInt
use geometry_header, only: Cell, Universe, Lattice, LatticeContainer, Surface
use geometry_header, only: Cell, Universe, Lattice, LatticeContainer
use material_header, only: Material
use mesh_header, only: RegularMesh
use plot_header, only: ObjectPlot
use set_header, only: SetInt
use surface_header, only: SurfaceContainer
use source_header, only: ExtSource
use tally_header, only: TallyObject, TallyMap, TallyResult
use trigger_header, only: KTrigger
@ -27,12 +28,12 @@ module global
! GEOMETRY-RELATED VARIABLES
! Main arrays
type(Cell), allocatable, target :: cells(:)
type(Universe), allocatable, target :: universes(:)
type(LatticeContainer), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(ObjectPlot), allocatable, target :: plots(:)
type(Cell), allocatable, target :: cells(:)
type(Universe), allocatable, target :: universes(:)
type(LatticeContainer), allocatable, target :: lattices(:)
type(SurfaceContainer), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(ObjectPlot), allocatable, target :: plots(:)
! Size of main arrays
integer :: n_cells ! # of cells

23
src/initialize.F90
View file

@ -571,15 +571,18 @@ contains
do i = 1, n_cells
! =======================================================================
! ADJUST SURFACE LIST FOR EACH CELL
! ADJUST REGION SPECIFICATION FOR EACH CELL
c => cells(i)
do j = 1, c%n_surfaces
id = c%surfaces(j)
if (id < OP_DIFFERENCE) then
do j = 1, size(c%region)
id = c%region(j)
! Make sure that only regions are checked. Since OP_UNION is the
! operator with the lowest integer value, anything below it must denote
! a half-space
if (id < OP_UNION) then
if (surface_dict%has_key(abs(id))) then
i_array = surface_dict%get_key(abs(id))
c%surfaces(j) = sign(i_array, id)
c%region(j) = sign(i_array, id)
else
call fatal_error("Could not find surface " // trim(to_str(abs(id)))&
&// " specified on cell " // trim(to_str(c%id)))
@ -587,6 +590,16 @@ contains
end if
end do
! Also adjust the indices in the reverse Polish notation
do j = 1, size(c%rpn)
id = c%rpn(j)
! Again, make sure that only regions are checked
if (id < OP_UNION) then
i_array = surface_dict%get_key(abs(id))
c%rpn(j) = sign(i_array, id)
end if
end do
! =======================================================================
! ADJUST UNIVERSE INDEX FOR EACH CELL

276
src/input_xml.F90
View file

@ -5,15 +5,17 @@ module input_xml
use dict_header, only: DictIntInt, ElemKeyValueCI
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
use geometry_header, only: Cell, Surface, Lattice, RectLattice, HexLattice
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice
use global
use list_header, only: ListChar, ListReal
use mesh_header, only: RegularMesh
use output, only: write_message
use plot_header
use random_lcg, only: prn
use surface_header
use stl_vector, only: VectorInt
use string, only: to_lower, to_str, str_to_int, str_to_real, &
starts_with, ends_with
starts_with, ends_with, tokenize
use tally_header, only: TallyObject, TallyFilter
use tally_initialize, only: add_tallies
use xml_interface
@ -987,13 +989,15 @@ contains
integer :: coeffs_reqd
integer, allocatable :: temp_int_array(:)
real(8) :: phi, theta, psi
real(8), allocatable :: coeffs(:)
logical :: file_exists
logical :: boundary_exists
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
class(Lattice), pointer :: lat => null()
character(MAX_LINE_LEN) :: region_spec
type(Cell), pointer :: c
class(Surface), pointer :: s
class(Lattice), pointer :: lat
type(Node), pointer :: doc => null()
type(Node), pointer :: node_cell => null()
type(Node), pointer :: node_surf => null()
@ -1002,6 +1006,8 @@ contains
type(NodeList), pointer :: node_surf_list => null()
type(NodeList), pointer :: node_rlat_list => null()
type(NodeList), pointer :: node_hlat_list => null()
type(VectorInt) :: tokens
type(VectorInt) :: rpn
! Display output message
call write_message("Reading geometry XML file...", 5)
@ -1112,17 +1118,42 @@ contains
call fatal_error("Cannot specify material and fill simultaneously")
end if
! Allocate array for surfaces and copy
! Check for region specification (also under deprecated name surfaces)
region_spec = ''
if (check_for_node(node_cell, "surfaces")) then
n = get_arraysize_integer(node_cell, "surfaces")
else
n = 0
call warning("The use of 'surfaces' is deprecated and will be &
&disallowed in a future release. Use 'region' instead. The &
&openmc-update-inputs utility can be used to automatically &
&update geometry.xml files.")
call get_node_value(node_cell, "surfaces", region_spec)
elseif (check_for_node(node_cell, "region")) then
call get_node_value(node_cell, "region", region_spec)
end if
c % n_surfaces = n
if (n > 0) then
allocate(c % surfaces(n))
call get_node_array(node_cell, "surfaces", c % surfaces)
if (len_trim(region_spec) > 0) then
! Create surfaces array from string
call tokenize(region_spec, tokens)
! Use shunting-yard algorithm to determine RPN for surface algorithm
call generate_rpn(c%id, tokens, rpn)
! Copy region spec and RPN form to cell arrays
allocate(c % region(tokens%size()))
allocate(c % rpn(rpn%size()))
c % region(:) = tokens%data(1:tokens%size())
c % rpn(:) = rpn%data(1:rpn%size())
call tokens%clear()
call rpn%clear()
end if
if (.not. allocated(c%region)) allocate(c%region(0))
if (.not. allocated(c%rpn)) allocate(c%rpn(0))
! Check if this is a simple cell
if (any(c%rpn == OP_COMPLEMENT) .or. any(c%rpn == OP_UNION)) then
c%simple = .false.
else
c%simple = .true.
end if
! Rotation matrix
@ -1222,71 +1253,71 @@ contains
allocate(surfaces(n_surfaces))
do i = 1, n_surfaces
s => surfaces(i)
! Get pointer to i-th surface node
call get_list_item(node_surf_list, i, node_surf)
! Copy data into cells
if (check_for_node(node_surf, "id")) then
call get_node_value(node_surf, "id", s % id)
else
call fatal_error("Must specify id of surface in geometry XML file.")
end if
! Check to make sure 'id' hasn't been used
if (surface_dict % has_key(s % id)) then
call fatal_error("Two or more surfaces use the same unique ID: " &
&// to_str(s % id))
end if
! Copy surface name
if (check_for_node(node_surf, "name")) then
call get_node_value(node_surf, "name", s % name)
end if
! Copy and interpret surface type
word = ''
if (check_for_node(node_surf, "type")) &
call get_node_value(node_surf, "type", word)
select case(to_lower(word))
case ('x-plane')
s % type = SURF_PX
coeffs_reqd = 1
allocate(SurfaceXPlane :: surfaces(i)%obj)
case ('y-plane')
s % type = SURF_PY
coeffs_reqd = 1
allocate(SurfaceYPlane :: surfaces(i)%obj)
case ('z-plane')
s % type = SURF_PZ
coeffs_reqd = 1
allocate(SurfaceZPlane :: surfaces(i)%obj)
case ('plane')
s % type = SURF_PLANE
coeffs_reqd = 4
allocate(SurfacePlane :: surfaces(i)%obj)
case ('x-cylinder')
s % type = SURF_CYL_X
coeffs_reqd = 3
allocate(SurfaceXCylinder :: surfaces(i)%obj)
case ('y-cylinder')
s % type = SURF_CYL_Y
coeffs_reqd = 3
allocate(SurfaceYCylinder :: surfaces(i)%obj)
case ('z-cylinder')
s % type = SURF_CYL_Z
coeffs_reqd = 3
allocate(SurfaceZCylinder :: surfaces(i)%obj)
case ('sphere')
s % type = SURF_SPHERE
coeffs_reqd = 4
allocate(SurfaceSphere :: surfaces(i)%obj)
case ('x-cone')
s % type = SURF_CONE_X
coeffs_reqd = 4
allocate(SurfaceXCone :: surfaces(i)%obj)
case ('y-cone')
s % type = SURF_CONE_Y
coeffs_reqd = 4
allocate(SurfaceYCone :: surfaces(i)%obj)
case ('z-cone')
s % type = SURF_CONE_Z
coeffs_reqd = 4
allocate(SurfaceZCone :: surfaces(i)%obj)
case default
call fatal_error("Invalid surface type: " // trim(word))
end select
s => surfaces(i)%obj
! Copy data into cells
if (check_for_node(node_surf, "id")) then
call get_node_value(node_surf, "id", s%id)
else
call fatal_error("Must specify id of surface in geometry XML file.")
end if
! Check to make sure 'id' hasn't been used
if (surface_dict % has_key(s%id)) then
call fatal_error("Two or more surfaces use the same unique ID: " &
&// to_str(s%id))
end if
! Copy surface name
if (check_for_node(node_surf, "name")) then
call get_node_value(node_surf, "name", s%name)
end if
! Check to make sure that the proper number of coefficients
! have been specified for the given type of surface. Then copy
! surface coordinates.
@ -1294,36 +1325,84 @@ contains
n = get_arraysize_double(node_surf, "coeffs")
if (n < coeffs_reqd) then
call fatal_error("Not enough coefficients specified for surface: " &
&// trim(to_str(s % id)))
&// trim(to_str(s%id)))
elseif (n > coeffs_reqd) then
call fatal_error("Too many coefficients specified for surface: " &
&// trim(to_str(s % id)))
else
allocate(s % coeffs(n))
call get_node_array(node_surf, "coeffs", s % coeffs)
&// trim(to_str(s%id)))
end if
allocate(coeffs(n))
call get_node_array(node_surf, "coeffs", coeffs)
select type(s)
type is (SurfaceXPlane)
s%x0 = coeffs(1)
type is (SurfaceYPlane)
s%y0 = coeffs(1)
type is (SurfaceZPlane)
s%z0 = coeffs(1)
type is (SurfacePlane)
s%A = coeffs(1)
s%B = coeffs(2)
s%C = coeffs(3)
s%D = coeffs(4)
type is (SurfaceXCylinder)
s%y0 = coeffs(1)
s%z0 = coeffs(2)
s%r = coeffs(3)
type is (SurfaceYCylinder)
s%x0 = coeffs(1)
s%z0 = coeffs(2)
s%r = coeffs(3)
type is (SurfaceZCylinder)
s%x0 = coeffs(1)
s%y0 = coeffs(2)
s%r = coeffs(3)
type is (SurfaceSphere)
s%x0 = coeffs(1)
s%y0 = coeffs(2)
s%z0 = coeffs(3)
s%r = coeffs(4)
type is (SurfaceXCone)
s%x0 = coeffs(1)
s%y0 = coeffs(2)
s%z0 = coeffs(3)
s%r2 = coeffs(4)
type is (SurfaceYCone)
s%x0 = coeffs(1)
s%y0 = coeffs(2)
s%z0 = coeffs(3)
s%r2 = coeffs(4)
type is (SurfaceZCone)
s%x0 = coeffs(1)
s%y0 = coeffs(2)
s%z0 = coeffs(3)
s%r2 = coeffs(4)
end select
! No longer need coefficients
deallocate(coeffs)
! Boundary conditions
word = ''
if (check_for_node(node_surf, "boundary")) &
call get_node_value(node_surf, "boundary", word)
select case (to_lower(word))
case ('transmission', 'transmit', '')
s % bc = BC_TRANSMIT
s%bc = BC_TRANSMIT
case ('vacuum')
s % bc = BC_VACUUM
s%bc = BC_VACUUM
boundary_exists = .true.
case ('reflective', 'reflect', 'reflecting')
s % bc = BC_REFLECT
s%bc = BC_REFLECT
boundary_exists = .true.
case default
call fatal_error("Unknown boundary condition '" // trim(word) // &
&"' specified on surface " // trim(to_str(s % id)))
&"' specified on surface " // trim(to_str(s%id)))
end select
! Add surface to dictionary
call surface_dict % add_key(s % id, i)
call surface_dict % add_key(s%id, i)
end do
! Check to make sure a boundary condition was applied to at least one
@ -4885,4 +4964,91 @@ contains
end subroutine expand_natural_element
!===============================================================================
! GENERATE_RPN implements the shunting-yard algorithm to generate a Reverse
! Polish notation (RPN) expression for the region specification of a cell given
! the infix notation.
!===============================================================================
subroutine generate_rpn(cell_id, tokens, output)
integer, intent(in) :: cell_id
type(VectorInt), intent(in) :: tokens ! infix notation
type(VectorInt), intent(inout) :: output ! RPN notation
integer :: i
integer :: token
integer :: op
type(VectorInt) :: stack
do i = 1, tokens%size()
token = tokens%data(i)
if (token < OP_UNION) then
! If token is not an operator, add it to output
call output%push_back(token)
elseif (token < OP_RIGHT_PAREN) then
! Regular operators union, intersection, complement
do while (stack%size() > 0)
op = stack%data(stack%size())
if (op < OP_RIGHT_PAREN .and. &
((token == OP_COMPLEMENT .and. token < op) .or. &
(token /= OP_COMPLEMENT .and. token <= op))) then
! While there is an operator, op, on top of the stack, if the token
! is left-associative and its precedence is less than or equal to
! that of op or if the token is right-associative and its precedence
! is less than that of op, move op to the output queue and push the
! token on to the stack. Note that only complement is
! right-associative.
call output%push_back(op)
call stack%pop_back()
else
exit
end if
end do
call stack%push_back(token)
elseif (token == OP_LEFT_PAREN) then
! If the token is a left parenthesis, push it onto the stack
call stack%push_back(token)
else
! If the token is a right parenthesis, move operators from the stack to
! the output queue until reaching the left parenthesis.
do
! If we run out of operators without finding a left parenthesis, it
! means there are mismatched parentheses.
if (stack%size() == 0) then
call fatal_error('Mimatched parentheses in region specification &
&for cell ' // trim(to_str(cell_id)) // '.')
end if
op = stack%data(stack%size())
if (op == OP_LEFT_PAREN) exit
call output%push_back(op)
call stack%pop_back()
end do
! Pop the left parenthesis.
call stack%pop_back()
end if
end do
! While there are operators on the stack, move them to the output queue
do while (stack%size() > 0)
op = stack%data(stack%size())
! If the operator is a parenthesis, it is mismatched
if (op >= OP_RIGHT_PAREN) then
call fatal_error('Mimatched parentheses in region specification &
&for cell ' // trim(to_str(cell_id)) // '.')
end if
call output%push_back(op)
call stack%pop_back()
end do
end subroutine generate_rpn
end module input_xml

14
src/output.F90
View file

@ -6,7 +6,7 @@ module output
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning
use geometry_header, only: Cell, Universe, Surface, Lattice, RectLattice, &
use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
HexLattice, BASE_UNIVERSE
use global
use math, only: t_percentile
@ -254,10 +254,9 @@ contains
type(Particle), intent(in) :: p
integer :: i ! index for coordinate levels
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
class(Lattice), pointer :: l => null()
type(Cell), pointer :: c
type(Universe), pointer :: u
class(Lattice), pointer :: l
! display type of particle
select case (p % type)
@ -304,8 +303,7 @@ contains
! Print surface
if (p % surface /= NONE) then
s => surfaces(abs(p % surface))
write(ou,*) ' Surface = ' // to_str(sign(s % id, p % surface))
write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%obj%id, p % surface))
end if
! Display weight, energy, grid index, and interpolation factor
@ -1388,7 +1386,7 @@ contains
univ, bin-1, offset, label)
case (FILTER_SURFACE)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(surfaces(i) % id)
label = to_str(surfaces(i)%obj%id)
case (FILTER_MESH)
m => meshes(t % filters(i_filter) % int_bins(1))
allocate(ijk(m % n_dimension))

2
src/relaxng/geometry.rnc
View file

@ -9,7 +9,7 @@ element geometry {
(element material { ( xsd:int | "void" ) } |
attribute material { ( xsd:int | "void" ) })
) &
(element surfaces { list { xsd:int* } } | attribute surfaces { list { xsd:int* } })? &
(element region { xsd:string } | attribute region { xsd:string })? &
(element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? &
(element translation { list { xsd:double+ } } | attribute translation { list { xsd:double+ } })?
}*

16
src/relaxng/geometry.rng
View file

@ -62,19 +62,11 @@
</choice>
<optional>
<choice>
<element name="surfaces">
<list>
<zeroOrMore>
<data type="int"/>
</zeroOrMore>
</list>
<element name="region">
<data type="string"/>
</element>
<attribute name="surfaces">
<list>
<zeroOrMore>
<data type="int"/>
</zeroOrMore>
</list>
<attribute name="region">
<data type="string"/>
</attribute>
</choice>
</optional>

351
src/stl_vector.F90 Normal file
View file

@ -0,0 +1,351 @@
module stl_vector
! This module provides derived types that are meant to mimic the
! std::vector<T> type in C++. The vector type has numerous advantages over
! simple arrays and linked lists in that storage can grow and shrink
! dynamically, yet it is still contiguous in memory. Vectors can be filled
! element-by-element with automatic memory allocation in amortized constant
! time. In the implementation here, we grow the vector by a factor of 1.5 each
! time the capacity is exceed.
!
! The member functions which have been implemented here are:
!
! capacity -- Returns the size of the storage space currently allocated for
! the vector
!
! clear -- Remove all elements from the vector, leaving it with a size of
! 0. Note that this doesn't imply that storage is deallocated.
!
! initialize -- Set the storage size of the vector and optionally fill it with
! a particular value.
!
! pop_back -- Remove the last element of the vector, reducing the size by one.
!
! push_back -- Add a new element at the end of the vector. This increases the
! size of the vector by one. Note that the underlying storage is
! reallocated only if the size exceeds the capacity.
!
! reserve -- Requests that the capacity of the vector be a certain size.
!
! resize -- Resize the vector so it contains n elements. If n is larger than
! the current size, an optional fill value can be used to set the
! extra elements.
!
! shrink_to_fit -- Request that the capacity be reduced to fit the size.
!
! size -- Returns the number of elements in the vector.
implicit none
private
real(8), parameter :: GROWTH_FACTOR = 1.5
type, public :: VectorInt
integer, private :: size_ = 0
integer, private :: capacity_ = 0
integer, allocatable :: data(:)
contains
procedure :: capacity => capacity_int
procedure :: clear => clear_int
generic :: initialize => &
initialize_fill_int
procedure, private :: initialize_fill_int
procedure :: pop_back => pop_back_int
procedure :: push_back => push_back_int
procedure :: reserve => reserve_int
procedure :: resize => resize_int
procedure :: shrink_to_fit => shrink_to_fit_int
procedure :: size => size_int
end type VectorInt
type, public :: VectorReal
integer, private :: size_ = 0
integer, private :: capacity_ = 0
real(8), allocatable :: data(:)
contains
procedure :: capacity => capacity_real
procedure :: clear => clear_real
generic :: initialize => &
initialize_fill_real
procedure, private :: initialize_fill_real
procedure :: pop_back => pop_back_real
procedure :: push_back => push_back_real
procedure :: reserve => reserve_real
procedure :: resize => resize_real
procedure :: shrink_to_fit => shrink_to_fit_real
procedure :: size => size_real
end type VectorReal
contains
!===============================================================================
! Implementation of VectorInt
!===============================================================================
pure function capacity_int(this) result(capacity)
class(VectorInt), intent(in) :: this
integer :: capacity
capacity = this%capacity_
end function capacity_int
subroutine clear_int(this)
class(VectorInt), intent(inout) :: this
! Since integer is trivially destructible, we only need to set size to zero
! and can leave capacity as is
this%size_ = 0
end subroutine clear_int
subroutine initialize_fill_int(this, n, val)
class(VectorInt), intent(inout) :: this
integer, intent(in) :: n
integer, optional, intent(in) :: val
integer :: val_
! If no value given, fill the vector with zeros
if (present(val)) then
val_ = val
else
val_ = 0
end if
if (allocated(this%data)) deallocate(this%data)
allocate(this%data(n), SOURCE=val_)
this%size_ = n
this%capacity_ = n
end subroutine initialize_fill_int
subroutine pop_back_int(this)
class(VectorInt), intent(inout) :: this
if (this%size_ > 0) this%size_ = this%size_ - 1
end subroutine pop_back_int
subroutine push_back_int(this, val)
class(VectorInt), intent(inout) :: this
integer, intent(in) :: val
integer :: capacity
integer, allocatable :: data(:)
if (this%capacity_ == this%size_) then
! Create new data array that is GROWTH_FACTOR larger. Note that
if (this%capacity_ == 0) then
capacity = 8
else
capacity = int(GROWTH_FACTOR*this%capacity_)
end if
allocate(data(capacity))
! Copy existing elements
if (this%size_ > 0) data(1:this%size_) = this%data
! Move allocation
call move_alloc(FROM=data, TO=this%data)
this%capacity_ = capacity
end if
! Increase size of vector by one and set new element
this%size_ = this%size_ + 1
this%data(this%size_) = val
end subroutine push_back_int
subroutine reserve_int(this, n)
class(VectorInt), intent(inout) :: this
integer, intent(in) :: n
integer, allocatable :: data(:)
if (n > this%capacity_) then
allocate(data(n))
! Copy existing elements
if (this%size_ > 0) data(1:this%size_) = this%data(1:this%size_)
! Move allocation
call move_alloc(FROM=data, TO=this%data)
this%capacity_ = n
end if
end subroutine reserve_int
subroutine resize_int(this, n, val)
class(VectorInt), intent(inout) :: this
integer, intent(in) :: n
integer, intent(in), optional :: val
if (n < this%size_) then
this%size_ = n
elseif (n > this%size_) then
! If requested size is greater than capacity, first reserve that many
! elements
if (n > this%capacity_) call this%reserve(n)
! Fill added elements with specified value and increase size
if (present(val)) this%data(this%size_ + 1 : n) = val
this%size_ = n
end if
end subroutine resize_int
subroutine shrink_to_fit_int(this)
class(VectorInt), intent(inout) :: this
integer, allocatable :: data(:)
if (this%capacity_ > this%size_) then
if (this%size_ > 0) then
allocate(data(this%size_))
data(:) = this%data(1:this%size_)
call move_alloc(FROM=data, TO=this%data)
this%capacity_ = this%size_
else
if (allocated(this%data)) deallocate(this%data)
end if
end if
end subroutine shrink_to_fit_int
pure function size_int(this) result(size)
class(VectorInt), intent(in) :: this
integer :: size
size = this%size_
end function size_int
!===============================================================================
! Implementation of VectorReal
!===============================================================================
pure function capacity_real(this) result(capacity)
class(VectorReal), intent(in) :: this
integer :: capacity
capacity = this%capacity_
end function capacity_real
subroutine clear_real(this)
class(VectorReal), intent(inout) :: this
! Since real is trivially destructible, we only need to set size to zero and
! can leave capacity as is
this%size_ = 0
end subroutine clear_real
subroutine initialize_fill_real(this, n, val)
class(VectorReal), intent(inout) :: this
integer, intent(in) :: n
real(8), optional, intent(in) :: val
real(8) :: val_
! If no value given, fill the vector with zeros
if (present(val)) then
val_ = val
else
val_ = 0
end if
if (allocated(this%data)) deallocate(this%data)
allocate(this%data(n), SOURCE=val_)
this%size_ = n
this%capacity_ = n
end subroutine initialize_fill_real
subroutine pop_back_real(this)
class(VectorReal), intent(inout) :: this
if (this%size_ > 0) this%size_ = this%size_ - 1
end subroutine pop_back_real
subroutine push_back_real(this, val)
class(VectorReal), intent(inout) :: this
real(8), intent(in) :: val
integer :: capacity
real(8), allocatable :: data(:)
if (this%capacity_ == this%size_) then
! Create new data array that is GROWTH_FACTOR larger. Note that
if (this%capacity_ == 0) then
capacity = 8
else
capacity = int(GROWTH_FACTOR*this%capacity_)
end if
allocate(data(capacity))
! Copy existing elements
if (this%size_ > 0) data(1:this%size_) = this%data
! Move allocation
call move_alloc(FROM=data, TO=this%data)
this%capacity_ = capacity
end if
! Increase size of vector by one and set new element
this%size_ = this%size_ + 1
this%data(this%size_) = val
end subroutine push_back_real
subroutine reserve_real(this, n)
class(VectorReal), intent(inout) :: this
integer, intent(in) :: n
real(8), allocatable :: data(:)
if (n > this%capacity_) then
allocate(data(n))
! Copy existing elements
if (this%size_ > 0) data(1:this%size_) = this%data(1:this%size_)
! Move allocation
call move_alloc(FROM=data, TO=this%data)
this%capacity_ = n
end if
end subroutine reserve_real
subroutine resize_real(this, n, val)
class(VectorReal), intent(inout) :: this
integer, intent(in) :: n
real(8), intent(in), optional :: val
if (n < this%size_) then
this%size_ = n
elseif (n > this%size_) then
! If requested size is greater than capacity, first reserve that many
! elements
if (n > this%capacity_) call this%reserve(n)
! Fill added elements with specified value and increase size
if (present(val)) this%data(this%size_ + 1 : n) = val
this%size_ = n
end if
end subroutine resize_real
subroutine shrink_to_fit_real(this)
class(VectorReal), intent(inout) :: this
real(8), allocatable :: data(:)
if (this%capacity_ > this%size_) then
if (this%size_ > 0) then
allocate(data(this%size_))
data(:) = this%data(1:this%size_)
call move_alloc(FROM=data, TO=this%data)
this%capacity_ = this%size_
else
if (allocated(this%data)) deallocate(this%data)
end if
end if
end subroutine shrink_to_fit_real
pure function size_real(this) result(size)
class(VectorReal), intent(in) :: this
integer :: size
size = this%size_
end function size_real
end module stl_vector

125
src/string.F90
View file

@ -1,8 +1,10 @@
module string
use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL
use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL, &
OP_LEFT_PAREN, OP_RIGHT_PAREN, OP_COMPLEMENT, OP_INTERSECTION, OP_UNION
use error, only: fatal_error, warning
use global, only: master
use stl_vector, only: VectorInt
implicit none
@ -63,63 +65,96 @@ contains
end subroutine split_string
!===============================================================================
! SPLIT_STRING_WL takes a string that includes logical expressions for a list of
! bounding surfaces in a cell and splits it into separate words. The characters
! (, ), :, and # count as separate words since they represent operators.
!
! Arguments:
! string = input line
! words = array of words
! n = number of words
! TOKENIZE takes a string that includes logical expressions for a list of
! bounding surfaces in a cell and splits it into separate tokens. The characters
! (, ), |, and ~ count as separate tokens since they represent operators.
!===============================================================================
subroutine split_string_wl(string, words, n)
subroutine tokenize(string, tokens)
character(*), intent(in) :: string
type(VectorInt), intent(inout) :: tokens
character(*), intent(in) :: string
character(*), intent(out) :: words(MAX_WORDS)
integer, intent(out) :: n
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: token
character(len=len_trim(string)) :: string_
character(1) :: chr ! current character
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: i_end ! ending index of word
! Remove leading blanks
string_ = adjustl(string)
i_start = 0
i_end = 0
n = 0
do i = 1, len_trim(string)
chr = string(i:i)
i = 1
do while (i <= len_trim(string_))
! Check for special characters
if (index('():#', chr) > 0) then
if (index('()|~ ', string_(i:i)) > 0) then
! If the special character appears immediately after a non-operator,
! create a token with the surface half-space
if (i_start > 0) then
i_end = i - 1
n = n + 1
words(n) = string(i_start:i_end)
call tokens%push_back(int(str_to_int(&
string_(i_start:i - 1)), 4))
end if
n = n + 1
words(n) = chr
select case (string_(i:i))
case ('(')
call tokens%push_back(OP_LEFT_PAREN)
case (')')
if (tokens%size() > 0) then
token = tokens%data(tokens%size())
if (token >= OP_UNION .and. token < OP_RIGHT_PAREN) then
call fatal_error("Right parentheses cannot follow an operator in &
&region specification: " // trim(string))
end if
end if
call tokens%push_back(OP_RIGHT_PAREN)
case ('|')
if (tokens%size() > 0) then
token = tokens%data(tokens%size())
if (.not. (token < OP_UNION .or. token == OP_RIGHT_PAREN)) then
call fatal_error("Union cannot follow an operator in region &
&specification: " // trim(string))
end if
end if
call tokens%push_back(OP_UNION)
case ('~')
call tokens%push_back(OP_COMPLEMENT)
case (' ')
! Find next non-space character
do while (string_(i+1:i+1) == ' ')
i = i + 1
end do
! If previous token is a halfspace or right parenthesis and next token
! is not a left parenthese or union operator, that implies that the
! whitespace is to be interpreted as an intersection operator
if (i_start > 0 .or. tokens%data(tokens%size()) == OP_RIGHT_PAREN) then
if (index(')|', string_(i+1:i+1)) == 0) then
call tokens%push_back(OP_INTERSECTION)
end if
end if
end select
i_start = 0
i_end = 0
cycle
else
! Check for invalid characters
if (index('-0123456789', string_(i:i)) == 0) then
call fatal_error("Invalid character '" // string_(i:i) // "' in &
&region specification.")
end if
! If we haven't yet reached the start of a word, start a new word
if (i_start == 0) i_start = i
end if
if ((i_start == 0) .and. (chr /= ' ')) then
i_start = i
end if
if (i_start > 0) then
if (chr == ' ') i_end = i - 1
if (i == len_trim(string)) i_end = i
if (i_end > 0) then
n = n + 1
words(n) = string(i_start:i_end)
! reset indices
i_start = 0
i_end = 0
end if
end if
i = i + 1
end do
end subroutine split_string_wl
! If we've reached the end and we're still in a word, create a token from it
! and add it to the list
if (i_start > 0) then
call tokens%push_back(int(str_to_int(&
string_(i_start:len_trim(string_))), 4))
end if
end subroutine tokenize
!===============================================================================
! CONCATENATE takes an array of words and concatenates them together in one

120
src/summary.F90
View file

@ -3,13 +3,14 @@ module summary
use ace_header, only: Reaction, UrrData, Nuclide
use constants
use endf, only: reaction_name
use geometry_header, only: Cell, Surface, Universe, Lattice, RectLattice, &
use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
&HexLattice
use global
use hdf5_interface
use material_header, only: Material
use mesh_header, only: RegularMesh
use output, only: time_stamp
use surface_header
use string, only: to_str
use tally_header, only: TallyObject
@ -106,14 +107,15 @@ contains
integer :: i, j, k, m
integer, allocatable :: lattice_universes(:,:,:)
integer, allocatable :: surface_ids(:)
integer(HID_T) :: geom_group
integer(HID_T) :: cells_group, cell_group
integer(HID_T) :: surfaces_group, surface_group
integer(HID_T) :: universes_group, univ_group
integer(HID_T) :: lattices_group, lattice_group
real(8), allocatable :: coeffs(:)
character(MAX_LINE_LEN) :: region_spec
type(Cell), pointer :: c
type(Surface), pointer :: s
class(Surface), pointer :: s
type(Universe), pointer :: u
class(Lattice), pointer :: lat
@ -174,15 +176,25 @@ contains
end select
! Write list of bounding surfaces
if (c%n_surfaces > 0) then
allocate(surface_ids(c%n_surfaces))
do j = 1, c%n_surfaces
k = c%surfaces(j)
surface_ids(j) = sign(surfaces(abs(k))%id, k)
end do
call write_dataset(cell_group, "surfaces", surface_ids)
deallocate(surface_ids)
end if
region_spec = ""
do j = 1, size(c%region)
k = c%region(j)
select case(k)
case (OP_LEFT_PAREN)
region_spec = trim(region_spec) // " ("
case (OP_RIGHT_PAREN)
region_spec = trim(region_spec) // " )"
case (OP_COMPLEMENT)
region_spec = trim(region_spec) // " ~"
case (OP_INTERSECTION)
case (OP_UNION)
region_spec = trim(region_spec) // " |"
case default
region_spec = trim(region_spec) // " " // to_str(&
sign(surfaces(abs(k))%obj%id, k))
end select
end do
call write_dataset(cell_group, "region", adjustl(region_spec))
call close_group(cell_group)
end do CELL_LOOP
@ -197,7 +209,7 @@ contains
! Write information on each surface
SURFACE_LOOP: do i = 1, n_surfaces
s => surfaces(i)
s => surfaces(i)%obj
surface_group = create_group(surfaces_group, "surface " // &
trim(to_str(s%id)))
@ -208,33 +220,65 @@ contains
call write_dataset(surface_group, "name", s%name)
! Write surface type
select case (s%type)
case (SURF_PX)
select type (s)
type is (SurfaceXPlane)
call write_dataset(surface_group, "type", "x-plane")
case (SURF_PY)
call write_dataset(surface_group, "type", "y-plane")
case (SURF_PZ)
call write_dataset(surface_group, "type", "z-plane")
case (SURF_PLANE)
call write_dataset(surface_group, "type", "plane")
case (SURF_CYL_X)
call write_dataset(surface_group, "type", "x-cylinder")
case (SURF_CYL_Y)
call write_dataset(surface_group, "type", "y-cylinder")
case (SURF_CYL_Z)
call write_dataset(surface_group, "type", "z-cylinder")
case (SURF_SPHERE)
call write_dataset(surface_group, "type", "sphere")
case (SURF_CONE_X)
call write_dataset(surface_group, "type", "x-cone")
case (SURF_CONE_Y)
call write_dataset(surface_group, "type", "y-cone")
case (SURF_CONE_Z)
call write_dataset(surface_group, "type", "z-cone")
end select
allocate(coeffs(1))
coeffs(1) = s%x0
! Write coefficients for surface
call write_dataset(surface_group, "coefficients", s%coeffs)
type is (SurfaceYPlane)
call write_dataset(surface_group, "type", "y-plane")
allocate(coeffs(1))
coeffs(1) = s%y0
type is (SurfaceZPlane)
call write_dataset(surface_group, "type", "z-plane")
allocate(coeffs(1))
coeffs(1) = s%z0
type is (SurfacePlane)
call write_dataset(surface_group, "type", "plane")
allocate(coeffs(4))
coeffs(:) = [s%A, s%B, s%C, s%D]
type is (SurfaceXCylinder)
call write_dataset(surface_group, "type", "x-cylinder")
allocate(coeffs(3))
coeffs(:) = [s%y0, s%z0, s%r]
type is (SurfaceYCylinder)
call write_dataset(surface_group, "type", "y-cylinder")
allocate(coeffs(3))
coeffs(:) = [s%x0, s%z0, s%r]
type is (SurfaceZCylinder)
call write_dataset(surface_group, "type", "z-cylinder")
allocate(coeffs(3))
coeffs(:) = [s%x0, s%y0, s%r]
type is (SurfaceSphere)
call write_dataset(surface_group, "type", "sphere")
allocate(coeffs(4))
coeffs(:) = [s%x0, s%y0, s%z0, s%r]
type is (SurfaceXCone)
call write_dataset(surface_group, "type", "x-cone")
allocate(coeffs(4))
coeffs(:) = [s%x0, s%y0, s%z0, s%r2]
type is (SurfaceYCone)
call write_dataset(surface_group, "type", "y-cone")
allocate(coeffs(4))
coeffs(:) = [s%x0, s%y0, s%z0, s%r2]
type is (SurfaceZCone)
call write_dataset(surface_group, "type", "z-cone")
allocate(coeffs(4))
coeffs(:) = [s%x0, s%y0, s%z0, s%r2]
end select
call write_dataset(surface_group, "coefficients", coeffs)
deallocate(coeffs)
! Write boundary condition
select case (s%bc)

956
src/surface_header.F90 Normal file
View file

@ -0,0 +1,956 @@
module surface_header
use constants, only: ONE, TWO, ZERO, INFINITY, FP_COINCIDENT
implicit none
!===============================================================================
! SURFACE type defines a first- or second-order surface that can be used to
! construct closed volumes (cells)
!===============================================================================
type, abstract :: Surface
integer :: id ! Unique ID
integer, allocatable :: &
neighbor_pos(:), & ! List of cells on positive side
neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
character(len=52) :: name = "" ! User-defined name
contains
procedure :: sense
procedure :: reflect
procedure(iEvaluate), deferred :: evaluate
procedure(iDistance), deferred :: distance
procedure(iNormal), deferred :: normal
end type Surface
abstract interface
pure function iEvaluate(this, xyz) result(f)
import Surface
class(Surface), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
end function iEvaluate
pure function iDistance(this, xyz, uvw, coincident) result(d)
import Surface
class(Surface), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
end function iDistance
pure function iNormal(this, xyz) result(uvw)
import Surface
class(Surface), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
end function iNormal
end interface
!===============================================================================
! SURFACECONTAINER allows us to store an array of different types of surfaces
!===============================================================================
type :: SurfaceContainer
class(Surface), allocatable :: obj
end type SurfaceContainer
!===============================================================================
! All the derived types below are extensions of the abstract Surface type. They
! inherent the reflect() and sense() type-bound procedures and must implement
! evaluate(), distance(), and normal()
!===============================================================================
type, extends(Surface) :: SurfaceXPlane
! x = x0
real(8) :: x0
contains
procedure :: evaluate => x_plane_evaluate
procedure :: distance => x_plane_distance
procedure :: normal => x_plane_normal
end type SurfaceXPlane
type, extends(Surface) :: SurfaceYPlane
! y = y0
real(8) :: y0
contains
procedure :: evaluate => y_plane_evaluate
procedure :: distance => y_plane_distance
procedure :: normal => y_plane_normal
end type SurfaceYPlane
type, extends(Surface) :: SurfaceZPlane
! z = z0
real(8) :: z0
contains
procedure :: evaluate => z_plane_evaluate
procedure :: distance => z_plane_distance
procedure :: normal => z_plane_normal
end type SurfaceZPlane
type, extends(Surface) :: SurfacePlane
! Ax + By + Cz = D
real(8) :: A
real(8) :: B
real(8) :: C
real(8) :: D
contains
procedure :: evaluate => plane_evaluate
procedure :: distance => plane_distance
procedure :: normal => plane_normal
end type SurfacePlane
type, extends(Surface) :: SurfaceXCylinder
! (y - y0)^2 + (z - z0)^2 = R^2
real(8) :: y0
real(8) :: z0
real(8) :: r
contains
procedure :: evaluate => x_cylinder_evaluate
procedure :: distance => x_cylinder_distance
procedure :: normal => x_cylinder_normal
end type SurfaceXCylinder
type, extends(Surface) :: SurfaceYCylinder
! (x - x0)^2 + (z - z0)^2 = R^2
real(8) :: x0
real(8) :: z0
real(8) :: r
contains
procedure :: evaluate => y_cylinder_evaluate
procedure :: distance => y_cylinder_distance
procedure :: normal => y_cylinder_normal
end type SurfaceYCylinder
type, extends(Surface) :: SurfaceZCylinder
! (x - x0)^2 + (y - y0)^2 = R^2
real(8) :: x0
real(8) :: y0
real(8) :: r
contains
procedure :: evaluate => z_cylinder_evaluate
procedure :: distance => z_cylinder_distance
procedure :: normal => z_cylinder_normal
end type SurfaceZCylinder
type, extends(Surface) :: SurfaceSphere
! (x - x0)^2 + (y - y0)^2 + (z - z0)^2 = R^2
real(8) :: x0
real(8) :: y0
real(8) :: z0
real(8) :: r
contains
procedure :: evaluate => sphere_evaluate
procedure :: distance => sphere_distance
procedure :: normal => sphere_normal
end type SurfaceSphere
type, extends(Surface) :: SurfaceXCone
! (y - y0)^2 + (z - z0)^2 = R^2*(x - x0)^2
real(8) :: x0
real(8) :: y0
real(8) :: z0
real(8) :: r2
contains
procedure :: evaluate => x_cone_evaluate
procedure :: distance => x_cone_distance
procedure :: normal => x_cone_normal
end type SurfaceXCone
type, extends(Surface) :: SurfaceYCone
! (x - x0)^2 + (z - z0)^2 = R^2*(y - y0)^2
real(8) :: x0
real(8) :: y0
real(8) :: z0
real(8) :: r2
contains
procedure :: evaluate => y_cone_evaluate
procedure :: distance => y_cone_distance
procedure :: normal => y_cone_normal
end type SurfaceYCone
type, extends(Surface) :: SurfaceZCone
! (x - x0)^2 + (y - y0)^2 = R^2*(z - z0)^2
real(8) :: x0
real(8) :: y0
real(8) :: z0
real(8) :: r2
contains
procedure :: evaluate => z_cone_evaluate
procedure :: distance => z_cone_distance
procedure :: normal => z_cone_normal
end type SurfaceZCone
contains
!===============================================================================
! SENSE determines whether a point is on the 'positive' or 'negative' side of a
! surface. This routine is crucial for determining what cell a particular point
! is in. The positive side is indicated by a returned value of .true. and the
! negative side is indicated by a returned value of .false.
!===============================================================================
pure function sense(this, xyz, uvw) result(s)
class(Surface), intent(in) :: this ! surface
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical :: s ! sense of particle
real(8) :: f ! surface function evaluated at point
! Evaluate the surface equation at the particle's coordinates to determine
! which side the particle is on
f = this%evaluate(xyz)
! Check which side of surface the point is on
if (abs(f) < FP_COINCIDENT) then
! Particle may be coincident with this surface. To determine the sense, we
! look at the direction of the particle relative to the surface normal (by
! default in the positive direction) via their dot product.
s = (dot_product(uvw, this%normal(xyz)) > ZERO)
else
s = (f > ZERO)
end if
end function sense
!===============================================================================
! REFLECT determines the direction a particle will travel if it is specularly
! reflected from the surface at a given position and direction. The position is
! needed because the reflection is performed using the surface normal, which
! depends on the position for second-order surfaces.
!===============================================================================
pure subroutine reflect(this, xyz, uvw)
class(Surface), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(inout) :: uvw(3)
real(8) :: projection
real(8) :: magnitude
real(8) :: n(3)
! Construct normal vector
n(:) = this%normal(xyz)
! Determine projection of direction onto normal and squared magnitude of
! normal
projection = n(1)*uvw(1) + n(2)*uvw(2) + n(3)*uvw(3)
magnitude = n(1)*n(1) + n(2)*n(2) + n(3)*n(3)
! Reflect direction according to normal
uvw(:) = uvw - TWO*projection/magnitude * n
end subroutine reflect
!===============================================================================
! SurfaceXPlane Implementation
!===============================================================================
pure function x_plane_evaluate(this, xyz) result(f)
class(SurfaceXPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
f = xyz(1) - this%x0
end function x_plane_evaluate
pure function x_plane_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceXPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: f
f = this%x0 - xyz(1)
if (coincident .or. abs(f) < FP_COINCIDENT .or. uvw(1) == ZERO) then
d = INFINITY
else
d = f/uvw(1)
if (d < ZERO) d = INFINITY
end if
end function x_plane_distance
pure function x_plane_normal(this, xyz) result(uvw)
class(SurfaceXPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(:) = [ONE, ZERO, ZERO]
end function x_plane_normal
!===============================================================================
! SurfaceYPlane Implementation
!===============================================================================
pure function y_plane_evaluate(this, xyz) result(f)
class(SurfaceYPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
f = xyz(2) - this%y0
end function y_plane_evaluate
pure function y_plane_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceYPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: f
f = this%y0 - xyz(2)
if (coincident .or. abs(f) < FP_COINCIDENT .or. uvw(2) == ZERO) then
d = INFINITY
else
d = f/uvw(2)
if (d < ZERO) d = INFINITY
end if
end function y_plane_distance
pure function y_plane_normal(this, xyz) result(uvw)
class(SurfaceYPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(:) = [ZERO, ONE, ZERO]
end function y_plane_normal
!===============================================================================
! SurfaceZPlane Implementation
!===============================================================================
pure function z_plane_evaluate(this, xyz) result(f)
class(SurfaceZPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
f = xyz(3) - this%z0
end function z_plane_evaluate
pure function z_plane_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceZPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: f
f = this%z0 - xyz(3)
if (coincident .or. abs(f) < FP_COINCIDENT .or. uvw(3) == ZERO) then
d = INFINITY
else
d = f/uvw(3)
if (d < ZERO) d = INFINITY
end if
end function z_plane_distance
pure function z_plane_normal(this, xyz) result(uvw)
class(SurfaceZPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(:) = [ZERO, ZERO, ONE]
end function z_plane_normal
!===============================================================================
! SurfacePlane Implementation
!===============================================================================
pure function plane_evaluate(this, xyz) result(f)
class(SurfacePlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
f = this%A*xyz(1) + this%B*xyz(2) + this%C*xyz(3) - this%D
end function plane_evaluate
pure function plane_distance(this, xyz, uvw, coincident) result(d)
class(SurfacePlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: f
real(8) :: projection
f = this%A*xyz(1) + this%B*xyz(2) + this%C*xyz(3) - this%D
projection = this%A*uvw(1) + this%B*uvw(2) + this%C*uvw(3)
if (coincident .or. abs(f) < FP_COINCIDENT .or. projection == ZERO) then
d = INFINITY
else
d = -f/projection
if (d < ZERO) d = INFINITY
end if
end function plane_distance
pure function plane_normal(this, xyz) result(uvw)
class(SurfacePlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(:) = [this%A, this%B, this%C]
end function plane_normal
!===============================================================================
! SurfaceXCylinder Implementation
!===============================================================================
pure function x_cylinder_evaluate(this, xyz) result(f)
class(SurfaceXCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
real(8) :: y, z
y = xyz(2) - this%y0
z = xyz(3) - this%z0
f = y*y + z*z - this%r*this%r
end function x_cylinder_evaluate
pure function x_cylinder_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceXCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: y, z, k, a, c, quad
a = ONE - uvw(1)*uvw(1) ! v^2 + w^2
if (a == ZERO) then
d = INFINITY
else
y = xyz(2) - this%y0
z = xyz(3) - this%z0
k = y*uvw(2) + z*uvw(3)
c = y*y + z*z - this%r*this%r
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cylinder, thus one distance is positive/negative
! and the other is zero. The sign of k determines if we are facing in or
! out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cylinder, thus one distance must be negative
! and one must be positive. The positive distance will be the one with
! negative sign on sqrt(quad)
d = (-k + sqrt(quad))/a
else
! particle is outside the cylinder, thus both distances are either
! positive or negative. If positive, the smaller distance is the one
! with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d < ZERO) d = INFINITY
end if
end if
end function x_cylinder_distance
pure function x_cylinder_normal(this, xyz) result(uvw)
class(SurfaceXCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(1) = ZERO
uvw(2) = TWO*(xyz(2) - this%y0)
uvw(3) = TWO*(xyz(3) - this%z0)
end function x_cylinder_normal
!===============================================================================
! SurfaceYCylinder Implementation
!===============================================================================
pure function y_cylinder_evaluate(this, xyz) result(f)
class(SurfaceYCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
real(8) :: x, z
x = xyz(1) - this%x0
z = xyz(3) - this%z0
f = x*x + z*z - this%r*this%r
end function y_cylinder_evaluate
pure function y_cylinder_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceYCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, z, k, a, c, quad
a = ONE - uvw(2)*uvw(2) ! u^2 + w^2
if (a == ZERO) then
d = INFINITY
else
x = xyz(1) - this%x0
z = xyz(3) - this%z0
k = x*uvw(1) + z*uvw(3)
c = x*x + z*z - this%r*this%r
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cylinder, thus one distance is positive/negative
! and the other is zero. The sign of k determines if we are facing in or
! out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cylinder, thus one distance must be negative
! and one must be positive. The positive distance will be the one with
! negative sign on sqrt(quad)
d = (-k + sqrt(quad))/a
else
! particle is outside the cylinder, thus both distances are either
! positive or negative. If positive, the smaller distance is the one
! with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d < ZERO) d = INFINITY
end if
end if
end function y_cylinder_distance
pure function y_cylinder_normal(this, xyz) result(uvw)
class(SurfaceYCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(1) = TWO*(xyz(1) - this%x0)
uvw(2) = ZERO
uvw(3) = TWO*(xyz(3) - this%z0)
end function y_cylinder_normal
!===============================================================================
! SurfaceZCylinder Implementation
!===============================================================================
pure function z_cylinder_evaluate(this, xyz) result(f)
class(SurfaceZCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
real(8) :: x, y
x = xyz(1) - this%x0
y = xyz(2) - this%y0
f = x*x + y*y - this%r*this%r
end function z_cylinder_evaluate
pure function z_cylinder_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceZCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, k, a, c, quad
a = ONE - uvw(3)*uvw(3) ! u^2 + v^2
if (a == ZERO) then
d = INFINITY
else
x = xyz(1) - this%x0
y = xyz(2) - this%y0
k = x*uvw(1) + y*uvw(2)
c = x*x + y*y - this%r*this%r
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cylinder, thus one distance is positive/negative
! and the other is zero. The sign of k determines if we are facing in or
! out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cylinder, thus one distance must be negative
! and one must be positive. The positive distance will be the one with
! negative sign on sqrt(quad)
d = (-k + sqrt(quad))/a
else
! particle is outside the cylinder, thus both distances are either
! positive or negative. If positive, the smaller distance is the one
! with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d < ZERO) d = INFINITY
end if
end if
end function z_cylinder_distance
pure function z_cylinder_normal(this, xyz) result(uvw)
class(SurfaceZCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(1) = TWO*(xyz(1) - this%x0)
uvw(2) = TWO*(xyz(2) - this%y0)
uvw(3) = ZERO
end function z_cylinder_normal
!===============================================================================
! SphereImplementation
!===============================================================================
pure function sphere_evaluate(this, xyz) result(f)
class(SurfaceSphere), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
real(8) :: x, y, z
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
f = x*x + y*y + z*z - this%r*this%r
end function sphere_evaluate
pure function sphere_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceSphere), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, z, k, c, quad
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
k = x*uvw(1) + y*uvw(2) + z*uvw(3)
c = x*x + y*y + z*z - this%r*this%r
quad = k*k - c
if (quad < ZERO) then
! no intersection with sphere
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the sphere, thus one distance is positive/negative and
! the other is zero. The sign of k determines if we are facing in or out
if (k >= ZERO) then
d = INFINITY
else
d = -k + sqrt(quad)
end if
elseif (c < ZERO) then
! particle is inside the sphere, thus one distance must be negative and
! one must be positive. The positive distance will be the one with
! negative sign on sqrt(quad)
d = -k + sqrt(quad)
else
! particle is outside the sphere, thus both distances are either positive
! or negative. If positive, the smaller distance is the one with positive
! sign on sqrt(quad)
d = -k - sqrt(quad)
if (d < ZERO) d = INFINITY
end if
end function sphere_distance
pure function sphere_normal(this, xyz) result(uvw)
class(SurfaceSphere), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(:) = TWO*(xyz - [this%x0, this%y0, this%z0])
end function sphere_normal
!===============================================================================
! SurfaceXCone Implementation
!===============================================================================
pure function x_cone_evaluate(this, xyz) result(f)
class(SurfaceXCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
real(8) :: x, y, z
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
f = y*y + z*z - this%r2*x*x
end function x_cone_evaluate
pure function x_cone_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceXCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, z, k, a, b, c, quad
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
a = uvw(2)*uvw(2) + uvw(3)*uvw(3) - this%r2*uvw(1)*uvw(1)
k = y*uvw(2) + z*uvw(3) - this%r2*x*uvw(1)
c = y*y + z*z - this%r2*x*x
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cone, thus one distance is positive/negative and the
! other is zero. The sign of k determines which distance is zero and which
! is not.
if (k >= ZERO) then
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
else
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
end function x_cone_distance
pure function x_cone_normal(this, xyz) result(uvw)
class(SurfaceXCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(1) = -TWO*this%r2*(xyz(1) - this%x0)
uvw(2) = TWO*(xyz(2) - this%y0)
uvw(3) = TWO*(xyz(3) - this%z0)
end function x_cone_normal
!===============================================================================
! SurfaceYCone Implementation
!===============================================================================
pure function y_cone_evaluate(this, xyz) result(f)
class(SurfaceYCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
real(8) :: x, y, z
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
f = x*x + z*z - this%r2*y*y
end function y_cone_evaluate
pure function y_cone_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceYCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, z, k, a, b, c, quad
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
a = uvw(1)*uvw(1) + uvw(3)*uvw(3) - this%r2*uvw(2)*uvw(2)
k = x*uvw(1) + z*uvw(3) - this%r2*y*uvw(2)
c = x*x + z*z - this%r2*y*y
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cone, thus one distance is positive/negative and the
! other is zero. The sign of k determines which distance is zero and which
! is not.
if (k >= ZERO) then
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
else
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
end function y_cone_distance
pure function y_cone_normal(this, xyz) result(uvw)
class(SurfaceYCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(1) = TWO*(xyz(1) - this%x0)
uvw(2) = -TWO*this%r2*(xyz(2) - this%y0)
uvw(3) = TWO*(xyz(3) - this%z0)
end function y_cone_normal
!===============================================================================
! SurfaceZConeImplementation
!===============================================================================
pure function z_cone_evaluate(this, xyz) result(f)
class(SurfaceZCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: f
real(8) :: x, y, z
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
f = x*x + y*y - this%r2*z*z
end function z_cone_evaluate
pure function z_cone_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceZCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, z, k, a, b, c, quad
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
a = uvw(1)*uvw(1) + uvw(2)*uvw(2) - this%r2*uvw(3)*uvw(3)
k = x*uvw(1) + y*uvw(2) - this%r2*z*uvw(3)
c = x*x + y*y - this%r2*z*z
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cone, thus one distance is positive/negative and the
! other is zero. The sign of k determines which distance is zero and which
! is not.
if (k >= ZERO) then
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
else
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
end function z_cone_distance
pure function z_cone_normal(this, xyz) result(uvw)
class(SurfaceZCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8) :: uvw(3)
uvw(1) = TWO*(xyz(1) - this%x0)
uvw(2) = TWO*(xyz(2) - this%y0)
uvw(3) = -TWO*this%r2*(xyz(3) - this%z0)
end function z_cone_normal
end module surface_header

571
tests/input_set.py Normal file
View file

@ -0,0 +1,571 @@
import openmc
class InputSet(object):
def __init__(self):
self.settings = openmc.SettingsFile()
self.materials = openmc.MaterialsFile()
self.geometry = openmc.GeometryFile()
self.tallies = None
self.plots = None
def export(self):
self.settings.export_to_xml()
self.materials.export_to_xml()
self.geometry.export_to_xml()
if self.tallies is not None: self.tallies.export_to_xml()
if self.plots is not None: self.plots.export_to_xml()
def build_default_materials_and_geometry(self):
# Define materials.
fuel = openmc.Material(name='Fuel', material_id=1)
fuel.set_density('g/cm3', 10.062)
fuel.add_nuclide("U-234", 4.9476e-6)
fuel.add_nuclide("U-235", 4.8218e-4)
fuel.add_nuclide("U-236", 9.0402e-5)
fuel.add_nuclide("U-238", 2.1504e-2)
fuel.add_nuclide("Np-237", 7.3733e-6)
fuel.add_nuclide("Pu-238", 1.5148e-6)
fuel.add_nuclide("Pu-239", 1.3955e-4)
fuel.add_nuclide("Pu-240", 3.4405e-5)
fuel.add_nuclide("Pu-241", 2.1439e-5)
fuel.add_nuclide("Pu-242", 3.7422e-6)
fuel.add_nuclide("Am-241", 4.5041e-7)
fuel.add_nuclide("Am-242m", 9.2301e-9)
fuel.add_nuclide("Am-243", 4.7878e-7)
fuel.add_nuclide("Cm-242", 1.0485e-7)
fuel.add_nuclide("Cm-243", 1.4268e-9)
fuel.add_nuclide("Cm-244", 8.8756e-8)
fuel.add_nuclide("Cm-245", 3.5285e-9)
fuel.add_nuclide("Mo-95", 2.6497e-5)
fuel.add_nuclide("Tc-99", 3.2772e-5)
fuel.add_nuclide("Ru-101", 3.0742e-5)
fuel.add_nuclide("Ru-103", 2.3505e-6)
fuel.add_nuclide("Ag-109", 2.0009e-6)
fuel.add_nuclide("Xe-135", 1.0801e-8)
fuel.add_nuclide("Cs-133", 3.4612e-5)
fuel.add_nuclide("Nd-143", 2.6078e-5)
fuel.add_nuclide("Nd-145", 1.9898e-5)
fuel.add_nuclide("Sm-147", 1.6128e-6)
fuel.add_nuclide("Sm-149", 1.1627e-7)
fuel.add_nuclide("Sm-150", 7.1727e-6)
fuel.add_nuclide("Sm-151", 5.4947e-7)
fuel.add_nuclide("Sm-152", 3.0221e-6)
fuel.add_nuclide("Eu-153", 2.6209e-6)
fuel.add_nuclide("Gd-155", 1.5369e-9)
fuel.add_nuclide("O-16", 4.5737e-2)
clad = openmc.Material(name='Cladding', material_id=2)
clad.set_density('g/cm3', 5.77)
clad.add_nuclide("Zr-90", 0.5145)
clad.add_nuclide("Zr-91", 0.1122)
clad.add_nuclide("Zr-92", 0.1715)
clad.add_nuclide("Zr-94", 0.1738)
clad.add_nuclide("Zr-96", 0.0280)
cold_water = openmc.Material(name='Cold borated water', material_id=3)
cold_water.set_density('atom/b-cm', 0.07416)
cold_water.add_nuclide("H-1", 2.0)
cold_water.add_nuclide("O-16", 1.0)
cold_water.add_nuclide("B-10", 6.490e-4)
cold_water.add_nuclide("B-11", 2.689e-3)
cold_water.add_s_alpha_beta('HH2O', '71t')
hot_water = openmc.Material(name='Hot borated water', material_id=4)
hot_water.set_density('atom/b-cm', 0.06614)
hot_water.add_nuclide("H-1", 2.0)
hot_water.add_nuclide("O-16", 1.0)
hot_water.add_nuclide("B-10", 6.490e-4)
hot_water.add_nuclide("B-11", 2.689e-3)
hot_water.add_s_alpha_beta('HH2O', '71t')
rpv_steel = openmc.Material(name='Reactor pressure vessel steel',
material_id=5)
rpv_steel.set_density('g/cm3', 7.9)
rpv_steel.add_nuclide("Fe-54", 0.05437098, 'wo')
rpv_steel.add_nuclide("Fe-56", 0.88500663, 'wo')
rpv_steel.add_nuclide("Fe-57", 0.0208008, 'wo')
rpv_steel.add_nuclide("Fe-58", 0.00282159, 'wo')
rpv_steel.add_nuclide("Ni-58", 0.0067198, 'wo')
rpv_steel.add_nuclide("Ni-60", 0.0026776, 'wo')
rpv_steel.add_nuclide("Ni-61", 0.0001183, 'wo')
rpv_steel.add_nuclide("Ni-62", 0.0003835, 'wo')
rpv_steel.add_nuclide("Ni-64", 0.0001008, 'wo')
rpv_steel.add_nuclide("Mn-55", 0.01, 'wo')
rpv_steel.add_nuclide("Mo-92", 0.000849, 'wo')
rpv_steel.add_nuclide("Mo-94", 0.0005418, 'wo')
rpv_steel.add_nuclide("Mo-95", 0.0009438, 'wo')
rpv_steel.add_nuclide("Mo-96", 0.0010002, 'wo')
rpv_steel.add_nuclide("Mo-97", 0.0005796, 'wo')
rpv_steel.add_nuclide("Mo-98", 0.0014814, 'wo')
rpv_steel.add_nuclide("Mo-100", 0.0006042, 'wo')
rpv_steel.add_nuclide("Si-28", 0.00367464, 'wo')
rpv_steel.add_nuclide("Si-29", 0.00019336, 'wo')
rpv_steel.add_nuclide("Si-30", 0.000132, 'wo')
rpv_steel.add_nuclide("Cr-50", 0.00010435, 'wo')
rpv_steel.add_nuclide("Cr-52", 0.002092475, 'wo')
rpv_steel.add_nuclide("Cr-53", 0.00024185, 'wo')
rpv_steel.add_nuclide("Cr-54", 6.1325e-05, 'wo')
rpv_steel.add_nuclide("C-Nat", 0.0025, 'wo')
rpv_steel.add_nuclide("Cu-63", 0.0013696, 'wo')
rpv_steel.add_nuclide("Cu-65", 0.0006304, 'wo')
lower_rad_ref = openmc.Material(name='Lower radial reflector',
material_id=6)
lower_rad_ref.set_density('g/cm3', 4.32)
lower_rad_ref.add_nuclide("H-1", 0.0095661, 'wo')
lower_rad_ref.add_nuclide("O-16", 0.0759107, 'wo')
lower_rad_ref.add_nuclide("B-10", 3.08409e-5, 'wo')
lower_rad_ref.add_nuclide("B-11", 1.40499e-4, 'wo')
lower_rad_ref.add_nuclide("Fe-54", 0.035620772088, 'wo')
lower_rad_ref.add_nuclide("Fe-56", 0.579805982228, 'wo')
lower_rad_ref.add_nuclide("Fe-57", 0.01362750048, 'wo')
lower_rad_ref.add_nuclide("Fe-58", 0.001848545204, 'wo')
lower_rad_ref.add_nuclide("Ni-58", 0.055298376566, 'wo')
lower_rad_ref.add_nuclide("Ni-60", 0.022034425592, 'wo')
lower_rad_ref.add_nuclide("Ni-61", 0.000973510811, 'wo')
lower_rad_ref.add_nuclide("Ni-62", 0.003155886695, 'wo')
lower_rad_ref.add_nuclide("Ni-64", 0.000829500336, 'wo')
lower_rad_ref.add_nuclide("Mn-55", 0.0182870, 'wo')
lower_rad_ref.add_nuclide("Si-28", 0.00839976771, 'wo')
lower_rad_ref.add_nuclide("Si-29", 0.00044199679, 'wo')
lower_rad_ref.add_nuclide("Si-30", 0.0003017355, 'wo')
lower_rad_ref.add_nuclide("Cr-50", 0.007251360806, 'wo')
lower_rad_ref.add_nuclide("Cr-52", 0.145407678031, 'wo')
lower_rad_ref.add_nuclide("Cr-53", 0.016806340306, 'wo')
lower_rad_ref.add_nuclide("Cr-54", 0.004261520857, 'wo')
lower_rad_ref.add_s_alpha_beta('HH2O', '71t')
upper_rad_ref = openmc.Material(name='Upper radial reflector /'
'Top plate region', material_id=7)
upper_rad_ref.set_density('g/cm3', 4.28)
upper_rad_ref.add_nuclide("H-1", 0.0086117, 'wo')
upper_rad_ref.add_nuclide("O-16", 0.0683369, 'wo')
upper_rad_ref.add_nuclide("B-10", 2.77638e-5, 'wo')
upper_rad_ref.add_nuclide("B-11", 1.26481e-4, 'wo')
upper_rad_ref.add_nuclide("Fe-54", 0.035953677186, 'wo')
upper_rad_ref.add_nuclide("Fe-56", 0.585224740891, 'wo')
upper_rad_ref.add_nuclide("Fe-57", 0.01375486056, 'wo')
upper_rad_ref.add_nuclide("Fe-58", 0.001865821363, 'wo')
upper_rad_ref.add_nuclide("Ni-58", 0.055815129186, 'wo')
upper_rad_ref.add_nuclide("Ni-60", 0.022240333032, 'wo')
upper_rad_ref.add_nuclide("Ni-61", 0.000982608081, 'wo')
upper_rad_ref.add_nuclide("Ni-62", 0.003185377845, 'wo')
upper_rad_ref.add_nuclide("Ni-64", 0.000837251856, 'wo')
upper_rad_ref.add_nuclide("Mn-55", 0.0184579, 'wo')
upper_rad_ref.add_nuclide("Si-28", 0.00847831314, 'wo')
upper_rad_ref.add_nuclide("Si-29", 0.00044612986, 'wo')
upper_rad_ref.add_nuclide("Si-30", 0.000304557, 'wo')
upper_rad_ref.add_nuclide("Cr-50", 0.00731912987, 'wo')
upper_rad_ref.add_nuclide("Cr-52", 0.146766614995, 'wo')
upper_rad_ref.add_nuclide("Cr-53", 0.01696340737, 'wo')
upper_rad_ref.add_nuclide("Cr-54", 0.004301347765, 'wo')
upper_rad_ref.add_s_alpha_beta('HH2O', '71t')
bot_plate = openmc.Material(name='Bottom plate region', material_id=8)
bot_plate.set_density('g/cm3', 7.184)
bot_plate.add_nuclide("H-1", 0.0011505, 'wo')
bot_plate.add_nuclide("O-16", 0.0091296, 'wo')
bot_plate.add_nuclide("B-10", 3.70915e-6, 'wo')
bot_plate.add_nuclide("B-11", 1.68974e-5, 'wo')
bot_plate.add_nuclide("Fe-54", 0.03855611055, 'wo')
bot_plate.add_nuclide("Fe-56", 0.627585036425, 'wo')
bot_plate.add_nuclide("Fe-57", 0.014750478, 'wo')
bot_plate.add_nuclide("Fe-58", 0.002000875025, 'wo')
bot_plate.add_nuclide("Ni-58", 0.059855207342, 'wo')
bot_plate.add_nuclide("Ni-60", 0.023850159704, 'wo')
bot_plate.add_nuclide("Ni-61", 0.001053732407, 'wo')
bot_plate.add_nuclide("Ni-62", 0.003415945715, 'wo')
bot_plate.add_nuclide("Ni-64", 0.000897854832, 'wo')
bot_plate.add_nuclide("Mn-55", 0.0197940, 'wo')
bot_plate.add_nuclide("Si-28", 0.00909197802, 'wo')
bot_plate.add_nuclide("Si-29", 0.00047842098, 'wo')
bot_plate.add_nuclide("Si-30", 0.000326601, 'wo')
bot_plate.add_nuclide("Cr-50", 0.007848910646, 'wo')
bot_plate.add_nuclide("Cr-52", 0.157390026871, 'wo')
bot_plate.add_nuclide("Cr-53", 0.018191270146, 'wo')
bot_plate.add_nuclide("Cr-54", 0.004612692337, 'wo')
bot_plate.add_s_alpha_beta('HH2O', '71t')
bot_nozzle = openmc.Material(name='Bottom nozzle region', material_id=9)
bot_nozzle.set_density('g/cm3', 2.53)
bot_nozzle.add_nuclide("H-1", 0.0245014, 'wo')
bot_nozzle.add_nuclide("O-16", 0.1944274, 'wo')
bot_nozzle.add_nuclide("B-10", 7.89917e-5, 'wo')
bot_nozzle.add_nuclide("B-11", 3.59854e-4, 'wo')
bot_nozzle.add_nuclide("Fe-54", 0.030411411144, 'wo')
bot_nozzle.add_nuclide("Fe-56", 0.495012237964, 'wo')
bot_nozzle.add_nuclide("Fe-57", 0.01163454624, 'wo')
bot_nozzle.add_nuclide("Fe-58", 0.001578204652, 'wo')
bot_nozzle.add_nuclide("Ni-58", 0.047211231662, 'wo')
bot_nozzle.add_nuclide("Ni-60", 0.018811987544, 'wo')
bot_nozzle.add_nuclide("Ni-61", 0.000831139127, 'wo')
bot_nozzle.add_nuclide("Ni-62", 0.002694352115, 'wo')
bot_nozzle.add_nuclide("Ni-64", 0.000708189552, 'wo')
bot_nozzle.add_nuclide("Mn-55", 0.0156126, 'wo')
bot_nozzle.add_nuclide("Si-28", 0.007171335558, 'wo')
bot_nozzle.add_nuclide("Si-29", 0.000377356542, 'wo')
bot_nozzle.add_nuclide("Si-30", 0.0002576079, 'wo')
bot_nozzle.add_nuclide("Cr-50", 0.006190885148, 'wo')
bot_nozzle.add_nuclide("Cr-52", 0.124142524198, 'wo')
bot_nozzle.add_nuclide("Cr-53", 0.014348496148, 'wo')
bot_nozzle.add_nuclide("Cr-54", 0.003638294506, 'wo')
bot_nozzle.add_s_alpha_beta('HH2O', '71t')
top_nozzle = openmc.Material(name='Top nozzle region', material_id=10)
top_nozzle.set_density('g/cm3', 1.746)
top_nozzle.add_nuclide("H-1", 0.0358870, 'wo')
top_nozzle.add_nuclide("O-16", 0.2847761, 'wo')
top_nozzle.add_nuclide("B-10", 1.15699e-4, 'wo')
top_nozzle.add_nuclide("B-11", 5.27075e-4, 'wo')
top_nozzle.add_nuclide("Fe-54", 0.02644016154, 'wo')
top_nozzle.add_nuclide("Fe-56", 0.43037146399, 'wo')
top_nozzle.add_nuclide("Fe-57", 0.0101152584, 'wo')
top_nozzle.add_nuclide("Fe-58", 0.00137211607, 'wo')
top_nozzle.add_nuclide("Ni-58", 0.04104621835, 'wo')
top_nozzle.add_nuclide("Ni-60", 0.0163554502, 'wo')
top_nozzle.add_nuclide("Ni-61", 0.000722605975, 'wo')
top_nozzle.add_nuclide("Ni-62", 0.002342513875, 'wo')
top_nozzle.add_nuclide("Ni-64", 0.0006157116, 'wo')
top_nozzle.add_nuclide("Mn-55", 0.0135739, 'wo')
top_nozzle.add_nuclide("Si-28", 0.006234853554, 'wo')
top_nozzle.add_nuclide("Si-29", 0.000328078746, 'wo')
top_nozzle.add_nuclide("Si-30", 0.0002239677, 'wo')
top_nozzle.add_nuclide("Cr-50", 0.005382452306, 'wo')
top_nozzle.add_nuclide("Cr-52", 0.107931450781, 'wo')
top_nozzle.add_nuclide("Cr-53", 0.012474806806, 'wo')
top_nozzle.add_nuclide("Cr-54", 0.003163190107, 'wo')
top_nozzle.add_s_alpha_beta('HH2O', '71t')
top_fa = openmc.Material(name='Top of fuel assemblies', material_id=11)
top_fa.set_density('g/cm3', 3.044)
top_fa.add_nuclide("H-1", 0.0162913, 'wo')
top_fa.add_nuclide("O-16", 0.1292776, 'wo')
top_fa.add_nuclide("B-10", 5.25228e-5, 'wo')
top_fa.add_nuclide("B-11", 2.39272e-4, 'wo')
top_fa.add_nuclide("Zr-90", 0.43313403903, 'wo')
top_fa.add_nuclide("Zr-91", 0.09549277374, 'wo')
top_fa.add_nuclide("Zr-92", 0.14759527104, 'wo')
top_fa.add_nuclide("Zr-94", 0.15280552077, 'wo')
top_fa.add_nuclide("Zr-96", 0.02511169542, 'wo')
top_fa.add_s_alpha_beta('HH2O', '71t')
bot_fa = openmc.Material(name='Bottom of fuel assemblies', material_id=12)
bot_fa.set_density('g/cm3', 1.762)
bot_fa.add_nuclide("H-1", 0.0292856, 'wo')
bot_fa.add_nuclide("O-16", 0.2323919, 'wo')
bot_fa.add_nuclide("B-10", 9.44159e-5, 'wo')
bot_fa.add_nuclide("B-11", 4.30120e-4, 'wo')
bot_fa.add_nuclide("Zr-90", 0.3741373658, 'wo')
bot_fa.add_nuclide("Zr-91", 0.0824858164, 'wo')
bot_fa.add_nuclide("Zr-92", 0.1274914944, 'wo')
bot_fa.add_nuclide("Zr-94", 0.1319920622, 'wo')
bot_fa.add_nuclide("Zr-96", 0.0216912612, 'wo')
bot_fa.add_s_alpha_beta('HH2O', '71t')
# Define the materials file.
self.materials.default_xs = '71c'
self.materials.add_materials((fuel, clad, cold_water, hot_water,
rpv_steel, lower_rad_ref, upper_rad_ref, bot_plate, bot_nozzle,
top_nozzle, top_fa, bot_fa))
# Define surfaces.
s1 = openmc.ZCylinder(R=0.41, surface_id=1)
s2 = openmc.ZCylinder(R=0.475, surface_id=2)
s3 = openmc.ZCylinder(R=0.56, surface_id=3)
s4 = openmc.ZCylinder(R=0.62, surface_id=4)
s5 = openmc.ZCylinder(R=187.6, surface_id=5)
s6 = openmc.ZCylinder(R=209.0, surface_id=6)
s7 = openmc.ZCylinder(R=229.0, surface_id=7)
s8 = openmc.ZCylinder(R=249.0, surface_id=8)
s8.boundary_type = 'vacuum'
s31 = openmc.ZPlane(z0=-229.0, surface_id=31)
s31.boundary_type = 'vacuum'
s32 = openmc.ZPlane(z0=-199.0, surface_id=32)
s33 = openmc.ZPlane(z0=-193.0, surface_id=33)
s34 = openmc.ZPlane(z0=-183.0, surface_id=34)
s35 = openmc.ZPlane(z0=0.0, surface_id=35)
s36 = openmc.ZPlane(z0=183.0, surface_id=36)
s37 = openmc.ZPlane(z0=203.0, surface_id=37)
s38 = openmc.ZPlane(z0=215.0, surface_id=38)
s39 = openmc.ZPlane(z0=223.0, surface_id=39)
s39.boundary_type = 'vacuum'
# Define pin cells.
fuel_cold = openmc.Universe(name='Fuel pin, cladding, cold water',
universe_id=1)
c21 = openmc.Cell(cell_id=21)
c21.region = -s1
c21.fill = fuel
c22 = openmc.Cell(cell_id=22)
c22.region = +s1 & -s2
c22.fill = clad
c23 = openmc.Cell(cell_id=23)
c23.region = +s2
c23.fill = cold_water
fuel_cold.add_cells((c21, c22, c23))
tube_cold = openmc.Universe(name='Instrumentation guide tube, '
'cold water', universe_id=2)
c24 = openmc.Cell(cell_id=24)
c24.region = -s3
c24.fill = cold_water
c25 = openmc.Cell(cell_id=25)
c25.region = +s3 & -s4
c25.fill = clad
c26 = openmc.Cell(cell_id=26)
c26.region = +s4
c26.fill = cold_water
tube_cold.add_cells((c24, c25, c26))
fuel_hot = openmc.Universe(name='Fuel pin, cladding, hot water',
universe_id=3)
c27 = openmc.Cell(cell_id=27)
c27.region = -s1
c27.fill = fuel
c28 = openmc.Cell(cell_id=28)
c28.region = +s1 & -s2
c28.fill = clad
c29 = openmc.Cell(cell_id=29)
c29.region = +s2
c29.fill = hot_water
fuel_hot.add_cells((c27, c28, c29))
tube_hot = openmc.Universe(name='Instrumentation guide tube, hot water',
universe_id=4)
c30 = openmc.Cell(cell_id=30)
c30.region = -s3
c30.fill = hot_water
c31 = openmc.Cell(cell_id=31)
c31.region = +s3 & -s4
c31.fill = clad
c32 = openmc.Cell(cell_id=32)
c32.region = +s4
c32.fill = hot_water
tube_hot.add_cells((c30, c31, c32))
# Define fuel lattices.
l100 = openmc.RectLattice(name='Fuel assembly (lower half)',
lattice_id=100)
l100.dimension = (17, 17)
l100.lower_left = (-10.71, -10.71)
l100.pitch = (1.26, 1.26)
l100.universes = [
[fuel_cold]*17,
[fuel_cold]*17,
[fuel_cold]*5 + [tube_cold] + [fuel_cold]*2 + [tube_cold]
+ [fuel_cold]*2 + [tube_cold] + [fuel_cold]*5,
[fuel_cold]*3 + [tube_cold] + [fuel_cold]*9 + [tube_cold]
+ [fuel_cold]*3,
[fuel_cold]*17,
[fuel_cold]*2 + [tube_cold] + [fuel_cold]*2 + [tube_cold]
+ [fuel_cold]*2 + [tube_cold] + [fuel_cold]*2 + [tube_cold]
+ [fuel_cold]*2 + [tube_cold] + [fuel_cold]*2,
[fuel_cold]*17,
[fuel_cold]*17,
[fuel_cold]*2 + [tube_cold] + [fuel_cold]*2 + [tube_cold]
+ [fuel_cold]*2 + [tube_cold] + [fuel_cold]*2 + [tube_cold]
+ [fuel_cold]*2 + [tube_cold] + [fuel_cold]*2,
[fuel_cold]*17,
[fuel_cold]*17,
[fuel_cold]*2 + [tube_cold] + [fuel_cold]*2 + [tube_cold]
+ [fuel_cold]*2 + [tube_cold] + [fuel_cold]*2 + [tube_cold]
+ [fuel_cold]*2 + [tube_cold] + [fuel_cold]*2,
[fuel_cold]*17,
[fuel_cold]*3 + [tube_cold] + [fuel_cold]*9 + [tube_cold]
+ [fuel_cold]*3,
[fuel_cold]*5 + [tube_cold] + [fuel_cold]*2 + [tube_cold]
+ [fuel_cold]*2 + [tube_cold] + [fuel_cold]*5,
[fuel_cold]*17,
[fuel_cold]*17 ]
l101 = openmc.RectLattice(name='Fuel assembly (upper half)',
lattice_id=101)
l101.dimension = (17, 17)
l101.lower_left = (-10.71, -10.71)
l101.pitch = (1.26, 1.26)
l101.universes = [
[fuel_hot]*17,
[fuel_hot]*17,
[fuel_hot]*5 + [tube_hot] + [fuel_hot]*2 + [tube_hot]
+ [fuel_hot]*2 + [tube_hot] + [fuel_hot]*5,
[fuel_hot]*3 + [tube_hot] + [fuel_hot]*9 + [tube_hot]
+ [fuel_hot]*3,
[fuel_hot]*17,
[fuel_hot]*2 + [tube_hot] + [fuel_hot]*2 + [tube_hot]
+ [fuel_hot]*2 + [tube_hot] + [fuel_hot]*2 + [tube_hot]
+ [fuel_hot]*2 + [tube_hot] + [fuel_hot]*2,
[fuel_hot]*17,
[fuel_hot]*17,
[fuel_hot]*2 + [tube_hot] + [fuel_hot]*2 + [tube_hot]
+ [fuel_hot]*2 + [tube_hot] + [fuel_hot]*2 + [tube_hot]
+ [fuel_hot]*2 + [tube_hot] + [fuel_hot]*2,
[fuel_hot]*17,
[fuel_hot]*17,
[fuel_hot]*2 + [tube_hot] + [fuel_hot]*2 + [tube_hot]
+ [fuel_hot]*2 + [tube_hot] + [fuel_hot]*2 + [tube_hot]
+ [fuel_hot]*2 + [tube_hot] + [fuel_hot]*2,
[fuel_hot]*17,
[fuel_hot]*3 + [tube_hot] + [fuel_hot]*9 + [tube_hot]
+ [fuel_hot]*3,
[fuel_hot]*5 + [tube_hot] + [fuel_hot]*2 + [tube_hot]
+ [fuel_hot]*2 + [tube_hot] + [fuel_hot]*5,
[fuel_hot]*17,
[fuel_hot]*17 ]
# Define assemblies.
fa_cw = openmc.Universe(name='Water assembly (cold)', universe_id=5)
c50 = openmc.Cell(cell_id=50)
c50.region = +s34 & -s35
c50.fill = cold_water
fa_cw.add_cells((c50, ))
fa_hw = openmc.Universe(name='Water assembly (hot)', universe_id=7)
c70 = openmc.Cell(cell_id=70)
c70.region = +s35 & -s36
c70.fill = hot_water
fa_hw.add_cells((c70, ))
fa_cold = openmc.Universe(name='Fuel assembly (cold)', universe_id=6)
c60 = openmc.Cell(cell_id=60)
c60.region = +s34 & -s35
c60.fill = l100
fa_cold.add_cells((c60, ))
fa_hot = openmc.Universe(name='Fuel assembly (hot)', universe_id=8)
c80 = openmc.Cell(cell_id=80)
c80.region = +s35 & -s36
c80.fill = l101
fa_hot.add_cells((c80, ))
# Define core lattices
l200 = openmc.RectLattice(name='Core lattice (lower half)',
lattice_id=200)
l200.dimension = (21, 21)
l200.lower_left = (-224.91, -224.91)
l200.pitch = (21.42, 21.42)
l200.universes = [
[fa_cw]*21,
[fa_cw]*21,
[fa_cw]*7 + [fa_cold]*7 + [fa_cw]*7,
[fa_cw]*5 + [fa_cold]*11 + [fa_cw]*5,
[fa_cw]*4 + [fa_cold]*13 + [fa_cw]*4,
[fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
[fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
[fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
[fa_cw]*4 + [fa_cold]*13 + [fa_cw]*4,
[fa_cw]*5 + [fa_cold]*11 + [fa_cw]*5,
[fa_cw]*7 + [fa_cold]*7 + [fa_cw]*7,
[fa_cw]*21,
[fa_cw]*21]
l201 = openmc.RectLattice(name='Core lattice (lower half)',
lattice_id=201)
l201.dimension = (21, 21)
l201.lower_left = (-224.91, -224.91)
l201.pitch = (21.42, 21.42)
l201.universes = [
[fa_hw]*21,
[fa_hw]*21,
[fa_hw]*7 + [fa_hot]*7 + [fa_hw]*7,
[fa_hw]*5 + [fa_hot]*11 + [fa_hw]*5,
[fa_hw]*4 + [fa_hot]*13 + [fa_hw]*4,
[fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
[fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
[fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
[fa_hw]*4 + [fa_hot]*13 + [fa_hw]*4,
[fa_hw]*5 + [fa_hot]*11 + [fa_hw]*5,
[fa_hw]*7 + [fa_hot]*7 + [fa_hw]*7,
[fa_hw]*21,
[fa_hw]*21]
# Define root universe.
root = openmc.Universe(universe_id=0, name='root universe')
c1 = openmc.Cell(cell_id=1)
c1.region = -s6 & +s34 & -s35
c1.fill = l200
c2 = openmc.Cell(cell_id=2)
c2.region = -s6 & +s35 & -s36
c2.fill = l201
c3 = openmc.Cell(cell_id=3)
c3.region = -s7 & +s31 & -s32
c3.fill = bot_plate
c4 = openmc.Cell(cell_id=4)
c4.region = -s5 & +s32 & -s33
c4.fill = bot_nozzle
c5 = openmc.Cell(cell_id=5)
c5.region = -s5 & +s33 & -s34
c5.fill = bot_fa
c6 = openmc.Cell(cell_id=6)
c6.region = -s5 & +s36 & -s37
c6.fill = top_fa
c7 = openmc.Cell(cell_id=7)
c7.region = -s5 & +s37 & -s38
c7.fill = top_nozzle
c8 = openmc.Cell(cell_id=8)
c8.region = -s7 & +s38 & -s39
c8.fill = upper_rad_ref
c9 = openmc.Cell(cell_id=9)
c9.region = +s6 & -s7 & +s32 & -s38
c9.fill = bot_nozzle
c10 = openmc.Cell(cell_id=10)
c10.region = +s7 & -s8 & +s31 & -s39
c10.fill = rpv_steel
c11 = openmc.Cell(cell_id=11)
c11.region = +s5 & -s6 & +s32 & -s34
c11.fill = lower_rad_ref
c12 = openmc.Cell(cell_id=12)
c12.region = +s5 & -s6 & +s36 & -s38
c12.fill = upper_rad_ref
root.add_cells((c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12))
# Define the geometry file.
geometry = openmc.Geometry()
geometry.root_universe = root
self.geometry.geometry = geometry
def build_default_settings(self):
self.settings.batches = 10
self.settings.inactive = 5
self.settings.particles = 100
self.settings.set_source_space('box', (-160, -160, -183, 160, 160, 183))
def build_defualt_plots(self):
plot = openmc.Plot()
plot.filename = 'mat'
plot.origin = (125, 125, 0)
plot.width = (250, 250)
plot.pixels = (3000, 3000)
plot.color = 'mat'
self.plots.add_plot(plot)

2
tests/test_basic/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

3
tests/test_basic/test_basic.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness

4
tests/test_cmfd_feed/geometry.xml
View file

@ -4,7 +4,7 @@
<!-- Definition of Cells -->
<cell id="1">
<universe>0</universe>
<surfaces>-1 2 -3 4 -5 6</surfaces>
<region>-1 2 -3 4 -5 6</region>
<material>1</material>
</cell>
@ -39,5 +39,5 @@
<coeffs>-1</coeffs>
<boundary>reflective</boundary>
</surface>
</geometry>

3
tests/test_cmfd_feed/test_cmfd_feed.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import CMFDTestHarness

4
tests/test_cmfd_nofeed/geometry.xml
View file

@ -4,7 +4,7 @@
<!-- Definition of Cells -->
<cell id="1">
<universe>0</universe>
<surfaces>-1 2 -3 4 -5 6</surfaces>
<region>-1 2 -3 4 -5 6</region>
<material>1</material>
</cell>
@ -39,5 +39,5 @@
<coeffs>-1</coeffs>
<boundary>reflective</boundary>
</surface>
</geometry>

3
tests/test_cmfd_nofeed/test_cmfd_nofeed.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import CMFDTestHarness

24
tests/test_complex_cell/geometry.xml Normal file
View file

@ -0,0 +1,24 @@
<?xml version="1.0"?>
<geometry>
<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="7" type="x-plane" coeffs="0" />
<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" />
<cell id="1" material="1" region="3 -4 13 -14" />
<cell id="2" material="2" region="2 -5 12 -15 ~(3 -4 13 -14)" />
<cell id="3" material="3" region="7 -6 11 -16 (-2 | 5 | -12 | 15)" />
<cell id="4" material="4" region="((1 -7) 11 -16) ~(2 -5 (12 -15))" />
</geometry>

23
tests/test_complex_cell/materials.xml Normal file
View file

@ -0,0 +1,23 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" ao="1.0" />
</material>
<material id="2">
<density value="4.5" units="g/cc" />
<nuclide name="U-238" ao="1.0" />
</material>
<material id="3">
<density value="2.0" units="g/cc" />
<element name="Zr" ao="1.0" />
</material>
<material id="4">
<density value="0.1" units="g/cc" />
<element name="N" ao="1.0" />
</material>
</materials>

11
tests/test_complex_cell/results_true.dat Normal file
View file

@ -0,0 +1,11 @@
k-combined:
2.651570E-01 2.116381E-03
tally 1:
2.639097E+00
1.394398E+00
2.743740E+00
1.506124E+00
1.041248E+00
2.177204E-01
1.087210E-01
2.365126E-03

7
tests/test_filter_energy/settings.xml → tests/test_complex_cell/settings.xml
View file

@ -4,15 +4,12 @@
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>100</particles>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>
-160 -160 -183
160 160 183
</parameters>
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>

6
tests/test_filter_cell/tallies.xml → tests/test_complex_cell/tallies.xml
View file

@ -1,9 +1,7 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<filter type="cell" bins="10 21 22 23" />
<filter type="cell" bins="1 2 3 4" />
<scores>total</scores>
</tally>
</tallies>
</tallies>

10
tests/test_complex_cell/test_complex_cell.py Executable file
View file

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

2
tests/test_confidence_intervals/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

3
tests/test_confidence_intervals/test_confidence_intervals.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness

2
tests/test_density_atombcm/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

3
tests/test_density_atombcm/test_density_atombcm.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness

2
tests/test_density_atomcm3/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

3
tests/test_density_atomcm3/test_density_atomcm3.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness

2
tests/test_density_kgm3/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

3
tests/test_density_kgm3/test_density_kgm3.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness

2
tests/test_density_sum/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

3
tests/test_density_sum/test_density_sum.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness

2
tests/test_eigenvalue_genperbatch/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

3
tests/test_eigenvalue_genperbatch/test_eigenvalue_genperbatch.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness

2
tests/test_eigenvalue_no_inactive/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

3
tests/test_eigenvalue_no_inactive/test_eigenvalue_no_inactive.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness

2
tests/test_energy_grid/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

3
tests/test_energy_grid/test_energy_grid.py
View file

@ -1,7 +1,8 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness

2
tests/test_entropy/geometry.xml
View file

@ -3,6 +3,6 @@
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
<cell id="1" material="1" region="-1" />
</geometry>

6
tests/test_entropy/test_entropy.py
View file

@ -3,9 +3,9 @@
import glob
import os
import sys
sys.path.insert(0, '..')
from testing_harness import *
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
from openmc.statepoint import StatePoint
class EntropyTestHarness(TestHarness):

181
tests/test_filter_cell/geometry.xml
View file

@ -1,181 +0,0 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="z-cylinder" coeffs="0. 0. 0.41" />
<surface id="2" type="z-cylinder" coeffs="0. 0. 0.475" />
<surface id="3" type="z-cylinder" coeffs="0. 0. 0.56" />
<surface id="4" type="z-cylinder" coeffs="0. 0. 0.62" />
<surface id="5" type="z-cylinder" coeffs="0. 0. 187.6" />
<surface id="6" type="z-cylinder" coeffs="0. 0. 209.0" />
<surface id="7" type="z-cylinder" coeffs="0. 0. 229.0" />
<surface id="8" type="z-cylinder" coeffs="0. 0. 249.0" boundary="vacuum" />
<surface id="31" type="z-plane" coeffs="-229.0" boundary="vacuum" />
<surface id="32" type="z-plane" coeffs="-199.0" />
<surface id="33" type="z-plane" coeffs="-193.0" />
<surface id="34" type="z-plane" coeffs="-183.0" />
<surface id="35" type="z-plane" coeffs="0.0" />
<surface id="36" type="z-plane" coeffs="183.0" />
<surface id="37" type="z-plane" coeffs="203.0" />
<surface id="38" type="z-plane" coeffs="215.0" />
<surface id="39" type="z-plane" coeffs="223.0" boundary="vacuum" />
<!-- All geometry on base universe -->
<cell id="1" fill="200" surfaces=" -6 34 -35" /> <!-- Lower core -->
<cell id="2" fill="201" surfaces=" -6 35 -36" /> <!-- Upper core -->
<cell id="3" material="8" surfaces=" -7 31 -32" /> <!-- Lower core plate region -->
<cell id="4" material="9" surfaces=" -5 32 -33" /> <!-- Bottom nozzle region -->
<cell id="5" material="12" surfaces=" -5 33 -34" /> <!-- Bottom FA region -->
<cell id="6" material="11" surfaces=" -5 36 -37" /> <!-- Top FA region -->
<cell id="7" material="10" surfaces=" -5 37 -38" /> <!-- Top nozzle region -->
<cell id="8" material="7" surfaces=" -7 38 -39" /> <!-- Upper plate region -->
<cell id="9" material="4" surfaces="6 -7 32 -38" /> <!-- Downcomer -->
<cell id="10" material="5" surfaces="7 -8 31 -39" /> <!-- RPV -->
<cell id="11" material="6" surfaces="5 -6 32 -34" /> <!-- Bottom of radial reflector -->
<cell id="12" material="7" surfaces="5 -6 36 -38" /> <!-- Top of radial reflector -->
<!-- Fuel pin, cladding, cold water -->
<cell id="21" universe="1" material="1" surfaces="-1" />
<cell id="22" universe="1" material="2" surfaces="1 -2" />
<cell id="23" universe="1" material="3" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="24" universe="2" material="3" surfaces="-3" />
<cell id="25" universe="2" material="2" surfaces="3 -4" />
<cell id="26" universe="2" material="3" surfaces="4" />
<!-- Fuel pin, cladding, hot water -->
<cell id="27" universe="3" material="1" surfaces="-1" />
<cell id="28" universe="3" material="2" surfaces="1 -2" />
<cell id="29" universe="3" material="4" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="30" universe="4" material="4" surfaces="-3" />
<cell id="31" universe="4" material="2" surfaces="3 -4" />
<cell id="32" universe="4" material="4" surfaces="4" />
<!-- cell for water assembly (cold) -->
<cell id="50" universe="5" material="4" surfaces="34 -35" />
<!-- containing cell for fuel assembly -->
<cell id="60" universe="6" fill="100" surfaces="34 -35" />
<!-- cell for water assembly (hot) -->
<cell id="70" universe="7" material="3" surfaces="35 -36" />
<!-- containing cell for fuel assembly -->
<cell id="80" universe="8" fill="101" surfaces="35 -36" />
<!-- Fuel Assembly (Lower Half) -->
<lattice id="100">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
</universes>
</lattice>
<!-- Fuel Assembly (Upper Half) -->
<lattice id="101">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
</universes>
</lattice>
<!-- Core Lattice (Lower Half) -->
<lattice id="200">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
</universes>
</lattice>
<!-- Core Lattice (Upper Half) -->
<lattice id="201">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
</universes>
</lattice>
</geometry>

1
tests/test_filter_cell/inputs_true.dat Normal file
View file

@ -0,0 +1 @@
caae173f01f7073d634a68a5c4ce97177423e13596a863976e1c40303dc8c05afed457d5a1aa0ae73627ec953143e4f9c1f45bdbd3b0cca76433062467d59777

272
tests/test_filter_cell/materials.xml
View file

@ -1,272 +0,0 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />
<nuclide name="U-234" ao="4.9476e-6" />
<nuclide name="U-235" ao="4.8218e-4" />
<nuclide name="U-236" ao="9.0402e-5" />
<nuclide name="U-238" ao="2.1504e-2" />
<nuclide name="Np-237" ao="7.3733e-6" />
<nuclide name="Pu-238" ao="1.5148e-6" />
<nuclide name="Pu-239" ao="1.3955e-4" />
<nuclide name="Pu-240" ao="3.4405e-5" />
<nuclide name="Pu-241" ao="2.1439e-5" />
<nuclide name="Pu-242" ao="3.7422e-6" />
<nuclide name="Am-241" ao="4.5041e-7" />
<nuclide name="Am-242m" ao="9.2301e-9" />
<nuclide name="Am-243" ao="4.7878e-7" />
<nuclide name="Cm-242" ao="1.0485e-7" />
<nuclide name="Cm-243" ao="1.4268e-9" />
<nuclide name="Cm-244" ao="8.8756e-8" />
<nuclide name="Cm-245" ao="3.5285e-9" />
<nuclide name="Mo-95" ao="2.6497e-5" />
<nuclide name="Tc-99" ao="3.2772e-5" />
<nuclide name="Ru-101" ao="3.0742e-5" />
<nuclide name="Ru-103" ao="2.3505e-6" />
<nuclide name="Ag-109" ao="2.0009e-6" />
<nuclide name="Xe-135" ao="1.0801e-8" />
<nuclide name="Cs-133" ao="3.4612e-5" />
<nuclide name="Nd-143" ao="2.6078e-5" />
<nuclide name="Nd-145" ao="1.9898e-5" />
<nuclide name="Sm-147" ao="1.6128e-6" />
<nuclide name="Sm-149" ao="1.1627e-7" />
<nuclide name="Sm-150" ao="7.1727e-6" />
<nuclide name="Sm-151" ao="5.4947e-7" />
<nuclide name="Sm-152" ao="3.0221e-6" />
<nuclide name="Eu-153" ao="2.6209e-6" />
<nuclide name="Gd-155" ao="1.5369e-9" />
<nuclide name="O-16" ao="4.5737e-2" />
</material>
<!-- Cladding composition -->
<material id="2">
<density value="5.77" units="g/cm3" />
<nuclide name="Zr-90" ao="0.5145" />
<nuclide name="Zr-91" ao="0.1122" />
<nuclide name="Zr-92" ao="0.1715" />
<nuclide name="Zr-94" ao="0.1738" />
<nuclide name="Zr-96" ao="0.0280" />
</material>
<!-- Cold borated water -->
<material id="3">
<density value="0.07416" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Hot borated water -->
<material id="4">
<density value="0.06614" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
<sab name="HH2O" xs="71t" />
</material>
<!-- RPV Composition -->
<material id="5">
<density value="7.9" units="g/cm3" />
<nuclide name="Fe-54" wo="0.05437098" />
<nuclide name="Fe-56" wo="0.88500663" />
<nuclide name="Fe-57" wo="0.0208008" />
<nuclide name="Fe-58" wo="0.00282159" />
<nuclide name="Ni-58" wo="0.0067198" />
<nuclide name="Ni-60" wo="0.0026776" />
<nuclide name="Ni-61" wo="0.0001183" />
<nuclide name="Ni-62" wo="0.0003835" />
<nuclide name="Ni-64" wo="0.0001008" />
<nuclide name="Mn-55" wo="0.01" />
<nuclide name="Mo-92" wo="0.000849" />
<nuclide name="Mo-94" wo="0.0005418" />
<nuclide name="Mo-95" wo="0.0009438" />
<nuclide name="Mo-96" wo="0.0010002" />
<nuclide name="Mo-97" wo="0.0005796" />
<nuclide name="Mo-98" wo="0.0014814" />
<nuclide name="Mo-100" wo="0.0006042" />
<nuclide name="Si-28" wo="0.00367464" />
<nuclide name="Si-29" wo="0.00019336" />
<nuclide name="Si-30" wo="0.000132" />
<nuclide name="Cr-50" wo="0.00010435" />
<nuclide name="Cr-52" wo="0.002092475" />
<nuclide name="Cr-53" wo="0.00024185" />
<nuclide name="Cr-54" wo="6.1325e-05" />
<nuclide name="C-Nat" wo="0.0025" />
<nuclide name="Cu-63" wo="0.0013696" />
<nuclide name="Cu-65" wo="0.0006304" />
</material>
<!-- Lower radial reflector -->
<material id="6">
<density value="4.32" units="g/cm3" />
<nuclide name="H-1" wo="0.0095661" />
<nuclide name="O-16" wo="0.0759107" />
<nuclide name="B-10" wo="3.08409e-5" />
<nuclide name="B-11" wo="1.40499e-4" />
<nuclide name="Fe-54" wo="0.035620772088" />
<nuclide name="Fe-56" wo="0.579805982228" />
<nuclide name="Fe-57" wo="0.01362750048" />
<nuclide name="Fe-58" wo="0.001848545204" />
<nuclide name="Ni-58" wo="0.055298376566" />
<nuclide name="Ni-60" wo="0.022034425592" />
<nuclide name="Ni-61" wo="0.000973510811" />
<nuclide name="Ni-62" wo="0.003155886695" />
<nuclide name="Ni-64" wo="0.000829500336" />
<nuclide name="Mn-55" wo="0.0182870" />
<nuclide name="Si-28" wo="0.00839976771" />
<nuclide name="Si-29" wo="0.00044199679" />
<nuclide name="Si-30" wo="0.0003017355" />
<nuclide name="Cr-50" wo="0.007251360806" />
<nuclide name="Cr-52" wo="0.145407678031" />
<nuclide name="Cr-53" wo="0.016806340306" />
<nuclide name="Cr-54" wo="0.004261520857" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Upper radial reflector / Top plate region -->
<material id="7">
<density value="4.28" units="g/cm3" />
<nuclide name="H-1" wo="0.0086117" />
<nuclide name="O-16" wo="0.0683369" />
<nuclide name="B-10" wo="2.77638e-5" />
<nuclide name="B-11" wo="1.26481e-4" />
<nuclide name="Fe-54" wo="0.035953677186" />
<nuclide name="Fe-56" wo="0.585224740891" />
<nuclide name="Fe-57" wo="0.01375486056" />
<nuclide name="Fe-58" wo="0.001865821363" />
<nuclide name="Ni-58" wo="0.055815129186" />
<nuclide name="Ni-60" wo="0.022240333032" />
<nuclide name="Ni-61" wo="0.000982608081" />
<nuclide name="Ni-62" wo="0.003185377845" />
<nuclide name="Ni-64" wo="0.000837251856" />
<nuclide name="Mn-55" wo="0.0184579" />
<nuclide name="Si-28" wo="0.00847831314" />
<nuclide name="Si-29" wo="0.00044612986" />
<nuclide name="Si-30" wo="0.000304557" />
<nuclide name="Cr-50" wo="0.00731912987" />
<nuclide name="Cr-52" wo="0.146766614995" />
<nuclide name="Cr-53" wo="0.01696340737" />
<nuclide name="Cr-54" wo="0.004301347765" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom plate region -->
<material id="8">
<density value="7.184" units="g/cm3" />
<nuclide name="H-1" wo="0.0011505" />
<nuclide name="O-16" wo="0.0091296" />
<nuclide name="B-10" wo="3.70915e-6" />
<nuclide name="B-11" wo="1.68974e-5" />
<nuclide name="Fe-54" wo="0.03855611055" />
<nuclide name="Fe-56" wo="0.627585036425" />
<nuclide name="Fe-57" wo="0.014750478" />
<nuclide name="Fe-58" wo="0.002000875025" />
<nuclide name="Ni-58" wo="0.059855207342" />
<nuclide name="Ni-60" wo="0.023850159704" />
<nuclide name="Ni-61" wo="0.001053732407" />
<nuclide name="Ni-62" wo="0.003415945715" />
<nuclide name="Ni-64" wo="0.000897854832" />
<nuclide name="Mn-55" wo="0.0197940" />
<nuclide name="Si-28" wo="0.00909197802" />
<nuclide name="Si-29" wo="0.00047842098" />
<nuclide name="Si-30" wo="0.000326601" />
<nuclide name="Cr-50" wo="0.007848910646" />
<nuclide name="Cr-52" wo="0.157390026871" />
<nuclide name="Cr-53" wo="0.018191270146" />
<nuclide name="Cr-54" wo="0.004612692337" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom nozzle region -->
<material id="9">
<density value="2.53" units="g/cm3" />
<nuclide name="H-1" wo="0.0245014" />
<nuclide name="O-16" wo="0.1944274" />
<nuclide name="B-10" wo="7.89917e-5" />
<nuclide name="B-11" wo="3.59854e-4" />
<nuclide name="Fe-54" wo="0.030411411144" />
<nuclide name="Fe-56" wo="0.495012237964" />
<nuclide name="Fe-57" wo="0.01163454624" />
<nuclide name="Fe-58" wo="0.001578204652" />
<nuclide name="Ni-58" wo="0.047211231662" />
<nuclide name="Ni-60" wo="0.018811987544" />
<nuclide name="Ni-61" wo="0.000831139127" />
<nuclide name="Ni-62" wo="0.002694352115" />
<nuclide name="Ni-64" wo="0.000708189552" />
<nuclide name="Mn-55" wo="0.0156126" />
<nuclide name="Si-28" wo="0.007171335558" />
<nuclide name="Si-29" wo="0.000377356542" />
<nuclide name="Si-30" wo="0.0002576079" />
<nuclide name="Cr-50" wo="0.006190885148" />
<nuclide name="Cr-52" wo="0.124142524198" />
<nuclide name="Cr-53" wo="0.014348496148" />
<nuclide name="Cr-54" wo="0.003638294506" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Top nozzle region -->
<material id="10">
<density value="1.746" units="g/cm3" />
<nuclide name="H-1" wo="0.0358870" />
<nuclide name="O-16" wo="0.2847761" />
<nuclide name="B-10" wo="1.15699e-4" />
<nuclide name="B-11" wo="5.27075e-4" />
<nuclide name="Fe-54" wo="0.02644016154" />
<nuclide name="Fe-56" wo="0.43037146399" />
<nuclide name="Fe-57" wo="0.0101152584" />
<nuclide name="Fe-58" wo="0.00137211607" />
<nuclide name="Ni-58" wo="0.04104621835" />
<nuclide name="Ni-60" wo="0.0163554502" />
<nuclide name="Ni-61" wo="0.000722605975" />
<nuclide name="Ni-62" wo="0.002342513875" />
<nuclide name="Ni-64" wo="0.0006157116" />
<nuclide name="Mn-55" wo="0.0135739" />
<nuclide name="Si-28" wo="0.006234853554" />
<nuclide name="Si-29" wo="0.000328078746" />
<nuclide name="Si-30" wo="0.0002239677" />
<nuclide name="Cr-50" wo="0.005382452306" />
<nuclide name="Cr-52" wo="0.107931450781" />
<nuclide name="Cr-53" wo="0.012474806806" />
<nuclide name="Cr-54" wo="0.003163190107" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Top of Fuel Assemblies -->
<material id="11">
<density value="3.044" units="g/cm3" />
<nuclide name="H-1" wo="0.0162913" />
<nuclide name="O-16" wo="0.1292776" />
<nuclide name="B-10" wo="5.25228e-5" />
<nuclide name="B-11" wo="2.39272e-4" />
<nuclide name="Zr-90" wo="0.43313403903" />
<nuclide name="Zr-91" wo="0.09549277374" />
<nuclide name="Zr-92" wo="0.14759527104" />
<nuclide name="Zr-94" wo="0.15280552077" />
<nuclide name="Zr-96" wo="0.02511169542" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom of Fuel Assemblies -->
<material id="12">
<density value="1.762" units="g/cm3" />
<nuclide name="H-1" wo="0.0292856" />
<nuclide name="O-16" wo="0.2323919" />
<nuclide name="B-10" wo="9.44159e-5" />
<nuclide name="B-11" wo="4.30120e-4" />
<nuclide name="Zr-90" wo="0.3741373658" />
<nuclide name="Zr-91" wo="0.0824858164" />
<nuclide name="Zr-92" wo="0.1274914944" />
<nuclide name="Zr-94" wo="0.1319920622" />
<nuclide name="Zr-96" wo="0.0216912612" />
<sab name="HH2O" xs="71t" />
</material>
</materials>

14
tests/test_filter_cell/results_true.dat
View file

@ -1,11 +1,11 @@
k-combined:
1.005983E+00 2.248579E-02
9.903196E-01 4.279617E-02
tally 1:
0.000000E+00
0.000000E+00
1.423676E+01
4.330937E+01
2.914798E+00
1.831649E+00
4.088282E+01
3.662539E+02
1.767552E+01
6.295417E+01
3.863588E+00
3.013300E+00
5.356594E+01
5.839391E+02

19
tests/test_filter_cell/settings.xml
View file

@ -1,19 +0,0 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>100</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>
-160 -160 -183
160 160 183
</parameters>
</space>
</source>
</settings>

25
tests/test_filter_cell/test_filter_cell.py
View file

@ -1,10 +1,29 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
from testing_harness import TestHarness
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterCellTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='cell', bins=(10, 21, 22, 23))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterCellTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterCellTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*', True)
harness = FilterCellTestHarness('statepoint.10.*', True)
harness.main()

181
tests/test_filter_cellborn/geometry.xml
View file

@ -1,181 +0,0 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="z-cylinder" coeffs="0. 0. 0.41" />
<surface id="2" type="z-cylinder" coeffs="0. 0. 0.475" />
<surface id="3" type="z-cylinder" coeffs="0. 0. 0.56" />
<surface id="4" type="z-cylinder" coeffs="0. 0. 0.62" />
<surface id="5" type="z-cylinder" coeffs="0. 0. 187.6" />
<surface id="6" type="z-cylinder" coeffs="0. 0. 209.0" />
<surface id="7" type="z-cylinder" coeffs="0. 0. 229.0" />
<surface id="8" type="z-cylinder" coeffs="0. 0. 249.0" boundary="vacuum" />
<surface id="31" type="z-plane" coeffs="-229.0" boundary="vacuum" />
<surface id="32" type="z-plane" coeffs="-199.0" />
<surface id="33" type="z-plane" coeffs="-193.0" />
<surface id="34" type="z-plane" coeffs="-183.0" />
<surface id="35" type="z-plane" coeffs="0.0" />
<surface id="36" type="z-plane" coeffs="183.0" />
<surface id="37" type="z-plane" coeffs="203.0" />
<surface id="38" type="z-plane" coeffs="215.0" />
<surface id="39" type="z-plane" coeffs="223.0" boundary="vacuum" />
<!-- All geometry on base universe -->
<cell id="1" fill="200" surfaces=" -6 34 -35" /> <!-- Lower core -->
<cell id="2" fill="201" surfaces=" -6 35 -36" /> <!-- Upper core -->
<cell id="3" material="8" surfaces=" -7 31 -32" /> <!-- Lower core plate region -->
<cell id="4" material="9" surfaces=" -5 32 -33" /> <!-- Bottom nozzle region -->
<cell id="5" material="12" surfaces=" -5 33 -34" /> <!-- Bottom FA region -->
<cell id="6" material="11" surfaces=" -5 36 -37" /> <!-- Top FA region -->
<cell id="7" material="10" surfaces=" -5 37 -38" /> <!-- Top nozzle region -->
<cell id="8" material="7" surfaces=" -7 38 -39" /> <!-- Upper plate region -->
<cell id="9" material="4" surfaces="6 -7 32 -38" /> <!-- Downcomer -->
<cell id="10" material="5" surfaces="7 -8 31 -39" /> <!-- RPV -->
<cell id="11" material="6" surfaces="5 -6 32 -34" /> <!-- Bottom of radial reflector -->
<cell id="12" material="7" surfaces="5 -6 36 -38" /> <!-- Top of radial reflector -->
<!-- Fuel pin, cladding, cold water -->
<cell id="21" universe="1" material="1" surfaces="-1" />
<cell id="22" universe="1" material="2" surfaces="1 -2" />
<cell id="23" universe="1" material="3" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="24" universe="2" material="3" surfaces="-3" />
<cell id="25" universe="2" material="2" surfaces="3 -4" />
<cell id="26" universe="2" material="3" surfaces="4" />
<!-- Fuel pin, cladding, hot water -->
<cell id="27" universe="3" material="1" surfaces="-1" />
<cell id="28" universe="3" material="2" surfaces="1 -2" />
<cell id="29" universe="3" material="4" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="30" universe="4" material="4" surfaces="-3" />
<cell id="31" universe="4" material="2" surfaces="3 -4" />
<cell id="32" universe="4" material="4" surfaces="4" />
<!-- cell for water assembly (cold) -->
<cell id="50" universe="5" material="4" surfaces="34 -35" />
<!-- containing cell for fuel assembly -->
<cell id="60" universe="6" fill="100" surfaces="34 -35" />
<!-- cell for water assembly (hot) -->
<cell id="70" universe="7" material="3" surfaces="35 -36" />
<!-- containing cell for fuel assembly -->
<cell id="80" universe="8" fill="101" surfaces="35 -36" />
<!-- Fuel Assembly (Lower Half) -->
<lattice id="100">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
</universes>
</lattice>
<!-- Fuel Assembly (Upper Half) -->
<lattice id="101">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
</universes>
</lattice>
<!-- Core Lattice (Lower Half) -->
<lattice id="200">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
</universes>
</lattice>
<!-- Core Lattice (Upper Half) -->
<lattice id="201">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
</universes>
</lattice>
</geometry>

1
tests/test_filter_cellborn/inputs_true.dat Normal file
View file

@ -0,0 +1 @@
2f24eb86cda981982a8db5bb110c72e9cef542c06e15b748c1c7e459f96b0d8ba0978b51dffc006e813cd2e2ae1fa0357336f8322ae263189841afde01f0327b

272
tests/test_filter_cellborn/materials.xml
View file

@ -1,272 +0,0 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />
<nuclide name="U-234" ao="4.9476e-6" />
<nuclide name="U-235" ao="4.8218e-4" />
<nuclide name="U-236" ao="9.0402e-5" />
<nuclide name="U-238" ao="2.1504e-2" />
<nuclide name="Np-237" ao="7.3733e-6" />
<nuclide name="Pu-238" ao="1.5148e-6" />
<nuclide name="Pu-239" ao="1.3955e-4" />
<nuclide name="Pu-240" ao="3.4405e-5" />
<nuclide name="Pu-241" ao="2.1439e-5" />
<nuclide name="Pu-242" ao="3.7422e-6" />
<nuclide name="Am-241" ao="4.5041e-7" />
<nuclide name="Am-242m" ao="9.2301e-9" />
<nuclide name="Am-243" ao="4.7878e-7" />
<nuclide name="Cm-242" ao="1.0485e-7" />
<nuclide name="Cm-243" ao="1.4268e-9" />
<nuclide name="Cm-244" ao="8.8756e-8" />
<nuclide name="Cm-245" ao="3.5285e-9" />
<nuclide name="Mo-95" ao="2.6497e-5" />
<nuclide name="Tc-99" ao="3.2772e-5" />
<nuclide name="Ru-101" ao="3.0742e-5" />
<nuclide name="Ru-103" ao="2.3505e-6" />
<nuclide name="Ag-109" ao="2.0009e-6" />
<nuclide name="Xe-135" ao="1.0801e-8" />
<nuclide name="Cs-133" ao="3.4612e-5" />
<nuclide name="Nd-143" ao="2.6078e-5" />
<nuclide name="Nd-145" ao="1.9898e-5" />
<nuclide name="Sm-147" ao="1.6128e-6" />
<nuclide name="Sm-149" ao="1.1627e-7" />
<nuclide name="Sm-150" ao="7.1727e-6" />
<nuclide name="Sm-151" ao="5.4947e-7" />
<nuclide name="Sm-152" ao="3.0221e-6" />
<nuclide name="Eu-153" ao="2.6209e-6" />
<nuclide name="Gd-155" ao="1.5369e-9" />
<nuclide name="O-16" ao="4.5737e-2" />
</material>
<!-- Cladding composition -->
<material id="2">
<density value="5.77" units="g/cm3" />
<nuclide name="Zr-90" ao="0.5145" />
<nuclide name="Zr-91" ao="0.1122" />
<nuclide name="Zr-92" ao="0.1715" />
<nuclide name="Zr-94" ao="0.1738" />
<nuclide name="Zr-96" ao="0.0280" />
</material>
<!-- Cold borated water -->
<material id="3">
<density value="0.07416" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Hot borated water -->
<material id="4">
<density value="0.06614" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
<sab name="HH2O" xs="71t" />
</material>
<!-- RPV Composition -->
<material id="5">
<density value="7.9" units="g/cm3" />
<nuclide name="Fe-54" wo="0.05437098" />
<nuclide name="Fe-56" wo="0.88500663" />
<nuclide name="Fe-57" wo="0.0208008" />
<nuclide name="Fe-58" wo="0.00282159" />
<nuclide name="Ni-58" wo="0.0067198" />
<nuclide name="Ni-60" wo="0.0026776" />
<nuclide name="Ni-61" wo="0.0001183" />
<nuclide name="Ni-62" wo="0.0003835" />
<nuclide name="Ni-64" wo="0.0001008" />
<nuclide name="Mn-55" wo="0.01" />
<nuclide name="Mo-92" wo="0.000849" />
<nuclide name="Mo-94" wo="0.0005418" />
<nuclide name="Mo-95" wo="0.0009438" />
<nuclide name="Mo-96" wo="0.0010002" />
<nuclide name="Mo-97" wo="0.0005796" />
<nuclide name="Mo-98" wo="0.0014814" />
<nuclide name="Mo-100" wo="0.0006042" />
<nuclide name="Si-28" wo="0.00367464" />
<nuclide name="Si-29" wo="0.00019336" />
<nuclide name="Si-30" wo="0.000132" />
<nuclide name="Cr-50" wo="0.00010435" />
<nuclide name="Cr-52" wo="0.002092475" />
<nuclide name="Cr-53" wo="0.00024185" />
<nuclide name="Cr-54" wo="6.1325e-05" />
<nuclide name="C-Nat" wo="0.0025" />
<nuclide name="Cu-63" wo="0.0013696" />
<nuclide name="Cu-65" wo="0.0006304" />
</material>
<!-- Lower radial reflector -->
<material id="6">
<density value="4.32" units="g/cm3" />
<nuclide name="H-1" wo="0.0095661" />
<nuclide name="O-16" wo="0.0759107" />
<nuclide name="B-10" wo="3.08409e-5" />
<nuclide name="B-11" wo="1.40499e-4" />
<nuclide name="Fe-54" wo="0.035620772088" />
<nuclide name="Fe-56" wo="0.579805982228" />
<nuclide name="Fe-57" wo="0.01362750048" />
<nuclide name="Fe-58" wo="0.001848545204" />
<nuclide name="Ni-58" wo="0.055298376566" />
<nuclide name="Ni-60" wo="0.022034425592" />
<nuclide name="Ni-61" wo="0.000973510811" />
<nuclide name="Ni-62" wo="0.003155886695" />
<nuclide name="Ni-64" wo="0.000829500336" />
<nuclide name="Mn-55" wo="0.0182870" />
<nuclide name="Si-28" wo="0.00839976771" />
<nuclide name="Si-29" wo="0.00044199679" />
<nuclide name="Si-30" wo="0.0003017355" />
<nuclide name="Cr-50" wo="0.007251360806" />
<nuclide name="Cr-52" wo="0.145407678031" />
<nuclide name="Cr-53" wo="0.016806340306" />
<nuclide name="Cr-54" wo="0.004261520857" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Upper radial reflector / Top plate region -->
<material id="7">
<density value="4.28" units="g/cm3" />
<nuclide name="H-1" wo="0.0086117" />
<nuclide name="O-16" wo="0.0683369" />
<nuclide name="B-10" wo="2.77638e-5" />
<nuclide name="B-11" wo="1.26481e-4" />
<nuclide name="Fe-54" wo="0.035953677186" />
<nuclide name="Fe-56" wo="0.585224740891" />
<nuclide name="Fe-57" wo="0.01375486056" />
<nuclide name="Fe-58" wo="0.001865821363" />
<nuclide name="Ni-58" wo="0.055815129186" />
<nuclide name="Ni-60" wo="0.022240333032" />
<nuclide name="Ni-61" wo="0.000982608081" />
<nuclide name="Ni-62" wo="0.003185377845" />
<nuclide name="Ni-64" wo="0.000837251856" />
<nuclide name="Mn-55" wo="0.0184579" />
<nuclide name="Si-28" wo="0.00847831314" />
<nuclide name="Si-29" wo="0.00044612986" />
<nuclide name="Si-30" wo="0.000304557" />
<nuclide name="Cr-50" wo="0.00731912987" />
<nuclide name="Cr-52" wo="0.146766614995" />
<nuclide name="Cr-53" wo="0.01696340737" />
<nuclide name="Cr-54" wo="0.004301347765" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom plate region -->
<material id="8">
<density value="7.184" units="g/cm3" />
<nuclide name="H-1" wo="0.0011505" />
<nuclide name="O-16" wo="0.0091296" />
<nuclide name="B-10" wo="3.70915e-6" />
<nuclide name="B-11" wo="1.68974e-5" />
<nuclide name="Fe-54" wo="0.03855611055" />
<nuclide name="Fe-56" wo="0.627585036425" />
<nuclide name="Fe-57" wo="0.014750478" />
<nuclide name="Fe-58" wo="0.002000875025" />
<nuclide name="Ni-58" wo="0.059855207342" />
<nuclide name="Ni-60" wo="0.023850159704" />
<nuclide name="Ni-61" wo="0.001053732407" />
<nuclide name="Ni-62" wo="0.003415945715" />
<nuclide name="Ni-64" wo="0.000897854832" />
<nuclide name="Mn-55" wo="0.0197940" />
<nuclide name="Si-28" wo="0.00909197802" />
<nuclide name="Si-29" wo="0.00047842098" />
<nuclide name="Si-30" wo="0.000326601" />
<nuclide name="Cr-50" wo="0.007848910646" />
<nuclide name="Cr-52" wo="0.157390026871" />
<nuclide name="Cr-53" wo="0.018191270146" />
<nuclide name="Cr-54" wo="0.004612692337" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom nozzle region -->
<material id="9">
<density value="2.53" units="g/cm3" />
<nuclide name="H-1" wo="0.0245014" />
<nuclide name="O-16" wo="0.1944274" />
<nuclide name="B-10" wo="7.89917e-5" />
<nuclide name="B-11" wo="3.59854e-4" />
<nuclide name="Fe-54" wo="0.030411411144" />
<nuclide name="Fe-56" wo="0.495012237964" />
<nuclide name="Fe-57" wo="0.01163454624" />
<nuclide name="Fe-58" wo="0.001578204652" />
<nuclide name="Ni-58" wo="0.047211231662" />
<nuclide name="Ni-60" wo="0.018811987544" />
<nuclide name="Ni-61" wo="0.000831139127" />
<nuclide name="Ni-62" wo="0.002694352115" />
<nuclide name="Ni-64" wo="0.000708189552" />
<nuclide name="Mn-55" wo="0.0156126" />
<nuclide name="Si-28" wo="0.007171335558" />
<nuclide name="Si-29" wo="0.000377356542" />
<nuclide name="Si-30" wo="0.0002576079" />
<nuclide name="Cr-50" wo="0.006190885148" />
<nuclide name="Cr-52" wo="0.124142524198" />
<nuclide name="Cr-53" wo="0.014348496148" />
<nuclide name="Cr-54" wo="0.003638294506" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Top nozzle region -->
<material id="10">
<density value="1.746" units="g/cm3" />
<nuclide name="H-1" wo="0.0358870" />
<nuclide name="O-16" wo="0.2847761" />
<nuclide name="B-10" wo="1.15699e-4" />
<nuclide name="B-11" wo="5.27075e-4" />
<nuclide name="Fe-54" wo="0.02644016154" />
<nuclide name="Fe-56" wo="0.43037146399" />
<nuclide name="Fe-57" wo="0.0101152584" />
<nuclide name="Fe-58" wo="0.00137211607" />
<nuclide name="Ni-58" wo="0.04104621835" />
<nuclide name="Ni-60" wo="0.0163554502" />
<nuclide name="Ni-61" wo="0.000722605975" />
<nuclide name="Ni-62" wo="0.002342513875" />
<nuclide name="Ni-64" wo="0.0006157116" />
<nuclide name="Mn-55" wo="0.0135739" />
<nuclide name="Si-28" wo="0.006234853554" />
<nuclide name="Si-29" wo="0.000328078746" />
<nuclide name="Si-30" wo="0.0002239677" />
<nuclide name="Cr-50" wo="0.005382452306" />
<nuclide name="Cr-52" wo="0.107931450781" />
<nuclide name="Cr-53" wo="0.012474806806" />
<nuclide name="Cr-54" wo="0.003163190107" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Top of Fuel Assemblies -->
<material id="11">
<density value="3.044" units="g/cm3" />
<nuclide name="H-1" wo="0.0162913" />
<nuclide name="O-16" wo="0.1292776" />
<nuclide name="B-10" wo="5.25228e-5" />
<nuclide name="B-11" wo="2.39272e-4" />
<nuclide name="Zr-90" wo="0.43313403903" />
<nuclide name="Zr-91" wo="0.09549277374" />
<nuclide name="Zr-92" wo="0.14759527104" />
<nuclide name="Zr-94" wo="0.15280552077" />
<nuclide name="Zr-96" wo="0.02511169542" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom of Fuel Assemblies -->
<material id="12">
<density value="1.762" units="g/cm3" />
<nuclide name="H-1" wo="0.0292856" />
<nuclide name="O-16" wo="0.2323919" />
<nuclide name="B-10" wo="9.44159e-5" />
<nuclide name="B-11" wo="4.30120e-4" />
<nuclide name="Zr-90" wo="0.3741373658" />
<nuclide name="Zr-91" wo="0.0824858164" />
<nuclide name="Zr-92" wo="0.1274914944" />
<nuclide name="Zr-94" wo="0.1319920622" />
<nuclide name="Zr-96" wo="0.0216912612" />
<sab name="HH2O" xs="71t" />
</material>
</materials>

6
tests/test_filter_cellborn/results_true.dat
View file

@ -1,10 +1,10 @@
k-combined:
1.005983E+00 2.248579E-02
9.903196E-01 4.279617E-02
tally 1:
0.000000E+00
0.000000E+00
7.007584E+01
1.050827E+03
8.921179E+01
1.601939E+03
0.000000E+00
0.000000E+00
0.000000E+00

19
tests/test_filter_cellborn/settings.xml
View file

@ -1,19 +0,0 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>100</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>
-160 -160 -183
160 160 183
</parameters>
</space>
</source>
</settings>

9
tests/test_filter_cellborn/tallies.xml
View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<filter type="cellborn" bins="10 21 22 23" />
<scores>total</scores>
</tally>
</tallies>

25
tests/test_filter_cellborn/test_filter_cellborn.py
View file

@ -1,10 +1,29 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
from testing_harness import TestHarness
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterCellbornTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='cellborn', bins=(10, 21, 22, 23))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterCellbornTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterCellbornTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*', True)
harness = FilterCellbornTestHarness('statepoint.10.*', True)
harness.main()

6
tests/test_filter_distribcell/case-1/geometry.xml
View file

@ -1,9 +1,9 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="5" surfaces="1 -2 3 -4" />
<cell id="201" universe="2" material="1" surfaces="-5" />
<cell id="202" universe="2" material="2" surfaces="5" />
<cell id="1" fill="5" region="1 -2 3 -4" />
<cell id="201" universe="2" material="1" region="-5" />
<cell id="202" universe="2" material="2" region="5" />
<cell id="301" universe="3" material="2"/>
<lattice id="5">

18
tests/test_filter_distribcell/case-2/geometry.xml
View file

@ -1,15 +1,15 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="5" surfaces="1 -2 3 -4" />
<cell id="201" universe="21" material="1" surfaces="-5" />
<cell id="202" universe="21" material="2" surfaces="5" />
<cell id="203" universe="22" material="1" surfaces="-5" />
<cell id="204" universe="22" material="2" surfaces="5" />
<cell id="205" universe="23" material="1" surfaces="-5" />
<cell id="206" universe="23" material="2" surfaces="5" />
<cell id="207" universe="20" material="1" surfaces="-5" />
<cell id="208" universe="20" material="2" surfaces="5" />
<cell id="1" fill="5" region="1 -2 3 -4" />
<cell id="201" universe="21" material="1" region="-5" />
<cell id="202" universe="21" material="2" region="5" />
<cell id="203" universe="22" material="1" region="-5" />
<cell id="204" universe="22" material="2" region="5" />
<cell id="205" universe="23" material="1" region="-5" />
<cell id="206" universe="23" material="2" region="5" />
<cell id="207" universe="20" material="1" region="-5" />
<cell id="208" universe="20" material="2" region="5" />
<lattice id="5">
<dimension>2 2</dimension>

56
tests/test_filter_distribcell/case-3/geometry.xml
View file

@ -21,50 +21,50 @@
<surface id="39" type="z-plane" coeffs="223.0" boundary="vacuum" />
<!-- All geometry on base universe -->
<cell id="1" fill="200" surfaces=" -6 34 -35" /> <!-- Lower core -->
<cell id="2" fill="201" surfaces=" -6 35 -36" /> <!-- Upper core -->
<cell id="3" material="8" surfaces=" -7 31 -32" /> <!-- Lower core plate region -->
<cell id="4" material="9" surfaces=" -5 32 -33" /> <!-- Bottom nozzle region -->
<cell id="5" material="12" surfaces=" -5 33 -34" /> <!-- Bottom FA region -->
<cell id="6" material="11" surfaces=" -5 36 -37" /> <!-- Top FA region -->
<cell id="7" material="10" surfaces=" -5 37 -38" /> <!-- Top nozzle region -->
<cell id="8" material="7" surfaces=" -7 38 -39" /> <!-- Upper plate region -->
<cell id="9" material="4" surfaces="6 -7 32 -38" /> <!-- Downcomer -->
<cell id="10" material="5" surfaces="7 -8 31 -39" /> <!-- RPV -->
<cell id="11" material="6" surfaces="5 -6 32 -34" /> <!-- Bottom of radial reflector -->
<cell id="12" material="7" surfaces="5 -6 36 -38" /> <!-- Top of radial reflector -->
<cell id="1" fill="200" region=" -6 34 -35" /> <!-- Lower core -->
<cell id="2" fill="201" region=" -6 35 -36" /> <!-- Upper core -->
<cell id="3" material="8" region=" -7 31 -32" /> <!-- Lower core plate region -->
<cell id="4" material="9" region=" -5 32 -33" /> <!-- Bottom nozzle region -->
<cell id="5" material="12" region=" -5 33 -34" /> <!-- Bottom FA region -->
<cell id="6" material="11" region=" -5 36 -37" /> <!-- Top FA region -->
<cell id="7" material="10" region=" -5 37 -38" /> <!-- Top nozzle region -->
<cell id="8" material="7" region=" -7 38 -39" /> <!-- Upper plate region -->
<cell id="9" material="4" region="6 -7 32 -38" /> <!-- Downcomer -->
<cell id="10" material="5" region="7 -8 31 -39" /> <!-- RPV -->
<cell id="11" material="6" region="5 -6 32 -34" /> <!-- Bottom of radial reflector -->
<cell id="12" material="7" region="5 -6 36 -38" /> <!-- Top of radial reflector -->
<!-- Fuel pin, cladding, cold water -->
<cell id="21" universe="1" material="1" surfaces="-1" />
<cell id="22" universe="1" material="2" surfaces="1 -2" />
<cell id="23" universe="1" material="3" surfaces="2" />
<cell id="21" universe="1" material="1" region="-1" />
<cell id="22" universe="1" material="2" region="1 -2" />
<cell id="23" universe="1" material="3" region="2" />
<!-- Instrumentation guide tube -->
<cell id="24" universe="2" material="3" surfaces="-3" />
<cell id="25" universe="2" material="2" surfaces="3 -4" />
<cell id="26" universe="2" material="3" surfaces="4" />
<cell id="24" universe="2" material="3" region="-3" />
<cell id="25" universe="2" material="2" region="3 -4" />
<cell id="26" universe="2" material="3" region="4" />
<!-- Fuel pin, cladding, hot water -->
<cell id="27" universe="3" material="1" surfaces="-1" />
<cell id="28" universe="3" material="2" surfaces="1 -2" />
<cell id="29" universe="3" material="4" surfaces="2" />
<cell id="27" universe="3" material="1" region="-1" />
<cell id="28" universe="3" material="2" region="1 -2" />
<cell id="29" universe="3" material="4" region="2" />
<!-- Instrumentation guide tube -->
<cell id="30" universe="4" material="4" surfaces="-3" />
<cell id="31" universe="4" material="2" surfaces="3 -4" />
<cell id="32" universe="4" material="4" surfaces="4" />
<cell id="30" universe="4" material="4" region="-3" />
<cell id="31" universe="4" material="2" region="3 -4" />
<cell id="32" universe="4" material="4" region="4" />
<!-- cell for water assembly (cold) -->
<cell id="50" universe="5" material="4" surfaces="34 -35" />
<cell id="50" universe="5" material="4" region="34 -35" />
<!-- containing cell for fuel assembly -->
<cell id="60" universe="6" fill="100" surfaces="34 -35" />
<cell id="60" universe="6" fill="100" region="34 -35" />
<!-- cell for water assembly (hot) -->
<cell id="70" universe="7" material="3" surfaces="35 -36" />
<cell id="70" universe="7" material="3" region="35 -36" />
<!-- containing cell for fuel assembly -->
<cell id="80" universe="8" fill="101" surfaces="35 -36" />
<cell id="80" universe="8" fill="101" region="35 -36" />
<!-- Fuel Assembly (Lower Half) -->
<lattice id="100">

8
tests/test_filter_distribcell/case-4/geometry.xml
View file

@ -1,9 +1,9 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell fill="5" id="1" name="Cell 1" surfaces="1 -2 3 -4" universe="0" />
<cell id="101" material="1" name="cell 2" surfaces="-5" universe="1" />
<cell id="102" material="2" name="cell 3" surfaces="5" universe="1" />
<cell id="500" material="2" name="cell 4" surfaces="" universe="3" />
<cell fill="5" id="1" name="Cell 1" region="1 -2 3 -4" universe="0" />
<cell id="101" material="1" name="cell 2" region="-5" universe="1" />
<cell id="102" material="2" name="cell 3" region="5" universe="1" />
<cell id="500" material="2" name="cell 4" region="" universe="3" />
<hex_lattice id="5" n_rings="2">
<pitch>1.0</pitch>
<outer>3</outer>

3
tests/test_filter_distribcell/test_filter_distribcell.py
View file

@ -4,8 +4,7 @@ import glob
import hashlib
import os
import sys
sys.path.insert(0, '..')
sys.path.insert(0, os.pardir)
from testing_harness import *

181
tests/test_filter_energy/geometry.xml
View file

@ -1,181 +0,0 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="z-cylinder" coeffs="0. 0. 0.41" />
<surface id="2" type="z-cylinder" coeffs="0. 0. 0.475" />
<surface id="3" type="z-cylinder" coeffs="0. 0. 0.56" />
<surface id="4" type="z-cylinder" coeffs="0. 0. 0.62" />
<surface id="5" type="z-cylinder" coeffs="0. 0. 187.6" />
<surface id="6" type="z-cylinder" coeffs="0. 0. 209.0" />
<surface id="7" type="z-cylinder" coeffs="0. 0. 229.0" />
<surface id="8" type="z-cylinder" coeffs="0. 0. 249.0" boundary="vacuum" />
<surface id="31" type="z-plane" coeffs="-229.0" boundary="vacuum" />
<surface id="32" type="z-plane" coeffs="-199.0" />
<surface id="33" type="z-plane" coeffs="-193.0" />
<surface id="34" type="z-plane" coeffs="-183.0" />
<surface id="35" type="z-plane" coeffs="0.0" />
<surface id="36" type="z-plane" coeffs="183.0" />
<surface id="37" type="z-plane" coeffs="203.0" />
<surface id="38" type="z-plane" coeffs="215.0" />
<surface id="39" type="z-plane" coeffs="223.0" boundary="vacuum" />
<!-- All geometry on base universe -->
<cell id="1" fill="200" surfaces=" -6 34 -35" /> <!-- Lower core -->
<cell id="2" fill="201" surfaces=" -6 35 -36" /> <!-- Upper core -->
<cell id="3" material="8" surfaces=" -7 31 -32" /> <!-- Lower core plate region -->
<cell id="4" material="9" surfaces=" -5 32 -33" /> <!-- Bottom nozzle region -->
<cell id="5" material="12" surfaces=" -5 33 -34" /> <!-- Bottom FA region -->
<cell id="6" material="11" surfaces=" -5 36 -37" /> <!-- Top FA region -->
<cell id="7" material="10" surfaces=" -5 37 -38" /> <!-- Top nozzle region -->
<cell id="8" material="7" surfaces=" -7 38 -39" /> <!-- Upper plate region -->
<cell id="9" material="4" surfaces="6 -7 32 -38" /> <!-- Downcomer -->
<cell id="10" material="5" surfaces="7 -8 31 -39" /> <!-- RPV -->
<cell id="11" material="6" surfaces="5 -6 32 -34" /> <!-- Bottom of radial reflector -->
<cell id="12" material="7" surfaces="5 -6 36 -38" /> <!-- Top of radial reflector -->
<!-- Fuel pin, cladding, cold water -->
<cell id="21" universe="1" material="1" surfaces="-1" />
<cell id="22" universe="1" material="2" surfaces="1 -2" />
<cell id="23" universe="1" material="3" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="24" universe="2" material="3" surfaces="-3" />
<cell id="25" universe="2" material="2" surfaces="3 -4" />
<cell id="26" universe="2" material="3" surfaces="4" />
<!-- Fuel pin, cladding, hot water -->
<cell id="27" universe="3" material="1" surfaces="-1" />
<cell id="28" universe="3" material="2" surfaces="1 -2" />
<cell id="29" universe="3" material="4" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="30" universe="4" material="4" surfaces="-3" />
<cell id="31" universe="4" material="2" surfaces="3 -4" />
<cell id="32" universe="4" material="4" surfaces="4" />
<!-- cell for water assembly (cold) -->
<cell id="50" universe="5" material="4" surfaces="34 -35" />
<!-- containing cell for fuel assembly -->
<cell id="60" universe="6" fill="100" surfaces="34 -35" />
<!-- cell for water assembly (hot) -->
<cell id="70" universe="7" material="3" surfaces="35 -36" />
<!-- containing cell for fuel assembly -->
<cell id="80" universe="8" fill="101" surfaces="35 -36" />
<!-- Fuel Assembly (Lower Half) -->
<lattice id="100">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
</universes>
</lattice>
<!-- Fuel Assembly (Upper Half) -->
<lattice id="101">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
</universes>
</lattice>
<!-- Core Lattice (Lower Half) -->
<lattice id="200">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
</universes>
</lattice>
<!-- Core Lattice (Upper Half) -->
<lattice id="201">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
</universes>
</lattice>
</geometry>

1
tests/test_filter_energy/inputs_true.dat Normal file
View file

@ -0,0 +1 @@
49835200052ee1a4c9583a7bb4e9430de7d01a6d7a8d4f63ae37be9bbac4fd8c7ed9135ecbc4ab4cb075fd4d77b322265783463dd07c127decceee8c9fe25bfd

272
tests/test_filter_energy/materials.xml
View file

@ -1,272 +0,0 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />
<nuclide name="U-234" ao="4.9476e-6" />
<nuclide name="U-235" ao="4.8218e-4" />
<nuclide name="U-236" ao="9.0402e-5" />
<nuclide name="U-238" ao="2.1504e-2" />
<nuclide name="Np-237" ao="7.3733e-6" />
<nuclide name="Pu-238" ao="1.5148e-6" />
<nuclide name="Pu-239" ao="1.3955e-4" />
<nuclide name="Pu-240" ao="3.4405e-5" />
<nuclide name="Pu-241" ao="2.1439e-5" />
<nuclide name="Pu-242" ao="3.7422e-6" />
<nuclide name="Am-241" ao="4.5041e-7" />
<nuclide name="Am-242m" ao="9.2301e-9" />
<nuclide name="Am-243" ao="4.7878e-7" />
<nuclide name="Cm-242" ao="1.0485e-7" />
<nuclide name="Cm-243" ao="1.4268e-9" />
<nuclide name="Cm-244" ao="8.8756e-8" />
<nuclide name="Cm-245" ao="3.5285e-9" />
<nuclide name="Mo-95" ao="2.6497e-5" />
<nuclide name="Tc-99" ao="3.2772e-5" />
<nuclide name="Ru-101" ao="3.0742e-5" />
<nuclide name="Ru-103" ao="2.3505e-6" />
<nuclide name="Ag-109" ao="2.0009e-6" />
<nuclide name="Xe-135" ao="1.0801e-8" />
<nuclide name="Cs-133" ao="3.4612e-5" />
<nuclide name="Nd-143" ao="2.6078e-5" />
<nuclide name="Nd-145" ao="1.9898e-5" />
<nuclide name="Sm-147" ao="1.6128e-6" />
<nuclide name="Sm-149" ao="1.1627e-7" />
<nuclide name="Sm-150" ao="7.1727e-6" />
<nuclide name="Sm-151" ao="5.4947e-7" />
<nuclide name="Sm-152" ao="3.0221e-6" />
<nuclide name="Eu-153" ao="2.6209e-6" />
<nuclide name="Gd-155" ao="1.5369e-9" />
<nuclide name="O-16" ao="4.5737e-2" />
</material>
<!-- Cladding composition -->
<material id="2">
<density value="5.77" units="g/cm3" />
<nuclide name="Zr-90" ao="0.5145" />
<nuclide name="Zr-91" ao="0.1122" />
<nuclide name="Zr-92" ao="0.1715" />
<nuclide name="Zr-94" ao="0.1738" />
<nuclide name="Zr-96" ao="0.0280" />
</material>
<!-- Cold borated water -->
<material id="3">
<density value="0.07416" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Hot borated water -->
<material id="4">
<density value="0.06614" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
<sab name="HH2O" xs="71t" />
</material>
<!-- RPV Composition -->
<material id="5">
<density value="7.9" units="g/cm3" />
<nuclide name="Fe-54" wo="0.05437098" />
<nuclide name="Fe-56" wo="0.88500663" />
<nuclide name="Fe-57" wo="0.0208008" />
<nuclide name="Fe-58" wo="0.00282159" />
<nuclide name="Ni-58" wo="0.0067198" />
<nuclide name="Ni-60" wo="0.0026776" />
<nuclide name="Ni-61" wo="0.0001183" />
<nuclide name="Ni-62" wo="0.0003835" />
<nuclide name="Ni-64" wo="0.0001008" />
<nuclide name="Mn-55" wo="0.01" />
<nuclide name="Mo-92" wo="0.000849" />
<nuclide name="Mo-94" wo="0.0005418" />
<nuclide name="Mo-95" wo="0.0009438" />
<nuclide name="Mo-96" wo="0.0010002" />
<nuclide name="Mo-97" wo="0.0005796" />
<nuclide name="Mo-98" wo="0.0014814" />
<nuclide name="Mo-100" wo="0.0006042" />
<nuclide name="Si-28" wo="0.00367464" />
<nuclide name="Si-29" wo="0.00019336" />
<nuclide name="Si-30" wo="0.000132" />
<nuclide name="Cr-50" wo="0.00010435" />
<nuclide name="Cr-52" wo="0.002092475" />
<nuclide name="Cr-53" wo="0.00024185" />
<nuclide name="Cr-54" wo="6.1325e-05" />
<nuclide name="C-Nat" wo="0.0025" />
<nuclide name="Cu-63" wo="0.0013696" />
<nuclide name="Cu-65" wo="0.0006304" />
</material>
<!-- Lower radial reflector -->
<material id="6">
<density value="4.32" units="g/cm3" />
<nuclide name="H-1" wo="0.0095661" />
<nuclide name="O-16" wo="0.0759107" />
<nuclide name="B-10" wo="3.08409e-5" />
<nuclide name="B-11" wo="1.40499e-4" />
<nuclide name="Fe-54" wo="0.035620772088" />
<nuclide name="Fe-56" wo="0.579805982228" />
<nuclide name="Fe-57" wo="0.01362750048" />
<nuclide name="Fe-58" wo="0.001848545204" />
<nuclide name="Ni-58" wo="0.055298376566" />
<nuclide name="Ni-60" wo="0.022034425592" />
<nuclide name="Ni-61" wo="0.000973510811" />
<nuclide name="Ni-62" wo="0.003155886695" />
<nuclide name="Ni-64" wo="0.000829500336" />
<nuclide name="Mn-55" wo="0.0182870" />
<nuclide name="Si-28" wo="0.00839976771" />
<nuclide name="Si-29" wo="0.00044199679" />
<nuclide name="Si-30" wo="0.0003017355" />
<nuclide name="Cr-50" wo="0.007251360806" />
<nuclide name="Cr-52" wo="0.145407678031" />
<nuclide name="Cr-53" wo="0.016806340306" />
<nuclide name="Cr-54" wo="0.004261520857" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Upper radial reflector / Top plate region -->
<material id="7">
<density value="4.28" units="g/cm3" />
<nuclide name="H-1" wo="0.0086117" />
<nuclide name="O-16" wo="0.0683369" />
<nuclide name="B-10" wo="2.77638e-5" />
<nuclide name="B-11" wo="1.26481e-4" />
<nuclide name="Fe-54" wo="0.035953677186" />
<nuclide name="Fe-56" wo="0.585224740891" />
<nuclide name="Fe-57" wo="0.01375486056" />
<nuclide name="Fe-58" wo="0.001865821363" />
<nuclide name="Ni-58" wo="0.055815129186" />
<nuclide name="Ni-60" wo="0.022240333032" />
<nuclide name="Ni-61" wo="0.000982608081" />
<nuclide name="Ni-62" wo="0.003185377845" />
<nuclide name="Ni-64" wo="0.000837251856" />
<nuclide name="Mn-55" wo="0.0184579" />
<nuclide name="Si-28" wo="0.00847831314" />
<nuclide name="Si-29" wo="0.00044612986" />
<nuclide name="Si-30" wo="0.000304557" />
<nuclide name="Cr-50" wo="0.00731912987" />
<nuclide name="Cr-52" wo="0.146766614995" />
<nuclide name="Cr-53" wo="0.01696340737" />
<nuclide name="Cr-54" wo="0.004301347765" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom plate region -->
<material id="8">
<density value="7.184" units="g/cm3" />
<nuclide name="H-1" wo="0.0011505" />
<nuclide name="O-16" wo="0.0091296" />
<nuclide name="B-10" wo="3.70915e-6" />
<nuclide name="B-11" wo="1.68974e-5" />
<nuclide name="Fe-54" wo="0.03855611055" />
<nuclide name="Fe-56" wo="0.627585036425" />
<nuclide name="Fe-57" wo="0.014750478" />
<nuclide name="Fe-58" wo="0.002000875025" />
<nuclide name="Ni-58" wo="0.059855207342" />
<nuclide name="Ni-60" wo="0.023850159704" />
<nuclide name="Ni-61" wo="0.001053732407" />
<nuclide name="Ni-62" wo="0.003415945715" />
<nuclide name="Ni-64" wo="0.000897854832" />
<nuclide name="Mn-55" wo="0.0197940" />
<nuclide name="Si-28" wo="0.00909197802" />
<nuclide name="Si-29" wo="0.00047842098" />
<nuclide name="Si-30" wo="0.000326601" />
<nuclide name="Cr-50" wo="0.007848910646" />
<nuclide name="Cr-52" wo="0.157390026871" />
<nuclide name="Cr-53" wo="0.018191270146" />
<nuclide name="Cr-54" wo="0.004612692337" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom nozzle region -->
<material id="9">
<density value="2.53" units="g/cm3" />
<nuclide name="H-1" wo="0.0245014" />
<nuclide name="O-16" wo="0.1944274" />
<nuclide name="B-10" wo="7.89917e-5" />
<nuclide name="B-11" wo="3.59854e-4" />
<nuclide name="Fe-54" wo="0.030411411144" />
<nuclide name="Fe-56" wo="0.495012237964" />
<nuclide name="Fe-57" wo="0.01163454624" />
<nuclide name="Fe-58" wo="0.001578204652" />
<nuclide name="Ni-58" wo="0.047211231662" />
<nuclide name="Ni-60" wo="0.018811987544" />
<nuclide name="Ni-61" wo="0.000831139127" />
<nuclide name="Ni-62" wo="0.002694352115" />
<nuclide name="Ni-64" wo="0.000708189552" />
<nuclide name="Mn-55" wo="0.0156126" />
<nuclide name="Si-28" wo="0.007171335558" />
<nuclide name="Si-29" wo="0.000377356542" />
<nuclide name="Si-30" wo="0.0002576079" />
<nuclide name="Cr-50" wo="0.006190885148" />
<nuclide name="Cr-52" wo="0.124142524198" />
<nuclide name="Cr-53" wo="0.014348496148" />
<nuclide name="Cr-54" wo="0.003638294506" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Top nozzle region -->
<material id="10">
<density value="1.746" units="g/cm3" />
<nuclide name="H-1" wo="0.0358870" />
<nuclide name="O-16" wo="0.2847761" />
<nuclide name="B-10" wo="1.15699e-4" />
<nuclide name="B-11" wo="5.27075e-4" />
<nuclide name="Fe-54" wo="0.02644016154" />
<nuclide name="Fe-56" wo="0.43037146399" />
<nuclide name="Fe-57" wo="0.0101152584" />
<nuclide name="Fe-58" wo="0.00137211607" />
<nuclide name="Ni-58" wo="0.04104621835" />
<nuclide name="Ni-60" wo="0.0163554502" />
<nuclide name="Ni-61" wo="0.000722605975" />
<nuclide name="Ni-62" wo="0.002342513875" />
<nuclide name="Ni-64" wo="0.0006157116" />
<nuclide name="Mn-55" wo="0.0135739" />
<nuclide name="Si-28" wo="0.006234853554" />
<nuclide name="Si-29" wo="0.000328078746" />
<nuclide name="Si-30" wo="0.0002239677" />
<nuclide name="Cr-50" wo="0.005382452306" />
<nuclide name="Cr-52" wo="0.107931450781" />
<nuclide name="Cr-53" wo="0.012474806806" />
<nuclide name="Cr-54" wo="0.003163190107" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Top of Fuel Assemblies -->
<material id="11">
<density value="3.044" units="g/cm3" />
<nuclide name="H-1" wo="0.0162913" />
<nuclide name="O-16" wo="0.1292776" />
<nuclide name="B-10" wo="5.25228e-5" />
<nuclide name="B-11" wo="2.39272e-4" />
<nuclide name="Zr-90" wo="0.43313403903" />
<nuclide name="Zr-91" wo="0.09549277374" />
<nuclide name="Zr-92" wo="0.14759527104" />
<nuclide name="Zr-94" wo="0.15280552077" />
<nuclide name="Zr-96" wo="0.02511169542" />
<sab name="HH2O" xs="71t" />
</material>
<!-- Bottom of Fuel Assemblies -->
<material id="12">
<density value="1.762" units="g/cm3" />
<nuclide name="H-1" wo="0.0292856" />
<nuclide name="O-16" wo="0.2323919" />
<nuclide name="B-10" wo="9.44159e-5" />
<nuclide name="B-11" wo="4.30120e-4" />
<nuclide name="Zr-90" wo="0.3741373658" />
<nuclide name="Zr-91" wo="0.0824858164" />
<nuclide name="Zr-92" wo="0.1274914944" />
<nuclide name="Zr-94" wo="0.1319920622" />
<nuclide name="Zr-96" wo="0.0216912612" />
<sab name="HH2O" xs="71t" />
</material>
</materials>

18
tests/test_filter_energy/results_true.dat
View file

@ -1,11 +1,11 @@
k-combined:
1.005983E+00 2.248579E-02
9.903196E-01 4.279617E-02
tally 1:
2.770022E+01
1.559167E+02
4.285040E+01
3.679177E+02
5.288354E+01
5.619308E+02
1.134904E+01
2.579927E+01
2.844008E+01
1.619630E+02
4.425619E+01
3.938244E+02
5.527425E+01
6.120383E+02
9.799897E+00
1.957877E+01

9
tests/test_filter_energy/tallies.xml
View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<filter type="energy" bins="0.0 0.253e-6 1.0e-3 1.0 20.0" />
<scores>total</scores>
</tally>
</tallies>

26
tests/test_filter_energy/test_filter_energy.py
View file

@ -1,10 +1,30 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, '..')
from testing_harness import TestHarness
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterEnergyTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='energy',
bins=(0.0, 0.253e-6, 1.0e-3, 1.0, 20.0))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterEnergyTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterEnergyTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*', True)
harness = FilterEnergyTestHarness('statepoint.10.*', True)
harness.main()

181
tests/test_filter_energyout/geometry.xml
View file

@ -1,181 +0,0 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="z-cylinder" coeffs="0. 0. 0.41" />
<surface id="2" type="z-cylinder" coeffs="0. 0. 0.475" />
<surface id="3" type="z-cylinder" coeffs="0. 0. 0.56" />
<surface id="4" type="z-cylinder" coeffs="0. 0. 0.62" />
<surface id="5" type="z-cylinder" coeffs="0. 0. 187.6" />
<surface id="6" type="z-cylinder" coeffs="0. 0. 209.0" />
<surface id="7" type="z-cylinder" coeffs="0. 0. 229.0" />
<surface id="8" type="z-cylinder" coeffs="0. 0. 249.0" boundary="vacuum" />
<surface id="31" type="z-plane" coeffs="-229.0" boundary="vacuum" />
<surface id="32" type="z-plane" coeffs="-199.0" />
<surface id="33" type="z-plane" coeffs="-193.0" />
<surface id="34" type="z-plane" coeffs="-183.0" />
<surface id="35" type="z-plane" coeffs="0.0" />
<surface id="36" type="z-plane" coeffs="183.0" />
<surface id="37" type="z-plane" coeffs="203.0" />
<surface id="38" type="z-plane" coeffs="215.0" />
<surface id="39" type="z-plane" coeffs="223.0" boundary="vacuum" />
<!-- All geometry on base universe -->
<cell id="1" fill="200" surfaces=" -6 34 -35" /> <!-- Lower core -->
<cell id="2" fill="201" surfaces=" -6 35 -36" /> <!-- Upper core -->
<cell id="3" material="8" surfaces=" -7 31 -32" /> <!-- Lower core plate region -->
<cell id="4" material="9" surfaces=" -5 32 -33" /> <!-- Bottom nozzle region -->
<cell id="5" material="12" surfaces=" -5 33 -34" /> <!-- Bottom FA region -->
<cell id="6" material="11" surfaces=" -5 36 -37" /> <!-- Top FA region -->
<cell id="7" material="10" surfaces=" -5 37 -38" /> <!-- Top nozzle region -->
<cell id="8" material="7" surfaces=" -7 38 -39" /> <!-- Upper plate region -->
<cell id="9" material="4" surfaces="6 -7 32 -38" /> <!-- Downcomer -->
<cell id="10" material="5" surfaces="7 -8 31 -39" /> <!-- RPV -->
<cell id="11" material="6" surfaces="5 -6 32 -34" /> <!-- Bottom of radial reflector -->
<cell id="12" material="7" surfaces="5 -6 36 -38" /> <!-- Top of radial reflector -->
<!-- Fuel pin, cladding, cold water -->
<cell id="21" universe="1" material="1" surfaces="-1" />
<cell id="22" universe="1" material="2" surfaces="1 -2" />
<cell id="23" universe="1" material="3" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="24" universe="2" material="3" surfaces="-3" />
<cell id="25" universe="2" material="2" surfaces="3 -4" />
<cell id="26" universe="2" material="3" surfaces="4" />
<!-- Fuel pin, cladding, hot water -->
<cell id="27" universe="3" material="1" surfaces="-1" />
<cell id="28" universe="3" material="2" surfaces="1 -2" />
<cell id="29" universe="3" material="4" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="30" universe="4" material="4" surfaces="-3" />
<cell id="31" universe="4" material="2" surfaces="3 -4" />
<cell id="32" universe="4" material="4" surfaces="4" />
<!-- cell for water assembly (cold) -->
<cell id="50" universe="5" material="4" surfaces="34 -35" />
<!-- containing cell for fuel assembly -->
<cell id="60" universe="6" fill="100" surfaces="34 -35" />
<!-- cell for water assembly (hot) -->
<cell id="70" universe="7" material="3" surfaces="35 -36" />
<!-- containing cell for fuel assembly -->
<cell id="80" universe="8" fill="101" surfaces="35 -36" />
<!-- Fuel Assembly (Lower Half) -->
<lattice id="100">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
</universes>
</lattice>
<!-- Fuel Assembly (Upper Half) -->
<lattice id="101">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
</universes>
</lattice>
<!-- Core Lattice (Lower Half) -->
<lattice id="200">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
</universes>
</lattice>
<!-- Core Lattice (Upper Half) -->
<lattice id="201">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
</universes>
</lattice>
</geometry>

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