Merge pull request #1527 from pshriwise/unstructured_mesh

Unstructured Mesh - Initial implementation
This commit is contained in:
Paul Romano 2020-03-31 08:54:56 -05:00 committed by GitHub
commit 9b499b4096
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
41 changed files with 29470 additions and 83 deletions

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

Binary file not shown.

After

Width:  |  Height:  |  Size: 97 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 60 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 118 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 129 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 194 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 129 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 110 KiB

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View file

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

View file

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

View file

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

View file

@ -764,7 +764,10 @@ class MeshFilter(Filter):
def mesh(self, mesh):
cv.check_type('filter mesh', mesh, openmc.MeshBase)
self._mesh = mesh
self.bins = list(mesh.indices)
if isinstance(mesh, openmc.UnstructuredMesh):
self.bins = list(range(len(mesh.volumes)))
else:
self.bins = list(mesh.indices)
def can_merge(self, other):
# Mesh filters cannot have more than one bin

View file

@ -80,6 +80,8 @@ class MeshBase(IDManagerMixin, metaclass=ABCMeta):
return RegularMesh.from_hdf5(group)
elif mesh_type == 'rectilinear':
return RectilinearMesh.from_hdf5(group)
elif mesh_type == 'unstructured':
return UnstructuredMesh.from_hdf5(group)
else:
raise ValueError('Unrecognized mesh type: "' + mesh_type + '"')
@ -105,8 +107,8 @@ class RegularMesh(MeshBase):
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.
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.
@ -586,3 +588,128 @@ class RectilinearMesh(MeshBase):
subelement.text = ' '.join(map(str, self.z_grid))
return element
class UnstructuredMesh(MeshBase):
"""A 3D unstructured mesh
Parameters
----------
filename : str
Location of the unstructured mesh file
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
filename : str
Name of the file containing the unstructured mesh
volumes : Iterable of float
Volumes of the unstructured mesh elements
total_volume : float
Volume of the unstructured mesh in total
centroids : Iterable of tuple
An iterable of element centroid coordinates, e.g. [(0.0, 0.0, 0.0),
(1.0, 1.0, 1.0), ...]
"""
def __init__(self, filename, mesh_id=None, name=''):
super().__init__(mesh_id, name)
self.filename = filename
self._volumes = []
self._centroids = []
@property
def filename(self):
return self._filename
@filename.setter
def filename(self, filename):
cv.check_type('Unstructured Mesh filename', filename, str)
self._filename = filename
@property
def volumes(self):
return self._volumes
@volumes.setter
def volumes(self, volumes):
cv.check_type("Unstructured mesh volumes", volumes, Iterable, Real)
self._volumes = volumes
@property
def total_volume(self):
return np.sum(self.volumes)
@property
def centroids(self):
return self._centroids
@centroids.setter
def centroids(self, centroids):
cv.check_type("Unstructured mesh centroids", centroids,
Iterable, Real)
self._centroids = centroids
def __repr__(self):
string = super().__repr__()
return string + '{: <16}=\t{}\n'.format('\tFilename', self.filename)
@classmethod
def from_hdf5(cls, group):
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
filename = group['filename'][()].decode()
mesh = cls(filename, mesh_id=mesh_id)
vol_data = group['volumes'][()]
centroids = group['centroids'][()]
mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],))
mesh.centroids = np.reshape(centroids, (vol_data.shape[0], 3))
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", "unstructured")
subelement = ET.SubElement(element, "filename")
subelement.text = self.filename
return element
@classmethod
def from_xml_element(cls, elem):
"""Generate unstructured mesh object from XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.UnstructuredMesh
UnstructuredMesh generated from an XML element
"""
mesh_id = int(get_text(elem, 'id'))
filename = get_text(elem, 'filename')
mesh = cls(filename, mesh_id)
return mesh

View file

