Address @smharper comments on #1169

include/openmc/hdf5_interface.h

@@ -223,6 +223,10 @@ read_attribute(hid_t obj_id, const char* name, std::vector<std::string>& 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<typename T> inline
 std::enable_if_t<std::is_scalar<std::decay_t<T>>::value>
 read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)

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

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<double> k_generation;

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.

include/openmc/tallies/tally.h

@@ -101,8 +101,6 @@ private:
 // Global variable declarations
 //==============================================================================
 
-extern "C" double total_weight;
-
 namespace model {
 
 extern std::vector<std::unique_ptr<Tally>> tallies;

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<int>& tallies);
 
 } // namespace openmc

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) {

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();
 }

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<PeriodicSurface*>(surf);
-    if (surf_p) {
-      auto other = dynamic_cast<PeriodicSurface*>(
-        model::surfaces[surf_p->i_periodic_]);
+    auto other = dynamic_cast<PeriodicSurface*>(
+      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)) {

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;

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.

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
 

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<int>& tallies)
 {