Merge pull request #1113 from paulromano/timer-bugfix

Assortment of small bug fixes

CMakeLists.txt

@@ -131,7 +131,7 @@ elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel)
   # Intel compiler options
   list(APPEND f90flags -fpp -std08 -assume byterecl -traceback)
   if(debug)
-    list(APPEND f90flags -g -warn -ftrapuv -fp-stack-check
+    list(APPEND f90flags -g -warn -fp-stack-check
       "-check all" -fpe0 -O0)
     list(APPEND ldflags -g)
   endif()

include/openmc/tallies/filter.h

@@ -13,19 +13,6 @@
 
 namespace openmc {
 
-//==============================================================================
-// Global variables
-//==============================================================================
-
-extern "C" int32_t n_filters;
-
-class FilterMatch;
-extern std::vector<FilterMatch> filter_matches;
-#pragma omp threadprivate(filter_matches)
-
-class Filter;
-extern std::vector<std::unique_ptr<Filter>> tally_filters;
-
 //==============================================================================
 //! Stores bins and weights for filtered tally events.
 //==============================================================================
@@ -37,6 +24,15 @@ public:
   std::vector<double> weights_;
 };
 
+} // namespace openmc
+
+// Without an explicit instantiation of vector<FilterMatch>, the Intel compiler
+// will complain about the threadprivate directive on filter_matches. Note that
+// this has to happen *outside* of the openmc namespace
+template class std::vector<openmc::FilterMatch>;
+
+namespace openmc {
+
 //==============================================================================
 //! Modifies tally score events.
 //==============================================================================
@@ -77,6 +73,17 @@ public:
   int n_bins_;
 };
 
+//==============================================================================
+// Global variables
+//==============================================================================
+
+extern "C" int32_t n_filters;
+
+extern std::vector<FilterMatch> filter_matches;
+#pragma omp threadprivate(filter_matches)
+
+extern std::vector<std::unique_ptr<Filter>> tally_filters;
+
 //==============================================================================
 
 extern "C" void free_memory_tally_c();

include/openmc/timer.h

@@ -31,7 +31,7 @@ public:
 private:
   bool running_ {false}; //!< is timer running?
   std::chrono::time_point<clock> start_; //!< starting point for clock
-  double elapsed_; //!< elasped time in [s]
+  double elapsed_ {0.0}; //!< elasped time in [s]
 };
 
 //==============================================================================

openmc/settings.py

@@ -227,7 +227,7 @@ class Settings(object):
         self._log_grid_bins = None
 
         self._dagmc = None
-        
+
     @property
     def run_mode(self):
         return self._run_mode
@@ -375,7 +375,7 @@ class Settings(object):
     @property
     def dagmc(self):
         return self._dagmc
-    
+
     @run_mode.setter
     def run_mode(self, run_mode):
         cv.check_value('run mode', run_mode, _RUN_MODES)
@@ -523,7 +523,7 @@ class Settings(object):
     def dagmc(self, dagmc):
         cv.check_type('dagmc geometry', dagmc, bool)
         self._dagmc = dagmc
-        
+
     @ptables.setter
     def ptables(self, ptables):
         cv.check_type('probability tables', ptables, bool)
@@ -961,8 +961,8 @@ class Settings(object):
     def _create_dagmc_subelement(self, root):
         if self._dagmc is not None:
             elem = ET.SubElement(root, "dagmc")