@ -9,11 +9,13 @@
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
#include <fmt/core.h> // for fmt
#include "xtensor/xbuilder.hpp"
#include "xtensor/xeval.hpp"
#include "xtensor/xmath.hpp"
#include "xtensor/xsort.hpp"
#include "xtensor/xtensor.hpp"
#include "xtensor/xview.hpp"
#include "openmc/capi.h"
#include "openmc/constants.h"
@ -21,6 +23,7 @@
#include "openmc/hdf5_interface.h"
#include "openmc/message_passing.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/tallies/filter.h"
#include "openmc/xml_interface.h"
@ -81,12 +84,30 @@ Mesh::Mesh(pugi::xml_node node)
}
}
//==============================================================================
// Structured Mesh implementation
//==============================================================================
std::string
StructuredMesh::bin_label(int bin) const {
std::vector<int> ijk(n_dimension_);
get_indices_from_bin(bin, ijk.data());
if (n_dimension_ > 2) {
return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]);
} else if (n_dimension_ > 1) {
return fmt::format("Mesh Index ({}, {})", ijk[0], ijk[1]);
} else {
return fmt::format("Mesh Index ({})", ijk[0]) ;
}
}
//==============================================================================
// RegularMesh implementation
//==============================================================================
RegularMesh::RegularMesh(pugi::xml_node node)
: Mesh {node}
: StructuredMesh {node}
{
// Determine number of dimensions for mesh
if (check_for_node(node, "dimension")) {
@ -792,11 +813,12 @@ void RegularMesh::to_hdf5(hid_t group) const
}
xt::xtensor<double, 1>
RegularMesh::count_sites(const Particle::Bank* bank, int64_t length,
bool* outside) const
RegularMesh::count_sites(const Particle::Bank* bank,
int64_t length,
bool* outside) const
{
// Determine shape of array for counts
std::size_t m = xt::prod(shape_)();
std::size_t m = this->n_bins();
std::vector<std::size_t> shape = {m};
// Create array of zeros
@ -851,7 +873,7 @@ RegularMesh::count_sites(const Particle::Bank* bank, int64_t length,
//==============================================================================
RectilinearMesh::RectilinearMesh(pugi::xml_node node)
: Mesh {node}
: StructuredMesh {node}
{
n_dimension_ = 3;
@ -1190,9 +1212,7 @@ void RectilinearMesh::get_indices_from_bin(int bin, int* ijk) const
int RectilinearMesh::n_bins() const
{
int n_bins = 1;
for (auto dim : shape_) n_bins *= dim;
return n_bins;
return xt::prod(shape_)();
}
int RectilinearMesh::n_surface_bins() const
@ -1383,6 +1403,10 @@ check_regular_mesh(int32_t index, RegularMesh** mesh)
// C API functions
//==============================================================================
RegularMesh* get_regular_mesh(int32_t index) {
return dynamic_cast<RegularMesh*>(model::meshes[index].get());
}
//! Extend the meshes array by n elements
extern "C" int
openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end)
@ -1450,7 +1474,6 @@ openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
mesh->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
mesh->n_dimension_ = mesh->shape_.size();
return 0;
}
@ -1502,6 +1525,550 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
return 0;
}
#ifdef DAGMC
UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
{
// unstructured always assumed to be 3D
n_dimension_ = 3;
// check the mesh type
if (check_for_node(node, "type")) {
auto temp = get_node_value(node, "type", true, true);
if (temp != "unstructured") {
fatal_error("Invalid mesh type: " + temp);
}
}
// get the filename of the unstructured mesh to load
if (check_for_node(node, "filename")) {
filename_ = get_node_value(node, "filename");
} else {
fatal_error("No filename supplied for unstructured mesh with ID: " +
std::to_string(id_));
}
// create MOAB instance
mbi_ = std::make_unique<moab::Core>();
// load unstructured mesh file
moab::ErrorCode rval = mbi_->load_file(filename_.c_str());
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to load the unstructured mesh file: " + filename_);
}
// set member range of tetrahedral entities
rval = mbi_->get_entities_by_dimension(0, n_dimension_, ehs_);
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to get all tetrahedral elements");
}
if (!ehs_.all_of_type(moab::MBTET)) {
warning("Non-tetrahedral elements found in unstructured "
"mesh file: " + filename_);
}
// make an entity set for all tetrahedra
// this is used for convenience later in output
rval = mbi_->create_meshset(moab::MESHSET_SET, tetset_);
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to create an entity set for the tetrahedral elements");
}
rval = mbi_->add_entities(tetset_, ehs_);
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to add tetrahedra to an entity set.");
}
// build acceleration data structures
compute_barycentric_data(ehs_);
build_kdtree(ehs_);
}
void
UnstructuredMesh::build_kdtree(const moab::Range& all_tets)
{
moab::Range all_tris;
int adj_dim = 2;
moab::ErrorCode rval = mbi_->get_adjacencies(all_tets,
adj_dim,
true,
all_tris,
moab::Interface::UNION);
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to get adjacent triangles for tets");
}
if (!all_tris.all_of_type(moab::MBTRI)) {
warning("Non-triangle elements found in tet adjacencies in "
"unstructured mesh file: " + filename_);
}
// combine into one range
moab::Range all_tets_and_tris;
all_tets_and_tris.merge(all_tets);
all_tets_and_tris.merge(all_tris);
// create a kd-tree instance
kdtree_ = std::make_unique<moab::AdaptiveKDTree>(mbi_.get());
// build the tree
rval = kdtree_->build_tree(all_tets_and_tris, &kdtree_root_);
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to construct KDTree for the "
"unstructured mesh file: " + filename_);
}
}
void
UnstructuredMesh::intersect_track(const moab::CartVect& start,
const moab::CartVect& dir,
double track_len,
std::vector<double>& hits) const {
hits.clear();
moab::ErrorCode rval;
std::vector<moab::EntityHandle> tris;
// get all intersections with triangles in the tet mesh
// (distances are relative to the start point, not the previous intersection)
rval = kdtree_->ray_intersect_triangles(kdtree_root_,
FP_COINCIDENT,
dir.array(),
start.array(),
tris,
hits,
0,
track_len);
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to compute intersections on unstructured mesh: " + filename_);
}
// remove duplicate intersection distances
std::unique(hits.begin(), hits.end());
// sorts by first component of std::pair by default
std::sort(hits.begin(), hits.end());
}
void
UnstructuredMesh::bins_crossed(const Particle* p,
std::vector<int>& bins,
std::vector<double>& lengths) const
{
moab::ErrorCode rval;
Position last_r{p->r_last_};
Position r{p->r()};
Direction u{p->u()};
u /= u.norm();
moab::CartVect r0(last_r.x, last_r.y, last_r.z);
moab::CartVect r1(r.x, r.y, r.z);
moab::CartVect dir(u.x, u.y, u.z);
double track_len = (r1 - r0).length();
r0 -= TINY_BIT * dir;
r1 += TINY_BIT * dir;
std::vector<double> hits;
intersect_track(r0, dir, track_len, hits);
bins.clear();
lengths.clear();
// if there are no intersections the track may lie entirely
// within a single tet. If this is the case, apply entire
// score to that tet and return.
if (hits.size() == 0) {
Position midpoint = last_r + u * (track_len * 0.5);
int bin = this->get_bin(midpoint);
if (bin != -1) {
bins.push_back(bin);
lengths.push_back(1.0);
}
return;
}
// for each segment in the set of tracks, try to look up a tet
// at the midpoint of the segment
Position current = last_r;
double last_dist = 0.0;
for (const auto& hit : hits) {
// get the segment length
double segment_length = hit - last_dist;
last_dist = hit;
// find the midpoint of this segment
Position midpoint = current + u * (segment_length * 0.5);
// try to find a tet for this position
int bin = this->get_bin(midpoint);
// determine the start point for this segment
current = last_r + u * hit;
if (bin == -1) {
continue;
}
bins.push_back(bin);
lengths.push_back(segment_length / track_len);
}
// tally remaining portion of track after last hit if
// the last segment of the track is in the mesh but doesn't
// reach the other side of the tet
if (hits.back() < track_len) {
Position segment_start = last_r + u * hits.back();
double segment_length = track_len - hits.back();
Position midpoint = segment_start + u * (segment_length * 0.5);
int bin = this->get_bin(midpoint);
if (bin != -1) {
bins.push_back(bin);
lengths.push_back(segment_length / track_len);
}
}
};
moab::EntityHandle
UnstructuredMesh::get_tet(const Position& r) const
{
moab::CartVect pos(r.x, r.y, r.z);
// find the leaf of the kd-tree for this position
moab::AdaptiveKDTreeIter kdtree_iter;
moab::ErrorCode rval = kdtree_->point_search(pos.array(), kdtree_iter);
if (rval != moab::MB_SUCCESS) { return 0; }
// retrieve the tet elements of this leaf
moab::EntityHandle leaf = kdtree_iter.handle();
moab::Range tets;
rval = mbi_->get_entities_by_dimension(leaf, 3, tets, false);
if (rval != moab::MB_SUCCESS) {
warning("MOAB error finding tets.");
}
// loop over the tets in this leaf, returning the containing tet if found
for (const auto& tet : tets) {
if (point_in_tet(pos, tet)) {
return tet;
}
}
// if no tet is found, return an invalid handle
return 0;
}
double UnstructuredMesh::tet_volume(moab::EntityHandle tet) const {
std::vector<moab::EntityHandle> conn;
moab::ErrorCode rval = mbi_->get_connectivity(&tet, 1, conn);
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to get tet connectivity");
}
moab::CartVect p[4];
rval = mbi_->get_coords(conn.data(), conn.size(), p[0].array());
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to get tet coords");
}
return 1.0 / 6.0 * (((p[1] - p[0]) * (p[2] - p[0])) % (p[3] - p[0]));
}
void UnstructuredMesh::surface_bins_crossed(const Particle* p, std::vector<int>& bins) const {
// TODO: Implement triangle crossings here
throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."};
}
int
UnstructuredMesh::get_bin(Position r) const {
moab::EntityHandle tet = get_tet(r);
if (tet == 0) {
return -1;
} else {
return get_bin_from_ent_handle(tet);
}
}
void
UnstructuredMesh::compute_barycentric_data(const moab::Range& tets) {
moab::ErrorCode rval;
baryc_data_.clear();
baryc_data_.resize(tets.size());
// compute the barycentric data for each tet element
// and store it as a 3x3 matrix
for (auto& tet : tets) {
std::vector<moab::EntityHandle> verts;
rval = mbi_->get_connectivity(&tet, 1, verts);
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to get connectivity of tet on umesh: " + filename_);
}
moab::CartVect p[4];
rval = mbi_->get_coords(verts.data(), verts.size(), p[0].array());
if (rval != moab::MB_SUCCESS) {
fatal_error("Failed to get coordinates of a tet in umesh: " + filename_);
}
moab::Matrix3 a(p[1] - p[0], p[2] - p[0], p[3] - p[0], true);
// invert now to avoid this cost later
a = a.transpose().inverse();
baryc_data_.at(get_bin_from_ent_handle(tet)) = a;
}
}
void
UnstructuredMesh::to_hdf5(hid_t group) const
{
hid_t mesh_group = create_group(group, fmt::format("mesh {}", id_));
write_dataset(mesh_group, "type", "unstructured");
write_dataset(mesh_group, "filename", filename_);
// write volume and centroid of each tet
std::vector<double> tet_vols;
xt::xtensor<double, 2> centroids({ehs_.size(), 3});
for (int i = 0; i < ehs_.size(); i++) {
const auto& eh = ehs_[i];
tet_vols.emplace_back(this->tet_volume(eh));
Position c = this->centroid(eh);
xt::view(centroids, i, xt::all()) = xt::xarray<double>({c.x, c.y, c.z});
}
write_dataset(mesh_group, "volumes", tet_vols);
write_dataset(mesh_group, "centroids", centroids);
close_group(mesh_group);
}
bool
UnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet) const {
moab::ErrorCode rval;
// get tet vertices
std::vector<moab::EntityHandle> verts;
rval = mbi_->get_connectivity(&tet, 1, verts);
if (rval != moab::MB_SUCCESS) {
warning("Failed to get vertices of tet in umesh: " + filename_);
return false;
}
// first vertex is used as a reference point for the barycentric data -
// retrieve its coordinates
moab::CartVect p_zero;
rval = mbi_->get_coords(verts.data(), 1, p_zero.array());
if (rval != moab::MB_SUCCESS) {
warning("Failed to get coordinates of a vertex in "
"unstructured mesh: " + filename_);
return false;
}
// look up barycentric data
int idx = get_bin_from_ent_handle(tet);
const moab::Matrix3& a_inv = baryc_data_[idx];
moab::CartVect bary_coords = a_inv * (r - p_zero);
return (bary_coords[0] >= 0.0 &&
bary_coords[1] >= 0.0 &&
bary_coords[2] >= 0.0 &&
bary_coords[0] + bary_coords[1] + bary_coords[2] <= 1.0);
}
int
UnstructuredMesh::get_bin_from_index(int idx) const {
if (idx >= n_bins()) {
fatal_error(fmt::format("Invalid bin index: {}", idx));
}
return ehs_[idx] - ehs_[0];
}
int
UnstructuredMesh::get_index(const Position& r,
bool* in_mesh) const {
int bin = get_bin(r);
*in_mesh = bin != -1;
return bin;
}
int UnstructuredMesh::get_index_from_bin(int bin) const {
return bin;
}
std::pair<std::vector<double>, std::vector<double>>
UnstructuredMesh::plot(Position plot_ll, Position plot_ur) const {
// TODO: Implement mesh lines
return {};
}
int
UnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const {
int bin = eh - ehs_[0];
if (bin >= n_bins()) {
fatal_error(fmt::format("Invalid bin: {}", bin));
}
return bin;
}
moab::EntityHandle
UnstructuredMesh::get_ent_handle_from_bin(int bin) const {
if (bin >= n_bins()) {
fatal_error(fmt::format("Invalid bin index: ", bin));
}
return ehs_[bin];
}
int UnstructuredMesh::n_bins() const {
return ehs_.size();
}
int UnstructuredMesh::n_surface_bins() const {
// collect all triangles in the set of tets for this mesh
moab::Range tris;
moab::ErrorCode rval;
rval = mbi_->get_entities_by_type(0, moab::MBTRI, tris);
if (rval != moab::MB_SUCCESS) {
warning("Failed to get all triangles in the mesh instance");
return -1;
}
return 2 * tris.size();
}
Position
UnstructuredMesh::centroid(moab::EntityHandle tet) const {
moab::ErrorCode rval;
// look up the tet connectivity
std::vector<moab::EntityHandle> conn;
rval = mbi_->get_connectivity(&tet, 1, conn);
if (rval != moab::MB_SUCCESS) {
warning("Failed to get connectivity of a mesh element.");
return {};
}
// get the coordinates
std::vector<moab::CartVect> coords(conn.size());
rval = mbi_->get_coords(conn.data(), conn.size(), coords[0].array());
if (rval != moab::MB_SUCCESS) {
warning("Failed to get the coordinates of a mesh element.");
return {};
}
// compute the centroid of the element vertices
moab::CartVect centroid(0.0, 0.0, 0.0);
for(const auto& coord : coords) {
centroid += coord;
}
centroid /= double(coords.size());
return {centroid[0], centroid[1], centroid[2]};
}
std::string
UnstructuredMesh::bin_label(int bin) const {
return fmt::format("Mesh Index ({})", bin);
};
std::pair<moab::Tag, moab::Tag>
UnstructuredMesh::get_score_tags(std::string score) const {
moab::ErrorCode rval;
// add a tag to the mesh
// all scores are treated as a single value
// with an uncertainty
moab::Tag value_tag;
// create the value tag if not present and get handle
double default_val = 0.0;
auto val_string = score + "_mean";
rval = mbi_->tag_get_handle(val_string.c_str(),
1,
moab::MB_TYPE_DOUBLE,
value_tag,
moab::MB_TAG_DENSE|moab::MB_TAG_CREAT,
&default_val);
if (rval != moab::MB_SUCCESS) {
auto msg = fmt::format("Could not create or retrieve the value tag for the score {}"
" on unstructured mesh {}", score, id_);
fatal_error(msg);
}
// create the std dev tag if not present and get handle
moab::Tag error_tag;
std::string err_string = score + "_std_dev";
rval = mbi_->tag_get_handle(err_string.c_str(),
1,
moab::MB_TYPE_DOUBLE,
error_tag,
moab::MB_TAG_DENSE|moab::MB_TAG_CREAT,
&default_val);
if (rval != moab::MB_SUCCESS) {
auto msg = fmt::format("Could not create or retrieve the error tag for the score {}"
" on unstructured mesh {}", score, id_);
fatal_error(msg);
}
// return the populated tag handles
return {value_tag, error_tag};
}
void
UnstructuredMesh::add_score(std::string score) const {
auto score_tags = this->get_score_tags(score);
}
void
UnstructuredMesh::set_score_data(const std::string& score,
std::vector<double> values,
std::vector<double> std_dev) const {
auto score_tags = this->get_score_tags(score);
// normalize tally values by element volume
for (int i = 0; i < ehs_.size(); i++) {
auto eh = this->get_ent_handle_from_bin(i);
double volume = this->tet_volume(eh);
values[i] /= volume;
std_dev[i] /= volume;
}
moab::ErrorCode rval;
// set the score value
rval = mbi_->tag_set_data(score_tags.first, ehs_, values.data());
if (rval != moab::MB_SUCCESS) {
auto msg = fmt::format("Failed to set the tally value for score '{}' "
"on unstructured mesh {}", score, id_);
warning(msg);
}
// set the error value
rval = mbi_->tag_set_data(score_tags.second, ehs_, std_dev.data());
if (rval != moab::MB_SUCCESS) {
auto msg = fmt::format("Failed to set the tally error for score '{}' "
"on unstructured mesh {}", score, id_);
warning(msg);
}
}
void
UnstructuredMesh::write(std::string base_filename) const {
// add extension to the base name
auto filename = base_filename + ".vtk";
write_message("Writing unstructured mesh " + filename + "...", 5);
filename = settings::path_output + filename;
// write the tetrahedral elements of the mesh only
// to avoid clutter from zero-value data on other
// elements during visualization
moab::ErrorCode rval;
rval = mbi_->write_mesh(filename.c_str(), &tetset_, 1);
if (rval != moab::MB_SUCCESS) {
auto msg = fmt::format("Failed to write unstructured mesh {}", id_);
warning(msg);
}
}
#endif
//==============================================================================
// Non-member functions
//==============================================================================
@ -1521,6 +2088,13 @@ void read_meshes(pugi::xml_node root)
model::meshes.push_back(std::make_unique<RegularMesh>(node));
} else if (mesh_type == "rectilinear") {
model::meshes.push_back(std::make_unique<RectilinearMesh>(node));
#ifdef DAGMC
} else if (mesh_type == "unstructured") {
model::meshes.push_back(std::make_unique<UnstructuredMesh>(node));
#else
} else if (mesh_type == "unstructured") {
fatal_error("Unstructured mesh support is disabled.");
#endif
} else {
fatal_error("Invalid mesh type: " + mesh_type);
}

