From 7325897ee9f4a0b80d215ae47ef51fc4005ea623 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Feb 2019 21:48:20 -0600 Subject: [PATCH] Address @smharper comments on #1169 --- include/openmc/hdf5_interface.h | 4 +++ include/openmc/particle.h | 4 +-- include/openmc/simulation.h | 1 + include/openmc/tallies/derivative.h | 12 ++++++++ include/openmc/tallies/tally.h | 2 -- include/openmc/tallies/tally_scoring.h | 31 +++++++++++++++++++ src/cell.cpp | 1 - src/geometry_aux.cpp | 3 -- src/particle.cpp | 41 ++++++++++++-------------- src/simulation.cpp | 2 +- src/tallies/derivative.cpp | 9 ------ src/tallies/tally.cpp | 4 +-- src/tallies/tally_scoring.cpp | 24 --------------- 13 files changed, 72 insertions(+), 66 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index f55b14171..dcf7ab8dd 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -223,6 +223,10 @@ read_attribute(hid_t obj_id, const char* name, std::vector& vec) // Templates/overloads for read_dataset and related methods //============================================================================== +// Template for scalars. We need to be careful that the compiler does not use +// this version of read_dataset for vectors, arrays, or other non-scalar types. +// enable_if_t allows us to conditionally remove the function from overload +// resolution when the type T doesn't meet a certain criterion. template inline std::enable_if_t>::value> read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index a81b3e77a..01c7eaefc 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -148,10 +148,10 @@ extern "C" { //! \param src Source site data void from_source(const Bank* src); - //! Transport the particle + //! Transport a particle from birth to death void transport(); - //! Cross a surface + //! Cross a surface and handle boundary conditions void cross_surface(); //! mark a particle as lost and create a particle restart file diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index af333f714..c7053c608 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -36,6 +36,7 @@ extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? extern "C" int restart_batch; //!< batch at which a restart job resumed extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? extern "C" int total_gen; //!< total number of generations simulated +extern "C" double total_weight; //!< Total source weight in a batch extern "C" int64_t work; //!< number of particles per process extern std::vector k_generation; diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index 21df00624..02532afc0 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -35,8 +35,20 @@ void apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, double atom_density, int score_bin, double& score); +//! Adjust diff tally flux derivatives for a particle scattering event. +// +//! Note that this subroutine will be called after absorption events in +//! addition to scattering events, but any flux derivatives scored after an +//! absorption will never be tallied. The paricle will be killed before any +//! further tallies are scored. +// +//! \param p The particle being tracked void score_collision_derivative(const Particle* p); +//! Adjust diff tally flux derivatives for a particle tracking event. +// +//! \param p The particle being tracked +//! \param distance The distance in [cm] traveled by the particle void score_track_derivative(const Particle* p, double distance); //! Set the flux derivatives on differential tallies to zero. diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 464ec9f3d..7767110c2 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -101,8 +101,6 @@ private: // Global variable declarations //============================================================================== -extern "C" double total_weight; - namespace model { extern std::vector> tallies; diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 12016956c..af03b1aa2 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -50,14 +50,45 @@ private: // Non-member functions //============================================================================== +//! Score tallies using a 1 / Sigma_t estimate of the flux. +// +//! This is triggered after every collision. It is invalid for tallies that +//! require post-collison information because it can score reactions that didn't +//! actually occur, and we don't a priori know what the outcome will be for +//! reactions that we didn't sample. It is assumed the material is not void +//! since collisions do not occur in voids. +// +//! \param p The particle being tracked void score_collision_tally(const Particle* p); +//! Score tallies based on a simple count of events (for continuous energy). +// +//! Analog tallies are triggered at every collision, not every event. +// +//! \param p The particle being tracked void score_analog_tally_ce(const Particle* p); +//! Score tallies based on a simple count of events (for multigroup). +// +//! Analog tallies are triggered at every collision, not every event. +// +//! \param p The particle being tracked void score_analog_tally_mg(const Particle* p); +//! Score tallies using a tracklength estimate of the flux. +// +//! This is triggered at every event (surface crossing, lattice crossing, or +//! collision) and thus cannot be done for tallies that require post-collision +//! information. +// +//! \param p The particle being tracked +//! \param distance The distance in [cm] traveled by the particle void score_tracklength_tally(const Particle* p, double distance); +//! Score surface or mesh-surface tallies for particle currents. +// +//! \param p The particle being tracked +//! \param tallies A vector of tallies to score to void score_surface_tally(const Particle* p, const std::vector& tallies); } // namespace openmc diff --git a/src/cell.cpp b/src/cell.cpp index af77094c9..6a962eced 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -753,7 +753,6 @@ extern "C" int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance) { if (index >= 0 && index < model::cells.size()) { - //TODO: off-by-one Cell& c {*model::cells[index]}; if (instance) { diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index f69a63d9b..d38b5f8f5 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -56,9 +56,6 @@ void read_geometry_xml() read_cells(root); read_lattices(root); - // ========================================================================== - // SETUP UNIVERSES - // Allocate universes, universe cell arrays, and assign base universe model::root_universe = find_root_universe(); } diff --git a/src/particle.cpp b/src/particle.cpp index 4c39087ce..db2421d78 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -149,8 +149,8 @@ Particle::transport() int n_event = 0; // Add paricle's starting weight to count for normalizing tallies later -#pragma omp atomic - total_weight += wgt; + #pragma omp atomic + simulation::total_weight += wgt; // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { @@ -375,8 +375,8 @@ Particle::transport() // 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."); + warning("Particle " + std::to_string(id) + + " underwent maximum number of events."); alive = false; } @@ -533,26 +533,23 @@ Particle::cross_surface() 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. + // Get a pointer to the partner periodic surface auto surf_p = dynamic_cast(surf); - if (surf_p) { - auto other = dynamic_cast( - model::surfaces[surf_p->i_periodic_]); + 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); + // 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); - } + // 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; @@ -620,7 +617,7 @@ Particle::cross_surface() 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 + // Couldn't find next cell anywhere! This probably means there is an actual // undefined region in the geometry. if (!find_cell(this, false)) { diff --git a/src/simulation.cpp b/src/simulation.cpp index 1fdda3061..096ea87ff 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -312,7 +312,7 @@ void initialize_batch() } // Reset total starting particle weight used for normalizing tallies - total_weight = 0.0; + simulation::total_weight = 0.0; // Determine if this batch is the first inactive or active batch. bool first_inactive = false; diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index c3bf8f1bf..b28b78d0a 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -569,8 +569,6 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, } } -//! Adjust diff tally flux derivatives for a particle tracking event. - void score_track_derivative(const Particle* p, double distance) { @@ -620,13 +618,6 @@ score_track_derivative(const Particle* p, double distance) } } -//! Adjust diff tally flux derivatives for a particle scattering event. -// -//! Note that this subroutine will be called after absorption events in -//! addition to scattering events, but any flux derivatives scored after an -//! absorption will never be tallied. The paricle will be killed before any -//! further tallies are scored. - void score_collision_derivative(const Particle* p) { // A void material cannot be perturbed so it will not affect flux derivatives. diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 4bfb15eff..93ba72a19 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -585,9 +585,9 @@ void reduce_tally_results() // We also need to determine the total starting weight of particles from the // last realization double weight_reduced; - MPI_Reduce(&total_weight, &weight_reduced, 1, MPI_DOUBLE, MPI_SUM, + MPI_Reduce(&simulation::total_weight, &weight_reduced, 1, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); - if (mpi::master) total_weight = weight_reduced; + if (mpi::master) simulation::total_weight = weight_reduced; } #endif diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index fd96d9d85..690f61f81 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1972,10 +1972,6 @@ score_all_nuclides(const Particle* p, int i_tally, double flux, } } -//! Score tallies based on a simple count of events (for continuous energy). -// -//! Analog tallies ar etriggered at every collision, not every event. - void score_analog_tally_ce(const Particle* p) { for (auto i_tally : model::active_analog_tallies) { @@ -2035,10 +2031,6 @@ void score_analog_tally_ce(const Particle* p) match.bins_present_ = false; } -//! Score tallies based on a simple count of events (for multigroup). -// -//! Analog tallies ar etriggered at every collision, not every event. - void score_analog_tally_mg(const Particle* p) { for (auto i_tally : model::active_analog_tallies) { @@ -2088,12 +2080,6 @@ void score_analog_tally_mg(const Particle* p) match.bins_present_ = false; } -//! Score tallies using a tracklength estimate of the flux. -// -//! This is triggered at every event (surface crossing, lattice crossing, or -//! collision) and thus cannot be done for tallies that require post-collision -//! information. - void score_tracklength_tally(const Particle* p, double distance) { @@ -2162,14 +2148,6 @@ score_tracklength_tally(const Particle* p, double distance) match.bins_present_ = false; } -//! Score tallies using a 1 / Sigma_t estimate of the flux. -// -//! This is triggered after every collision. It is invalid for tallies that -//! require post-collison information because it can score reactions that didn't -//! actually occur, and we don't a priori know what the outcome will be for -//! reactions that we didn't sample. It is assumed the material is not void -//! since collisions do not occur in voids. - void score_collision_tally(const Particle* p) { // Determine the collision estimate of the flux @@ -2238,8 +2216,6 @@ void score_collision_tally(const Particle* p) match.bins_present_ = false; } -//! Score surface or mesh-surface tallies for particle currents. - void score_surface_tally(const Particle* p, const std::vector& tallies) {