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

Merge pull request #1232 from paulromano/rotate-lattice

Browse source 
Allow cells filled with lattices to be translated/rotated
This commit is contained in:
Patrick Shriwise 2019-06-24 15:47:47 -05:00 • committed by GitHub
commit dccc9105cc
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
29 changed files with 710 additions and 348 deletions
Download patch file Download diff file Expand all files Collapse all files

13
docs/source/examples/hexagonal.rst Normal file
View file

@ -0,0 +1,13 @@
.. _notebook_hexagonal:
===========================
Modeling Hexagonal Lattices
===========================
.. only:: html
.. notebook:: ../../../examples/jupyter/hexagonal-lattice.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

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

@ -8,9 +8,9 @@ The following series of `Jupyter <https://jupyter.org/>`_ Notebooks provide
examples for how to use various features of OpenMC by leveraging the
:ref:`pythonapi`.
-----------
Basic Usage
-----------
-------------
General Usage
-------------
.. toctree::
:maxdepth: 1
@ -21,12 +21,21 @@ Basic Usage
tally-arithmetic
expansion-filters
search
triso
candu
nuclear-data
nuclear-data-resonance-covariance
cad-geom
--------
Geometry
--------
.. toctree::
:maxdepth: 1
hexagonal
triso
candu
------------------------------------
Multi-Group Cross Section Generation
------------------------------------

4
docs/source/pythonapi/model.rst
View file

@ -12,8 +12,8 @@ Convenience Functions
openmc.model.borated_water
openmc.model.cylinder_from_points
openmc.model.get_hexagonal_prism
openmc.model.get_rectangular_prism
openmc.model.hexagonal_prism
openmc.model.rectangular_prism
openmc.model.subdivide
TRISO Fuel Modeling

2
examples/jupyter/expansion-filters.ipynb
View file