View file

@ -24,6 +24,7 @@
#include "openmc/simulation.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/filter_mesh.h"
#include "openmc/tallies/tally.h"
#include "openmc/timer.h"
@ -164,6 +165,11 @@ openmc_statepoint_write(const char* filename, bool* write_source)
tally_ids.push_back(tally->id_);
write_attribute(tallies_group, "ids", tally_ids);
#ifdef DAGMC
// write unstructured mesh tallies to VTK if possible
write_unstructured_mesh_results();
#endif
// Write all tally information except results
for (const auto& tally : model::tallies) {
hid_t tally_group = create_group(tallies_group,
@ -677,6 +683,82 @@ void read_source_bank(hid_t group_id)
H5Tclose(banktype);
}
#ifdef DAGMC
void write_unstructured_mesh_results() {
for (auto& tally : model::tallies) {
for (auto filter_idx : tally->filters()) {
auto& filter = model::tally_filters[filter_idx];
if (filter->type() != "mesh") continue;
// check if the filter uses an unstructured mesh
auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
auto mesh_idx = mesh_filter->mesh();
auto umesh = dynamic_cast<UnstructuredMesh*>(model::meshes[mesh_idx].get());
if (!umesh) continue;
// if this tally has more than one filter, print
// warning and skip writing the mesh
if (tally->filters().size() > 1) {
warning(fmt::format("Skipping unstructured mesh writing for tally "
"{}. More than one filter is present on the tally.",
tally->id_));
break;
}
int n_realizations = tally->n_realizations_;
// write each score/nuclide combination for this tally
for (int i_score = 0; i_score < tally->scores_.size(); i_score++) {
for (int i_nuc = 0; i_nuc < tally->nuclides_.size(); i_nuc++) {
// index for this nuclide and score
int nuc_score_idx = i_score + i_nuc*tally->scores_.size();
// construct result vectors
std::vector<double> mean_vec, std_dev_vec;
for (int j = 0; j < tally->results_.shape()[0]; j++) {
// mean
double mean = tally->results_(j, nuc_score_idx, TallyResult::SUM) / n_realizations;
mean_vec.push_back(mean);
// std. dev.
double sum_sq = tally->results_(j , nuc_score_idx, TallyResult::SUM_SQ);
if (n_realizations > 1) {
double std_dev = sum_sq/n_realizations - mean*mean;
std_dev = std::sqrt(std_dev / (n_realizations - 1));
std_dev_vec.push_back(std_dev);
} else {
std_dev_vec.push_back(0.0);
}
}
// generate a name for the value
std::string nuclide_name = "total"; // start with total by default
if (tally->nuclides_[i_nuc] > -1) {
nuclide_name = data::nuclides[tally->nuclides_[i_nuc]]->name_;
}
std::string score_name = tally->score_name(i_score);
auto score_str = fmt::format("{}_{}", score_name, nuclide_name);
umesh->set_score_data(score_str, mean_vec, std_dev_vec);
}
}
// Generate a file name based on the tally id
// and the current batch number
int w = std::to_string(settings::n_max_batches).size();
std::string filename = fmt::format("tally_{0}.{1:0{2}}",
tally->id_,
simulation::current_batch,
w);
// Write the unstructured mesh and data to file
umesh->write(filename);
}
}
}
#endif
void write_tally_results_nr(hid_t file_id)
{
// ==========================================================================

View file

@ -55,18 +55,7 @@ std::string
MeshFilter::text_label(int bin) const
{
auto& mesh = *model::meshes[mesh_];
int n_dim = mesh.n_dimension_;
std::vector<int> ijk(n_dim);
mesh.get_indices_from_bin(bin, ijk.data());
if (n_dim > 2) {
return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]);
} else if (n_dim > 1) {
return fmt::format("Mesh Index ({}, {})", ijk[0], ijk[1]);
} else {
return fmt::format("Mesh Index ({})", ijk[0]) ;
}
return mesh.bin_label(bin);
}
void

