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

Merge pull request #1246 from smharper/rect_mesh

Browse source 
Add rectilinear tally meshes
This commit is contained in:
Paul Romano 2019-06-12 22:40:46 -05:00 • committed by GitHub
commit c0512b3213
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
60 changed files with 1531 additions and 644 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -312,21 +312,30 @@ a separate element with the tag name ``<mesh>``. 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 ``<upper_right>`` or ``<width>`` must be specified, but not both

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

@ -124,7 +124,8 @@ Constructing Tallies
openmc.ZernikeFilter
openmc.ZernikeRadialFilter
openmc.ParticleFilter
openmc.Mesh
openmc.RegularMesh
openmc.RectilinearMesh
openmc.Trigger
openmc.TallyDerivative
openmc.Tally

81
examples/jupyter/mdgxs-part-ii.ipynb
View file

@ -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
}

45
examples/jupyter/mg-mode-part-i.ipynb
View file

@ -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,

9
examples/jupyter/mg-mode-part-ii.ipynb
View file

@ -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,

9
examples/jupyter/mg-mode-part-iii.ipynb
View file

@ -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,

9
examples/jupyter/mgxs-part-iii.ipynb
View file

@ -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,

5
examples/jupyter/pandas-dataframes.ipynb
View file

@ -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,

2
examples/jupyter/post-processing.ipynb
View file

@ -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",

74
examples/jupyter/tally-arithmetic.ipynb
View file

@ -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,

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

@ -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()

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

@ -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]

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

@ -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]

2
examples/python/pincell_depletion/restart_depletion.py
View file

@ -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]

2
examples/python/pincell_depletion/run_depletion.py
View file

@ -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]

3
examples/python/pincell_multigroup/build-xml.py
View file

@ -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]

166
include/openmc/mesh.h
View file

