From c98142dd239c249330356793c93c09af5dccec76 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 20 Aug 2018 22:45:28 -0400 Subject: [PATCH] Move cross_lattice to C++ --- src/geometry.F90 | 82 +++++------------------------------------------- src/geometry.cpp | 64 ++++++++++++++++++++++++++++++++++++- src/geometry.h | 11 +++++-- src/particle.h | 10 +++++- 4 files changed, 89 insertions(+), 78 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index dcaabb9dda..09f7066bef 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -41,11 +41,18 @@ module geometry function find_cell_c(p, n_search_cells, search_cells) & bind(C, name="find_cell") result(found) import Particle, C_INT, C_BOOL - type(Particle), intent(in) :: p + type(Particle), intent(inout) :: p integer(C_INT), intent(in), value :: n_search_cells integer(C_INT), intent(in), optional :: search_cells(n_search_cells) logical(C_BOOL) :: found end function find_cell_c + + subroutine cross_lattice(p, lattice_translation) & + bind(C, name="cross_lattice") + import Particle, C_INT + type(Particle), intent(inout) :: p + integer(C_INT), intent(in) :: lattice_translation(3) + end subroutine cross_lattice end interface contains @@ -77,79 +84,6 @@ contains end subroutine find_cell -!=============================================================================== -! CROSS_LATTICE moves a particle into a new lattice element -!=============================================================================== - - subroutine cross_lattice(p, lattice_translation) - - type(Particle), intent(inout) :: p - integer, intent(in) :: lattice_translation(3) - integer :: j - integer :: i_xyz(3) ! indices in lattice - logical :: found ! particle found in cell? - class(Lattice), pointer :: lat - - j = p % n_coord - lat => lattices(p % coord(j) % lattice) % obj - - if (verbosity >= 10 .or. trace) then - call write_message(" Crossing lattice " // trim(to_str(lat % id())) & - &// ". Current position (" // trim(to_str(p % coord(j) % lattice_x)) & - &// "," // trim(to_str(p % coord(j) % lattice_y)) // "," & - &// trim(to_str(p % coord(j) % lattice_z)) // ")") - end if - - ! Set the lattice indices. - p % coord(j) % lattice_x = p % coord(j) % lattice_x + lattice_translation(1) - p % coord(j) % lattice_y = p % coord(j) % lattice_y + lattice_translation(2) - p % coord(j) % lattice_z = p % coord(j) % lattice_z + lattice_translation(3) - i_xyz(1) = p % coord(j) % lattice_x - i_xyz(2) = p % coord(j) % lattice_y - i_xyz(3) = p % coord(j) % lattice_z - - ! Set the new coordinate position. - p % coord(j) % xyz = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz) - - OUTSIDE_LAT: if (.not. lat % are_valid_indices(i_xyz)) then - ! The particle is outside the lattice. Search for it from base coord - p % n_coord = 1 - call find_cell(p, found) - if (.not. found) then - if (p % alive) then ! Particle may have been killed in find_cell - call particle_mark_as_lost(p, "Could not locate particle " & - // trim(to_str(p % id)) // " after crossing a lattice boundary.") - return - end if - end if - - else OUTSIDE_LAT - - ! Find cell in next lattice element - p % coord(j) % universe = & - lat % get([i_xyz(1), i_xyz(2), i_xyz(3)]) + 1 - - call find_cell(p, found) - if (.not. found) then - ! In some circumstances, a particle crossing the corner of a cell may - ! not be able to be found in the next universe. In this scenario we cut - ! off all lower-level coordinates and search from universe zero - - ! Remove lower coordinates - p % n_coord = 1 - - ! Search for particle - call find_cell(p, found) - if (.not. found) then - call particle_mark_as_lost(p, "Could not locate particle " // & - trim(to_str(p % id)) // " after crossing a lattice boundary.") - return - end if - end if - end if OUTSIDE_LAT - - end subroutine cross_lattice - !=============================================================================== ! DISTANCE_TO_BOUNDARY calculates the distance to the nearest boundary for a ! particle 'p' traveling in a certain direction. For a cell in a subuniverse diff --git a/src/geometry.cpp b/src/geometry.cpp index e400a2245d..85cefb7e2d 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -1,5 +1,6 @@ #include "geometry.h" +#include #include #include "cell.h" @@ -230,7 +231,7 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) { << lat.id << " but the lattice has no defined outer " "universe."; warning(err_msg); - found = false; + return false; } } @@ -243,4 +244,65 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) { return found; } +//============================================================================== + +extern "C" void +cross_lattice(Particle* p, int lattice_translation[3]) +{ + Lattice& lat {*lattices_c[p->coord[p->n_coord-1].lattice-1]}; + + if (openmc_verbosity >= 10 || openmc_trace) { + std::stringstream msg; + msg << " Crossing lattice " << lat.id << ". Current position (" + << p->coord[p->n_coord-1].lattice_x << "," + << p->coord[p->n_coord-1].lattice_y << "," + << p->coord[p->n_coord-1].lattice_z << ")"; + write_message(msg, 1); + } + + // Set the lattice indices. + p->coord[p->n_coord-1].lattice_x += lattice_translation[0]; + p->coord[p->n_coord-1].lattice_y += lattice_translation[1]; + p->coord[p->n_coord-1].lattice_z += lattice_translation[2]; + std::array i_xyz {p->coord[p->n_coord-1].lattice_x, + p->coord[p->n_coord-1].lattice_y, + p->coord[p->n_coord-1].lattice_z}; + + // Set the new coordinate position. + auto r = lat.get_local_position(p->coord[p->n_coord-2].xyz, i_xyz); + p->coord[p->n_coord-1].xyz[0] = r.x; + p->coord[p->n_coord-1].xyz[1] = r.y; + p->coord[p->n_coord-1].xyz[2] = r.z; + + if (!lat.are_valid_indices(i_xyz)) { + // The particle is outside the lattice. Search for it from the base coords. + p->n_coord = 1; + bool found = find_cell(p, 0, nullptr); + if (!found && p->alive) { + std::stringstream err_msg; + err_msg << "Could not locate particle " << p->id + << " after crossing a lattice boundary"; + p->mark_as_lost(err_msg); + } + + } else { + // Find cell in next lattice element. + p->coord[p->n_coord-1].universe = lat[i_xyz] + 1; + bool found = find_cell(p, 0, nullptr); + + if (!found) { + // A particle crossing the corner of a lattice tile may not be found. In + // this case, search for it from the base coords. + p->n_coord = 1; + bool found = find_cell(p, 0, nullptr); + if (!found && p->alive) { + std::stringstream err_msg; + err_msg << "Could not locate particle " << p->id + << " after crossing a lattice boundary"; + p->mark_as_lost(err_msg); + } + } + } +} + } // namespace openmc diff --git a/src/geometry.h b/src/geometry.h index b73ae314c4..47162c73e1 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -13,19 +13,26 @@ extern "C" int openmc_root_universe; extern std::vector overlap_check_count; //============================================================================== -//! Check for overlapping cells at the particle's position. +//! Check for overlapping cells at a particle's position. //============================================================================== extern "C" bool check_cell_overlap(Particle* p); //============================================================================== -//! Locate the particle in the geometry tree and set its geometry data fields. +//! Locate a particle in the geometry tree and set its geometry data fields. //============================================================================== extern "C" bool find_cell(Particle* p, int n_search_cells, int* search_cells); +//============================================================================== +//! Move a particle into a new lattice tile. +//============================================================================== + +extern "C" void +cross_lattice(Particle* p, int lattice_translation[3]); + } // namespace openmc #endif // OPENMC_GEOMETRY_H diff --git a/src/particle.h b/src/particle.h index b533971710..e03af3d69b 100644 --- a/src/particle.h +++ b/src/particle.h @@ -4,8 +4,10 @@ //! \file particle.h //! \brief Particle type -#include #include +#include +#include +#include #include "openmc.h" @@ -152,6 +154,12 @@ extern "C" { //! \param message A warning message to display void mark_as_lost(const char* message); + void mark_as_lost(const std::string& message) + {mark_as_lost(message.c_str());} + + void mark_as_lost(const std::stringstream& message) + {mark_as_lost(message.str());} + //! create a particle restart HDF5 file void write_restart(); };