View file

@ -833,6 +833,14 @@ void Tally::accumulate()
}
}
std::string
Tally::score_name(int score_idx) const {
if (score_idx < 0 || score_idx >= scores_.size()) {
fatal_error("Index in scores array is out of bounds.");
}
return reaction_name(scores_[score_idx]);
}
//==============================================================================
// Non-member functions
//==============================================================================

View file

@ -0,0 +1,92 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" name="fuel" region="1 -2 3 -4 5 -6" universe="1" />
<cell id="2" material="2" name="clad" region="(-1 | 2 | -3 | 4 | -5 | 6) (7 -8 9 -10 11 -12)" universe="1" />
<cell id="3" material="3" name="water" region="(-7 | 8 | -9 | 10 | -11 | 12) (13 -14 15 -16 17 -18)" universe="1" />
<surface coeffs="-5.0" id="1" name="minimum x" type="x-plane" />
<surface coeffs="5.0" id="2" name="maximum x" type="x-plane" />
<surface coeffs="-5.0" id="3" name="minimum y" type="y-plane" />
<surface coeffs="5.0" id="4" name="maximum y" type="y-plane" />
<surface coeffs="-5.0" id="5" name="minimum z" type="z-plane" />
<surface coeffs="5.0" id="6" name="maximum z" type="z-plane" />
<surface coeffs="-6.0" id="7" name="minimum x" type="x-plane" />
<surface coeffs="6.0" id="8" name="maximum x" type="x-plane" />
<surface coeffs="-6.0" id="9" name="minimum y" type="y-plane" />
<surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
<surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
<surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
<surface boundary="vacuum" coeffs="-15" id="13" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="15" id="14" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-15" id="15" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="15" id="16" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-15" id="17" name="minimum z" type="z-plane" />
<surface boundary="vacuum" coeffs="15" id="18" name="maximum z" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="1" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="2" name="zircaloy">
<density units="g/cc" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="3" name="water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>
</space>
<angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>15000000.0 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1">
<dimension>10 10 10</dimension>
<lower_left>-10.0 -10.0 -10.0</lower_left>
<upper_right>10.0 10.0 10.0</upper_right>
</mesh>
<mesh id="2" type="unstructured">
<filename>test_mesh_tets_w_holes.h5m</filename>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="2" type="mesh">
<bins>2</bins>
</filter>
<tally id="1" name="regular mesh tally">
<filters>1</filters>
<scores>flux</scores>
<estimator>collision</estimator>
</tally>
<tally id="2" name="unstructured mesh tally">
<filters>2</filters>
<scores>flux</scores>
<estimator>collision</estimator>
</tally>
</tallies>