-            element.text = str(self._dagmc).lower()
-        
+            elem.text = str(self._dagmc).lower()
+
     def export_to_xml(self, path='settings.xml'):
         """Export simulation settings to an XML file.
 
@@ -1011,7 +1011,7 @@ class Settings(object):
         self._create_create_fission_neutrons_subelement(root_element)
         self._create_log_grid_bins_subelement(root_element)
         self._create_dagmc_subelement(root_element)
-        
+
         # Clean the indentation in the file to be user-readable
         clean_indentation(root_element)
 

src/mesh.cpp

@@ -466,42 +466,38 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
     // ========================================================================
     // Compute the length of the track segment in the each mesh cell and return
 
-    double distance;
-    int j;
     if (ijk0 == ijk1) {
       // The track ends in this cell.  Use the particle end location rather
       // than the mesh surface.
-      distance = (r1 - r0).norm();
-    } else {
-      // The track exits this cell.  Determine the distance to the closest mesh
-      // surface.
-      xt::xtensor<double, 1> d = xt::zeros<double>({n});
-      for (int j = 0; j < n; ++j) {
-        if (std::fabs(u[j]) < FP_PRECISION) {
-          d(j) = INFTY;
-        } else if (u[j] > 0) {
-          double xyz_cross = lower_left_[j] + ijk0(j) * width_[j];
-          d(j) = (xyz_cross - r0[j]) / u[j];
-        } else {
-          double xyz_cross = lower_left_[j] + (ijk0(j) - 1) * width_[j];
-          d(j) = (xyz_cross - r0[j]) / u[j];
-        }
+      double distance = (r1 - r0).norm();
+      bins.push_back(get_bin_from_indices(ijk0.data()));
+      lengths.push_back(distance / total_distance);
+      break;
+    }
+
+    // The track exits this cell.  Determine the distance to the closest mesh
+    // surface.
+    xt::xtensor<double, 1> d = xt::zeros<double>({n});
+    for (int k = 0; k < n; ++k) {
+      if (std::fabs(u[k]) < FP_PRECISION) {
+        d(k) = INFTY;
+      } else if (u[k] > 0) {
+        double xyz_cross = lower_left_[k] + ijk0(k) * width_[k];
+        d(k) = (xyz_cross - r0[k]) / u[k];
+      } else {
+        double xyz_cross = lower_left_[k] + (ijk0(k) - 1) * width_[k];
+        d(k) = (xyz_cross - r0[k]) / u[k];
       }
-      j = xt::argmin(d)(0);
-      distance = d(j);
     }
 
     // Assign the next tally bin and the score.
-    int bin = get_bin_from_indices(ijk0.data());
-    bins.push_back(bin);
+    auto j = xt::argmin(d)(0);
+    double distance = d(j);
+    bins.push_back(get_bin_from_indices(ijk0.data()));
     lengths.push_back(distance / total_distance);
 
-    // If the particle track ends in that bin, then we are done.
-    if (ijk0 == ijk1) break;
-
-    // Translate the starting coordintes by the distance to that face. This
-    // should be the xyz that we computed the distance to in the last
-    // iteration of the filter loop.
+    // Translate the starting coordintes by the distance to the oncoming mesh
+    // surface.
     r0 += distance * u;
 
     // Increment the indices into the next mesh cell.

src/mesh_header.F90

@@ -140,21 +140,21 @@ module mesh_header
 
     pure function mesh_get_bin(ptr, xyz) result(bin) bind(C)
       import C_PTR, C_DOUBLE, C_INT
-      type(C_PTR), value :: ptr
+      type(C_PTR), intent(in), value :: ptr
       real(C_DOUBLE), intent(in) :: xyz(*)
       integer(C_INT) :: bin
     end function
 
     pure function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C)
       import C_PTR, C_INT
-      type(C_PTR), value :: ptr
+      type(C_PTR), intent(in), value :: ptr
       integer(C_INT), intent(in) :: ijk(*)
       integer(C_INT) :: bin
     end function
 
     pure subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C)
       import C_PTR, C_DOUBLE, C_INT, C_BOOL
-      type(C_PTR), value :: ptr
+      type(C_PTR), intent(in), value :: ptr
       real(C_DOUBLE), intent(in) :: xyz(*)
       integer(C_INT), intent(out) :: ijk(*)
       logical(C_BOOL), intent(out) :: in_mesh
@@ -162,8 +162,8 @@ module mesh_header
 
     pure subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C)
       import C_PTR, C_INT
-      type(C_PTR), value :: ptr
-      integer(C_INT), value :: bin
+      type(C_PTR), intent(in), value :: ptr
+      integer(C_INT), intent(in), value :: bin
       integer(C_INT), intent(out) :: ijk(*)
     end subroutine
 

src/reaction_header.F90

@@ -74,22 +74,22 @@ module reaction_header
 
     pure function reaction_product_decay_rate(ptr, product) result(rate) bind(C)
       import C_PTR, C_INT, C_DOUBLE
-      type(C_PTR), value :: ptr
-      integer(C_INT), value :: product
+      type(C_PTR), intent(in), value :: ptr
+      integer(C_INT), intent(in), value :: product
       real(C_DOUBLE) :: rate
     end function
 
     pure function reaction_product_emission_mode(ptr, product) result(m) bind(C)
       import C_PTR, C_INT
-      type(C_PTR), value :: ptr
-      integer(C_INT), value :: product
+      type(C_PTR), intent(in), value :: ptr
+      integer(C_INT), intent(in), value :: product
       integer(C_INT) :: m
     end function
 
     pure function reaction_product_particle(ptr, product) result(particle) bind(C)
       import C_PTR, C_INT
-      type(C_PTR), value :: ptr
-      integer(C_INT), value :: product
+      type(C_PTR), intent(in), value :: ptr
+      integer(C_INT), intent(in), value :: product
       integer(C_INT) :: particle
     end function
 
@@ -104,15 +104,15 @@ module reaction_header
 
     pure function reaction_product_yield(ptr, product, E) result(val) bind(C)
       import C_PTR, C_INT, C_DOUBLE
-      type(C_PTR), value :: ptr
-      integer(C_INT), value :: product
-      real(C_DOUBLE), value :: E
+      type(C_PTR), intent(in), value :: ptr
+      integer(C_INT), intent(in), value :: product
+      real(C_DOUBLE), intent(in), value :: E
       real(C_DOUBLE) :: val
     end function
 
     pure function reaction_products_size(ptr) result(sz) bind(C)
       import C_PTR, C_INT
-      type(C_PTR), value :: ptr
+      type(C_PTR), intent(in), value :: ptr
       integer(C_INT) :: sz
     end function
 

src/settings.cpp

@@ -205,7 +205,7 @@ void read_settings_xml()
 
   // Parse settings.xml file
   xml_document doc;
-  auto result = doc.load_file("settings.xml");
+  auto result = doc.load_file(filename.c_str());
   if (!result) {
     fatal_error("Error processing settings.xml file.");
   }

src/tallies/filter_energy.cpp

@@ -180,6 +180,7 @@ openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies)
   for (int i = 0; i < n; i++) filt->bins_[i] = energies[i];
   filt->n_bins_ = n - 1;
   filter_update_n_bins(index);
+  return 0;
 }
 
 }// namespace openmc

src/timer.cpp

@@ -28,7 +28,7 @@ double Timer::elapsed()
 {
   if (running_) {
     std::chrono::duration<double> diff = clock::now() - start_;
-    return elapsed_ += diff.count();
+    return elapsed_ + diff.count();
   } else {
     return elapsed_;
   }