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

Merge pull request #1814 from pshriwise/mft-pr

Browse source 
Mesh Filter Translations
This commit is contained in:
Paul Romano 2021-04-29 17:02:31 -05:00 • committed by GitHub
commit 24318e1416
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
14 changed files with 799 additions and 79 deletions
Download patch file Download diff file Expand all files Collapse all files

2
include/openmc/capi.h
View file

@ -74,6 +74,8 @@ extern "C" {
int openmc_material_filter_set_bins(int32_t index, size_t n, const int32_t* bins);
int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh);
int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_mesh_filter_get_translation(int32_t index, double translation[3]);
int openmc_mesh_filter_set_translation(int32_t index, double translation[3]);
int openmc_mesh_get_id(int32_t index, int32_t* id);
int openmc_mesh_set_id(int32_t index, int32_t id);
int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);

64
include/openmc/mesh.h
View file

@ -71,18 +71,28 @@ public:
//! Determine which bins were crossed by a particle
//
//! \param[in] p Particle to check
//! \param[in] r0 Previous position of the particle
//! \param[in] r1 Current position of the particle
//! \param[in] u Particle direction
//! \param[out] bins Bins that were crossed
//! \param[out] lengths Fraction of tracklength in each bin
virtual void bins_crossed(const Particle& p, std::vector<int>& bins,
virtual void bins_crossed(Position r0,
Position r1,
const Direction& u,
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[in] r0 Previous position of the particle
//! \param[in] r1 Current position of the particle
//! \param[in] u Particle direction
//! \param[out] bins Surface bins that were crossed
virtual void
surface_bins_crossed(const Particle& p, std::vector<int>& bins) const = 0;
surface_bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins) const = 0;
//! Get bin at a given position in space
//
@ -137,7 +147,10 @@ public:
int n_surface_bins() const override;
void bins_crossed(const Particle& p, std::vector<int>& bins,
void bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins,
std::vector<double>& lengths) const override;
//! Count number of bank sites in each mesh bin / energy bin
@ -219,9 +232,11 @@ public:
RegularMesh(pugi::xml_node node);
// Overriden methods
void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
const override;
void
surface_bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins) const override;
int get_index_in_direction(double r, int i) const override;
@ -259,9 +274,11 @@ public:
RectilinearMesh(pugi::xml_node node);
// Overriden methods
void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
const override;
void
surface_bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins) const override;
int get_index_in_direction(double r, int i) const override;
@ -324,8 +341,9 @@ public:
std::string bin_label(int bin) const override;
void surface_bins_crossed(const Particle& p,
std::vector<int>& bins) const override;
void surface_bins_crossed(Position r0,
Position r1,
std::vector<int>& bins) const;
void to_hdf5(hid_t group) const override;
@ -348,10 +366,18 @@ public:
MOABMesh(pugi::xml_node);
MOABMesh(const std::string& filename);
void bins_crossed(const Particle& p,
void bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins,
std::vector<double>& lengths) const override;
void
surface_bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins) const override;
int get_bin(Position r) const;
int n_bins() const override;
@ -496,10 +522,18 @@ public:
LibMesh(const std::string& filename);
// Methods
void bins_crossed(const Particle& p,
void bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins,
std::vector<double>& lengths) const override;
void
surface_bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins) const override;
int get_bin(Position r) const override;
int n_bins() const override;

12
include/openmc/tallies/filter_mesh.h
View file

@ -4,6 +4,7 @@
#include <cstdint>
#include "openmc/tallies/filter.h"
#include "openmc/position.h"
namespace openmc {
@ -42,11 +43,22 @@ public:
virtual void set_mesh(int32_t mesh);
virtual void set_translation(const Position& translation);
virtual void set_translation(const double translation[3]);
virtual const Position& translation() const {return translation_;}
virtual bool translated() const {return translated_;}
protected:
//----------------------------------------------------------------------------
// Data members
int32_t mesh_;
bool translated_ {false};
Position translation_ {0.0, 0.0, 0.0};
};
} // namespace openmc

31
openmc/filter.py
View file