View file

@ -0,0 +1,92 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="4" material="4" name="fuel" region="19 -20 21 -22 23 -24" universe="2" />
<cell id="5" material="5" name="clad" region="(-19 | 20 | -21 | 22 | -23 | 24) (25 -26 27 -28 29 -30)" universe="2" />
<cell id="6" material="6" name="water" region="(-25 | 26 | -27 | 28 | -29 | 30) (31 -32 33 -34 35 -36)" universe="2" />
<surface coeffs="-5.0" id="19" name="minimum x" type="x-plane" />
<surface coeffs="5.0" id="20" name="maximum x" type="x-plane" />
<surface coeffs="-5.0" id="21" name="minimum y" type="y-plane" />
<surface coeffs="5.0" id="22" name="maximum y" type="y-plane" />
<surface coeffs="-5.0" id="23" name="minimum z" type="z-plane" />
<surface coeffs="5.0" id="24" name="maximum z" type="z-plane" />
<surface coeffs="-6.0" id="25" name="minimum x" type="x-plane" />
<surface coeffs="6.0" id="26" name="maximum x" type="x-plane" />
<surface coeffs="-6.0" id="27" name="minimum y" type="y-plane" />
<surface coeffs="6.0" id="28" name="maximum y" type="y-plane" />
<surface coeffs="-6.0" id="29" name="minimum z" type="z-plane" />
<surface coeffs="6.0" id="30" name="maximum z" type="z-plane" />
<surface boundary="vacuum" coeffs="-15" id="31" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="15" id="32" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-15" id="33" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="15" id="34" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-15" id="35" name="minimum z" type="z-plane" />
<surface boundary="vacuum" coeffs="15" id="36" name="maximum z" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="4" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="5" name="zircaloy">
<density units="g/cc" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="6" name="water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>
</space>
<angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>15000000.0 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="3">
<dimension>10 10 10</dimension>
<lower_left>-10.0 -10.0 -10.0</lower_left>
<upper_right>10.0 10.0 10.0</upper_right>
</mesh>
<mesh id="4" type="unstructured">
<filename>test_mesh_tets.h5m</filename>
</mesh>
<filter id="3" type="mesh">
<bins>3</bins>
</filter>
<filter id="4" type="mesh">
<bins>4</bins>
</filter>
<tally id="3" name="regular mesh tally">
<filters>3</filters>
<scores>flux</scores>
<estimator>collision</estimator>
</tally>
<tally id="4" name="unstructured mesh tally">
<filters>4</filters>
<scores>flux</scores>
<estimator>collision</estimator>
</tally>
</tallies>

View file

@ -0,0 +1,92 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="7" material="7" name="fuel" region="37 -38 39 -40 41 -42" universe="3" />
<cell id="8" material="8" name="clad" region="(-37 | 38 | -39 | 40 | -41 | 42) (43 -44 45 -46 47 -48)" universe="3" />
<cell id="9" material="9" name="water" region="(-43 | 44 | -45 | 46 | -47 | 48) (49 -50 51 -52 53 -54)" universe="3" />
<surface coeffs="-5.0" id="37" name="minimum x" type="x-plane" />
<surface coeffs="5.0" id="38" name="maximum x" type="x-plane" />
<surface coeffs="-5.0" id="39" name="minimum y" type="y-plane" />
<surface coeffs="5.0" id="40" name="maximum y" type="y-plane" />
<surface coeffs="-5.0" id="41" name="minimum z" type="z-plane" />
<surface coeffs="5.0" id="42" name="maximum z" type="z-plane" />
<surface coeffs="-6.0" id="43" name="minimum x" type="x-plane" />
<surface coeffs="6.0" id="44" name="maximum x" type="x-plane" />
<surface coeffs="-6.0" id="45" name="minimum y" type="y-plane" />
<surface coeffs="6.0" id="46" name="maximum y" type="y-plane" />
<surface coeffs="-6.0" id="47" name="minimum z" type="z-plane" />
<surface coeffs="6.0" id="48" name="maximum z" type="z-plane" />
<surface boundary="vacuum" coeffs="-10" id="49" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="10" id="50" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-10" id="51" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="10" id="52" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-10" id="53" name="minimum z" type="z-plane" />
<surface boundary="vacuum" coeffs="10" id="54" name="maximum z" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="7" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="8" name="zircaloy">
<density units="g/cc" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="9" name="water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>
</space>
<angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>15000000.0 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="5">
<dimension>10 10 10</dimension>
<lower_left>-10.0 -10.0 -10.0</lower_left>
<upper_right>10.0 10.0 10.0</upper_right>
</mesh>
<mesh id="6" type="unstructured">
<filename>test_mesh_tets_w_holes.h5m</filename>
</mesh>
<filter id="5" type="mesh">
<bins>5</bins>
</filter>
<filter id="6" type="mesh">
<bins>6</bins>
</filter>
<tally id="5" name="regular mesh tally">
<filters>5</filters>
<scores>flux</scores>
<estimator>collision</estimator>
</tally>
<tally id="6" name="unstructured mesh tally">
<filters>6</filters>
<scores>flux</scores>
<estimator>collision</estimator>
</tally>
</tallies>

View file

