crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1701 from smharper/bc_update

Browse source 
Refactor and generalize the boundary condition implementation
This commit is contained in:
Paul Romano 2020-12-22 15:41:11 -06:00 • committed by GitHub
commit 10715497c1
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
16 changed files with 779 additions and 442 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -225,6 +225,7 @@ target_compile_options(faddeeva PRIVATE ${cxxflags})
list(APPEND libopenmc_SOURCES
src/bank.cpp
src/boundary_condition.cpp
src/bremsstrahlung.cpp
src/dagmc.cpp
src/cell.cpp

24
docs/source/usersguide/geometry.rst
View file

@ -162,22 +162,26 @@ surface. To specify a vacuum boundary condition, simply change the
Reflective and periodic boundary conditions can be set with the strings
'reflective' and 'periodic'. Vacuum and reflective boundary conditions can be
applied to any type of surface. Periodic boundary conditions can be applied to
pairs of planar surfaces. For axis-aligned planes, matching periodic surfaces
can be determined automatically. For non-axis-aligned planes, it is necessary to
specify pairs explicitly using the :attr:`Surface.periodic_surface` attribute as
in the following example::
pairs of planar surfaces. If there are only two periodic surfaces they will be
matched automatically. Otherwise it is necessary to specify pairs explicitly
using the :attr:`Surface.periodic_surface` attribute as in the following
example::
p1 = openmc.Plane(a=0.3, b=5.0, d=1.0, boundary_type='periodic')
p2 = openmc.Plane(a=0.3, b=5.0, d=-1.0, boundary_type='periodic')
p1.periodic_surface = p2
Rotationally-periodic boundary conditions can be specified for a pair of
:class:`XPlane` and :class:`YPlane`; in that case, the
:attr:`Surface.periodic_surface` attribute must be specified manually as well.
Both rotational and translational periodic boundary conditions are specified in
the same fashion. If both planes have the same normal vector, a translational
periodicity is assumed; rotational periodicity is assumed otherwise. Currently,
only rotations about the :math:`z`-axis are supported.
.. caution:: When using rotationally-periodic boundary conditions, your geometry
must be defined in the first quadrant, i.e., above the y-plane and
to the right of the x-plane.
For a rotational periodic BC, the normal vectors of each surface must point
inwards---towards the valid geometry. For example, a :class:`XPlane` and
:class:`YPlane` would be valid for a 90-degree periodic rotation if the geometry
lies in the first quadrant of the Cartesian grid. If the geometry instead lies
in the fourth quadrant, the :class:`YPlane` must be replaced by a
:class:`Plane` with the normal vector pointing in the :math:`-y` direction.
.. _usersguide_cells:

118
include/openmc/boundary_condition.h Normal file
View file

@ -0,0 +1,118 @@
#ifndef OPENMC_BOUNDARY_CONDITION_H
#define OPENMC_BOUNDARY_CONDITION_H
#include "openmc/position.h"
namespace openmc {
// Forward declare some types used in function arguments.
class Particle;
class Surface;
//==============================================================================
//! A class that tells particles what to do after they strike an outer boundary.
//==============================================================================
class BoundaryCondition {
public:
//! Perform tracking operations for a particle that strikes the boundary.
//! \param p The particle that struck the boundary. This class is not meant
//! to directly modify anything about the particle, but it will do so
//! indirectly by calling the particle's appropriate cross_*_bc function.
//! \param surf The specific surface on the boundary the particle struck.
virtual void
handle_particle(Particle& p, const Surface& surf) const = 0;
//! Return a string classification of this BC.
virtual std::string type() const = 0;
};
//==============================================================================
//! A BC that kills particles, indicating they left the problem.
//==============================================================================
class VacuumBC : public BoundaryCondition {
public:
void
handle_particle(Particle& p, const Surface& surf) const override;
std::string type() const override {return "vacuum";}
};
//==============================================================================
//! A BC that returns particles via specular reflection.
//==============================================================================
class ReflectiveBC : public BoundaryCondition {
public:
void
handle_particle(Particle& p, const Surface& surf) const override;
std::string type() const override {return "reflective";}
};
//==============================================================================
//! A BC that returns particles via diffuse reflection.
//==============================================================================
class WhiteBC : public BoundaryCondition {
public:
void
handle_particle(Particle& p, const Surface& surf) const override;
std::string type() const override {return "white";}
};
//==============================================================================
//! A BC that moves particles to another part of the problem.
//==============================================================================
class PeriodicBC : public BoundaryCondition {
public:
PeriodicBC(int i_surf, int j_surf)
: i_surf_(i_surf), j_surf_(j_surf)
{};
std::string type() const override {return "periodic";}
protected:
int i_surf_;
int j_surf_;
};
//==============================================================================
//! A BC that moves particles to another part of the problem without rotation.
//==============================================================================
class TranslationalPeriodicBC : public PeriodicBC {
public:
TranslationalPeriodicBC(int i_surf, int j_surf);
void
handle_particle(Particle& p, const Surface& surf) const override;
protected:
//! Vector along which incident particles will be moved
Position translation_;
};
//==============================================================================
//! A BC that rotates particles about a global axis.
//
//! Currently only rotations about the z-axis are supported.
//==============================================================================
class RotationalPeriodicBC : public PeriodicBC {
public:
RotationalPeriodicBC(int i_surf, int j_surf);
void
handle_particle(Particle& p, const Surface& surf) const override;
protected:
//! Angle about the axis by which particle coordinates will be rotated
double angle_;
};
} // namespace openmc
#endif // OPENMC_BOUNDARY_CONDITION_H

27
include/openmc/particle.h
View file

@ -40,6 +40,9 @@ constexpr double CACHE_INVALID {-1.0};
// Class declarations
//==============================================================================
// Forward declare the Surface class for use in function arguments.
class Surface;
class LocalCoord {
public:
void rotate(const std::vector<double>& rotation);
@ -237,6 +240,30 @@ public:
//! Cross a surface and handle boundary conditions
void cross_surface();
//! Cross a vacuum boundary condition.
//
//! \param surf The surface (with the vacuum boundary condition) that the
//! particle struck.
void cross_vacuum_bc(const Surface& surf);
//! Cross a reflective boundary condition.
//
//! \param surf The surface (with the reflective boundary condition) that the
//! particle struck.
//! \param new_u The direction of the particle after reflection.
void cross_reflective_bc(const Surface& surf, Direction new_u);
//! Cross a periodic boundary condition.
//
//! \param surf The surface (with the periodic boundary condition) that the
//! particle struck.
//! \param new_r The position of the particle after translation/rotation.
//! \param new_u The direction of the particle after translation/rotation.
//! \param new_surface The signed index of the surface that the particle will
//! reside on after translation/rotation.
void cross_periodic_bc(const Surface& surf, Position new_r, Direction new_u,
int new_surface);
//! mark a particle as lost and create a particle restart file
//! \param message A warning message to display
void mark_as_lost(const char* message);

60
include/openmc/surface.h
View file

@ -10,6 +10,7 @@
#include "hdf5.h"
#include "pugixml.hpp"
#include "openmc/boundary_condition.h"
#include "openmc/constants.h"
#include "openmc/particle.h"
#include "openmc/position.h"
@ -84,18 +85,9 @@ class Surface
{
public:
// Types of available boundary conditions on a surface
enum class BoundaryType {
TRANSMIT,
VACUUM,
REFLECT,
PERIODIC,
WHITE
};
int id_; //!< Unique ID
BoundaryType bc_; //!< Boundary condition
std::string name_; //!< User-defined name
int id_; //!< Unique ID
std::string name_; //!< User-defined name
std::shared_ptr<BoundaryCondition> bc_ {nullptr}; //!< Boundary condition
explicit Surface(pugi::xml_node surf_node);
Surface();
@ -141,7 +133,6 @@ public:
//! Write all information needed to reconstruct the surface to an HDF5 group.
//! \param group_id An HDF5 group id.
//TODO: this probably needs to include i_periodic for PeriodicSurface
virtual void to_hdf5(hid_t group_id) const = 0;
//! Get the BoundingBox for this surface.
@ -180,33 +171,6 @@ public:
int32_t dag_index_; //!< DagMC index of surface
};
#endif
//==============================================================================
//! A `Surface` that supports periodic boundary conditions.
//!
//! Translational periodicity is supported for the `XPlane`, `YPlane`, `ZPlane`,
//! and `Plane` types. Rotational periodicity is supported for
//! `XPlane`-`YPlane` pairs.
//==============================================================================
class PeriodicSurface : public CSGSurface
{
public:
int i_periodic_{C_NONE}; //!< Index of corresponding periodic surface
explicit PeriodicSurface(pugi::xml_node surf_node);
//! Translate a particle onto this surface from a periodic partner surface.
//! \param other A pointer to the partner surface in this periodic BC.
//! \param r A point on the partner surface that will be translated onto
//! this surface.
//! \param u A direction that will be rotated for systems with rotational
//! periodicity.
//! \return true if this surface and its partner make a rotationally-periodic
//! boundary condition.
virtual bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const = 0;
};
//==============================================================================
//! A plane perpendicular to the x-axis.
@ -214,7 +178,7 @@ public:
//! The plane is described by the equation \f$x - x_0 = 0\f$
//==============================================================================
class SurfaceXPlane : public PeriodicSurface
class SurfaceXPlane : public CSGSurface
{
public:
explicit SurfaceXPlane(pugi::xml_node surf_node);
@ -222,8 +186,6 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box(bool pos_side) const;
double x0_;
@ -235,7 +197,7 @@ public:
//! The plane is described by the equation \f$y - y_0 = 0\f$
//==============================================================================
class SurfaceYPlane : public PeriodicSurface
class SurfaceYPlane : public CSGSurface
{
public:
explicit SurfaceYPlane(pugi::xml_node surf_node);
@ -243,8 +205,6 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box(bool pos_side) const;
double y0_;
@ -256,7 +216,7 @@ public:
//! The plane is described by the equation \f$z - z_0 = 0\f$
//==============================================================================
class SurfaceZPlane : public PeriodicSurface
class SurfaceZPlane : public CSGSurface
{
public:
explicit SurfaceZPlane(pugi::xml_node surf_node);
@ -264,8 +224,6 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box(bool pos_side) const;
double z0_;
@ -277,7 +235,7 @@ public:
//! The plane is described by the equation \f$A x + B y + C z - D = 0\f$
//==============================================================================
class SurfacePlane : public PeriodicSurface
class SurfacePlane : public CSGSurface
{
public:
explicit SurfacePlane(pugi::xml_node surf_node);
@ -285,8 +243,6 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
double A_, B_, C_, D_;
};

261
src/boundary_condition.cpp Normal file
View file

@ -0,0 +1,261 @@
#include "openmc/boundary_condition.h"
#include <exception>
#include <fmt/core.h>
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/surface.h"
namespace openmc {
//==============================================================================
// VacuumBC implementation
//==============================================================================
void
VacuumBC::handle_particle(Particle& p, const Surface& surf) const
{
p.cross_vacuum_bc(surf);
}
//==============================================================================
// ReflectiveBC implementation
//==============================================================================
void
ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
{
Direction u = surf.reflect(p.r(), p.u(), &p);
u /= u.norm();
p.cross_reflective_bc(surf, u);
}
//==============================================================================
// WhiteBC implementation
//==============================================================================
void
WhiteBC::handle_particle(Particle& p, const Surface& surf) const
{
Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed());
u /= u.norm();
p.cross_reflective_bc(surf, u);
}
//==============================================================================
// TranslationalPeriodicBC implementation
//==============================================================================
TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf)
: PeriodicBC(i_surf, j_surf)
{
Surface& surf1 {*model::surfaces[i_surf_]};
Surface& surf2 {*model::surfaces[j_surf_]};
// Make sure the first surface has an appropriate type.
if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
} else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
} else if (const auto* ptr = dynamic_cast<const SurfaceZPlane*>(&surf1)) {
} else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf1)) {
} else {
throw std::invalid_argument(fmt::format("Surface {} is an invalid type for "
"translational periodic BCs. Only planes are supported for these BCs.",
surf1.id_));
}
// Make sure the second surface has an appropriate type.
if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
} else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
} else if (const auto* ptr = dynamic_cast<const SurfaceZPlane*>(&surf2)) {
} else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf2)) {
} else {
throw std::invalid_argument(fmt::format("Surface {} is an invalid type for "
"translational periodic BCs. Only planes are supported for these BCs.",
surf2.id_));
}
// Compute the distance from the first surface to the origin. Check the
// surface evaluate function to decide if the distance is positive, negative,
// or zero.
Position origin {0, 0, 0};
Direction u = surf1.normal(origin);
double d1;
double e1 = surf1.evaluate(origin);
if (e1 > FP_COINCIDENT) {
d1 = -surf1.distance(origin, -u, false);
} else if (e1 < -FP_COINCIDENT) {
d1 = surf1.distance(origin, u, false);
} else {
d1 = 0.0;
}
// Compute the distance from the second surface to the origin.
double d2;
double e2 = surf2.evaluate(origin);
if (e2 > FP_COINCIDENT) {
d2 = -surf2.distance(origin, -u, false);
} else if (e2 < -FP_COINCIDENT) {
d2 = surf2.distance(origin, u, false);
} else {
d2 = 0.0;
}
// Set the translation vector; it's length is the difference in the two
// distances.
translation_ = u * (d2 - d1);
}
void
TranslationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const
{
// TODO: off-by-one on surface indices throughout this function.
int i_particle_surf = std::abs(p.surface_) - 1;
// Figure out which of the two BC surfaces were struck then find the
// particle's new location and surface.
Position new_r;
int new_surface;
if (i_particle_surf == i_surf_) {
new_r = p.r() + translation_;
new_surface = p.surface_ > 0 ? j_surf_ + 1 : -(j_surf_ + 1);
} else if (i_particle_surf == j_surf_) {
new_r = p.r() - translation_;
new_surface = p.surface_ > 0 ? i_surf_ + 1 : -(i_surf_ + 1);
} else {
throw std::runtime_error("Called BoundaryCondition::handle_particle after "
"hitting a surface, but that surface is not recognized by the BC.");
}
// Pass the new location and surface to the particle.
p.cross_periodic_bc(surf, new_r, p.u(), new_surface);
}
//==============================================================================
// RotationalPeriodicBC implementation
//==============================================================================
RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
: PeriodicBC(i_surf, j_surf)
{
Surface& surf1 {*model::surfaces[i_surf_]};
Surface& surf2 {*model::surfaces[j_surf_]};
// Check the type of the first surface
bool surf1_is_xyplane;
if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
surf1_is_xyplane = true;
} else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
surf1_is_xyplane = true;
} else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf1)) {
surf1_is_xyplane = false;
} else {
throw std::invalid_argument(fmt::format("Surface {} is an invalid type for "
"rotational periodic BCs. Only x-planes, y-planes, or general planes "
"(that are perpendicular to z) are supported for these BCs.", surf1.id_));
}
// Check the type of the second surface
bool surf2_is_xyplane;
if (const auto* ptr = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
surf2_is_xyplane = true;
} else if (const auto* ptr = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
surf2_is_xyplane = true;
} else if (const auto* ptr = dynamic_cast<const SurfacePlane*>(&surf2)) {
surf2_is_xyplane = false;
} else {
throw std::invalid_argument(fmt::format("Surface {} is an invalid type for "
"rotational periodic BCs. Only x-planes, y-planes, or general planes "
"(that are perpendicular to z) are supported for these BCs.", surf2.id_));
}
// Compute the surface normal vectors and make sure they are perpendicular
// to the z-axis
Direction norm1 = surf1.normal({0, 0, 0});
Direction norm2 = surf2.normal({0, 0, 0});
if (std::abs(norm1.z) > FP_PRECISION) {
throw std::invalid_argument(fmt::format("Rotational periodic BCs are only "
"supported for rotations about the z-axis, but surface {} is not "
"perpendicular to the z-axis.", surf1.id_));
}
if (std::abs(norm2.z) > FP_PRECISION) {
throw std::invalid_argument(fmt::format("Rotational periodic BCs are only "
"supported for rotations about the z-axis, but surface {} is not "
"perpendicular to the z-axis.", surf2.id_));
}
// Make sure both surfaces intersect the origin
if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) {
throw std::invalid_argument(fmt::format("Rotational periodic BCs are only "
"supported for rotations about the origin, but surface {} does not "
"intersect the origin.", surf1.id_));
}
if (std::abs(surf2.evaluate({0, 0, 0})) > FP_COINCIDENT) {
throw std::invalid_argument(fmt::format("Rotational periodic BCs are only "
"supported for rotations about the origin, but surface {} does not "
"intersect the origin.", surf2.id_));
}
// Compute the BC rotation angle. Here it is assumed that both surface
// normal vectors point inwards---towards the valid geometry region.
// Consequently, the rotation angle is not the difference between the two
// normals, but is instead the difference between one normal and one
// anti-normal. (An incident ray on one surface must be an outgoing ray on
// the other surface after rotation hence the anti-normal.)
double theta1 = std::atan2(norm1.y, norm1.x);
double theta2 = std::atan2(norm2.y, norm2.x) + PI;
angle_ = theta2 - theta1;
// Warn the user if the angle does not evenly divide a circle
double rem = std::abs(std::remainder((2 * PI / angle_), 1.0));
if (rem > FP_REL_PRECISION && rem < 1 - FP_REL_PRECISION) {
warning(fmt::format("Rotational periodic BC specified with a rotation "
"angle of {} degrees which does not evenly divide 360 degrees.",
angle_ * 180 / PI));
}
}
void
RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const
{
// TODO: off-by-one on surface indices throughout this function.
int i_particle_surf = std::abs(p.surface_) - 1;
// Figure out which of the two BC surfaces were struck to figure out if a
// forward or backward rotation is required. Specify the other surface as
// the particle's new surface.
double theta;
int new_surface;
if (i_particle_surf == i_surf_) {
theta = angle_;
new_surface = p.surface_ > 0 ? -(j_surf_ + 1) : j_surf_ + 1;
} else if (i_particle_surf == j_surf_) {
theta = -angle_;
new_surface = p.surface_ > 0 ? -(i_surf_ + 1) : i_surf_ + 1;
} else {
throw std::runtime_error("Called BoundaryCondition::handle_particle after "
"hitting a surface, but that surface is not recognized by the BC.");
}
// Rotate the particle's position and direction about the z-axis.
Position r = p.r();
Direction u = p.u();
double cos_theta = std::cos(theta);
double sin_theta = std::sin(theta);
Position new_r = {
cos_theta*r.x - sin_theta*r.y,
sin_theta*r.x + cos_theta*r.y,
r.z};
Direction new_u = {
cos_theta*u.x - sin_theta*u.y,
sin_theta*u.x + cos_theta*u.y,
u.z};
// Pass the new location, direction, and surface to the particle.
p.cross_periodic_bc(surf, new_r, new_u, new_surface);
}
} // namespace openmc

