diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 23a0201ca9..0e1cd16f4e 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -312,21 +312,30 @@ a separate element with the tag name ````. This element has the following attributes/sub-elements: :type: - The type of structured mesh. The only valid option is "regular". + The type of structured mesh. This can be either "regular" or "rectilinear". :dimension: - The number of mesh cells in each direction. + The number of mesh cells in each direction. (For regular mesh only.) :lower_left: The lower-left corner of the structured mesh. If only two coordinates are - given, it is assumed that the mesh is an x-y mesh. + given, it is assumed that the mesh is an x-y mesh. (For regular mesh only.) :upper_right: The upper-right corner of the structured mesh. If only two coordinates are - given, it is assumed that the mesh is an x-y mesh. + given, it is assumed that the mesh is an x-y mesh. (For regular mesh only.) :width: - The width of mesh cells in each direction. + The width of mesh cells in each direction. (For regular mesh only.) + + :x_grid: + The mesh divisions along the x-axis. (For rectilinear mesh only.) + + :y_grid: + The mesh divisions along the y-axis. (For rectilinear mesh only.) + + :z_grid: + The mesh divisions along the z-axis. (For rectilinear mesh only.) .. note:: One of ```` or ```` must be specified, but not both diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 8cf65a997f..2c0ac8797d 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -124,7 +124,8 @@ Constructing Tallies openmc.ZernikeFilter openmc.ZernikeRadialFilter openmc.ParticleFilter - openmc.Mesh + openmc.RegularMesh + openmc.RectilinearMesh openmc.Trigger openmc.TallyDerivative openmc.Tally diff --git a/examples/jupyter/mdgxs-part-ii.ipynb b/examples/jupyter/mdgxs-part-ii.ipynb index 4c2cc40486..23ee4d2514 100644 --- a/examples/jupyter/mdgxs-part-ii.ipynb +++ b/examples/jupyter/mdgxs-part-ii.ipynb @@ -22,9 +22,7 @@ { "cell_type": "code", "execution_count": 1, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", @@ -82,9 +80,7 @@ { "cell_type": "code", "execution_count": 3, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Create a materials collection and export to XML\n", @@ -102,9 +98,7 @@ { "cell_type": "code", "execution_count": 4, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Create cylinders for the fuel and clad\n", @@ -167,9 +161,7 @@ { "cell_type": "code", "execution_count": 6, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Create a Universe to encapsulate a control rod guide tube\n", @@ -204,9 +196,7 @@ { "cell_type": "code", "execution_count": 7, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Create fuel assembly Lattice\n", @@ -225,9 +215,7 @@ { "cell_type": "code", "execution_count": 8, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Create array indices for guide tube locations in lattice\n", @@ -254,9 +242,7 @@ { "cell_type": "code", "execution_count": 9, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Create root Cell\n", @@ -280,9 +266,7 @@ { "cell_type": "code", "execution_count": 10, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Create Geometry and export to XML\n", @@ -300,9 +284,7 @@ { "cell_type": "code", "execution_count": 11, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# OpenMC simulation parameters\n", @@ -336,9 +318,7 @@ { "cell_type": "code", "execution_count": 12, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -383,9 +363,7 @@ { "cell_type": "code", "execution_count": 13, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Instantiate a 20-group EnergyGroups object\n", @@ -407,14 +385,11 @@ { "cell_type": "code", "execution_count": 14, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Instantiate a tally mesh \n", - "mesh = openmc.Mesh(mesh_id=1)\n", - "mesh.type = 'regular'\n", + "mesh = openmc.RegularMesh(mesh_id=1)\n", "mesh.dimension = [17, 17, 1]\n", "mesh.lower_left = [-10.71, -10.71, -10000.]\n", "mesh.width = [1.26, 1.26, 20000.]\n", @@ -464,9 +439,7 @@ { "cell_type": "code", "execution_count": 15, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -635,9 +608,7 @@ { "cell_type": "code", "execution_count": 16, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Load the last statepoint file\n", @@ -654,9 +625,7 @@ { "cell_type": "code", "execution_count": 17, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Initialize MGXS Library with OpenMC statepoint data\n", @@ -687,9 +656,7 @@ { "cell_type": "code", "execution_count": 18, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stderr", @@ -923,9 +890,7 @@ { "cell_type": "code", "execution_count": 19, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -1110,9 +1075,7 @@ { "cell_type": "code", "execution_count": 20, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stderr", @@ -1189,7 +1152,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python [default]", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -1203,9 +1166,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.2" + "version": "3.6.7" } }, "nbformat": 4, - "nbformat_minor": 0 + "nbformat_minor": 1 } diff --git a/examples/jupyter/mg-mode-part-i.ipynb b/examples/jupyter/mg-mode-part-i.ipynb index 8e52787879..23227b1ed3 100644 --- a/examples/jupyter/mg-mode-part-i.ipynb +++ b/examples/jupyter/mg-mode-part-i.ipynb @@ -86,9 +86,7 @@ { "cell_type": "code", "execution_count": 3, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# The scattering matrix is ordered with incoming groups as rows and outgoing groups as columns\n", @@ -158,9 +156,7 @@ { "cell_type": "code", "execution_count": 5, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# For every cross section data set in the library, assign an openmc.Macroscopic object to a material\n", @@ -207,9 +203,7 @@ { "cell_type": "code", "execution_count": 7, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Create the surface used for each pin\n", @@ -249,9 +243,7 @@ { "cell_type": "code", "execution_count": 8, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "lattices = {}\n", @@ -352,9 +344,7 @@ { "cell_type": "code", "execution_count": 9, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "lattices['Core'] = openmc.RectLattice(name='3x3 core lattice')\n", @@ -396,9 +386,7 @@ { "cell_type": "code", "execution_count": 10, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -426,9 +414,7 @@ { "cell_type": "code", "execution_count": 11, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Create Geometry and set root Universe\n", @@ -456,8 +442,7 @@ "tallies_file = openmc.Tallies()\n", "\n", "# Instantiate a tally Mesh\n", - "mesh = openmc.Mesh()\n", - "mesh.type = 'regular'\n", + "mesh = openmc.RegularMesh()\n", "mesh.dimension = [17 * 2, 17 * 2]\n", "mesh.lower_left = [-32.13, -10.71]\n", "mesh.upper_right = [+10.71, +32.13]\n", @@ -489,9 +474,7 @@ { "cell_type": "code", "execution_count": 13, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# OpenMC simulation parameters\n", @@ -533,9 +516,7 @@ { "cell_type": "code", "execution_count": 14, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -649,9 +630,7 @@ { "cell_type": "code", "execution_count": 15, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -714,7 +693,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.0" + "version": "3.6.7" } }, "nbformat": 4, diff --git a/examples/jupyter/mg-mode-part-ii.ipynb b/examples/jupyter/mg-mode-part-ii.ipynb index e2df6fd3db..85a8d73675 100644 --- a/examples/jupyter/mg-mode-part-ii.ipynb +++ b/examples/jupyter/mg-mode-part-ii.ipynb @@ -564,8 +564,7 @@ ], "source": [ "# Instantiate a tally Mesh\n", - "mesh = openmc.Mesh()\n", - "mesh.type = 'regular'\n", + "mesh = openmc.RegularMesh()\n", "mesh.dimension = [17, 17]\n", "mesh.lower_left = [-10.71, -10.71]\n", "mesh.upper_right = [+10.71, +10.71]\n", @@ -1518,9 +1517,9 @@ ], "metadata": { "kernelspec": { - "display_name": "openmc", + "display_name": "Python 3", "language": "python", - "name": "openmc" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -1532,7 +1531,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.6.7" } }, "nbformat": 4, diff --git a/examples/jupyter/mg-mode-part-iii.ipynb b/examples/jupyter/mg-mode-part-iii.ipynb index 7487b842c6..2cf9981889 100644 --- a/examples/jupyter/mg-mode-part-iii.ipynb +++ b/examples/jupyter/mg-mode-part-iii.ipynb @@ -545,8 +545,7 @@ "outputs": [], "source": [ "# Instantiate a tally Mesh\n", - "mesh = openmc.Mesh()\n", - "mesh.type = 'regular'\n", + "mesh = openmc.RegularMesh()\n", "mesh.dimension = [10, 10]\n", "mesh.lower_left = [0., 0.]\n", "mesh.upper_right = [length, length]\n", @@ -1417,9 +1416,9 @@ ], "metadata": { "kernelspec": { - "display_name": "openmc", + "display_name": "Python 3", "language": "python", - "name": "openmc" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -1431,7 +1430,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.6.7" } }, "nbformat": 4, diff --git a/examples/jupyter/mgxs-part-iii.ipynb b/examples/jupyter/mgxs-part-iii.ipynb index d3e90d0946..5ffb66d3a0 100644 --- a/examples/jupyter/mgxs-part-iii.ipynb +++ b/examples/jupyter/mgxs-part-iii.ipynb @@ -586,8 +586,7 @@ "outputs": [], "source": [ "# Instantiate a tally Mesh\n", - "mesh = openmc.Mesh(mesh_id=1)\n", - "mesh.type = 'regular'\n", + "mesh = openmc.RegularMesh(mesh_id=1)\n", "mesh.dimension = [17, 17]\n", "mesh.lower_left = [-10.71, -10.71]\n", "mesh.upper_right = [+10.71, +10.71]\n", @@ -1527,9 +1526,9 @@ "metadata": { "anaconda-cloud": {}, "kernelspec": { - "display_name": "openmc", + "display_name": "Python 3", "language": "python", - "name": "openmc" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -1541,7 +1540,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.6.7" } }, "nbformat": 4, diff --git a/examples/jupyter/pandas-dataframes.ipynb b/examples/jupyter/pandas-dataframes.ipynb index b8b10311c4..80a0b880f9 100644 --- a/examples/jupyter/pandas-dataframes.ipynb +++ b/examples/jupyter/pandas-dataframes.ipynb @@ -341,8 +341,7 @@ "outputs": [], "source": [ "# Instantiate a tally Mesh\n", - "mesh = openmc.Mesh(mesh_id=1)\n", - "mesh.type = 'regular'\n", + "mesh = openmc.RegularMesh(mesh_id=1)\n", "mesh.dimension = [17, 17]\n", "mesh.lower_left = [-10.71, -10.71]\n", "mesh.width = [1.26, 1.26]\n", @@ -1986,7 +1985,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.6" + "version": "3.6.7" } }, "nbformat": 4, diff --git a/examples/jupyter/post-processing.ipynb b/examples/jupyter/post-processing.ipynb index 1223f67f33..a6d134b54b 100644 --- a/examples/jupyter/post-processing.ipynb +++ b/examples/jupyter/post-processing.ipynb @@ -319,7 +319,7 @@ "outputs": [], "source": [ "# Create mesh which will be used for tally\n", - "mesh = openmc.Mesh()\n", + "mesh = openmc.RegularMesh()\n", "mesh.dimension = [100, 100]\n", "mesh.lower_left = [-0.63, -0.63]\n", "mesh.upper_right = [0.63, 0.63]\n", diff --git a/examples/jupyter/tally-arithmetic.ipynb b/examples/jupyter/tally-arithmetic.ipynb index 2bccf73f2a..35c096ba4d 100644 --- a/examples/jupyter/tally-arithmetic.ipynb +++ b/examples/jupyter/tally-arithmetic.ipynb @@ -284,9 +284,7 @@ { "cell_type": "code", "execution_count": 11, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Run openmc in plotting mode\n", @@ -296,9 +294,7 @@ { "cell_type": "code", "execution_count": 12, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -376,8 +372,7 @@ "tallies_file.append(tally)\n", "\n", "# Instantiate a tally mesh\n", - "mesh = openmc.Mesh(mesh_id=1)\n", - "mesh.type = 'regular'\n", + "mesh = openmc.RegularMesh(mesh_id=1)\n", "mesh.dimension = [1, 1, 1]\n", "mesh.lower_left = [-0.63, -0.63, -100.]\n", "mesh.width = [1.26, 1.26, 200.]\n", @@ -485,9 +480,7 @@ { "cell_type": "code", "execution_count": 20, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stderr", @@ -518,7 +511,6 @@ "cell_type": "code", "execution_count": 21, "metadata": { - "collapsed": false, "scrolled": true }, "outputs": [ @@ -674,9 +666,7 @@ { "cell_type": "code", "execution_count": 23, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -740,9 +730,7 @@ { "cell_type": "code", "execution_count": 24, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -807,9 +795,7 @@ { "cell_type": "code", "execution_count": 25, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -873,9 +859,7 @@ { "cell_type": "code", "execution_count": 26, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -939,9 +923,7 @@ { "cell_type": "code", "execution_count": 27, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -1004,9 +986,7 @@ { "cell_type": "code", "execution_count": 28, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -1066,9 +1046,7 @@ { "cell_type": "code", "execution_count": 29, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -1128,9 +1106,7 @@ { "cell_type": "code", "execution_count": 30, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -1210,9 +1186,7 @@ { "cell_type": "code", "execution_count": 32, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -1358,9 +1332,7 @@ { "cell_type": "code", "execution_count": 33, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -1390,9 +1362,7 @@ { "cell_type": "code", "execution_count": 34, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -1414,9 +1384,7 @@ { "cell_type": "code", "execution_count": 35, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -1445,9 +1413,7 @@ { "cell_type": "code", "execution_count": 36, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -1539,9 +1505,7 @@ { "cell_type": "code", "execution_count": 37, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -1709,7 +1673,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.0" + "version": "3.6.7" } }, "nbformat": 4, diff --git a/examples/python/lattice/nested/build-xml.py b/examples/python/lattice/nested/build-xml.py index 37483a9397..06641f5ac5 100644 --- a/examples/python/lattice/nested/build-xml.py +++ b/examples/python/lattice/nested/build-xml.py @@ -150,8 +150,7 @@ plot_file.export_to_xml() ############################################################################### # Instantiate a tally mesh -mesh = openmc.Mesh(mesh_id=1) -mesh.type = 'regular' +mesh = openmc.RegularMesh(mesh_id=1) mesh.dimension = [4, 4] mesh.lower_left = [-2, -2] mesh.width = [1, 1] @@ -164,6 +163,7 @@ tally = openmc.Tally(tally_id=1) tally.filters = [mesh_filter] tally.scores = ['total'] -# Instantiate a Tallies collection, register Tally/Mesh, and export to XML +# Instantiate a Tallies collection, register Tally/RegularMesh, and export to +# XML tallies_file = openmc.Tallies([tally]) tallies_file.export_to_xml() diff --git a/examples/python/lattice/simple/build-xml.py b/examples/python/lattice/simple/build-xml.py index 1cdfa24f3b..17544b87f5 100644 --- a/examples/python/lattice/simple/build-xml.py +++ b/examples/python/lattice/simple/build-xml.py @@ -143,8 +143,7 @@ plot_file.export_to_xml() ############################################################################### # Instantiate a tally mesh -mesh = openmc.Mesh(mesh_id=1) -mesh.type = 'regular' +mesh = openmc.RegularMesh(mesh_id=1) mesh.dimension = [4, 4] mesh.lower_left = [-2, -2] mesh.width = [1, 1] diff --git a/examples/python/pincell/build-xml.py b/examples/python/pincell/build-xml.py index c728d5f36c..072cc911f0 100644 --- a/examples/python/pincell/build-xml.py +++ b/examples/python/pincell/build-xml.py @@ -106,7 +106,7 @@ bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) settings_file.source = openmc.source.Source(space=uniform_dist) -entropy_mesh = openmc.Mesh() +entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] entropy_mesh.dimension = [10, 10, 1] @@ -119,8 +119,7 @@ settings_file.export_to_xml() ############################################################################### # Instantiate a tally mesh -mesh = openmc.Mesh() -mesh.type = 'regular' +mesh = openmc.RegularMesh() mesh.dimension = [100, 100, 1] mesh.lower_left = [-0.62992, -0.62992, -1.e50] mesh.upper_right = [0.62992, 0.62992, 1.e50] diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py index f9324a8cb8..d5a0e2234a 100644 --- a/examples/python/pincell_depletion/restart_depletion.py +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -47,7 +47,7 @@ bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) settings_file.source = openmc.source.Source(space=uniform_dist) -entropy_mesh = openmc.Mesh() +entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] entropy_mesh.dimension = [10, 10, 1] diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py index c1ad9d20aa..5a018cda10 100644 --- a/examples/python/pincell_depletion/run_depletion.py +++ b/examples/python/pincell_depletion/run_depletion.py @@ -119,7 +119,7 @@ bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) settings_file.source = openmc.source.Source(space=uniform_dist) -entropy_mesh = openmc.Mesh() +entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] entropy_mesh.dimension = [10, 10, 1] diff --git a/examples/python/pincell_multigroup/build-xml.py b/examples/python/pincell_multigroup/build-xml.py index 995cd87db5..bb9fc78606 100644 --- a/examples/python/pincell_multigroup/build-xml.py +++ b/examples/python/pincell_multigroup/build-xml.py @@ -155,8 +155,7 @@ settings_file.export_to_xml() ############################################################################### # Instantiate a tally mesh -mesh = openmc.Mesh(mesh_id=1) -mesh.type = 'regular' +mesh = openmc.RegularMesh(mesh_id=1) mesh.dimension = [100, 100, 1] mesh.lower_left = [-0.63, -0.63, -1.e50] mesh.upper_right = [0.63, 0.63, 1.e50] diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index f53bba224d..6d59826aaf 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -10,7 +10,7 @@ #include "hdf5.h" #include "pugixml.hpp" -#include "xtensor/xarray.hpp" +#include "xtensor/xtensor.hpp" #include "openmc/particle.h" #include "openmc/position.h" @@ -21,24 +21,21 @@ namespace openmc { // Global variables //============================================================================== -class RegularMesh; +class Mesh; namespace model { -extern std::vector> meshes; +extern std::vector> meshes; extern std::unordered_map mesh_map; } // namespace model -//============================================================================== -//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes -//============================================================================== - -class RegularMesh { +class Mesh +{ public: // Constructors - RegularMesh() = default; - RegularMesh(pugi::xml_node node); + Mesh() = default; + Mesh(pugi::xml_node node); // Methods @@ -47,39 +44,108 @@ public: //! \param[in] p Particle to check //! \param[out] bins Bins that were crossed //! \param[out] lengths Fraction of tracklength in each bin - void bins_crossed(const Particle* p, std::vector& bins, - std::vector& lengths) const; + virtual void bins_crossed(const Particle* p, std::vector& bins, + std::vector& lengths) const = 0; //! Determine which surface bins were crossed by a particle // //! \param[in] p Particle to check //! \param[out] bins Surface bins that were crossed - void surface_bins_crossed(const Particle* p, std::vector& bins) const; + virtual void + surface_bins_crossed(const Particle* p, std::vector& bins) const = 0; //! Get bin at a given position in space // //! \param[in] r Position to get bin for //! \return Mesh bin - int get_bin(Position r) const; + virtual int get_bin(Position r) const = 0; //! Get bin given mesh indices // //! \param[in] Array of mesh indices //! \return Mesh bin - int get_bin_from_indices(const int* ijk) const; + virtual int get_bin_from_indices(const int* ijk) const = 0; //! Get mesh indices given a position // //! \param[in] r Position to get indices for //! \param[out] ijk Array of mesh indices //! \param[out] in_mesh Whether position is in mesh - void get_indices(Position r, int* ijk, bool* in_mesh) const; + virtual void get_indices(Position r, int* ijk, bool* in_mesh) const = 0; //! Get mesh indices corresponding to a mesh bin // //! \param[in] bin Mesh bin //! \param[out] ijk Mesh indices - void get_indices_from_bin(int bin, int* ijk) const; + virtual void get_indices_from_bin(int bin, int* ijk) const = 0; + + //! Get the number of mesh cells. + virtual int n_bins() const = 0; + + //! Get the number of mesh cell surfaces. + virtual int n_surface_bins() const = 0; + + //! Find the mesh lines that intersect an axis-aligned slice plot + // + //! \param[in] plot_ll The lower-left coordinates of the slice plot. + //! \param[in] plot_ur The upper-right coordinates of the slice plot. + //! \return A pair of vectors indicating where the mesh lines lie along each + //! of the plot's axes. For example an xy-slice plot will get back a vector + //! of x-coordinates and another of y-coordinates. These vectors may be + //! empty for low-dimensional meshes. + virtual std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const = 0; + + //! Write mesh data to an HDF5 group + // + //! \param[in] group HDF5 group + virtual void to_hdf5(hid_t group) const = 0; + + // Data members + + int id_ {-1}; //!< User-specified ID + int n_dimension_; //!< Number of dimensions + xt::xtensor lower_left_; //!< Lower-left coordinates of mesh + xt::xtensor upper_right_; //!< Upper-right coordinates of mesh +}; + +//============================================================================== +//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes +//============================================================================== + +class RegularMesh : public Mesh +{ +public: + // Constructors + RegularMesh() = default; + RegularMesh(pugi::xml_node node); + + // Overriden methods + + void bins_crossed(const Particle* p, std::vector& bins, + std::vector& lengths) const override; + + void surface_bins_crossed(const Particle* p, std::vector& bins) + const override; + + int get_bin(Position r) const override; + + int get_bin_from_indices(const int* ijk) const override; + + void get_indices(Position r, int* ijk, bool* in_mesh) const override; + + void get_indices_from_bin(int bin, int* ijk) const override; + + int n_bins() const override; + + int n_surface_bins() const override; + + std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const override; + + void to_hdf5(hid_t group) const override; + + // New methods //! Check where a line segment intersects the mesh and if it intersects at all // @@ -89,26 +155,19 @@ public: //! \return Whether the line segment connecting r0 and r1 intersects mesh bool intersects(Position& r0, Position r1, int* ijk) const; - //! Write mesh data to an HDF5 group - // - //! \param[in] group HDF5 group - void to_hdf5(hid_t group) const; - //! Count number of bank sites in each mesh bin / energy bin // //! \param[in] bank Array of bank sites //! \param[out] Whether any bank sites are outside the mesh //! \return Array indicating number of sites in each mesh/energy bin - xt::xarray count_sites(const std::vector& bank, + xt::xtensor count_sites(const std::vector& bank, bool* outside) const; - int id_ {-1}; //!< User-specified ID - int n_dimension_; //!< Number of dimensions + // Data members + double volume_frac_; //!< Volume fraction of each mesh element - xt::xarray shape_; //!< Number of mesh elements in each dimension - xt::xarray lower_left_; //!< Lower-left coordinates of mesh - xt::xarray upper_right_; //!< Upper-right coordinates of mesh - xt::xarray width_; //!< Width of each mesh element + xt::xtensor shape_; //!< Number of mesh elements in each dimension + xt::xtensor width_; //!< Width of each mesh element private: bool intersects_1d(Position& r0, Position r1, int* ijk) const; @@ -116,6 +175,55 @@ private: bool intersects_3d(Position& r0, Position r1, int* ijk) const; }; +class RectilinearMesh : public Mesh +{ +public: + // Constructors + RectilinearMesh(pugi::xml_node node); + + // Overriden methods + + void bins_crossed(const Particle* p, std::vector& bins, + std::vector& lengths) const override; + + void surface_bins_crossed(const Particle* p, std::vector& bins) + const override; + + int get_bin(Position r) const override; + + int get_bin_from_indices(const int* ijk) const override; + + void get_indices(Position r, int* ijk, bool* in_mesh) const override; + + void get_indices_from_bin(int bin, int* ijk) const override; + + int n_bins() const override; + + int n_surface_bins() const override; + + std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const override; + + void to_hdf5(hid_t group) const override; + + // New methods + + //! Check where a line segment intersects the mesh and if it intersects at all + // + //! \param[in,out] r0 In: starting position, out: intersection point + //! \param[in] r1 Ending position + //! \param[out] ijk Indices of the mesh bin containing the intersection point + //! \return Whether the line segment connecting r0 and r1 intersects mesh + bool intersects(Position& r0, Position r1, int* ijk) const; + + // Data members + + xt::xtensor shape_; //!< Number of mesh elements in each dimension + +private: + std::vector> grid_; +}; + //============================================================================== // Non-member functions //============================================================================== diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 4b76c3dc59..77127db0a4 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -59,9 +59,6 @@ extern std::string path_source; extern std::string path_sourcepoint; //!< path to a source file extern std::string path_statepoint; //!< path to a statepoint file -extern int32_t index_entropy_mesh; //!< Index of entropy mesh in global mesh array -extern int32_t index_ufs_mesh; //!< Index of UFS mesh in global mesh array - extern "C" int32_t n_batches; //!< number of (inactive+active) batches extern "C" int32_t n_inactive; //!< number of inactive batches extern "C" int32_t gen_per_batch; //!< number of generations per batch diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index fef6ebbcd1..8efac4c7f8 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -4,6 +4,7 @@ #ifndef OPENMC_SIMULATION_H #define OPENMC_SIMULATION_H +#include "openmc/mesh.h" #include "openmc/particle.h" #include @@ -39,6 +40,9 @@ extern "C" int total_gen; //!< total number of generations simulated extern double total_weight; //!< Total source weight in a batch extern int64_t work_per_rank; //!< number of particles per MPI rank +extern const RegularMesh* entropy_mesh; +extern const RegularMesh* ufs_mesh; + extern std::vector k_generation; extern std::vector work_index; diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 05ec85d010..c0716bcef6 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -11,7 +11,7 @@ from . import _dll from .core import _FortranObjectWithID from .error import _error_handler from .material import Material -from .mesh import Mesh +from .mesh import RegularMesh __all__ = ['Filter', 'AzimuthalFilter', 'CellFilter', @@ -255,7 +255,7 @@ class MeshFilter(Filter): def mesh(self): index_mesh = c_int32() _dll.openmc_mesh_filter_get_mesh(self._index, index_mesh) - return Mesh(index=index_mesh.value) + return RegularMesh(index=index_mesh.value) @mesh.setter def mesh(self, mesh): @@ -274,7 +274,7 @@ class MeshSurfaceFilter(Filter): def mesh(self): index_mesh = c_int32() _dll.openmc_meshsurface_filter_get_mesh(self._index, index_mesh) - return Mesh(index=index_mesh.value) + return RegularMesh(index=index_mesh.value) @mesh.setter def mesh(self, mesh): diff --git a/openmc/capi/mesh.py b/openmc/capi/mesh.py index 70a4345b09..a1161f48f3 100644 --- a/openmc/capi/mesh.py +++ b/openmc/capi/mesh.py @@ -11,7 +11,7 @@ from .core import _FortranObjectWithID from .error import _error_handler from .material import Material -__all__ = ['Mesh', 'meshes'] +__all__ = ['RegularMesh', 'meshes'] # Mesh functions _dll.openmc_extend_meshes.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)] @@ -45,8 +45,8 @@ _dll.n_meshes.argtypes = [] _dll.n_meshes.restype = c_int -class Mesh(_FortranObjectWithID): - """Mesh stored internally. +class RegularMesh(_FortranObjectWithID): + """RegularMesh stored internally. This class exposes a mesh that is stored internally in the OpenMC library. To obtain a view of a mesh with a given ID, use the @@ -170,11 +170,11 @@ class _MeshMapping(Mapping): except (AllocationError, InvalidIDError) as e: # __contains__ expects a KeyError to work correctly raise KeyError(str(e)) - return Mesh(index=index.value) + return RegularMesh(index=index.value) def __iter__(self): for i in range(len(self)): - yield Mesh(index=i).id + yield RegularMesh(index=i).id def __len__(self): return _dll.n_meshes() diff --git a/openmc/cmfd.py b/openmc/cmfd.py index d1ea77fa06..244ceadcf0 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -2605,16 +2605,16 @@ class CMFDRun(object): def _create_cmfd_tally(self): """Creates all tallies in-memory that are used to solve CMFD problem""" # Create Mesh object based on CMFDMesh, stored internally - cmfd_mesh = openmc.capi.Mesh() - # Store id of Mesh object + cmfd_mesh = openmc.capi.RegularMesh() + # Store id of mesh object self._mesh_id = cmfd_mesh.id - # Set dimension and parameters of Mesh object + # Set dimension and parameters of mesh object cmfd_mesh.dimension = self._mesh.dimension cmfd_mesh.set_parameters(lower_left=self._mesh.lower_left, upper_right=self._mesh.upper_right, width=self._mesh.width) - # Create Mesh Filter object, stored internally + # Create mesh Filter object, stored internally mesh_filter = openmc.capi.MeshFilter() # Set mesh for Mesh Filter mesh_filter.mesh = cmfd_mesh diff --git a/openmc/filter.py b/openmc/filter.py index 49d3ff5de1..dbdaf2e19e 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -161,8 +161,8 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): Keyword arguments ----------------- meshes : dict - Dictionary mapping integer IDs to openmc.Mesh objects. Only used - for openmc.MeshFilter objects. + Dictionary mapping integer IDs to openmc.MeshBase objects. Only + used for openmc.MeshFilter objects. """ @@ -587,15 +587,15 @@ class MeshFilter(Filter): Parameters ---------- - mesh : openmc.Mesh - The Mesh object that events will be tallied onto + mesh : openmc.MeshBase + The mesh object that events will be tallied onto filter_id : int Unique identifier for the filter Attributes ---------- - mesh : openmc.Mesh - The Mesh object that events will be tallied onto + mesh : openmc.MeshBase + The mesh object that events will be tallied onto id : int Unique identifier for the filter bins : list of tuple @@ -610,6 +610,11 @@ class MeshFilter(Filter): self.mesh = mesh self.id = filter_id + def __hash__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id) + return hash(string) + def __repr__(self): string = type(self).__name__ + '\n' string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id) @@ -641,7 +646,7 @@ class MeshFilter(Filter): @mesh.setter def mesh(self, mesh): - cv.check_type('filter mesh', mesh, openmc.Mesh) + cv.check_type('filter mesh', mesh, openmc.MeshBase) self._mesh = mesh self.bins = list(mesh.indices) @@ -683,7 +688,7 @@ class MeshFilter(Filter): # Initialize dictionary to build Pandas Multi-index column filter_dict = {} - # Append Mesh ID as outermost index of multi-index + # Append mesh ID as outermost index of multi-index mesh_key = 'mesh {}'.format(self.mesh.id) # Find mesh dimensions - use 3D indices for simplicity @@ -745,17 +750,17 @@ class MeshSurfaceFilter(MeshFilter): Parameters ---------- - mesh : openmc.Mesh - The Mesh object that events will be tallied onto + mesh : openmc.MeshBase + The mesh object that events will be tallied onto filter_id : int Unique identifier for the filter Attributes ---------- bins : Integral - The Mesh ID - mesh : openmc.Mesh - The Mesh object that events will be tallied onto + The mesh ID + mesh : openmc.MeshBase + The mesh object that events will be tallied onto id : int Unique identifier for the filter bins : list of tuple @@ -770,11 +775,11 @@ class MeshSurfaceFilter(MeshFilter): @MeshFilter.mesh.setter def mesh(self, mesh): - cv.check_type('filter mesh', mesh, openmc.Mesh) + cv.check_type('filter mesh', mesh, openmc.MeshBase) self._mesh = mesh # Take the product of mesh indices and current names - n_dim = len(mesh.dimension) + n_dim = mesh.n_dimension self.bins = [mesh_tuple + (surf,) for mesh_tuple, surf in product(mesh.indices, _CURRENT_NAMES[:4*n_dim])] @@ -812,7 +817,7 @@ class MeshSurfaceFilter(MeshFilter): # Initialize dictionary to build Pandas Multi-index column filter_dict = {} - # Append Mesh ID as outermost index of multi-index + # Append mesh ID as outermost index of multi-index mesh_key = 'mesh {}'.format(self.mesh.id) # Find mesh dimensions - use 3D indices for simplicity diff --git a/openmc/mesh.py b/openmc/mesh.py index 2489d5faa7..4186b56a88 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1,7 +1,9 @@ +from abc import ABCMeta from collections.abc import Iterable from numbers import Real, Integral from xml.etree import ElementTree as ET import sys +import warnings import numpy as np @@ -11,8 +13,8 @@ from openmc._xml import get_text from openmc.mixin import EqualityMixin, IDManagerMixin -class Mesh(IDManagerMixin): - """A structured Cartesian mesh in one, two, or three dimensions +class MeshBase(IDManagerMixin, metaclass=ABCMeta): + """A mesh that partitions geometry for tallying purposes. Parameters ---------- @@ -27,21 +29,6 @@ class Mesh(IDManagerMixin): Unique identifier for the mesh name : str Name of the mesh - type : str - Type of the mesh - dimension : Iterable of int - The number of mesh cells in each direction. - lower_left : Iterable of float - The lower-left corner of the structured mesh. If only two coordinate are - given, it is assumed that the mesh is an x-y mesh. - upper_right : Iterable of float - The upper-right corner of the structrued mesh. If only two coordinate - are given, it is assumed that the mesh is an x-y mesh. - width : Iterable of float - The width of mesh cells in each direction. - indices : list of tuple - A list of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, - 1), ...] """ @@ -52,24 +39,95 @@ class Mesh(IDManagerMixin): # Initialize Mesh class attributes self.id = mesh_id self.name = name - self._type = 'regular' + + @property + def name(self): + return self._name + + @name.setter + def name(self, name): + if name is not None: + cv.check_type('name for mesh ID="{0}"'.format(self._id), + name, str) + self._name = name + else: + self._name = '' + + @classmethod + def from_hdf5(cls, group): + """Create mesh from HDF5 group + + Parameters + ---------- + group : h5py.Group + Group in HDF5 file + + Returns + ------- + openmc.MeshBase + Instance of a MeshBase subclass + + """ + + mesh_type = group['type'][()].decode() + if mesh_type == 'regular': + return RegularMesh.from_hdf5(group) + elif mesh_type == 'rectilinear': + return RectilinearMesh.from_hdf5(group) + else: + raise ValueError('Unrecognized mesh type: "' + mesh_type + '"') + + +class RegularMesh(MeshBase): + """A regular Cartesian mesh in one, two, or three dimensions + + Parameters + ---------- + mesh_id : int + Unique identifier for the mesh + name : str + Name of the mesh + + Attributes + ---------- + id : int + Unique identifier for the mesh + name : str + Name of the mesh + dimension : Iterable of int + The number of mesh cells in each direction. + n_dimension : int + Number of mesh dimensions. + lower_left : Iterable of float + The lower-left corner of the structured mesh. If only two coordinate are + given, it is assumed that the mesh is an x-y mesh. + upper_right : Iterable of float + The upper-right corner of the structrued mesh. If only two coordinate + are given, it is assumed that the mesh is an x-y mesh. + width : Iterable of float + The width of mesh cells in each direction. + indices : Iterable of tuple + An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), + (2, 1, 1), ...] + + """ + + def __init__(self, mesh_id=None, name=''): + super().__init__(mesh_id, name) + self._dimension = None self._lower_left = None self._upper_right = None self._width = None - @property - def name(self): - return self._name - - @property - def type(self): - return self._type - @property def dimension(self): return self._dimension + @property + def n_dimension(self): + return len(self._dimension) + @property def lower_left(self): return self._lower_left @@ -104,23 +162,6 @@ class Mesh(IDManagerMixin): nx, = self.dimension return ((x,) for x in range(1, nx + 1)) - @name.setter - def name(self, name): - if name is not None: - cv.check_type('name for mesh ID="{0}"'.format(self._id), - name, str) - self._name = name - else: - self._name = '' - - @type.setter - def type(self, meshtype): - cv.check_type('type for mesh ID="{0}"'.format(self._id), - meshtype, str) - cv.check_value('type for mesh ID="{0}"'.format(self._id), - meshtype, ['regular']) - self._type = meshtype - @dimension.setter def dimension(self, dimension): cv.check_type('mesh dimension', dimension, Iterable, Integral) @@ -146,7 +187,7 @@ class Mesh(IDManagerMixin): self._width = width def __repr__(self): - string = 'Mesh\n' + string = 'RegularMesh\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type) @@ -158,24 +199,10 @@ class Mesh(IDManagerMixin): @classmethod def from_hdf5(cls, group): - """Create mesh from HDF5 group - - Parameters - ---------- - group : h5py.Group - Group in HDF5 file - - Returns - ------- - openmc.Mesh - Mesh instance - - """ mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) # Read and assign mesh properties mesh = cls(mesh_id) - mesh.type = group['type'][()].decode() mesh.dimension = group['dimension'][()] mesh.lower_left = group['lower_left'][()] mesh.upper_right = group['upper_right'][()] @@ -201,8 +228,8 @@ class Mesh(IDManagerMixin): Returns ------- - openmc.Mesh - Mesh instance + openmc.RegularMesh + RegularMesh instance """ cv.check_type('rectangular lattice', lattice, openmc.RectLattice) @@ -229,7 +256,6 @@ class Mesh(IDManagerMixin): element = ET.Element("mesh") element.set("id", str(self._id)) - element.set("type", self._type) if self._dimension is not None: subelement = ET.SubElement(element, "dimension") @@ -416,3 +442,125 @@ class Mesh(IDManagerMixin): root_cell.fill = lattice return root_cell, cells + + +def Mesh(*args, **kwargs): + warnings.warn("Mesh has been renamed RegularMesh. Future versions of " + "OpenMC will not accept the name Mesh.") + return RegularMesh(*args, **kwargs) + + +class RectilinearMesh(MeshBase): + """A 3D rectilinear Cartesian mesh + + Parameters + ---------- + mesh_id : int + Unique identifier for the mesh + name : str + Name of the mesh + + Attributes + ---------- + id : int + Unique identifier for the mesh + name : str + Name of the mesh + n_dimension : int + Number of mesh dimensions (always 3 for a RectilinearMesh). + x_grid : Iterable of float + Mesh boundary points along the x-axis. + y_grid : Iterable of float + Mesh boundary points along the y-axis. + z_grid : Iterable of float + Mesh boundary points along the z-axis. + indices : Iterable of tuple + An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), + (2, 1, 1), ...] + + """ + + def __init__(self, mesh_id=None, name=''): + super().__init__(mesh_id, name) + + self._x_grid = None + self._y_grid = None + self._z_grid = None + + @property + def n_dimension(self): + return 3 + + @property + def x_grid(self): + return self._x_grid + + @property + def y_grid(self): + return self._y_grid + + @property + def z_grid(self): + return self._z_grid + + @property + def indices(self): + nx = len(self.x_grid) - 1 + ny = len(self.y_grid) - 1 + nz = len(self.z_grid) - 1 + return ((x, y, z) + for z in range(1, nz + 1) + for y in range(1, ny + 1) + for x in range(1, nx + 1)) + + @x_grid.setter + def x_grid(self, grid): + cv.check_type('mesh x_grid', grid, Iterable, Real) + self._x_grid = grid + + @y_grid.setter + def y_grid(self, grid): + cv.check_type('mesh y_grid', grid, Iterable, Real) + self._y_grid = grid + + @z_grid.setter + def z_grid(self, grid): + cv.check_type('mesh z_grid', grid, Iterable, Real) + self._z_grid = grid + + @classmethod + def from_hdf5(cls, group): + mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) + + # Read and assign mesh properties + mesh = cls(mesh_id) + mesh.x_grid = group['x_grid'][()] + mesh.y_grid = group['y_grid'][()] + mesh.z_grid = group['z_grid'][()] + + return mesh + + def to_xml_element(self): + """Return XML representation of the mesh + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing mesh data + + """ + + element = ET.Element("mesh") + element.set("id", str(self._id)) + element.set("type", "rectilinear") + + subelement = ET.SubElement(element, "x_grid") + subelement.text = ' '.join(map(str, self.x_grid)) + + subelement = ET.SubElement(element, "y_grid") + subelement.text = ' '.join(map(str, self.y_grid)) + + subelement = ET.SubElement(element, "z_grid") + subelement.text = ' '.join(map(str, self.z_grid)) + + return element diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 7b75d90fff..fbf08f85f7 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -51,7 +51,7 @@ class Library(object): The types of cross sections in the library (e.g., ['total', 'scatter']) domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization - domains : Iterable of openmc.Material, openmc.Cell, openmc.Universe or openmc.Mesh + domains : Iterable of openmc.Material, openmc.Cell, openmc.Universe or openmc.RegularMesh The spatial domain(s) for which MGXS in the Library are computed correction : {'P0', None} Apply the P0 correction to scattering matrices if set to 'P0' @@ -324,7 +324,7 @@ class Library(object): cv.check_type('domain', domains, Iterable, openmc.Universe) all_domains = self.geometry.get_all_universes().values() elif self.domain_type == 'mesh': - cv.check_type('domain', domains, Iterable, openmc.Mesh) + cv.check_type('domain', domains, Iterable, openmc.RegularMesh) # The mesh and geometry are independent, so set all_domains # to the input domains @@ -606,7 +606,7 @@ class Library(object): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh or Integral + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh or Integral The material, cell, or universe object of interest (or its ID) mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', chi', 'chi-prompt', 'inverse-velocity', 'prompt-nu-fission', 'prompt-nu-fission matrix', 'delayed-nu-fission', 'delayed-nu-fission matrix', 'chi-delayed', 'beta'} The type of multi-group cross section object to return @@ -631,7 +631,7 @@ class Library(object): elif self.domain_type == 'universe': cv.check_type('domain', domain, (openmc.Universe, Integral)) elif self.domain_type == 'mesh': - cv.check_type('domain', domain, (openmc.Mesh, Integral)) + cv.check_type('domain', domain, (openmc.RegularMesh, Integral)) # Check that requested domain is included in library if isinstance(domain, Integral): @@ -916,7 +916,7 @@ class Library(object): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization xsdata_name : str Name to apply to the "xsdata" entry produced by this method @@ -930,7 +930,7 @@ class Library(object): subdomain : iterable of int This parameter is not used unless using a mesh domain. In that case, the subdomain is an [i,j,k] index (1-based indexing) of the - mesh cell of interest in the openmc.Mesh object. Note: + mesh cell of interest in the openmc.RegularMesh object. Note: this parameter currently only supports subdomains within a mesh, and not the subdomains of a distribcell. @@ -952,7 +952,7 @@ class Library(object): """ cv.check_type('domain', domain, (openmc.Material, openmc.Cell, - openmc.Universe, openmc.Mesh)) + openmc.Universe, openmc.RegularMesh)) cv.check_type('xsdata_name', xsdata_name, str) cv.check_type('nuclide', nuclide, str) cv.check_value('xs_type', xs_type, ['macro', 'micro']) diff --git a/openmc/mgxs/mdgxs.py b/openmc/mgxs/mdgxs.py index 2d278d62d7..14bc524ee9 100644 --- a/openmc/mgxs/mdgxs.py +++ b/openmc/mgxs/mdgxs.py @@ -40,7 +40,7 @@ class MDGXS(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -68,7 +68,7 @@ class MDGXS(MGXS): Reaction type (e.g., 'chi-delayed', 'beta', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -238,7 +238,7 @@ class MDGXS(MGXS): mdgxs_type : {'delayed-nu-fission', 'chi-delayed', 'beta', 'decay-rate', 'delayed-nu-fission matrix'} The type of multi-delayed-group cross section object to return domain : openmc.Material or openmc.Cell or openmc.Universe or - openmc.Mesh + openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -917,7 +917,7 @@ class ChiDelayed(MDGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -945,7 +945,7 @@ class ChiDelayed(MDGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -1429,7 +1429,7 @@ class DelayedNuFissionXS(MDGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -1457,7 +1457,7 @@ class DelayedNuFissionXS(MDGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -1563,7 +1563,7 @@ class Beta(MDGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -1591,7 +1591,7 @@ class Beta(MDGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -1746,7 +1746,7 @@ class DecayRate(MDGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -1774,7 +1774,7 @@ class DecayRate(MDGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -1921,7 +1921,7 @@ class MatrixMDGXS(MDGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -1949,7 +1949,7 @@ class MatrixMDGXS(MDGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -2512,7 +2512,7 @@ class DelayedNuFissionMatrixXS(MatrixMDGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -2540,7 +2540,7 @@ class DelayedNuFissionMatrixXS(MatrixMDGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 01a03f6506..d3d39e0f67 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -56,7 +56,7 @@ _DOMAIN_TO_FILTER = {'cell': openmc.CellFilter, _DOMAINS = (openmc.Cell, openmc.Universe, openmc.Material, - openmc.Mesh) + openmc.RegularMesh) # Supported ScatterMatrixXS angular distribution types MU_TREATMENTS = ('legendre', 'histogram') @@ -124,7 +124,7 @@ class MGXS(metaclass=ABCMeta): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -150,7 +150,7 @@ class MGXS(metaclass=ABCMeta): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -606,7 +606,7 @@ class MGXS(metaclass=ABCMeta): self._domain_type = 'cell' elif isinstance(domain, openmc.Universe): self._domain_type = 'universe' - elif isinstance(domain, openmc.Mesh): + elif isinstance(domain, openmc.RegularMesh): self._domain_type = 'mesh' @domain_type.setter @@ -675,7 +675,7 @@ class MGXS(metaclass=ABCMeta): ---------- mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', 'chi', 'chi-prompt', 'inverse-velocity', 'prompt-nu-fission', 'prompt-nu-fission matrix'} The type of multi-group cross section object to return - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -1985,7 +1985,7 @@ class MatrixMGXS(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -2011,7 +2011,7 @@ class MatrixMGXS(MGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -2478,7 +2478,7 @@ class TotalXS(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -2504,7 +2504,7 @@ class TotalXS(MGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -2609,7 +2609,7 @@ class TransportXS(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -2640,7 +2640,7 @@ class TransportXS(MGXS): If True, the cross section data will include neutron multiplication by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -2834,7 +2834,7 @@ class AbsorptionXS(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -2860,7 +2860,7 @@ class AbsorptionXS(MGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -2961,7 +2961,7 @@ class CaptureXS(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -2987,7 +2987,7 @@ class CaptureXS(MGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -3103,7 +3103,7 @@ class FissionXS(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -3140,7 +3140,7 @@ class FissionXS(MGXS): If true, computes cross sections which only includes prompt neutrons by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -3282,7 +3282,7 @@ class KappaFissionXS(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -3308,7 +3308,7 @@ class KappaFissionXS(MGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -3408,7 +3408,7 @@ class ScatterXS(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -3439,7 +3439,7 @@ class ScatterXS(MGXS): If True, the cross section data will include neutron multiplication by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -3602,7 +3602,7 @@ class ScatterMatrixXS(MatrixMGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -3656,7 +3656,7 @@ class ScatterMatrixXS(MatrixMGXS): If True, the cross section data will include neutron multiplication by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -4715,7 +4715,7 @@ class MultiplicityMatrixXS(MatrixMGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -4741,7 +4741,7 @@ class MultiplicityMatrixXS(MatrixMGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -4880,7 +4880,7 @@ class ScatterProbabilityMatrix(MatrixMGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -4906,7 +4906,7 @@ class ScatterProbabilityMatrix(MatrixMGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -5042,7 +5042,7 @@ class NuFissionMatrixXS(MatrixMGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -5073,7 +5073,7 @@ class NuFissionMatrixXS(MatrixMGXS): If true, computes cross sections which only includes prompt neutrons by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -5198,7 +5198,7 @@ class Chi(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -5229,7 +5229,7 @@ class Chi(MGXS): If true, computes cross sections which only includes prompt neutrons by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization @@ -5780,7 +5780,7 @@ class InverseVelocity(MGXS): Parameters ---------- - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -5806,7 +5806,7 @@ class InverseVelocity(MGXS): Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain - domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization diff --git a/openmc/plots.py b/openmc/plots.py index 00d96fd7ee..5a52822989 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -214,8 +214,8 @@ class Plot(IDManagerMixin): level : int Universe depth to plot at meshlines : dict - Dictionary defining type, id, linewidth and color of a regular mesh - to be plotted on top of a plot + Dictionary defining type, id, linewidth and color of a mesh to be + plotted on top of a plot """ diff --git a/openmc/settings.py b/openmc/settings.py index 4fe87d908c..784c2ac6e6 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -9,7 +9,7 @@ import numpy as np from openmc._xml import clean_indentation, get_text import openmc.checkvalue as cv -from openmc import VolumeCalculation, Source, Mesh +from openmc import VolumeCalculation, Source, RegularMesh _RUN_MODES = ['eigenvalue', 'fixed source', 'plot', 'volume', 'particle restart'] _RES_SCAT_METHODS = ['dbrc', 'rvs'] @@ -45,7 +45,7 @@ class Settings(object): secondary bremsstrahlung photons ('ttb'). energy_mode : {'continuous-energy', 'multi-group'} Set whether the calculation should be continuous-energy or multi-group. - entropy_mesh : openmc.Mesh + entropy_mesh : openmc.RegularMesh Mesh to be used to calculate Shannon entropy. If the mesh dimensions are not specified. OpenMC assigns a mesh such that 20 source sites per mesh cell are to be expected on average. @@ -145,7 +145,7 @@ class Settings(object): Maximum number of batches simulated. If this is set, the number of batches specified via ``batches`` is interpreted as the minimum number of batches - ufs_mesh : openmc.Mesh + ufs_mesh : openmc.RegularMesh Mesh to be used for redistributing source sites via the uniform fision site (UFS) method. verbosity : int @@ -550,7 +550,7 @@ class Settings(object): @entropy_mesh.setter def entropy_mesh(self, entropy): - cv.check_type('entropy mesh', entropy, Mesh) + cv.check_type('entropy mesh', entropy, RegularMesh) if entropy.dimension: cv.check_length('entropy mesh dimension', entropy.dimension, 3) cv.check_length('entropy mesh lower-left corner', entropy.lower_left, 3) @@ -640,7 +640,7 @@ class Settings(object): @ufs_mesh.setter def ufs_mesh(self, ufs_mesh): - cv.check_type('UFS mesh', ufs_mesh, Mesh) + cv.check_type('UFS mesh', ufs_mesh, RegularMesh) cv.check_length('UFS mesh dimension', ufs_mesh.dimension, 3) cv.check_length('UFS mesh lower-left corner', ufs_mesh.lower_left, 3) cv.check_length('UFS mesh upper-right corner', ufs_mesh.upper_right, 3) @@ -1061,7 +1061,7 @@ class Settings(object): path = "./mesh[@id='{}']".format(int(text)) elem = root.find(path) if elem is not None: - self.entropy_mesh = Mesh.from_xml_element(elem) + self.entropy_mesh = RegularMesh.from_xml_element(elem) def _trigger_from_xml_element(self, root): elem = root.find('trigger') @@ -1126,7 +1126,7 @@ class Settings(object): path = "./mesh[@id='{}']".format(int(text)) elem = root.find(path) if elem is not None: - self.ufs_mesh = Mesh.from_xml_element(elem) + self.ufs_mesh = RegularMesh.from_xml_element(elem) def _resonance_scattering_from_xml_element(self, root): elem = root.find('resonance_scattering') diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 3dc6535841..c9a8f9f9a8 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -72,7 +72,7 @@ class StatePoint(object): k_generation : numpy.ndarray Estimate of k-effective for each batch/generation meshes : dict - Dictionary whose keys are mesh IDs and whose values are Mesh objects + Dictionary whose keys are mesh IDs and whose values are MeshBase objects n_batches : int Number of batches n_inactive : int @@ -292,9 +292,9 @@ class StatePoint(object): if not self._meshes_read: mesh_group = self._f['tallies/meshes'] - # Iterate over all Meshes + # Iterate over all meshes for group in mesh_group.values(): - mesh = openmc.Mesh.from_hdf5(group) + mesh = openmc.MeshBase.from_hdf5(group) self._meshes[mesh.id] = mesh self._meshes_read = True diff --git a/openmc/trigger.py b/openmc/trigger.py index 71f9c0b92b..98557aab58 100644 --- a/openmc/trigger.py +++ b/openmc/trigger.py @@ -31,7 +31,6 @@ class Trigger(object): """ def __init__(self, trigger_type, threshold): - # Initialize Mesh class attributes self.trigger_type = trigger_type self.threshold = threshold self._scores = [] diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 04fdabab1b..a0b7f82c5f 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -538,13 +538,10 @@ int openmc_get_keff(double* k_combined) void shannon_entropy() { - // Get pointer to entropy mesh - auto& m = model::meshes[settings::index_entropy_mesh]; - // Get source weight in each mesh bin bool sites_outside; - xt::xtensor p = m->count_sites(simulation::fission_bank, - &sites_outside); + xt::xtensor p = simulation::entropy_mesh->count_sites( + simulation::fission_bank, &sites_outside); // display warning message if there were sites outside entropy box if (sites_outside) { @@ -570,21 +567,19 @@ void shannon_entropy() void ufs_count_sites() { - auto &m = model::meshes[settings::index_ufs_mesh]; - if (simulation::current_batch == 1 && simulation::current_gen == 1) { // On the first generation, just assume that the source is already evenly // distributed so that effectively the production of fission sites is not // biased auto s = xt::view(simulation::source_frac, xt::all()); - s = m->volume_frac_; + s = simulation::ufs_mesh->volume_frac_; } else { // count number of source sites in each ufs mesh cell bool sites_outside; - simulation::source_frac = m->count_sites(simulation::source_bank, - &sites_outside); + simulation::source_frac = simulation::ufs_mesh->count_sites( + simulation::source_bank, &sites_outside); // Check for sites outside of the mesh if (mpi::master && sites_outside) { @@ -593,7 +588,7 @@ void ufs_count_sites() #ifdef OPENMC_MPI // Send source fraction to all processors - int n_bins = xt::prod(m->shape_)(); + int n_bins = xt::prod(simulation::ufs_mesh->shape_)(); MPI_Bcast(simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm); #endif @@ -611,17 +606,16 @@ void ufs_count_sites() double ufs_get_weight(const Particle* p) { - auto& m = model::meshes[settings::index_ufs_mesh]; - // Determine indices on ufs mesh for current location - int mesh_bin = m->get_bin(p->r()); + int mesh_bin = simulation::ufs_mesh->get_bin(p->r()); if (mesh_bin < 0) { p->write_restart(); fatal_error("Source site outside UFS mesh!"); } if (simulation::source_frac(mesh_bin) != 0.0) { - return m->volume_frac_ / simulation::source_frac(mesh_bin); + return simulation::ufs_mesh->volume_frac_ + / simulation::source_frac(mesh_bin); } else { return 1.0; } diff --git a/src/finalize.cpp b/src/finalize.cpp index 34c4b9781b..fa6b7628cc 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -71,8 +71,6 @@ int openmc_finalize() settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0}; settings::entropy_on = false; settings::gen_per_batch = 1; - settings::index_entropy_mesh = -1; - settings::index_ufs_mesh = -1; settings::legendre_to_tabular = true; settings::legendre_to_tabular_points = -1; settings::n_particles = -1; @@ -114,6 +112,9 @@ int openmc_finalize() simulation::satisfy_triggers = false; simulation::total_gen = 0; + simulation::entropy_mesh = nullptr; + simulation::ufs_mesh = nullptr; + data::energy_max = {INFTY, INFTY}; data::energy_min = {0.0, 0.0}; model::root_universe = -1; diff --git a/src/mesh.cpp b/src/mesh.cpp index ddfa272a6d..315381706e 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -31,7 +31,7 @@ namespace openmc { namespace model { -std::vector> meshes; +std::vector> meshes; std::unordered_map mesh_map; } // namespace model @@ -63,10 +63,10 @@ inline bool check_intersection_point(double x1, double x0, double y1, } //============================================================================== -// RegularMesh implementation +// Mesh implementation //============================================================================== -RegularMesh::RegularMesh(pugi::xml_node node) +Mesh::Mesh(pugi::xml_node node) { // Copy mesh id if (check_for_node(node, "id")) { @@ -78,17 +78,15 @@ RegularMesh::RegularMesh(pugi::xml_node node) std::to_string(id_)); } } +} - // Read mesh type - if (check_for_node(node, "type")) { - auto temp = get_node_value(node, "type", true, true); - if (temp == "regular") { - // TODO: move elsewhere - } else { - fatal_error("Invalid mesh type: " + temp); - } - } +//============================================================================== +// RegularMesh implementation +//============================================================================== +RegularMesh::RegularMesh(pugi::xml_node node) + : Mesh {node} +{ // Determine number of dimensions for mesh if (check_for_node(node, "dimension")) { shape_ = get_node_xarray(node, "dimension"); @@ -230,6 +228,18 @@ void RegularMesh::get_indices_from_bin(int bin, int* ijk) const } } +int RegularMesh::n_bins() const +{ + int n_bins = 1; + for (auto dim : shape_) n_bins *= dim; + return n_bins; +} + +int RegularMesh::n_surface_bins() const +{ + return 4 * n_dimension_ * n_bins(); +} + bool RegularMesh::intersects(Position& r0, Position r1, int* ijk) const { switch(n_dimension_) { @@ -719,6 +729,42 @@ void RegularMesh::surface_bins_crossed(const Particle* p, } } +std::pair, std::vector> +RegularMesh::plot(Position plot_ll, Position plot_ur) const +{ + // Figure out which axes lie in the plane of the plot. + std::array axes {-1, -1}; + if (plot_ur.z == plot_ll.z) { + axes[0] = 0; + if (n_dimension_ > 1) axes[1] = 1; + } else if (plot_ur.y == plot_ll.y) { + axes[0] = 0; + if (n_dimension_ > 2) axes[1] = 2; + } else if (plot_ur.x == plot_ll.x) { + if (n_dimension_ > 1) axes[0] = 1; + if (n_dimension_ > 2) axes[1] = 2; + } else { + fatal_error("Can only plot mesh lines on an axis-aligned plot"); + } + + // Get the coordinates of the mesh lines along both of the axes. + std::array, 2> axis_lines; + for (int i_ax = 0; i_ax < 2; ++i_ax) { + int axis = axes[i_ax]; + if (axis == -1) continue; + auto& lines {axis_lines[i_ax]}; + + double coord = lower_left_[axis]; + for (int i = 0; i < shape_[axis] + 1; ++i) { + if (coord >= plot_ll[axis] && coord <= plot_ur[axis]) + lines.push_back(coord); + coord += width_[axis]; + } + } + + return {axis_lines[0], axis_lines[1]}; +} + void RegularMesh::to_hdf5(hid_t group) const { hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); @@ -732,7 +778,7 @@ void RegularMesh::to_hdf5(hid_t group) const close_group(mesh_group); } -xt::xarray +xt::xtensor RegularMesh::count_sites(const std::vector& bank, bool* outside) const { @@ -783,6 +829,529 @@ RegularMesh::count_sites(const std::vector& bank, return counts; } +//============================================================================== +// RectilinearMesh implementation +//============================================================================== + +RectilinearMesh::RectilinearMesh(pugi::xml_node node) + : Mesh {node} +{ + n_dimension_ = 3; + + grid_.resize(3); + grid_[0] = get_node_array(node, "x_grid"); + grid_[1] = get_node_array(node, "y_grid"); + grid_[2] = get_node_array(node, "z_grid"); + + shape_ = {static_cast(grid_[0].size()) - 1, + static_cast(grid_[1].size()) - 1, + static_cast(grid_[2].size()) - 1}; + + for (const auto& g : grid_) { + if (g.size() < 2) fatal_error("x-, y-, and z- grids for rectilinear meshes " + "must each have at least 2 points"); + for (int i = 1; i < g.size(); ++i) { + if (g[i] <= g[i-1]) fatal_error("Values in for x-, y-, and z- grids for " + "rectilinear meshes must be sorted and unique."); + } + } + + lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()}; + upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()}; +} + +void RectilinearMesh::bins_crossed(const Particle* p, std::vector& bins, + std::vector& lengths) const +{ + // ======================================================================== + // Determine where the track intersects the mesh and if it intersects at all. + + // Copy the starting and ending coordinates of the particle. + Position last_r {p->r_last_}; + Position r {p->r()}; + Direction u {p->u()}; + + // Compute the length of the entire track. + double total_distance = (r - last_r).norm(); + + // While determining if this track intersects the mesh, offset the starting + // and ending coords by a bit. This avoid finite-precision errors that can + // occur when the mesh surfaces coincide with lattice or geometric surfaces. + Position r0 = last_r + TINY_BIT*u; + Position r1 = r - TINY_BIT*u; + + // Determine the mesh indices for the starting and ending coords. + int ijk0[3], ijk1[3]; + bool start_in_mesh; + get_indices(r0, ijk0, &start_in_mesh); + bool end_in_mesh; + get_indices(r1, ijk1, &end_in_mesh); + + // Reset coordinates and check for a mesh intersection if necessary. + if (start_in_mesh) { + // The initial coords lie in the mesh, use those coords for tallying. + r0 = last_r; + } else { + // The initial coords do not lie in the mesh. Check to see if the particle + // eventually intersects the mesh and compute the relevant coords and + // indices. + if (!intersects(r0, r1, ijk0)) return; + } + r1 = r; + + // ======================================================================== + // Find which mesh cells are traversed and the length of each traversal. + + while (true) { + if (std::equal(ijk0, ijk0+3, ijk1)) { + // The track ends in this cell. Use the particle end location rather + // than the mesh surface and stop iterating. + double distance = (r1 - r0).norm(); + bins.push_back(get_bin_from_indices(ijk0)); + lengths.push_back(distance / total_distance); + break; + } + + // The track exits this cell. Determine the distance to each mesh surface. + double d[3]; + for (int k = 0; k < 3; ++k) { + if (std::fabs(u[k]) < FP_PRECISION) { + d[k] = INFTY; + } else if (u[k] > 0) { + double xyz_cross = grid_[k][ijk0[k]]; + d[k] = (xyz_cross - r0[k]) / u[k]; + } else { + double xyz_cross = grid_[k][ijk0[k] - 1]; + d[k] = (xyz_cross - r0[k]) / u[k]; + } + } + + // Pick the closest mesh surface and append this traversal to the output. + auto j = std::min_element(d, d+3) - d; + double distance = d[j]; + bins.push_back(get_bin_from_indices(ijk0)); + lengths.push_back(distance / total_distance); + + // Translate to the oncoming mesh surface. + r0 += distance * u; + + // Increment the indices into the next mesh cell. + if (u[j] > 0.0) { + ++ijk0[j]; + } else { + --ijk0[j]; + } + + // If the next indices are invalid, then the track has left the mesh and + // we are done. + bool in_mesh = true; + for (int i = 0; i < 3; ++i) { + if (ijk0[i] < 1 || ijk0[i] > shape_[i]) { + in_mesh = false; + break; + } + } + if (!in_mesh) break; + } +} + +void RectilinearMesh::surface_bins_crossed(const Particle* p, + std::vector& bins) const +{ + // ======================================================================== + // Determine if the track intersects the tally mesh. + + // Copy the starting and ending coordinates of the particle. + Position r0 {p->r_last_current_}; + Position r1 {p->r()}; + Direction u {p->u()}; + + // Determine indices for starting and ending location. + int ijk0[3], ijk1[3]; + bool start_in_mesh; + get_indices(r0, ijk0, &start_in_mesh); + bool end_in_mesh; + get_indices(r1, ijk1, &end_in_mesh); + + // If the starting coordinates do not lie in the mesh, compute the coords and + // mesh indices of the first intersection, and add the bin for this first + // intersection. Return if the particle does not intersect the mesh at all. + if (!start_in_mesh) { + // Compute the incoming intersection coordinates and indices. + if (!intersects(r0, r1, ijk0)) return; + + // Determine which surface the particle entered. + double min_dist = INFTY; + int i_surf; + for (int i = 0; i < 3; ++i) { + if (u[i] > 0.0 && ijk0[i] == 1) { + double d = std::abs(r0[i] - grid_[i][0]); + if (d < min_dist) { + min_dist = d; + i_surf = 4*i + 2; + } + } else if (u[i] < 0.0 && ijk0[i] == shape_[i]) { + double d = std::abs(r0[i] - grid_[i][shape_[i]]); + if (d < min_dist) { + min_dist = d; + i_surf = 4*i + 4; + } + } // u[i] == 0 intentionally skipped + } + + // Add the incoming current bin. + int i_mesh = get_bin_from_indices(ijk0); + int i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); + } + + // If the ending coordinates do not lie in the mesh, compute the coords and + // mesh indices of the last intersection, and add the bin for this last + // intersection. + if (!end_in_mesh) { + // Compute the outgoing intersection coordinates and indices. + intersects(r1, r0, ijk1); + + // Determine which surface the particle exited. + double min_dist = INFTY; + int i_surf; + for (int i = 0; i < 3; ++i) { + if (u[i] > 0.0 && ijk1[i] == shape_[i]) { + double d = std::abs(r1[i] - grid_[i][shape_[i]]); + if (d < min_dist) { + min_dist = d; + i_surf = 4*i + 3; + } + } else if (u[i] < 0.0 && ijk1[i] == 1) { + double d = std::abs(r1[i] - grid_[i][0]); + if (d < min_dist) { + min_dist = d; + i_surf = 4*i + 1; + } + } // u[i] == 0 intentionally skipped + } + + // Add the outgoing current bin. + int i_mesh = get_bin_from_indices(ijk1); + int i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); + } + + // ======================================================================== + // Find which mesh surfaces are crossed. + + // Calculate number of surface crossings + int n_cross = 0; + for (int i = 0; i < 3; ++i) n_cross += std::abs(ijk1[i] - ijk0[i]); + if (n_cross == 0) return; + + // Bounding coordinates + Position xyz_cross; + for (int i = 0; i < 3; ++i) { + if (u[i] > 0.0) { + xyz_cross[i] = grid_[i][ijk0[i]]; + } else { + xyz_cross[i] = grid_[i][ijk0[i] - 1]; + } + } + + for (int j = 0; j < n_cross; ++j) { + // Set the distances to infinity + Position d {INFTY, INFTY, INFTY}; + + // Determine closest bounding surface. We need to treat + // special case where the cosine of the angle is zero since this would + // result in a divide-by-zero. + double distance = INFTY; + for (int i = 0; i < 3; ++i) { + if (u[i] == 0) { + d[i] = INFTY; + } else { + d[i] = (xyz_cross[i] - r0[i])/u[i]; + } + distance = std::min(distance, d[i]); + } + + // Loop over the dimensions + for (int i = 0; i < 3; ++i) { + // Check whether distance is the shortest distance + if (distance == d[i]) { + + // Check whether particle is moving in positive i direction + if (u[i] > 0) { + + // Outward current on i max surface + int i_surf = 4*i + 3; + int i_mesh = get_bin_from_indices(ijk0); + int i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); + + // Advance position + ++ijk0[i]; + xyz_cross[i] = grid_[i][ijk0[i]]; + + // Inward current on i min surface + i_surf = 4*i + 2; + i_mesh = get_bin_from_indices(ijk0); + i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); + + } else { + // The particle is moving in the negative i direction + + // Outward current on i min surface + int i_surf = 4*i + 1; + int i_mesh = get_bin_from_indices(ijk0); + int i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); + + // Advance position + --ijk0[i]; + xyz_cross[i] = grid_[i][ijk0[i] - 1]; + + // Inward current on i min surface + i_surf = 4*i + 4; + i_mesh = get_bin_from_indices(ijk0); + i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); + } + } + } + + // Calculate new coordinates + r0 += distance * u; + } +} + +int RectilinearMesh::get_bin(Position r) const +{ + // Determine indices + int ijk[3]; + bool in_mesh; + get_indices(r, ijk, &in_mesh); + if (!in_mesh) return -1; + + // Convert indices to bin + return get_bin_from_indices(ijk); +} + +int RectilinearMesh::get_bin_from_indices(const int* ijk) const +{ + return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0] - 1; +} + +void RectilinearMesh::get_indices(Position r, int* ijk, bool* in_mesh) const +{ + *in_mesh = true; + + for (int i = 0; i < 3; ++i) { + if (r[i] < grid_[i].front() || r[i] > grid_[i].back()) { + ijk[i] = -1; + *in_mesh = false; + } else { + ijk[i] = lower_bound_index(grid_[i].begin(), grid_[i].end(), r[i]) + 1; + } + } +} + +void RectilinearMesh::get_indices_from_bin(int bin, int* ijk) const +{ + ijk[0] = bin % shape_[0] + 1; + ijk[1] = (bin % (shape_[0] * shape_[1])) / shape_[0] + 1; + ijk[2] = bin / (shape_[0] * shape_[1]) + 1; +} + +int RectilinearMesh::n_bins() const +{ + int n_bins = 1; + for (auto dim : shape_) n_bins *= dim; + return n_bins; +} + +int RectilinearMesh::n_surface_bins() const +{ + return 4 * n_dimension_ * n_bins(); +} + +std::pair, std::vector> +RectilinearMesh::plot(Position plot_ll, Position plot_ur) const +{ + // Figure out which axes lie in the plane of the plot. + std::array axes {-1, -1}; + if (plot_ur.z == plot_ll.z) { + axes = {0, 1}; + } else if (plot_ur.y == plot_ll.y) { + axes = {0, 2}; + } else if (plot_ur.x == plot_ll.x) { + axes = {1, 2}; + } else { + fatal_error("Can only plot mesh lines on an axis-aligned plot"); + } + + // Get the coordinates of the mesh lines along both of the axes. + std::array, 2> axis_lines; + for (int i_ax = 0; i_ax < 2; ++i_ax) { + int axis = axes[i_ax]; + std::vector& lines {axis_lines[i_ax]}; + + for (auto coord : grid_[axis]) { + if (coord >= plot_ll[axis] && coord <= plot_ur[axis]) + lines.push_back(coord); + } + } + + return {axis_lines[0], axis_lines[1]}; +} + +void RectilinearMesh::to_hdf5(hid_t group) const +{ + hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); + + write_dataset(mesh_group, "type", "rectilinear"); + write_dataset(mesh_group, "x_grid", grid_[0]); + write_dataset(mesh_group, "y_grid", grid_[1]); + write_dataset(mesh_group, "z_grid", grid_[2]); + + close_group(mesh_group); +} + +bool RectilinearMesh::intersects(Position& r0, Position r1, int* ijk) const +{ + // Copy coordinates of starting point + double x0 = r0.x; + double y0 = r0.y; + double z0 = r0.z; + + // Copy coordinates of ending point + double x1 = r1.x; + double y1 = r1.y; + double z1 = r1.z; + + // Copy coordinates of mesh lower_left + double xm0 = grid_[0].front(); + double ym0 = grid_[1].front(); + double zm0 = grid_[2].front(); + + // Copy coordinates of mesh upper_right + double xm1 = grid_[0].back(); + double ym1 = grid_[1].back(); + double zm1 = grid_[2].back(); + + double min_dist = INFTY; + + // Check if line intersects left surface -- calculate the intersection point + // (y,z) + if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) { + double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0); + double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0); + if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) { + if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) { + ijk[0] = 1; + ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1; + ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1; + } + } + } + + // Check if line intersects back surface -- calculate the intersection point + // (x,z) + if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) { + double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0); + double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0); + if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) { + if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) { + ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1; + ijk[1] = 1; + ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1; + } + } + } + + // Check if line intersects bottom surface -- calculate the intersection + // point (x,y) + if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) { + double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0); + double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0); + if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) { + if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) { + ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1; + ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1; + ijk[2] = 1; + } + } + } + + // Check if line intersects right surface -- calculate the intersection point + // (y,z) + if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) { + double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0); + double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0); + if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) { + if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) { + ijk[0] = shape_[0]; + ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1; + ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1; + } + } + } + + // Check if line intersects front surface -- calculate the intersection point + // (x,z) + if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) { + double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0); + double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0); + if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) { + if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) { + ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1; + ijk[1] = shape_[1]; + ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1; + } + } + } + + // Check if line intersects top surface -- calculate the intersection point + // (x,y) + if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) { + double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0); + double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0); + if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) { + if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) { + ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1; + ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1; + ijk[2] = shape_[2]; + } + } + } + + return min_dist < INFTY; +} + +//============================================================================== +// Helper functions for the C API +//============================================================================== + +int +check_mesh(int32_t index) +{ + if (index < 0 || index >= model::meshes.size()) { + set_errmsg("Index in meshes array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + return 0; +} + +int +check_regular_mesh(int32_t index, RegularMesh** mesh) +{ + if (int err = check_mesh(index)) return err; + *mesh = dynamic_cast(model::meshes[index].get()); + if (!*mesh) { + set_errmsg("This function is only valid for regular meshes."); + return OPENMC_E_INVALID_TYPE; + } + return 0; +} + //============================================================================== // C API functions //============================================================================== @@ -817,10 +1386,7 @@ openmc_get_mesh_index(int32_t id, int32_t* index) extern "C" int openmc_mesh_get_id(int32_t index, int32_t* id) { - if (index < 0 || index >= model::meshes.size()) { - set_errmsg("Index in meshes array is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } + if (int err = check_mesh(index)) return err; *id = model::meshes[index]->id_; return 0; } @@ -829,10 +1395,7 @@ openmc_mesh_get_id(int32_t index, int32_t* id) extern "C" int openmc_mesh_set_id(int32_t index, int32_t id) { - if (index < 0 || index >= model::meshes.size()) { - set_errmsg("Index in meshes array is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } + if (int err = check_mesh(index)) return err; model::meshes[index]->id_ = id; model::mesh_map[id] = index; return 0; @@ -842,12 +1405,10 @@ openmc_mesh_set_id(int32_t index, int32_t id) extern "C" int openmc_mesh_get_dimension(int32_t index, int** dims, int* n) { - if (index < 0 || index >= model::meshes.size()) { - set_errmsg("Index in meshes array is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } - *dims = model::meshes[index]->shape_.data(); - *n = model::meshes[index]->n_dimension_; + RegularMesh* mesh; + if (int err = check_regular_mesh(index, &mesh)) return err; + *dims = mesh->shape_.data(); + *n = mesh->n_dimension_; return 0; } @@ -855,16 +1416,13 @@ openmc_mesh_get_dimension(int32_t index, int** dims, int* n) extern "C" int openmc_mesh_set_dimension(int32_t index, int n, const int* dims) { - if (index < 0 || index >= model::meshes.size()) { - set_errmsg("Index in meshes array is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } + RegularMesh* mesh; + if (int err = check_regular_mesh(index, &mesh)) return err; // Copy dimension std::vector shape = {static_cast(n)}; - auto& m = model::meshes[index]; - m->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape); - m->n_dimension_ = m->shape_.size(); + mesh->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape); + mesh->n_dimension_ = mesh->shape_.size(); return 0; } @@ -873,12 +1431,9 @@ openmc_mesh_set_dimension(int32_t index, int n, const int* dims) extern "C" int openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n) { - if (index < 0 || index >= model::meshes.size()) { - set_errmsg("Index in meshes array is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } + RegularMesh* m; + if (int err = check_regular_mesh(index, &m)) return err; - auto& m = model::meshes[index]; if (m->lower_left_.dimension() == 0) { set_errmsg("Mesh parameters have not been set."); return OPENMC_E_ALLOCATE; @@ -896,12 +1451,9 @@ extern "C" int openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, const double* width) { - if (index < 0 || index >= model::meshes.size()) { - set_errmsg("Index in meshes array is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } + RegularMesh* m; + if (int err = check_regular_mesh(index, &m)) return err; - auto& m = model::meshes[index]; std::vector shape = {static_cast(n)}; if (ll && ur) { m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape); @@ -930,8 +1482,21 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, void read_meshes(pugi::xml_node root) { for (auto node : root.children("mesh")) { + std::string mesh_type; + if (check_for_node(node, "type")) { + mesh_type = get_node_value(node, "type", true, true); + } else { + mesh_type = "regular"; + } + // Read mesh and add to vector - model::meshes.push_back(std::make_unique(node)); + if (mesh_type == "regular") { + model::meshes.push_back(std::make_unique(node)); + } else if (mesh_type == "rectilinear") { + model::meshes.push_back(std::make_unique(node)); + } else { + fatal_error("Invalid mesh type: " + mesh_type); + } // Map ID to position in vector model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1; diff --git a/src/plot.cpp b/src/plot.cpp index 57466eca1b..6bc82e1641 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -19,6 +19,7 @@ #include "openmc/progress_bar.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/string_utils.h" namespace openmc { @@ -497,20 +498,38 @@ Plot::set_meshlines(pugi::xml_node plot_node) // Set mesh based on type if ("ufs" == meshtype) { - if (settings::index_ufs_mesh < 0) { + if (!simulation::ufs_mesh) { std::stringstream err_msg; err_msg << "No UFS mesh for meshlines on plot " << id_; fatal_error(err_msg); } else { - index_meshlines_mesh_ = settings::index_ufs_mesh; + for (int i = 0; i < model::meshes.size(); ++i) { + if (const auto* m + = dynamic_cast(model::meshes[i].get())) { + if (m == simulation::ufs_mesh) { + index_meshlines_mesh_ = i; + } + } + } + if (index_meshlines_mesh_ == -1) + fatal_error("Could not find the UFS mesh for meshlines plot"); } } else if ("entropy" == meshtype) { - if (settings::index_entropy_mesh < 0) { + if (!simulation::entropy_mesh) { std::stringstream err_msg; - err_msg <<"No entropy mesh for meshlines on plot " << id_; + err_msg << "No entropy mesh for meshlines on plot " << id_; fatal_error(err_msg); } else { - index_meshlines_mesh_ = settings::index_entropy_mesh; + for (int i = 0; i < model::meshes.size(); ++i) { + if (const auto* m + = dynamic_cast(model::meshes[i].get())) { + if (m == simulation::entropy_mesh) { + index_meshlines_mesh_ = i; + } + } + } + if (index_meshlines_mesh_ == -1) + fatal_error("Could not find the entropy mesh for meshlines plot"); } } else if ("tally" == meshtype) { // Ensure that there is a mesh id if the type is tally @@ -679,19 +698,19 @@ void draw_mesh_lines(Plot pl, ImageData& data) RGBColor rgb; rgb = pl.meshlines_color_; - int outer, inner; + int ax1, ax2; switch(pl.basis_) { case PlotBasis::xy : - outer = 0; - inner = 1; + ax1 = 0; + ax2 = 1; break; case PlotBasis::xz : - outer = 0; - inner = 2; + ax1 = 0; + ax2 = 2; break; case PlotBasis::yz : - outer = 1; - inner = 2; + ax1 = 1; + ax2 = 2; break; default: UNREACHABLE(); @@ -700,70 +719,73 @@ void draw_mesh_lines(Plot pl, ImageData& data) Position ll_plot {pl.origin_}; Position ur_plot {pl.origin_}; - ll_plot[outer] -= pl.width_[0] / 2.; - ll_plot[inner] -= pl.width_[1] / 2.; - ur_plot[outer] += pl.width_[0] / 2.; - ur_plot[inner] += pl.width_[1] / 2.; + ll_plot[ax1] -= pl.width_[0] / 2.; + ll_plot[ax2] -= pl.width_[1] / 2.; + ur_plot[ax1] += pl.width_[0] / 2.; + ur_plot[ax2] += pl.width_[1] / 2.; Position width = ur_plot - ll_plot; - auto& m = model::meshes[pl.index_meshlines_mesh_]; + // Find the (axis-aligned) lines of the mesh that intersect this plot. + auto axis_lines = model::meshes[pl.index_meshlines_mesh_] + ->plot(ll_plot, ur_plot); - int ijk_ll[3], ijk_ur[3]; - bool in_mesh; - m->get_indices(ll_plot, &(ijk_ll[0]), &in_mesh); - m->get_indices(ur_plot, &(ijk_ur[0]), &in_mesh); + // Find the bounds along the second axis (accounting for low-D meshes). + int ax2_min, ax2_max; + if (axis_lines.second.size() > 0) { + double frac = (axis_lines.second.back() - ll_plot[ax2]) / width[ax2]; + ax2_min = (1.0 - frac) * pl.pixels_[1]; + if (ax2_min < 0) ax2_min = 0; + frac = (axis_lines.second.front() - ll_plot[ax2]) / width[ax2]; + ax2_max = (1.0 - frac) * pl.pixels_[1]; + if (ax2_max > pl.pixels_[1]) ax2_max = pl.pixels_[1]; + } else { + ax2_min = 0; + ax2_max = pl.pixels_[1]; + } - // Fortran/C++ index correction - ijk_ur[0]++; ijk_ur[1]++; ijk_ur[2]++; - - Position r_ll, r_ur; - // sweep through all meshbins on this plane and draw borders - for (int i = ijk_ll[outer]; i <= ijk_ur[outer]; i++) { - for (int j = ijk_ll[inner]; j <= ijk_ur[inner]; j++) { - // check if we're in the mesh for this ijk - if (i > 0 && i <= m->shape_[outer] && j >0 && j <= m->shape_[inner] ) { - int outrange[3], inrange[3]; - // get xyz's of lower left and upper right of this mesh cell - r_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1); - r_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1); - r_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i; - r_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j; - - // map the xyz ranges to pixel ranges - double frac = (r_ll[outer] - ll_plot[outer]) / width[outer]; - outrange[0] = int(frac * double(pl.pixels_[0])); - frac = (r_ur[outer] - ll_plot[outer]) / width[outer]; - outrange[1] = int(frac * double(pl.pixels_[0])); - - frac = (r_ur[inner] - ll_plot[inner]) / width[inner]; - inrange[0] = int((1. - frac) * (double)pl.pixels_[1]); - frac = (r_ll[inner] - ll_plot[inner]) / width[inner]; - inrange[1] = int((1. - frac) * (double)pl.pixels_[1]); - - // draw lines - for (int out_ = outrange[0]; out_ <= outrange[1]; out_++) { - for (int plus = 0; plus <= pl.meshlines_width_; plus++) { - data(out_, inrange[0] + plus) = rgb; - data(out_, inrange[1] + plus) = rgb; - data(out_, inrange[0] - plus) = rgb; - data(out_, inrange[1] - plus) = rgb; - } - } - - for (int in_ = inrange[0]; in_ <= inrange[1]; in_++) { - for (int plus = 0; plus <= pl.meshlines_width_; plus++) { - data(outrange[0] + plus, in_) = rgb; - data(outrange[1] + plus, in_) = rgb; - data(outrange[0] - plus, in_) = rgb; - data(outrange[1] - plus, in_) = rgb; - } - } - - } // end if(in mesh) + // Iterate across the first axis and draw lines. + for (auto ax1_val : axis_lines.first) { + double frac = (ax1_val - ll_plot[ax1]) / width[ax1]; + int ax1_ind = frac * pl.pixels_[0]; + for (int ax2_ind = ax2_min; ax2_ind < ax2_max; ++ax2_ind) { + for (int plus = 0; plus <= pl.meshlines_width_; plus++) { + if (ax1_ind+plus >= 0 && ax1_ind+plus < pl.pixels_[0]) + data(ax1_ind+plus, ax2_ind) = rgb; + if (ax1_ind-plus >= 0 && ax1_ind-plus < pl.pixels_[0]) + data(ax1_ind-plus, ax2_ind) = rgb; + } } - } // end outer loops -} // end draw_mesh_lines + } + + // Find the bounds along the first axis. + int ax1_min, ax1_max; + if (axis_lines.first.size() > 0) { + double frac = (axis_lines.first.front() - ll_plot[ax1]) / width[ax1]; + ax1_min = frac * pl.pixels_[0]; + if (ax1_min < 0) ax1_min = 0; + frac = (axis_lines.first.back() - ll_plot[ax1]) / width[ax1]; + ax1_max = frac * pl.pixels_[0]; + if (ax1_max > pl.pixels_[0]) ax1_max = pl.pixels_[0]; + } else { + ax1_min = 0; + ax1_max = pl.pixels_[0]; + } + + // Iterate across the second axis and draw lines. + for (auto ax2_val : axis_lines.second) { + double frac = (ax2_val - ll_plot[ax2]) / width[ax2]; + int ax2_ind = (1.0 - frac) * pl.pixels_[1]; + for (int ax1_ind = ax1_min; ax1_ind < ax1_max; ++ax1_ind) { + for (int plus = 0; plus <= pl.meshlines_width_; plus++) { + if (ax2_ind+plus >= 0 && ax2_ind+plus < pl.pixels_[1]) + data(ax1_ind, ax2_ind+plus) = rgb; + if (ax2_ind-plus >= 0 && ax2_ind-plus < pl.pixels_[1]) + data(ax1_ind, ax2_ind-plus) = rgb; + } + } + } +} //============================================================================== // CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 204284c480..ec511481f6 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -1,17 +1,30 @@ element tallies { element mesh { (element id { xsd:int } | attribute id { xsd:int }) & - (element type { ( "regular" ) } | - attribute type { ( "regular" ) }) & - (element dimension { list { xsd:positiveInteger+ } } | - attribute dimension { list { xsd:positiveInteger+ } }) & - (element lower_left { list { xsd:double+ } } | - attribute lower_left { list { xsd:double+ } }) & ( - (element upper_right { list { xsd:double+ } } | - attribute upper_right { list { xsd:double+ } }) | - (element width { list { xsd:double+ } } | - attribute width { list { xsd:double+ } }) + ( + (element type { ( "regular" ) } | + attribute type { ( "regular" ) }) & + (element dimension { list { xsd:positiveInteger+ } } | + attribute dimension { list { xsd:positiveInteger+ } }) & + (element lower_left { list { xsd:double+ } } | + attribute lower_left { list { xsd:double+ } }) & + ( + (element upper_right { list { xsd:double+ } } | + attribute upper_right { list { xsd:double+ } }) | + (element width { list { xsd:double+ } } | + attribute width { list { xsd:double+ } }) + ) + ) | ( + (element type { ( "rectilinear" ) } | + attribute type { ( "rectilinear" ) }) & + (element x_grid { list { xsd:double+ } } | + attribute x_grid { list { xsd:double+ } }) & + (element y_grid { list { xsd:double+ } } | + attribute y_grid { list { xsd:double+ } }) & + (element z_grid { list { xsd:double+ } } | + attribute z_grid { list { xsd:double+ } }) + ) ) }* & diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index 15b6f5b249..98a48eeb85 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -13,78 +13,140 @@ - - regular - - - regular - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + regular + + + regular + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + rectilinear + + + rectilinear + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/settings.cpp b/src/settings.cpp index 32032ed2a7..19ec9eabda 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -73,9 +73,6 @@ std::string path_source; std::string path_sourcepoint; std::string path_statepoint; -int32_t index_entropy_mesh {-1}; -int32_t index_ufs_mesh {-1}; - int32_t n_batches; int32_t n_inactive {0}; int32_t gen_per_batch {1}; @@ -481,6 +478,7 @@ void read_settings_xml() read_meshes(root); // Shannon Entropy mesh + int32_t index_entropy_mesh = -1; if (check_for_node(root, "entropy_mesh")) { int temp = std::stoi(get_node_value(root, "entropy_mesh")); if (model::mesh_map.find(temp) == model::mesh_map.end()) { @@ -508,17 +506,21 @@ void read_settings_xml() } if (index_entropy_mesh >= 0) { - auto& m = *model::meshes[index_entropy_mesh]; - if (m.shape_.dimension() == 0) { + auto* m = dynamic_cast( + model::meshes[index_entropy_mesh].get()); + if (!m) fatal_error("Only regular meshes can be used as an entropy mesh"); + simulation::entropy_mesh = m; + + if (m->shape_.dimension() == 0) { // If the user did not specify how many mesh cells are to be used in // each direction, we automatically determine an appropriate number of // cells int n = std::ceil(std::pow(n_particles / 20.0, 1.0/3.0)); - m.shape_ = {n, n, n}; - m.n_dimension_ = 3; + m->shape_ = {n, n, n}; + m->n_dimension_ = 3; // Calculate width - m.width_ = (m.upper_right_ - m.lower_left_) / m.shape_; + m->width_ = (m->upper_right_ - m->lower_left_) / m->shape_; } // Turn on Shannon entropy calculation @@ -526,6 +528,7 @@ void read_settings_xml() } // Uniform fission source weighting mesh + int32_t i_ufs_mesh = -1; if (check_for_node(root, "ufs_mesh")) { auto temp = std::stoi(get_node_value(root, "ufs_mesh")); if (model::mesh_map.find(temp) == model::mesh_map.end()) { @@ -534,7 +537,7 @@ void read_settings_xml() "does not exist."; fatal_error(msg); } - index_ufs_mesh = model::mesh_map.at(temp); + i_ufs_mesh = model::mesh_map.at(temp); } else if (check_for_node(root, "uniform_fs")) { warning("Specifying a UFS mesh via the element " @@ -546,14 +549,18 @@ void read_settings_xml() model::meshes.push_back(std::make_unique(node_ufs)); // Set entropy mesh index - index_ufs_mesh = model::meshes.size() - 1; + i_ufs_mesh = model::meshes.size() - 1; // Assign ID and set mapping model::meshes.back()->id_ = 10001; model::mesh_map[10001] = index_entropy_mesh; } - if (index_ufs_mesh >= 0) { + if (i_ufs_mesh >= 0) { + auto* m = dynamic_cast(model::meshes[i_ufs_mesh].get()); + if (!m) fatal_error("Only regular meshes can be used as a UFS mesh"); + simulation::ufs_mesh = m; + // Turn on uniform fission source weighting ufs_on = true; } diff --git a/src/simulation.cpp b/src/simulation.cpp index fc7a788517..f94ba232ac 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -247,6 +247,9 @@ int total_gen {0}; double total_weight; int64_t work_per_rank; +const RegularMesh* entropy_mesh {nullptr}; +const RegularMesh* ufs_mesh {nullptr}; + std::vector k_generation; std::vector work_index; diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index f69a966523..5dce9c618a 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -74,8 +74,7 @@ void MeshFilter::set_mesh(int32_t mesh) { mesh_ = mesh; - n_bins_ = 1; - for (auto dim : model::meshes[mesh_]->shape_) n_bins_ *= dim; + n_bins_ = model::meshes[mesh_]->n_bins(); } //============================================================================== diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp index 57f10b09c3..edb0492e9f 100644 --- a/src/tallies/filter_meshsurface.cpp +++ b/src/tallies/filter_meshsurface.cpp @@ -75,8 +75,7 @@ void MeshSurfaceFilter::set_mesh(int32_t mesh) { mesh_ = mesh; - n_bins_ = 4 * model::meshes[mesh_]->n_dimension_; - for (auto dim : model::meshes[mesh_]->shape_) n_bins_ *= dim; + n_bins_ = model::meshes[mesh_]->n_surface_bins(); } //============================================================================== diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat index 511b298482..6755a64c82 100644 --- a/tests/regression_tests/filter_mesh/inputs_true.dat +++ b/tests/regression_tests/filter_mesh/inputs_true.dat @@ -309,45 +309,56 @@ - + 17 -182.07 182.07 - + 17 17 -182.07 -182.07 182.07 182.07 - + 17 17 17 -182.07 -182.07 -183.0 182.07 182.07 183.0 + + -182.07 -160.65 -139.23 -117.81 -96.39 -74.97 -53.55000000000001 -32.129999999999995 -10.710000000000008 10.70999999999998 32.129999999999995 53.54999999999998 74.96999999999997 96.38999999999999 117.81 139.22999999999996 160.64999999999998 182.07 + -182.07 -160.65 -139.23 -117.81 -96.39 -74.97 -53.55000000000001 -32.129999999999995 -10.710000000000008 10.70999999999998 32.129999999999995 53.54999999999998 74.96999999999997 96.38999999999999 117.81 139.22999999999996 160.64999999999998 182.07 + 1.0 1.683624003879018 2.8345897864376153 4.772383405596668 8.034899257376447 13.52774925846868 22.77564337001445 38.34561988154435 64.55960607618856 108.69410247084474 182.99999999999991 + 1 - + 1 2 - + 2 3 - + 3 + + 4 + + + 4 + 1 total - 4 + 5 current @@ -355,7 +366,7 @@ total - 5 + 6 current @@ -363,7 +374,15 @@ total - 6 + 7 + current + + + 4 + total + + + 8 current diff --git a/tests/regression_tests/filter_mesh/results_true.dat b/tests/regression_tests/filter_mesh/results_true.dat index 227d0d608a..05efffe2fb 100644 --- a/tests/regression_tests/filter_mesh/results_true.dat +++ b/tests/regression_tests/filter_mesh/results_true.dat @@ -1 +1 @@ -cc306dfd3e3669827cb2f3bdc38005c0616f6faecb9277f3b0f869fdecf826ae662c65d67a402705956f100d2e13c7c11332364fa41ccc2e215ab100709b627f \ No newline at end of file +35f04a6f062ef64116ef4eb0e9b803cd44cff7e185e2b53c9174afad8a26ca1a436ca9b800d6a228e006a9129f4536d7dce289d7a11cd56c6949d71d6a201b31 \ No newline at end of file diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py index 8a7f7a6fca..ed2e9b35b4 100644 --- a/tests/regression_tests/filter_mesh/test.py +++ b/tests/regression_tests/filter_mesh/test.py @@ -1,3 +1,5 @@ +import numpy as np + import openmc from tests.testing_harness import HashedPyAPITestHarness @@ -8,31 +10,35 @@ class FilterMeshTestHarness(HashedPyAPITestHarness): super().__init__(*args, **kwargs) # Initialize Meshes - mesh_1d = openmc.Mesh(mesh_id=1) - mesh_1d.type = 'regular' + mesh_1d = openmc.RegularMesh(mesh_id=1) mesh_1d.dimension = [17] mesh_1d.lower_left = [-182.07] mesh_1d.upper_right = [182.07] - mesh_2d = openmc.Mesh(mesh_id=2) - mesh_2d.type = 'regular' + mesh_2d = openmc.RegularMesh(mesh_id=2) mesh_2d.dimension = [17, 17] mesh_2d.lower_left = [-182.07, -182.07] mesh_2d.upper_right = [182.07, 182.07] - mesh_3d = openmc.Mesh(mesh_id=3) - mesh_3d.type = 'regular' + mesh_3d = openmc.RegularMesh(mesh_id=3) mesh_3d.dimension = [17, 17, 17] mesh_3d.lower_left = [-182.07, -182.07, -183.00] mesh_3d.upper_right = [182.07, 182.07, 183.00] + recti_mesh = openmc.RectilinearMesh(mesh_id=4) + recti_mesh.x_grid = np.linspace(-182.07, 182.07, 18) + recti_mesh.y_grid = np.linspace(-182.07, 182.07, 18) + recti_mesh.z_grid = np.logspace(0, np.log10(183), 11) + # Initialize the filters mesh_1d_filter = openmc.MeshFilter(mesh_1d) mesh_2d_filter = openmc.MeshFilter(mesh_2d) mesh_3d_filter = openmc.MeshFilter(mesh_3d) + recti_mesh_filter = openmc.MeshFilter(recti_mesh) meshsurf_1d_filter = openmc.MeshSurfaceFilter(mesh_1d) meshsurf_2d_filter = openmc.MeshSurfaceFilter(mesh_2d) meshsurf_3d_filter = openmc.MeshSurfaceFilter(mesh_3d) + recti_meshsurf_filter = openmc.MeshSurfaceFilter(recti_mesh) # Initialized the tallies tally = openmc.Tally(name='tally 1') @@ -65,6 +71,16 @@ class FilterMeshTestHarness(HashedPyAPITestHarness): tally.scores = ['current'] self._model.tallies.append(tally) + tally = openmc.Tally(name='tally 7') + tally.filters = [recti_mesh_filter] + tally.scores = ['total'] + self._model.tallies.append(tally) + + tally = openmc.Tally(name='tally 8') + tally.filters = [recti_meshsurf_filter] + tally.scores = ['current'] + self._model.tallies.append(tally) + def test_filter_mesh(): harness = FilterMeshTestHarness('statepoint.10.h5') diff --git a/tests/regression_tests/mg_tallies/inputs_true.dat b/tests/regression_tests/mg_tallies/inputs_true.dat index b75573b0a7..5629587222 100644 --- a/tests/regression_tests/mg_tallies/inputs_true.dat +++ b/tests/regression_tests/mg_tallies/inputs_true.dat @@ -33,7 +33,7 @@ - + 10 1 1 0.0 0.0 0.0 929.45 1000 1000 diff --git a/tests/regression_tests/mg_tallies/test.py b/tests/regression_tests/mg_tallies/test.py index b0a8987108..e1c19e6b7e 100644 --- a/tests/regression_tests/mg_tallies/test.py +++ b/tests/regression_tests/mg_tallies/test.py @@ -62,8 +62,7 @@ def test_mg_tallies(): model = slab_mg() # Instantiate a tally mesh - mesh = openmc.Mesh(mesh_id=1) - mesh.type = 'regular' + mesh = openmc.RegularMesh(mesh_id=1) mesh.dimension = [10, 1, 1] mesh.lower_left = [0.0, 0.0, 0.0] mesh.upper_right = [929.45, 1000, 1000] diff --git a/tests/regression_tests/mgxs_library_mesh/inputs_true.dat b/tests/regression_tests/mgxs_library_mesh/inputs_true.dat index 299da0713c..4f138b24f6 100644 --- a/tests/regression_tests/mgxs_library_mesh/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_mesh/inputs_true.dat @@ -309,7 +309,7 @@ - + 2 2 -100.0 -100.0 100.0 100.0 diff --git a/tests/regression_tests/mgxs_library_mesh/test.py b/tests/regression_tests/mgxs_library_mesh/test.py index a9d24da934..b72fdf2f1e 100644 --- a/tests/regression_tests/mgxs_library_mesh/test.py +++ b/tests/regression_tests/mgxs_library_mesh/test.py @@ -27,8 +27,7 @@ class MGXSTestHarness(PyAPITestHarness): self.mgxs_lib.domain_type = 'mesh' # Instantiate a tally mesh - mesh = openmc.Mesh(mesh_id=1) - mesh.type = 'regular' + mesh = openmc.RegularMesh(mesh_id=1) mesh.dimension = [2, 2] mesh.lower_left = [-100., -100.] mesh.width = [100., 100.] diff --git a/tests/regression_tests/plot/plots.xml b/tests/regression_tests/plot/plots.xml index ecfe69125b..ce63da1442 100644 --- a/tests/regression_tests/plot/plots.xml +++ b/tests/regression_tests/plot/plots.xml @@ -14,6 +14,7 @@ 25 25 200 200 + diff --git a/tests/regression_tests/plot/results_true.dat b/tests/regression_tests/plot/results_true.dat index 18251e8d0c..db6beb9cdb 100644 --- a/tests/regression_tests/plot/results_true.dat +++ b/tests/regression_tests/plot/results_true.dat @@ -1 +1 @@ -6b3eb36488b995b42233e6b7c0164cf6c92a1823b3c91985b97fd076f86c0aad93fbdb74e70026f5f12c61a6aee52a09588171a7fd839511ee7d9efc3952beb2 \ No newline at end of file +fd61f6a5de632f06a39e2a938480ba3dc314810655e684291e4b255bb8c142a16f10fec414c587ecd727fa559d3760d6e5b43f53fe86e90238a69d9c5698db8e \ No newline at end of file diff --git a/tests/regression_tests/plot/tallies.xml b/tests/regression_tests/plot/tallies.xml new file mode 100644 index 0000000000..b7e678ca0f --- /dev/null +++ b/tests/regression_tests/plot/tallies.xml @@ -0,0 +1,20 @@ + + + + + -10 10 + -10 10 + -10 0 5 7.5 8.75 10 + + + + mesh + 2 + + + + 1 + total + + + diff --git a/tests/regression_tests/tallies/inputs_true.dat b/tests/regression_tests/tallies/inputs_true.dat index 3ec8d8d4ce..9409ca4a71 100644 --- a/tests/regression_tests/tallies/inputs_true.dat +++ b/tests/regression_tests/tallies/inputs_true.dat @@ -309,7 +309,7 @@ - + 2 2 -182.07 -182.07 182.07 182.07 diff --git a/tests/regression_tests/tallies/test.py b/tests/regression_tests/tallies/test.py index 543d0acd21..d4857b3ec7 100644 --- a/tests/regression_tests/tallies/test.py +++ b/tests/regression_tests/tallies/test.py @@ -1,6 +1,6 @@ from openmc.filter import * from openmc.filter_expansion import * -from openmc import Mesh, Tally +from openmc import RegularMesh, Tally from tests.testing_harness import HashedPyAPITestHarness @@ -28,7 +28,7 @@ def test_tallies(): azimuthal_tally2.scores = ['flux'] azimuthal_tally2.estimator = 'analog' - mesh_2x2 = Mesh(mesh_id=1) + mesh_2x2 = RegularMesh(mesh_id=1) mesh_2x2.lower_left = [-182.07, -182.07] mesh_2x2.upper_right = [182.07, 182.07] mesh_2x2.dimension = [2, 2] diff --git a/tests/regression_tests/tally_arithmetic/inputs_true.dat b/tests/regression_tests/tally_arithmetic/inputs_true.dat index 3605005add..842b0df623 100644 --- a/tests/regression_tests/tally_arithmetic/inputs_true.dat +++ b/tests/regression_tests/tally_arithmetic/inputs_true.dat @@ -309,7 +309,7 @@ - + 2 2 2 -160.0 -160.0 -183.0 160.0 160.0 183.0 diff --git a/tests/regression_tests/tally_arithmetic/test.py b/tests/regression_tests/tally_arithmetic/test.py index 3adf0c5bef..236f5cc384 100644 --- a/tests/regression_tests/tally_arithmetic/test.py +++ b/tests/regression_tests/tally_arithmetic/test.py @@ -10,8 +10,7 @@ class TallyArithmeticTestHarness(PyAPITestHarness): super().__init__(*args, **kwargs) # Initialize Mesh - mesh = openmc.Mesh(mesh_id=1) - mesh.type = 'regular' + mesh = openmc.RegularMesh(mesh_id=1) mesh.dimension = [2, 2, 2] mesh.lower_left = [-160.0, -160.0, -183.0] mesh.upper_right = [160.0, 160.0, 183.0] diff --git a/tests/regression_tests/tally_slice_merge/inputs_true.dat b/tests/regression_tests/tally_slice_merge/inputs_true.dat index ccd3655f9e..6f421b1d95 100644 --- a/tests/regression_tests/tally_slice_merge/inputs_true.dat +++ b/tests/regression_tests/tally_slice_merge/inputs_true.dat @@ -310,7 +310,7 @@ - + 2 2 -50.0 -50.0 50.0 50.0 diff --git a/tests/regression_tests/tally_slice_merge/test.py b/tests/regression_tests/tally_slice_merge/test.py index 20ab57a077..090e314483 100644 --- a/tests/regression_tests/tally_slice_merge/test.py +++ b/tests/regression_tests/tally_slice_merge/test.py @@ -24,8 +24,7 @@ class TallySliceMergeTestHarness(PyAPITestHarness): cell_27 = openmc.CellFilter(27) distribcell_filter = openmc.DistribcellFilter(21) - mesh = openmc.Mesh(name='mesh') - mesh.type = 'regular' + mesh = openmc.RegularMesh(name='mesh') mesh.dimension = [2, 2] mesh.lower_left = [-50., -50.] mesh.upper_right = [+50., +50.] diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 418fe1d3f1..17219c82e4 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -336,11 +336,11 @@ def test_find_material(capi_init): def test_mesh(capi_init): - mesh = openmc.capi.Mesh() + mesh = openmc.capi.RegularMesh() mesh.dimension = (2, 3, 4) assert mesh.dimension == (2, 3, 4) with pytest.raises(exc.AllocationError): - mesh2 = openmc.capi.Mesh(mesh.id) + mesh2 = openmc.capi.RegularMesh(mesh.id) # Make sure each combination of parameters works ll = (0., 0., 0.) @@ -360,7 +360,7 @@ def test_mesh(capi_init): assert isinstance(meshes, Mapping) assert len(meshes) == 1 for mesh_id, mesh in meshes.items(): - assert isinstance(mesh, openmc.capi.Mesh) + assert isinstance(mesh, openmc.capi.RegularMesh) assert mesh_id == mesh.id mf = openmc.capi.MeshFilter(mesh) diff --git a/tests/unit_tests/test_mesh_from_lattice.py b/tests/unit_tests/test_mesh_from_lattice.py index feac873260..0a01387cde 100644 --- a/tests/unit_tests/test_mesh_from_lattice.py +++ b/tests/unit_tests/test_mesh_from_lattice.py @@ -90,12 +90,12 @@ def test_mesh2d(rlat2): shape = np.array(rlat2.shape) width = shape*rlat2.pitch - mesh1 = openmc.Mesh.from_rect_lattice(rlat2) + mesh1 = openmc.RegularMesh.from_rect_lattice(rlat2) assert np.array_equal(mesh1.dimension, (3, 3)) assert np.array_equal(mesh1.lower_left, rlat2.lower_left) assert np.array_equal(mesh1.upper_right, rlat2.lower_left + width) - mesh2 = openmc.Mesh.from_rect_lattice(rlat2, division=3) + mesh2 = openmc.RegularMesh.from_rect_lattice(rlat2, division=3) assert np.array_equal(mesh2.dimension, (9, 9)) assert np.array_equal(mesh2.lower_left, rlat2.lower_left) assert np.array_equal(mesh2.upper_right, rlat2.lower_left + width) @@ -105,12 +105,12 @@ def test_mesh3d(rlat3): shape = np.array(rlat3.shape) width = shape*rlat3.pitch - mesh1 = openmc.Mesh.from_rect_lattice(rlat3) + mesh1 = openmc.RegularMesh.from_rect_lattice(rlat3) assert np.array_equal(mesh1.dimension, (3, 3, 2)) assert np.array_equal(mesh1.lower_left, rlat3.lower_left) assert np.array_equal(mesh1.upper_right, rlat3.lower_left + width) - mesh2 = openmc.Mesh.from_rect_lattice(rlat3, division=3) + mesh2 = openmc.RegularMesh.from_rect_lattice(rlat3, division=3) assert np.array_equal(mesh2.dimension, (9, 9, 6)) assert np.array_equal(mesh2.lower_left, rlat3.lower_left) assert np.array_equal(mesh2.upper_right, rlat3.lower_left + width) diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index f99a32a52c..9c33f6a2d4 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -25,7 +25,7 @@ def test_export_to_xml(run_in_tmpdir): s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5, 'energy_photon': 1000.0, 'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5} - mesh = openmc.Mesh() + mesh = openmc.RegularMesh() mesh.lower_left = (-10., -10., -10.) mesh.upper_right = (10., 10., 10.) mesh.dimension = (5, 5, 5) @@ -79,7 +79,7 @@ def test_export_to_xml(run_in_tmpdir): assert s.cutoff == {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5, 'energy_photon': 1000.0, 'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5} - assert isinstance(s.entropy_mesh, openmc.Mesh) + assert isinstance(s.entropy_mesh, openmc.RegularMesh) assert s.entropy_mesh.lower_left == [-10., -10., -10.] assert s.entropy_mesh.upper_right == [10., 10., 10.] assert s.entropy_mesh.dimension == [5, 5, 5] @@ -92,7 +92,7 @@ def test_export_to_xml(run_in_tmpdir): 'multipole': True, 'range': [200., 1000.]} assert s.trace == [10, 1, 20] assert s.track == [1, 1, 1, 2, 1, 1] - assert isinstance(s.ufs_mesh, openmc.Mesh) + assert isinstance(s.ufs_mesh, openmc.RegularMesh) assert s.ufs_mesh.lower_left == [-10., -10., -10.] assert s.ufs_mesh.upper_right == [10., 10., 10.] assert s.ufs_mesh.dimension == [5, 5, 5]