@ -0,0 +1,92 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="10" material="10" name="fuel" region="55 -56 57 -58 59 -60" universe="4" />
<cell id="11" material="11" name="clad" region="(-55 | 56 | -57 | 58 | -59 | 60) (61 -62 63 -64 65 -66)" universe="4" />
<cell id="12" material="12" name="water" region="(-61 | 62 | -63 | 64 | -65 | 66) (67 -68 69 -70 71 -72)" universe="4" />
<surface coeffs="-5.0" id="55" name="minimum x" type="x-plane" />
<surface coeffs="5.0" id="56" name="maximum x" type="x-plane" />
<surface coeffs="-5.0" id="57" name="minimum y" type="y-plane" />
<surface coeffs="5.0" id="58" name="maximum y" type="y-plane" />
<surface coeffs="-5.0" id="59" name="minimum z" type="z-plane" />
<surface coeffs="5.0" id="60" name="maximum z" type="z-plane" />
<surface coeffs="-6.0" id="61" name="minimum x" type="x-plane" />
<surface coeffs="6.0" id="62" name="maximum x" type="x-plane" />
<surface coeffs="-6.0" id="63" name="minimum y" type="y-plane" />
<surface coeffs="6.0" id="64" name="maximum y" type="y-plane" />
<surface coeffs="-6.0" id="65" name="minimum z" type="z-plane" />
<surface coeffs="6.0" id="66" name="maximum z" type="z-plane" />
<surface boundary="vacuum" coeffs="-10" id="67" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="10" id="68" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-10" id="69" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="10" id="70" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-10" id="71" name="minimum z" type="z-plane" />
<surface boundary="vacuum" coeffs="10" id="72" name="maximum z" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="10" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="11" name="zircaloy">
<density units="g/cc" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="12" name="water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>
</space>
<angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>15000000.0 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="7">
<dimension>10 10 10</dimension>
<lower_left>-10.0 -10.0 -10.0</lower_left>
<upper_right>10.0 10.0 10.0</upper_right>
</mesh>
<mesh id="8" type="unstructured">
<filename>test_mesh_tets.h5m</filename>
</mesh>
<filter id="7" type="mesh">
<bins>7</bins>
</filter>
<filter id="8" type="mesh">
<bins>8</bins>
</filter>
<tally id="7" name="regular mesh tally">
<filters>7</filters>
<scores>flux</scores>
<estimator>collision</estimator>
</tally>
<tally id="8" name="unstructured mesh tally">
<filters>8</filters>
<scores>flux</scores>
<estimator>collision</estimator>
</tally>
</tallies>

View file

@ -0,0 +1,92 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="13" material="13" name="fuel" region="73 -74 75 -76 77 -78" universe="5" />
<cell id="14" material="14" name="clad" region="(-73 | 74 | -75 | 76 | -77 | 78) (79 -80 81 -82 83 -84)" universe="5" />
<cell id="15" material="15" name="water" region="(-79 | 80 | -81 | 82 | -83 | 84) (85 -86 87 -88 89 -90)" universe="5" />
<surface coeffs="-5.0" id="73" name="minimum x" type="x-plane" />
<surface coeffs="5.0" id="74" name="maximum x" type="x-plane" />
<surface coeffs="-5.0" id="75" name="minimum y" type="y-plane" />
<surface coeffs="5.0" id="76" name="maximum y" type="y-plane" />
<surface coeffs="-5.0" id="77" name="minimum z" type="z-plane" />
<surface coeffs="5.0" id="78" name="maximum z" type="z-plane" />
<surface coeffs="-6.0" id="79" name="minimum x" type="x-plane" />
<surface coeffs="6.0" id="80" name="maximum x" type="x-plane" />
<surface coeffs="-6.0" id="81" name="minimum y" type="y-plane" />
<surface coeffs="6.0" id="82" name="maximum y" type="y-plane" />
<surface coeffs="-6.0" id="83" name="minimum z" type="z-plane" />
<surface coeffs="6.0" id="84" name="maximum z" type="z-plane" />
<surface boundary="vacuum" coeffs="-15" id="85" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="15" id="86" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-15" id="87" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="15" id="88" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-15" id="89" name="minimum z" type="z-plane" />
<surface boundary="vacuum" coeffs="15" id="90" name="maximum z" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="13" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="14" name="zircaloy">
<density units="g/cc" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="15" name="water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>
</space>
<angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>15000000.0 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="9">
<dimension>10 10 10</dimension>
<lower_left>-10.0 -10.0 -10.0</lower_left>
<upper_right>10.0 10.0 10.0</upper_right>
</mesh>
<mesh id="10" type="unstructured">
<filename>test_mesh_tets_w_holes.h5m</filename>
</mesh>
<filter id="9" type="mesh">
<bins>9</bins>
</filter>
<filter id="10" type="mesh">
<bins>10</bins>
</filter>
<tally id="9" name="regular mesh tally">
<filters>9</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="10" name="unstructured mesh tally">
<filters>10</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
</tallies>

View file

@ -0,0 +1,92 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="16" material="16" name="fuel" region="91 -92 93 -94 95 -96" universe="6" />
<cell id="17" material="17" name="clad" region="(-91 | 92 | -93 | 94 | -95 | 96) (97 -98 99 -100 101 -102)" universe="6" />
<cell id="18" material="18" name="water" region="(-97 | 98 | -99 | 100 | -101 | 102) (103 -104 105 -106 107 -108)" universe="6" />
<surface coeffs="-5.0" id="91" name="minimum x" type="x-plane" />
<surface coeffs="5.0" id="92" name="maximum x" type="x-plane" />
<surface coeffs="-5.0" id="93" name="minimum y" type="y-plane" />
<surface coeffs="5.0" id="94" name="maximum y" type="y-plane" />
<surface coeffs="-5.0" id="95" name="minimum z" type="z-plane" />
<surface coeffs="5.0" id="96" name="maximum z" type="z-plane" />
<surface coeffs="-6.0" id="97" name="minimum x" type="x-plane" />
<surface coeffs="6.0" id="98" name="maximum x" type="x-plane" />
<surface coeffs="-6.0" id="99" name="minimum y" type="y-plane" />
<surface coeffs="6.0" id="100" name="maximum y" type="y-plane" />
<surface coeffs="-6.0" id="101" name="minimum z" type="z-plane" />
<surface coeffs="6.0" id="102" name="maximum z" type="z-plane" />
<surface boundary="vacuum" coeffs="-15" id="103" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="15" id="104" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-15" id="105" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="15" id="106" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-15" id="107" name="minimum z" type="z-plane" />
<surface boundary="vacuum" coeffs="15" id="108" name="maximum z" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="16" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="17" name="zircaloy">
<density units="g/cc" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="18" name="water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>
</space>
<angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>15000000.0 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="11">
<dimension>10 10 10</dimension>
<lower_left>-10.0 -10.0 -10.0</lower_left>
<upper_right>10.0 10.0 10.0</upper_right>
</mesh>
<mesh id="12" type="unstructured">
<filename>test_mesh_tets.h5m</filename>
</mesh>
<filter id="11" type="mesh">
<bins>11</bins>
</filter>
<filter id="12" type="mesh">
<bins>12</bins>
</filter>
<tally id="11" name="regular mesh tally">
<filters>11</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="12" name="unstructured mesh tally">
<filters>12</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
</tallies>

View file