@ -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<std::unique_ptr<RegularMesh>> meshes;
extern std::vector<std::unique_ptr<Mesh>> meshes;
extern std::unordered_map<int32_t, int32_t> 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<int>& bins,
std::vector<double>& lengths) const;
virtual void bins_crossed(const Particle* p, std::vector<int>& bins,
std::vector<double>& 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<int>& bins) const;
virtual void
surface_bins_crossed(const Particle* p, std::vector<int>& 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<double>, std::vector<double>>
plot(Position plot_ll, Position plot_ur) const = 0;
//! Write mesh data to an HDF5 group
//
//! \param[in] group HDF5 group
virtual void to_hdf5(hid_t group) const = 0;
// Data members
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions
xt::xtensor<double, 1> lower_left_; //!< Lower-left coordinates of mesh
xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
};
//==============================================================================
//! 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<int>& bins,
std::vector<double>& lengths) const override;
void surface_bins_crossed(const Particle* p, std::vector<int>& 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<double>, std::vector<double>>
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<double> count_sites(const std::vector<Particle::Bank>& bank,
xt::xtensor<double, 1> count_sites(const std::vector<Particle::Bank>& 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<int> shape_; //!< Number of mesh elements in each dimension
xt::xarray<double> lower_left_; //!< Lower-left coordinates of mesh
xt::xarray<double> upper_right_; //!< Upper-right coordinates of mesh
xt::xarray<double> width_; //!< Width of each mesh element
xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
xt::xtensor<double, 1> 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<int>& bins,
std::vector<double>& lengths) const override;
void surface_bins_crossed(const Particle* p, std::vector<int>& 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<double>, std::vector<double>>
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<int, 1> shape_; //!< Number of mesh elements in each dimension
private:
std::vector<std::vector<double>> grid_;
};
//==============================================================================
// Non-member functions
//==============================================================================

3
include/openmc/settings.h
View file

@ -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

4
include/openmc/simulation.h
View file

@ -4,6 +4,7 @@
#ifndef OPENMC_SIMULATION_H
#define OPENMC_SIMULATION_H
#include "openmc/mesh.h"
#include "openmc/particle.h"
#include <cstdint>
@ -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<double> k_generation;
extern std::vector<int64_t> work_index;

6
openmc/capi/filter.py
View file

@ -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):

10
openmc/capi/mesh.py
View file

@ -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()

8
openmc/cmfd.py
View file

@ -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

37
openmc/filter.py
View file

@ -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

270
openmc/mesh.py
View file

@ -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

14
openmc/mgxs/library.py
View file

@ -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'])

30
openmc/mgxs/mdgxs.py
View file

@ -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

66
openmc/mgxs/mgxs.py
View file

@ -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

4
openmc/plots.py
View file

@ -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
"""

14
openmc/settings.py
View file

@ -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')

6
openmc/statepoint.py
View file

@ -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

1
openmc/trigger.py
View file

@ -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 = []

24
src/eigenvalue.cpp
View file

@ -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<double, 1> p = m->count_sites(simulation::fission_bank,
&sites_outside);
xt::xtensor<double, 1> 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;
}

5
src/finalize.cpp
View file

@ -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;

655
src/mesh.cpp
View file

@ -31,7 +31,7 @@ namespace openmc {
namespace model {
std::vector<std::unique_ptr<RegularMesh>> meshes;
std::vector<std::unique_ptr<Mesh>> meshes;
std::unordered_map<int32_t, int32_t> 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<int>(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<double>, std::vector<double>>
RegularMesh::plot(Position plot_ll, Position plot_ur) const
{
// Figure out which axes lie in the plane of the plot.
std::array<int, 2> 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<std::vector<double>, 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<double>
xt::xtensor<double, 1>
RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
bool* outside) const
{
@ -783,6 +829,529 @@ RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
return counts;
}
//==============================================================================
// RectilinearMesh implementation
//==============================================================================
RectilinearMesh::RectilinearMesh(pugi::xml_node node)
: Mesh {node}
{
n_dimension_ = 3;
grid_.resize(3);
grid_[0] = get_node_array<double>(node, "x_grid");
grid_[1] = get_node_array<double>(node, "y_grid");
grid_[2] = get_node_array<double>(node, "z_grid");
shape_ = {static_cast<int>(grid_[0].size()) - 1,
static_cast<int>(grid_[1].size()) - 1,
static_cast<int>(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<int>& bins,
std::vector<double>& 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<int>& 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<double>, std::vector<double>>
RectilinearMesh::plot(Position plot_ll, Position plot_ur) const
{
// Figure out which axes lie in the plane of the plot.
std::array<int, 2> 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<std::vector<double>, 2> axis_lines;
for (int i_ax = 0; i_ax < 2; ++i_ax) {
int axis = axes[i_ax];
std::vector<double>& 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<RegularMesh*>(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<std::size_t> shape = {static_cast<std::size_t>(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<std::size_t> shape = {static_cast<std::size_t>(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<RegularMesh>(node));
if (mesh_type == "regular") {
model::meshes.push_back(std::make_unique<RegularMesh>(node));
} else if (mesh_type == "rectilinear") {
model::meshes.push_back(std::make_unique<RectilinearMesh>(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;

162
src/plot.cpp
View file

@ -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<const RegularMesh*>(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<const RegularMesh*>(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

33
src/relaxng/tallies.rnc
View file

@ -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+ } })
)
)
}* &

206
src/relaxng/tallies.rng
View file

@ -13,78 +13,140 @@
</attribute>
</choice>
<choice>
<element name="type">
<value>regular</value>
</element>
<attribute name="type">
<value>regular</value>
</attribute>
</choice>
<choice>
<element name="dimension">
<list>
<oneOrMore>
<data type="positiveInteger"/>
</oneOrMore>
</list>
</element>
<attribute name="dimension">
<list>
<oneOrMore>
<data type="positiveInteger"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="lower_left">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="lower_left">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<choice>
<element name="upper_right">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="upper_right">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="width">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="width">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<interleave>
<choice>
<element name="type">
<value>regular</value>
</element>
<attribute name="type">
<value>regular</value>
</attribute>
</choice>
<choice>
<element name="dimension">
<list>
<oneOrMore>
<data type="positiveInteger"/>
</oneOrMore>
</list>
</element>
<attribute name="dimension">
<list>
<oneOrMore>
<data type="positiveInteger"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="lower_left">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="lower_left">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<choice>
<element name="upper_right">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="upper_right">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="width">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="width">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</choice>
</interleave>
<interleave>
<choice>
<element name="type">
<value>rectilinear</value>
</element>
<attribute name="type">
<value>rectilinear</value>
</attribute>
</choice>
<choice>
<element name="x_grid">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="x_grid">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="y_grid">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="y_grid">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="z_grid">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="z_grid">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</interleave>
</choice>
</interleave>
</element>

29
src/settings.cpp
View file

@ -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<RegularMesh*>(
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 <uniform_fs> element "
@ -546,14 +549,18 @@ void read_settings_xml()
model::meshes.push_back(std::make_unique<RegularMesh>(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<RegularMesh*>(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;
}

3
src/simulation.cpp
View file

@ -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<double> k_generation;
std::vector<int64_t> work_index;

3
src/tallies/filter_mesh.cpp
View file

@ -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();
}
//==============================================================================

3
src/tallies/filter_meshsurface.cpp
View file

@ -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();
}
//==============================================================================

37
tests/regression_tests/filter_mesh/inputs_true.dat
View file

@ -309,45 +309,56 @@
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" type="regular">
<mesh id="1">
<dimension>17</dimension>
<lower_left>-182.07</lower_left>
<upper_right>182.07</upper_right>
</mesh>
<mesh id="2" type="regular">
<mesh id="2">
<dimension>17 17</dimension>
<lower_left>-182.07 -182.07</lower_left>
<upper_right>182.07 182.07</upper_right>
</mesh>
<mesh id="3" type="regular">
<mesh id="3">
<dimension>17 17 17</dimension>
<lower_left>-182.07 -182.07 -183.0</lower_left>
<upper_right>182.07 182.07 183.0</upper_right>
</mesh>
<mesh id="4" type="rectilinear">
<x_grid>-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</x_grid>
<y_grid>-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</y_grid>
<z_grid>1.0 1.683624003879018 2.8345897864376153 4.772383405596668 8.034899257376447 13.52774925846868 22.77564337001445 38.34561988154435 64.55960607618856 108.69410247084474 182.99999999999991</z_grid>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="4" type="meshsurface">
<filter id="5" type="meshsurface">
<bins>1</bins>
</filter>
<filter id="2" type="mesh">
<bins>2</bins>
</filter>
<filter id="5" type="meshsurface">
<filter id="6" type="meshsurface">
<bins>2</bins>
</filter>
<filter id="3" type="mesh">
<bins>3</bins>
</filter>
<filter id="6" type="meshsurface">
<filter id="7" type="meshsurface">
<bins>3</bins>
</filter>
<filter id="4" type="mesh">
<bins>4</bins>
</filter>
<filter id="8" type="meshsurface">
<bins>4</bins>
</filter>
<tally id="1" name="tally 1">
<filters>1</filters>
<scores>total</scores>
</tally>
<tally id="2" name="tally 2">
<filters>4</filters>
<filters>5</filters>
<scores>current</scores>
</tally>
<tally id="3" name="tally 3">
@ -355,7 +366,7 @@
<scores>total</scores>
</tally>
<tally id="4" name="tally 4">
<filters>5</filters>
<filters>6</filters>
<scores>current</scores>
</tally>
<tally id="5" name="tally 5">
@ -363,7 +374,15 @@
<scores>total</scores>
</tally>
<tally id="6" name="tally 6">
<filters>6</filters>
<filters>7</filters>
<scores>current</scores>
</tally>
<tally id="7" name="tally 7">
<filters>4</filters>
<scores>total</scores>
</tally>
<tally id="8" name="tally 8">
<filters>8</filters>
<scores>current</scores>
</tally>
</tallies>

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

@ -1 +1 @@
cc306dfd3e3669827cb2f3bdc38005c0616f6faecb9277f3b0f869fdecf826ae662c65d67a402705956f100d2e13c7c11332364fa41ccc2e215ab100709b627f
35f04a6f062ef64116ef4eb0e9b803cd44cff7e185e2b53c9174afad8a26ca1a436ca9b800d6a228e006a9129f4536d7dce289d7a11cd56c6949d71d6a201b31

28
tests/regression_tests/filter_mesh/test.py
View file

@ -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')

2
tests/regression_tests/mg_tallies/inputs_true.dat
View file

@ -33,7 +33,7 @@
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" type="regular">
<mesh id="1">
<dimension>10 1 1</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>929.45 1000 1000</upper_right>

3
tests/regression_tests/mg_tallies/test.py
View file

@ -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]

2
tests/regression_tests/mgxs_library_mesh/inputs_true.dat
View file

@ -309,7 +309,7 @@
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" type="regular">
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-100.0 -100.0</lower_left>
<width>100.0 100.0</width>

3
tests/regression_tests/mgxs_library_mesh/test.py
View file

@ -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.]

1
tests/regression_tests/plot/plots.xml
View file

@ -14,6 +14,7 @@
<width>25 25</width>
<pixels>200 200</pixels>
<mask components="1 3" background="255 255 255" />
<meshlines meshtype="tally" id="2" linewidth="0" />
</plot>
<plot id="3" basis="yz" color_by="material">

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

@ -1 +1 @@
6b3eb36488b995b42233e6b7c0164cf6c92a1823b3c91985b97fd076f86c0aad93fbdb74e70026f5f12c61a6aee52a09588171a7fd839511ee7d9efc3952beb2
fd61f6a5de632f06a39e2a938480ba3dc314810655e684291e4b255bb8c142a16f10fec414c587ecd727fa559d3760d6e5b43f53fe86e90238a69d9c5698db8e

20
tests/regression_tests/plot/tallies.xml Normal file
View file

@ -0,0 +1,20 @@
<?xml version="1.0"?>
<tallies>
<mesh id="2" type="rectilinear">
<x_grid>-10 10</x_grid>
<y_grid>-10 10</y_grid>
<z_grid>-10 0 5 7.5 8.75 10</z_grid>
</mesh>
<filter id="1">
<type>mesh</type>
<bins>2</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>total</scores>
</tally>
</tallies>

2
tests/regression_tests/tallies/inputs_true.dat
View file

@ -309,7 +309,7 @@
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" type="regular">
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-182.07 -182.07</lower_left>
<upper_right>182.07 182.07</upper_right>

4
tests/regression_tests/tallies/test.py
View file

@ -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]

2
tests/regression_tests/tally_arithmetic/inputs_true.dat
View file

@ -309,7 +309,7 @@
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" type="regular">
<mesh id="1">
<dimension>2 2 2</dimension>
<lower_left>-160.0 -160.0 -183.0</lower_left>
<upper_right>160.0 160.0 183.0</upper_right>

3
tests/regression_tests/tally_arithmetic/test.py
View file

@ -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]

2
tests/regression_tests/tally_slice_merge/inputs_true.dat
View file

@ -310,7 +310,7 @@
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<!--mesh-->
<mesh id="1" type="regular">
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-50.0 -50.0</lower_left>
<upper_right>50.0 50.0</upper_right>

3
tests/regression_tests/tally_slice_merge/test.py
View file

@ -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.]

6
tests/unit_tests/test_capi.py
View file

@ -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)

8
tests/unit_tests/test_mesh_from_lattice.py
View file

@ -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)

6
tests/unit_tests/test_settings.py
View file

@ -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]