@ -62,7 +62,7 @@
"source": [
"# Define surfaces used to construct regions\n",
"zmin, zmax = -10., 10.\n",
"box = openmc.model.get_rectangular_prism(10., 10., boundary_type='reflective')\n",
"box = openmc.model.rectangular_prism(10., 10., boundary_type='reflective')\n",
"bottom = openmc.ZPlane(z0=zmin, boundary_type='vacuum')\n",
"boron_lower = openmc.ZPlane(z0=-0.5)\n",
"boron_upper = openmc.ZPlane(z0=0.5)\n",

372
examples/jupyter/hexagonal-lattice.ipynb Normal file
View file

@ -0,0 +1,372 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example, we will create a hexagonal lattice and show how the orientation can be changed via the cell rotation property. Let's first just set up some materials and universes that we will use to fill the lattice."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import openmc"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"fuel = openmc.Material(name='fuel')\n",
"fuel.add_nuclide('U235', 1.0)\n",
"fuel.set_density('g/cm3', 10.0)\n",
"\n",
"fuel2 = openmc.Material(name='fuel2')\n",
"fuel2.add_nuclide('U238', 1.0)\n",
"fuel2.set_density('g/cm3', 10.0)\n",
"\n",
"water = openmc.Material(name='water')\n",
"water.add_nuclide('H1', 2.0)\n",
"water.add_nuclide('O16', 1.0)\n",
"water.set_density('g/cm3', 1.0)\n",
"\n",
"mats = openmc.Materials((fuel, fuel2, water))\n",
"mats.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With our three materials, we will set up two universes that represent pin-cells: one with a small pin and one with a big pin. Since we will be using these universes in a lattice, it's always a good idea to have an \"outer\" universe as well that is applied outside the defined lattice."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"r_pin = openmc.ZCylinder(r=0.25)\n",
"fuel_cell = openmc.Cell(fill=fuel, region=-r_pin)\n",
"water_cell = openmc.Cell(fill=water, region=+r_pin)\n",
"pin_universe = openmc.Universe(cells=(fuel_cell, water_cell))\n",
"\n",
"r_big_pin = openmc.ZCylinder(r=0.5)\n",
"fuel2_cell = openmc.Cell(fill=fuel2, region=-r_big_pin)\n",
"water2_cell = openmc.Cell(fill=water, region=+r_big_pin)\n",
"big_pin_universe = openmc.Universe(cells=(fuel2_cell, water2_cell))\n",
"\n",
"all_water_cell = openmc.Cell(fill=water)\n",
"outer_universe = openmc.Universe(cells=(all_water_cell,))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's create a hexagonal lattice using the `HexLattice` class:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"lat = openmc.HexLattice()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We need to set the `center` of the lattice, the `pitch`, an `outer` universe (which is applied to all lattice elements outside of those that are defined), and a list of `universes`. Let's start with the easy ones first. Note that for a 2D lattice, we only need to specify a single number for the pitch."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"lat.center = (0., 0.)\n",
"lat.pitch = (1.25,)\n",
"lat.outer = outer_universe"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we need to set the `universes` property on our lattice. It needs to be set to a list of lists of Universes, where each list of Universes corresponds to a ring of the lattice. The rings are ordered from outermost to innermost, and within each ring the indexing starts at the \"top\". To help visualize the proper indices, we can use the `show_indices()` helper method."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" (0, 0)\n",
" (0,11) (0, 1)\n",
"(0,10) (1, 0) (0, 2)\n",
" (1, 5) (1, 1)\n",
"(0, 9) (2, 0) (0, 3)\n",
" (1, 4) (1, 2)\n",
"(0, 8) (1, 3) (0, 4)\n",
" (0, 7) (0, 5)\n",
" (0, 6)\n"
]
}
],
"source": [
"print(lat.show_indices(num_rings=3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's set up a lattice where the first element in each ring is the big pin universe and all other elements are regular pin universes. From the diagram above, we see that the outer ring has 12 elements, the middle ring has 6, and the innermost degenerate ring has a single element."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"HexLattice\n",
"\tID =\t4\n",
"\tName =\t\n",
"\tOrientation =\ty\n",
"\t# Rings =\t3\n",
"\t# Axial =\tNone\n",
"\tCenter =\t(0.0, 0.0)\n",
"\tPitch =\t(1.25,)\n",
"\tOuter =\t3\n",
"\tUniverses \n",
" 2\n",
" 1 1\n",
"1 2 1\n",
" 1 1\n",
"1 2 1\n",
" 1 1\n",
"1 1 1\n",
" 1 1\n",
" 1\n"
]
}
],
"source": [
"outer_ring = [big_pin_universe] + [pin_universe]*11\n",
"middle_ring = [big_pin_universe] + [pin_universe]*5\n",
"inner_ring = [big_pin_universe]\n",
"lat.universes = [outer_ring, middle_ring, inner_ring]\n",
"print(lat)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's put our lattice inside a circular cell that will serve as the top-level cell for our geometry."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"outer_surface = openmc.ZCylinder(r=4.0, boundary_type='vacuum')\n",
"main_cell = openmc.Cell(fill=lat, region=-outer_surface)\n",
"geom = openmc.Geometry([main_cell])\n",
"geom.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's create a plot to see what our geometry looks like."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQAgMAAAD90d5fAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///8AAP///wCAgACerKf2AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+MGCAIxFif9eE4AAAhVSURBVHja7Z3LdeM4EEXlBUJgPgyBC9E8xxvtlQSj6BBmMc5HoWjZR2OLlkVK9Xn1IeRRd626DdDXVQ8ECiAIbjYmK/3ZOttVJutnVgGxEualv7M2m1F6wrpcRtOTtq3ASKWUnrW0iAmMNMpLL1qbAukVy2A0GiRB/KIxMmTRGfGAqcFKCNgLwoi2MIwRCxigelx7lBFxBXYk4grO8LticMTvioXhdcXkiNcVG8PnitERnyuNFeLowcBea27t+o44XHE4YnfFLLtHeg/D2opd0bLGyyH7p9mk9zFs8SpeiEV6Z7RM8XLK/mm49O5oWeLljpYhXoFo4fEKRAuPVyBacLxC0ULjFYoWGq9QtNB4xRhY/xWUBBMlGC0sXlEGEi8tWq/v77/C8SoqQ6XojViRZP8J+TcqCuCI6kpQkrcJ8k9MFEWS9y+LiSJL8nqByPHSRIGipcYrJMn+AlHalywKKElMlCYHIosiX/t6hfzyi4LdJcE7pWRBJFEUSfZXSKD76rMggijaWPL+jjYvQZSSB+FFafIgvCh9HqT3Ql7nkF9OCJBD4BBO+ZIJ4ZRvMiGc8n0mhBFFTYJtkPZxkJILoZVvciG08tpVRkjvg5i6FRoCzLBsEEr5kg2hlG+yIVsfxDD8MhD9KiOEUB6ZWeMp0dnax0AKcBWepp7tvnk1wFVwwj3Z9jEQ5Cp4EvRlHt3x6dyXtS4IODHlIEWouzsdL/+8lWRWRNht82r4qsPpdLqJ13e0PooO/JVbHLL7+E2Xv/dm2WNeBEAEp8eP3/T78p/9om0tigjDIafTLF7LpahFkQoRGtew/E2LRbVzkSBKa4McjEUUpPAVdxrkyF+7bMONBjkaiybb1ofw9c7NVITAbbgGROoeQ5C2NqQI9SLCL9pwFUizFmQLQiLdygIiVItBehBy05/v56mK1tXPIXIWsWimb4vRd1QhLQgZ5+Le5hHybTKHFLHeMIv7TSIxaJLM2rAM6W+H+OUgL1/bQS34/Pd+//MuFVYcmbVhBXI1Ywa5gKBXWHPhT7NDjFm9D2Kbnywg8NNe43TubG1NSLkvou9jbfY7EjdNx0J29CihQMireAjTJSkrEiPVW3bcvTgwnasCIQeXLQfZMWORvEo0kGPxlrtNRmboliE7Oj8SIb+tkNEGOfkhJwbyQkPuqyurqcxV7QMhQwRyoCFlXUj3ZJBmXciWSSNSW1dFSE9Xz7rjewaS2neJkKMVsrNBMseTCfLCVL//sWtknPqVF7p61hjPQxKzlQlSiJ8n5l1Tv0JCaPNkkD8T4snqJ0hjqL+XGhdnWyPkTdQ9CWKfM8KQyOwXhkTm8V8QwJHAisSnIZDQ2goKCa0SWSDu9S4QElu5myDq/DoOaXVIaDX1z4OE1uq/IGV9SPcnQeLC/yzIgSpSnzPCkHi3AkO+/2t69nuBNFqVyFPsybYAZDy5n8fjkN3Ju7PAAPHvkTBA/Ls9LJCruTJIK8SVC1shrqzeCnHNT2pBLLV907mPxMtSWZv92iADfY8pkIG7M2nI6F2RoDsyGuJfW8Ehg3+ViIwXCdn517uO1O8jIbmLahIkbXmQgTA5m7KaOrKDMQ/JWrJlILkr3H8hfyFEdR+E+n2PuxmrdCtVOsgqXX2VQavK8FsnkaiSEjGWm9wx5k5TG0ttqXHx9hOnDnupcWVB3ly6/5yJaY0pdpXFgirLHlUWcBKWoopWJWNRDYQcjEVGSJWFzudZF36etfoMSJWHND/jcdPzQDIeZlZ5LGt9wIwWOSBMMx3TICMv7g4aGa3bF/CiJaQB4sX+tSO4EUOHDMIvSttSErTngpS1Idb9Xf9DiHAfj0IRu7HPsEXRV2TcbDn9tVKRYbPlwHeugzCCMEWtbQNsL603iRtgDVt5ZQi3qGbclOwrMm6vBiDc9uqGhtDK24uMu9FTIUMEcqAhZV2I8V0HiJ8AcXsCv+Tigxhf1/EVVYTgr1D1njve+jKYq8j6Wlvv6YWtL+j1jvGEfQswc2S8QApVnRvIhWkPOcabX//0FPEQKe8SZiTSi6zhU54lM7/3G4EkHPPM2+b5IM16DPhN/8js9wopcsXIPL4DIaEViSvEcFjF0vS1lRaEhFaJrhDLASJo0ZdtMEjaUSjNWhD05JjzgJtxPE1ZC9LVhjzPCUvPc+rVE51EVjTIgSoyHtxW5Qg6/DC9hdkg8LGAcNHZNjBklActwwGHDV9zJw+/R/5Kw3mQgxASWZJbSBGcFrZ02I7PXCWzbx8BWSXpvmWskXQ7z7KNQko+pLuDrKB864YMyJNFDoI2rxF5RjrZPQNVXk2xv815DncvL+AAkAJduVMG3Jl1BARTftSSoKsRuoPKGyAbN+SkjbgKpMmF0Mf7l1xIR0Ig5XFI+zgIpDwOoRmQ8jCE+6BHyYR0DAQRBYYwkkCiwJBNAAJ3KyykyYPwn/Ap+sVoV9+xEEB5dNBidYeVB6rxDFCUkCSQKFhK1AmQtAxPkCQvt5cYWWm3/Em4kgPpREiSKKIkWaLIjBxRtE8nlgxIp0BSRFEkyRFFY2SIon9itMQhmiQpoqiS1PmAbZ1P8Vb5qHCVzyPX+dBziUH0BpwgCiRJnc+I1/kgepVPu1f5SH0kXmi0QvGCoxWIFx6tQLzwaG38/ZeF4Y2XJVpu6Q2y++NlY/ikN8nujZcxWi7pbbI7XTE74nDF7ohDeqvsZ7NCPAyrKy5HrK74GDZXnI7YXPEyLK64HbG44mfgrgQcwV2JMNAerA1BsB7M02uZAxZlINp3YYgesHCwkIBlMLQW1qZAZFm6HIZISWMI4qeIrlBSGUzEEmM1GdHG2mzG5v5+WQFxi1kJ8WnFJcZ/bkPEpRmBEaUAAAAldEVYdGRhdGU6Y3JlYXRlADIwMTktMDYtMDhUMDc6NDk6MjItMDU6MDA/X9I4AAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE5LTA2LTA4VDA3OjQ5OjIyLTA1OjAwTgJqhAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = openmc.Plot.from_geometry(geom)\n",
"p.color_by = 'material'\n",
"p.colors = colors = {\n",
" water: 'blue',\n",
" fuel: 'olive',\n",
" fuel2: 'yellow'\n",
"}\n",
"p.to_ipython_image()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this point, if we wanted to simulate the model, we would need to create an instance of `openmc.Settings`, export it to XML, and run."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Lattice orientation\n",
"\n",
"Now let's say we want our hexagonal lattice orientated such that two sides of the lattice are parallel to the x-axis. This can be achieved by two means: either we can rotate the cell that contains the lattice, or we can can change the `HexLattice.orientation` attribute. By default, the `orientation` is set to \"y\", indicating that two sides of the lattice are parallel to the y-axis, but we can also change it to \"x\" to make them parallel to the x-axis."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQAgMAAAD90d5fAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///8AAP+AgAD//wDoUCWoAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+MGCAIxFif9eE4AAAeGSURBVHja7Z1NkuMoEIWrFhyB++gILCxNhDZezUaX8CnqCF607tNH6WWFpmxXWfDEX0KSdmucO4dxfmQ+hARG8PZGMtVfzdB+RbLeMgFEI8x7v7GOm6F6jxlehu69dhBgsFJUHzS2jEUYbJT3PmodC6RPGAdDpyAM4qsUg0OWNKM+YclkMSTsPYdR28LyGHUJy1C9XvtcRk0o2YHUhJLPKA+FEEh5KBRGaSikQEpDoTHKQiEGUhaKpkIKerDMXsu2rn0gBaEUBEIPhSx7ifQlDGorLsoWNV8Fsl+MJn0Zg5YvVQqhSF+YLVK+CmW/WL70xdmi5Ks4W4R8VWQrP18V2crPV0W2svNVla3cfFVlKzdfVdnKzVcdI6//qpQkT5TKbOXlq5aRk6/qbOXkS9VD0o24WpIcUeoZaVEYJEmLojggKVE0ByQlCgcjJQqLJClRFA8kLormgcRF4WHERWGSJC6K4oLERNFckJgoXIyYKGySxERRfJCwKJoPEhaFjxERRQLCqHtYecUJCSmvOSEh5TkZIVFYJQmJIgJRvBC/8poX4lc+UHiY3c/TL/fz/BH4IQUyz47XYXa9TvNMgAR0H2fX6+RCv5jz2f9Ln/LKX3QCL18f7apf6vDL/0uf8jqULcfLpeZ2ZBNALTtkQ65OLS+XmtuRzQBNQPosyASRRSAe5SO6214AcqvD2f/brhQyu5ERIcpfcAIvALnVIbt5aX7I4TGQnh+yaV7vLSDdIyCqBQSblw6Uq7lONsqLQPpMCKVb2TSvULGaDhIh4cchcIpdfeR+0mPzCkPgzog3rcidcQNRwXJ4j4eaDzFJoA3rYLkB0jGC05gk0LzCkK98uTUFCcZItgDSNzJxCPMD/WqdNES1ghhpiG4FOUhDWjGcNiwBadaC7TYsAlHtIEYWottBDrKQdgyrDe8G0vAyWS8UEYhqCTG7g+iWkIMkpCXjfqHsBhK6FmEANcDT+/QRLX63LgaZYEwA45IRhhBTaAgRgwwwuhlhMAQjrCE4GOoi1yIOBmGsOMCwbgoO60wEAsPaAUa9I0DDsx8xyOyBWF4mD2SOQLTvqwFmAUbwAjMRWNyyQwqyepl8kDWy8Skg3kyOXsjdy+CFnL2uUpDVy+yDrJFN/xvIe2tI9wyQM3gBSF7r6qL9IxPEBCF8F2MS8gFeAHIOFs+D9D4I9sII8Xsyb9F/5Vi6+mvnFYBMnjvjWvNbZJHieZARWsu0vf3+ihTPg+A06Qg1nbYPEgFHMci4eSSymRevkN5AIFdI39xekGeENB1f36x7QZ4QotpDzAvygjwR5J9l+W1/Pi6fzvfL8mcL0WSGQzl+fbQpl6+RcqBCTovrdXGhx+vnSsg1EMvrzekKvX0NoVAhR/CyuNDvOrgqkSHfThfX6R16AmgR5Mfpj5fvwO5VXxZfvoiQI3g5uZHd6/BZAzmBl8WN7F6HpQayuF7uNf+O7ARQFsha8884hMJYa37zcgLoApH9GAlyBC8AserwmYCM4TEBQhY3MgJkDo9uTuAFIFYdljgEx2mTBQWIVfNrZPkQHHHao97F9YIQqw5u89pAYOw8NIHMkVkAgNjp+ayBOJMmXBD/9M8NYqfn4gUh9td/KJBJDnJ+QV6Qv7UJN7kYRboVmQ5SpKsXuWmJ3H5FHiRkHok29m8fgpzcmtdAVit+TNUUCjg9Qs1Pbh1+7BmHDie35m0GQUe35m2GcyIDU5Ehtsxkgci0h8gETuFUlKJByibVqJACe0FekJaQ/fwT9IKQIPv57/cFoUF04CvWhRghiMSSEpHFMSLLfHgXLKkIZPVSufTKDxFZRCayHG70Qu5eqAv7RJYoPnyxJRtkP6tsHw45gxeAEJZXa983f+Vq9AjkA7wA5BwsvoEobybZuvqHv1Axzqyvhoi85PLA13WYXzzaz3taMhDdkiH6FqBqCdnf65/7ee93P+9iy0B0O4bwm/6qHcTIQvazI8Z+dinZ0c4xqhXESEP2s8PSfna9ktkkLAxh3O5MBcuN4AWGcZSN20S2oBPZTE9kW0CZDQ51JgSmoJ5wP0gVKIezZQDBuTfXDED2s9uoyOasO9rLVvFDzAYisomxyHbMIhtLy2yR7Ycwb/at/EVxmhScjpTGJbQBu8hW8jKb4uue1fzb+yteiPFC9nOug8gxGDIHeihOiAlARA5Z2c+ZNDJH+Cg+iAlCRI5VEjkgSuaoK8UFMRGIyPFjIgepyRwJp3ggJgoROaZP5MBBmaMTFQfEJCAix1mKHMwpc8SoqoekJBE69lXkAFuZo3hFDhUWOR5Z5qBnVQdJN2AGUbIkkTlGXOZAdJGj3UUOqa/JV262qvKVna2KfOVnqyJf+dl6K++/KIzSfFGyVSw9QfbyfNEYZdKTZC/NFzFbRdLTZC8MhRxIQSj0QAqkp8p+NSqkhEENpSgQaihlDFoohYHQQillUEIpDoQSSjkjP5SKQPJDqWHk9mBdFSSvByvptcgJq2XkaG+qIemEVScrJ2EcjFQL61ggcVkMDyNKYWNExGcRPUFhZQQyxpirm3naWMfNeNteLw0QiGmEuJgqEuM/rS4lqiAX++MAAAAldEVYdGRhdGU6Y3JlYXRlADIwMTktMDYtMDhUMDc6NDk6MjItMDU6MDA/X9I4AAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE5LTA2LTA4VDA3OjQ5OjIyLTA1OjAwTgJqhAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Change the orientation of the lattice and re-export the geometry\n",
"lat.orientation = 'x'\n",
"geom.export_to_xml()\n",
"\n",
"# Run OpenMC in plotting mode\n",
"p.to_ipython_image()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When we change the orientation to 'x', you can see that the first universe in each ring starts to the right along the x-axis. As before, the universes are defined in a clockwise fashion around each ring. To see the proper indices for a hexagonal lattice in this orientation, we can again call `show_indices` but pass an extra orientation argument:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" (0, 8) (0, 9) (0,10)\n",
"\n",
" (0, 7) (1, 4) (1, 5) (0,11)\n",
"\n",
"(0, 6) (1, 3) (2, 0) (1, 0) (0, 0)\n",
"\n",
" (0, 5) (1, 2) (1, 1) (0, 1)\n",
"\n",
" (0, 4) (0, 3) (0, 2)\n"
]
}
],
"source": [
"print(lat.show_indices(3, orientation='x'))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hexagonal prisms\n",
"\n",
"OpenMC also contains a convenience function that can create a hexagonal prism representing the interior region of six surfaces defining a hexagon. This can be useful as a bounding surface of a hexagonal lattice. For example, if we wanted the outer boundary of our geometry to be hexagonal, we could change the `region` of the main cell:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQBAMAAABykSv/AAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAAElBMVEX///+TUVD//wCAgAA4vPIAAP/pte6jAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+MGCAIxFif9eE4AAAn2SURBVHja7d3tceM4DAbgsANJHVh0BVIFl7E7WPffymV346wlEh8UCVDGAP9u3jkvHhP+JCN/fHh5eXl5eXl5eXl5eXl5eXl5eXl5eXl5eXl5eXl5WatRtPQcQRYyOKS0ZB16syW8IHpLYgYi7VCbLSuQIA8ZHFJS8g6d2VJYEJ0lMQPRcKjMlhVI0IEMDuGWjkN+trILMl/QniYinnssSQ4yz6hkIuOcRBoCdDLjjWISKFaH/GkUlnzHoASMZR0BahSUMONUMuguyI8DaJUdpxJdyIy1solzkgmLJR0BazQj2cQz6kjjQRGybTSRTEXxrAjBHftWp8J4L1GDpJ3MRLyRELEkJFCOV0k2vmCOvWTQWZBsoy+SQ/FWogOZ8VaAeMbvhlkFEhiOpwSMZ9SxlQwKELjRP5LpcDwrQHiO361OFfGrRByCdzIT8YWIpSGBuSDVJT1b3MlqCJFZEisQ5nNWW4jEbJmBwK+GcV2xvr7iBcuTWPbFHX5ZX38X3GpkxNs7QvbFHYSsKyqJrHjVg4wA5NkJNF3MeAEg7WcLgqwruiSRGa9akABA/nWSv8/Z8QJAWs/WQUhc8RXTh4wAZMU7jex4BSCNZyschKwrumIcSNslgSCR6LQgXlQgYz9I29myAgk9IS1nywxk7AlpOFuhL6TdksAQhdeRlpCxL6TdbB2HRHasAQkIJKL3eIM3jS1nC4NQn0dwJ+PzSEvIiEEi2ij1wYr+hNhwtgIKoT6z407yM3vLJSEgEW30KSHizd0gBRlxyF/JPKOtEvF2OTP/oAqkeQlBwhkgLWZrPAOkxZJYgYRzQOpnywwkd6uTLOQiMlv6DpmzmvpPvr9LYLa6OCReE/s42r8DDp0czT+UJJBJC3JpC+nmyEhaQhQdqaThZKk6EknNbHV6oD+r3ZL0dbTbhAudHc32RTeQqQfk0gbS3bGTtIB0cmwlDSarm2MjOTpbnR/oz6pfklMsSIPjs+eYrAazZQZyhuesHeTQbIUzQo4sCQLBT5XOET9VSsTJrddCwJf1iO3rzN87Hofj9NZrX9whSMR3qGZ8h4qKM7deCQHfMB7YFEwbJbccIUj5bEGQ7D4sO2ZvAr/ElZARgLD3m5cDcfbWN5Di2QoAhH8CYD0QZ299CyldEgaEf7gkC1m4cR1kPCukdLYgyMqHLMVx/tarIIEDWfBOCMjKjXeQstkyAxnPCymarXBmSMmSmIGMLMiMt7IUx6xnraLZsgIJMCTyITPe6cKNEwh/tniQhej0QBwbQ0YYMrMhC9EpO04g7NkKPMiMt7IciBlv40uWBIVEvJMZd1Jx7taPQ0YMwj02uhyKc7eeQrizhUNmvBMqjujdkIsPQwIBmdFGyThijlycgfBma6QgypWB8JbECiS8A4QzW2Yg4ztAOLNlBZKbrMq90OsnHt/w//2W64ieLQHH7YZKyDgr6fDk+9XJjYgRydfdcPvvACSIOBDJ3xiUXP/EGQk1W0HGAUquvPi/YoiUA2j1yo1TSSmkzvGvk+z0vMRZyUtcCAlijpxkE2ckm7hstgQdGUlRfCtaEklH0mphfCuA7CZrqoMkndyIeCNJ7ob9A35gQ9o47o+v+pVKnu29xi+S60v8yEowSHvHH8ZLr/tGmfEjJ2EuSFvHvtV9o1Sckww8SJ3junc8W/0EGt1K0viRPuBhSPMnrEemlU+g0VfoLb0bnlDWbMk6XiS5Rsk4lTAma6qD5Dv5mZ5b3vkdfwJ3Q/owGWiIyIL8PI0CTipOlgSCSD7SN61ATiJOH+/0c1YLyANsBXJS8a89ZDgD5AF2Wgtp/OR7xzvNxg/c+eA9ATeGPPBO70fiXxwI/LK+lm4EcCAPuNMCyFACWYu3Zq5gJ/86PRQ/OJARgDD30tYyyB3ptAAyEguSvThQHpLbpr3CnTx4ECj+tYcMbMiR/eaukBGA8M8qFEEeSKclkLEV5CVeTgEJACRmO6Uhd6LTg/EjgQxtIPF0kLEesnaBjOiCvBNkcIhDNCBmHux2nn7NvCCaeYti501jACARnawTvo2HIG/3wcrMR107Xz6AkHf7OqjiC7q4+28O5H4kTiAfbSFzDnIHO+FAar4yDW0h/b7EVobcy2MmxMxGj52ttxfIVCcBWvnuRHwz1Mz2tJ0DA6G15A40Kn2EYwOZ6iTJnbppVPhQjZljTnYOnm0hb3wU0MzhTDvHZYOg5BOP2x5gNnOk3M4h/yAlIZaMiFPHUAx50z+EMfOnSXb+WCwISD7xWObP98z8QaUdSHgHCGOy7EDM/GG+HUigIMi2DyOOpZsUGQhrsuxcTsTOBV4CBvnZM1zwTs5wyR0Ugm+7n+siSGYuS2XnQmE8yBtcus3MxfTsQAILsuCdEBC/cmYJxMxFWe1cJtcMxMylpO1c3DucFVI4WXYugG/mJwns/EgEBHm7n+0w80Mqdn7aJkCQd/uxIQSiW7UQMz/IZecn0sxAzPyMoJ0fdjzHbNVP1jmet6qfs3aQ9/452tBdUvuGMQd5559sNvMj2nZ+1nw7W2/8Q/PwV6f6jpoFwfYX1B11kNBNkjiqJiuzwasFGdtCMocHOjnqJit7wqaPoxYSukgy/2jlZPU59HTJ/KO1jh5HHHOOekhutoSHK/cvVk+WHYjA0ewjkHqHHYj+E7DMZHV4TZRZEBoS8X2CmREvKpCAQyK+czOz4o1EaLIoCL7Vxj2UqgHB3wFzNwUBaC4WmiziQ8mB/eacc8UgjRYEhUS8U/75gkUBgv2V5Yp3GtnxikBaOSog64qumDYk9IU0mywEEolOC+JFAfLRF9LOwbuwphSk4YLYgbCuECoFaemwAwn9IE0nC4QovI60hXCuELrgnS4HIW0djCuECr1pbLwgRyHU23h9CHiFUPwepz5YxZxTcrLoS51Cn2WZ8aIGCQCE+sweWTH4mb35ZFGXOoUaJQ+lxnS9ZCHIBeoifuo0rvix0iQWnayWF9akSnZB7ED0TgbKTpYhSNCHiExWwfFZIr4QsfCC8GfrgkOJWHyy2M9bF3zJiFj8OeuDe6nTCwElZlP4Zf1vjQzJhYASj7LNPyHlYF3qlIAS8UUHEmgJEb80SjsGMQj4h6/ZyZiIRimH3ILQx2cvOJSIRb/Q2lbAJbt7dN8qEYt+VUpARrTRPTSJJzSWhKCXOk0dW0kmnrBY0oFejGfMFTtO7wZRSDJb/yTjiLZKxKljEIXA194CGn3GFygHY1kHeFZzBGtCHX8luVgYEvKtII3+lhCx2AmtQsiINvrVKhFLHWLEa1QqaUd+SdqX+ILYgSjNlrzDDiRoOAaHFJQGRMOhsSQqC2IHojBbOg47kCDtGBxSWNIQLYf0kqgtiB2I8GzpOby8vLy8vLy8vLy8vLy8vLy8vLy8vLy8vLy8vLy8lOp/nb6TXPV6EWgAAAAldEVYdGRhdGU6Y3JlYXRlADIwMTktMDYtMDhUMDc6NDk6MjItMDU6MDA/X9I4AAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE5LTA2LTA4VDA3OjQ5OjIyLTA1OjAwTgJqhAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"main_cell.region = openmc.model.hexagonal_prism(\n",
" edge_length=3*lat.pitch[0],\n",
" orientation='x',\n",
" boundary_type='vacuum'\n",
")\n",
"geom.export_to_xml()\n",
"\n",
"# Run OpenMC in plotting mode\n",
"p.color_by = 'cell'\n",
"p.to_ipython_image()"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.0"
}
},
"nbformat": 4,
"nbformat_minor": 1
}