@ -716,6 +716,9 @@ class MeshFilter(Filter):
The mesh object that events will be tallied onto
id : int
Unique identifier for the filter
translation : Iterable of float
This array specifies a vector that is used to translate (shift)
the mesh for this filter
bins : list of tuple
A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1),
...]
@ -727,6 +730,7 @@ class MeshFilter(Filter):
def __init__(self, mesh, filter_id=None):
self.mesh = mesh
self.id = filter_id
self._translation = None
def __hash__(self):
string = type(self).__name__ + '\n'
@ -737,6 +741,7 @@ class MeshFilter(Filter):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id)
string += '{: <16}=\t{}\n'.format('\tID', self.id)
string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation)
return string
@classmethod
@ -754,8 +759,13 @@ class MeshFilter(Filter):
mesh_obj = kwargs['meshes'][mesh_id]
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
out = cls(mesh_obj, filter_id=filter_id)
translation = group.get('translation')
if translation:
out.translation = translation[()]
return out
@property
@ -774,6 +784,16 @@ class MeshFilter(Filter):
else:
self.bins = list(mesh.indices)
@property
def translation(self):
return self._translation
@translation.setter
def translation(self, t):
cv.check_type('mesh filter translation', t, Iterable, Real)
cv.check_length('mesh filter translation', t, 3)
self._translation = np.asarray(t)
def can_merge(self, other):
# Mesh filters cannot have more than one bin
return False
@ -854,6 +874,8 @@ class MeshFilter(Filter):
"""
element = super().to_xml_element()
element[0].text = str(self.mesh.id)
if self.translation is not None:
element.set('translation', ' '.join(map(str, self.translation)))
return element
@ -873,6 +895,9 @@ class MeshSurfaceFilter(MeshFilter):
The mesh ID
mesh : openmc.MeshBase
The mesh object that events will be tallied onto
translation : Iterable of float
This array specifies a vector that is used to translate (shift)
the mesh for this filter
id : int
Unique identifier for the filter
bins : list of tuple
@ -971,9 +996,9 @@ class CollisionFilter(Filter):
Parameters
----------
bins : Iterable of int
A list or iterable of the number of collisions, as integer values.
The events whose post-scattering collision number equals one of
the provided values will be counted.
A list or iterable of the number of collisions, as integer values.
The events whose post-scattering collision number equals one of
the provided values will be counted.
filter_id : int
Unique identifier for the filter

26
openmc/lib/filter.py
View file

@ -84,6 +84,12 @@ _dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler
_dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
_dll.openmc_meshsurface_filter_set_mesh.restype = c_int
_dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
_dll.openmc_mesh_filter_get_translation.argtypes = [c_int32, POINTER(c_double*3)]
_dll.openmc_mesh_filter_get_translation.restype = c_int
_dll.openmc_mesh_filter_get_translation.errcheck = _error_handler
_dll.openmc_mesh_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)]
_dll.openmc_mesh_filter_set_translation.restype = c_int
_dll.openmc_mesh_filter_set_translation.errcheck = _error_handler
_dll.openmc_new_filter.argtypes = [c_char_p, POINTER(c_int32)]
_dll.openmc_new_filter.restype = c_int
_dll.openmc_new_filter.errcheck = _error_handler
@ -325,6 +331,16 @@ class MeshFilter(Filter):
def mesh(self, mesh):
_dll.openmc_mesh_filter_set_mesh(self._index, mesh._index)
@property
def translation(self):
translation = (c_double*3)()
_dll.openmc_mesh_filter_get_translation(self._index, translation)
return tuple(translation)
@translation.setter
def translation(self, translation):
_dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation))
class MeshSurfaceFilter(Filter):
filter_type = 'meshsurface'
@ -344,6 +360,16 @@ class MeshSurfaceFilter(Filter):
def mesh(self, mesh):
_dll.openmc_meshsurface_filter_set_mesh(self._index, mesh._index)
@property
def translation(self):
translation = (c_double*3)()
_dll.openmc_mesh_filter_get_translation(self._index, translation)
return tuple(translation)
@translation.setter
def translation(self, translation):
_dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation))
class MuFilter(Filter):
filter_type = 'mu'

8
openmc/mesh.py
View file

@ -471,6 +471,8 @@ class RectilinearMesh(MeshBase):
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 (always 3 for a RectilinearMesh).
x_grid : Iterable of float
@ -492,6 +494,12 @@ class RectilinearMesh(MeshBase):
self._y_grid = None
self._z_grid = None
@property
def dimension(self):
return (len(self.x_grid) - 1,
len(self.y_grid) - 1,
len(self.z_grid) - 1)
@property
def n_dimension(self):
return 3

123
src/mesh.cpp
View file

@ -190,7 +190,8 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) {
}
void
UnstructuredMesh::surface_bins_crossed(const Particle& p,
UnstructuredMesh::surface_bins_crossed(Position r0,
Position r1,
std::vector<int>& bins) const {
fatal_error("Unstructured mesh surface tallies are not implemented.");
}
@ -581,26 +582,24 @@ bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
return min_dist < INFTY;
}
void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
std::vector<double>& lengths) const
void StructuredMesh::bins_crossed(Position r0,
Position r1,
const Direction& u,
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();
double total_distance = (r1 - r0).norm();
if (total_distance == 0.0) return;
// 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;
Position last_r = r0 + TINY_BIT*u;
Position r = r1 - TINY_BIT*u;
// Determine the mesh indices for the starting and ending coords. Here, we
// use arrays for ijk0 and ijk1 instead of std::vector because we obtain a
@ -611,21 +610,21 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
int n = n_dimension_;
int ijk0[3], ijk1[3];
bool start_in_mesh;
get_indices(r0, ijk0, &start_in_mesh);
get_indices(last_r, ijk0, &start_in_mesh);
bool end_in_mesh;
get_indices(r1, ijk1, &end_in_mesh);
get_indices(r, 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;
last_r = r0;
} 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;
if (!intersects(last_r, r, ijk0)) return;
}
r1 = r;
r = r1;
// The TINY_BIT offsets above mean that the preceding logic cannot always find
// the correct ijk0 and ijk1 indices. For tracks shorter than 2*TINY_BIT, just
@ -644,7 +643,7 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
if (std::equal(ijk0, ijk0 + n, 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();
double distance = (last_r - r).norm();
bins.push_back(get_bin_from_indices(ijk0));
lengths.push_back(distance / total_distance);
break;
@ -657,10 +656,10 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
d[k] = INFTY;
} else if (u[k] > 0) {
double xyz_cross = positive_grid_boundary(ijk0, k);
d[k] = (xyz_cross - r0[k]) / u[k];
d[k] = (xyz_cross - last_r[k]) / u[k];
} else {
double xyz_cross = negative_grid_boundary(ijk0, k);
d[k] = (xyz_cross - r0[k]) / u[k];
d[k] = (xyz_cross - last_r[k]) / u[k];
}
}
@ -671,7 +670,7 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
lengths.push_back(distance / total_distance);
// Translate to the oncoming mesh surface.
r0 += distance * u;
last_r += distance * u;
// Increment the indices into the next mesh cell.
if (u[j] > 0.0) {
@ -797,17 +796,12 @@ double RegularMesh::negative_grid_boundary(int* ijk, int i) const
return lower_left_[i] + (ijk[i] - 1) * width_[i];
}
void RegularMesh::surface_bins_crossed(const Particle& p,
std::vector<int>& bins) const
void
RegularMesh::surface_bins_crossed(Position r0,
Position r1,
const Direction& u,
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 n = n_dimension_;
std::vector<int> ijk0(n), ijk1(n);
@ -1107,17 +1101,11 @@ int RectilinearMesh::set_grid()
return 0;
}
void RectilinearMesh::surface_bins_crossed(const Particle& p,
void RectilinearMesh::surface_bins_crossed(Position r0,
Position r1,
const Direction& u,
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;
@ -1444,7 +1432,7 @@ openmc_get_mesh_index(int32_t id, int32_t* index)
return 0;
}
// Return the ID of a mesh
//! Return the ID of a mesh
extern "C" int
openmc_mesh_get_id(int32_t index, int32_t* id)
{
@ -1701,26 +1689,25 @@ MOABMesh::intersect_track(const moab::CartVect& start,
}
void
MOABMesh::bins_crossed(const Particle& p,
MOABMesh::bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins,
std::vector<double>& lengths) const
{
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 start(r0.x, r0.y, r0.z);
moab::CartVect end(r1.x, r1.y, r1.z);
moab::CartVect dir(u.x, u.y, u.z);
dir.normalize();
double track_len = (r1 - r0).length();
double track_len = (end - start).length();
if (track_len == 0.0) return;
r0 -= TINY_BIT * dir;
r1 += TINY_BIT * dir;
start -= TINY_BIT * dir;
end += TINY_BIT * dir;
std::vector<double> hits;
intersect_track(r0, dir, track_len, hits);
intersect_track(start, dir, track_len, hits);
bins.clear();
lengths.clear();
@ -1729,7 +1716,7 @@ MOABMesh::bins_crossed(const Particle& p,
// 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);
Position midpoint = r0 + u * (track_len * 0.5);
int bin = this->get_bin(midpoint);
if (bin != -1) {
bins.push_back(bin);
@ -1740,7 +1727,7 @@ MOABMesh::bins_crossed(const Particle& p,
// for each segment in the set of tracks, try to look up a tet
// at the midpoint of the segment
Position current = last_r;
Position current = r0;
double last_dist = 0.0;
for (const auto& hit : hits) {
// get the segment length
@ -1752,7 +1739,7 @@ MOABMesh::bins_crossed(const Particle& p,
int bin = this->get_bin(midpoint);
// determine the start point for this segment
current = last_r + u * hit;
current = r0 + u * hit;
if (bin == -1) {
continue;
@ -1767,7 +1754,7 @@ MOABMesh::bins_crossed(const Particle& p,
// 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();
Position segment_start = r0 + 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);
@ -1833,6 +1820,16 @@ double MOABMesh::tet_volume(moab::EntityHandle tet) const
return 1.0 / 6.0 * (((p[1] - p[0]) * (p[2] - p[0])) % (p[3] - p[0]));
}
void MOABMesh::surface_bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins) const
{
// TODO: Implement triangle crossings here
throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."};
}
int MOABMesh::get_bin(Position r) const
{
moab::EntityHandle tet = get_tet(r);
@ -2212,8 +2209,18 @@ int LibMesh::n_bins() const
return m_->n_elem();
}
int LibMesh::n_surface_bins() const
void
LibMesh::surface_bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins) const
{
// TODO: Implement triangle crossings here
throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."};
}
int
LibMesh::n_surface_bins() const {
int n_bins = 0;
for (int i = 0; i < this->n_bins(); i++) {
const libMesh::Elem& e = get_element_from_bin(i);
@ -2300,7 +2307,9 @@ void LibMesh::write(const std::string& filename) const
}
void
LibMesh::bins_crossed(const Particle& p,
LibMesh::bins_crossed(Position r0,
Position r1,
const Direction& u,
std::vector<int>& bins,
std::vector<double>& lengths) const
{

80
src/tallies/filter_mesh.cpp
View file

@ -28,20 +28,36 @@ MeshFilter::from_xml(pugi::xml_node node)
fatal_error(fmt::format(
"Could not find mesh {} specified on tally filter.", id));
}
if (check_for_node(node, "translation")) {
set_translation(get_node_array<double>(node, "translation"));
}
}
void
MeshFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match)
const
{
Position last_r = p.r_last_;
Position r = p.r();
Position u = p.u();
// apply translation if present
if (translated_) {
last_r -= translation();
r -= translation();
}
if (estimator != TallyEstimator::TRACKLENGTH) {
auto bin = model::meshes[mesh_]->get_bin(p.r());
auto bin = model::meshes[mesh_]->get_bin(r);
if (bin >= 0) {
match.bins_.push_back(bin);
match.weights_.push_back(1.0);
}
} else {
model::meshes[mesh_]->bins_crossed(p, match.bins_, match.weights_);
model::meshes[mesh_]->bins_crossed(last_r, r, u, match.bins_, match.weights_);
}
}
@ -50,13 +66,17 @@ MeshFilter::to_statepoint(hid_t filter_group) const
{
Filter::to_statepoint(filter_group);
write_dataset(filter_group, "bins", model::meshes[mesh_]->id_);
if (translated_) {
write_dataset(filter_group, "translation", translation_);
}
}
std::string
MeshFilter::text_label(int bin) const
{
auto& mesh = *model::meshes.at(mesh_);
return mesh.bin_label(bin);
std::string label = mesh.bin_label(bin);
return label;
}
void
@ -66,6 +86,17 @@ MeshFilter::set_mesh(int32_t mesh)
n_bins_ = model::meshes[mesh_]->n_bins();
}
void MeshFilter::set_translation(const Position& translation)
{
translated_ = true;
translation_ = translation;
}
void MeshFilter::set_translation(const double translation[3])
{
this->set_translation({translation[0], translation[1], translation[2]});
}
//==============================================================================
// C-API functions
//==============================================================================
@ -123,4 +154,47 @@ openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh)
return 0;
}
extern "C" int
openmc_mesh_filter_get_translation(int32_t index, double translation[3])
{
// Make sure this is a valid index to an allocated filter
if (int err = verify_filter(index)) return err;
// Check the filter type
const auto& filter = model::tally_filters[index];
if (filter->type() != "mesh" && filter->type() != "meshsurface") {
set_errmsg("Tried to get a translation from a non-mesh-based filter.");
return OPENMC_E_INVALID_TYPE;
}
// Get translation from the mesh filter and set value
auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
const auto& t = mesh_filter->translation();
for (int i = 0; i < 3; i++) { translation[i] = t[i]; }
return 0;
}
extern "C" int
openmc_mesh_filter_set_translation(int32_t index, double translation[3])
{
// Make sure this is a valid index to an allocated filter
if (int err = verify_filter(index)) return err;
const auto& filter = model::tally_filters[index];
// Check the filter type
if (filter->type() != "mesh" && filter->type() != "meshsurface") {
set_errmsg("Tried to set mesh on a non-mesh-based filter.");
return OPENMC_E_INVALID_TYPE;
}
// Get a pointer to the filter and downcast
auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
// Set the translation
mesh_filter->set_translation(translation);
return 0;
}
} // namespace openmc

10
src/tallies/filter_meshsurface.cpp
View file

@ -11,7 +11,15 @@ void
MeshSurfaceFilter::get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const
{
model::meshes[mesh_]->surface_bins_crossed(p, match.bins_);
Position r0 = p.r_last_current_;
Position r1 = p.r();
if (translated_) {
r0 -= translation();
r1 -= translation();
}
Direction u = p.u();
model::meshes[mesh_]->surface_bins_crossed(r0, r1, u, match.bins_);
for (auto b : match.bins_) match.weights_.push_back(1.0);
}

0
tests/regression_tests/filter_translations/__init__.py Normal file
View file

84
tests/regression_tests/filter_translations/inputs_true.dat Normal file
View file

@ -0,0 +1,84 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2 3 -4 10 -9" universe="1" />
<cell id="2" material="2" region="~(1 -2 3 -4) (5 -6 7 -8) 10 -9" 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 boundary="reflective" coeffs="-10.0" id="5" name="minimum x" type="x-plane" />
<surface boundary="reflective" coeffs="10.0" id="6" name="maximum x" type="x-plane" />
<surface boundary="reflective" coeffs="-10.0" id="7" name="minimum y" type="y-plane" />
<surface boundary="reflective" coeffs="10.0" id="8" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="10.0" id="9" type="z-plane" />
<surface boundary="vacuum" coeffs="-10.0" id="10" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material depletable="true" id="1">
<density units="g/cm3" value="10.0" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="2">
<density units="g/cm3" value="1.0" />
<nuclide ao="1.0" name="Zr90" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>5</batches>
<inactive>0</inactive>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1">
<dimension>3 4 5</dimension>
<lower_left>-9 -9 -9</lower_left>
<upper_right>9 9 9</upper_right>
</mesh>
<mesh id="2" type="rectilinear">
<x_grid>-9.0 0.0 9.0</x_grid>
<y_grid>-9.0 -3.0 3.0 9.0</y_grid>
<z_grid>-9.0 -4.5 0.0 4.5 9.0</z_grid>
</mesh>
<mesh id="3">
<dimension>3 4 5</dimension>
<lower_left>-19 -4 -9</lower_left>
<upper_right>-1 14 9</upper_right>
</mesh>
<mesh id="4" type="rectilinear">
<x_grid>-19.0 -10.0 -1.0</x_grid>
<y_grid>-4.0 2.0 8.0 14.0</y_grid>
<z_grid>-9.0 -4.5 0.0 4.5 9.0</z_grid>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="2" type="mesh">
<bins>2</bins>
</filter>
<filter id="3" translation="10 -5 0" type="mesh">
<bins>3</bins>
</filter>
<filter id="4" translation="10 -5 0" type="mesh">
<bins>4</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>total</scores>
</tally>
<tally id="2">
<filters>2</filters>
<scores>total</scores>
</tally>
<tally id="3">
<filters>3</filters>
<scores>total</scores>
</tally>
<tally id="4">
<filters>4</filters>
<scores>total</scores>
</tally>
</tallies>

342
tests/regression_tests/filter_translations/results_true.dat Normal file
View file

@ -0,0 +1,342 @@
k-combined:
7.656044E-01 5.168921E-02
tally 1:
4.195134E-02
3.639654E-04
7.345058E-02
1.082591E-03
5.325577E-02
5.772001E-04
1.390109E-01
3.954675E-03
3.829062E-01
3.001419E-02
1.408217E-01
4.113283E-03
1.412139E-01
4.218445E-03
3.567122E-01
2.663446E-02
1.393013E-01
3.978946E-03
5.117767E-02
5.394412E-04
8.111129E-02
1.362452E-03
5.594852E-02
6.333549E-04
6.298061E-02
8.073054E-04
9.082908E-02
1.668351E-03
6.783427E-02
9.250122E-04
1.999222E-01
8.028088E-03
5.909662E-01
7.061368E-02
1.896895E-01
7.242074E-03
1.877543E-01
7.145844E-03
6.189100E-01
7.667032E-02
1.847604E-01
6.905108E-03
6.017074E-02
7.311068E-04
1.109867E-01
2.474797E-03
7.322493E-02
1.112742E-03
6.859899E-02
9.576740E-04
1.062208E-01
2.360724E-03
6.762553E-02
9.247454E-04
2.367276E-01
1.130443E-02
9.674870E-01
2.606562E-01
2.113618E-01
8.988266E-03
2.079097E-01
8.765092E-03
9.535566E-01
2.216359E-01
2.131094E-01
9.143533E-03
6.644946E-02
9.024189E-04
1.064224E-01
2.284439E-03
8.143154E-02
1.362849E-03
6.209959E-02
7.875797E-04
9.702493E-02
1.902685E-03
6.416842E-02
8.369934E-04
1.727520E-01
6.027045E-03
5.966710E-01
7.197863E-02
1.889340E-01
7.213441E-03
1.914452E-01
7.420265E-03
5.726219E-01
6.723631E-02
1.646287E-01
5.489910E-03
5.663491E-02
6.494568E-04
9.923269E-02
2.053123E-03
6.809958E-02
9.386140E-04
4.850339E-02
4.856214E-04
8.169408E-02
1.354206E-03
5.474566E-02
6.128204E-04
1.208954E-01
3.074917E-03
3.393101E-01
2.398955E-02
1.358915E-01
3.758521E-03
1.347810E-01
3.718991E-03
3.404279E-01
2.380238E-02
1.315829E-01
3.585161E-03
4.754269E-02
4.593490E-04
7.591439E-02
1.170776E-03
5.121227E-02
5.346435E-04
tally 2:
2.194240E-01
9.811952E-03
2.718491E-01
1.501408E-02
5.823038E-01
6.967280E-02
5.955012E-01
7.246929E-02
2.496521E-01
1.268521E-02
2.623766E-01
1.397567E-02
3.475667E-01
2.453188E-02
3.290723E-01
2.187833E-02
1.123843E+00
2.718539E-01
1.060797E+00
2.384769E-01
3.373095E-01
2.280788E-02
3.563621E-01
2.594087E-02
3.442434E-01
2.400708E-02
3.319540E-01
2.206590E-02
1.051427E+00
2.459968E-01
1.065112E+00
2.465202E-01
3.182627E-01
2.049613E-02
3.724575E-01
2.824300E-02
2.274581E-01
1.052009E-02
2.558629E-01
1.318940E-02
5.434276E-01
6.081708E-02
5.623197E-01
6.408372E-02
2.320052E-01
1.096050E-02
2.380175E-01
1.168927E-02
tally 3:
4.195134E-02
3.639654E-04
7.345058E-02
1.082591E-03
5.325577E-02
5.772001E-04
1.390109E-01
3.954675E-03
3.829062E-01
3.001419E-02
1.408217E-01
4.113283E-03
1.412139E-01
4.218445E-03
3.567122E-01
2.663446E-02
1.393013E-01
3.978946E-03
5.117767E-02
5.394412E-04
8.111129E-02
1.362452E-03
5.594852E-02
6.333549E-04
6.298061E-02
8.073054E-04
9.082908E-02
1.668351E-03
6.783427E-02
9.250122E-04
1.999222E-01
8.028088E-03
5.909662E-01
7.061368E-02
1.896895E-01
7.242074E-03
1.877543E-01
7.145844E-03
6.189100E-01
7.667032E-02
1.847604E-01
6.905108E-03
6.017074E-02
7.311068E-04
1.109867E-01
2.474797E-03
7.322493E-02
1.112742E-03
6.859899E-02
9.576740E-04
1.062208E-01
2.360724E-03
6.762553E-02
9.247454E-04
2.367276E-01
1.130443E-02
9.674870E-01
2.606562E-01
2.113618E-01
8.988266E-03
2.079097E-01
8.765092E-03
9.535566E-01
2.216359E-01
2.131094E-01
9.143533E-03
6.644946E-02
9.024189E-04
1.064224E-01
2.284439E-03
8.143154E-02
1.362849E-03
6.209959E-02
7.875797E-04
9.702493E-02
1.902685E-03
6.416842E-02
8.369934E-04
1.727520E-01
6.027045E-03
5.966710E-01
7.197863E-02
1.889340E-01
7.213441E-03
1.914452E-01
7.420265E-03
5.726219E-01
6.723631E-02
1.646287E-01
5.489910E-03
5.663491E-02
6.494568E-04
9.923269E-02
2.053123E-03
6.809958E-02
9.386140E-04
4.850339E-02
4.856214E-04
8.169408E-02
1.354206E-03
5.474566E-02
6.128204E-04
1.208954E-01
3.074917E-03
3.393101E-01
2.398955E-02
1.358915E-01
3.758521E-03
1.347810E-01
3.718991E-03
3.404279E-01
2.380238E-02
1.315829E-01
3.585161E-03
4.754269E-02
4.593490E-04
7.591439E-02
1.170776E-03
5.121227E-02
5.346435E-04
tally 4:
2.194240E-01
9.811952E-03
2.718491E-01
1.501408E-02
5.823038E-01
6.967280E-02
5.955012E-01
7.246929E-02
2.496521E-01
1.268521E-02
2.623766E-01
1.397567E-02
3.475667E-01
2.453188E-02
3.290723E-01
2.187833E-02
1.123843E+00
2.718539E-01
1.060797E+00
2.384769E-01
3.373095E-01
2.280788E-02
3.563621E-01
2.594087E-02
3.442434E-01
2.400708E-02
3.319540E-01
2.206590E-02
1.051427E+00
2.459968E-01
1.065112E+00
2.465202E-01
3.182627E-01
2.049613E-02
3.724575E-01
2.824300E-02
2.274581E-01
1.052009E-02
2.558629E-01
1.318940E-02
5.434276E-01
6.081708E-02
5.623197E-01
6.408372E-02
2.320052E-01
1.096050E-02
2.380175E-01
1.168927E-02

90
tests/regression_tests/filter_translations/test.py Normal file
View file

@ -0,0 +1,90 @@
import numpy as np
import openmc
import pytest
from tests.testing_harness import PyAPITestHarness
@pytest.fixture
def model():
model = openmc.model.Model()
fuel = openmc.Material()
fuel.set_density('g/cm3', 10.0)
fuel.add_nuclide('U235', 1.0)
zr = openmc.Material()
zr.set_density('g/cm3', 1.0)
zr.add_nuclide('Zr90', 1.0)
model.materials.extend([fuel, zr])
box1 = openmc.model.rectangular_prism(10.0, 10.0)
box2 = openmc.model.rectangular_prism(20.0, 20.0, boundary_type='reflective')
top = openmc.ZPlane(z0=10.0, boundary_type='vacuum')
bottom = openmc.ZPlane(z0=-10.0, boundary_type='vacuum')
cell1 = openmc.Cell(fill=fuel, region=box1 & +bottom & -top)
cell2 = openmc.Cell(fill=zr, region=~box1 & box2 & +bottom & -top)
model.geometry = openmc.Geometry([cell1, cell2])
model.settings.batches = 5
model.settings.inactive = 0
model.settings.particles = 1000
translation = np.array((10, -5, 0))
llc = np.array([-9, -9, -9])
urc = np.array([9, 9, 9])
mesh_dims = (3, 4, 5)
filters = []
# un-translated meshes
reg_mesh = openmc.RegularMesh()
reg_mesh.dimension = mesh_dims
reg_mesh.lower_left = llc
reg_mesh.upper_right = urc
filters.append(openmc.MeshFilter(reg_mesh))
recti_mesh = openmc.RectilinearMesh()
recti_mesh.x_grid = np.linspace(llc[0], urc[0], mesh_dims[0])
recti_mesh.y_grid = np.linspace(llc[1], urc[1], mesh_dims[1])
recti_mesh.z_grid = np.linspace(llc[2], urc[2], mesh_dims[2])
filters.append(openmc.MeshFilter(recti_mesh))
llc = np.array(llc - translation)
urc = np.array(urc - translation)
# translated meshes
translated_reg_mesh = openmc.RegularMesh()
translated_reg_mesh.dimension = mesh_dims
translated_reg_mesh.lower_left = llc
translated_reg_mesh.upper_right = urc
filters.append(openmc.MeshFilter(translated_reg_mesh))
filters[-1].translation = translation
translated_recti_mesh = openmc.RectilinearMesh()
translated_recti_mesh.x_grid = np.linspace(llc[0], urc[0], mesh_dims[0])
translated_recti_mesh.y_grid = np.linspace(llc[1], urc[1], mesh_dims[1])
translated_recti_mesh.z_grid = np.linspace(llc[2], urc[2], mesh_dims[2])
filters.append(openmc.MeshFilter(translated_recti_mesh))
filters[-1].translation = translation
# Create tallies
for f in filters:
tally = openmc.Tally()
tally.filters = [f]
tally.scores = ['total']
model.tallies.append(tally)
return model
def test_filter_mesh_translations(model):
harness = PyAPITestHarness('statepoint.5.h5', model)
harness.main()

6
tests/unit_tests/test_lib.py
View file

@ -470,11 +470,17 @@ def test_regular_mesh(lib_init):
assert isinstance(mesh, openmc.lib.RegularMesh)
assert mesh_id == mesh.id
translation = (1.0, 2.0, 3.0)
mf = openmc.lib.MeshFilter(mesh)
assert mf.mesh == mesh
mf.translation = translation
assert mf.translation == translation
msf = openmc.lib.MeshSurfaceFilter(mesh)
assert msf.mesh == mesh
msf.translation = translation
assert msf.translation == translation
def test_rectilinear_mesh(lib_init):