diff --git a/CMakeLists.txt b/CMakeLists.txt index a59e6f7e3a..5e99b6151e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/docs/source/usersguide/geometry.rst b/docs/source/usersguide/geometry.rst index c929503df1..18bdf94b5f 100644 --- a/docs/source/usersguide/geometry.rst +++ b/docs/source/usersguide/geometry.rst @@ -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: diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h new file mode 100644 index 0000000000..7eb93a61d4 --- /dev/null +++ b/include/openmc/boundary_condition.h @@ -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 diff --git a/include/openmc/particle.h b/include/openmc/particle.h index e551261a75..f70faa611d 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -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& 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); diff --git a/include/openmc/surface.h b/include/openmc/surface.h index cef5f4f880..1e16c4919b 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -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 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_; }; diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp new file mode 100644 index 0000000000..fccae191b2 --- /dev/null +++ b/src/boundary_condition.cpp @@ -0,0 +1,261 @@ +#include "openmc/boundary_condition.h" + +#include + +#include + +#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(&surf1)) { + } else if (const auto* ptr = dynamic_cast(&surf1)) { + } else if (const auto* ptr = dynamic_cast(&surf1)) { + } else if (const auto* ptr = dynamic_cast(&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(&surf2)) { + } else if (const auto* ptr = dynamic_cast(&surf2)) { + } else if (const auto* ptr = dynamic_cast(&surf2)) { + } else if (const auto* ptr = dynamic_cast(&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(&surf1)) { + surf1_is_xyplane = true; + } else if (const auto* ptr = dynamic_cast(&surf1)) { + surf1_is_xyplane = true; + } else if (const auto* ptr = dynamic_cast(&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(&surf2)) { + surf2_is_xyplane = true; + } else if (const auto* ptr = dynamic_cast(&surf2)) { + surf2_is_xyplane = true; + } else if (const auto* ptr = dynamic_cast(&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 diff --git a/src/dagmc.cpp b/src/dagmc.cpp index f0bd787c47..5c614e610f 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -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(); } else if (bc_value == "reflective" || bc_value == "reflect" || bc_value == "reflecting") { - s->bc_ = Surface::BoundaryType::REFLECT; + s->bc_ = std::make_shared(); + } 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(); } // add to global array and map diff --git a/src/particle.cpp b/src/particle.cpp index 0b50671659..4453701be4 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -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(surf); - auto other = dynamic_cast( - 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) { diff --git a/src/surface.cpp b/src/surface.cpp index f2faa3efee..e0f2787b6e 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -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(); } else if (surf_bc == "reflective" || surf_bc == "reflect" || surf_bc == "reflecting") { - bc_ = BoundaryType::REFLECT; + bc_ = std::make_shared(); } else if (surf_bc == "white") { - bc_ = BoundaryType::WHITE; + bc_ = std::make_shared(); } 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> 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(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(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 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(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(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(i_surf, j_surf); + surf2.bc_ = surf1.bc_; + } else { + surf1.bc_ = std::make_shared(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; } diff --git a/tests/regression_tests/periodic/inputs_true.dat b/tests/regression_tests/periodic/inputs_true.dat index 6f613a1602..c1f6d15631 100644 --- a/tests/regression_tests/periodic/inputs_true.dat +++ b/tests/regression_tests/periodic/inputs_true.dat @@ -1,14 +1,14 @@ - - - - - - - - + + + + + + + + @@ -31,7 +31,7 @@ 0 - -5.0 -5.0 -5.0 5.0 5.0 5.0 + 0 0 0 5 5 0 diff --git a/tests/regression_tests/periodic/results_true.dat b/tests/regression_tests/periodic/results_true.dat index 5a9820dd35..e2994ee2c4 100644 --- a/tests/regression_tests/periodic/results_true.dat +++ b/tests/regression_tests/periodic/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.040109E+00 6.527466E-02 +1.638526E+00 6.505124E-04 diff --git a/tests/regression_tests/periodic/test.py b/tests/regression_tests/periodic/test.py index 879bd84809..f9d26392f2 100644 --- a/tests/regression_tests/periodic/test.py +++ b/tests/regression_tests/periodic/test.py @@ -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() diff --git a/tests/regression_tests/periodic_6fold/__init__.py b/tests/regression_tests/periodic_6fold/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/periodic_6fold/inputs_true.dat b/tests/regression_tests/periodic_6fold/inputs_true.dat new file mode 100644 index 0000000000..46c2a6ff55 --- /dev/null +++ b/tests/regression_tests/periodic_6fold/inputs_true.dat @@ -0,0 +1,34 @@ + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 4 + 0 + + + 0 0 0 5 5 0 + + + diff --git a/tests/regression_tests/periodic_6fold/results_true.dat b/tests/regression_tests/periodic_6fold/results_true.dat new file mode 100644 index 0000000000..2f0874c002 --- /dev/null +++ b/tests/regression_tests/periodic_6fold/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.858773E+00 1.550950E-02 diff --git a/tests/regression_tests/periodic_6fold/test.py b/tests/regression_tests/periodic_6fold/test.py new file mode 100644 index 0000000000..1f8efdbf71 --- /dev/null +++ b/tests/regression_tests/periodic_6fold/test.py @@ -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()