From 0f214eb1309ab6cf5f907ab9b1087fb401c02e70 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 17 Oct 2020 20:13:38 -0600 Subject: [PATCH 01/21] Move BC tracking software into separate functions --- include/openmc/particle.h | 12 ++ src/particle.cpp | 273 ++++++++++++++++++++------------------ 2 files changed, 156 insertions(+), 129 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index e551261a75..9168dd6e0a 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,15 @@ public: //! Cross a surface and handle boundary conditions void cross_surface(); + //! Cross a vacuum boundary condition. + void cross_vacuum_bc(const Surface& surf); + + //! Cross a reflective boundary condition. + void cross_reflective_bc(const Surface& surf); + + //! Cross a periodic boundary condition. + void cross_periodic_bc(const Surface& surf); + //! 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/src/particle.cpp b/src/particle.cpp index 0b50671659..2097ba0da3 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -417,147 +417,21 @@ Particle::cross_surface() // ======================================================================= // PARTICLE LEAKS OUT OF PROBLEM - // Kill particle - alive_ = false; + cross_vacuum_bc(*surf); - // 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; - } + cross_reflective_bc(*surf); - // 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 + cross_periodic_bc(*surf); - // 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_); - } return; } @@ -613,6 +487,147 @@ 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) +{ + // 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_); + } +} + +void +Particle::cross_periodic_bc(const Surface& surf) +{ + // 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_); + } +} + void Particle::mark_as_lost(const char* message) { From 05ea9ea585baca72ff80326713a31352bc748f83 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 17 Oct 2020 20:34:41 -0600 Subject: [PATCH 02/21] Move reflective/white BCs to a new class --- CMakeLists.txt | 1 + include/openmc/particle.h | 2 +- include/openmc/surface.h | 2 ++ src/boundary_condition.cpp | 29 +++++++++++++++++++++++++++++ src/particle.cpp | 14 ++++---------- src/surface.cpp | 2 ++ 6 files changed, 39 insertions(+), 11 deletions(-) create mode 100644 src/boundary_condition.cpp 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/include/openmc/particle.h b/include/openmc/particle.h index 9168dd6e0a..48570c4ea9 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -244,7 +244,7 @@ public: void cross_vacuum_bc(const Surface& surf); //! Cross a reflective boundary condition. - void cross_reflective_bc(const Surface& surf); + void cross_reflective_bc(const Surface& surf, Direction new_u); //! Cross a periodic boundary condition. void cross_periodic_bc(const Surface& surf); diff --git a/include/openmc/surface.h b/include/openmc/surface.h index cef5f4f880..010d9dfa28 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" @@ -95,6 +96,7 @@ public: int id_; //!< Unique ID BoundaryType bc_; //!< Boundary condition + std::shared_ptr new_bc_ {nullptr}; std::string name_; //!< User-defined name explicit Surface(pugi::xml_node surf_node); diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp new file mode 100644 index 0000000000..d9c463ca09 --- /dev/null +++ b/src/boundary_condition.cpp @@ -0,0 +1,29 @@ +#include "openmc/boundary_condition.h" +#include "openmc/surface.h" + +namespace openmc { + +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); +} + +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); +} + +//void +//PeriodicBC::handle_particle(Particle& p, const Surface& surf) const +//{ +//} + +} // namespace openmc diff --git a/src/particle.cpp b/src/particle.cpp index 2097ba0da3..ba68ce2eb8 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -414,8 +414,6 @@ Particle::cross_surface() } if (surf->bc_ == Surface::BoundaryType::VACUUM && (settings::run_mode != RunMode::PLOTTING)) { - // ======================================================================= - // PARTICLE LEAKS OUT OF PROBLEM cross_vacuum_bc(*surf); @@ -425,7 +423,7 @@ Particle::cross_surface() surf->bc_ == Surface::BoundaryType::WHITE) && (settings::run_mode != RunMode::PLOTTING)) { - cross_reflective_bc(*surf); + surf->new_bc_->handle_particle(*this, *surf); return; @@ -515,7 +513,7 @@ Particle::cross_vacuum_bc(const Surface& surf) } void -Particle::cross_reflective_bc(const Surface& surf) +Particle::cross_reflective_bc(const Surface& surf, Direction new_u) { // Do not handle reflective boundary conditions on lower universes if (n_coord_ != 1) { @@ -542,12 +540,8 @@ Particle::cross_reflective_bc(const Surface& surf) 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(); + // Set the new particle direction + this->u() = new_u; // Reassign particle's cell and surface coord_[0].cell = cell_last_[n_coord_last_ - 1]; diff --git a/src/surface.cpp b/src/surface.cpp index f2faa3efee..6296e128b4 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -134,8 +134,10 @@ Surface::Surface(pugi::xml_node surf_node) } else if (surf_bc == "reflective" || surf_bc == "reflect" || surf_bc == "reflecting") { bc_ = BoundaryType::REFLECT; + new_bc_ = std::make_shared(); } else if (surf_bc == "white") { bc_ = BoundaryType::WHITE; + new_bc_ = std::make_shared(); } else if (surf_bc == "periodic") { bc_ = BoundaryType::PERIODIC; } else { From 1d07f498612987e6740d1b237db8c2afacc2f990 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 17 Oct 2020 20:59:58 -0600 Subject: [PATCH 03/21] Move all BCs to the new BoundaryCondtion class --- include/openmc/particle.h | 3 ++- src/boundary_condition.cpp | 31 +++++++++++++++++++++++++++---- src/particle.cpp | 29 +++++++---------------------- src/surface.cpp | 2 ++ 4 files changed, 38 insertions(+), 27 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 48570c4ea9..f122455fc7 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -247,7 +247,8 @@ public: void cross_reflective_bc(const Surface& surf, Direction new_u); //! Cross a periodic boundary condition. - void cross_periodic_bc(const Surface& surf); + 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 diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index d9c463ca09..334f40b311 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -3,6 +3,12 @@ namespace openmc { +void +VacuumBC::handle_particle(Particle& p, const Surface& surf) const +{ + p.cross_vacuum_bc(surf); +} + void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const { @@ -21,9 +27,26 @@ WhiteBC::handle_particle(Particle& p, const Surface& surf) const p.cross_reflective_bc(surf, u); } -//void -//PeriodicBC::handle_particle(Particle& p, const Surface& surf) const -//{ -//} +void +PeriodicBC::handle_particle(Particle& p, const Surface& surf) const +{ + // 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()); + + // Compute the new particle location and direction + Position r {p.r()}; + Direction u {p.u()}; + bool rotational = other->periodic_translate(surf_p, r, u); + + // Pick the particle's new surface + // TODO: off-by-one + int new_surface = rotational ? + surf_p->i_periodic_ + 1 : + ((p.surface_ > 0) ? surf_p->i_periodic_ + 1 : -(surf_p->i_periodic_ + 1)); + + p.cross_periodic_bc(surf, r, u, new_surface); +} } // namespace openmc diff --git a/src/particle.cpp b/src/particle.cpp index ba68ce2eb8..a7532f8c98 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -413,23 +413,9 @@ Particle::cross_surface() write_message(1, " Crossing surface {}", surf->id_); } - if (surf->bc_ == Surface::BoundaryType::VACUUM && (settings::run_mode != RunMode::PLOTTING)) { - - cross_vacuum_bc(*surf); - - return; - - } else if ((surf->bc_ == Surface::BoundaryType::REFLECT || - surf->bc_ == Surface::BoundaryType::WHITE) - && (settings::run_mode != RunMode::PLOTTING)) { - + // Handle any applicable boundary conditions. + if (surf->new_bc_ && settings::run_mode != RunMode::PLOTTING) { surf->new_bc_->handle_particle(*this, *surf); - - return; - - } else if (surf->bc_ == Surface::BoundaryType::PERIODIC && settings::run_mode != RunMode::PLOTTING) { - cross_periodic_bc(*surf); - return; } @@ -570,7 +556,8 @@ Particle::cross_reflective_bc(const Surface& surf, Direction new_u) } void -Particle::cross_periodic_bc(const Surface& surf) +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) { @@ -596,13 +583,11 @@ Particle::cross_periodic_bc(const Surface& surf) 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()); + r() = new_r; + u() = new_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)); + surface_ = new_surface; // Figure out what cell particle is in now n_coord_ = 1; diff --git a/src/surface.cpp b/src/surface.cpp index 6296e128b4..6957c85235 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -130,6 +130,7 @@ Surface::Surface(pugi::xml_node surf_node) } else if (surf_bc == "vacuum") { bc_ = BoundaryType::VACUUM; + new_bc_ = std::make_shared(); } else if (surf_bc == "reflective" || surf_bc == "reflect" || surf_bc == "reflecting") { @@ -140,6 +141,7 @@ Surface::Surface(pugi::xml_node surf_node) new_bc_ = std::make_shared(); } else if (surf_bc == "periodic") { bc_ = BoundaryType::PERIODIC; + new_bc_ = std::make_shared(); } else { fatal_error(fmt::format("Unknown boundary condition \"{}\" specified " "on surface {}", surf_bc, id_)); From 58a0aa3297e6b722dd27325773eb0fa6ccfedeae Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 17 Oct 2020 22:16:28 -0600 Subject: [PATCH 04/21] Add surface indices to the PeriodicBC class --- src/boundary_condition.cpp | 12 ++- src/surface.cpp | 167 +++++++++++++++---------------------- 2 files changed, 75 insertions(+), 104 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 334f40b311..d577f6bbf9 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -30,10 +30,16 @@ WhiteBC::handle_particle(Particle& p, const Surface& surf) const void PeriodicBC::handle_particle(Particle& p, const Surface& surf) const { - // Get a pointer to the partner periodic surface + // Get a pointer of the first surface downcast to PeriodicSurface auto surf_p = dynamic_cast(&surf); - auto other = dynamic_cast( - model::surfaces[surf_p->i_periodic_].get()); + + // Get a pointer to the partner PeriodicSurface + auto other = + surf.id_ == model::surfaces[i_surf_]->id_ ? + dynamic_cast(model::surfaces[j_surf_].get()) : + dynamic_cast(model::surfaces[i_surf_].get()); + //auto other = dynamic_cast( + // model::surfaces[surf_p->i_periodic_].get()); // Compute the new particle location and direction Position r {p.r()}; diff --git a/src/surface.cpp b/src/surface.cpp index 6957c85235..3006ea34ba 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -141,7 +141,6 @@ Surface::Surface(pugi::xml_node surf_node) new_bc_ = std::make_shared(); } else if (surf_bc == "periodic") { bc_ = BoundaryType::PERIODIC; - new_bc_ = std::make_shared(); } else { fatal_error(fmt::format("Unknown boundary condition \"{}\" specified " "on surface {}", surf_bc, id_)); @@ -1118,15 +1117,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::vector> periodic_pairs; { pugi::xml_node surf_node; int i_surf; @@ -1134,6 +1135,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)); @@ -1173,10 +1176,24 @@ 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")); + periodic_pairs.push_back({model::surfaces.back()->id_, i_periodic}); + } else { + periodic_pairs.push_back({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); @@ -1188,106 +1205,54 @@ 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; - } + // Remove duplicate permutations from the pairs of periodic surfaces (using a + // functor to compare the pairs) + struct compare_pairs { + bool operator()(const std::pair p1, const std::pair p2) + const { + return p1.first == p2.second; } + }; + compare_pairs compare_f; + auto last = std::unique(periodic_pairs.begin(), periodic_pairs.end(), + compare_f); + periodic_pairs.erase(last, periodic_pairs.end()); + + // Resolve unpaired periodic surfaces + struct is_unresolved_pair { + bool operator()(const std::pair p) const + {return p.second == -1;} + }; + is_unresolved_pair unresolved_f; + auto first_unresolved = std::find_if(periodic_pairs.begin(), + periodic_pairs.end(), unresolved_f); + if (first_unresolved != periodic_pairs.end()) { + auto second_unresolved = std::find_if(first_unresolved+1, + periodic_pairs.end(), unresolved_f); + 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."); + } + auto third_unresolved = std::find_if(second_unresolved+1, + periodic_pairs.end(), unresolved_f); + 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."); + } + first_unresolved->second = second_unresolved->first; + 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); - - // Also try downcasting to the SurfacePlane type (which must be handled - // differently). - SurfacePlane* surf_p = dynamic_cast(surf); - - 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_)); - } - } + // 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]}; + surf1.new_bc_ = std::make_shared(i_surf, j_surf); + surf2.new_bc_ = surf1.new_bc_; + dynamic_cast(&surf1)->i_periodic_ = j_surf; + dynamic_cast(&surf2)->i_periodic_ = i_surf; } // Check to make sure a boundary condition was applied to at least one From 16d4c585c45316cd3ed21f1f74973eb0b403ad58 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 18 Oct 2020 11:23:33 -0600 Subject: [PATCH 05/21] Move periodic translation from Surface class to BC --- include/openmc/boundary_condition.h | 83 +++++++++++++++++++++++++++++ src/boundary_condition.cpp | 83 +++++++++++++++++++++++++++-- src/surface.cpp | 15 +++++- 3 files changed, 176 insertions(+), 5 deletions(-) create mode 100644 include/openmc/boundary_condition.h diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h new file mode 100644 index 0000000000..cbc1575aac --- /dev/null +++ b/include/openmc/boundary_condition.h @@ -0,0 +1,83 @@ +#ifndef OPENMC_BOUNDARY_CONDITION_H +#define OPENMC_BOUNDARY_CONDITION_H + +#include "openmc/position.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +//class BoundaryCondition; +// +//namespace model { +// extern std::vector> boundary_conditions; +//} // namespace model + +class Particle; +class Surface; + +class BoundaryCondition +{ +public: + virtual void + handle_particle(Particle& p, const Surface& surf) const = 0; +}; + +class VacuumBC : public BoundaryCondition +{ +public: + void + handle_particle(Particle& p, const Surface& surf) const override; +}; + +class ReflectiveBC : public BoundaryCondition +{ +public: + void + handle_particle(Particle& p, const Surface& surf) const override; +}; + +class WhiteBC : public BoundaryCondition +{ +public: + void + handle_particle(Particle& p, const Surface& surf) const override; +}; + +class PeriodicBC : public BoundaryCondition +{ +public: + PeriodicBC(int i_surf, int j_surf) + : i_surf_(i_surf), j_surf_(j_surf) + {}; + +protected: + int i_surf_; + int j_surf_; +}; + +class TranslationalPeriodicBC : public PeriodicBC +{ +public: + TranslationalPeriodicBC(int i_surf, int j_surf); + + void + handle_particle(Particle& p, const Surface& surf) const override; + +protected: + Position translation_; +}; + +class RotationalPeriodicBC : public PeriodicBC +{ +public: + RotationalPeriodicBC(int i_surf, int j_surf); + + void + handle_particle(Particle& p, const Surface& surf) const override; +}; + +} // namespace openmc +#endif // OPENMC_BOUNDARY_CONDITION_H diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index d577f6bbf9..1c19eba595 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -1,4 +1,9 @@ #include "openmc/boundary_condition.h" + +#include + +#include + #include "openmc/surface.h" namespace openmc { @@ -27,8 +32,82 @@ WhiteBC::handle_particle(Particle& p, const Surface& surf) const p.cross_reflective_bc(surf, u); } +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_]}; + + // The following blocks will resolve the type of each surface and compute the + // appropriate translation vector + + // Check for a pair of x-planes + if (const auto* xplane1 = dynamic_cast(&surf1)) { + if (const auto* xplane2 = dynamic_cast(&surf2)) { + translation_ = {xplane2->x0_ - xplane1->x0_, 0, 0}; + } else { + throw std::invalid_argument(fmt::format("Invalid pair of periodic " + "surfaces ({} and {}). For a translational periodic BC, both surfaces " + "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_)); + } + + // Check for a pair of y-planes + } else if (const auto* yplane1 = dynamic_cast(&surf1)) { + if (const auto* yplane2 = dynamic_cast(&surf2)) { + translation_ = {0, yplane2->y0_ - yplane1->y0_, 0}; + } else { + throw std::invalid_argument(fmt::format("Invalid pair of periodic " + "surfaces ({} and {}). For a translational periodic BC, both surfaces " + "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_)); + } + + // Check for a pair of z-planes + } else if (const auto* zplane1 = dynamic_cast(&surf1)) { + if (const auto* zplane2 = dynamic_cast(&surf2)) { + translation_ = {0, 0, zplane2->z0_ - zplane1->z0_}; + } else { + throw std::invalid_argument(fmt::format("Invalid pair of periodic " + "surfaces ({} and {}). For a translational periodic BC, both surfaces " + "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_)); + } + } + + // TODO: Check for a pair of general planes + +} + void -PeriodicBC::handle_particle(Particle& p, const Surface& surf) const +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::RotationalPeriodicBC(int i_surf, int j_surf) + : PeriodicBC(i_surf, j_surf) +{ +} + +void +RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const { // Get a pointer of the first surface downcast to PeriodicSurface auto surf_p = dynamic_cast(&surf); @@ -38,8 +117,6 @@ PeriodicBC::handle_particle(Particle& p, const Surface& surf) const surf.id_ == model::surfaces[i_surf_]->id_ ? dynamic_cast(model::surfaces[j_surf_].get()) : dynamic_cast(model::surfaces[i_surf_].get()); - //auto other = dynamic_cast( - // model::surfaces[surf_p->i_periodic_].get()); // Compute the new particle location and direction Position r {p.r()}; diff --git a/src/surface.cpp b/src/surface.cpp index 3006ea34ba..5540cf860a 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1249,8 +1249,19 @@ void read_surfaces(pugi::xml_node node) int j_surf = model::surface_map[periodic_pair.second]; Surface& surf1 {*model::surfaces[i_surf]}; Surface& surf2 {*model::surfaces[j_surf]}; - surf1.new_bc_ = std::make_shared(i_surf, j_surf); - surf2.new_bc_ = surf1.new_bc_; + + Direction norm1 = surf1.normal({0, 0, 0}); + Direction norm2 = surf2.normal({0, 0, 0}); + double dot_prod = norm1.dot(norm2); + + if (std::abs(1.0 - dot_prod) < FP_PRECISION) { + surf1.new_bc_ = std::make_shared(i_surf, j_surf); + surf2.new_bc_ = surf1.new_bc_; + } else { + surf1.new_bc_ = std::make_shared(i_surf, j_surf); + surf2.new_bc_ = surf1.new_bc_; + } + dynamic_cast(&surf1)->i_periodic_ = j_surf; dynamic_cast(&surf2)->i_periodic_ = i_surf; } From 8e6842dfc52dce988b1121039a1491e1e269de11 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 18 Oct 2020 12:05:31 -0600 Subject: [PATCH 06/21] Move rotational periodic BCs to the new class Remove the now unnecessary Surface::periodic_translate function. --- include/openmc/boundary_condition.h | 3 ++ include/openmc/surface.h | 12 ----- src/boundary_condition.cpp | 63 +++++++++++++++++-------- src/surface.cpp | 71 ----------------------------- 4 files changed, 48 insertions(+), 101 deletions(-) diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index cbc1575aac..3dfc2fc747 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -77,6 +77,9 @@ public: void handle_particle(Particle& p, const Surface& surf) const override; + +protected: + double angle_; }; } // namespace openmc diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 010d9dfa28..4de7fd6ac8 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -196,18 +196,6 @@ 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; - }; //============================================================================== diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 1c19eba595..5f1b0d68a5 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -4,6 +4,7 @@ #include +#include "openmc/constants.h" #include "openmc/surface.h" namespace openmc { @@ -104,32 +105,58 @@ TranslationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const 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_]}; + + // The following blocks will resolve the type of each surface and compute the + // appropriate rotation angle + + if (const auto* xplane = dynamic_cast(&surf1)) { + if (const auto* yplane = dynamic_cast(&surf2)) { + angle_ = 0.5 * PI; + } + } else if (const auto* yplane = dynamic_cast(&surf1)) { + if (const auto* xplane = dynamic_cast(&surf2)) { + angle_ = -0.5 * PI; + } + } } void RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const { - // Get a pointer of the first surface downcast to PeriodicSurface - auto surf_p = dynamic_cast(&surf); + // TODO: off-by-one on surface indices throughout this function. + int i_particle_surf = std::abs(p.surface_) - 1; - // Get a pointer to the partner PeriodicSurface - auto other = - surf.id_ == model::surfaces[i_surf_]->id_ ? - dynamic_cast(model::surfaces[j_surf_].get()) : - dynamic_cast(model::surfaces[i_surf_].get()); + // Figure out which of the two BC surfaces were struck to figure out if a + // forward or backward rotation is required. Also specify the 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."); + } - // Compute the new particle location and direction - Position r {p.r()}; - Direction u {p.u()}; - bool rotational = other->periodic_translate(surf_p, r, u); + // Rotate the particle's position and direction about the z-axis. + Position r = p.r(); + Direction u = p.u(); + Position new_r = { + std::cos(theta)*r[0] - std::sin(theta)*r[1], + std::sin(theta)*r[0] + std::cos(theta)*r[1], + r[2]}; + Direction new_u = { + std::cos(theta)*u[0] - std::sin(theta)*u[1], + std::sin(theta)*u[0] + std::cos(theta)*u[1], + u[2]}; - // Pick the particle's new surface - // TODO: off-by-one - int new_surface = rotational ? - surf_p->i_periodic_ + 1 : - ((p.surface_ > 0) ? surf_p->i_periodic_ + 1 : -(surf_p->i_periodic_ + 1)); - - p.cross_periodic_bc(surf, r, u, new_surface); + // 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/surface.cpp b/src/surface.cpp index 5540cf860a..9b270472cc 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -352,29 +352,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 { @@ -417,30 +394,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 { @@ -483,14 +436,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 { @@ -544,22 +489,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 //============================================================================== From 2fa449de5ad92110cdfc93537a50dc9a3d33567f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 18 Oct 2020 12:19:46 -0600 Subject: [PATCH 07/21] Remove the PeriodicSurface class --- include/openmc/surface.h | 32 ++++---------------------------- src/boundary_condition.cpp | 20 ++++++++++++++++++++ src/particle.cpp | 5 ----- src/surface.cpp | 22 ++++------------------ 4 files changed, 28 insertions(+), 51 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 4de7fd6ac8..14fbe047b0 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -143,7 +143,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. @@ -182,21 +181,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); -}; //============================================================================== //! A plane perpendicular to the x-axis. @@ -204,7 +188,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); @@ -212,8 +196,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_; @@ -225,7 +207,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); @@ -233,8 +215,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_; @@ -246,7 +226,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); @@ -254,8 +234,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_; @@ -267,7 +245,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); @@ -275,8 +253,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 index 5f1b0d68a5..a1bd8721a0 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -9,12 +9,20 @@ 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 { @@ -24,6 +32,10 @@ ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const p.cross_reflective_bc(surf, u); } +//============================================================================== +// WhiteBC implementation +//============================================================================== + void WhiteBC::handle_particle(Particle& p, const Surface& surf) const { @@ -33,6 +45,10 @@ WhiteBC::handle_particle(Particle& p, const Surface& surf) const p.cross_reflective_bc(surf, u); } +//============================================================================== +// TranslationalPeriodicBC implementation +//============================================================================== + TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) : PeriodicBC(i_surf, j_surf) { @@ -102,6 +118,10 @@ TranslationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const 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) { diff --git a/src/particle.cpp b/src/particle.cpp index a7532f8c98..b2e4013502 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -577,11 +577,6 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r, 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. r() = new_r; u() = new_u; diff --git a/src/surface.cpp b/src/surface.cpp index 9b270472cc..5eba9d2100 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -293,17 +293,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. @@ -325,7 +314,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_); } @@ -367,7 +356,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_); } @@ -409,7 +398,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_); } @@ -451,7 +440,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_); } @@ -1190,9 +1179,6 @@ void read_surfaces(pugi::xml_node node) surf1.new_bc_ = std::make_shared(i_surf, j_surf); surf2.new_bc_ = surf1.new_bc_; } - - dynamic_cast(&surf1)->i_periodic_ = j_surf; - dynamic_cast(&surf2)->i_periodic_ = i_surf; } // Check to make sure a boundary condition was applied to at least one From 593a4addcfe07fdb403d120a971608ddbe51922f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 18 Oct 2020 12:33:36 -0600 Subject: [PATCH 08/21] Remove the old Surface::bc_ attribute --- include/openmc/boundary_condition.h | 16 ++++++---- include/openmc/surface.h | 16 ++-------- src/particle.cpp | 4 +-- src/surface.cpp | 48 +++++++++-------------------- 4 files changed, 29 insertions(+), 55 deletions(-) diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index 3dfc2fc747..7524856c46 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -9,12 +9,6 @@ namespace openmc { // Global variables //============================================================================== -//class BoundaryCondition; -// -//namespace model { -// extern std::vector> boundary_conditions; -//} // namespace model - class Particle; class Surface; @@ -23,6 +17,8 @@ class BoundaryCondition public: virtual void handle_particle(Particle& p, const Surface& surf) const = 0; + + virtual std::string type() const = 0; }; class VacuumBC : public BoundaryCondition @@ -30,6 +26,8 @@ class VacuumBC : public BoundaryCondition public: void handle_particle(Particle& p, const Surface& surf) const override; + + std::string type() const override {return "vacuum";} }; class ReflectiveBC : public BoundaryCondition @@ -37,6 +35,8 @@ class ReflectiveBC : public BoundaryCondition public: void handle_particle(Particle& p, const Surface& surf) const override; + + std::string type() const override {return "reflective";} }; class WhiteBC : public BoundaryCondition @@ -44,6 +44,8 @@ class WhiteBC : public BoundaryCondition public: void handle_particle(Particle& p, const Surface& surf) const override; + + std::string type() const override {return "white";} }; class PeriodicBC : public BoundaryCondition @@ -53,6 +55,8 @@ public: : i_surf_(i_surf), j_surf_(j_surf) {}; + std::string type() const override {return "periodic";} + protected: int i_surf_; int j_surf_; diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 14fbe047b0..1e16c4919b 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -85,19 +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::shared_ptr new_bc_ {nullptr}; - 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(); diff --git a/src/particle.cpp b/src/particle.cpp index b2e4013502..b2902238d7 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -414,8 +414,8 @@ Particle::cross_surface() } // Handle any applicable boundary conditions. - if (surf->new_bc_ && settings::run_mode != RunMode::PLOTTING) { - surf->new_bc_->handle_particle(*this, *surf); + if (surf->bc_ && settings::run_mode != RunMode::PLOTTING) { + surf->bc_->handle_particle(*this, *surf); return; } diff --git a/src/surface.cpp b/src/surface.cpp index 5eba9d2100..af3b0aa371 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -126,28 +126,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; - new_bc_ = std::make_shared(); - + bc_ = std::make_shared(); } else if (surf_bc == "reflective" || surf_bc == "reflect" || surf_bc == "reflecting") { - bc_ = BoundaryType::REFLECT; - new_bc_ = std::make_shared(); + bc_ = std::make_shared(); } else if (surf_bc == "white") { - bc_ = BoundaryType::WHITE; - new_bc_ = std::make_shared(); + 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; } } @@ -212,22 +204,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()) { @@ -1173,11 +1153,11 @@ void read_surfaces(pugi::xml_node node) double dot_prod = norm1.dot(norm2); if (std::abs(1.0 - dot_prod) < FP_PRECISION) { - surf1.new_bc_ = std::make_shared(i_surf, j_surf); - surf2.new_bc_ = surf1.new_bc_; + surf1.bc_ = std::make_shared(i_surf, j_surf); + surf2.bc_ = surf1.bc_; } else { - surf1.new_bc_ = std::make_shared(i_surf, j_surf); - surf2.new_bc_ = surf1.new_bc_; + surf1.bc_ = std::make_shared(i_surf, j_surf); + surf2.bc_ = surf1.bc_; } } @@ -1185,7 +1165,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; } From fe0cffebd820e018461d8bdf83f51eedcec12749 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 18 Oct 2020 14:05:21 -0600 Subject: [PATCH 09/21] Reimplement SurfacePlane translational BCs --- src/boundary_condition.cpp | 78 +++++++++++++++++++++++--------------- 1 file changed, 48 insertions(+), 30 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index a1bd8721a0..e671ad554c 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -58,39 +58,57 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) // The following blocks will resolve the type of each surface and compute the // appropriate translation vector - // Check for a pair of x-planes - if (const auto* xplane1 = dynamic_cast(&surf1)) { - if (const auto* xplane2 = dynamic_cast(&surf2)) { - translation_ = {xplane2->x0_ - xplane1->x0_, 0, 0}; - } else { - throw std::invalid_argument(fmt::format("Invalid pair of periodic " - "surfaces ({} and {}). For a translational periodic BC, both surfaces " - "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_)); - } - - // Check for a pair of y-planes - } else if (const auto* yplane1 = dynamic_cast(&surf1)) { - if (const auto* yplane2 = dynamic_cast(&surf2)) { - translation_ = {0, yplane2->y0_ - yplane1->y0_, 0}; - } else { - throw std::invalid_argument(fmt::format("Invalid pair of periodic " - "surfaces ({} and {}). For a translational periodic BC, both surfaces " - "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_)); - } - - // Check for a pair of z-planes - } else if (const auto* zplane1 = dynamic_cast(&surf1)) { - if (const auto* zplane2 = dynamic_cast(&surf2)) { - translation_ = {0, 0, zplane2->z0_ - zplane1->z0_}; - } else { - throw std::invalid_argument(fmt::format("Invalid pair of periodic " - "surfaces ({} and {}). For a translational periodic BC, both surfaces " - "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_)); - } + // 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_)); } - // TODO: Check for a pair of general planes + // 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 From b64a3d3cb7203d791836a3527c6a9ab2237857f0 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 18 Oct 2020 15:46:38 -0600 Subject: [PATCH 10/21] Add general plane support for rotational BCs --- include/openmc/boundary_condition.h | 1 + src/boundary_condition.cpp | 107 +++++++++++++++++++++++----- 2 files changed, 89 insertions(+), 19 deletions(-) diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index 7524856c46..caf83b6117 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -84,6 +84,7 @@ public: protected: double angle_; + bool aligned_normals_; }; } // namespace openmc diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index e671ad554c..67496187f6 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -5,6 +5,7 @@ #include #include "openmc/constants.h" +#include "openmc/error.h" #include "openmc/surface.h" namespace openmc { @@ -55,9 +56,6 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) Surface& surf1 {*model::surfaces[i_surf_]}; Surface& surf2 {*model::surfaces[j_surf_]}; - // The following blocks will resolve the type of each surface and compute the - // appropriate translation vector - // Make sure the first surface has an appropriate type. if (const auto* ptr = dynamic_cast(&surf1)) { } else if (const auto* ptr = dynamic_cast(&surf1)) { @@ -146,18 +144,81 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf) Surface& surf1 {*model::surfaces[i_surf_]}; Surface& surf2 {*model::surfaces[j_surf_]}; - // The following blocks will resolve the type of each surface and compute the - // appropriate rotation angle - - if (const auto* xplane = dynamic_cast(&surf1)) { - if (const auto* yplane = dynamic_cast(&surf2)) { - angle_ = 0.5 * PI; - } - } else if (const auto* yplane = dynamic_cast(&surf1)) { - if (const auto* xplane = dynamic_cast(&surf2)) { - angle_ = -0.5 * PI; - } + // 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 angle between the two surfaces; this is the BC rotation angle + double theta1 = std::atan2(norm1.y, norm1.x); + double theta2 = std::atan2(norm2.y, norm2.x); + 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)); + } + + // Guess whether or not the normal vectors of the two planes are aligned, i.e. + // if an arc passing from one surface through the geometry to the other + // surface will pass through the same "side" of both surfaces. If the user + // specified an x-plane and a y-plane then the geometry likely lies in the + // first quadrant which means the normals are not aligned. Otherwise, assume + // the opposite. + aligned_normals_ = !(surf1_is_xyplane && surf2_is_xyplane); } void @@ -167,20 +228,28 @@ RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const 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. Also specify the new surface. + // 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; + 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."); } + // If the normal vectors of the two surfaces are not aligned, then the logic + // must be reversed for rotation and picking a new surface halfspace. + if (not aligned_normals_) { + theta = -theta; + new_surface = -new_surface; + } + // Rotate the particle's position and direction about the z-axis. Position r = p.r(); Direction u = p.u(); From b403f877c24779070dae9c0e1109a739937f9e03 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 18 Oct 2020 16:20:44 -0600 Subject: [PATCH 11/21] Increase periodic BC test coverage --- .../regression_tests/periodic/inputs_true.dat | 26 ++++++++++----- .../periodic/results_true.dat | 2 +- tests/regression_tests/periodic/test.py | 33 +++++++++++++------ 3 files changed, 41 insertions(+), 20 deletions(-) diff --git a/tests/regression_tests/periodic/inputs_true.dat b/tests/regression_tests/periodic/inputs_true.dat index 6f613a1602..ab25560f81 100644 --- a/tests/regression_tests/periodic/inputs_true.dat +++ b/tests/regression_tests/periodic/inputs_true.dat @@ -1,14 +1,14 @@ - - - - - - - - + + + + + + + + @@ -31,7 +31,15 @@ 0 - -5.0 -5.0 -5.0 5.0 5.0 5.0 + 0 0 0 5 5 0 + + + + 2.5 2.5 0 + 6 11 + 300 300 + + 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..22d38ef271 100644 --- a/tests/regression_tests/periodic/test.py +++ b/tests/regression_tests/periodic/test.py @@ -21,33 +21,46 @@ 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() geometry.root_universe = root_universe geometry.export_to_xml() + plots = openmc.Plots() + + plot = openmc.Plot() + plot.basis = 'xz' + plot.origin = (2.5, 2.5, 0) + plot.width = (6, 11) + plot.pixels = (300, 300) + plots.append(plot) + + plots.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( - *outside_cyl.region.bounding_box)) + (0, 0, 0), (5, 5, 0))) settings.export_to_xml() From f060b4eb72a2952e2a1b6b7529953611b58d6857 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 18 Oct 2020 16:39:19 -0600 Subject: [PATCH 12/21] Add docstrings to boundary_condition.h --- include/openmc/boundary_condition.h | 45 ++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 4 deletions(-) diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index caf83b6117..d40c645695 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -5,22 +5,33 @@ namespace openmc { -//============================================================================== -// Global variables -//============================================================================== - +// 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: @@ -30,6 +41,10 @@ public: std::string type() const override {return "vacuum";} }; +//============================================================================== +//! A BC that returns particles via specular reflection. +//============================================================================== + class ReflectiveBC : public BoundaryCondition { public: @@ -39,6 +54,10 @@ public: std::string type() const override {return "reflective";} }; +//============================================================================== +//! A BC that returns particles via diffuse reflection. +//============================================================================== + class WhiteBC : public BoundaryCondition { public: @@ -48,6 +67,10 @@ public: std::string type() const override {return "white";} }; +//============================================================================== +//! A BC that moves particles to another part of the problem. +//============================================================================== + class PeriodicBC : public BoundaryCondition { public: @@ -62,6 +85,10 @@ protected: int j_surf_; }; +//============================================================================== +//! A BC that moves particles to another part of the problem without rotation. +//============================================================================== + class TranslationalPeriodicBC : public PeriodicBC { public: @@ -71,9 +98,16 @@ public: 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: @@ -83,7 +117,10 @@ public: handle_particle(Particle& p, const Surface& surf) const override; protected: + //! Angle about the axis by which particle coordinates will be rotated double angle_; + + //! Whether or not the boundary normal vector reverses after a rotation bool aligned_normals_; }; From fb8bae8e9e98d103907426455c3bb07eca08412a Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 18 Oct 2020 18:43:04 -0600 Subject: [PATCH 13/21] Update BCs for DAGMC surfaces --- src/dagmc.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) 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 From 1af80b1d39abec972fb25818a2cc8674c73b490c Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 26 Oct 2020 21:31:26 -0600 Subject: [PATCH 14/21] Apply suggestions from code review Co-authored-by: Paul Romano --- include/openmc/boundary_condition.h | 23 ++++++++--------------- src/boundary_condition.cpp | 18 ++++++++++-------- 2 files changed, 18 insertions(+), 23 deletions(-) diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h index d40c645695..c03e19e50f 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -13,13 +13,12 @@ class Surface; //! A class that tells particles what to do after they strike an outer boundary. //============================================================================== -class BoundaryCondition -{ +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. + //! 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; @@ -32,8 +31,7 @@ public: //! A BC that kills particles, indicating they left the problem. //============================================================================== -class VacuumBC : public BoundaryCondition -{ +class VacuumBC : public BoundaryCondition { public: void handle_particle(Particle& p, const Surface& surf) const override; @@ -45,8 +43,7 @@ public: //! A BC that returns particles via specular reflection. //============================================================================== -class ReflectiveBC : public BoundaryCondition -{ +class ReflectiveBC : public BoundaryCondition { public: void handle_particle(Particle& p, const Surface& surf) const override; @@ -58,8 +55,7 @@ public: //! A BC that returns particles via diffuse reflection. //============================================================================== -class WhiteBC : public BoundaryCondition -{ +class WhiteBC : public BoundaryCondition { public: void handle_particle(Particle& p, const Surface& surf) const override; @@ -71,8 +67,7 @@ public: //! A BC that moves particles to another part of the problem. //============================================================================== -class PeriodicBC : public BoundaryCondition -{ +class PeriodicBC : public BoundaryCondition { public: PeriodicBC(int i_surf, int j_surf) : i_surf_(i_surf), j_surf_(j_surf) @@ -89,8 +84,7 @@ protected: //! A BC that moves particles to another part of the problem without rotation. //============================================================================== -class TranslationalPeriodicBC : public PeriodicBC -{ +class TranslationalPeriodicBC : public PeriodicBC { public: TranslationalPeriodicBC(int i_surf, int j_surf); @@ -108,8 +102,7 @@ protected: //! Currently only rotations about the z-axis are supported. //============================================================================== -class RotationalPeriodicBC : public PeriodicBC -{ +class RotationalPeriodicBC : public PeriodicBC { public: RotationalPeriodicBC(int i_surf, int j_surf); diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 67496187f6..bc01c21994 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -190,12 +190,12 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf) // 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 " + "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 " + "supported for rotations about the origin, but surface {} does not " "intersect the origin.", surf2.id_)); } @@ -253,14 +253,16 @@ RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const // 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 = { - std::cos(theta)*r[0] - std::sin(theta)*r[1], - std::sin(theta)*r[0] + std::cos(theta)*r[1], - r[2]}; + cos_theta*r.x - sin_theta*r.y, + sin_theta*r.x + cos_theta*r.y, + r.z}; Direction new_u = { - std::cos(theta)*u[0] - std::sin(theta)*u[1], - std::sin(theta)*u[0] + std::cos(theta)*u[1], - u[2]}; + 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); From 5fb806bcb5ac222cd4e574b42809b69be7ca08e2 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 26 Oct 2020 22:15:20 -0600 Subject: [PATCH 15/21] Add comments to periodic BC implementation --- include/openmc/particle.h | 14 ++++++++++++++ src/surface.cpp | 4 ++++ 2 files changed, 18 insertions(+) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index f122455fc7..f70faa611d 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -241,12 +241,26 @@ public: 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); diff --git a/src/surface.cpp b/src/surface.cpp index af3b0aa371..d3d0cb7655 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1148,10 +1148,14 @@ void read_surfaces(pugi::xml_node node) Surface& surf1 {*model::surfaces[i_surf]}; Surface& surf2 {*model::surfaces[j_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 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_; From 803ae8d02d7964e14dcd8ff553ad2620eb61301b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 26 Oct 2020 22:31:27 -0600 Subject: [PATCH 16/21] Simplify pairing logic for periodic surfaces --- src/surface.cpp | 56 ++++++++++++++++++++++++------------------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/src/surface.cpp b/src/surface.cpp index d3d0cb7655..e0f2787b6e 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -1025,7 +1026,7 @@ void read_surfaces(pugi::xml_node node) // Loop over XML surface elements and populate the array. Keep track of // periodic surfaces. model::surfaces.reserve(n_surfaces); - std::vector> periodic_pairs; + std::set> periodic_pairs; { pugi::xml_node surf_node; int i_surf; @@ -1082,9 +1083,11 @@ void read_surfaces(pugi::xml_node node) if (check_for_node(surf_node, "periodic_surface_id")) { int i_periodic = std::stoi(get_node_value(surf_node, "periodic_surface_id")); - periodic_pairs.push_back({model::surfaces.back()->id_, i_periodic}); + 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.push_back({model::surfaces.back()->id_, -1}); + periodic_pairs.insert({model::surfaces.back()->id_, -1}); } } } @@ -1103,41 +1106,38 @@ void read_surfaces(pugi::xml_node node) } } - // Remove duplicate permutations from the pairs of periodic surfaces (using a - // functor to compare the pairs) - struct compare_pairs { - bool operator()(const std::pair p1, const std::pair p2) - const { - return p1.first == p2.second; - } - }; - compare_pairs compare_f; - auto last = std::unique(periodic_pairs.begin(), periodic_pairs.end(), - compare_f); - periodic_pairs.erase(last, periodic_pairs.end()); - - // Resolve unpaired periodic surfaces - struct is_unresolved_pair { - bool operator()(const std::pair p) const - {return p.second == -1;} - }; - is_unresolved_pair unresolved_f; + // 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(), unresolved_f); + periodic_pairs.end(), is_unresolved_pair); if (first_unresolved != periodic_pairs.end()) { - auto second_unresolved = std::find_if(first_unresolved+1, - periodic_pairs.end(), unresolved_f); + // 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."); } - auto third_unresolved = std::find_if(second_unresolved+1, - periodic_pairs.end(), unresolved_f); + + // 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."); } - first_unresolved->second = second_unresolved->first; + + // 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); } From c75790d3ac0d743eaefc0b142aa60bb7995cf79a Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 26 Oct 2020 22:32:57 -0600 Subject: [PATCH 17/21] Remove unneeded plot from periodic BC test --- tests/regression_tests/periodic/test.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/tests/regression_tests/periodic/test.py b/tests/regression_tests/periodic/test.py index 22d38ef271..f9d26392f2 100644 --- a/tests/regression_tests/periodic/test.py +++ b/tests/regression_tests/periodic/test.py @@ -43,17 +43,6 @@ class PeriodicTest(PyAPITestHarness): geometry.root_universe = root_universe geometry.export_to_xml() - plots = openmc.Plots() - - plot = openmc.Plot() - plot.basis = 'xz' - plot.origin = (2.5, 2.5, 0) - plot.width = (6, 11) - plot.pixels = (300, 300) - plots.append(plot) - - plots.export_to_xml() - # Define settings settings = openmc.Settings() settings.particles = 1000 From a8ecb795a337ac311f0de6202e12305438314489 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 14 Nov 2020 14:00:16 -0700 Subject: [PATCH 18/21] Assume un-aligned normals for rotational BCs --- docs/source/usersguide/geometry.rst | 24 +++++++++++-------- include/openmc/boundary_condition.h | 3 --- src/boundary_condition.cpp | 23 ++++-------------- src/particle.cpp | 3 ++- .../regression_tests/periodic/inputs_true.dat | 8 ------- 5 files changed, 20 insertions(+), 41 deletions(-) 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 index c03e19e50f..7eb93a61d4 100644 --- a/include/openmc/boundary_condition.h +++ b/include/openmc/boundary_condition.h @@ -112,9 +112,6 @@ public: protected: //! Angle about the axis by which particle coordinates will be rotated double angle_; - - //! Whether or not the boundary normal vector reverses after a rotation - bool aligned_normals_; }; } // namespace openmc diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index bc01c21994..2b88e4c7ca 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -211,14 +211,6 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf) "angle of {} degrees which does not evenly divide 360 degrees.", angle_ * 180 / PI)); } - - // Guess whether or not the normal vectors of the two planes are aligned, i.e. - // if an arc passing from one surface through the geometry to the other - // surface will pass through the same "side" of both surfaces. If the user - // specified an x-plane and a y-plane then the geometry likely lies in the - // first quadrant which means the normals are not aligned. Otherwise, assume - // the opposite. - aligned_normals_ = !(surf1_is_xyplane && surf2_is_xyplane); } void @@ -233,23 +225,16 @@ RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const 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); + 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."); } - // If the normal vectors of the two surfaces are not aligned, then the logic - // must be reversed for rotation and picking a new surface halfspace. - if (not aligned_normals_) { - theta = -theta; - new_surface = -new_surface; - } - // Rotate the particle's position and direction about the z-axis. Position r = p.r(); Direction u = p.u(); diff --git a/src/particle.cpp b/src/particle.cpp index b2902238d7..4453701be4 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -589,7 +589,8 @@ Particle::cross_periodic_bc(const Surface& surf, Position new_r, if (!find_cell(*this, true)) { this->mark_as_lost("Couldn't find particle after hitting periodic " - "boundary on surface " + std::to_string(surf.id_) + "."); + "boundary on surface " + std::to_string(surf.id_) + ". The normal vector " + "of one periodic surface may need to be reversed."); return; } diff --git a/tests/regression_tests/periodic/inputs_true.dat b/tests/regression_tests/periodic/inputs_true.dat index ab25560f81..c1f6d15631 100644 --- a/tests/regression_tests/periodic/inputs_true.dat +++ b/tests/regression_tests/periodic/inputs_true.dat @@ -35,11 +35,3 @@ - - - - 2.5 2.5 0 - 6 11 - 300 300 - - From f8c7a74d6cc4f2222c20cc6acc6238b53ac8545e Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 18 Dec 2020 23:20:04 -0700 Subject: [PATCH 19/21] Fix non-90-degree rotational periodic BCs --- src/boundary_condition.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 2b88e4c7ca..fccae191b2 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -199,9 +199,14 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf) "intersect the origin.", surf2.id_)); } - // Compute the angle between the two surfaces; this is the BC rotation angle + // 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); + double theta2 = std::atan2(norm2.y, norm2.x) + PI; angle_ = theta2 - theta1; // Warn the user if the angle does not evenly divide a circle @@ -225,10 +230,10 @@ RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const double theta; int new_surface; if (i_particle_surf == i_surf_) { - theta = -angle_; + theta = angle_; new_surface = p.surface_ > 0 ? -(j_surf_ + 1) : j_surf_ + 1; } else if (i_particle_surf == j_surf_) { - theta = angle_; + theta = -angle_; new_surface = p.surface_ > 0 ? -(i_surf_ + 1) : i_surf_ + 1; } else { throw std::runtime_error("Called BoundaryCondition::handle_particle after " From 6ed023c790a42b79edda376e55b47e0206c7e734 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 18 Dec 2020 23:20:50 -0700 Subject: [PATCH 20/21] Add a test for non-90-degree periodic BCs --- .../periodic_6fold/inputs_true.dat | 34 ++++++++++ .../periodic_6fold/results_true.dat | 2 + tests/regression_tests/periodic_6fold/test.py | 62 +++++++++++++++++++ 3 files changed, 98 insertions(+) create mode 100644 tests/regression_tests/periodic_6fold/inputs_true.dat create mode 100644 tests/regression_tests/periodic_6fold/results_true.dat create mode 100644 tests/regression_tests/periodic_6fold/test.py 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() From b03b98317c5b72211805d015bbbb3e47484ce671 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 19 Dec 2020 08:35:17 -0700 Subject: [PATCH 21/21] Add missing test __init__.py --- tests/regression_tests/periodic_6fold/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 tests/regression_tests/periodic_6fold/__init__.py 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