11
src/dagmc.cpp
View file

@ -297,12 +297,13 @@ void load_dagmc_geometry()
std::string bc_value = DMD.get_surface_property("boundary", surf_handle);
to_lower(bc_value);
if (bc_value.empty() || bc_value == "transmit" || bc_value == "transmission") {
// set to transmission by default
s->bc_ = Surface::BoundaryType::TRANSMIT;
// Leave the bc_ a nullptr
} else if (bc_value == "vacuum") {
s->bc_ = Surface::BoundaryType::VACUUM;
s->bc_ = std::make_shared<VacuumBC>();
} else if (bc_value == "reflective" || bc_value == "reflect" || bc_value == "reflecting") {
s->bc_ = Surface::BoundaryType::REFLECT;
s->bc_ = std::make_shared<ReflectiveBC>();
} else if (bc_value == "white") {
fatal_error("White boundary condition not supported in DAGMC.");
} else if (bc_value == "periodic") {
fatal_error("Periodic boundary condition not supported in DAGMC.");
} else {
@ -318,7 +319,7 @@ void load_dagmc_geometry()
// if this surface belongs to the graveyard
if (graveyard && parent_vols.find(graveyard) != parent_vols.end()) {
// set graveyard surface BC's to vacuum
s->bc_ = Surface::BoundaryType::VACUUM;
s->bc_ = std::make_shared<VacuumBC>();
}
// add to global array and map

280
src/particle.cpp
View file

@ -413,151 +413,9 @@ Particle::cross_surface()
write_message(1, " Crossing surface {}", surf->id_);
}
if (surf->bc_ == Surface::BoundaryType::VACUUM && (settings::run_mode != RunMode::PLOTTING)) {
// =======================================================================
// PARTICLE LEAKS OUT OF PROBLEM
// Kill particle
alive_ = false;
// Score any surface current tallies -- note that the particle is moved
// forward slightly so that if the mesh boundary is on the surface, it is
// still processed
if (!model::active_meshsurf_tallies.empty()) {
// TODO: Find a better solution to score surface currents than
// physically moving the particle forward slightly
this->r() += TINY_BIT * this->u();
score_surface_tally(*this, model::active_meshsurf_tallies);
}
// Score to global leakage tally
keff_tally_leakage_ += wgt_;
// Display message
if (settings::verbosity >= 10 || trace_) {
write_message(1, " Leaked out of surface {}", surf->id_);
}
return;
} else if ((surf->bc_ == Surface::BoundaryType::REFLECT ||
surf->bc_ == Surface::BoundaryType::WHITE)
&& (settings::run_mode != RunMode::PLOTTING)) {
// =======================================================================
// PARTICLE REFLECTS FROM SURFACE
// Do not handle reflective boundary conditions on lower universes
if (n_coord_ != 1) {
this->mark_as_lost("Cannot reflect particle " + std::to_string(id_) +
" off surface in a lower universe.");
return;
}
// Score surface currents since reflection causes the direction of the
// particle to change. For surface filters, we need to score the tallies
// twice, once before the particle's surface attribute has changed and
// once after. For mesh surface filters, we need to artificially move
// the particle slightly back in case the surface crossing is coincident
// with a mesh boundary
if (!model::active_surface_tallies.empty()) {
score_surface_tally(*this, model::active_surface_tallies);
}
if (!model::active_meshsurf_tallies.empty()) {
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
score_surface_tally(*this, model::active_meshsurf_tallies);
this->r() = r;
}
Direction u = (surf->bc_ == Surface::BoundaryType::REFLECT) ?
surf->reflect(this->r(), this->u(), this) :
surf->diffuse_reflect(this->r(), this->u(), this->current_seed());
// Make sure new particle direction is normalized
this->u() = u / u.norm();
// Reassign particle's cell and surface
coord_[0].cell = cell_last_[n_coord_last_ - 1];
surface_ = -surface_;
// If a reflective surface is coincident with a lattice or universe
// boundary, it is necessary to redetermine the particle's coordinates in
// the lower universes.
// (unless we're using a dagmc model, which has exactly one universe)
if (!settings::dagmc) {
n_coord_ = 1;
if (!find_cell(*this, true)) {
this->mark_as_lost("Couldn't find particle after reflecting from surface "
+ std::to_string(surf->id_) + ".");
return;
}
}
// Set previous coordinate going slightly past surface crossing
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || trace_) {
write_message(1, " Reflected from surface {}", surf->id_);
}
return;
} else if (surf->bc_ == Surface::BoundaryType::PERIODIC && settings::run_mode != RunMode::PLOTTING) {
// =======================================================================
// PERIODIC BOUNDARY
// Do not handle periodic boundary conditions on lower universes
if (n_coord_ != 1) {
this->mark_as_lost("Cannot transfer particle " + std::to_string(id_) +
" across surface in a lower universe. Boundary conditions must be "
"applied to root universe.");
return;
}
// Score surface currents since reflection causes the direction of the
// particle to change -- artificially move the particle slightly back in
// case the surface crossing is coincident with a mesh boundary
if (!model::active_meshsurf_tallies.empty()) {
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
score_surface_tally(*this, model::active_meshsurf_tallies);
this->r() = r;
}
// Get a pointer to the partner periodic surface
auto surf_p = dynamic_cast<PeriodicSurface*>(surf);
auto other = dynamic_cast<PeriodicSurface*>(
model::surfaces[surf_p->i_periodic_].get());
// Adjust the particle's location and direction.
bool rotational = other->periodic_translate(surf_p, this->r(), this->u());
// Reassign particle's surface
// TODO: off-by-one
surface_ = rotational ?
surf_p->i_periodic_ + 1 :
((surface_ > 0) ? surf_p->i_periodic_ + 1 : -(surf_p->i_periodic_ + 1));
// Figure out what cell particle is in now
n_coord_ = 1;
if (!find_cell(*this, true)) {
this->mark_as_lost("Couldn't find particle after hitting periodic "
"boundary on surface " + std::to_string(surf->id_) + ".");
return;
}
// Set previous coordinate going slightly past surface crossing
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || trace_) {
write_message(1, " Hit periodic boundary on surface {}", surf->id_);
}
// Handle any applicable boundary conditions.
if (surf->bc_ && settings::run_mode != RunMode::PLOTTING) {
surf->bc_->handle_particle(*this, *surf);
return;
}
@ -613,6 +471,138 @@ Particle::cross_surface()
}
}
void
Particle::cross_vacuum_bc(const Surface& surf)
{
// Kill the particle
alive_ = false;
// Score any surface current tallies -- note that the particle is moved
// forward slightly so that if the mesh boundary is on the surface, it is
// still processed
if (!model::active_meshsurf_tallies.empty()) {
// TODO: Find a better solution to score surface currents than
// physically moving the particle forward slightly
this->r() += TINY_BIT * this->u();
score_surface_tally(*this, model::active_meshsurf_tallies);
}
// Score to global leakage tally
keff_tally_leakage_ += wgt_;
// Display message
if (settings::verbosity >= 10 || trace_) {
write_message(1, " Leaked out of surface {}", surf.id_);
}
}
void
Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
{
// Do not handle reflective boundary conditions on lower universes
if (n_coord_ != 1) {
this->mark_as_lost("Cannot reflect particle " + std::to_string(id_) +
" off surface in a lower universe.");
return;
}
// Score surface currents since reflection causes the direction of the
// particle to change. For surface filters, we need to score the tallies
// twice, once before the particle's surface attribute has changed and
// once after. For mesh surface filters, we need to artificially move
// the particle slightly back in case the surface crossing is coincident
// with a mesh boundary
if (!model::active_surface_tallies.empty()) {
score_surface_tally(*this, model::active_surface_tallies);
}
if (!model::active_meshsurf_tallies.empty()) {
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
score_surface_tally(*this, model::active_meshsurf_tallies);
this->r() = r;
}
// Set the new particle direction
this->u() = new_u;
// Reassign particle's cell and surface
coord_[0].cell = cell_last_[n_coord_last_ - 1];
surface_ = -surface_;
// If a reflective surface is coincident with a lattice or universe
// boundary, it is necessary to redetermine the particle's coordinates in
// the lower universes.
// (unless we're using a dagmc model, which has exactly one universe)
if (!settings::dagmc) {
n_coord_ = 1;
if (!find_cell(*this, true)) {
this->mark_as_lost("Couldn't find particle after reflecting from surface "
+ std::to_string(surf.id_) + ".");
return;
}
}
// Set previous coordinate going slightly past surface crossing
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || trace_) {
write_message(1, " Reflected from surface {}", surf.id_);
}
}
void
Particle::cross_periodic_bc(const Surface& surf, Position new_r,
Direction new_u, int new_surface)
{
// Do not handle periodic boundary conditions on lower universes
if (n_coord_ != 1) {
this->mark_as_lost("Cannot transfer particle " + std::to_string(id_) +
" across surface in a lower universe. Boundary conditions must be "
"applied to root universe.");
return;
}
// Score surface currents since reflection causes the direction of the
// particle to change -- artificially move the particle slightly back in
// case the surface crossing is coincident with a mesh boundary
if (!model::active_meshsurf_tallies.empty()) {
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
score_surface_tally(*this, model::active_meshsurf_tallies);
this->r() = r;
}
// Adjust the particle's location and direction.
r() = new_r;
u() = new_u;
// Reassign particle's surface
surface_ = new_surface;
// Figure out what cell particle is in now
n_coord_ = 1;
if (!find_cell(*this, true)) {
this->mark_as_lost("Couldn't find particle after hitting periodic "
"boundary on surface " + std::to_string(surf.id_) + ". The normal vector "
"of one periodic surface may need to be reversed.");
return;
}
// Set previous coordinate going slightly past surface crossing
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || trace_) {
write_message(1, " Hit periodic boundary on surface {}", surf.id_);
}
}
void
Particle::mark_as_lost(const char* message)
{

299
src/surface.cpp
View file

@ -3,6 +3,7 @@
#include <array>
#include <cmath>
#include <utility>
#include <set>
#include <fmt/core.h>
#include <gsl/gsl>
@ -126,25 +127,20 @@ Surface::Surface(pugi::xml_node surf_node)
std::string surf_bc = get_node_value(surf_node, "boundary", true, true);
if (surf_bc == "transmission" || surf_bc == "transmit" ||surf_bc.empty()) {
bc_ = BoundaryType::TRANSMIT;
// Leave the bc_ a nullptr
} else if (surf_bc == "vacuum") {
bc_ = BoundaryType::VACUUM;
bc_ = std::make_shared<VacuumBC>();
} else if (surf_bc == "reflective" || surf_bc == "reflect"
|| surf_bc == "reflecting") {
bc_ = BoundaryType::REFLECT;
bc_ = std::make_shared<ReflectiveBC>();
} else if (surf_bc == "white") {
bc_ = BoundaryType::WHITE;
bc_ = std::make_shared<WhiteBC>();
} else if (surf_bc == "periodic") {
bc_ = BoundaryType::PERIODIC;
// periodic BC's are handled separately
} else {
fatal_error(fmt::format("Unknown boundary condition \"{}\" specified "
"on surface {}", surf_bc, id_));
}
} else {
bc_ = BoundaryType::TRANSMIT;
}
}
@ -209,22 +205,10 @@ CSGSurface::to_hdf5(hid_t group_id) const
hid_t surf_group = create_group(group_id, group_name);
switch(bc_) {
case BoundaryType::TRANSMIT :
write_string(surf_group, "boundary_type", "transmission", false);
break;
case BoundaryType::VACUUM :
write_string(surf_group, "boundary_type", "vacuum", false);
break;
case BoundaryType::REFLECT :
write_string(surf_group, "boundary_type", "reflective", false);
break;
case BoundaryType::WHITE :
write_string(surf_group, "boundary_type", "white", false);
break;
case BoundaryType::PERIODIC :
write_string(surf_group, "boundary_type", "periodic", false);
break;
if (bc_) {
write_string(surf_group, "boundary_type", bc_->type(), false);
} else {
write_string(surf_group, "boundary_type", "transmission", false);
}
if (!name_.empty()) {
@ -290,17 +274,6 @@ Direction DAGSurface::reflect(Position r, Direction u, Particle* p) const
void DAGSurface::to_hdf5(hid_t group_id) const {}
#endif
//==============================================================================
// PeriodicSurface implementation
//==============================================================================
PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node)
: CSGSurface {surf_node}
{
if (check_for_node(surf_node, "periodic_surface_id")) {
i_periodic_ = std::stoi(get_node_value(surf_node, "periodic_surface_id"));
}
}
//==============================================================================
// Generic functions for x-, y-, and z-, planes.
@ -322,7 +295,7 @@ axis_aligned_plane_distance(Position r, Direction u, bool coincident, double off
//==============================================================================
SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node)
: PeriodicSurface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, x0_);
}
@ -349,29 +322,6 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other,
Position& r, Direction& u) const
{
Direction other_n = other->normal(r);
if (other_n.x == 1 && other_n.y == 0 && other_n.z == 0) {
r.x = x0_;
return false;
} else {
// Assume the partner is an YPlane (the only supported partner). Use the
// evaluate function to find y0, then adjust position/Direction for
// rotational symmetry.
double y0_ = -other->evaluate({0., 0., 0.});
r.y = r.x - x0_ + y0_;
r.x = x0_;
double ux = u.x;
u.x = -u.y;
u.y = ux;
return true;
}
}
BoundingBox
SurfaceXPlane::bounding_box(bool pos_side) const
{
@ -387,7 +337,7 @@ SurfaceXPlane::bounding_box(bool pos_side) const
//==============================================================================
SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node)
: PeriodicSurface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, y0_);
}
@ -414,30 +364,6 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfaceYPlane::periodic_translate(const PeriodicSurface* other,
Position& r, Direction& u) const
{
Direction other_n = other->normal(r);
if (other_n.x == 0 && other_n.y == 1 && other_n.z == 0) {
// The periodic partner is also aligned along y. Just change the y coord.
r.y = y0_;
return false;
} else {
// Assume the partner is an XPlane (the only supported partner). Use the
// evaluate function to find x0, then adjust position/Direction for rotational
// symmetry.
double x0_ = -other->evaluate({0., 0., 0.});
r.x = r.y - y0_ + x0_;
r.y = y0_;
double ux = u.x;
u.x = u.y;
u.y = -ux;
return true;
}
}
BoundingBox
SurfaceYPlane::bounding_box(bool pos_side) const
{
@ -453,7 +379,7 @@ SurfaceYPlane::bounding_box(bool pos_side) const
//==============================================================================
SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node)
: PeriodicSurface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, z0_);
}
@ -480,14 +406,6 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfaceZPlane::periodic_translate(const PeriodicSurface* other,
Position& r, Direction& u) const
{
// Assume the other plane is aligned along z. Just change the z coord.
r.z = z0_;
return false;
}
BoundingBox
SurfaceZPlane::bounding_box(bool pos_side) const
{
@ -503,7 +421,7 @@ SurfaceZPlane::bounding_box(bool pos_side) const
//==============================================================================
SurfacePlane::SurfacePlane(pugi::xml_node surf_node)
: PeriodicSurface(surf_node)
: CSGSurface(surf_node)
{
read_coeffs(surf_node, id_, A_, B_, C_, D_);
}
@ -541,22 +459,6 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const
{
// This function assumes the other plane shares this plane's normal direction.
// Determine the distance to intersection.
double d = evaluate(r) / (A_*A_ + B_*B_ + C_*C_);
// Move the particle that distance along the normal vector.
r.x -= d * A_;
r.y -= d * B_;
r.z -= d * C_;
return false;
}
//==============================================================================
// Generic functions for x-, y-, and z-, cylinders
//==============================================================================
@ -1114,15 +1016,17 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
void read_surfaces(pugi::xml_node node)
{
// Count the number of surfaces.
// Count the number of surfaces
int n_surfaces = 0;
for (pugi::xml_node surf_node : node.children("surface")) {n_surfaces++;}
if (n_surfaces == 0) {
fatal_error("No surfaces found in geometry.xml!");
}
// Loop over XML surface elements and populate the array.
// Loop over XML surface elements and populate the array. Keep track of
// periodic surfaces.
model::surfaces.reserve(n_surfaces);
std::set<std::pair<int, int>> periodic_pairs;
{
pugi::xml_node surf_node;
int i_surf;
@ -1130,6 +1034,8 @@ void read_surfaces(pugi::xml_node node)
surf_node = surf_node.next_sibling("surface"), i_surf++) {
std::string surf_type = get_node_value(surf_node, "type", true, true);
// Allocate and initialize the new surface
if (surf_type == "x-plane") {
model::surfaces.push_back(std::make_unique<SurfaceXPlane>(surf_node));
@ -1169,10 +1075,26 @@ void read_surfaces(pugi::xml_node node)
} else {
fatal_error(fmt::format("Invalid surface type, \"{}\"", surf_type));
}
// Check for a periodic surface
if (check_for_node(surf_node, "boundary")) {
std::string surf_bc = get_node_value(surf_node, "boundary", true, true);
if (surf_bc == "periodic") {
if (check_for_node(surf_node, "periodic_surface_id")) {
int i_periodic = std::stoi(get_node_value(surf_node,
"periodic_surface_id"));
int lo_id = std::min(model::surfaces.back()->id_, i_periodic);
int hi_id = std::max(model::surfaces.back()->id_, i_periodic);
periodic_pairs.insert({lo_id, hi_id});
} else {
periodic_pairs.insert({model::surfaces.back()->id_, -1});
}
}
}
}
}
// Fill the surface map.
// Fill the surface map
for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
int id = model::surfaces[i_surf]->id_;
auto in_map = model::surface_map.find(id);
@ -1184,105 +1106,62 @@ void read_surfaces(pugi::xml_node node)
}
}
// Find the global bounding box (of periodic BC surfaces).
double xmin {INFTY}, xmax {-INFTY}, ymin {INFTY}, ymax {-INFTY},
zmin {INFTY}, zmax {-INFTY};
int i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax;
for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
if (model::surfaces[i_surf]->bc_ == Surface::BoundaryType::PERIODIC) {
// Downcast to the PeriodicSurface type.
Surface* surf_base = model::surfaces[i_surf].get();
auto surf = dynamic_cast<PeriodicSurface*>(surf_base);
// Make sure this surface inherits from PeriodicSurface.
if (!surf) {
fatal_error(fmt::format(
"Periodic boundary condition not supported for surface {}. Periodic "
"BCs are only supported for planar surfaces.", surf_base->id_));
}
// See if this surface makes part of the global bounding box.
auto bb = surf->bounding_box(true) & surf->bounding_box(false);
if (bb.xmin > -INFTY && bb.xmin < xmin) {
xmin = bb.xmin;
i_xmin = i_surf;
}
if (bb.xmax < INFTY && bb.xmax > xmax) {
xmax = bb.xmax;
i_xmax = i_surf;
}
if (bb.ymin > -INFTY && bb.ymin < ymin) {
ymin = bb.ymin;
i_ymin = i_surf;
}
if (bb.ymax < INFTY && bb.ymax > ymax) {
ymax = bb.ymax;
i_ymax = i_surf;
}
if (bb.zmin > -INFTY && bb.zmin < zmin) {
zmin = bb.zmin;
i_zmin = i_surf;
}
if (bb.zmax < INFTY && bb.zmax > zmax) {
zmax = bb.zmax;
i_zmax = i_surf;
}
// Resolve unpaired periodic surfaces. A lambda function is used with
// std::find_if to identify the unpaired surfaces.
auto is_unresolved_pair =
[](const std::pair<int, int> p){return p.second == -1;};
auto first_unresolved = std::find_if(periodic_pairs.begin(),
periodic_pairs.end(), is_unresolved_pair);
if (first_unresolved != periodic_pairs.end()) {
// Found one unpaired surface; search for a second one
auto next_elem = first_unresolved;
next_elem++;
auto second_unresolved = std::find_if(next_elem, periodic_pairs.end(),
is_unresolved_pair);
if (second_unresolved == periodic_pairs.end()) {
fatal_error("Found only one periodic surface without a specified partner."
" Please specify the partner for each periodic surface.");
}
// Make sure there isn't a third unpaired surface
next_elem = second_unresolved;
next_elem++;
auto third_unresolved = std::find_if(next_elem,
periodic_pairs.end(), is_unresolved_pair);
if (third_unresolved != periodic_pairs.end()) {
fatal_error("Found at least three periodic surfaces without a specified "
"partner. Please specify the partner for each periodic surface.");
}
// Add the completed pair and remove the old, unpaired entries
int lo_id = std::min(first_unresolved->first, second_unresolved->first);
int hi_id = std::max(first_unresolved->first, second_unresolved->first);
periodic_pairs.insert({lo_id, hi_id});
periodic_pairs.erase(first_unresolved);
periodic_pairs.erase(second_unresolved);
}
// Set i_periodic for periodic BC surfaces.
for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
if (model::surfaces[i_surf]->bc_ == Surface::BoundaryType::PERIODIC) {
// Downcast to the PeriodicSurface type.
Surface* surf_base = model::surfaces[i_surf].get();
auto surf = dynamic_cast<PeriodicSurface*>(surf_base);
// Assign the periodic boundary conditions
for (auto periodic_pair : periodic_pairs) {
int i_surf = model::surface_map[periodic_pair.first];
int j_surf = model::surface_map[periodic_pair.second];
Surface& surf1 {*model::surfaces[i_surf]};
Surface& surf2 {*model::surfaces[j_surf]};
// Also try downcasting to the SurfacePlane type (which must be handled
// differently).
SurfacePlane* surf_p = dynamic_cast<SurfacePlane*>(surf);
// Compute the dot product of the surface normals
Direction norm1 = surf1.normal({0, 0, 0});
Direction norm2 = surf2.normal({0, 0, 0});
double dot_prod = norm1.dot(norm2);
if (!surf_p) {
// This is not a SurfacePlane.
if (surf->i_periodic_ == C_NONE) {
// The user did not specify the matching periodic surface. See if we
// can find the partnered surface from the bounding box information.
if (i_surf == i_xmin) {
surf->i_periodic_ = i_xmax;
} else if (i_surf == i_xmax) {
surf->i_periodic_ = i_xmin;
} else if (i_surf == i_ymin) {
surf->i_periodic_ = i_ymax;
} else if (i_surf == i_ymax) {
surf->i_periodic_ = i_ymin;
} else if (i_surf == i_zmin) {
surf->i_periodic_ = i_zmax;
} else if (i_surf == i_zmax) {
surf->i_periodic_ = i_zmin;
} else {
fatal_error("Periodic boundary condition applied to interior "
"surface");
}
} else {
// Convert the surface id to an index.
surf->i_periodic_ = model::surface_map[surf->i_periodic_];
}
} else {
// This is a SurfacePlane. We won't try to find it's partner if the
// user didn't specify one.
if (surf->i_periodic_ == C_NONE) {
fatal_error(fmt::format("No matching periodic surface specified for "
"periodic boundary condition on surface {}", surf->id_));
} else {
// Convert the surface id to an index.
surf->i_periodic_ = model::surface_map[surf->i_periodic_];
}
}
// Make sure the opposite surface is also periodic.
if (model::surfaces[surf->i_periodic_]->bc_ != Surface::BoundaryType::PERIODIC) {
fatal_error(fmt::format("Could not find matching surface for periodic "
"boundary condition on surface {}", surf->id_));
}
// If the dot product is 1 (to within floating point precision) then the
// planes are parallel which indicates a translational periodic boundary
// condition. Otherwise, it is a rotational periodic BC.
if (std::abs(1.0 - dot_prod) < FP_PRECISION) {
surf1.bc_ = std::make_shared<TranslationalPeriodicBC>(i_surf, j_surf);
surf2.bc_ = surf1.bc_;
} else {
surf1.bc_ = std::make_shared<RotationalPeriodicBC>(i_surf, j_surf);
surf2.bc_ = surf1.bc_;
}
}
@ -1290,7 +1169,7 @@ void read_surfaces(pugi::xml_node node)
// surface
bool boundary_exists = false;
for (const auto& surf : model::surfaces) {
if (surf->bc_ != Surface::BoundaryType::TRANSMIT) {
if (surf->bc_) {
boundary_exists = true;
break;
}

18
tests/regression_tests/periodic/inputs_true.dat
View file

@ -1,14 +1,14 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2 3 -4 5 -6 7" universe="0" />
<cell id="2" material="2" region="5 -6 -7" universe="0" />
<surface boundary="periodic" coeffs="-5.0" id="1" periodic_surface_id="2" type="x-plane" />
<surface boundary="periodic" coeffs="5.0" id="2" periodic_surface_id="1" type="x-plane" />
<surface boundary="periodic" coeffs="-5.0" id="3" type="y-plane" />
<surface boundary="periodic" coeffs="5.0" id="4" type="y-plane" />
<surface boundary="reflective" coeffs="-5.0" id="5" type="z-plane" />
<surface boundary="reflective" coeffs="5.0" id="6" type="z-plane" />
<surface coeffs="-2.5 2.5 2.0" id="7" type="z-cylinder" />
<cell id="2" material="2" region="3 5 -6 -7" universe="0" />
<surface boundary="periodic" coeffs="0.0" id="1" periodic_surface_id="3" type="x-plane" />
<surface boundary="reflective" coeffs="5.0" id="2" type="x-plane" />
<surface boundary="periodic" coeffs="0.0" id="3" periodic_surface_id="1" type="y-plane" />
<surface boundary="reflective" coeffs="5.0" id="4" type="y-plane" />
<surface boundary="periodic" coeffs="-5.0" id="5" type="z-plane" />
<surface boundary="periodic" coeffs="0 0 1 5.0" id="6" type="plane" />
<surface coeffs="2.5 0.0 2.0" id="7" type="z-cylinder" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
@ -31,7 +31,7 @@
<inactive>0</inactive>
<source strength="1.0">
<space type="box">
<parameters>-5.0 -5.0 -5.0 5.0 5.0 5.0</parameters>
<parameters>0 0 0 5 5 0</parameters>
</space>
</source>
</settings>

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

@ -1,2 +1,2 @@
k-combined:
1.040109E+00 6.527466E-02
1.638526E+00 6.505124E-04

22
tests/regression_tests/periodic/test.py
View file

@ -21,20 +21,22 @@ class PeriodicTest(PyAPITestHarness):
materials.export_to_xml()
# Define geometry
x_min = openmc.XPlane(surface_id=1, x0=-5., boundary_type='periodic')
x_max = openmc.XPlane(surface_id=2, x0=5., boundary_type='periodic')
x_max.periodic_surface = x_min
x_min = openmc.XPlane(surface_id=1, x0=0., boundary_type='periodic')
x_max = openmc.XPlane(surface_id=2, x0=5., boundary_type='reflective')
y_min = openmc.YPlane(surface_id=3, y0=-5., boundary_type='periodic')
y_max = openmc.YPlane(surface_id=4, y0=5., boundary_type='periodic')
y_min = openmc.YPlane(surface_id=3, y0=0., boundary_type='periodic')
y_max = openmc.YPlane(surface_id=4, y0=5., boundary_type='reflective')
y_min.periodic_surface = x_min
z_min = openmc.ZPlane(surface_id=5, z0=-5., boundary_type='reflective')
z_max = openmc.ZPlane(surface_id=6, z0=5., boundary_type='reflective')
z_cyl = openmc.ZCylinder(surface_id=7, x0=-2.5, y0=2.5, r=2.0)
z_min = openmc.ZPlane(surface_id=5, z0=-5., boundary_type='periodic')
z_max = openmc.Plane(surface_id=6, a=0, b=0, c=1, d=5.,
boundary_type='periodic')
z_cyl = openmc.ZCylinder(surface_id=7, x0=2.5, y0=0., r=2.0)
outside_cyl = openmc.Cell(1, fill=water, region=(
+x_min & -x_max & +y_min & -y_max & +z_min & -z_max & +z_cyl))
inside_cyl = openmc.Cell(2, fill=fuel, region=+z_min & -z_max & -z_cyl)
inside_cyl = openmc.Cell(2, fill=fuel, region=(
+y_min & +z_min & -z_max & -z_cyl))
root_universe = openmc.Universe(0, cells=(outside_cyl, inside_cyl))
geometry = openmc.Geometry()
@ -47,7 +49,7 @@ class PeriodicTest(PyAPITestHarness):
settings.batches = 4
settings.inactive = 0
settings.source = openmc.Source(space=openmc.stats.Box(
*outside_cyl.region.bounding_box))
(0, 0, 0), (5, 5, 0)))
settings.export_to_xml()

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

34
tests/regression_tests/periodic_6fold/inputs_true.dat Normal file
View file

@ -0,0 +1,34 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="9 10 -11 12" universe="0" />
<cell id="2" material="2" region="9 10 -12" universe="0" />
<surface boundary="periodic" coeffs="0.4999999999999999 0.8660254037844387 0.0 0.0" id="9" type="plane" />
<surface boundary="periodic" coeffs="0.4999999999999999 -0.8660254037844387 0.0 0.0" id="10" type="plane" />
<surface boundary="reflective" coeffs="5.0" id="11" type="x-plane" />
<surface coeffs="2.598076211353316 1.4999999999999998 2.0" id="12" type="z-cylinder" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material id="1">
<density units="g/cc" value="1.0" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
<sab name="c_H_in_H2O" />
</material>
<material depletable="true" id="2">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>4</batches>
<inactive>0</inactive>
<source strength="1.0">
<space type="box">
<parameters>0 0 0 5 5 0</parameters>
</space>
</source>
</settings>

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

@ -0,0 +1,2 @@
k-combined:
1.858773E+00 1.550950E-02

62
tests/regression_tests/periodic_6fold/test.py Normal file
View file

@ -0,0 +1,62 @@
import openmc
import numpy as np
from tests.testing_harness import PyAPITestHarness
class Periodic6FoldTest(PyAPITestHarness):
def _build_inputs(self):
# Define materials
water = openmc.Material(1)
water.add_nuclide('H1', 2.0)
water.add_nuclide('O16', 1.0)
water.add_s_alpha_beta('c_H_in_H2O')
water.set_density('g/cc', 1.0)
fuel = openmc.Material(2)
fuel.add_nuclide('U235', 1.0)
fuel.set_density('g/cc', 4.5)
materials = openmc.Materials((water, fuel))
materials.default_temperature = '294K'
materials.export_to_xml()
# Define the geometry. Note that this geometry is somewhat non-sensical
# (it essentially defines a circle of half-cylinders), but it is
# designed so that periodic and reflective BCs will give different
# answers.
theta1 = (-1/6 + 1/2) * np.pi
theta2 = (1/6 - 1/2) * np.pi
plane1 = openmc.Plane(a=np.cos(theta1), b=np.sin(theta1),
boundary_type='periodic')
plane2 = openmc.Plane(a=np.cos(theta2), b=np.sin(theta2),
boundary_type='periodic')
x_max = openmc.XPlane(x0=5., boundary_type='reflective')
z_cyl = openmc.ZCylinder(x0=3*np.cos(np.pi/6), y0=3*np.sin(np.pi/6),
r=2.0)
outside_cyl = openmc.Cell(1, fill=water, region=(
+plane1 & +plane2 & -x_max & +z_cyl))
inside_cyl = openmc.Cell(2, fill=fuel, region=(
+plane1 & +plane2 & -z_cyl))
root_universe = openmc.Universe(0, cells=(outside_cyl, inside_cyl))
geometry = openmc.Geometry()
geometry.root_universe = root_universe
geometry.export_to_xml()
# Define settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 4
settings.inactive = 0
settings.source = openmc.Source(space=openmc.stats.Box(
(0, 0, 0), (5, 5, 0)))
settings.export_to_xml()
def test_periodic():
harness = Periodic6FoldTest('statepoint.4.h5')
harness.main()