4
examples/jupyter/pincell.ipynb
View file

@ -903,8 +903,8 @@
}
],
"source": [
"box = openmc.get_rectangular_prism(width=pitch, height=pitch,\n",
" boundary_type='reflective')\n",
"box = openmc.rectangular_prism(width=pitch, height=pitch,\n",
" boundary_type='reflective')\n",
"type(box)"
]
},

12
include/openmc/particle.h
View file

@ -9,6 +9,7 @@
#include <memory> // for unique_ptr
#include <sstream>
#include <string>
#include <vector>
#include "openmc/constants.h"
#include "openmc/position.h"
@ -38,7 +39,13 @@ constexpr double CACHE_INVALID {-1.0};
// Class declarations
//==============================================================================
struct LocalCoord {
class LocalCoord {
public:
void rotate(const std::vector<double>& rotation);
//! clear data from a single coordinate level
void reset();
Position r; //!< particle position
Direction u; //!< particle direction
int cell {-1};
@ -48,9 +55,6 @@ struct LocalCoord {
int lattice_y {-1};
int lattice_z {-1};
bool rotated {false}; //!< Is the level rotated?
//! clear data from a single coordinate level
void reset();
};
//==============================================================================

10
include/openmc/position.h
View file

@ -3,6 +3,7 @@
#include <array>
#include <cmath> // for sqrt
#include <iostream>
#include <stdexcept> // for out_of_range
#include <vector>
@ -55,13 +56,16 @@ struct Position {
//! Dot product of two vectors
//! \param[in] other Vector to take dot product with
//! \result Resulting dot product
inline double dot(Position other) {
inline double dot(Position other) const {
return x*other.x + y*other.y + z*other.z;
}
inline double norm() {
inline double norm() const {
return std::sqrt(x*x + y*y + z*z);
}
//! Rotate the position based on a rotation matrix
Position rotate(const std::vector<double>& rotation) const;
// Data members
double x = 0.;
double y = 0.;
@ -91,6 +95,8 @@ inline bool operator==(Position a, Position b)
inline bool operator!=(Position a, Position b)
{return a.x != b.x || a.y != b.y || a.z != b.z;}
std::ostream& operator<<(std::ostream& os, Position a);
//==============================================================================
//! Type representing a vector direction in Cartesian coordinates
//==============================================================================

2
openmc/__init__.py
View file

@ -32,6 +32,6 @@ from openmc.polynomial import *
from . import examples
# Import a few convencience functions that used to be here
from openmc.model import get_rectangular_prism, get_hexagonal_prism
from openmc.model import rectangular_prism, hexagonal_prism
__version__ = '0.10.0'

4
openmc/cell.py
View file

@ -227,10 +227,6 @@ class Cell(IDManagerMixin):
@rotation.setter
def rotation(self, rotation):
if not isinstance(self.fill, openmc.Universe):
raise TypeError('Cell rotation can only be applied if the cell '
'is filled with a Universe.')
cv.check_type('cell rotation', rotation, Iterable, Real)
cv.check_length('cell rotation', rotation, 3)
self._rotation = np.asarray(rotation)

3
openmc/mgxs/__init__.py
View file

@ -7,6 +7,7 @@ from openmc.mgxs.mdgxs import *
GROUP_STRUCTURES = {}
"""Dictionary of commonly used energy group structures:
- "CASMO-X" (where X is 2, 4, 8, 16, 25, 40 or 70) from the CASMO_ lattice
physics code
- "XMAS-172_" designed for LWR analysis ([SAR1990]_, [SAN2004]_)
@ -19,7 +20,7 @@ GROUP_STRUCTURES = {}
.. _XMAS-172: https://www-nds.iaea.org/wimsd/energy.htm
.. _SHEM-361: https://www.polymtl.ca/merlin/libraries.htm
.. _activation: https://fispact.ukaea.uk/wiki/Keyword:GETXS
.. _CCFE-709: https://fispact.ukaea.uk/wiki/CCFE-709_group_structure
.. _CCFE-709: https://fispact.ukaea.uk/wiki/CCFE-709_group_structure
.. _UKAEA-1102: https://fispact.ukaea.uk/wiki/UKAEA-1102_group_structure
.. [SAR1990] Sartori, E., OECD/NEA Data Bank: Standard Energy Group Structures
of Cross Section Libraries for Reactor Shielding, Reactor Cell and Fusion

31
openmc/model/funcs.py
View file

@ -3,6 +3,7 @@ from collections.abc import Iterable
from math import sqrt
from numbers import Real
from functools import partial
from warnings import warn
from openmc import XPlane, YPlane, Plane, ZCylinder, Quadric
from openmc.checkvalue import check_type, check_value
@ -101,10 +102,14 @@ def borated_water(boron_ppm, temperature=293., pressure=0.1013, temp_unit='K',
return out
def get_rectangular_prism(width, height, axis='z', origin=(0., 0.),
boundary_type='transmission', corner_radius=0.):
def rectangular_prism(width, height, axis='z', origin=(0., 0.),
boundary_type='transmission', corner_radius=0.):
"""Get an infinite rectangular prism from four planar surfaces.
.. versionchanged:: 0.11
This function was renamed from `get_rectangular_prism` to
`rectangular_prism`.
Parameters
----------
width: float
@ -202,10 +207,21 @@ def get_rectangular_prism(width, height, axis='z', origin=(0., 0.),
return prism
def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
boundary_type='transmission', corner_radius=0.):
def get_rectangular_prism(*args, **kwargs):
warn("get_rectangular_prism(...) has been renamed rectangular_prism(...). "
"Future versions of OpenMC will not accept get_rectangular_prism.",
FutureWarning)
return rectangular_prism(*args, **kwargs)
def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
boundary_type='transmission', corner_radius=0.):
"""Create a hexagon region from six surface planes.
.. versionchanged:: 0.11
This function was renamed from `get_hexagonal_prism` to
`hexagonal_prism`.
Parameters
----------
edge_length : float
@ -346,6 +362,13 @@ def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
return prism
def get_hexagonal_prism(*args, **kwargs):
warn("get_hexagonal_prism(...) has been renamed hexagonal_prism(...). "
"Future versions of OpenMC will not accept get_hexagonal_prism.",
FutureWarning)
return hexagonal_prism(*args, **kwargs)
def cylinder_from_points(p1, p2, r, **kwargs):
"""Return cylinder defined by two points passing through its center.

110
src/geometry.cpp
View file

@ -173,32 +173,19 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
//! Found a lower universe, update this coord level then search the next.
// Set the lower coordinate level universe.
p->coord_[p->n_coord_].universe = c.fill_;
auto& coord {p->coord_[p->n_coord_]};
coord.universe = c.fill_;
// Set the position and direction.
p->coord_[p->n_coord_].r = p->coord_[p->n_coord_-1].r;
p->coord_[p->n_coord_].u = p->coord_[p->n_coord_-1].u;
coord.r = p->r_local();
coord.u = p->u_local();
// Apply translation.
p->coord_[p->n_coord_].r -= c.translation_;
coord.r -= c.translation_;
// Apply rotation.
if (!c.rotation_.empty()) {
Position r = p->coord_[p->n_coord_].r;
p->coord_[p->n_coord_].r.x = r.x*c.rotation_[3] + r.y*c.rotation_[4]
+ r.z*c.rotation_[5];
p->coord_[p->n_coord_].r.y = r.x*c.rotation_[6] + r.y*c.rotation_[7]
+ r.z*c.rotation_[8];
p->coord_[p->n_coord_].r.z = r.x*c.rotation_[9] + r.y*c.rotation_[10]
+ r.z*c.rotation_[11];
Direction u = p->coord_[p->n_coord_].u;
p->coord_[p->n_coord_].u.x = u.x*c.rotation_[3] + u.y*c.rotation_[4]
+ u.z*c.rotation_[5];
p->coord_[p->n_coord_].u.y = u.x*c.rotation_[6] + u.y*c.rotation_[7]
+ u.z*c.rotation_[8];
p->coord_[p->n_coord_].u.z = u.x*c.rotation_[9] + u.y*c.rotation_[10]
+ u.z*c.rotation_[11];
p->coord_[p->n_coord_].rotated = true;
coord.rotate(c.rotation_);
}
// Update the coordinate level and recurse.
@ -211,26 +198,37 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
Lattice& lat {*model::lattices[c.fill_]};
// Determine lattice indices.
auto i_xyz = lat.get_indices(p->r_local(), p->u_local());
// Set the position and direction.
auto& coord {p->coord_[p->n_coord_]};
coord.r = p->r_local();
coord.u = p->u_local();
// Store lower level coordinates.
Position r = lat.get_local_position(p->r_local(), i_xyz);
p->coord_[p->n_coord_].r = r;
p->coord_[p->n_coord_].u = p->u_local();
// Apply translation.
coord.r -= c.translation_;
// Apply rotation.
if (!c.rotation_.empty()) {
coord.rotate(c.rotation_);
}
// Determine lattice indices.
auto i_xyz = lat.get_indices(coord.r, coord.u);
// Get local position in appropriate lattice cell
coord.r = lat.get_local_position(coord.r, i_xyz);
// Set lattice indices.
p->coord_[p->n_coord_].lattice = c.fill_;
p->coord_[p->n_coord_].lattice_x = i_xyz[0];
p->coord_[p->n_coord_].lattice_y = i_xyz[1];
p->coord_[p->n_coord_].lattice_z = i_xyz[2];
coord.lattice = c.fill_;
coord.lattice_x = i_xyz[0];
coord.lattice_y = i_xyz[1];
coord.lattice_z = i_xyz[2];
// Set the lower coordinate level universe.
if (lat.are_valid_indices(i_xyz)) {
p->coord_[p->n_coord_].universe = lat[i_xyz];
coord.universe = lat[i_xyz];
} else {
if (lat.outer_ != NO_OUTER_UNIVERSE) {
p->coord_[p->n_coord_].universe = lat.outer_;
coord.universe = lat.outer_;
} else {
std::stringstream err_msg;
err_msg << "Particle " << p->id_ << " is outside lattice "
@ -297,27 +295,32 @@ find_cell(Particle* p, bool use_neighbor_lists)
void
cross_lattice(Particle* p, const BoundaryInfo& boundary)
{
auto& lat {*model::lattices[p->coord_[p->n_coord_-1].lattice]};
auto& coord {p->coord_[p->n_coord_ - 1]};
auto& lat {*model::lattices[coord.lattice]};
if (settings::verbosity >= 10 || simulation::trace) {
std::stringstream msg;
msg << " Crossing lattice " << lat.id_ << ". Current position ("
<< p->coord_[p->n_coord_-1].lattice_x << ","
<< p->coord_[p->n_coord_-1].lattice_y << ","
<< p->coord_[p->n_coord_-1].lattice_z << ")";
<< coord.lattice_x << "," << coord.lattice_y << ","
<< coord.lattice_z << "). r=" << p->r();
write_message(msg, 1);
}
// Set the lattice indices.
p->coord_[p->n_coord_-1].lattice_x += boundary.lattice_translation[0];
p->coord_[p->n_coord_-1].lattice_y += boundary.lattice_translation[1];
p->coord_[p->n_coord_-1].lattice_z += boundary.lattice_translation[2];
std::array<int, 3> i_xyz {p->coord_[p->n_coord_-1].lattice_x,
p->coord_[p->n_coord_-1].lattice_y,
p->coord_[p->n_coord_-1].lattice_z};
coord.lattice_x += boundary.lattice_translation[0];
coord.lattice_y += boundary.lattice_translation[1];
coord.lattice_z += boundary.lattice_translation[2];
std::array<int, 3> i_xyz {coord.lattice_x, coord.lattice_y, coord.lattice_z};
// Set the new coordinate position.
p->r_local() = lat.get_local_position(p->coord_[p->n_coord_-2].r, i_xyz);
const auto& upper_coord {p->coord_[p->n_coord_ - 2]};
const auto& cell {model::cells[upper_coord.cell]};
Position r = upper_coord.r;
r -= cell->translation_;
if (!cell->rotation_.empty()) {
r = r.rotate(cell->rotation_);
}
p->r_local() = lat.get_local_position(r, i_xyz);
if (!lat.are_valid_indices(i_xyz)) {
// The particle is outside the lattice. Search for it from the base coords.
@ -362,9 +365,10 @@ BoundaryInfo distance_to_boundary(Particle* p)
// Loop over each coordinate level.
for (int i = 0; i < p->n_coord_; i++) {
Position r {p->coord_[i].r};
Direction u {p->coord_[i].u};
Cell& c {*model::cells[p->coord_[i].cell]};
const auto& coord {p->coord_[i]};
Position r {coord.r};
Direction u {coord.u};
Cell& c {*model::cells[coord.cell]};
// Find the oncoming surface in this cell and the distance to it.
auto surface_distance = c.distance(r, u, p->surface_);
@ -372,10 +376,9 @@ BoundaryInfo distance_to_boundary(Particle* p)
level_surf_cross = surface_distance.second;
// Find the distance to the next lattice tile crossing.
if (p->coord_[i].lattice != C_NONE) {
auto& lat {*model::lattices[p->coord_[i].lattice]};
std::array<int, 3> i_xyz {p->coord_[i].lattice_x, p->coord_[i].lattice_y,
p->coord_[i].lattice_z};
if (coord.lattice != C_NONE) {
auto& lat {*model::lattices[coord.lattice]};
std::array<int, 3> i_xyz {coord.lattice_x, coord.lattice_y, coord.lattice_z};
//TODO: refactor so both lattice use the same position argument (which
//also means the lat.type attribute can be removed)
std::pair<double, std::array<int, 3>> lattice_distance;
@ -384,8 +387,13 @@ BoundaryInfo distance_to_boundary(Particle* p)
lattice_distance = lat.distance(r, u, i_xyz);
break;
case LatticeType::hex:
Position r_hex {p->coord_[i-1].r.x, p->coord_[i-1].r.y,
p->coord_[i].r.z};
auto& cell_above {model::cells[p->coord_[i-1].cell]};
Position r_hex {p->coord_[i-1].r};
r_hex -= cell_above->translation_;
if (coord.rotated) {
r_hex = r_hex.rotate(cell_above->rotation_);
}
r_hex.z = coord.r.z;
lattice_distance = lat.distance(r_hex, u, i_xyz);
break;
}

13
src/particle.cpp
View file

@ -33,6 +33,14 @@ namespace openmc {
// LocalCoord implementation
//==============================================================================
void
LocalCoord::rotate(const std::vector<double>& rotation)
{
this->r = this->r.rotate(rotation);
this->u = this->u.rotate(rotation);
this->rotated = true;
}
void
LocalCoord::reset()
{
@ -41,6 +49,7 @@ LocalCoord::reset()
lattice = C_NONE;
lattice_x = 0;
lattice_y = 0;
lattice_z = 0;
rotated = false;
}
@ -337,9 +346,7 @@ Particle::transport()
// If next level is rotated, apply rotation matrix
const auto& m {model::cells[coord_[j].cell]->rotation_};
const auto& u {coord_[j].u};
coord_[j + 1].u.x = m[3]*u.x + m[4]*u.y + m[5]*u.z;
coord_[j + 1].u.y = m[6]*u.x + m[7]*u.y + m[8]*u.z;
coord_[j + 1].u.z = m[9]*u.x + m[10]*u.y + m[11]*u.z;
coord_[j + 1].u = u.rotate(m);
} else {
// Otherwise, copy this level's direction
coord_[j+1].u = coord_[j].u;

17
src/position.cpp
View file

@ -84,4 +84,21 @@ Position::operator-() const
return {-x, -y, -z};
}
Position
Position::rotate(const std::vector<double>& rotation) const
{
return {
x*rotation[3] + y*rotation[4] + z*rotation[5],
x*rotation[6] + y*rotation[7] + z*rotation[8],
x*rotation[9] + y*rotation[10] + z*rotation[11]
};
}
std::ostream&
operator<<(std::ostream& os, Position r)
{
os << "(" << r.x << ", " << r.y << ", " << r.z << ")";
return os;
}
} // namespace openmc

8
tests/regression_tests/lattice_hex_coincident/test.py
View file

@ -103,10 +103,10 @@ class HexLatticeCoincidentTestHarness(PyAPITestHarness):
inf_mat_univ = openmc.Universe(cells=[inf_mat,])
# a hex surface for the core to go inside of
hexprism = openmc.model.get_hexagonal_prism(edge_length=edge_length,
origin=(0.0, 0.0),
boundary_type = 'reflective',
orientation='x')
hexprism = openmc.model.hexagonal_prism(edge_length=edge_length,
origin=(0.0, 0.0),
boundary_type = 'reflective',
orientation='x')
pincell_only_lattice = openmc.HexLattice(name="regular fuel assembly")
pincell_only_lattice.center = (0., 0.)

8
tests/regression_tests/lattice_hex_x/test.py
View file

@ -139,10 +139,10 @@ class HexLatticeOXTestHarness(PyAPITestHarness):
# a hex surface for the core to go inside of
hexprism = openmc.model.get_hexagonal_prism(edge_length=edge_length,
origin=(0.0, 0.0),
boundary_type='reflective',
orientation='x')
hexprism = openmc.model.hexagonal_prism(edge_length=edge_length,
origin=(0.0, 0.0),
boundary_type='reflective',
orientation='x')
region = hexprism & +fuel_bottom & -fuel_top
inf_mat = openmc.Cell(cell_id=12)

156
tests/regression_tests/lattice_mixed/geometry.xml
View file

@ -1,156 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<!--
This geometry approximately describes a Russian-style VVER-1000 fuel assembly.
-->
<geometry>
<surface id="001" type="z-plane" coeffs="-10000"/>
<!-- Universe 011: Fuel pellet -->
<surface id="011" type="z-cylinder" coeffs="0.0 0.0 0.06"/>
<surface id="012" type="z-cylinder" coeffs="0.0 0.0 0.378"/>
<surface id="013" type="z-cylinder" coeffs="0.0 0.0 0.3865"/>
<surface id="014" type="z-cylinder" coeffs="0.0 0.0 0.455"/>
<surface id="015" type="sphere" coeffs="0.0 0.0 -2.5639 2.0"/>
<surface id="016" type="sphere" coeffs="0.0 0.0 2.5639 2.0"/>
<cell id="011" universe="011" material="void" region="-011"/>
<cell id="012" universe="011" material="13" region="011 -012 015 016"/>
<cell id="013" universe="011" material="void" region="011 -012 -015"/>
<cell id="014" universe="011" material="void" region="011 -012 -016"/>
<cell id="015" universe="011" material="void" region="012 -013"/>
<cell id="016" universe="011" material="21" region="013 -014"/>
<cell id="017" universe="011" material="01" region="014"/>
<!-- Universe 051: Central guide tube -->
<surface id="051" type="z-cylinder" coeffs="0.0 0.0 0.44"/>
<surface id="052" type="z-cylinder" coeffs="0.0 0.0 0.515"/>
<cell id="051" universe="051" material="01" region="-051"/>
<cell id="052" universe="051" material="21" region="051 -052"/>
<cell id="053" universe="051" material="01" region="052"/>
<!-- Universe 055: Cluster tube -->
<surface id="155" type="z-cylinder" coeffs="0.0 0.0 0.55"/>
<surface id="156" type="z-cylinder" coeffs="0.0 0.0 0.63"/>
<cell id="155" universe="055" material="01" region="-155"/>
<cell id="156" universe="055" material="22" region="155 -156"/>
<cell id="157" universe="055" material="01" region="156"/>
<!-- Lattice 101: Fuel assembly pins -->
<cell id="100" universe="100" material="01"/>
<hex_lattice id="101" outer="100" n_rings="11"
center="0.0 0.0" pitch="1.265">
<universes>
011
011 011
011 011 011
011 011 011 011
011 011 011 011 011
011 011 011 011 011 011
011 011 011 011 011 011 011
011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 055 011 011 011 011 011
011 011 011 055 011 011 055 011 011 011
011 011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011 011
011 011 055 011 011 055 011 055 011 011
011 011 011 011 055 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 055 011 011 011
011 011 055 011 011 051 011 011 055 011 011
011 011 011 055 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 055 011 011 011 011
011 011 055 011 055 011 011 055 011 011
011 011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011 011
011 011 011 055 011 011 055 011 011 011
011 011 011 011 011 055 011 011 011 011 011
011 011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011 011
011 011 011 011 011 011 011 011
011 011 011 011 011 011 011
011 011 011 011 011 011
011 011 011 011 011
011 011 011 011
011 011 011
011 011
011
</universes>
</hex_lattice>
<cell id="102" universe="102" fill="101"/>
<lattice id="103" outer="100" dimension="1 1 294"
lower_left="-25.0 -25.0 -177.0" pitch="50.0 50.0 1.2">
<universes>
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102 102 102 102 102 102 102
102 102 102 102
</universes>
</lattice>
<!-- Universe 0 -->
<surface id="1001" type="plane" coeffs="1.0 0.0 0.0 11.8"
boundary="reflective"/>
<surface id="1002" type="plane" coeffs="0.5 0.866 0.0 11.8"
boundary="reflective"/>
<surface id="1003" type="plane" coeffs="-0.5 0.866 0.0 11.8"
boundary="reflective"/>
<surface id="1004" type="plane" coeffs="-1.0 0.0 0.0 11.8"
boundary="reflective"/>
<surface id="1005" type="plane" coeffs="-0.5 -0.866 0.0 11.8"
boundary="reflective"/>
<surface id="1006" type="plane" coeffs="0.5 -0.866 0.0 11.8"
boundary="reflective"/>
<surface id="1007" type="z-plane" coeffs="-177.4" boundary="vacuum"/>
<surface id="1008" type="z-plane" coeffs="177.4" boundary="vacuum"/>
<cell id="1001" universe="0" fill="103"
region="-1001 -1002 -1003 -1004 -1005 -1006 1007 -1008"/>
</geometry>

42
tests/regression_tests/lattice_mixed/materials.xml
View file

@ -1,42 +0,0 @@
<?xml version="1.0"?>
<materials>
<!-- Material 01: Water -->
<material id="01">
<density value="0.998" units="g/cc"/>
<nuclide name="H1" ao="2.0"/>
<nuclide name="O16" ao="1.0"/>
<nuclide name="B10" ao="0.00085"/>
</material>
<!-- Material 13: Fuel (3.30%) -->
<material id="13">
<density value="10.371" units="g/cc"/>
<nuclide name="U235" ao="0.033408078"/>
<nuclide name="U238" ao="0.966591922"/>
<nuclide name="O16" ao="2.0"/>
</material>
<!-- Material 21: Zircaloy -->
<material id="21">
<density value="6.44" units="g/cc" />
<nuclide name="Zr90" wo="0.501807"/>
<nuclide name="Zr91" wo="0.110712"/>
<nuclide name="Zr92" wo="0.170991"/>
<nuclide name="Zr94" wo="0.177057"/>
<nuclide name="Zr96" wo="0.029133"/>
</material>
<!-- Material 22: Steel -->
<material id="22">
<density value="7.9" units="g/cc"/>
<nuclide name="Mn55" wo="2"/>
<nuclide name="Ni58" wo="6.83"/>
<nuclide name="Ni60" wo="2.61"/>
<nuclide name="Cr53" wo="1.71"/>
<nuclide name="Fe54" wo="3.98257"/>
<nuclide name="Fe56" wo="63.03447"/>
<nuclide name="Fe57" wo="1.4763"/>
</material>
</materials>

32
tests/regression_tests/lattice_mixed/plots.xml
View file

@ -1,32 +0,0 @@
<?xml version="1.0"?>
<plots>
<plot id="1" type="slice" basis="xy" color_by="cell">
<filename>xy_cell</filename>
<origin>0 0 0</origin>
<width>30 30</width>
<pixels>500 500</pixels>
</plot>
<plot id="2" type="slice" basis="xy" color_by="material">
<filename>xy_material</filename>
<origin>0 0 0</origin>
<width>30 30</width>
<pixels>500 500</pixels>
</plot>
<plot id="3" type="slice" basis="yz" color_by="cell">
<filename>yz_cell</filename>
<origin>0 0 0</origin>
<width>50 400</width>
<pixels>500 4000</pixels>
</plot>
<plot id="4" type="slice" basis="yz" color_by="material">
<filename>yz_material</filename>
<origin>0 0 0</origin>
<width>5 5</width>
<pixels>500 500</pixels>
</plot>
</plots>

2
tests/regression_tests/lattice_mixed/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
9.877395E-01 1.744626E-02

14
tests/regression_tests/lattice_mixed/settings.xml
View file

@ -1,14 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<settings>
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>5</inactive>
<particles>500</particles>
<source>
<space type="box">
<parameters>-8.0 -8.0 -176
8.0 8.0 176</parameters>
</space>
</source>
</settings>

6
tests/regression_tests/lattice_mixed/test.py
View file

@ -1,6 +0,0 @@
from tests.testing_harness import TestHarness
def test_lattice_mixed():
harness = TestHarness('statepoint.10.h5')
harness.main()

0
tests/regression_tests/lattice_mixed/__init__.py → tests/regression_tests/lattice_rotated/__init__.py
View file

71
tests/regression_tests/lattice_rotated/inputs_true.dat Normal file
View file

@ -0,0 +1,71 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="-1" universe="1" />
<cell id="2" material="3" region="1" universe="1" />
<cell id="3" material="2" region="-2" universe="2" />
<cell id="4" material="3" region="2" universe="2" />
<cell id="5" material="3" universe="30" />
<cell fill="4" id="6" region="-3" rotation="0.0 0.0 45.0" translation="-4.0 0.0 0.0" universe="5" />
<cell fill="3" id="7" region="-4" rotation="0.0 0.0 30.0" translation="4.0 0.0 0.0" universe="5" />
<cell id="8" material="3" region="-5 3 4" universe="5" />
<hex_lattice id="3" n_rings="3">
<pitch>1.25</pitch>
<outer>30</outer>
<center>0.0 0.0</center>
<universes>
2
1 1
1 2 1
1 1
1 2 1
1 1
1 1 1
1 1
1</universes>
</hex_lattice>
<lattice id="4">
<pitch>1.25 1.25</pitch>
<outer>30</outer>
<dimension>4 4</dimension>
<lower_left>-2.5 -2.5</lower_left>
<universes>
2 2 2 2
1 1 1 1
1 1 1 1
1 1 1 1 </universes>
</lattice>
<surface coeffs="0.0 0.0 0.25" id="1" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.5" id="2" type="z-cylinder" />
<surface coeffs="-4.0 0.0 4.0" id="3" type="z-cylinder" />
<surface coeffs="4.0 0.0 4.0" id="4" type="z-cylinder" />
<surface boundary="vacuum" coeffs="0.0 0.0 8.0" id="5" type="z-cylinder" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="1">
<density units="g/cm3" value="10.0" />
<nuclide ao="1.0" name="U235" />
</material>
<material depletable="true" id="2">
<density units="g/cm3" value="10.0" />
<nuclide ao="1.0" name="U238" />
</material>
<material id="3">
<density units="g/cm3" value="1.0" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
<sab name="c_H_in_H2O" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>5</batches>
<inactive>0</inactive>
<source strength="1.0">
<space type="point">
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>

2
tests/regression_tests/lattice_rotated/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
4.538090E-01 2.452457E-02

85
tests/regression_tests/lattice_rotated/test.py Normal file
View file

@ -0,0 +1,85 @@
import numpy as np
import openmc
from tests.testing_harness import PyAPITestHarness
def rotated_lattice_model():
model = openmc.model.Model()
# Create some materials
fuel1 = openmc.Material()
fuel1.set_density('g/cm3', 10.0)
fuel1.add_nuclide('U235', 1.0)
fuel2 = openmc.Material()
fuel2.set_density('g/cm3', 10.0)
fuel2.add_nuclide('U238', 1.0)
water = openmc.Material()
water.set_density('g/cm3', 1.0)
water.add_nuclide('H1', 2.0)
water.add_nuclide('O16', 1.0)
water.add_s_alpha_beta('c_H_in_H2O')
model.materials.extend([fuel1, fuel2, water])
# Create universes for lattices
r_pin = openmc.ZCylinder(r=0.25)
fuel_cell = openmc.Cell(fill=fuel1, region=-r_pin)
water_cell = openmc.Cell(fill=water, region=+r_pin)
pin_universe = openmc.Universe(cells=(fuel_cell, water_cell))
r_big_pin = openmc.ZCylinder(r=0.5)
fuel2_cell = openmc.Cell(fill=fuel2, region=-r_big_pin)
water2_cell = openmc.Cell(fill=water, region=+r_big_pin)
big_pin_universe = openmc.Universe(cells=(fuel2_cell, water2_cell))
all_water_cell = openmc.Cell(fill=water)
outer_universe = openmc.Universe(30, cells=(all_water_cell,))
# Create hexagonal lattice
pitch = 1.25
hexlat = openmc.HexLattice()
hexlat.center = (0., 0.)
hexlat.pitch = [pitch]
hexlat.outer = outer_universe
outer_ring = [big_pin_universe] + [pin_universe]*11
middle_ring = [big_pin_universe] + [pin_universe]*5
inner_ring = [big_pin_universe]
hexlat.universes = [outer_ring, middle_ring, inner_ring]
# Create rectangular lattice
rectlat = openmc.RectLattice()
rectlat.center = (0., 0.)
rectlat.pitch = (pitch, pitch)
rectlat.lower_left = (-2*pitch, -2*pitch)
rectlat.outer = outer_universe
rectlat.universes = np.full((4, 4), pin_universe)
rectlat.universes[0] = big_pin_universe
# Create cell filled with translated/rotated rectangular lattice on left
left_cyl = openmc.ZCylinder(x0=-4.0, r=4.0)
left_cell = openmc.Cell(fill=rectlat, region=-left_cyl)
left_cell.translation = (-4.0, 0.0, 0.0)
left_cell.rotation = (0.0, 0.0, 45.0)
# Create cell filled with translated/rotated hexagonal lattice on right
right_cyl = openmc.ZCylinder(x0=4.0, r=4.0)
right_cell = openmc.Cell(fill=hexlat, region=-right_cyl)
right_cell.translation = (4.0, 0.0, 0.0)
right_cell.rotation = (0.0, 0.0, 30.0)
# Finish up with the geometry
outer_cyl = openmc.ZCylinder(r=8.0, boundary_type='vacuum')
main_cell = openmc.Cell(fill=water, region=-outer_cyl & +left_cyl & +right_cyl)
model.geometry = openmc.Geometry([main_cell, left_cell, right_cell])
model.settings.batches = 5
model.settings.inactive = 0
model.settings.particles = 1000
model.settings.source = openmc.Source(space=openmc.stats.Point())
model.settings.export_to_xml()
return model
def test():
model = rotated_lattice_model()
harness = PyAPITestHarness('statepoint.5.h5', model)
harness.main()

2
tests/unit_tests/test_filters.py
View file

@ -11,7 +11,7 @@ def box_model():
m.add_nuclide('U235', 1.0)
m.set_density('g/cm3', 1.0)
box = openmc.model.get_rectangular_prism(10., 10., boundary_type='vacuum')
box = openmc.model.rectangular_prism(10., 10., boundary_type='vacuum')
c = openmc.Cell(fill=m, region=box)
model.geometry.root_universe = openmc.Universe(cells=[c])

14
tests/unit_tests/test_geometry.py
View file

@ -204,9 +204,9 @@ def test_get_by_name():
def test_hex_prism():
hex_prism = openmc.model.get_hexagonal_prism(edge_length=5.0,
origin=(0.0, 0.0),
orientation='y')
hex_prism = openmc.model.hexagonal_prism(edge_length=5.0,
origin=(0.0, 0.0),
orientation='y')
# clear checks
assert (0.0, 0.0, 0.0) in hex_prism
assert (10.0, 10.0, 10.0) not in hex_prism
@ -214,10 +214,10 @@ def test_hex_prism():
assert (0.0, 5.01, 0.0) not in hex_prism
assert (0.0, 4.99, 0.0) in hex_prism
rounded_hex_prism = openmc.model.get_hexagonal_prism(edge_length=5.0,
origin=(0.0, 0.0),
orientation='y',
corner_radius=1.0)
rounded_hex_prism = openmc.model.hexagonal_prism(edge_length=5.0,
origin=(0.0, 0.0),
orientation='y',
corner_radius=1.0)
# clear checks
assert (0.0, 0.0, 0.0) in rounded_hex_prism