diff --git a/CMakeLists.txt b/CMakeLists.txt index 5a204fb66..6bc51e812 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -337,7 +337,6 @@ add_library(libopenmc SHARED src/stl_vector.F90 src/string.F90 src/surface_header.F90 - src/tracking.F90 src/track_output.F90 src/vector_header.F90 src/xml_interface.F90 diff --git a/include/openmc/cell.h b/include/openmc/cell.h index c90cb7f69..ba202c343 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -195,5 +195,13 @@ public: }; #endif +//============================================================================== +// Non-member functions +//============================================================================== + +#ifdef DAGMC +int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed); +#endif + } // namespace openmc #endif // OPENMC_CELL_H diff --git a/include/openmc/error.h b/include/openmc/error.h index ee7168a4e..206c9b5b2 100644 --- a/include/openmc/error.h +++ b/include/openmc/error.h @@ -49,7 +49,7 @@ void warning(const std::stringstream& message) warning(message.str()); } -void write_message(const std::string& message, int level); +void write_message(const std::string& message, int level=0); inline void write_message(const std::stringstream& message, int level) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 84e628fe2..a81b3e77a 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -148,6 +148,12 @@ extern "C" { //! \param src Source site data void from_source(const Bank* src); + //! Transport the particle + void transport(); + + //! Cross a surface + void cross_surface(); + //! mark a particle as lost and create a particle restart file //! \param message A warning message to display void mark_as_lost(const char* message); diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 5d480e0bf..0c2b399d1 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -14,7 +14,7 @@ namespace openmc { //============================================================================== //! Sample a nuclide and reaction and then calls the appropriate routine -extern "C" void collision(Particle* p); +void collision(Particle* p); //! Samples an incident neutron reaction void sample_neutron_reaction(Particle* p); diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index 8febad19a..511a0658b 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -12,7 +12,7 @@ namespace openmc { //! \brief samples particle behavior after a collision event. //! \param p Particle to operate on -extern "C" void +void collision_mg(Particle* p); //! \brief samples a reaction type. diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index 3237fb3ba..21df00624 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -35,6 +35,13 @@ void apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, double atom_density, int score_bin, double& score); +void score_collision_derivative(const Particle* p); + +void score_track_derivative(const Particle* p, double distance); + +//! Set the flux derivatives on differential tallies to zero. +void zero_flux_derivs(); + } // namespace openmc //============================================================================== diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 8aa659325..12016956c 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -50,6 +50,16 @@ private: // Non-member functions //============================================================================== +void score_collision_tally(const Particle* p); + +void score_analog_tally_ce(const Particle* p); + +void score_analog_tally_mg(const Particle* p); + +void score_tracklength_tally(const Particle* p, double distance); + +void score_surface_tally(const Particle* p, const std::vector& tallies); + } // namespace openmc #endif // OPENMC_TALLIES_TALLY_SCORING_H diff --git a/src/cell.cpp b/src/cell.cpp index 642059506..56635309f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -805,12 +805,12 @@ extern "C" { } #ifdef DAGMC - int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed ) + int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed) { moab::EntityHandle surf = - surf_xed->dagmc_ptr_->entity_by_id(2,surf_xed->id_); + surf_xed->dagmc_ptr_->entity_by_id(2, surf_xed->id_); moab::EntityHandle vol = - cur_cell->dagmc_ptr_->entity_by_id(3,cur_cell->id_); + cur_cell->dagmc_ptr_->entity_by_id(3, cur_cell->id_); moab::EntityHandle new_vol; cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol); diff --git a/src/geometry.F90 b/src/geometry.F90 index 7456644e6..02ababae6 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -50,34 +50,10 @@ module geometry integer(C_INT), intent(out) :: lattice_translation(3) integer(C_INT), intent(out) :: next_level end subroutine distance_to_boundary - -#ifdef DAGMC - - function next_cell_c(current_cell, surface_crossed) & - bind(C, name="next_cell") result(new_cell) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: current_cell - type(C_PTR), intent(in), value :: surface_crossed - integer(C_INT32_T) :: new_cell - end function next_cell_c - -#endif - end interface contains -#ifdef DAGMC - - function next_cell(c, s) result(new_cell) - type(Cell), intent(in) :: c - type(Surface), intent(in) :: s - integer :: new_cell - new_cell = next_cell_c(c%ptr, s%ptr) - end function next_cell - -#endif - !=============================================================================== ! FIND_CELL determines what cell a source particle is in within a particular ! universe. If the base universe is passed, the particle should be found as long diff --git a/src/material.cpp b/src/material.cpp index 510cea966..41c9d0a01 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1135,11 +1135,6 @@ extern "C" { return model::materials[i_mat - 1]->density_gpcc_; } - void material_calculate_xs(const Particle* p) - { - model::materials[p->material - 1]->calculate_xs(*p); - } - void free_memory_material() { for (Material *mat : model::materials) {delete mat;} diff --git a/src/material_header.F90 b/src/material_header.F90 index cda79ec60..dc1e68944 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -2,12 +2,9 @@ module material_header use, intrinsic :: ISO_C_BINDING - use particle_header, only: Particle - implicit none private - public :: material_calculate_xs public :: material_id public :: material_nuclide public :: material_nuclide_size @@ -21,11 +18,6 @@ module material_header integer(C_INT32_T) :: id end function - subroutine material_calculate_xs(p) bind(C) - import Particle - type(Particle), intent(in) :: p - end subroutine - function material_nuclide(i_mat, idx) bind(C) result(nuc) import C_INT32_T, C_INT integer(C_INT32_T), value :: i_mat diff --git a/src/particle.cpp b/src/particle.cpp index e5a38a983..26c4e8a7f 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -1,16 +1,29 @@ #include "openmc/particle.h" -#include +#include // copy, min +#include // log, abs, copysign #include #include "openmc/bank.h" #include "openmc/capi.h" +#include "openmc/cell.h" #include "openmc/constants.h" #include "openmc/error.h" +#include "openmc/geometry.h" #include "openmc/hdf5_interface.h" +#include "openmc/material.h" +#include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" +#include "openmc/nuclide.h" +#include "openmc/physics.h" +#include "openmc/physics_mg.h" +#include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/surface.h" #include "openmc/simulation.h" +#include "openmc/tallies/derivative.h" +#include "openmc/tallies/tally.h" +#include "openmc/tallies/tally_scoring.h" namespace openmc { @@ -119,6 +132,500 @@ Particle::from_source(const Bank* src) last_E = E; } +extern "C" void initialize_particle_track(); +extern "C" void write_particle_track(const Particle*); +extern "C" void add_particle_track(); +extern "C" void finalize_particle_track(const Particle*); + +void +Particle::transport() +{ + // Display message if high verbosity or trace is on + if (settings::verbosity >= 9 || simulation::trace) { + write_message("Simulating Particle " + std::to_string(id)); + } + + // Initialize number of events to zero + int n_event = 0; + + // Add paricle's starting weight to count for normalizing tallies later +#pragma omp atomic + total_weight += wgt; + + // Force calculation of cross-sections by setting last energy to zero + if (settings::run_CE) { + for (int i = 0; i < data::nuclides.size(); ++i) { + simulation::micro_xs[i].last_E = 0.0; + } + } + + // Prepare to write out particle track. + if (write_track) initialize_particle_track(); + + // Every particle starts with no accumulated flux derivative. + if (!model::active_tallies.empty()) zero_flux_derivs(); + + while (true) { + // Set the random number stream + if (type == static_cast(ParticleType::neutron)) { + prn_set_stream(STREAM_TRACKING); + } else { + prn_set_stream(STREAM_PHOTON); + } + + // Store pre-collision particle properties + last_wgt = wgt; + last_E = E; + std::copy(coord[0].uvw, coord[0].uvw + 3, last_uvw); + std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz); + + // If the cell hasn't been determined based on the particle's location, + // initiate a search for the current cell. This generally happens at the + // beginning of the history and again for any secondary particles + if (coord[n_coord - 1].cell == C_NONE) { + if (!find_cell(this, false)) { + this->mark_as_lost("Could not find the cell containing particle " + + std::to_string(id)); + return; + } + + // set birth cell attribute + if (cell_born == C_NONE) cell_born = coord[n_coord - 1].cell; + } + + // Write particle track. + if (write_track) write_particle_track(this); + + if (settings::check_overlaps) check_cell_overlap(this); + + // Calculate microscopic and macroscopic cross sections + if (material != MATERIAL_VOID) { + if (settings::run_CE) { + if (material != last_material || sqrtkT != last_sqrtkT) { + // If the material is the same as the last material and the + // temperature hasn't changed, we don't need to lookup cross + // sections again. + model::materials[material - 1]->calculate_xs(*this); + } + } else { + // Get the MG data + calculate_xs_c(material, g, sqrtkT, coord[n_coord-1].uvw, + simulation::material_xs.total, simulation::material_xs.absorption, + simulation::material_xs.nu_fission); + + // Finally, update the particle group while we have already checked + // for if multi-group + last_g = g; + } + } else { + simulation::material_xs.total = 0.0; + simulation::material_xs.absorption = 0.0; + simulation::material_xs.fission = 0.0; + simulation::material_xs.nu_fission = 0.0; + } + + // Find the distance to the nearest boundary + double d_boundary; + int surface_crossed; + int lattice_translation[3]; + int next_level; + distance_to_boundary(this, &d_boundary, &surface_crossed, + lattice_translation, &next_level); + + // Sample a distance to collision + double d_collision; + if (type == static_cast(ParticleType::electron) || + type == static_cast(ParticleType::positron)) { + d_collision = 0.0; + } else if (simulation::material_xs.total == 0.0) { + d_collision = INFINITY; + } else { + d_collision = -std::log(prn()) / simulation::material_xs.total; + } + + // Select smaller of the two distances + double distance = std::min(d_boundary, d_collision); + + // Advance particle + for (int j = 0; j < n_coord; ++j) { + // TODO: use Position + coord[j].xyz[0] += distance * coord[j].uvw[0]; + coord[j].xyz[1] += distance * coord[j].uvw[1]; + coord[j].xyz[2] += distance * coord[j].uvw[2]; + } + + // Score track-length tallies + if (!model::active_tracklength_tallies.empty()) { + score_tracklength_tally(this, distance); + } + + // Score track-length estimate of k-eff + if (settings::run_mode == RUN_MODE_EIGENVALUE && + type == static_cast(ParticleType::neutron)) { + global_tally_tracklength += wgt * distance * simulation::material_xs.nu_fission; + } + + // Score flux derivative accumulators for differential tallies. + if (!model::active_tallies.empty()) { + score_track_derivative(this, distance); + } + + if (d_collision > d_boundary) { + // ==================================================================== + // PARTICLE CROSSES SURFACE + + if (next_level > 0) n_coord = next_level; + + // Saving previous cell data + for (int j = 0; j < n_coord; ++j) { + last_cell[j] = coord[j].cell; + } + last_n_coord = n_coord; + + if (lattice_translation[0] != 0 || lattice_translation[1] != 0 || + lattice_translation[2] != 0) { + // Particle crosses lattice boundary + surface = ERROR_INT; + cross_lattice(this, lattice_translation); + event = EVENT_LATTICE; + } else { + // Particle crosses surface + surface = surface_crossed; + this->cross_surface(); + event = EVENT_SURFACE; + } + // Score cell to cell partial currents + if (!model::active_surface_tallies.empty()) { + score_surface_tally(this, model::active_surface_tallies); + } + } else { + // ==================================================================== + // PARTICLE HAS COLLISION + + // Score collision estimate of keff + if (settings::run_mode == RUN_MODE_EIGENVALUE && + type == static_cast(ParticleType::neutron)) { + global_tally_collision += wgt * simulation::material_xs.nu_fission + / simulation::material_xs.total; + } + + // Score surface current tallies -- this has to be done before the collision + // since the direction of the particle will change and we need to use the + // pre-collision direction to figure out what mesh surfaces were crossed + + if (!model::active_meshsurf_tallies.empty()) + score_surface_tally(this, model::active_meshsurf_tallies); + + // Clear surface component + surface = ERROR_INT; + + if (settings::run_CE) { + collision(this); + } else { + collision_mg(this); + } + + // Score collision estimator tallies -- this is done after a collision + // has occurred rather than before because we need information on the + // outgoing energy for any tallies with an outgoing energy filter + if (!model::active_collision_tallies.empty()) score_collision_tally(this); + if (!model::active_analog_tallies.empty()) { + if (settings::run_CE) { + score_analog_tally_ce(this); + } else { + score_analog_tally_mg(this); + } + } + + // Reset banked weight during collision + n_bank = 0; + wgt_bank = 0.0; + for (int& v : n_delayed_bank) v = 0; + + // Reset fission logical + fission = false; + + // Save coordinates for tallying purposes + std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz_current); + + // Set last material to none since cross sections will need to be + // re-evaluated + last_material = F90_NONE; + + // Set all uvws to base level -- right now, after a collision, only the + // base level uvws are changed + for (int j = 0; j < n_coord - 1; ++j) { + if (coord[j + 1].rotated) { + // If next level is rotated, apply rotation matrix + // FIXME + // coord[j + 1].uvw = matmul(cells(coord[j].cell + 1) % & + // rotation_matrix, coord[j].uvw) + } else { + // Otherwise, copy this level's direction + std::copy(coord[j].uvw, coord[j].uvw + 3, coord[j + 1].uvw); + } + } + + // Score flux derivative accumulators for differential tallies. + if (!model::active_tallies.empty()) score_collision_derivative(this); + } + + // If particle has too many events, display warning and kill it + ++n_event; + if (n_event == MAX_EVENTS) { + if (mpi::master) warning("Particle " + std::to_string(id) + + " underwent maximum number of events."); + alive = false; + } + + // Check for secondary particles if this particle is dead + if (!alive) { + // If no secondary particles, break out of event loop + if (n_secondary == 0) break; + + this->from_source(&secondary_bank[n_secondary - 1]); + --n_secondary; + n_event = 0; + + // Enter new particle in particle track file + if (write_track) add_particle_track(); + } + } + + // Finish particle track output. + if (write_track) { + write_particle_track(this); + finalize_particle_track(this); + } +} + +void +Particle::cross_surface() +{ + int i_surface = std::abs(surface); + // TODO: off-by-one + const auto& surf {model::surfaces[i_surface - 1]}; + if (settings::verbosity >= 10 || simulation::trace) { + write_message(" Crossing surface " + std::to_string(surf->id_)); + } + + if (surf->bc_ == BC_VACUUM && (settings::run_mode != RUN_MODE_PLOTTING)) { + // ======================================================================= + // PARTICLE LEAKS OUT OF PROBLEM + + // Kill particle + alive = false; + + // Score any surface current tallies -- note that the particle is moved + // forward slightly so that if the mesh boundary is on the surface, it is + // still processed + + if (!model::active_meshsurf_tallies.empty()) { + // TODO: Find a better solution to score surface currents than + // physically moving the particle forward slightly + + // TODO: Use Position + coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0]; + coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1]; + coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2]; + score_surface_tally(this, model::active_meshsurf_tallies); + } + + // Score to global leakage tally + global_tally_leakage += wgt; + + // Display message + if (settings::verbosity >= 10 || simulation::trace) { + write_message(" Leaked out of surface " + std::to_string(surf->id_)); + } + return; + + } else if (surf->bc_ == BC_REFLECT && (settings::run_mode != RUN_MODE_PLOTTING)) { + // ======================================================================= + // PARTICLE REFLECTS FROM SURFACE + + // Do not handle reflective boundary conditions on lower universes + if (n_coord != 1) { + this->mark_as_lost("Cannot reflect particle " + std::to_string(id) + + " off surface in a lower universe."); + return; + } + + // Score surface currents since reflection causes the direction of the + // particle to change. For surface filters, we need to score the tallies + // twice, once before the particle's surface attribute has changed and + // once after. For mesh surface filters, we need to artificially move + // the particle slightly back in case the surface crossing is coincident + // with a mesh boundary + + if (!model::active_surface_tallies.empty()) { + score_surface_tally(this, model::active_surface_tallies); + } + + + if (!model::active_meshsurf_tallies.empty()) { + Position r {coord[0].xyz}; + coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0]; + coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1]; + coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2]; + score_surface_tally(this, model::active_meshsurf_tallies); + std::copy(&r.x, &r.x + 3, coord[0].xyz); + } + + // Reflect particle off surface + Direction u = surf->reflect(coord[0].xyz, coord[0].uvw); + + // Make sure new particle direction is normalized + double norm = u.norm(); + coord[0].uvw[0] = u.x/norm; + coord[0].uvw[1] = u.y/norm; + coord[0].uvw[2] = u.z/norm; + + // Reassign particle's cell and surface + coord[0].cell = last_cell[last_n_coord - 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. + + 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 + last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT*coord[0].uvw[0]; + last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT*coord[0].uvw[1]; + last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT*coord[0].uvw[2]; + + // Diagnostic message + if (settings::verbosity >= 10 || simulation::trace) { + write_message(" Reflected from surface " + std::to_string(surf->id_)); + } + return; + + } else if (surf->bc_ == BC_PERIODIC && settings::run_mode != RUN_MODE_PLOTTING) { + // ======================================================================= + // PERIODIC BOUNDARY + + // Do not handle periodic boundary conditions on lower universes + if (n_coord != 1) { + this->mark_as_lost("Cannot transfer particle " + std::to_string(id) + + " across surface in a lower universe. Boundary conditions must be " + "applied to root universe."); + return; + } + + // Score surface currents since reflection causes the direction of the + // particle to change -- artificially move the particle slightly back in + // case the surface crossing is coincident with a mesh boundary + if (!model::active_meshsurf_tallies.empty()) { + Position r {coord[0].xyz}; + coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0]; + coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1]; + coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2]; + score_surface_tally(this, model::active_meshsurf_tallies); + std::copy(&r.x, &r.x + 3, coord[0].xyz); + } + + // Get a pointer to the partner periodic surface. Offset the index to + // correct for C vs. Fortran indexing. + auto surf_p = dynamic_cast(surf); + if (surf_p) { + auto other = dynamic_cast( + model::surfaces[surf_p->i_periodic_]); + + // Adjust the particle's location and direction. + Position r {coord[0].xyz}; + Direction u {coord[0].uvw}; + bool rotational = other->periodic_translate(surf_p, r, u); + std::copy(&r.x, &r.x + 3, coord[0].xyz); + std::copy(&u.x, &u.x + 3, coord[0].uvw); + + // Reassign particle's surface + // TODO: off-by-one + surface = rotational ? + surf_p->i_periodic_ + 1 : + std::copysign(surf_p->i_periodic_ + 1, surface); + } + + // Figure out what cell particle is in now + n_coord = 1; + + if (!find_cell(this, true)) { + this->mark_as_lost("Couldn't find particle after hitting periodic " + "boundary on surface " + std::to_string(surf->id_) + "."); + return; + } + + // Set previous coordinate going slightly past surface crossing + last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT * coord[0].uvw[0]; + last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT * coord[0].uvw[1]; + last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT * coord[0].uvw[2]; + + // Diagnostic message + if (settings::verbosity >= 10 || simulation::trace) { + write_message(" Hit periodic boundary on surface " + + std::to_string(surf->id_)); + } + return; + } + + // ========================================================================== + // SEARCH NEIGHBOR LISTS FOR NEXT CELL + +#ifdef DAGMC + if (settings::dagmc) { + int32_t i_cell = next_cell(model::cells[last_cell[0]], + model::surfaces[std::abs(surface)]); + // save material and temp + last_material = material; + last_sqrtkT = sqrtKT; + // set new cell value + coord[0].cell = i_cell + cell_instance = 0; + material = model::cells[i_cell]->material[0]; + sqrtKT = model::cells[i_cell]->sqrtKT[0]; + return + } +#endif + + if (find_cell(this, true)) return; + + // ========================================================================== + // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS + + // Remove lower coordinate levels and assignment of surface + surface = ERROR_INT; + n_coord = 1; + bool found = find_cell(this, false); + + if (settings::run_mode != RUN_MODE_PLOTTING && (!found)) { + // If a cell is still not found, there are two possible causes: 1) there is + // a void in the model, and 2) the particle hit a surface at a tangent. If + // the particle is really traveling tangent to a surface, if we move it + // forward a tiny bit it should fix the problem. + + n_coord = 1; + coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0]; + coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1]; + coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2]; + + // Couldn't find next cell anywhere// This probably means there is an actual + // undefined region in the geometry. + + if (!find_cell(this, false)) { + this->mark_as_lost("After particle " + std::to_string(id) + + " crossed surface " + std::to_string(surf->id_) + + " it could not be located in any cell and it did not leak."); + return; + } + } +} + void Particle::mark_as_lost(const char* message) { diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index dd1c89669..31ee5ff66 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -74,7 +74,8 @@ contains call set_particle_seed(particle_seed) ! Transport neutron - call transport(p) + ! FIXME + !call transport(p) ! Write output if particle made it call print_particle(p) diff --git a/src/simulation.cpp b/src/simulation.cpp index 7c7bbc4e6..acc488353 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -29,11 +29,9 @@ namespace openmc { // data/functions from Fortran side extern "C" void accumulate_tallies(); extern "C" void allocate_tally_results(); -extern "C" void init_tally_routines(); extern "C" void setup_active_tallies(); extern "C" void simulation_init_f(); extern "C" void simulation_finalize_f(); -extern "C" void transport(Particle* p); extern "C" void write_tallies(); } // namespace openmc @@ -76,9 +74,6 @@ int openmc_simulation_init() simulation::filter_matches.resize(model::tally_filters.size()); } - // Set up tally procedure pointers - init_tally_routines(); - // Allocate source bank, and for eigenvalue simulations also allocate the // fission bank allocate_banks(); @@ -213,7 +208,7 @@ int openmc_next_batch(int* status) initialize_history(&p, simulation::current_work); // transport particle - transport(&p); + p.transport(); } // Accumulate time for transport diff --git a/src/surface.cpp b/src/surface.cpp index 3ffbb3277..8b9b20ea1 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -340,7 +340,7 @@ bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const { Direction other_n = other->normal(r); - if (other_n.x == 1 and other_n.y == 0 and other_n.z == 0) { + if (other_n.x == 1 && other_n.y == 0 && other_n.z == 0) { r.x = x0_; return false; } else { @@ -397,11 +397,11 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, +bool SurfaceYPlane::periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const { Direction other_n = other->normal(r); - if (other_n.x == 0 and other_n.y == 1 and other_n.z == 0) { + 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; diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 600cd2876..c3bf8f1bf 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -571,7 +571,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, //! Adjust diff tally flux derivatives for a particle tracking event. -extern "C" void +void score_track_derivative(const Particle* p, double distance) { // A void material cannot be perturbed so it will not affect flux derivatives. @@ -627,8 +627,7 @@ score_track_derivative(const Particle* p, double distance) //! absorption will never be tallied. The paricle will be killed before any //! further tallies are scored. -extern "C" void -score_collision_derivative(const Particle* p) +void score_collision_derivative(const Particle* p) { // A void material cannot be perturbed so it will not affect flux derivatives. if (p->material == MATERIAL_VOID) return; @@ -698,9 +697,7 @@ score_collision_derivative(const Particle* p) } } -//! Set the flux derivatives on differential tallies to zero. -extern "C" void -zero_flux_derivs() +void zero_flux_derivs() { for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.; } diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 96b2fff06..003eb438d 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -11,87 +11,16 @@ module tally use message_passing use mgxs_interface use nuclide_header - use particle_header, only: LocalCoord, Particle use settings use simulation_header - use string, only: to_str use tally_derivative_header use tally_filter use tally_header implicit none - procedure(score_analog_tally_), pointer :: score_analog_tally => null() - - abstract interface - subroutine score_analog_tally_(p) - import Particle - type(Particle), intent(in) :: p - end subroutine score_analog_tally_ - end interface - - interface - subroutine score_analog_tally_ce(p) bind(C) - import Particle - type(Particle), intent(in) :: p - end subroutine - - subroutine score_analog_tally_mg(p) bind(C) - import Particle - type(Particle), intent(in) :: p - end subroutine - - subroutine score_tracklength_tally(p, distance) bind(C) - import Particle, C_DOUBLE - type(Particle) :: p - real(C_DOUBLE), value :: distance - end subroutine - - subroutine score_collision_tally(p) bind(C) - import Particle - type(Particle) :: p - end subroutine - - subroutine score_meshsurface_tally(p) bind(C) - import Particle - type(Particle) :: p - end subroutine - - subroutine score_surface_tally(p) bind(C) - import Particle - type(Particle) :: p - end subroutine - - subroutine score_track_derivative(p, distance) bind(C) - import Particle, C_DOUBLE - type(Particle) :: p - real(C_DOUBLE), value :: distance - end subroutine - - subroutine score_collision_derivative(p) bind(C) - import Particle - type(Particle) :: p - end subroutine - - subroutine zero_flux_derivs() bind(C) - end subroutine - end interface - contains -!=============================================================================== -! INIT_TALLY_ROUTINES Sets the procedure pointers needed for minimizing code -! with the CE and MG modes. -!=============================================================================== - - subroutine init_tally_routines() bind(C) - if (run_CE) then - score_analog_tally => score_analog_tally_ce - else - score_analog_tally => score_analog_tally_mg - end if - end subroutine init_tally_routines - !=============================================================================== ! ACCUMULATE_TALLIES accumulates the sum of the contributions from each history ! within the batch to a new random variable diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 8c509a0aa..fd96d9d85 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1976,8 +1976,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux, // //! Analog tallies ar etriggered at every collision, not every event. -extern "C" void -score_analog_tally_ce(const Particle* p) +void score_analog_tally_ce(const Particle* p) { for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2040,8 +2039,7 @@ score_analog_tally_ce(const Particle* p) // //! Analog tallies ar etriggered at every collision, not every event. -extern "C" void -score_analog_tally_mg(const Particle* p) +void score_analog_tally_mg(const Particle* p) { for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2096,7 +2094,7 @@ score_analog_tally_mg(const Particle* p) //! collision) and thus cannot be done for tallies that require post-collision //! information. -extern "C" void +void score_tracklength_tally(const Particle* p, double distance) { // Determine the tracklength estimate of the flux @@ -2172,8 +2170,7 @@ score_tracklength_tally(const Particle* p, double distance) //! reactions that we didn't sample. It is assumed the material is not void //! since collisions do not occur in voids. -extern "C" void -score_collision_tally(const Particle* p) +void score_collision_tally(const Particle* p) { // Determine the collision estimate of the flux double flux; @@ -2243,8 +2240,8 @@ score_collision_tally(const Particle* p) //! Score surface or mesh-surface tallies for particle currents. -static void -score_surface_tally_inner(const Particle* p, const std::vector& tallies) +void +score_surface_tally(const Particle* p, const std::vector& tallies) { // No collision, so no weight change when survival biasing double flux = p->wgt; @@ -2289,20 +2286,4 @@ score_surface_tally_inner(const Particle* p, const std::vector& tallies) match.bins_present_ = false; } -//! Score mesh-surface tallies for particle currents. - -extern "C" void -score_meshsurface_tally(const Particle* p) -{ - score_surface_tally_inner(p, model::active_meshsurf_tallies); -} - -//! Score surface tallies for particle currents. - -extern "C" void -score_surface_tally(const Particle* p) -{ - score_surface_tally_inner(p, model::active_surface_tallies); -} - } // namespace openmc diff --git a/src/track_output.F90 b/src/track_output.F90 index 95af9e6b0..cf1625ca3 100644 --- a/src/track_output.F90 +++ b/src/track_output.F90 @@ -5,6 +5,8 @@ module track_output + use, intrinsic :: ISO_C_BINDING + use constants use hdf5_interface use particle_header, only: Particle @@ -34,7 +36,7 @@ contains ! information !=============================================================================== - subroutine initialize_particle_track() + subroutine initialize_particle_track() bind(C) allocate(tracks(1)) end subroutine initialize_particle_track @@ -42,7 +44,7 @@ contains ! WRITE_PARTICLE_TRACK copies particle position to an array. !=============================================================================== - subroutine write_particle_track(p) + subroutine write_particle_track(p) bind(C) type(Particle), intent(in) :: p real(8), allocatable :: new_coords(:, :) @@ -71,7 +73,7 @@ contains ! secondary particle !=============================================================================== - subroutine add_particle_track() + subroutine add_particle_track() bind(C) type(TrackCoordinates), allocatable :: new_tracks(:) integer :: i @@ -92,7 +94,7 @@ contains ! FINALIZE_PARTICLE_TRACK writes the particle track array to disk. !=============================================================================== - subroutine finalize_particle_track(p) + subroutine finalize_particle_track(p) bind(C) type(Particle), intent(in) :: p integer :: i diff --git a/src/tracking.F90 b/src/tracking.F90 deleted file mode 100644 index d3c85bd79..000000000 --- a/src/tracking.F90 +++ /dev/null @@ -1,539 +0,0 @@ -module tracking - - use, intrinsic :: ISO_C_BINDING - - use constants - use error, only: warning, write_message - use geometry_header, only: cells - use geometry, only: find_cell, distance_to_boundary, cross_lattice,& - check_cell_overlap -#ifdef DAGMC - use geometry, only: next_cell -#endif - - use material_header, only: material_calculate_xs - use message_passing - use mgxs_interface - use nuclide_header - use particle_header - use random_lcg, only: prn, prn_set_stream - use settings - use simulation_header - use string, only: to_str - use surface_header - use tally_header - use tally, only: score_analog_tally, score_tracklength_tally, & - score_collision_tally, score_surface_tally, & - score_meshsurface_tally, & - score_track_derivative, zero_flux_derivs, & - score_collision_derivative - use track_output, only: initialize_particle_track, write_particle_track, & - add_particle_track, finalize_particle_track - - implicit none - - interface - subroutine collision(p) bind(C) - import Particle - type(Particle), intent(inout) :: p - end subroutine - - subroutine collision_mg(p) bind(C) - import Particle, C_DOUBLE - type(Particle), intent(inout) :: p - end subroutine collision_mg - end interface - -contains - -!=============================================================================== -! TRANSPORT encompasses the main logic for moving a particle through geometry. -!=============================================================================== - - subroutine transport(p) bind(C) - - type(Particle), intent(inout) :: p - - integer :: j ! coordinate level - integer :: next_level ! next coordinate level to check - integer :: surface_crossed ! surface which particle is on - integer :: lattice_translation(3) ! in-lattice translation vector - integer :: n_event ! number of collisions/crossings - real(8) :: d_boundary ! distance to nearest boundary - real(8) :: d_collision ! sampled distance to collision - real(8) :: distance ! distance particle travels - logical :: found_cell ! found cell which particle is in? - - ! Display message if high verbosity or trace is on - if (verbosity >= 9 .or. trace) then - call write_message("Simulating Particle " // trim(to_str(p % id))) - end if - - ! Initialize number of events to zero - n_event = 0 - - ! Add paricle's starting weight to count for normalizing tallies later -!$omp atomic - total_weight = total_weight + p % wgt - - ! Force calculation of cross-sections by setting last energy to zero - if (run_CE) then - micro_xs % last_E = ZERO - end if - - ! Prepare to write out particle track. - if (p % write_track) then - call initialize_particle_track() - endif - - ! Every particle starts with no accumulated flux derivative. - if (active_tallies_size() > 0) call zero_flux_derivs() - - EVENT_LOOP: do - ! Set the random number stream - if (p % type == NEUTRON) then - call prn_set_stream(STREAM_TRACKING) - else - call prn_set_stream(STREAM_PHOTON) - end if - - ! Store pre-collision particle properties - p % last_wgt = p % wgt - p % last_E = p % E - p % last_uvw = p % coord(1) % uvw - p % last_xyz = p % coord(1) % xyz - - ! If the cell hasn't been determined based on the particle's location, - ! initiate a search for the current cell. This generally happens at the - ! beginning of the history and again for any secondary particles - if (p % coord(p % n_coord) % cell == C_NONE) then - call find_cell(p, found_cell) - if (.not. found_cell) then - call particle_mark_as_lost(p, "Could not find the cell containing" & - // " particle " // trim(to_str(p %id))) - return - end if - - ! set birth cell attribute - if (p % cell_born == C_NONE) p % cell_born = p % coord(p % n_coord) % cell - end if - - ! Write particle track. - if (p % write_track) call write_particle_track(p) - - if (check_overlaps) call check_cell_overlap(p) - - ! Calculate microscopic and macroscopic cross sections - if (p % material /= MATERIAL_VOID) then - if (run_CE) then - if (p % material /= p % last_material .or. & - p % sqrtkT /= p % last_sqrtkT) then - ! If the material is the same as the last material and the - ! temperature hasn't changed, we don't need to lookup cross - ! sections again. - call material_calculate_xs(p) - end if - else - ! Get the MG data - call calculate_xs_c(p % material, p % g, p % sqrtkT, & - p % coord(p % n_coord) % uvw, material_xs % total, & - material_xs % absorption, material_xs % nu_fission) - - - ! Finally, update the particle group while we have already checked - ! for if multi-group - p % last_g = p % g - end if - else - material_xs % total = ZERO - material_xs % absorption = ZERO - material_xs % fission = ZERO - material_xs % nu_fission = ZERO - end if - - ! Find the distance to the nearest boundary - call distance_to_boundary(p, d_boundary, surface_crossed, & - lattice_translation, next_level) - - ! Sample a distance to collision - if (p % type == ELECTRON .or. p % type == POSITRON) then - d_collision = ZERO - else if (material_xs % total == ZERO) then - d_collision = INFINITY - else - d_collision = -log(prn()) / material_xs % total - end if - - ! Select smaller of the two distances - distance = min(d_boundary, d_collision) - - ! Advance particle - do j = 1, p % n_coord - p % coord(j) % xyz = p % coord(j) % xyz + distance * p % coord(j) % uvw - end do - - ! Score track-length tallies - if (active_tracklength_tallies_size() > 0) then - call score_tracklength_tally(p, distance) - end if - - ! Score track-length estimate of k-eff - if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then - global_tally_tracklength = global_tally_tracklength + p % wgt * & - distance * material_xs % nu_fission - end if - - ! Score flux derivative accumulators for differential tallies. - if (active_tallies_size() > 0) call score_track_derivative(p, distance) - - if (d_collision > d_boundary) then - ! ==================================================================== - ! PARTICLE CROSSES SURFACE - - if (next_level > 0) p % n_coord = next_level - - ! Saving previous cell data - do j = 1, p % n_coord - p % last_cell(j) = p % coord(j) % cell - end do - p % last_n_coord = p % n_coord - - if (any(lattice_translation /= 0)) then - ! Particle crosses lattice boundary - p % surface = ERROR_INT - call cross_lattice(p, lattice_translation) - p % event = EVENT_LATTICE - else - ! Particle crosses surface - p % surface = surface_crossed - - call cross_surface(p) - p % event = EVENT_SURFACE - end if - ! Score cell to cell partial currents - if (active_surface_tallies_size() > 0) call score_surface_tally(p) - else - ! ==================================================================== - ! PARTICLE HAS COLLISION - - ! Score collision estimate of keff - if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then - global_tally_collision = global_tally_collision + p % wgt * & - material_xs % nu_fission / material_xs % total - end if - - ! Score surface current tallies -- this has to be done before the collision - ! since the direction of the particle will change and we need to use the - ! pre-collision direction to figure out what mesh surfaces were crossed - - if (active_meshsurf_tallies_size() > 0) call score_meshsurface_tally(p) - - ! Clear surface component - p % surface = ERROR_INT - - if (run_CE) then - call collision(p) - else - call collision_mg(p) - end if - - ! Score collision estimator tallies -- this is done after a collision - ! has occurred rather than before because we need information on the - ! outgoing energy for any tallies with an outgoing energy filter - if (active_collision_tallies_size() > 0) call score_collision_tally(p) - if (active_analog_tallies_size() > 0) call score_analog_tally(p) - - ! Reset banked weight during collision - p % n_bank = 0 - p % wgt_bank = ZERO - p % n_delayed_bank(:) = 0 - - ! Reset fission logical - p % fission = .false. - - ! Save coordinates for tallying purposes - p % last_xyz_current = p % coord(1) % xyz - - ! Set last material to none since cross sections will need to be - ! re-evaluated - p % last_material = NONE - - ! Set all uvws to base level -- right now, after a collision, only the - ! base level uvws are changed - do j = 1, p % n_coord - 1 - if (p % coord(j + 1) % rotated) then - ! If next level is rotated, apply rotation matrix - p % coord(j + 1) % uvw = matmul(cells(p % coord(j) % cell + 1) % & - rotation_matrix, p % coord(j) % uvw) - else - ! Otherwise, copy this level's direction - p % coord(j + 1) % uvw = p % coord(j) % uvw - end if - end do - - ! Score flux derivative accumulators for differential tallies. - if (active_tallies_size() > 0) call score_collision_derivative(p) - end if - - ! If particle has too many events, display warning and kill it - n_event = n_event + 1 - if (n_event == MAX_EVENTS) then - if (master) call warning("Particle " // trim(to_str(p%id)) & - &// " underwent maximum number of events.") - p % alive = .false. - end if - - ! Check for secondary particles if this particle is dead - if (.not. p % alive) then - if (p % n_secondary > 0) then - call particle_from_source(p, p % secondary_bank(p % n_secondary)) - p % n_secondary = p % n_secondary - 1 - n_event = 0 - - ! Enter new particle in particle track file - if (p % write_track) call add_particle_track() - else - exit EVENT_LOOP - end if - end if - end do EVENT_LOOP - - ! Finish particle track output. - if (p % write_track) then - call write_particle_track(p) - call finalize_particle_track(p) - endif - - end subroutine transport - -!=============================================================================== -! CROSS_SURFACE handles all surface crossings, whether the particle leaks out of -! the geometry, is reflected, or crosses into a new lattice or cell -!=============================================================================== - - subroutine cross_surface(p) - type(Particle), intent(inout) :: p - - real(8) :: u ! x-component of direction - real(8) :: v ! y-component of direction - real(8) :: w ! z-component of direction - real(8) :: norm ! "norm" of surface normal - real(8) :: xyz(3) ! Saved global coordinate - integer :: i_surface ! index in surfaces -#ifdef DAGMC - integer :: i_cell ! index of new cell -#endif - logical :: rotational ! if rotational periodic BC applied - logical :: found ! particle found in universe? - class(Surface), pointer :: surf - class(Surface), pointer :: surf2 ! periodic partner surface - - i_surface = abs(p % surface) - surf => surfaces(i_surface) - if (verbosity >= 10 .or. trace) then - call write_message(" Crossing surface " // trim(to_str(surf % id()))) - end if - - if (surf % bc() == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then - ! ======================================================================= - ! PARTICLE LEAKS OUT OF PROBLEM - - ! Kill particle - p % 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 (active_meshsurf_tallies_size() > 0) then - ! TODO: Find a better solution to score surface currents than - ! physically moving the particle forward slightly - - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - call score_meshsurface_tally(p) - end if - - ! Score to global leakage tally - global_tally_leakage = global_tally_leakage + p % wgt - - ! Display message - if (verbosity >= 10 .or. trace) then - call write_message(" Leaked out of surface " & - &// trim(to_str(surf % id()))) - end if - return - - elseif (surf % bc() == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then - ! ======================================================================= - ! PARTICLE REFLECTS FROM SURFACE - - ! Do not handle reflective boundary conditions on lower universes - if (p % n_coord /= 1) then - call particle_mark_as_lost(p, "Cannot reflect particle " & - // trim(to_str(p % id)) // " off surface in a lower universe.") - return - end if - - ! 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 (active_surface_tallies_size() > 0) call score_surface_tally(p) - - - if (active_meshsurf_tallies_size() > 0) then - xyz = p % coord(1) % xyz - p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_meshsurface_tally(p) - p % coord(1) % xyz = xyz - end if - - ! Reflect particle off surface - call surf % reflect(p%coord(1)%xyz, p%coord(1)%uvw) - - ! Make sure new particle direction is normalized - u = p%coord(1)%uvw(1) - v = p%coord(1)%uvw(2) - w = p%coord(1)%uvw(3) - norm = sqrt(u*u + v*v + w*w) - p%coord(1)%uvw(:) = [u, v, w] / norm - - ! Reassign particle's cell and surface - p % coord(1) % cell = p % last_cell(p % last_n_coord) - p % surface = -p % 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. - - p % n_coord = 1 - call find_cell(p, found, .true.) - if (.not. found) then - call particle_mark_as_lost(p, "Couldn't find particle after reflecting& - & from surface " // trim(to_str(surf % id())) // ".") - return - end if - - ! Set previous coordinate going slightly past surface crossing - p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - - ! Diagnostic message - if (verbosity >= 10 .or. trace) then - call write_message(" Reflected from surface " & - &// trim(to_str(surf%id()))) - end if - return - elseif (surf % bc() == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then - ! ======================================================================= - ! PERIODIC BOUNDARY - - ! Do not handle periodic boundary conditions on lower universes - if (p % n_coord /= 1) then - call particle_mark_as_lost(p, "Cannot transfer particle " & - // trim(to_str(p % id)) // " across surface in a lower universe.& - & Boundary conditions must be applied to universe 0.") - return - end if - - ! 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 (active_meshsurf_tallies_size() > 0) then - xyz = p % coord(1) % xyz - p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_meshsurface_tally(p) - p % coord(1) % xyz = xyz - end if - - ! Get a pointer to the partner periodic surface. Offset the index to - ! correct for C vs. Fortran indexing. - surf2 => surfaces(surf % i_periodic() + 1) - - ! Adjust the particle's location and direction. - rotational = surf2 % periodic_translate(surf, p % coord(1) % xyz, & - p % coord(1) % uvw) - - ! Reassign particle's surface - if (rotational) then - p % surface = surf % i_periodic() + 1 - else - p % surface = sign(surf % i_periodic() + 1, p % surface) - end if - - ! Figure out what cell particle is in now - p % n_coord = 1 - call find_cell(p, found, .true.) - if (.not. found) then - call particle_mark_as_lost(p, "Couldn't find particle after hitting & - &periodic boundary on surface " // trim(to_str(surf % id())) & - // ".") - return - end if - - ! Set previous coordinate going slightly past surface crossing - p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - - ! Diagnostic message - if (verbosity >= 10 .or. trace) then - call write_message(" Hit periodic boundary on surface " & - // trim(to_str(surf%id()))) - end if - return - end if - - ! ========================================================================== - ! SEARCH NEIGHBOR LISTS FOR NEXT CELL - -#ifdef DAGMC - if (dagmc) then - i_cell = next_cell(cells(p % last_cell(1) + 1), surfaces(abs(p % surface))) - ! save material and temp - p % last_material = p % material - p % last_sqrtkT = p % sqrtKT - ! set new cell value - p % coord(1) % cell = i_cell-1 ! decrement for C++ indexing - p % cell_instance = 1 - p % material = cells(i_cell) % material(1) - p % sqrtKT = cells(i_cell) % sqrtKT(0) - return - end if -#endif - - call find_cell(p, found, .true.) - if (found) return - - ! ========================================================================== - ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS - - ! Remove lower coordinate levels and assignment of surface - p % surface = ERROR_INT - p % n_coord = 1 - call find_cell(p, found) - - if (run_mode /= MODE_PLOTTING .and. (.not. found)) then - ! If a cell is still not found, there are two possible causes: 1) there is - ! a void in the model, and 2) the particle hit a surface at a tangent. If - ! the particle is really traveling tangent to a surface, if we move it - ! forward a tiny bit it should fix the problem. - - p % n_coord = 1 - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - call find_cell(p, found) - - ! Couldn't find next cell anywhere! This probably means there is an actual - ! undefined region in the geometry. - - if (.not. found) then - call particle_mark_as_lost(p, "After particle " // trim(to_str(p % id)) & - // " crossed surface " // trim(to_str(surf % id())) & - // " it could not be located in any cell and it did not leak.") - return - end if - end if - - end subroutine cross_surface - -end module tracking