diff --git a/CMakeLists.txt b/CMakeLists.txt index 4a9843ff2..fded5876c 100644 --- a/CMakeLists.txt +++ b/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() diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index af84b43eb..24ce5cb89 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -13,19 +13,6 @@ namespace openmc { -//============================================================================== -// Global variables -//============================================================================== - -extern "C" int32_t n_filters; - -class FilterMatch; -extern std::vector filter_matches; -#pragma omp threadprivate(filter_matches) - -class Filter; -extern std::vector> tally_filters; - //============================================================================== //! Stores bins and weights for filtered tally events. //============================================================================== @@ -37,6 +24,15 @@ public: std::vector weights_; }; +} // namespace openmc + +// Without an explicit instantiation of vector, 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; + +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 filter_matches; +#pragma omp threadprivate(filter_matches) + +extern std::vector> tally_filters; + //============================================================================== extern "C" void free_memory_tally_c(); diff --git a/include/openmc/timer.h b/include/openmc/timer.h index 1b4d71a5c..48b3a9e75 100644 --- a/include/openmc/timer.h +++ b/include/openmc/timer.h @@ -31,7 +31,7 @@ public: private: bool running_ {false}; //!< is timer running? std::chrono::time_point start_; //!< starting point for clock - double elapsed_; //!< elasped time in [s] + double elapsed_ {0.0}; //!< elasped time in [s] }; //============================================================================== diff --git a/openmc/settings.py b/openmc/settings.py index af64af01b..e0f3634e7 100644 --- a/openmc/settings.py +++ b/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) diff --git a/src/mesh.cpp b/src/mesh.cpp index a9b716ecd..2c3e508fe 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -466,42 +466,38 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& 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 d = xt::zeros({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 d = xt::zeros({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. diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index d7e26fc6b..a0c14077a 100644 --- a/src/mesh_header.F90 +++ b/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 diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index 9d6add0f0..0762846ef 100644 --- a/src/reaction_header.F90 +++ b/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 diff --git a/src/settings.cpp b/src/settings.cpp index 2394597c7..d4cea83bc 100644 --- a/src/settings.cpp +++ b/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."); } diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp index d6dd8b9ac..4c983b267 100644 --- a/src/tallies/filter_energy.cpp +++ b/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 diff --git a/src/timer.cpp b/src/timer.cpp index 46a84a6f1..f1c211446 100644 --- a/src/timer.cpp +++ b/src/timer.cpp @@ -28,7 +28,7 @@ double Timer::elapsed() { if (running_) { std::chrono::duration diff = clock::now() - start_; - return elapsed_ += diff.count(); + return elapsed_ + diff.count(); } else { return elapsed_; }