@ -0,0 +1,92 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="19" material="19" name="fuel" region="109 -110 111 -112 113 -114" universe="7" />
<cell id="20" material="20" name="clad" region="(-109 | 110 | -111 | 112 | -113 | 114) (115 -116 117 -118 119 -120)" universe="7" />
<cell id="21" material="21" name="water" region="(-115 | 116 | -117 | 118 | -119 | 120) (121 -122 123 -124 125 -126)" universe="7" />
<surface coeffs="-5.0" id="109" name="minimum x" type="x-plane" />
<surface coeffs="5.0" id="110" name="maximum x" type="x-plane" />
<surface coeffs="-5.0" id="111" name="minimum y" type="y-plane" />
<surface coeffs="5.0" id="112" name="maximum y" type="y-plane" />
<surface coeffs="-5.0" id="113" name="minimum z" type="z-plane" />
<surface coeffs="5.0" id="114" name="maximum z" type="z-plane" />
<surface coeffs="-6.0" id="115" name="minimum x" type="x-plane" />
<surface coeffs="6.0" id="116" name="maximum x" type="x-plane" />
<surface coeffs="-6.0" id="117" name="minimum y" type="y-plane" />
<surface coeffs="6.0" id="118" name="maximum y" type="y-plane" />
<surface coeffs="-6.0" id="119" name="minimum z" type="z-plane" />
<surface coeffs="6.0" id="120" name="maximum z" type="z-plane" />
<surface boundary="vacuum" coeffs="-10" id="121" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="10" id="122" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-10" id="123" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="10" id="124" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-10" id="125" name="minimum z" type="z-plane" />
<surface boundary="vacuum" coeffs="10" id="126" name="maximum z" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="19" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="20" name="zircaloy">
<density units="g/cc" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="21" name="water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>
</space>
<angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>15000000.0 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="13">
<dimension>10 10 10</dimension>
<lower_left>-10.0 -10.0 -10.0</lower_left>
<upper_right>10.0 10.0 10.0</upper_right>
</mesh>
<mesh id="14" type="unstructured">
<filename>test_mesh_tets_w_holes.h5m</filename>
</mesh>
<filter id="13" type="mesh">
<bins>13</bins>
</filter>
<filter id="14" type="mesh">
<bins>14</bins>
</filter>
<tally id="13" name="regular mesh tally">
<filters>13</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="14" name="unstructured mesh tally">
<filters>14</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
</tallies>

View file

@ -0,0 +1,92 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="22" material="22" name="fuel" region="127 -128 129 -130 131 -132" universe="8" />
<cell id="23" material="23" name="clad" region="(-127 | 128 | -129 | 130 | -131 | 132) (133 -134 135 -136 137 -138)" universe="8" />
<cell id="24" material="24" name="water" region="(-133 | 134 | -135 | 136 | -137 | 138) (139 -140 141 -142 143 -144)" universe="8" />
<surface coeffs="-5.0" id="127" name="minimum x" type="x-plane" />
<surface coeffs="5.0" id="128" name="maximum x" type="x-plane" />
<surface coeffs="-5.0" id="129" name="minimum y" type="y-plane" />
<surface coeffs="5.0" id="130" name="maximum y" type="y-plane" />
<surface coeffs="-5.0" id="131" name="minimum z" type="z-plane" />
<surface coeffs="5.0" id="132" name="maximum z" type="z-plane" />
<surface coeffs="-6.0" id="133" name="minimum x" type="x-plane" />
<surface coeffs="6.0" id="134" name="maximum x" type="x-plane" />
<surface coeffs="-6.0" id="135" name="minimum y" type="y-plane" />
<surface coeffs="6.0" id="136" name="maximum y" type="y-plane" />
<surface coeffs="-6.0" id="137" name="minimum z" type="z-plane" />
<surface coeffs="6.0" id="138" name="maximum z" type="z-plane" />
<surface boundary="vacuum" coeffs="-10" id="139" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="10" id="140" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-10" id="141" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="10" id="142" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-10" id="143" name="minimum z" type="z-plane" />
<surface boundary="vacuum" coeffs="10" id="144" name="maximum z" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="22" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="23" name="zircaloy">
<density units="g/cc" value="5.77" />
<nuclide ao="0.5145" name="Zr90" />
<nuclide ao="0.1122" name="Zr91" />
<nuclide ao="0.1715" name="Zr92" />
<nuclide ao="0.1738" name="Zr94" />
<nuclide ao="0.028" name="Zr96" />
</material>
<material id="24" name="water">
<density units="atom/b-cm" value="0.07416" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100</particles>
<batches>10</batches>
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>
</space>
<angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>15000000.0 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="15">
<dimension>10 10 10</dimension>
<lower_left>-10.0 -10.0 -10.0</lower_left>
<upper_right>10.0 10.0 10.0</upper_right>
</mesh>
<mesh id="16" type="unstructured">
<filename>test_mesh_tets.h5m</filename>
</mesh>
<filter id="15" type="mesh">
<bins>15</bins>
</filter>
<filter id="16" type="mesh">
<bins>16</bins>
</filter>
<tally id="15" name="regular mesh tally">
<filters>15</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="16" name="unstructured mesh tally">
<filters>16</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
</tallies>

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,237 @@
import glob
from itertools import product
import os
import openmc
import openmc.lib
import numpy as np
import pytest
from tests.testing_harness import PyAPITestHarness
pytestmark = pytest.mark.skipif(
not openmc.lib._dagmc_enabled(),
reason="Mesh library is not enabled.")
TETS_PER_VOXEL = 12
class UnstructuredMeshTest(PyAPITestHarness):
def __init__(self, statepoint_name, model, inputs_true, holes):
super().__init__(statepoint_name, model, inputs_true)
self.holes = holes # holes in the test mesh
def _compare_results(self):
with openmc.StatePoint(self._sp_name) as sp:
# loop over the tallies and get data
for tally in sp.tallies.values():
# find the regular and unstructured meshes
if tally.contains_filter(openmc.MeshFilter):
flt = tally.find_filter(openmc.MeshFilter)
if isinstance(flt.mesh, openmc.RegularMesh):
reg_mesh_data, reg_mesh_std_dev = self.get_mesh_tally_data(tally)
if self.holes:
reg_mesh_data = np.delete(reg_mesh_data, self.holes)
reg_mesh_std_dev = np.delete(reg_mesh_std_dev, self.holes)
else:
umesh_tally = tally
unstructured_data, unstructured_std_dev = self.get_mesh_tally_data(tally, True)
# we expect these results to be the same to within at least ten
# decimal places
decimals = 10 if umesh_tally.estimator == 'collision' else 8
np.testing.assert_array_almost_equal(unstructured_data,
reg_mesh_data,
decimals)
@staticmethod
def get_mesh_tally_data(tally, structured=False):
data = tally.get_reshaped_data(value='mean')
std_dev = tally.get_reshaped_data(value='std_dev')
if structured:
data.shape = (data.size // TETS_PER_VOXEL, TETS_PER_VOXEL)
std_dev.shape = (std_dev.size // TETS_PER_VOXEL, TETS_PER_VOXEL)
else:
data.shape = (data.size, 1)
std_dev.shape = (std_dev.size, 1)
return np.sum(data, axis=1), np.sum(std_dev, axis=1)
def _cleanup(self):
super()._cleanup()
output = glob.glob('tally*.vtk')
for f in output:
if os.path.exists(f):
os.remove(f)
param_values = (['collision', 'tracklength'], # estimators
[True, False], # geometry outside of the mesh
[(333, 90, 77), None]) # location of holes in the mesh
test_cases = []
for i, (estimator, ext_geom, holes) in enumerate(product(*param_values)):
test_cases.append({'estimator' : estimator,
'external_geom' : ext_geom,
'holes' : holes,
'inputs_true' : 'inputs_true{}.dat'.format(i)})
@pytest.mark.parametrize("test_opts", test_cases)
def test_unstructured_mesh(test_opts):
### Materials ###
materials = openmc.Materials()
fuel_mat = openmc.Material(name="fuel")
fuel_mat.add_nuclide("U235", 1.0)
fuel_mat.set_density('g/cc', 4.5)
materials.append(fuel_mat)
zirc_mat = openmc.Material(name="zircaloy")
zirc_mat.add_element("Zr", 1.0)
zirc_mat.set_density("g/cc", 5.77)
materials.append(zirc_mat)
water_mat = openmc.Material(name="water")
water_mat.add_nuclide("H1", 2.0)
water_mat.add_nuclide("O16", 1.0)
water_mat.set_density("atom/b-cm", 0.07416)
materials.append(water_mat)
materials.export_to_xml()
### Geometry ###
fuel_min_x = openmc.XPlane(-5.0, name="minimum x")
fuel_max_x = openmc.XPlane(5.0, name="maximum x")
fuel_min_y = openmc.YPlane(-5.0, name="minimum y")
fuel_max_y = openmc.YPlane(5.0, name="maximum y")
fuel_min_z = openmc.ZPlane(-5.0, name="minimum z")
fuel_max_z = openmc.ZPlane(5.0, name="maximum z")
fuel_cell = openmc.Cell(name="fuel")
fuel_cell.region = +fuel_min_x & -fuel_max_x & \
+fuel_min_y & -fuel_max_y & \
+fuel_min_z & -fuel_max_z
fuel_cell.fill = fuel_mat
clad_min_x = openmc.XPlane(-6.0, name="minimum x")
clad_max_x = openmc.XPlane(6.0, name="maximum x")
clad_min_y = openmc.YPlane(-6.0, name="minimum y")
clad_max_y = openmc.YPlane(6.0, name="maximum y")
clad_min_z = openmc.ZPlane(-6.0, name="minimum z")
clad_max_z = openmc.ZPlane(6.0, name="maximum z")
clad_cell = openmc.Cell(name="clad")
clad_cell.region = (-fuel_min_x | +fuel_max_x |
-fuel_min_y | +fuel_max_y |
-fuel_min_z | +fuel_max_z) & \
(+clad_min_x & -clad_max_x &
+clad_min_y & -clad_max_y &
+clad_min_z & -clad_max_z)
clad_cell.fill = zirc_mat
if test_opts['external_geom']:
bounds = (15, 15, 15)
else:
bounds = (10, 10, 10)
water_min_x = openmc.XPlane(x0=-bounds[0],
name="minimum x",
boundary_type='vacuum')
water_max_x = openmc.XPlane(x0=bounds[0],
name="maximum x",
boundary_type='vacuum')
water_min_y = openmc.YPlane(y0=-bounds[1],
name="minimum y",
boundary_type='vacuum')
water_max_y = openmc.YPlane(y0=bounds[1],
name="maximum y",
boundary_type='vacuum')
water_min_z = openmc.ZPlane(z0=-bounds[2],
name="minimum z",
boundary_type='vacuum')
water_max_z = openmc.ZPlane(z0=bounds[2],
name="maximum z",
boundary_type='vacuum')
water_cell = openmc.Cell(name="water")
water_cell.region = (-clad_min_x | +clad_max_x |
-clad_min_y | +clad_max_y |
-clad_min_z | +clad_max_z) & \
(+water_min_x & -water_max_x &
+water_min_y & -water_max_y &
+water_min_z & -water_max_z)
water_cell.fill = water_mat
# create a containing universe
geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell])
### Tallies ###
# create meshes
regular_mesh = openmc.RegularMesh()
regular_mesh.dimension = (10, 10, 10)
regular_mesh.lower_left = (-10.0, -10.0, -10.0)
regular_mesh.upper_right = (10.0, 10.0, 10.0)
regular_mesh_filter = openmc.MeshFilter(mesh=regular_mesh)
if test_opts['holes']:
mesh_filename = "test_mesh_tets_w_holes.h5m"
else:
mesh_filename = "test_mesh_tets.h5m"
uscd_mesh = openmc.UnstructuredMesh(mesh_filename)
uscd_mesh.mesh_lib = 'moab'
uscd_filter = openmc.MeshFilter(mesh=uscd_mesh)
# create tallies
tallies = openmc.Tallies()
regular_mesh_tally = openmc.Tally(name="regular mesh tally")
regular_mesh_tally.filters = [regular_mesh_filter]
regular_mesh_tally.scores = ['flux']
regular_mesh_tally.estimator = test_opts['estimator']
tallies.append(regular_mesh_tally)
uscd_tally = openmc.Tally(name="unstructured mesh tally")
uscd_tally.filters = [uscd_filter]
uscd_tally.scores = ['flux']
uscd_tally.estimator = test_opts['estimator']
tallies.append(uscd_tally)
### Settings ###
settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.particles = 100
settings.batches = 10
# source setup
r = openmc.stats.Uniform(a=0.0, b=0.0)
theta = openmc.stats.Discrete(x=[0.0], p=[1.0])
phi = openmc.stats.Discrete(x=[0.0], p=[1.0])
space = openmc.stats.SphericalIndependent(r, theta, phi)
angle = openmc.stats.Monodirectional((-1.0, 0.0, 0.0))
energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0])
source = openmc.Source(space=space, energy=energy, angle=angle)
settings.source = source
model = openmc.model.Model(geometry=geometry,
materials=materials,
tallies=tallies,
settings=settings)
harness = UnstructuredMeshTest('statepoint.10.h5',
model,
test_opts['inputs_true'],
test_opts['holes'])
harness.main()

View file

@ -276,7 +276,7 @@ class ParticleRestartTestHarness(TestHarness):
class PyAPITestHarness(TestHarness):
def __init__(self, statepoint_name, model=None):
def __init__(self, statepoint_name, model=None, inputs_true=None):
super().__init__(statepoint_name)
if model is None:
self._model = pwr_core()
@ -284,6 +284,7 @@ class PyAPITestHarness(TestHarness):
self._model = model
self._model.plots = []
self.inputs_true = "inputs_true.dat" if not inputs_true else inputs_true
def main(self):
"""Accept commandline arguments and either run or update tests."""
@ -342,15 +343,15 @@ class PyAPITestHarness(TestHarness):
def _overwrite_inputs(self):
"""Overwrite inputs_true.dat with inputs_test.dat"""
shutil.copyfile('inputs_test.dat', 'inputs_true.dat')
shutil.copyfile('inputs_test.dat', self.inputs_true)
def _compare_inputs(self):
"""Make sure the current inputs agree with the _true standard."""
compare = filecmp.cmp('inputs_test.dat', 'inputs_true.dat')
compare = filecmp.cmp('inputs_test.dat', self.inputs_true)
if not compare:
expected = open('inputs_true.dat', 'r').readlines()
expected = open(self.inputs_true, 'r').readlines()
actual = open('inputs_test.dat', 'r').readlines()
diff = unified_diff(expected, actual, 'inputs_true.dat',
diff = unified_diff(expected, actual, self.inputs_true,
'inputs_test.dat')
print('Input differences:')
print(''.join(colorize(diff)))