diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst index d7da3261e..85e5c8ced 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -465,13 +465,15 @@ Functions :return: Return status (negative if an error occurred) :rtype: int -.. c:function:: int openmc_statepoint_write(const char filename[]) +.. c:function:: int openmc_statepoint_write(const char filename[], const bool* write_source) Write a statepoint file :param filename: Name of file to create. If a null pointer is passed, a filename is assigned automatically. :type filename: const char[] + :param write_source: Whether to include the source bank + :type write_source: const bool* :return: Return status (negative if an error occurs) :rtype: int diff --git a/docs/source/devguide/styleguide.rst b/docs/source/devguide/styleguide.rst index 2b366c191..52c6bc3e0 100644 --- a/docs/source/devguide/styleguide.rst +++ b/docs/source/devguide/styleguide.rst @@ -237,6 +237,26 @@ are written differently by convention (e.g., ``E`` for energy). Data members of classes (but not structs) additionally have trailing underscores (e.g., ``a_class_member_``). +The following conventions are used for variables with short names: + +- ``d`` stands for "distance" +- ``E`` stands for "energy" +- ``p`` stands for "particle" +- ``r`` stands for "position" +- ``rx`` stands for "reaction" +- ``u`` stands for "direction" +- ``xs`` stands for "cross section" + +All classes and non-member functions should be declared within the ``openmc`` +namespace. Global variables must be declared in a namespace nested within the +``openmc`` namespace. The following sub-namespaces are in use: + +- ``openmc::data``: Fundamental nuclear data (cross sections, multigroup data, + decay constants, etc.) +- ``openmc::model``: Variables related to geometry, materials, and tallies +- ``openmc::settings``: Global settings / options +- ``openmc::simulation``: Variables used only during a simulation + Accessors and mutators (get and set functions) may be named like variables. These often correspond to actual member variables, but this is not required. For example, ``int count()`` and ``void set_count(int count)``. diff --git a/include/openmc/capi.h b/include/openmc/capi.h index fd664b2a9..fdea48e9e 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -135,17 +135,10 @@ extern "C" { // Global variables extern char openmc_err_msg[256]; - extern int32_t n_cells; - extern int32_t n_lattices; - extern int32_t n_materials; extern int n_nuclides; - extern int32_t n_plots; extern int32_t n_realizations; extern int32_t n_sab_tables; - extern int32_t n_sources; - extern int32_t n_surfaces; extern int32_t n_tallies; - extern int32_t n_universes; // Variables that are shared by necessity (can be removed from public header // later) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 5a889023b..f2c8f0f71 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -39,16 +39,21 @@ constexpr int32_t OP_UNION {std::numeric_limits::max() - 4}; // Global variables //============================================================================== +class Cell; +class Universe; + +namespace model { + extern "C" int32_t n_cells; -class Cell; extern std::vector cells; extern std::unordered_map cell_map; -class Universe; extern std::vector universes; extern std::unordered_map universe_map; +} // namespace model + //============================================================================== //! A geometry primitive that fills all space and contains cells. //============================================================================== @@ -145,13 +150,13 @@ public: virtual ~Cell() {} }; - + class CSGCell : public Cell { public: CSGCell(); - + explicit CSGCell(pugi::xml_node cell_node); bool @@ -163,7 +168,7 @@ public: void to_hdf5(hid_t group_id) const; - + protected: bool contains_simple(Position r, Direction u, int32_t on_surface) const; bool contains_complex(Position r, Direction u, int32_t on_surface) const; @@ -183,6 +188,6 @@ public: }; #endif - + } // namespace openmc #endif // OPENMC_CELL_H diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index afbd9cbc6..51d97d737 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -10,12 +10,24 @@ namespace openmc { +//============================================================================== +// Global variables +//============================================================================== + +namespace model { + extern moab::DagMC* DAG; +} // namespace model + +//============================================================================== +// Non-member functions +//============================================================================== + extern "C" void load_dagmc_geometry(); extern "C" void free_memory_dagmc(); -} +} // namespace openmc #endif // DAGMC diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index b7a75becb..6d918eb40 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -18,6 +18,8 @@ namespace openmc { // Global variables //============================================================================== +namespace simulation { + extern double keff_generation; //!< Single-generation k on each processor extern std::array k_sum; //!< Used to reduce sum and sum_sq extern std::vector entropy; //!< Shannon entropy at each generation @@ -26,6 +28,8 @@ extern xt::xtensor source_frac; //!< Source fraction for UFS extern "C" int64_t n_bank; #pragma omp threadprivate(n_bank) +} // namespace simulation + //============================================================================== // Non-member functions //============================================================================== diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index f045209ee..2f0853e0c 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -13,10 +13,14 @@ namespace openmc { // Global variables //============================================================================== -extern "C" int openmc_root_universe; +namespace model { + +extern "C" int root_universe; extern std::vector overlap_check_count; +} // namespace model + //============================================================================== //! Check for overlapping cells at a particle's position. //============================================================================== diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index cacb87abe..5a5eff0b4 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -31,10 +31,14 @@ enum class LatticeType { //============================================================================== class Lattice; -extern std::vector lattices; +namespace model { + +extern std::vector lattices; extern std::unordered_map lattice_map; +} // namespace model + //============================================================================== //! \class Lattice //! \brief Abstract type for ordered array of universes. diff --git a/include/openmc/material.h b/include/openmc/material.h index b6d7e19cd..fcab023a5 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -14,9 +14,14 @@ namespace openmc { //============================================================================== class Material; + +namespace model { + extern std::vector materials; extern std::unordered_map material_map; +} // namespace model + //============================================================================== //! A substance with constituent nuclides and thermal scattering data //============================================================================== diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 3d31fb489..7fcd58a06 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -17,6 +17,19 @@ namespace openmc { +//============================================================================== +// Global variables +//============================================================================== + +class RegularMesh; + +namespace model { + +extern std::vector> meshes; +extern std::unordered_map mesh_map; + +} // namespace model + //============================================================================== //! Tessellation of n-dimensional Euclidean space by congruent squares or cubes //============================================================================== @@ -117,15 +130,6 @@ extern "C" void read_meshes(pugi::xml_node* root); //! \param[in] group HDF5 group extern "C" void meshes_to_hdf5(hid_t group); -//============================================================================== -// Global variables -//============================================================================== - -extern std::vector> meshes; - -extern std::unordered_map mesh_map; - - } // namespace openmc #endif // OPENMC_MESH_H diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h index c81375511..17be1a18e 100644 --- a/include/openmc/mgxs_interface.h +++ b/include/openmc/mgxs_interface.h @@ -15,6 +15,8 @@ namespace openmc { // Global variables //============================================================================== +namespace data { + extern std::vector nuclides_MG; extern std::vector macro_xs; extern "C" int num_energy_groups; @@ -22,6 +24,8 @@ extern std::vector energy_bins; extern std::vector energy_bin_avg; extern std::vector rev_energy_bins; +} // namespace data + //============================================================================== // Mgxs data loading interface methods //============================================================================== diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index e3d43450d..db30a1a50 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -14,11 +14,15 @@ namespace openmc { // Global variables //============================================================================== +namespace data { + // Minimum/maximum transport energy for each particle type. Order corresponds to // that of the ParticleType enum extern std::array energy_min; extern std::array energy_max; +} // namespace data + //=============================================================================== //! Cached microscopic cross sections for a particular nuclide at the current //! energy diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index 14b7e6cbd..e5648c4a5 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -10,34 +10,29 @@ namespace openmc { -//TODO: Remove energy_bin_avg and material_xs parameters when they reside on +//TODO: Remove material_xs parameters when they reside on // the C-side this should happen after materials, physics, input, and tallies // are brought over //! \brief samples particle behavior after a collision event. //! \param p Particle to operate on -//! \param energy_bin_avg Average energy within each energy bin //! \param material_xs The cross section cache for the current material extern "C" void -collision_mg(Particle* p, const double* energy_bin_avg, - const MaterialMacroXS* material_xs); +collision_mg(Particle* p, const MaterialMacroXS* material_xs); //! \brief samples a reaction type. //! //! Note that there is special logic when suvival biasing is turned on since //! fission and disappearance are treated implicitly. //! \param p Particle to operate on -//! \param energy_bin_avg Average energy within each energy bin //! \param material_xs The cross section cache for the current material void -sample_reaction(Particle* p, const double* energy_bin_avg, - const MaterialMacroXS* material_xs); +sample_reaction(Particle* p, const MaterialMacroXS* material_xs); //! \brief Samples the scattering event //! \param p Particle to operate on -//! \param energy_bin_avg Average energy within each energy bin void -scatter(Particle* p, const double* energy_bin_avg); +scatter(Particle* p); //! \brief Determines the average total, prompt and delayed neutrons produced //! from fission and creates the appropriate bank sites. diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 64c092ce8..1096d2484 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -18,14 +18,14 @@ namespace openmc { // Global variables //=============================================================================== -extern int PLOT_LEVEL_LOWEST; //!< lower bound on plot universe level +class Plot; +namespace model { + +extern std::vector plots; //!< Plot instance container extern std::unordered_map plot_map; //!< map of plot ids to index -extern "C" int32_t n_plots; //!< number of plots in openmc run - -class Plot; -extern std::vector plots; //!< Plot instance container +} // namespace model //=============================================================================== // RGBColor holds color information for plotted objects @@ -124,8 +124,7 @@ void draw_mesh_lines(Plot pl, ImageData& data); //! Write a ppm image to file using a plot object's image data //! \param[in] plot object //! \param[out] image data associated with the plot object -void output_ppm(Plot pl, - const ImageData& data); +void output_ppm(Plot pl, const ImageData& data); //! Get the rgb color for a given particle position in a plot //! \param[in] particle with position for current pixel @@ -141,7 +140,7 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id); //! \param[out] dataset pointer to voxesl data //! \param[out] pointer to memory space of voxel data void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, - hid_t* dset, hid_t* memspace); + hid_t* dset, hid_t* memspace); //! Write a section of the voxel data to hdf5 //! \param[in] voxel slice @@ -149,7 +148,8 @@ void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, //! \param[out] dataset pointer to voxesl data //! \param[out] pointer to data to write void voxel_write_slice(int x, hid_t dspace, hid_t dset, - hid_t memspace, void* buf); + hid_t memspace, void* buf); + //! Close voxel file entities //! \param[in] data space to close //! \param[in] dataset to close diff --git a/include/openmc/source.h b/include/openmc/source.h index efc6f12a6..45f4896b4 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -16,6 +16,18 @@ namespace openmc { +//============================================================================== +// Global variables +//============================================================================== + +class SourceDistribution; + +namespace model { + +extern std::vector external_sources; + +} // namespace model + //============================================================================== //! External source distribution //============================================================================== @@ -40,12 +52,6 @@ private: UPtrDist energy_; //!< Energy distribution }; -//============================================================================== -// Global variables -//============================================================================== - -extern std::vector external_sources; - //============================================================================== // Functions //============================================================================== diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 23c3b61ec..d0ada79ea 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -32,13 +32,15 @@ extern "C" const int BC_PERIODIC; // Global variables //============================================================================== -extern "C" int32_t n_surfaces; - class Surface; -extern std::vector surfaces; +namespace model { + +extern std::vector surfaces; extern std::map surface_map; +} // namespace model + //============================================================================== //! Coordinates for an axis-aligned cube that bounds a geometric object. //============================================================================== diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 24ce5cb89..fd3be0aaf 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -29,7 +29,7 @@ public: // 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; +extern template class std::vector; namespace openmc { @@ -77,13 +77,20 @@ public: // Global variables //============================================================================== -extern "C" int32_t n_filters; +namespace simulation { extern std::vector filter_matches; #pragma omp threadprivate(filter_matches) +} // namespace simulation + +namespace model { + +extern "C" int32_t n_filters; extern std::vector> tally_filters; +} // namespace model + //============================================================================== extern "C" void free_memory_tally_c(); diff --git a/include/openmc/timer.h b/include/openmc/timer.h index 48b3a9e75..8fbcee46b 100644 --- a/include/openmc/timer.h +++ b/include/openmc/timer.h @@ -5,6 +5,27 @@ namespace openmc { +//============================================================================== +// Global variables +//============================================================================== + +class Timer; + +namespace simulation { + +extern Timer time_active; +extern Timer time_bank; +extern Timer time_bank_sample; +extern Timer time_bank_sendrecv; +extern Timer time_finalize; +extern Timer time_inactive; +extern Timer time_initialize; +extern Timer time_tallies; +extern Timer time_total; +extern Timer time_transport; + +} // namespace simulation + //============================================================================== //! Class for measuring time elapsed //============================================================================== @@ -34,21 +55,6 @@ private: double elapsed_ {0.0}; //!< elasped time in [s] }; -//============================================================================== -// Global variables -//============================================================================== - -extern Timer time_active; -extern Timer time_bank; -extern Timer time_bank_sample; -extern Timer time_bank_sendrecv; -extern Timer time_finalize; -extern Timer time_inactive; -extern Timer time_initialize; -extern Timer time_tallies; -extern Timer time_total; -extern Timer time_transport; - //============================================================================== // Non-member functions //============================================================================== diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 5be31ce86..bedde6e16 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -52,7 +52,7 @@ _dll.openmc_simulation_init.restype = c_int _dll.openmc_simulation_init.errcheck = _error_handler _dll.openmc_simulation_finalize.restype = c_int _dll.openmc_simulation_finalize.errcheck = _error_handler -_dll.openmc_statepoint_write.argtypes = [POINTER(c_char_p), POINTER(c_bool)] +_dll.openmc_statepoint_write.argtypes = [c_char_p, POINTER(c_bool)] _dll.openmc_statepoint_write.restype = c_int _dll.openmc_statepoint_write.errcheck = _error_handler diff --git a/src/cell.cpp b/src/cell.cpp index 65582d780..c3d806aef 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -21,6 +21,8 @@ namespace openmc { // Global variables //============================================================================== +namespace model { + int32_t n_cells {0}; std::vector cells; @@ -29,6 +31,8 @@ std::unordered_map cell_map; std::vector universes; std::unordered_map universe_map; +} // namespace model + //============================================================================== //! Convert region specification string to integer tokens. //! @@ -197,7 +201,7 @@ Universe::to_hdf5(hid_t universes_group) const // Write the contained cells. if (cells_.size() > 0) { std::vector cell_ids; - for (auto i_cell : cells_) cell_ids.push_back(cells[i_cell]->id_); + for (auto i_cell : cells_) cell_ids.push_back(model::cells[i_cell]->id_); write_dataset(group, "cells", cell_ids); } @@ -209,7 +213,7 @@ Universe::to_hdf5(hid_t universes_group) const //============================================================================== CSGCell::CSGCell() {} // empty constructor - + CSGCell::CSGCell(pugi::xml_node cell_node) { if (check_for_node(cell_node, "id")) { @@ -314,7 +318,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // Convert user IDs to surface indices. for (auto& r : region_) { if (r < OP_UNION) { - r = copysign(surface_map[abs(r)] + 1, r); + r = copysign(model::surface_map[abs(r)] + 1, r); } } @@ -418,7 +422,7 @@ CSGCell::distance(Position r, Direction u, int32_t on_surface) const // Calculate the distance to this surface. // Note the off-by-one indexing bool coincident {token == on_surface}; - double d {surfaces[abs(token)-1]->distance(r, u, coincident)}; + double d {model::surfaces[abs(token)-1]->distance(r, u, coincident)}; // Check if this distance is the new minimum. if (d < min_dist) { @@ -446,7 +450,7 @@ CSGCell::to_hdf5(hid_t cell_group) const write_string(group, "name", name_, false); } - write_dataset(group, "universe", universes[universe_]->id_); + write_dataset(group, "universe", model::universes[universe_]->id_); // Write the region specification. if (!region_.empty()) { @@ -464,7 +468,7 @@ CSGCell::to_hdf5(hid_t cell_group) const } else { // Note the off-by-one indexing region_spec << " " - << copysign(surfaces[abs(token)-1]->id_, token); + << copysign(model::surfaces[abs(token)-1]->id_, token); } } write_string(group, "region", region_spec.str(), false); @@ -476,7 +480,7 @@ CSGCell::to_hdf5(hid_t cell_group) const std::vector mat_ids; for (auto i_mat : material_) { if (i_mat != MATERIAL_VOID) { - mat_ids.push_back(materials[i_mat]->id_); + mat_ids.push_back(model::materials[i_mat]->id_); } else { mat_ids.push_back(MATERIAL_VOID); } @@ -494,7 +498,7 @@ CSGCell::to_hdf5(hid_t cell_group) const } else if (type_ == FILL_UNIVERSE) { write_dataset(group, "fill_type", "universe"); - write_dataset(group, "fill", universes[fill_]->id_); + write_dataset(group, "fill", model::universes[fill_]->id_); if (translation_ != Position(0, 0, 0)) { write_dataset(group, "translation", translation_); } @@ -505,7 +509,7 @@ CSGCell::to_hdf5(hid_t cell_group) const } else if (type_ == FILL_LATTICE) { write_dataset(group, "fill_type", "lattice"); - write_dataset(group, "lattice", lattices[fill_]->id_); + write_dataset(group, "lattice", model::lattices[fill_]->id_); } close_group(group); @@ -527,7 +531,7 @@ CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const return false; } else { // Note the off-by-one indexing - bool sense = surfaces[abs(token)-1]->sense(r, u); + bool sense = model::surfaces[abs(token)-1]->sense(r, u); if (sense != (token > 0)) {return false;} } } @@ -569,7 +573,7 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const stack[i_stack] = false; } else { // Note the off-by-one indexing - bool sense = surfaces[abs(token)-1]->sense(r, u); + bool sense = model::surfaces[abs(token)-1]->sense(r, u); stack[i_stack] = (sense == (token > 0)); } } @@ -609,10 +613,10 @@ DAGCell::distance(Position r, Direction u, int32_t on_surface) const } else { // indicate that particle is lost surf_idx = -1; } - + return {dist, surf_idx}; } - + bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const { moab::ErrorCode rval; @@ -622,7 +626,7 @@ bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const double pnt[3] = {r.x, r.y, r.z}; double dir[3] = {u.x, u.y, u.z}; rval = dagmc_ptr_->point_in_volume(vol, pnt, result, dir); - MB_CHK_ERR_CONT(rval); + MB_CHK_ERR_CONT(rval); return result; } @@ -638,23 +642,23 @@ extern "C" void read_cells(pugi::xml_node* node) { // Count the number of cells. - for (pugi::xml_node cell_node: node->children("cell")) {n_cells++;} - if (n_cells == 0) { + for (pugi::xml_node cell_node: node->children("cell")) {model::n_cells++;} + if (model::n_cells == 0) { fatal_error("No cells found in geometry.xml!"); } // Loop over XML cell elements and populate the array. - cells.reserve(n_cells); - for (pugi::xml_node cell_node: node->children("cell")) { - cells.push_back(new CSGCell(cell_node)); + model::cells.reserve(model::n_cells); + for (pugi::xml_node cell_node : node->children("cell")) { + model::cells.push_back(new CSGCell(cell_node)); } // Fill the cell map. - for (int i = 0; i < cells.size(); i++) { - int32_t id = cells[i]->id_; - auto search = cell_map.find(id); - if (search == cell_map.end()) { - cell_map[id] = i; + for (int i = 0; i < model::cells.size(); i++) { + int32_t id = model::cells[i]->id_; + auto search = model::cell_map.find(id); + if (search == model::cell_map.end()) { + model::cell_map[id] = i; } else { std::stringstream err_msg; err_msg << "Two or more cells use the same unique ID: " << id; @@ -663,22 +667,24 @@ read_cells(pugi::xml_node* node) } // Populate the Universe vector and map. - for (int i = 0; i < cells.size(); i++) { - int32_t uid = cells[i]->universe_; - auto it = universe_map.find(uid); - if (it == universe_map.end()) { - universes.push_back(new Universe()); - universes.back()->id_ = uid; - universes.back()->cells_.push_back(i); - universe_map[uid] = universes.size() - 1; + for (int i = 0; i < model::cells.size(); i++) { + int32_t uid = model::cells[i]->universe_; + auto it = model::universe_map.find(uid); + if (it == model::universe_map.end()) { + model::universes.push_back(new Universe()); + model::universes.back()->id_ = uid; + model::universes.back()->cells_.push_back(i); + model::universe_map[uid] = model::universes.size() - 1; } else { - universes[it->second]->cells_.push_back(i); + model::universes[it->second]->cells_.push_back(i); } } - universes.shrink_to_fit(); + model::universes.shrink_to_fit(); // Allocate the cell overlap count if necessary. - if (settings::check_overlaps) overlap_check_count.resize(n_cells, 0); + if (settings::check_overlaps) { + model::overlap_check_count.resize(model::cells.size(), 0); + } } //============================================================================== @@ -688,9 +694,9 @@ read_cells(pugi::xml_node* node) extern "C" int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n) { - if (index >= 1 && index <= cells.size()) { + if (index >= 1 && index <= model::cells.size()) { //TODO: off-by-one - Cell& c {*cells[index - 1]}; + Cell& c {*model::cells[index - 1]}; *type = c.type_; if (c.type_ == FILL_MATERIAL) { *indices = c.material_.data(); @@ -710,9 +716,9 @@ extern "C" int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices) { - if (index >= 1 && index <= cells.size()) { + if (index >= 1 && index <= model::cells.size()) { //TODO: off-by-one - Cell& c {*cells[index - 1]}; + Cell& c {*model::cells[index - 1]}; if (type == FILL_MATERIAL) { c.type_ = FILL_MATERIAL; c.material_.clear(); @@ -720,7 +726,7 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n, int i_mat = indices[i]; if (i_mat == MATERIAL_VOID) { c.material_.push_back(MATERIAL_VOID); - } else if (i_mat >= 1 && i_mat <= materials.size()) { + } else if (i_mat >= 1 && i_mat <= model::materials.size()) { //TODO: off-by-one c.material_.push_back(i_mat - 1); } else { @@ -745,9 +751,9 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n, extern "C" int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance) { - if (index >= 1 && index <= cells.size()) { + if (index >= 1 && index <= model::cells.size()) { //TODO: off-by-one - Cell& c {*cells[index - 1]}; + Cell& c {*model::cells[index - 1]}; if (instance) { if (*instance >= 0 && *instance < c.sqrtkT_.size()) { @@ -775,7 +781,7 @@ openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance) //============================================================================== extern "C" { - Cell* cell_pointer(int32_t cell_ind) {return cells[cell_ind];} + Cell* cell_pointer(int32_t cell_ind) {return model::cells[cell_ind];} int32_t cell_id(Cell* c) {return c->id_;} @@ -785,9 +791,9 @@ extern "C" { c->id_ = id; // Find the index of this cell and update the cell map. - for (int i = 0; i < cells.size(); i++) { - if (cells[i] == c) { - cell_map[id] = i; + for (int i = 0; i < model::cells.size(); i++) { + if (model::cells[i] == c) { + model::cell_map[id] = i; break; } } @@ -830,17 +836,17 @@ extern "C" { void extend_cells_c(int32_t n) { - cells.reserve(cells.size() + n); + model::cells.reserve(model::cells.size() + n); for (int32_t i = 0; i < n; i++) { - cells.push_back(new CSGCell()); + model::cells.push_back(new CSGCell()); } - n_cells = cells.size(); + model::n_cells = model::cells.size(); } - int32_t universe_id(int i_univ) {return universes[i_univ]->id_;} + int32_t universe_id(int i_univ) {return model::universes[i_univ]->id_;} void universes_to_hdf5(hid_t universes_group) - {for (Universe* u : universes) u->to_hdf5(universes_group);} + {for (Universe* u : model::universes) u->to_hdf5(universes_group);} } diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp index 46371601b..0af584c1b 100644 --- a/src/cmfd_execute.cpp +++ b/src/cmfd_execute.cpp @@ -24,7 +24,7 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies, openmc_source_bank(&source_bank, &n); // Get source counts in each mesh bin / energy bin - auto& m = meshes.at(settings::index_cmfd_mesh); + auto& m = model::meshes.at(settings::index_cmfd_mesh); xt::xarray counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside); // Copy data from the xarray into the source counts array diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 7c1127867..e6ccee24d 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -1,5 +1,6 @@ - #include "openmc/dagmc.h" + +#include "openmc/cell.h" #include "openmc/error.h" #include "openmc/string_utils.h" #include "openmc/settings.h" @@ -13,20 +14,24 @@ namespace openmc { +namespace model { + moab::DagMC* DAG; +} // namespace model + void load_dagmc_geometry() { - if (!DAG) { - DAG = new moab::DagMC(); + if (!model::DAG) { + model::DAG = new moab::DagMC(); } int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC - moab::ErrorCode rval = DAG->load_file("dagmc.h5m"); + moab::ErrorCode rval = model::DAG->load_file("dagmc.h5m"); MB_CHK_ERR_CONT(rval); - rval = DAG->init_OBBTree(); + rval = model::DAG->init_OBBTree(); MB_CHK_ERR_CONT(rval); std::vector prop_keywords; @@ -34,48 +39,45 @@ void load_dagmc_geometry() prop_keywords.push_back("boundary"); std::map ph; - DAG->parse_properties(prop_keywords, ph, ":"); + model::DAG->parse_properties(prop_keywords, ph, ":"); MB_CHK_ERR_CONT(rval); // initialize cell objects - n_cells = DAG->num_entities(3); + model::n_cells = model::DAG->num_entities(3); - // Allocate the cell overlap count if necessary. - if (settings::check_overlaps) overlap_check_count.resize(n_cells, 0); - - for (int i = 0; i < n_cells; i++) { - moab::EntityHandle vol_handle = DAG->entity_by_index(3, i+1); + for (int i = 0; i < model::n_cells; i++) { + moab::EntityHandle vol_handle = model::DAG->entity_by_index(3, i+1); // set cell ids using global IDs DAGCell* c = new DAGCell(); - c->id_ = DAG->id_by_index(3, i+1); - c->dagmc_ptr_ = DAG; + c->id_ = model::DAG->id_by_index(3, i+1); + c->dagmc_ptr_ = model::DAG; c->universe_ = dagmc_univ_id; // set to zero for now c->fill_ = C_NONE; // no fill, single universe - cells.push_back(c); - cell_map[c->id_] = i; + model::cells.push_back(c); + model::cell_map[c->id_] = i; // Populate the Universe vector and dict - auto it = universe_map.find(dagmc_univ_id); - if (it == universe_map.end()) { - universes.push_back(new Universe()); - universes.back()-> id_ = dagmc_univ_id; - universes.back()->cells_.push_back(i); - universe_map[dagmc_univ_id] = universes.size() - 1; + auto it = model::universe_map.find(dagmc_univ_id); + if (it == model::universe_map.end()) { + model::universes.push_back(new Universe()); + model::universes.back()-> id_ = dagmc_univ_id; + model::universes.back()->cells_.push_back(i); + model::universe_map[dagmc_univ_id] = model::universes.size() - 1; } else { - universes[it->second]->cells_.push_back(i); + model::universes[it->second]->cells_.push_back(i); } - if (DAG->is_implicit_complement(vol_handle)) { + if (model::DAG->is_implicit_complement(vol_handle)) { // assuming implicit complement is void for now c->material_.push_back(MATERIAL_VOID); continue; } - if (DAG->has_prop(vol_handle, "mat")){ + if (model::DAG->has_prop(vol_handle, "mat")){ std::string mat_value; - rval = DAG->prop_value(vol_handle, "mat", mat_value); + rval = model::DAG->prop_value(vol_handle, "mat", mat_value); MB_CHK_ERR_CONT(rval); to_lower(mat_value); @@ -91,21 +93,26 @@ void load_dagmc_geometry() } } + // Allocate the cell overlap count if necessary. + if (settings::check_overlaps) { + model::overlap_check_count.resize(model::cells.size(), 0); + } + // initialize surface objects - n_surfaces = DAG->num_entities(2); - surfaces.resize(n_surfaces); + int n_surfaces = model::DAG->num_entities(2); + model::surfaces.resize(n_surfaces); for (int i = 0; i < n_surfaces; i++) { - moab::EntityHandle surf_handle = DAG->entity_by_index(2, i+1); + moab::EntityHandle surf_handle = model::DAG->entity_by_index(2, i+1); // set cell ids using global IDs DAGSurface* s = new DAGSurface(); - s->id_ = DAG->id_by_index(2, i+1); - s->dagmc_ptr_ = DAG; + s->id_ = model::DAG->id_by_index(2, i+1); + s->dagmc_ptr_ = model::DAG; - if (DAG->has_prop(surf_handle, "boundary")) { + if (model::DAG->has_prop(surf_handle, "boundary")) { std::string bc_value; - rval = DAG->prop_value(surf_handle, "boundary", bc_value); + rval = model::DAG->prop_value(surf_handle, "boundary", bc_value); MB_CHK_ERR_CONT(rval); to_lower(bc_value); @@ -128,8 +135,8 @@ void load_dagmc_geometry() } // add to global array and map - surfaces[i] = s; - surface_map[s->id_] = s->id_; + model::surfaces[i] = s; + model::surface_map[s->id_] = s->id_; } return; @@ -137,7 +144,7 @@ void load_dagmc_geometry() void free_memory_dagmc() { - delete DAG; + delete model::DAG; } } diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 927b014b5..64e0c4e42 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -30,11 +30,15 @@ namespace openmc { // Global variables //============================================================================== +namespace simulation { + double keff_generation; std::array k_sum; std::vector entropy; xt::xtensor source_frac; +} // namespace simulation + //============================================================================== // Non-member functions //============================================================================== @@ -44,15 +48,15 @@ void calculate_generation_keff() auto gt = global_tallies(); // Get keff for this generation by subtracting off the starting value - keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - keff_generation; + simulation::keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - simulation::keff_generation; double keff_reduced; #ifdef OPENMC_MPI // Combine values across all processors - MPI_Allreduce(&keff_generation, &keff_reduced, 1, MPI_DOUBLE, + MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE, MPI_SUM, mpi::intracomm); #else - keff_reduced = keff_generation; + keff_reduced = simulation::keff_generation; #endif // Normalize single batch estimate of k @@ -63,7 +67,7 @@ void calculate_generation_keff() void synchronize_bank() { - time_bank.start(); + simulation::time_bank.start(); // Get pointers to source/fission bank Bank* source_bank; @@ -83,21 +87,21 @@ void synchronize_bank() #ifdef OPENMC_MPI int64_t start = 0; - MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm); + MPI_Exscan(&simulation::n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm); // While we would expect the value of start on rank 0 to be 0, the MPI // standard says that the receive buffer on rank 0 is undefined and not // significant if (mpi::rank == 0) start = 0; - int64_t finish = start + n_bank; + int64_t finish = start + simulation::n_bank; int64_t total = finish; MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm); #else int64_t start = 0; - int64_t finish = n_bank; - int64_t total = n_bank; + int64_t finish = simulation::n_bank; + int64_t total = simulation::n_bank; #endif // If there are not that many particles per generation, it's possible that no @@ -105,7 +109,7 @@ void synchronize_bank() // extra logic to treat this circumstance, we really want to ensure the user // runs enough particles to avoid this in the first place. - if (n_bank == 0) { + if (simulation::n_bank == 0) { fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank)); } @@ -127,7 +131,7 @@ void synchronize_bank() } double p_sample = static_cast(sites_needed) / total; - time_bank_sample.start(); + simulation::time_bank_sample.start(); // ========================================================================== // SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES @@ -136,7 +140,7 @@ void synchronize_bank() int64_t index_temp = 0; std::vector temp_sites(3*simulation::work); - for (int64_t i = 0; i < n_bank; ++i) { + for (int64_t i = 0; i < simulation::n_bank; ++i) { // If there are less than n_particles particles banked, automatically add // int(n_particles/total) sites to temp_sites. For example, if you need // 1000 and 300 were banked, this would add 3 source sites per banked site @@ -191,7 +195,7 @@ void synchronize_bank() // fission bank sites_needed = settings::n_particles - finish; for (int i = 0; i < sites_needed; ++i) { - temp_sites[index_temp] = fission_bank[n_bank - sites_needed + i]; + temp_sites[index_temp] = fission_bank[simulation::n_bank - sites_needed + i]; ++index_temp; } } @@ -200,8 +204,8 @@ void synchronize_bank() finish = simulation::work_index[mpi::rank + 1]; } - time_bank_sample.stop(); - time_bank_sendrecv.start(); + simulation::time_bank_sample.stop(); + simulation::time_bank_sendrecv.start(); #ifdef OPENMC_MPI // ========================================================================== @@ -299,8 +303,8 @@ void synchronize_bank() std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, source_bank); #endif - time_bank_sendrecv.stop(); - time_bank.stop(); + simulation::time_bank_sendrecv.stop(); + simulation::time_bank.stop(); } void calculate_average_keff() @@ -320,11 +324,11 @@ void calculate_average_keff() simulation::keff = simulation::k_generation[i]; } else { // Sample mean of keff - k_sum[0] += simulation::k_generation[i]; - k_sum[1] += std::pow(simulation::k_generation[i], 2); + simulation::k_sum[0] += simulation::k_generation[i]; + simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); // Determine mean - simulation::keff = k_sum[0] / n; + simulation::keff = simulation::k_sum[0] / n; if (n > 1) { double t_value; @@ -337,7 +341,7 @@ void calculate_average_keff() } // Standard deviation of the sample mean of k - simulation::keff_std = t_value * std::sqrt((k_sum[1]/n - + simulation::keff_std = t_value * std::sqrt((simulation::k_sum[1]/n - std::pow(simulation::keff, 2)) / (n - 1)); } } @@ -493,7 +497,7 @@ int openmc_get_keff(double* k_combined) void shannon_entropy() { // Get pointer to entropy mesh - auto& m = meshes[settings::index_entropy_mesh]; + auto& m = model::meshes[settings::index_entropy_mesh]; // Get pointer to fission bank Bank* fission_bank; @@ -503,7 +507,7 @@ void shannon_entropy() // Get source weight in each mesh bin bool sites_outside; xt::xtensor p = m->count_sites( - n_bank, fission_bank, 0, nullptr, &sites_outside); + simulation::n_bank, fission_bank, 0, nullptr, &sites_outside); // display warning message if there were sites outside entropy box if (sites_outside) { @@ -523,20 +527,20 @@ void shannon_entropy() } // Add value to vector - entropy.push_back(H); + simulation::entropy.push_back(H); } } void ufs_count_sites() { - auto &m = meshes[settings::index_ufs_mesh]; + auto &m = model::meshes[settings::index_ufs_mesh]; if (simulation::current_batch == 1 && simulation::current_gen == 1) { // On the first generation, just assume that the source is already evenly // distributed so that effectively the production of fission sites is not // biased - auto s = xt::view(source_frac, xt::all()); + auto s = xt::view(simulation::source_frac, xt::all()); s = m->volume_frac_; } else { @@ -547,7 +551,7 @@ void ufs_count_sites() // count number of source sites in each ufs mesh cell bool sites_outside; - source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr, + simulation::source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr, &sites_outside); // Check for sites outside of the mesh @@ -558,12 +562,12 @@ void ufs_count_sites() #ifdef OPENMC_MPI // Send source fraction to all processors int n_bins = xt::prod(m->shape_)(); - MPI_Bcast(source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm); + MPI_Bcast(simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm); #endif // Normalize to total weight to get fraction of source in each cell - double total = xt::sum(source_frac)(); - source_frac /= total; + double total = xt::sum(simulation::source_frac)(); + simulation::source_frac /= total; // Since the total starting weight is not equal to n_particles, we need to // renormalize the weight of the source sites @@ -575,7 +579,7 @@ void ufs_count_sites() double ufs_get_weight(const Particle* p) { - auto& m = meshes[settings::index_ufs_mesh]; + auto& m = model::meshes[settings::index_ufs_mesh]; // Determine indices on ufs mesh for current location // TODO: off by one @@ -585,8 +589,8 @@ double ufs_get_weight(const Particle* p) fatal_error("Source site outside UFS mesh!"); } - if (source_frac(mesh_bin) != 0.0) { - return m->volume_frac_ / source_frac(mesh_bin); + if (simulation::source_frac(mesh_bin) != 0.0) { + return m->volume_frac_ / simulation::source_frac(mesh_bin); } else { return 1.0; } @@ -598,7 +602,7 @@ extern "C" void write_eigenvalue_hdf5(hid_t group) write_dataset(group, "generations_per_batch", settings::gen_per_batch); write_dataset(group, "k_generation", simulation::k_generation); if (settings::entropy_on) { - write_dataset(group, "entropy", entropy); + write_dataset(group, "entropy", simulation::entropy); } write_dataset(group, "k_col_abs", simulation::k_col_abs); write_dataset(group, "k_col_tra", simulation::k_col_tra); @@ -615,7 +619,7 @@ extern "C" void read_eigenvalue_hdf5(hid_t group) simulation::k_generation.resize(n); read_dataset(group, "k_generation", simulation::k_generation); if (settings::entropy_on) { - read_dataset(group, "entropy", entropy); + read_dataset(group, "entropy", simulation::entropy); } read_dataset(group, "k_col_abs", simulation::k_col_abs); read_dataset(group, "k_col_tra", simulation::k_col_tra); @@ -628,14 +632,14 @@ extern "C" void read_eigenvalue_hdf5(hid_t group) extern "C" double entropy_c(int i) { - return entropy.at(i - 1); + return simulation::entropy.at(i - 1); } extern "C" void entropy_clear() { - entropy.clear(); + simulation::entropy.clear(); } -extern "C" void k_sum_reset() { k_sum.fill(0.0); } +extern "C" void k_sum_reset() { simulation::k_sum.fill(0.0); } } // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 94e3d1a7f..5885470c1 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -75,10 +75,10 @@ int openmc_finalize() simulation::satisfy_triggers = false; simulation::total_gen = 0; - energy_max = {INFTY, INFTY}; - energy_min = {0.0, 0.0}; + data::energy_max = {INFTY, INFTY}; + data::energy_min = {0.0, 0.0}; n_tallies = 0; - openmc_root_universe = -1; + model::root_universe = -1; openmc_set_seed(DEFAULT_SEED); // Deallocate arrays @@ -113,7 +113,7 @@ int openmc_reset() simulation::k_col_abs = 0.0; simulation::k_col_tra = 0.0; simulation::k_abs_tra = 0.0; - k_sum = {0.0, 0.0}; + simulation::k_sum = {0.0, 0.0}; // Reset timers reset_timers(); diff --git a/src/geometry.cpp b/src/geometry.cpp index cff5657a5..fe3d52c68 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -14,8 +14,21 @@ namespace openmc { +//============================================================================== +// Global variables +//============================================================================== + + +namespace model { + +int root_universe {-1}; + std::vector overlap_check_count; +} // namespace model + +//============================================================================== +// Non-member functions //============================================================================== extern "C" bool @@ -24,21 +37,21 @@ check_cell_overlap(Particle* p) { // Loop through each coordinate level for (int j = 0; j < n_coord; j++) { - Universe& univ = *universes[p->coord[j].universe]; + Universe& univ = *model::universes[p->coord[j].universe]; int n = univ.cells_.size(); // Loop through each cell on this level for (auto index_cell : univ.cells_) { - Cell& c = *cells[index_cell]; + Cell& c = *model::cells[index_cell]; if (c.contains(p->coord[j].xyz, p->coord[j].uvw, p->surface)) { if (index_cell != p->coord[j].cell) { std::stringstream err_msg; err_msg << "Overlapping cells detected: " << c.id_ << ", " - << cells[p->coord[j].cell]->id_ << " on universe " + << model::cells[p->coord[j].cell]->id_ << " on universe " << univ.id_; fatal_error(err_msg); } - ++overlap_check_count[index_cell]; + ++model::overlap_check_count[index_cell]; } } } @@ -57,8 +70,8 @@ find_cell(Particle* p, int search_surf) { // Determine universe (if not yet set, use root universe) int i_universe = p->coord[p->n_coord-1].universe; if (i_universe == C_NONE) { - p->coord[p->n_coord-1].universe = openmc_root_universe; - i_universe = openmc_root_universe; + p->coord[p->n_coord-1].universe = model::root_universe; + i_universe = model::root_universe; } // If a surface was indicated, only search cells from the neighbor list of @@ -66,12 +79,12 @@ find_cell(Particle* p, int search_surf) { // the positive or negative side of the surface should be searched. const std::vector* search_cells; if (search_surf > 0) { - search_cells = &surfaces[search_surf-1]->neighbor_pos_; + search_cells = &model::surfaces[search_surf-1]->neighbor_pos_; } else if (search_surf < 0) { - search_cells = &surfaces[-search_surf-1]->neighbor_neg_; + search_cells = &model::surfaces[-search_surf-1]->neighbor_neg_; } else { // No surface was indicated, search all cells in the universe. - search_cells = &universes[i_universe]->cells_; + search_cells = &model::universes[i_universe]->cells_; } // Find which cell of this universe the particle is in. @@ -81,17 +94,17 @@ find_cell(Particle* p, int search_surf) { i_cell = (*search_cells)[i]; // Make sure the search cell is in the same universe. - if (cells[i_cell]->universe_ != i_universe) continue; + if (model::cells[i_cell]->universe_ != i_universe) continue; Position r {p->coord[p->n_coord-1].xyz}; Direction u {p->coord[p->n_coord-1].uvw}; int32_t surf = p->surface; - if (cells[i_cell]->contains(r, u, surf)) { + if (model::cells[i_cell]->contains(r, u, surf)) { p->coord[p->n_coord-1].cell = i_cell; if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; - msg << " Entering cell " << cells[i_cell]->id_; + msg << " Entering cell " << model::cells[i_cell]->id_; write_message(msg, 1); } found = true; @@ -100,7 +113,7 @@ find_cell(Particle* p, int search_surf) { } if (found) { - Cell& c {*cells[i_cell]}; + Cell& c {*model::cells[i_cell]}; if (c.type_ == FILL_MATERIAL) { //======================================================================= //! Found a material cell which means this is the lowest coord level. @@ -109,11 +122,11 @@ find_cell(Particle* p, int search_surf) { if (c.material_.size() > 1 || c.sqrtkT_.size() > 1) { int offset = 0; for (int i = 0; i < p->n_coord; i++) { - Cell& c_i {*cells[p->coord[i].cell]}; + Cell& c_i {*model::cells[p->coord[i].cell]}; if (c_i.type_ == FILL_UNIVERSE) { offset += c_i.offset_[c.distribcell_index_]; } else if (c_i.type_ == FILL_LATTICE) { - Lattice& lat {*lattices[p->coord[i+1].lattice-1]}; + Lattice& lat {*model::lattices[p->coord[i+1].lattice-1]}; int i_xyz[3] {p->coord[i+1].lattice_x, p->coord[i+1].lattice_y, p->coord[i+1].lattice_z}; @@ -198,7 +211,7 @@ find_cell(Particle* p, int search_surf) { //======================================================================== //! Found a lower lattice, update this coord level then search the next. - Lattice& lat {*lattices[c.fill_]}; + Lattice& lat {*model::lattices[c.fill_]}; // Determine lattice indices. Position r {p->coord[p->n_coord-1].xyz}; @@ -251,7 +264,7 @@ find_cell(Particle* p, int search_surf) { extern "C" void cross_lattice(Particle* p, int lattice_translation[3]) { - Lattice& lat {*lattices[p->coord[p->n_coord-1].lattice-1]}; + Lattice& lat {*model::lattices[p->coord[p->n_coord-1].lattice-1]}; if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; @@ -326,7 +339,7 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, for (int i = 0; i < p->n_coord; i++) { Position r {p->coord[i].xyz}; Direction u {p->coord[i].uvw}; - Cell& c {*cells[p->coord[i].cell]}; + Cell& c {*model::cells[p->coord[i].cell]}; // Find the oncoming surface in this cell and the distance to it. auto surface_distance = c.distance(r, u, p->surface); @@ -335,7 +348,7 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, // Find the distance to the next lattice tile crossing. if (p->coord[i].lattice != F90_NONE) { - Lattice& lat {*lattices[p->coord[i].lattice-1]}; + Lattice& lat {*model::lattices[p->coord[i].lattice-1]}; std::array i_xyz {p->coord[i].lattice_x, p->coord[i].lattice_y, p->coord[i].lattice_z}; //TODO: refactor so both lattice use the same position argument (which @@ -377,7 +390,7 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, *surface_crossed = level_surf_cross; } else { Position r_hit = r + d_surf * u; - Surface& surf {*surfaces[std::abs(level_surf_cross)-1]}; + Surface& surf {*model::surfaces[std::abs(level_surf_cross)-1]}; Direction norm = surf.normal(r_hit); if (u.dot(norm) > 0) { *surface_crossed = std::abs(level_surf_cross); diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index e3934480b..c326744eb 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -24,15 +24,15 @@ void adjust_indices() { // Adjust material/fill idices. - for (Cell* c : cells) { + for (Cell* c : model::cells) { if (c->fill_ != C_NONE) { int32_t id = c->fill_; - auto search_univ = universe_map.find(id); - auto search_lat = lattice_map.find(id); - if (search_univ != universe_map.end()) { + auto search_univ = model::universe_map.find(id); + auto search_lat = model::lattice_map.find(id); + if (search_univ != model::universe_map.end()) { c->type_ = FILL_UNIVERSE; c->fill_ = search_univ->second; - } else if (search_lat != lattice_map.end()) { + } else if (search_lat != model::lattice_map.end()) { c->type_ = FILL_LATTICE; c->fill_ = search_lat->second; } else { @@ -46,8 +46,8 @@ adjust_indices() for (auto it = c->material_.begin(); it != c->material_.end(); it++) { int32_t mid = *it; if (mid != MATERIAL_VOID) { - auto search = material_map.find(mid); - if (search != material_map.end()) { + auto search = model::material_map.find(mid); + if (search != model::material_map.end()) { *it = search->second; } else { std::stringstream err_msg; @@ -61,9 +61,9 @@ adjust_indices() } // Change cell.universe values from IDs to indices. - for (Cell* c : cells) { - auto search = universe_map.find(c->universe_); - if (search != universe_map.end()) { + for (Cell* c : model::cells) { + auto search = model::universe_map.find(c->universe_); + if (search != model::universe_map.end()) { c->universe_ = search->second; } else { std::stringstream err_msg; @@ -74,7 +74,7 @@ adjust_indices() } // Change all lattice universe values from IDs to indices. - for (Lattice* l : lattices) { + for (Lattice* l : model::lattices) { l->adjust_indices(); } } @@ -84,7 +84,7 @@ adjust_indices() void assign_temperatures() { - for (Cell* c : cells) { + for (Cell* c : model::cells) { // Ignore non-material cells and cells with defined temperature. if (c->material_.size() == 0) continue; if (c->sqrtkT_.size() > 0) continue; @@ -96,9 +96,9 @@ assign_temperatures() c->sqrtkT_.push_back(0); } else { - if (materials[i_mat]->temperature_ >= 0) { + if (model::materials[i_mat]->temperature_ >= 0) { // This material has a default temperature; use that value. - auto T = materials[i_mat]->temperature_; + auto T = model::materials[i_mat]->temperature_; c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * T)); } else { // Use the global default temperature. @@ -117,12 +117,12 @@ find_root_universe() { // Find all the universes listed as a cell fill. std::unordered_set fill_univ_ids; - for (Cell* c : cells) { + for (Cell* c : model::cells) { fill_univ_ids.insert(c->fill_); } // Find all the universes contained in a lattice. - for (Lattice* lat : lattices) { + for (Lattice* lat : model::lattices) { for (auto it = lat->begin(); it != lat->end(); ++it) { fill_univ_ids.insert(*it); } @@ -134,8 +134,8 @@ find_root_universe() // Figure out which universe is not in the set. This is the root universe. bool root_found {false}; int32_t root_univ; - for (int32_t i = 0; i < universes.size(); i++) { - auto search = fill_univ_ids.find(universes[i]->id_); + for (int32_t i = 0; i < model::universes.size(); i++) { + auto search = fill_univ_ids.find(model::universes[i]->id_); if (search == fill_univ_ids.end()) { if (root_found) { fatal_error("Two or more universes are not used as fill universes, so " @@ -160,21 +160,21 @@ neighbor_lists() { write_message("Building neighboring cells lists for each surface...", 6); - for (int i = 0; i < cells.size(); i++) { - for (auto token : cells[i]->region_) { + for (int i = 0; i < model::cells.size(); i++) { + for (auto token : model::cells[i]->region_) { // Skip operator tokens. if (std::abs(token) >= OP_UNION) continue; // This token is a surface index. Add the cell to the surface's list. if (token > 0) { - surfaces[std::abs(token)-1]->neighbor_pos_.push_back(i); + model::surfaces[std::abs(token)-1]->neighbor_pos_.push_back(i); } else { - surfaces[std::abs(token)-1]->neighbor_neg_.push_back(i); + model::surfaces[std::abs(token)-1]->neighbor_neg_.push_back(i); } } } - for (Surface* surf : surfaces) { + for (Surface* surf : model::surfaces) { surf->neighbor_pos_.shrink_to_fit(); surf->neighbor_neg_.shrink_to_fit(); } @@ -187,7 +187,7 @@ prepare_distribcell() { // Find all cells listed in a DistribcellFilter. std::unordered_set distribcells; - for (auto& filt : tally_filters) { + for (auto& filt : model::tally_filters) { auto* distrib_filt = dynamic_cast(filt.get()); if (distrib_filt) { distribcells.insert(distrib_filt->cell_); @@ -196,8 +196,8 @@ prepare_distribcell() // Find all cells with distributed materials or temperatures. Make sure that // the number of materials/temperatures matches the number of cell instances. - for (int i = 0; i < cells.size(); i++) { - Cell& c {*cells[i]}; + for (int i = 0; i < model::cells.size(); i++) { + Cell& c {*model::cells[i]}; if (c.material_.size() > 1) { if (c.material_.size() != c.n_instances_) { @@ -228,10 +228,10 @@ prepare_distribcell() // unique distribcell array index. int distribcell_index = 0; std::vector target_univ_ids; - for (Universe* u : universes) { + for (Universe* u : model::universes) { for (auto cell_indx : u->cells_) { if (distribcells.find(cell_indx) != distribcells.end()) { - cells[cell_indx]->distribcell_index_ = distribcell_index; + model::cells[cell_indx]->distribcell_index_ = distribcell_index; target_univ_ids.push_back(u->id_); ++distribcell_index; } @@ -240,22 +240,22 @@ prepare_distribcell() // Allocate the cell and lattice offset tables. int n_maps = target_univ_ids.size(); - for (Cell* c : cells) { + for (Cell* c : model::cells) { if (c->type_ != FILL_MATERIAL) { c->offset_.resize(n_maps, C_NONE); } } - for (Lattice* lat : lattices) { + for (Lattice* lat : model::lattices) { lat->allocate_offset_table(n_maps); } // Fill the cell and lattice offset tables. for (int map = 0; map < target_univ_ids.size(); map++) { auto target_univ_id = target_univ_ids[map]; - for (Universe* univ : universes) { + for (Universe* univ : model::universes) { int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation. for (int32_t cell_indx : univ->cells_) { - Cell& c = *cells[cell_indx]; + Cell& c = *model::cells[cell_indx]; if (c.type_ == FILL_UNIVERSE) { c.offset_[map] = offset; @@ -263,7 +263,7 @@ prepare_distribcell() offset += count_universe_instances(search_univ, target_univ_id); } else if (c.type_ == FILL_LATTICE) { - Lattice& lat = *lattices[c.fill_]; + Lattice& lat = *model::lattices[c.fill_]; offset = lat.fill_offset_table(offset, target_univ_id, map); } } @@ -276,8 +276,8 @@ prepare_distribcell() void count_cell_instances(int32_t univ_indx) { - for (int32_t cell_indx : universes[univ_indx]->cells_) { - Cell& c = *cells[cell_indx]; + for (int32_t cell_indx : model::universes[univ_indx]->cells_) { + Cell& c = *model::cells[cell_indx]; ++c.n_instances_; if (c.type_ == FILL_UNIVERSE) { @@ -286,7 +286,7 @@ count_cell_instances(int32_t univ_indx) } else if (c.type_ == FILL_LATTICE) { // This cell contains a lattice. Recurse into the lattice universes. - Lattice& lat = *lattices[c.fill_]; + Lattice& lat = *model::lattices[c.fill_]; for (auto it = lat.begin(); it != lat.end(); ++it) { count_cell_instances(*it); } @@ -300,20 +300,20 @@ int count_universe_instances(int32_t search_univ, int32_t target_univ_id) { // If this is the target, it can't contain itself. - if (universes[search_univ]->id_ == target_univ_id) { + if (model::universes[search_univ]->id_ == target_univ_id) { return 1; } int count {0}; - for (int32_t cell_indx : universes[search_univ]->cells_) { - Cell& c = *cells[cell_indx]; + for (int32_t cell_indx : model::universes[search_univ]->cells_) { + Cell& c = *model::cells[cell_indx]; if (c.type_ == FILL_UNIVERSE) { int32_t next_univ = c.fill_; count += count_universe_instances(next_univ, target_univ_id); } else if (c.type_ == FILL_LATTICE) { - Lattice& lat = *lattices[c.fill_]; + Lattice& lat = *model::lattices[c.fill_]; for (auto it = lat.begin(); it != lat.end(); ++it) { int32_t next_univ = *it; count += count_universe_instances(next_univ, target_univ_id); @@ -338,7 +338,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, // write to the path and return. for (int32_t cell_indx : search_univ.cells_) { if ((cell_indx == target_cell) && (offset == target_offset)) { - Cell& c = *cells[cell_indx]; + Cell& c = *model::cells[cell_indx]; path << "c" << c.id_; return path.str(); } @@ -350,7 +350,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, std::vector::const_reverse_iterator cell_it {search_univ.cells_.crbegin()}; for (; cell_it != search_univ.cells_.crend(); ++cell_it) { - Cell& c = *cells[*cell_it]; + Cell& c = *model::cells[*cell_it]; // Material cells don't contain other cells so ignore them. if (c.type_ != FILL_MATERIAL) { @@ -358,7 +358,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, if (c.type_ == FILL_UNIVERSE) { temp_offset = offset + c.offset_[map]; } else { - Lattice& lat = *lattices[c.fill_]; + Lattice& lat = *model::lattices[c.fill_]; int32_t indx = lat.universes_.size()*map + lat.begin().indx_; temp_offset = offset + lat.offsets_[indx]; } @@ -370,18 +370,18 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, } // Add the cell to the path string. - Cell& c = *cells[*cell_it]; + Cell& c = *model::cells[*cell_it]; path << "c" << c.id_ << "->"; if (c.type_ == FILL_UNIVERSE) { // Recurse into the fill cell. offset += c.offset_[map]; path << distribcell_path_inner(target_cell, map, target_offset, - *universes[c.fill_], offset); + *model::universes[c.fill_], offset); return path.str(); } else { // Recurse into the lattice cell. - Lattice& lat = *lattices[c.fill_]; + Lattice& lat = *model::lattices[c.fill_]; path << "l" << lat.id_; for (ReverseLatticeIter it = lat.rbegin(); it != lat.rend(); ++it) { int32_t indx = lat.universes_.size()*map + it.indx_; @@ -390,7 +390,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, offset = temp_offset; path << "(" << lat.index_to_string(it.indx_) << ")->"; path << distribcell_path_inner(target_cell, map, target_offset, - *universes[*it], offset); + *model::universes[*it], offset); return path.str(); } } @@ -401,7 +401,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, std::string distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset) { - auto& root_univ = *universes[openmc_root_universe]; + auto& root_univ = *model::universes[model::root_universe]; return distribcell_path_inner(target_cell, map, target_offset, root_univ, 0); } @@ -412,13 +412,13 @@ maximum_levels(int32_t univ) { int levels_below {0}; - for (int32_t cell_indx : universes[univ]->cells_) { - Cell& c = *cells[cell_indx]; + for (int32_t cell_indx : model::universes[univ]->cells_) { + Cell& c = *model::cells[cell_indx]; if (c.type_ == FILL_UNIVERSE) { int32_t next_univ = c.fill_; levels_below = std::max(levels_below, maximum_levels(next_univ)); } else if (c.type_ == FILL_LATTICE) { - Lattice& lat = *lattices[c.fill_]; + Lattice& lat = *model::lattices[c.fill_]; for (auto it = lat.begin(); it != lat.end(); ++it) { int32_t next_univ = *it; levels_below = std::max(levels_below, maximum_levels(next_univ)); @@ -435,20 +435,20 @@ maximum_levels(int32_t univ) void free_memory_geometry_c() { - for (Cell* c : cells) {delete c;} - cells.clear(); - cell_map.clear(); - n_cells = 0; + for (Cell* c : model::cells) {delete c;} + model::cells.clear(); + model::cell_map.clear(); + model::n_cells = 0; - for (Universe* u : universes) {delete u;} - universes.clear(); - universe_map.clear(); + for (Universe* u : model::universes) {delete u;} + model::universes.clear(); + model::universe_map.clear(); - for (Lattice* lat : lattices) {delete lat;} - lattices.clear(); - lattice_map.clear(); + for (Lattice* lat : model::lattices) {delete lat;} + model::lattices.clear(); + model::lattice_map.clear(); - overlap_check_count.clear(); + model::overlap_check_count.clear(); } } // namespace openmc diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 9af1fcbaa..df2819382 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -143,7 +143,7 @@ module geometry_header end type Cell ! array index of the root universe - integer(C_INT), bind(C, name='openmc_root_universe') :: root_universe = -1 + integer(C_INT), bind(C) :: root_universe integer(C_INT32_T), bind(C) :: n_cells ! # of cells integer(C_INT32_T), bind(C) :: n_universes ! # of universes diff --git a/src/initialize.cpp b/src/initialize.cpp index e04ac2922..3a4a3a480 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -54,8 +54,8 @@ int openmc_init(int argc, char* argv[], const void* intracomm) if (err) return err; // Start total and initialization timer - time_total.start(); - time_initialize.start(); + simulation::time_total.start(); + simulation::time_initialize.start(); #ifdef _OPENMP // If OMP_SCHEDULE is not set, default to a static schedule @@ -79,7 +79,7 @@ int openmc_init(int argc, char* argv[], const void* intracomm) if (settings::particle_restart_run) settings::run_mode = RUN_MODE_PARTICLE; // Stop initialization timer - time_initialize.stop(); + simulation::time_initialize.stop(); return 0; } @@ -210,7 +210,7 @@ parse_command_line(int argc, char* argv[]) } omp_set_num_threads(simulation::n_threads); #else - if (openmc_master) + if (mpi::master) warning("Ignoring number of threads specified on command line."); #endif diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e2bf633b0..1da1bf9df 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -347,8 +347,8 @@ contains call read_surfaces(root % ptr) ! Allocate surfaces array - allocate(surfaces(n_surfaces)) - do i = 1, n_surfaces + allocate(surfaces(surfaces_size())) + do i = 1, size(surfaces) surfaces(i) % ptr = surface_pointer(i - 1); if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true. @@ -505,9 +505,9 @@ contains integer :: i ! Allocate surfaces array - allocate(surfaces(n_surfaces)) + allocate(surfaces(surfaces_size())) - do i = 1, n_surfaces + do i = 1, size(surfaces) surfaces(i) % ptr = surface_pointer(i - 1); end do diff --git a/src/lattice.cpp b/src/lattice.cpp index 14b679903..a8966f898 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -18,10 +18,13 @@ namespace openmc { // Global variables //============================================================================== -std::vector lattices; +namespace model { +std::vector lattices; std::unordered_map lattice_map; +} + //============================================================================== // Lattice implementation //============================================================================== @@ -65,8 +68,8 @@ Lattice::adjust_indices() // Adjust the indices for the universes array. for (LatticeIter it = begin(); it != end(); ++it) { int uid = *it; - auto search = universe_map.find(uid); - if (search != universe_map.end()) { + auto search = model::universe_map.find(uid); + if (search != model::universe_map.end()) { *it = search->second; } else { std::stringstream err_msg; @@ -78,8 +81,8 @@ Lattice::adjust_indices() // Adjust the index for the outer universe. if (outer_ != NO_OUTER_UNIVERSE) { - auto search = universe_map.find(outer_); - if (search != universe_map.end()) { + auto search = model::universe_map.find(outer_); + if (search != model::universe_map.end()) { outer_ = search->second; } else { std::stringstream err_msg; @@ -118,7 +121,7 @@ Lattice::to_hdf5(hid_t lattices_group) const } if (outer_ != NO_OUTER_UNIVERSE) { - int32_t outer_id = universes[outer_]->id_; + int32_t outer_id = model::universes[outer_]->id_; write_dataset(lat_group, "outer", outer_id); } else { write_dataset(lat_group, "outer", outer_); @@ -370,7 +373,7 @@ RectLattice::to_hdf5_inner(hid_t lat_group) const for (int j = 0; j < nx; j++) { int indx1 = nx*ny*m + nx*k + j; int indx2 = nx*ny*m + nx*(ny-k-1) + j; - out[indx2] = universes[universes_[indx1]]->id_; + out[indx2] = model::universes[universes_[indx1]]->id_; } } } @@ -387,7 +390,7 @@ RectLattice::to_hdf5_inner(hid_t lat_group) const for (int j = 0; j < nx; j++) { int indx1 = nx*k + j; int indx2 = nx*(ny-k-1) + j; - out[indx2] = universes[universes_[indx1]]->id_; + out[indx2] = model::universes[universes_[indx1]]->id_; } } @@ -847,7 +850,7 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const // This array position is never used; put a -1 to indicate this. out[indx] = -1; } else { - out[indx] = universes[universes_[indx]]->id_; + out[indx] = model::universes[universes_[indx]]->id_; } } } @@ -865,18 +868,18 @@ extern "C" void read_lattices(pugi::xml_node *node) { for (pugi::xml_node lat_node : node->children("lattice")) { - lattices.push_back(new RectLattice(lat_node)); + model::lattices.push_back(new RectLattice(lat_node)); } for (pugi::xml_node lat_node : node->children("hex_lattice")) { - lattices.push_back(new HexLattice(lat_node)); + model::lattices.push_back(new HexLattice(lat_node)); } // Fill the lattice map. - for (int i_lat = 0; i_lat < lattices.size(); i_lat++) { - int id = lattices[i_lat]->id_; - auto in_map = lattice_map.find(id); - if (in_map == lattice_map.end()) { - lattice_map[id] = i_lat; + for (int i_lat = 0; i_lat < model::lattices.size(); i_lat++) { + int id = model::lattices[i_lat]->id_; + auto in_map = model::lattice_map.find(id); + if (in_map == model::lattice_map.end()) { + model::lattice_map[id] = i_lat; } else { std::stringstream err_msg; err_msg << "Two or more lattices use the same unique ID: " << id; @@ -890,7 +893,7 @@ read_lattices(pugi::xml_node *node) //============================================================================== extern "C" { - Lattice* lattice_pointer(int lat_ind) {return lattices[lat_ind];} + Lattice* lattice_pointer(int lat_ind) {return model::lattices[lat_ind];} int32_t lattice_id(Lattice *lat) {return lat->id_;} } diff --git a/src/main.cpp b/src/main.cpp index f4e2d3f8e..b757c67ff 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,13 +1,15 @@ #ifdef OPENMC_MPI -#include "mpi.h" +#include #endif #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/error.h" +#include "openmc/message_passing.h" #include "openmc/settings.h" int main(int argc, char* argv[]) { + using namespace openmc; int err; // Initialize run -- when run with MPI, pass communicator @@ -21,30 +23,30 @@ int main(int argc, char* argv[]) { // This happens for the -h and -v flags return 0; } else if (err) { - openmc::fatal_error(openmc_err_msg); + fatal_error(openmc_err_msg); } // start problem based on mode - switch (openmc::settings::run_mode) { - case openmc::RUN_MODE_FIXEDSOURCE: - case openmc::RUN_MODE_EIGENVALUE: + switch (settings::run_mode) { + case RUN_MODE_FIXEDSOURCE: + case RUN_MODE_EIGENVALUE: err = openmc_run(); break; - case openmc::RUN_MODE_PLOTTING: + case RUN_MODE_PLOTTING: err = openmc_plot_geometry(); break; - case openmc::RUN_MODE_PARTICLE: - if (openmc_master) err = openmc_particle_restart(); + case RUN_MODE_PARTICLE: + if (mpi::master) err = openmc_particle_restart(); break; - case openmc::RUN_MODE_VOLUME: + case RUN_MODE_VOLUME: err = openmc_calculate_volumes(); break; } - if (err) openmc::fatal_error(openmc_err_msg); + if (err) fatal_error(openmc_err_msg); // Finalize and free up memory err = openmc_finalize(); - if (err) openmc::fatal_error(openmc_err_msg); + if (err) fatal_error(openmc_err_msg); // If MPI is in use and enabled, terminate it #ifdef OPENMC_MPI diff --git a/src/material.cpp b/src/material.cpp index ecc6eed2b..7f8cfc50b 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -13,9 +13,13 @@ namespace openmc { // Global variables //============================================================================== +namespace model { + std::vector materials; std::unordered_map material_map; +} // namespace model + //============================================================================== // Material implementation //============================================================================== @@ -46,16 +50,16 @@ read_materials(pugi::xml_node* node) { // Loop over XML material elements and populate the array. for (pugi::xml_node material_node : node->children("material")) { - materials.push_back(new Material(material_node)); + model::materials.push_back(new Material(material_node)); } - materials.shrink_to_fit(); + model::materials.shrink_to_fit(); // Populate the material map. - for (int i = 0; i < materials.size(); i++) { - int32_t mid = materials[i]->id_; - auto search = material_map.find(mid); - if (search == material_map.end()) { - material_map[mid] = i; + for (int i = 0; i < model::materials.size(); i++) { + int32_t mid = model::materials[i]->id_; + auto search = model::material_map.find(mid); + if (search == model::material_map.end()) { + model::material_map[mid] = i; } else { std::stringstream err_msg; err_msg << "Two or more materials use the same unique ID: " << mid; @@ -71,8 +75,8 @@ read_materials(pugi::xml_node* node) extern "C" int openmc_material_get_volume(int32_t index, double* volume) { - if (index >= 1 && index <= materials.size()) { - Material* m = materials[index - 1]; + if (index >= 1 && index <= model::materials.size()) { + Material* m = model::materials[index - 1]; if (m->volume_ >= 0.0) { *volume = m->volume_; return 0; @@ -91,8 +95,8 @@ openmc_material_get_volume(int32_t index, double* volume) extern "C" int openmc_material_set_volume(int32_t index, double volume) { - if (index >= 1 && index <= materials.size()) { - Material* m = materials[index - 1]; + if (index >= 1 && index <= model::materials.size()) { + Material* m = model::materials[index - 1]; if (volume >= 0.0) { m->volume_ = volume; return 0; @@ -111,7 +115,7 @@ openmc_material_set_volume(int32_t index, double volume) //============================================================================== extern "C" { - Material* material_pointer(int32_t indx) {return materials[indx];} + Material* material_pointer(int32_t indx) {return model::materials[indx];} int32_t material_id(Material* mat) {return mat->id_;} @@ -119,7 +123,7 @@ extern "C" { { mat->id_ = id; //TODO: off-by-one - material_map[id] = index - 1; + model::material_map[id] = index - 1; } bool material_fissionable(Material* mat) {return mat->fissionable;} @@ -131,17 +135,17 @@ extern "C" { void extend_materials_c(int32_t n) { - materials.reserve(materials.size() + n); + model::materials.reserve(model::materials.size() + n); for (int32_t i = 0; i < n; i++) { - materials.push_back(new Material()); + model::materials.push_back(new Material()); } } void free_memory_material_c() { - for (Material *mat : materials) {delete mat;} - materials.clear(); - material_map.clear(); + for (Material *mat : model::materials) {delete mat;} + model::materials.clear(); + model::material_map.clear(); } } diff --git a/src/mesh.cpp b/src/mesh.cpp index 2c3e508fe..6ec298a78 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -29,10 +29,13 @@ namespace openmc { // Global variables //============================================================================== -std::vector> meshes; +namespace model { +std::vector> meshes; std::unordered_map mesh_map; +} // namespace model + //============================================================================== // RegularMesh implementation //============================================================================== @@ -44,7 +47,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) id_ = std::stoi(get_node_value(node, "id")); // Check to make sure 'id' hasn't been used - if (mesh_map.find(id_) != mesh_map.end()) { + if (model::mesh_map.find(id_) != model::mesh_map.end()) { fatal_error("Two or more meshes use the same unique ID: " + std::to_string(id_)); } @@ -724,11 +727,11 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, extern "C" int openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end) { - if (index_start) *index_start = meshes.size(); + if (index_start) *index_start = model::meshes.size(); for (int i = 0; i < n; ++i) { - meshes.emplace_back(new RegularMesh{}); + model::meshes.emplace_back(new RegularMesh{}); } - if (index_end) *index_end = meshes.size() - 1; + if (index_end) *index_end = model::meshes.size() - 1; return 0; } @@ -737,8 +740,8 @@ openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end) extern "C" int openmc_get_mesh_index(int32_t id, int32_t* index) { - auto pair = mesh_map.find(id); - if (pair == mesh_map.end()) { + auto pair = model::mesh_map.find(id); + if (pair == model::mesh_map.end()) { set_errmsg("No mesh exists with ID=" + std::to_string(id) + "."); return OPENMC_E_INVALID_ID; } @@ -750,11 +753,11 @@ openmc_get_mesh_index(int32_t id, int32_t* index) extern "C" int openmc_mesh_get_id(int32_t index, int32_t* id) { - if (index < 0 || index >= meshes.size()) { + if (index < 0 || index >= model::meshes.size()) { set_errmsg("Index in meshes array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - *id = meshes[index]->id_; + *id = model::meshes[index]->id_; return 0; } @@ -762,12 +765,12 @@ openmc_mesh_get_id(int32_t index, int32_t* id) extern "C" int openmc_mesh_set_id(int32_t index, int32_t id) { - if (index < 0 || index >= meshes.size()) { + if (index < 0 || index >= model::meshes.size()) { set_errmsg("Index in meshes array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - meshes[index]->id_ = id; - mesh_map[id] = index; + model::meshes[index]->id_ = id; + model::mesh_map[id] = index; return 0; } @@ -775,12 +778,12 @@ openmc_mesh_set_id(int32_t index, int32_t id) extern "C" int openmc_mesh_get_dimension(int32_t index, int** dims, int* n) { - if (index < 0 || index >= meshes.size()) { + if (index < 0 || index >= model::meshes.size()) { set_errmsg("Index in meshes array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - *dims = meshes[index]->shape_.data(); - *n = meshes[index]->n_dimension_; + *dims = model::meshes[index]->shape_.data(); + *n = model::meshes[index]->n_dimension_; return 0; } @@ -788,14 +791,14 @@ openmc_mesh_get_dimension(int32_t index, int** dims, int* n) extern "C" int openmc_mesh_set_dimension(int32_t index, int n, const int* dims) { - if (index < 0 || index >= meshes.size()) { + if (index < 0 || index >= model::meshes.size()) { set_errmsg("Index in meshes array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } // Copy dimension std::vector shape = {static_cast(n)}; - auto& m = meshes[index]; + auto& m = model::meshes[index]; m->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape); m->n_dimension_ = m->shape_.size(); @@ -806,12 +809,12 @@ openmc_mesh_set_dimension(int32_t index, int n, const int* dims) extern "C" int openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n) { - if (index < 0 || index >= meshes.size()) { + if (index < 0 || index >= model::meshes.size()) { set_errmsg("Index in meshes array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - auto& m = meshes[index]; + auto& m = model::meshes[index]; if (m->lower_left_.dimension() == 0) { set_errmsg("Mesh parameters have not been set."); return OPENMC_E_ALLOCATE; @@ -829,12 +832,12 @@ extern "C" int openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, const double* width) { - if (index < 0 || index >= meshes.size()) { + if (index < 0 || index >= model::meshes.size()) { set_errmsg("Index in meshes array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - auto& m = meshes[index]; + auto& m = model::meshes[index]; std::vector shape = {static_cast(n)}; if (ll && ur) { m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape); @@ -864,10 +867,10 @@ void read_meshes(pugi::xml_node* root) { for (auto node : root->children("mesh")) { // Read mesh and add to vector - meshes.emplace_back(new RegularMesh{node}); + model::meshes.emplace_back(new RegularMesh{node}); // Map ID to position in vector - mesh_map[meshes.back()->id_] = meshes.size() - 1; + model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1; } } @@ -875,13 +878,13 @@ void meshes_to_hdf5(hid_t group) { // Write number of meshes hid_t meshes_group = create_group(group, "meshes"); - int32_t n_meshes = meshes.size(); + int32_t n_meshes = model::meshes.size(); write_attribute(meshes_group, "n_meshes", n_meshes); if (n_meshes > 0) { // Write IDs of meshes std::vector ids; - for (const auto& m : meshes) { + for (const auto& m : model::meshes) { m->to_hdf5(meshes_group); ids.push_back(m->id_); } @@ -896,9 +899,9 @@ void meshes_to_hdf5(hid_t group) //============================================================================== extern "C" { - int n_meshes() { return meshes.size(); } + int n_meshes() { return model::meshes.size(); } - RegularMesh* mesh_ptr(int i) { return meshes.at(i).get(); } + RegularMesh* mesh_ptr(int i) { return model::meshes.at(i).get(); } int32_t mesh_id(RegularMesh* m) { return m->id_; } @@ -936,8 +939,8 @@ extern "C" { void free_memory_mesh() { - meshes.clear(); - mesh_map.clear(); + model::meshes.clear(); + model::mesh_map.clear(); } } diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 8c14f324e..85107cbf9 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -22,10 +22,18 @@ namespace openmc { +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + // Storage for the MGXS data std::vector nuclides_MG; std::vector macro_xs; +} // namespace data + //============================================================================== // Mgxs base-class methods //============================================================================== diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 80aa66621..710d76a52 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -13,10 +13,14 @@ namespace openmc { // Global variable definitions //============================================================================== +namespace data { + std::vector energy_bins; std::vector energy_bin_avg; std::vector rev_energy_bins; +} // namesapce data + //============================================================================== // Mgxs data loading interface methods //============================================================================== @@ -44,7 +48,7 @@ add_mgxs_c(hid_t file_id, const char* name, int energy_groups, Mgxs mg(xs_grp, energy_groups, delayed_groups, temperature, tolerance, max_order, legendre_to_tabular, legendre_to_tabular_points, method); - nuclides_MG.push_back(mg); + data::nuclides_MG.push_back(mg); close_group(xs_grp); } @@ -55,7 +59,7 @@ query_fissionable_c(int n_nuclides, const int i_nuclides[]) { bool result = false; for (int n = 0; n < n_nuclides; n++) { - if (nuclides_MG[i_nuclides[n] - 1].fissionable) result = true; + if (data::nuclides_MG[i_nuclides[n] - 1].fissionable) result = true; } return result; } @@ -79,16 +83,16 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[], // material std::vector mgxs_ptr(n_nuclides); for (int n = 0; n < n_nuclides; n++) { - mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1]; + mgxs_ptr[n] = &data::nuclides_MG[i_nuclides[n] - 1]; } Mgxs macro(mat_name, temperature, mgxs_ptr, atom_densities_vec, tolerance, method); - macro_xs.emplace_back(macro); + data::macro_xs.emplace_back(macro); } else { // Preserve the ordering of materials by including a blank entry Mgxs macro; - macro_xs.emplace_back(macro); + data::macro_xs.emplace_back(macro); } } @@ -97,18 +101,18 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[], void read_mg_cross_sections_header_c(hid_t file_id) { ensure_exists(file_id, "energy_groups", true); - read_attribute(file_id, "energy_groups", num_energy_groups); + read_attribute(file_id, "energy_groups", data::num_energy_groups); ensure_exists(file_id, "group structure", true); - read_attribute(file_id, "group structure", rev_energy_bins); + read_attribute(file_id, "group structure", data::rev_energy_bins); // Reverse energy bins - std::copy(rev_energy_bins.crbegin(), rev_energy_bins.crend(), - std::back_inserter(energy_bins)); + std::copy(data::rev_energy_bins.crbegin(), data::rev_energy_bins.crend(), + std::back_inserter(data::energy_bins)); // Create average energies - for (int i = 0; i < energy_bins.size() - 1; ++i) { - energy_bin_avg.push_back(0.5*(energy_bins[i] + energy_bins[i+1])); + for (int i = 0; i < data::energy_bins.size() - 1; ++i) { + data::energy_bin_avg.push_back(0.5*(data::energy_bins[i] + data::energy_bins[i+1])); } // Add entries into libraries for MG data @@ -134,7 +138,7 @@ void calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs) { - macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs, + data::macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs, nu_fiss_xs); } @@ -159,7 +163,7 @@ get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg) } else { dg_c_p = dg; } - return nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p); + return data::nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p); } //============================================================================== @@ -183,7 +187,7 @@ get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg) } else { dg_c_p = dg; } - return macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p); + return data::macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p); } //============================================================================== @@ -192,7 +196,7 @@ void set_nuclide_angle_index_c(int index, const double uvw[3]) { // Update the values - nuclides_MG[index - 1].set_angle_index(uvw); + data::nuclides_MG[index - 1].set_angle_index(uvw); } //============================================================================== @@ -201,7 +205,7 @@ void set_macro_angle_index_c(int index, const double uvw[3]) { // Update the values - macro_xs[index - 1].set_angle_index(uvw); + data::macro_xs[index - 1].set_angle_index(uvw); } //============================================================================== @@ -210,7 +214,7 @@ void set_nuclide_temperature_index_c(int index, double sqrtkT) { // Update the values - nuclides_MG[index - 1].set_temperature_index(sqrtkT); + data::nuclides_MG[index - 1].set_temperature_index(sqrtkT); } //============================================================================== @@ -225,7 +229,7 @@ get_name_c(int index, int name_len, char* name) std::strcpy(name, str.c_str()); // Now get the data and copy to the C-string - str = nuclides_MG[index - 1].name; + str = data::nuclides_MG[index - 1].name; std::strcpy(name, str.c_str()); // Finally, remove the null terminator @@ -237,7 +241,7 @@ get_name_c(int index, int name_len, char* name) double get_awr_c(int index) { - return nuclides_MG[index - 1].awr; + return data::nuclides_MG[index - 1].awr; } } // namespace openmc diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 79e2289ae..18e7e5a7e 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -6,9 +6,13 @@ namespace openmc { // Global variables //============================================================================== +namespace data { + std::array energy_min {0.0, 0.0}; std::array energy_max {INFTY, INFTY}; +} // namespace data + //============================================================================== // Fortran compatibility functions //============================================================================== @@ -16,8 +20,8 @@ std::array energy_max {INFTY, INFTY}; extern "C" void set_particle_energy_bounds(int particle, double E_min, double E_max) { - energy_min[particle - 1] = E_min; - energy_max[particle - 1] = E_max; + data::energy_min[particle - 1] = E_min; + data::energy_max[particle - 1] = E_max; } } // namespace openmc diff --git a/src/output.cpp b/src/output.cpp index 79800e21a..1b197a29e 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -51,7 +51,7 @@ std::string time_stamp() << ":" << now->tm_min << ":" << now->tm_sec; return ts.str(); } - + //============================================================================== //=============================================================================== @@ -62,7 +62,7 @@ void print_plot() { header("PLOTTING SUMMARY", 5); - for (auto pl : plots) { + for (auto pl : model::plots) { // Plot id std::cout << "Plot ID: " << pl.id_ << "\n"; // Plot filename @@ -74,7 +74,7 @@ void print_plot() { if (PlotType::slice == pl.type_) { std::cout << "Plot Type: Slice" << "\n"; } else if (PlotType::voxel == pl.type_) { - std::cout << "Plot Type: Voxel" << "\n"; + std::cout << "Plot Type: Voxel" << "\n"; } // Plot parameters @@ -98,9 +98,9 @@ void print_plot() { if (PlotColorBy::cells == pl.color_by_) { std::cout << "Coloring: Cells" << "\n"; } else if (PlotColorBy::mats == pl.color_by_) { - std::cout << "Coloring: Materials" << "\n"; + std::cout << "Coloring: Materials" << "\n"; } - + if (PlotType::slice == pl.type_) { switch(pl.basis_) { case PlotBasis::xy: @@ -122,29 +122,29 @@ void print_plot() { } std::cout << "\n"; - + } } - + void print_overlap_check() { #ifdef OPENMC_MPI - std::vector temp(overlap_check_count); - int err = MPI_Reduce(temp.data(), overlap_check_count.data(), - overlap_check_count.size(), MPI_INT64_T, MPI_SUM, 0, + std::vector temp(model::overlap_check_count); + int err = MPI_Reduce(temp.data(), model::overlap_check_count.data(), + model::overlap_check_count.size(), MPI_INT64_T, MPI_SUM, 0, mpi::intracomm); #endif - if (openmc_master) { + if (mpi::master) { header("cell overlap check summary", 1); std::cout << " Cell ID No. Overlap Checks\n"; std::vector sparse_cell_ids; - for (int i = 0; i < n_cells; i++) { - std::cout << " " << std::setw(8) << cells[i]->id_ << std::setw(17) - << overlap_check_count[i] << "\n"; - if (overlap_check_count[i] < 10) { - sparse_cell_ids.push_back(cells[i]->id_); + for (int i = 0; i < model::cells.size(); i++) { + std::cout << " " << std::setw(8) << model::cells[i]->id_ << std::setw(17) + << model::overlap_check_count[i] << "\n"; + if (model::overlap_check_count[i] < 10) { + sparse_cell_ids.push_back(model::cells[i]->id_); } } diff --git a/src/particle.cpp b/src/particle.cpp index c8a5a1b56..f87ff3066 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -113,7 +113,7 @@ Particle::from_source(const Bank* src) } else { g = static_cast(src->E); last_g = static_cast(src->E); - E = energy_bin_avg[g - 1]; + E = data::energy_bin_avg[g - 1]; } last_E = E; } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 40e5ba1ae..a4795f82a 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -21,14 +21,13 @@ namespace openmc { void -collision_mg(Particle* p, const double* energy_bin_avg, - const MaterialMacroXS* material_xs) +collision_mg(Particle* p, const MaterialMacroXS* material_xs) { // Add to the collision counter for the particle p->n_collision++; // Sample the reaction type - sample_reaction(p, energy_bin_avg, material_xs); + sample_reaction(p, material_xs); // Display information about collision if ((settings::verbosity >= 10) || (simulation::trace)) { @@ -39,21 +38,20 @@ collision_mg(Particle* p, const double* energy_bin_avg, } void -sample_reaction(Particle* p, const double* energy_bin_avg, - const MaterialMacroXS* material_xs) +sample_reaction(Particle* p, const MaterialMacroXS* material_xs) { // Create fission bank sites. Note that while a fission reaction is sampled, // it never actually "happens", i.e. the weight of the particle does not // change when sampling fission sites. The following block handles all // absorption (including fission) - if (materials[p->material - 1]->fissionable) { + if (model::materials[p->material - 1]->fissionable) { if (settings::run_mode == RUN_MODE_EIGENVALUE) { Bank* result_bank; int64_t result_bank_size; // Get pointer to fission bank from Fortran side openmc_fission_bank(&result_bank, &result_bank_size); - create_fission_sites(p, result_bank, &n_bank, result_bank_size, + create_fission_sites(p, result_bank, &simulation::n_bank, result_bank_size, material_xs); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { @@ -72,7 +70,7 @@ sample_reaction(Particle* p, const double* energy_bin_avg, if (!p->alive) return; // Sample a scattering event to determine the energy of the exiting neutron - scatter(p, energy_bin_avg); + scatter(p); // Play Russian roulette if survival biasing is turned on if (settings::survival_biasing) { @@ -82,14 +80,14 @@ sample_reaction(Particle* p, const double* energy_bin_avg, } void -scatter(Particle* p, const double* energy_bin_avg) +scatter(Particle* p) { // Adjust indices for Fortran to C++ indexing // TODO: Remove when no longer needed int gin = p->last_g - 1; int gout = p->g - 1; int i_mat = p->material - 1; - macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt); + data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt); // Adjust return value for fortran indexing // TODO: Remove when no longer needed @@ -99,7 +97,7 @@ scatter(Particle* p, const double* energy_bin_avg) rotate_angle_c(p->coord[0].uvw, p->mu, nullptr); // Update energy value for downstream compatability (in tallying) - p->E = energy_bin_avg[gout]; + p->E = data::energy_bin_avg[gout]; // Set event component p->event = EVENT_SCATTER; @@ -180,7 +178,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, // the energy in the fission bank int dg; int gout; - macro_xs[p->material - 1].sample_fission_energy(p->g - 1, dg, gout); + data::macro_xs[p->material - 1].sample_fission_energy(p->g - 1, dg, gout); bank_array[i].E = static_cast(gout + 1); bank_array[i].delayed_group = dg + 1; diff --git a/src/plot.cpp b/src/plot.cpp index 9a8aba139..70dcd11b5 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -9,6 +9,7 @@ #include "openmc/geometry.h" #include "openmc/cell.h" #include "openmc/material.h" +#include "openmc/message_passing.h" #include "openmc/string_utils.h" #include "openmc/mesh.h" #include "openmc/output.h" @@ -19,19 +20,23 @@ namespace openmc { +//============================================================================== +// Constants +//============================================================================== + +const RGBColor WHITE {255, 255, 255}; +constexpr int PLOT_LEVEL_LOWEST {-1}; //!< lower bound on plot universe level + //============================================================================== // Global variables //============================================================================== -int PLOT_LEVEL_LOWEST = -1; - -std::unordered_map plot_map; - -int n_plots; +namespace model { std::vector plots; +std::unordered_map plot_map; -const RGBColor WHITE = {255, 255, 255}; +} // namespace model //============================================================================== // RUN_PLOT controls the logic for making one or many plots @@ -42,7 +47,7 @@ int openmc_plot_geometry() { int err; - for (auto pl : plots) { + for (auto pl : model::plots) { std::stringstream ss; ss << "Processing plot " << pl.id_ << ": " << pl.path_plot_ << "..."; @@ -63,11 +68,10 @@ int openmc_plot_geometry() void read_plots(pugi::xml_node* plots_node) { - n_plots = 0; for (auto node : plots_node->children("plot")) { Plot pl(node); - plots.push_back(pl); - plot_map[pl.id_] = n_plots++; + model::plots.push_back(pl); + model::plot_map[pl.id_] = model::plots.size() - 1; } } @@ -122,7 +126,7 @@ void create_ppm(Plot pl) p.initialize(); std::copy(xyz, xyz+3, p.coord[0].xyz); std::copy(dir, dir+3, p.coord[0].uvw); - p.coord[0].universe = openmc_root_universe; + p.coord[0].universe = model::root_universe; #pragma omp for for (int y = 0; y < height; y++) { @@ -155,7 +159,7 @@ Plot::set_id(pugi::xml_node plot_node) } // Check to make sure 'id' hasn't been used - if (plot_map.find(id_) != plot_map.end()) { + if (model::plot_map.find(id_) != model::plot_map.end()) { std::stringstream err_msg; err_msg << "Two or more plots use the same unique ID: " << id_; fatal_error(err_msg.str()); @@ -243,7 +247,7 @@ Plot::set_bg_color(pugi::xml_node plot_node) if (check_for_node(plot_node, "background")) { std::vector bg_rgb = get_node_array(plot_node, "background"); if (PlotType::voxel == type_) { - if (openmc_master) { + if (mpi::master) { std::stringstream err_msg; err_msg << "Background color ignored in voxel plot " << id_; @@ -361,30 +365,27 @@ Plot::set_default_colors(pugi::xml_node plot_node) } if ("cell" == pl_color_by) { color_by_ = PlotColorBy::cells; - colors_.resize(n_cells); - for (int i = 0; i < n_cells; i++) { - colors_[i] = random_color(); - } - + colors_.resize(model::cells.size()); } else if("material" == pl_color_by) { color_by_ = PlotColorBy::mats; - colors_.resize(n_materials); - for (int i = 0; i < materials.size(); i++) { - colors_[i] = random_color(); - } + colors_.resize(model::materials.size()); } else { std::stringstream err_msg; err_msg << "Unsupported plot color type '" << pl_color_by << "' in plot " << id_; fatal_error(err_msg); } + + for (auto& c : colors_) { + c = random_color(); + } } void Plot::set_user_colors(pugi::xml_node plot_node) { if (!plot_node.select_nodes("color").empty() && PlotType::voxel == type_) { - if (openmc_master) { + if (mpi::master) { std::stringstream err_msg; err_msg << "Color specifications ignored in voxel plot " << id_; @@ -412,8 +413,8 @@ Plot::set_user_colors(pugi::xml_node plot_node) } // Add RGB if (PlotColorBy::cells == color_by_) { - if (cell_map.find(col_id) != cell_map.end()) { - col_id = cell_map[col_id]; + if (model::cell_map.find(col_id) != model::cell_map.end()) { + col_id = model::cell_map[col_id]; colors_[col_id] = user_rgb; } else { std::stringstream err_msg; @@ -422,8 +423,8 @@ Plot::set_user_colors(pugi::xml_node plot_node) fatal_error(err_msg); } } else if (PlotColorBy::mats == color_by_) { - if (material_map.find(col_id) != material_map.end()) { - col_id = material_map[col_id]; + if (model::material_map.find(col_id) != model::material_map.end()) { + col_id = model::material_map[col_id]; colors_[col_id] = user_rgb; } else { std::stringstream err_msg; @@ -552,7 +553,7 @@ Plot::set_mask(pugi::xml_node plot_node) if (!mask_nodes.empty()) { if (PlotType::voxel == type_) { - if (openmc_master) { + if (mpi::master) { std::stringstream wrn_msg; wrn_msg << "Mask ignored in voxel plot " << id_; warning(wrn_msg); @@ -575,8 +576,8 @@ Plot::set_mask(pugi::xml_node plot_node) // in the cell and material arrays for (auto& col_id : iarray) { if (PlotColorBy::cells == color_by_) { - if (cell_map.find(col_id) != cell_map.end()) { - col_id = cell_map[col_id]; + if (model::cell_map.find(col_id) != model::cell_map.end()) { + col_id = model::cell_map[col_id]; } else { std::stringstream err_msg; @@ -585,8 +586,8 @@ Plot::set_mask(pugi::xml_node plot_node) fatal_error(err_msg); } } else if (PlotColorBy::mats == color_by_) { - if (material_map.find(col_id) != material_map.end()) { - col_id = material_map[col_id]; + if (model::material_map.find(col_id) != model::material_map.end()) { + col_id = model::material_map[col_id]; } else { std::stringstream err_msg; @@ -660,7 +661,7 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id) } else { if (PlotColorBy::mats == pl.color_by_) { // Assign color based on material - Cell* c = cells[p.coord[j].cell]; + Cell* c = model::cells[p.coord[j].cell]; if (c->type_ == FILL_UNIVERSE) { // If we stopped on a middle universe level, treat as if not found rgb = pl.not_found_; @@ -671,12 +672,12 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id) id = -1; } else { rgb = pl.colors_[p.material - 1]; - id = materials[p.material - 1]->id_; + id = model::materials[p.material - 1]->id_; } } else if (PlotColorBy::cells == pl.color_by_) { // Assign color based on cell rgb = pl.colors_[p.coord[j].cell]; - id = cells[p.coord[j].cell]->id_; + id = model::cells[p.coord[j].cell]->id_; } } // endif found_cell } @@ -759,7 +760,7 @@ void draw_mesh_lines(Plot pl, ImageData& data) width[1] = xyz_ur_plot[1] - xyz_ll_plot[1]; width[2] = xyz_ur_plot[2] - xyz_ll_plot[2]; - auto& m = meshes[pl.index_meshlines_mesh_]; + auto& m = model::meshes[pl.index_meshlines_mesh_]; int ijk_ll[3], ijk_ur[3]; bool in_mesh; @@ -851,7 +852,7 @@ void create_voxel(Plot pl) p.initialize(); std::copy(ll.begin(), ll.begin()+ll.size(), p.coord[0].xyz); std::copy(dir, dir+3, p.coord[0].uvw); - p.coord[0].universe = openmc_root_universe; + p.coord[0].universe = model::root_universe; // Open binary plot file for writing std::ofstream of; diff --git a/src/settings.cpp b/src/settings.cpp index 459d7115d..36657c6c6 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -17,6 +17,7 @@ #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/mesh.h" +#include "openmc/message_passing.h" #include "openmc/output.h" #include "openmc/random_lcg.h" #include "openmc/simulation.h" @@ -150,7 +151,7 @@ void get_run_parameters(pugi::xml_node node_base) // Preallocate space for keff and entropy by generation int m = settings::n_max_batches * settings::gen_per_batch; simulation::k_generation.reserve(m); - entropy.reserve(m); + simulation::entropy.reserve(m); // Get the trigger information for keff if (check_for_node(node_base, "keff_trigger")) { @@ -231,7 +232,7 @@ void read_settings_xml() // To this point, we haven't displayed any output since we didn't know what // the verbosity is. Now that we checked for it, show the title if necessary - if (openmc_master) { + if (mpi::master) { if (verbosity >= 2) title(); } write_message("Reading settings XML file...", 5); @@ -398,7 +399,7 @@ void read_settings_xml() omp_set_num_threads(simulation::n_threads); } #else - if (openmc_master) warning("OpenMC was not compiled with OpenMP support; " + if (mpi::master) warning("OpenMC was not compiled with OpenMP support; " "ignoring number of threads."); #endif } @@ -415,17 +416,17 @@ void read_settings_xml() // Get point to list of elements and make sure there is at least one for (pugi::xml_node node : root.children("source")) { - external_sources.emplace_back(node); + model::external_sources.emplace_back(node); } // If no source specified, default to isotropic point source at origin with Watt spectrum - if (external_sources.empty()) { + if (model::external_sources.empty()) { SourceDistribution source { UPtrSpace{new SpatialPoint({0.0, 0.0, 0.0})}, UPtrAngle{new Isotropic()}, UPtrDist{new Watt(0.988, 2.249e-6)} }; - external_sources.push_back(std::move(source)); + model::external_sources.push_back(std::move(source)); } // Check if we want to write out source @@ -506,12 +507,12 @@ void read_settings_xml() // Shannon Entropy mesh if (check_for_node(root, "entropy_mesh")) { int temp = std::stoi(get_node_value(root, "entropy_mesh")); - if (mesh_map.find(temp) == mesh_map.end()) { + if (model::mesh_map.find(temp) == model::mesh_map.end()) { std::stringstream msg; msg << "Mesh " << temp << " specified for Shannon entropy does not exist."; fatal_error(msg); } - index_entropy_mesh = mesh_map.at(temp); + index_entropy_mesh = model::mesh_map.at(temp); } else if (check_for_node(root, "entropy")) { warning("Specifying a Shannon entropy mesh via the element " @@ -520,18 +521,18 @@ void read_settings_xml() // Read entropy mesh from auto node_entropy = root.child("entropy"); - meshes.emplace_back(new RegularMesh{node_entropy}); + model::meshes.emplace_back(new RegularMesh{node_entropy}); // Set entropy mesh index - index_entropy_mesh = meshes.size() - 1; + index_entropy_mesh = model::meshes.size() - 1; // Assign ID and set mapping - meshes.back()->id_ = 10000; - mesh_map[10000] = index_entropy_mesh; + model::meshes.back()->id_ = 10000; + model::mesh_map[10000] = index_entropy_mesh; } if (index_entropy_mesh >= 0) { - auto& m = *meshes[index_entropy_mesh]; + auto& m = *model::meshes[index_entropy_mesh]; if (m.shape_.dimension() == 0) { // If the user did not specify how many mesh cells are to be used in // each direction, we automatically determine an appropriate number of @@ -551,13 +552,13 @@ void read_settings_xml() // Uniform fission source weighting mesh if (check_for_node(root, "ufs_mesh")) { auto temp = std::stoi(get_node_value(root, "ufs_mesh")); - if (mesh_map.find(temp) == mesh_map.end()) { + if (model::mesh_map.find(temp) == model::mesh_map.end()) { std::stringstream msg; msg << "Mesh " << temp << " specified for uniform fission site method " "does not exist."; fatal_error(msg); } - index_ufs_mesh = mesh_map.at(temp); + index_ufs_mesh = model::mesh_map.at(temp); } else if (check_for_node(root, "uniform_fs")) { warning("Specifying a UFS mesh via the element " @@ -566,14 +567,14 @@ void read_settings_xml() // Read entropy mesh from auto node_ufs = root.child("uniform_fs"); - meshes.emplace_back(new RegularMesh{node_ufs}); + model::meshes.emplace_back(new RegularMesh{node_ufs}); // Set entropy mesh index - index_ufs_mesh = meshes.size() - 1; + index_ufs_mesh = model::meshes.size() - 1; // Assign ID and set mapping - meshes.back()->id_ = 10001; - mesh_map[10001] = index_entropy_mesh; + model::meshes.back()->id_ = 10001; + model::mesh_map[10001] = index_entropy_mesh; } if (index_ufs_mesh >= 0) { diff --git a/src/simulation.cpp b/src/simulation.cpp index 84ebb0192..50bf12521 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -25,10 +25,10 @@ extern "C" bool cmfd_on; extern "C" void accumulate_tallies(); extern "C" void allocate_banks(); +extern "C" void allocate_tally_results(); extern "C" void check_triggers(); extern "C" void cmfd_init_batch(); extern "C" void cmfd_tally_init(); -extern "C" void configure_tallies(); extern "C" void execute_cmfd(); extern "C" void init_tally_routines(); extern "C" void join_bank_from_threads(); @@ -81,7 +81,7 @@ int openmc_simulation_init() // Allocate array for matching filter bins #pragma omp parallel { - filter_matches.resize(n_filters); + simulation::filter_matches.resize(model::tally_filters.size()); } // Set up tally procedure pointers @@ -92,7 +92,7 @@ int openmc_simulation_init() allocate_banks(); // Allocate tally results arrays if they're not allocated yet - configure_tallies(); + allocate_tally_results(); // Activate the CMFD tallies cmfd_tally_init(); @@ -104,8 +104,9 @@ int openmc_simulation_init() // will potentially populate k_generation and entropy) simulation::current_batch = 0; simulation::k_generation.clear(); - entropy.clear(); + simulation::entropy.clear(); simulation::need_depletion_rx = false; + openmc_reset(); // If this is a restart run, load the state point data and binary source // file @@ -139,12 +140,12 @@ int openmc_simulation_finalize() if (!simulation::initialized) return 0; // Stop active batch timer and start finalization timer - time_active.stop(); - time_finalize.start(); + simulation::time_active.stop(); + simulation::time_finalize.start(); #pragma omp parallel { - filter_matches.clear(); + simulation::filter_matches.clear(); } // Deallocate Fortran variables, set tallies to inactive @@ -166,8 +167,8 @@ int openmc_simulation_finalize() } // Stop timers and show timing statistics - time_finalize.stop(); - time_total.stop(); + simulation::time_finalize.stop(); + simulation::time_total.stop(); if (mpi::master) { if (settings::verbosity >= 6) print_runtime(); if (settings::verbosity >= 4) print_results(); @@ -200,7 +201,7 @@ int openmc_next_batch(int* status) initialize_generation(); // Start timer for transport - time_transport.start(); + simulation::time_transport.start(); // ==================================================================== // LOOP OVER PARTICLES @@ -218,7 +219,7 @@ int openmc_next_batch(int* status) } // Accumulate time for transport - time_transport.stop(); + simulation::time_transport.stop(); finalize_generation(); } @@ -305,10 +306,10 @@ void initialize_batch() // Manage active/inactive timers and activate tallies if necessary. if (first_inactive) { - time_inactive.start(); + simulation::time_inactive.start(); } else if (first_active) { - time_inactive.stop(); - time_active.start(); + simulation::time_inactive.stop(); + simulation::time_active.start(); for (int i = 1; i <= n_tallies; ++i) { // TODO: change one-based index openmc_tally_set_active(i, true); @@ -327,9 +328,9 @@ void initialize_batch() void finalize_batch() { // Reduce tallies onto master process and accumulate - time_tallies.start(); + simulation::time_tallies.start(); accumulate_tallies(); - time_tallies.stop(); + simulation::time_tallies.stop(); // Reset global tally results if (simulation::current_batch <= settings::n_inactive) { @@ -384,13 +385,13 @@ void initialize_generation() { if (settings::run_mode == RUN_MODE_EIGENVALUE) { // Reset number of fission bank sites - n_bank = 0; + simulation::n_bank = 0; // Count source sites if using uniform fission source weighting if (settings::ufs_on) ufs_count_sites(); // Store current value of tracklength k - keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE); + simulation::keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE); } } @@ -564,4 +565,10 @@ extern "C" void k_generation_clear() { simulation::k_generation.clear(); } extern "C" void k_generation_reserve(int i) { simulation::k_generation.reserve(i); } extern "C" int64_t work_index(int rank) { return simulation::work_index[rank]; } +// This function was moved here to get around a bug on macOS whereby an invalid +// pointer is returned for the threadprivate filter_matches +extern "C" FilterMatch* filter_match_pointer(int indx) { + return &simulation::filter_matches[indx]; +} + } // namespace openmc diff --git a/src/source.cpp b/src/source.cpp index 315ff2cdd..6d071c0c7 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -27,8 +27,12 @@ namespace openmc { // Global variables //============================================================================== +namespace model { + std::vector external_sources; +} + //============================================================================== // SourceDistribution implementation //============================================================================== @@ -167,9 +171,9 @@ Bank SourceDistribution::sample() const if (space_box) { if (space_box->only_fissionable()) { // Determine material - auto c = cells[cell_index - 1]; + auto c = model::cells[cell_index - 1]; int32_t mat_index = c->material_[instance]; - auto m = materials[mat_index]; + auto m = model::materials[mat_index]; if (mat_index == MATERIAL_VOID) { found = false; @@ -207,10 +211,10 @@ Bank SourceDistribution::sample() const auto energy_ptr = dynamic_cast(energy_.get()); if (energy_ptr) { auto energies = xt::adapt(energy_ptr->x()); - if (xt::any(energies > energy_max[p-1])) { + if (xt::any(energies > data::energy_max[p-1])) { fatal_error("Source energy above range of energies of at least " "one cross section table"); - } else if (xt::any(energies < energy_min[p-1])) { + } else if (xt::any(energies < data::energy_min[p-1])) { fatal_error("Source energy below range of energies of at least " "one cross section table"); } @@ -221,7 +225,7 @@ Bank SourceDistribution::sample() const site.E = energy_->sample(); // Resample if energy falls outside minimum or maximum particle energy - if (site.E < energy_max[p-1] && site.E > energy_min[p-1]) break; + if (site.E < data::energy_max[p-1] && site.E > data::energy_min[p-1]) break; } // Set delayed group @@ -299,28 +303,28 @@ Bank sample_external_source() // Determine total source strength double total_strength = 0.0; - for (auto& s : external_sources) + for (auto& s : model::external_sources) total_strength += s.strength(); // Sample from among multiple source distributions int i = 0; - if (external_sources.size() > 1) { + if (model::external_sources.size() > 1) { double xi = prn()*total_strength; double c = 0.0; - for (; i < external_sources.size(); ++i) { - c += external_sources[i].strength(); + for (; i < model::external_sources.size(); ++i) { + c += model::external_sources[i].strength(); if (xi < c) break; } } // Sample source site from i-th source distribution - Bank site {external_sources[i].sample()}; + Bank site {model::external_sources[i].sample()}; // If running in MG, convert site % E to group if (!settings::run_CE) { - site.E = lower_bound_index(rev_energy_bins.begin(), rev_energy_bins.end(), - site.E); - site.E = num_energy_groups - site.E; + site.E = lower_bound_index(data::rev_energy_bins.begin(), + data::rev_energy_bins.end(), site.E); + site.E = data::num_energy_groups - site.E; } // Set the random number generator back to the tracking stream. @@ -335,13 +339,13 @@ Bank sample_external_source() extern "C" void free_memory_source() { - external_sources.clear(); + model::external_sources.clear(); } extern "C" double total_source_strength() { double strength = 0.0; - for (const auto& s : external_sources) { + for (const auto& s : model::external_sources) { strength += s.strength(); } return strength; diff --git a/src/state_point.F90 b/src/state_point.F90 index 1286bb63d..2692856dc 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -55,7 +55,7 @@ contains !=============================================================================== function openmc_statepoint_write(filename, write_source) result(err) bind(C) - type(C_PTR), intent(in), optional :: filename + type(C_PTR), value :: filename logical(C_BOOL), intent(in), optional :: write_source integer(C_INT) :: err @@ -91,7 +91,7 @@ contains err = 0 ! Set the filename - if (present(filename)) then + if (c_associated(filename)) then call c_f_pointer(filename, string, [MAX_FILE_LEN]) filename_ = to_f_string(string) else diff --git a/src/state_point.cpp b/src/state_point.cpp index 620b5f71a..fdf700c44 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -113,7 +113,7 @@ write_source_bank(hid_t group_id, Bank* source_bank) #else - if (openmc_master) { + if (mpi::master) { // Create dataset big enough to hold all source sites hsize_t dims[] {static_cast(settings::n_particles)}; hid_t dspace = H5Screate_simple(1, dims, nullptr); @@ -320,11 +320,11 @@ void restart_set_keff() { if (simulation::restart_batch > settings::n_inactive) { for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) { - k_sum[0] += simulation::k_generation[i]; - k_sum[1] += std::pow(simulation::k_generation[i], 2); + simulation::k_sum[0] += simulation::k_generation[i]; + simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); } int n = settings::gen_per_batch*n_realizations; - simulation::keff = k_sum[0] / n; + simulation::keff = simulation::k_sum[0] / n; } else { simulation::keff = simulation::k_generation.back(); } diff --git a/src/summary.cpp b/src/summary.cpp index d2f083657..398f304fb 100644 --- a/src/summary.cpp +++ b/src/summary.cpp @@ -10,7 +10,7 @@ extern "C" void write_geometry(hid_t file_id) { auto geom_group = create_group(file_id, "geometry"); - + #ifdef DAGMC if (settings::dagmc) { write_attribute(geom_group, "dagmc", 1); @@ -18,25 +18,25 @@ write_geometry(hid_t file_id) { } #endif - write_attribute(geom_group, "n_cells", cells.size()); - write_attribute(geom_group, "n_surfaces", surfaces.size()); - write_attribute(geom_group, "n_universes", universes.size()); - write_attribute(geom_group, "n_lattices", lattices.size()); + write_attribute(geom_group, "n_cells", model::cells.size()); + write_attribute(geom_group, "n_surfaces", model::surfaces.size()); + write_attribute(geom_group, "n_universes", model::universes.size()); + write_attribute(geom_group, "n_lattices", model::lattices.size()); auto cells_group = create_group(geom_group, "cells"); - for (Cell* c : cells) c->to_hdf5(cells_group); + for (Cell* c : model::cells) c->to_hdf5(cells_group); close_group(cells_group); auto surfaces_group = create_group(geom_group, "surfaces"); - for (Surface* surf : surfaces) surf->to_hdf5(surfaces_group); + for (Surface* surf : model::surfaces) surf->to_hdf5(surfaces_group); close_group(surfaces_group); auto universes_group = create_group(geom_group, "universes"); - for (Universe* u : universes) u->to_hdf5(universes_group); + for (Universe* u : model::universes) u->to_hdf5(universes_group); close_group(universes_group); auto lattices_group = create_group(geom_group, "lattices"); - for (Lattice* lat : lattices) lat->to_hdf5(lattices_group); + for (Lattice* lat : model::lattices) lat->to_hdf5(lattices_group); close_group(lattices_group); close_group(geom_group); diff --git a/src/surface.cpp b/src/surface.cpp index 66217a58b..3ffbb3277 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -25,12 +25,13 @@ extern "C" const int BC_PERIODIC {3}; // Global variables //============================================================================== -int32_t n_surfaces; +namespace model { std::vector surfaces; - std::map surface_map; +} // namespace model + //============================================================================== // Helper functions for reading the "coeffs" node of an XML surface element //============================================================================== @@ -1063,13 +1064,14 @@ extern "C" void read_surfaces(pugi::xml_node* node) { // Count the number of surfaces. - for (pugi::xml_node surf_node: node->children("surface")) {n_surfaces++;} + int n_surfaces = 0; + for (pugi::xml_node surf_node : node->children("surface")) {n_surfaces++;} if (n_surfaces == 0) { fatal_error("No surfaces found in geometry.xml!"); } // Loop over XML surface elements and populate the array. - surfaces.reserve(n_surfaces); + model::surfaces.reserve(n_surfaces); { pugi::xml_node surf_node; int i_surf; @@ -1078,40 +1080,40 @@ read_surfaces(pugi::xml_node* node) std::string surf_type = get_node_value(surf_node, "type", true, true); if (surf_type == "x-plane") { - surfaces.push_back(new SurfaceXPlane(surf_node)); + model::surfaces.push_back(new SurfaceXPlane(surf_node)); } else if (surf_type == "y-plane") { - surfaces.push_back(new SurfaceYPlane(surf_node)); + model::surfaces.push_back(new SurfaceYPlane(surf_node)); } else if (surf_type == "z-plane") { - surfaces.push_back(new SurfaceZPlane(surf_node)); + model::surfaces.push_back(new SurfaceZPlane(surf_node)); } else if (surf_type == "plane") { - surfaces.push_back(new SurfacePlane(surf_node)); + model::surfaces.push_back(new SurfacePlane(surf_node)); } else if (surf_type == "x-cylinder") { - surfaces.push_back(new SurfaceXCylinder(surf_node)); + model::surfaces.push_back(new SurfaceXCylinder(surf_node)); } else if (surf_type == "y-cylinder") { - surfaces.push_back(new SurfaceYCylinder(surf_node)); + model::surfaces.push_back(new SurfaceYCylinder(surf_node)); } else if (surf_type == "z-cylinder") { - surfaces.push_back(new SurfaceZCylinder(surf_node)); + model::surfaces.push_back(new SurfaceZCylinder(surf_node)); } else if (surf_type == "sphere") { - surfaces.push_back(new SurfaceSphere(surf_node)); + model::surfaces.push_back(new SurfaceSphere(surf_node)); } else if (surf_type == "x-cone") { - surfaces.push_back(new SurfaceXCone(surf_node)); + model::surfaces.push_back(new SurfaceXCone(surf_node)); } else if (surf_type == "y-cone") { - surfaces.push_back(new SurfaceYCone(surf_node)); + model::surfaces.push_back(new SurfaceYCone(surf_node)); } else if (surf_type == "z-cone") { - surfaces.push_back(new SurfaceZCone(surf_node)); + model::surfaces.push_back(new SurfaceZCone(surf_node)); } else if (surf_type == "quadric") { - surfaces.push_back(new SurfaceQuadric(surf_node)); + model::surfaces.push_back(new SurfaceQuadric(surf_node)); } else { std::stringstream err_msg; @@ -1122,11 +1124,11 @@ read_surfaces(pugi::xml_node* node) } // Fill the surface map. - for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { - int id = surfaces[i_surf]->id_; - auto in_map = surface_map.find(id); - if (in_map == surface_map.end()) { - surface_map[id] = i_surf; + for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) { + int id = model::surfaces[i_surf]->id_; + auto in_map = model::surface_map.find(id); + if (in_map == model::surface_map.end()) { + model::surface_map[id] = i_surf; } else { std::stringstream err_msg; err_msg << "Two or more surfaces use the same unique ID: " << id; @@ -1138,10 +1140,10 @@ read_surfaces(pugi::xml_node* node) double xmin {INFTY}, xmax {-INFTY}, ymin {INFTY}, ymax {-INFTY}, zmin {INFTY}, zmax {-INFTY}; int i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax; - for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { - if (surfaces[i_surf]->bc_ == BC_PERIODIC) { + for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) { + if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) { // Downcast to the PeriodicSurface type. - Surface* surf_base = surfaces[i_surf]; + Surface* surf_base = model::surfaces[i_surf]; PeriodicSurface* surf = dynamic_cast(surf_base); // Make sure this surface inherits from PeriodicSurface. @@ -1183,10 +1185,10 @@ read_surfaces(pugi::xml_node* node) } // Set i_periodic for periodic BC surfaces. - for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { - if (surfaces[i_surf]->bc_ == BC_PERIODIC) { + for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) { + if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) { // Downcast to the PeriodicSurface type. - Surface* surf_base = surfaces[i_surf]; + Surface* surf_base = model::surfaces[i_surf]; PeriodicSurface* surf = dynamic_cast(surf_base); // Also try downcasting to the SurfacePlane type (which must be handled @@ -1216,7 +1218,7 @@ read_surfaces(pugi::xml_node* node) } } else { // Convert the surface id to an index. - surf->i_periodic_ = surface_map[surf->i_periodic_]; + surf->i_periodic_ = model::surface_map[surf->i_periodic_]; } } else { // This is a SurfacePlane. We won't try to find it's partner if the @@ -1228,12 +1230,12 @@ read_surfaces(pugi::xml_node* node) fatal_error(err_msg); } else { // Convert the surface id to an index. - surf->i_periodic_ = surface_map[surf->i_periodic_]; + surf->i_periodic_ = model::surface_map[surf->i_periodic_]; } } // Make sure the opposite surface is also periodic. - if (surfaces[surf->i_periodic_]->bc_ != BC_PERIODIC) { + if (model::surfaces[surf->i_periodic_]->bc_ != BC_PERIODIC) { std::stringstream err_msg; err_msg << "Could not find matching surface for periodic boundary " "condition on surface " << surf->id_; @@ -1248,7 +1250,7 @@ read_surfaces(pugi::xml_node* node) //============================================================================== extern "C" { - Surface* surface_pointer(int surf_ind) {return surfaces[surf_ind];} + Surface* surface_pointer(int surf_ind) {return model::surfaces[surf_ind];} int surface_id(Surface* surf) {return surf->id_;} @@ -1288,11 +1290,12 @@ extern "C" { void free_memory_surfaces_c() { - for (Surface* surf : surfaces) {delete surf;} - surfaces.clear(); - n_surfaces = 0; - surface_map.clear(); + for (Surface* surf : model::surfaces) {delete surf;} + model::surfaces.clear(); + model::surface_map.clear(); } + + int surfaces_size() { return model::surfaces.size(); } } } // namespace openmc diff --git a/src/surface_header.F90 b/src/surface_header.F90 index a2ef9c056..17f9515f3 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -78,10 +78,15 @@ module surface_header end type Surface - integer(C_INT32_T), bind(C) :: n_surfaces ! # of surfaces - type(Surface), allocatable, target :: surfaces(:) + interface + function surfaces_size() result(sz) bind(C) + import C_INT + integer(C_INT) :: sz + end function + end interface + contains pure function surface_id(this) result(id) diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index ef4ab345d..fd5e38f7e 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -26,6 +26,8 @@ #include "openmc/tallies/filter_universe.h" #include "openmc/tallies/filter_zernike.h" +// explicit template instantiation definition +template class std::vector; namespace openmc { @@ -33,8 +35,13 @@ namespace openmc { // Global variables //============================================================================== +namespace simulation { std::vector filter_matches; +} // namespace simulation + +namespace model { std::vector> tally_filters; +} // namespace model //============================================================================== // Non-member functions @@ -45,10 +52,10 @@ free_memory_tally_c() { #pragma omp parallel { - filter_matches.clear(); + simulation::filter_matches.clear(); } - tally_filters.clear(); + model::tally_filters.clear(); } //============================================================================== @@ -56,8 +63,7 @@ free_memory_tally_c() //============================================================================== extern "C" { - FilterMatch* filter_match_pointer(int indx) - {return &filter_matches[indx];} + // filter_match_point moved to simulation.cpp void filter_match_bins_push_back(FilterMatch* match, int val) @@ -96,53 +102,53 @@ extern "C" { { std::string type_ {type}; if (type_ == "azimuthal") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "cell") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "cellborn") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "cellfrom") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "distribcell") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "delayedgroup") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "energyfunction") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "energy") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "energyout") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "legendre") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "material") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "mesh") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "meshsurface") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "mu") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "particle") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "polar") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "surface") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "spatiallegendre") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "sphericalharmonics") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "universe") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "zernike") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else if (type_ == "zernikeradial") { - tally_filters.push_back(std::make_unique()); + model::tally_filters.push_back(std::make_unique()); } else { return nullptr; } - return tally_filters.back().get(); + return model::tally_filters.back().get(); } void filter_from_xml(Filter* filt, pugi::xml_node* node) diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp index 9e0a6721e..12dc38df1 100644 --- a/src/tallies/filter_cell.cpp +++ b/src/tallies/filter_cell.cpp @@ -21,8 +21,8 @@ CellFilter::initialize() { // Convert cell IDs to indices of the global array. for (auto& c : cells_) { - auto search = cell_map.find(c); - if (search != cell_map.end()) { + auto search = model::cell_map.find(c); + if (search != model::cell_map.end()) { c = search->second; } else { std::stringstream err_msg; @@ -56,7 +56,7 @@ CellFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); std::vector cell_ids; - for (auto c : cells_) cell_ids.push_back(cells[c]->id_); + for (auto c : cells_) cell_ids.push_back(model::cells[c]->id_); write_dataset(filter_group, "bins", cell_ids); } @@ -64,7 +64,7 @@ std::string CellFilter::text_label(int bin) const { //TODO: off-by-one - return "Cell " + std::to_string(cells[cells_[bin-1]]->id_); + return "Cell " + std::to_string(model::cells[cells_[bin-1]]->id_); } //============================================================================== diff --git a/src/tallies/filter_cellborn.cpp b/src/tallies/filter_cellborn.cpp index 4f86ad34c..b173b2e39 100644 --- a/src/tallies/filter_cellborn.cpp +++ b/src/tallies/filter_cellborn.cpp @@ -20,7 +20,7 @@ std::string CellbornFilter::text_label(int bin) const { //TODO: off-by-one - return "Birth Cell " + std::to_string(cells[cells_[bin-1]]->id_); + return "Birth Cell " + std::to_string(model::cells[cells_[bin-1]]->id_); } } // namespace openmc diff --git a/src/tallies/filter_cellfrom.cpp b/src/tallies/filter_cellfrom.cpp index b4b1a57ef..1262ebecd 100644 --- a/src/tallies/filter_cellfrom.cpp +++ b/src/tallies/filter_cellfrom.cpp @@ -22,7 +22,7 @@ std::string CellFromFilter::text_label(int bin) const { //TODO: off-by-one - return "Cell from " + std::to_string(cells[cells_[bin-1]]->id_); + return "Cell from " + std::to_string(model::cells[cells_[bin-1]]->id_); } } // namespace openmc diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp index 39459672f..ee4f08dd4 100644 --- a/src/tallies/filter_distribcell.cpp +++ b/src/tallies/filter_distribcell.cpp @@ -22,10 +22,10 @@ void DistribcellFilter::initialize() { // Convert the cell ID to an index of the global array. - auto search = cell_map.find(cell_); - if (search != cell_map.end()) { + auto search = model::cell_map.find(cell_); + if (search != model::cell_map.end()) { cell_ = search->second; - n_bins_ = cells[cell_]->n_instances_; + n_bins_ = model::cells[cell_]->n_instances_; } else { std::stringstream err_msg; err_msg << "Could not find cell " << cell_ @@ -39,13 +39,13 @@ DistribcellFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { int offset = 0; - auto distribcell_index = cells[cell_]->distribcell_index_; + auto distribcell_index = model::cells[cell_]->distribcell_index_; for (int i = 0; i < p->n_coord; i++) { - auto& c {*cells[p->coord[i].cell]}; + auto& c {*model::cells[p->coord[i].cell]}; if (c.type_ == FILL_UNIVERSE) { offset += c.offset_[distribcell_index]; } else if (c.type_ == FILL_LATTICE) { - auto& lat {*lattices[p->coord[i+1].lattice-1]}; + auto& lat {*model::lattices[p->coord[i+1].lattice-1]}; int i_xyz[3] {p->coord[i+1].lattice_x, p->coord[i+1].lattice_y, p->coord[i+1].lattice_z}; @@ -66,13 +66,13 @@ void DistribcellFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); - write_dataset(filter_group, "bins", cells[cell_]->id_); + write_dataset(filter_group, "bins", model::cells[cell_]->id_); } std::string DistribcellFilter::text_label(int bin) const { - auto map = cells[cell_]->distribcell_index_; + auto map = model::cells[cell_]->distribcell_index_; //TODO: off-by-one auto path = distribcell_path(cell_, map, bin-1); return "Distributed Cell " + path; diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp index 4c983b267..65339b519 100644 --- a/src/tallies/filter_energy.cpp +++ b/src/tallies/filter_energy.cpp @@ -25,10 +25,10 @@ EnergyFilter::from_xml(pugi::xml_node node) // (after flipping for the different ordering of the library and tallying // systems). if (!settings::run_CE) { - if (n_bins_ == num_energy_groups) { + if (n_bins_ == data::num_energy_groups) { matches_transport_groups_ = true; for (auto i = 0; i < n_bins_ + 1; i++) { - if (rev_energy_bins[i] != bins_[i]) { + if (data::rev_energy_bins[i] != bins_[i]) { matches_transport_groups_ = false; break; } @@ -44,10 +44,10 @@ const if (p->g != F90_NONE && matches_transport_groups_) { if (estimator == ESTIMATOR_TRACKLENGTH) { //TODO: off-by-one - match.bins_.push_back(num_energy_groups - p->g + 1); + match.bins_.push_back(data::num_energy_groups - p->g + 1); } else { //TODO: off-by-one - match.bins_.push_back(num_energy_groups - p->last_g + 1); + match.bins_.push_back(data::num_energy_groups - p->last_g + 1); } match.weights_.push_back(1.0); @@ -90,7 +90,7 @@ EnergyoutFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { if (p->g != F90_NONE && matches_transport_groups_) { - match.bins_.push_back(num_energy_groups - p->g + 1); + match.bins_.push_back(data::num_energy_groups - p->g + 1); match.weights_.push_back(1.0); } else { diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp index d36175e9b..be462b915 100644 --- a/src/tallies/filter_material.cpp +++ b/src/tallies/filter_material.cpp @@ -21,8 +21,8 @@ MaterialFilter::initialize() { // Convert material IDs to indices of the global array. for (auto& m : materials_) { - auto search = material_map.find(m); - if (search != material_map.end()) { + auto search = model::material_map.find(m); + if (search != model::material_map.end()) { m = search->second; } else { std::stringstream err_msg; @@ -55,7 +55,7 @@ MaterialFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); std::vector material_ids; - for (auto c : materials_) material_ids.push_back(materials[c]->id_); + for (auto c : materials_) material_ids.push_back(model::materials[c]->id_); write_dataset(filter_group, "bins", material_ids); } @@ -63,7 +63,7 @@ std::string MaterialFilter::text_label(int bin) const { //TODO: off-by-one - return "Material " + std::to_string(materials[materials_[bin-1]]->id_); + return "Material " + std::to_string(model::materials[materials_[bin-1]]->id_); } //============================================================================== diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index 5ff54bcf5..7f2b0d60c 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -20,8 +20,8 @@ MeshFilter::from_xml(pugi::xml_node node) } auto id = bins_[0]; - auto search = mesh_map.find(id); - if (search != mesh_map.end()) { + auto search = model::mesh_map.find(id); + if (search != model::mesh_map.end()) { set_mesh(search->second); } else{ std::stringstream err_msg; @@ -35,14 +35,14 @@ MeshFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { if (estimator != ESTIMATOR_TRACKLENGTH) { - auto bin = meshes[mesh_]->get_bin(p->coord[0].xyz); + auto bin = model::meshes[mesh_]->get_bin(p->coord[0].xyz); if (bin >= 0) { //TODO: off-by-one match.bins_.push_back(bin); match.weights_.push_back(1.0); } } else { - meshes[mesh_]->bins_crossed(p, match.bins_, match.weights_); + model::meshes[mesh_]->bins_crossed(p, match.bins_, match.weights_); } } @@ -50,13 +50,13 @@ void MeshFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); - write_dataset(filter_group, "bins", meshes[mesh_]->id_); + write_dataset(filter_group, "bins", model::meshes[mesh_]->id_); } std::string MeshFilter::text_label(int bin) const { - auto& mesh = *meshes[mesh_]; + auto& mesh = *model::meshes[mesh_]; int n_dim = mesh.n_dimension_; int ijk[n_dim]; @@ -76,7 +76,7 @@ MeshFilter::set_mesh(int32_t mesh) { mesh_ = mesh; n_bins_ = 1; - for (auto dim : meshes[mesh_]->shape_) n_bins_ *= dim; + for (auto dim : model::meshes[mesh_]->shape_) n_bins_ *= dim; } //============================================================================== @@ -123,7 +123,7 @@ openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh) } // Check the mesh index. - if (index_mesh < 0 || index_mesh >= meshes.size()) { + if (index_mesh < 0 || index_mesh >= model::meshes.size()) { set_errmsg("Index in 'meshes' array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp index af9bc5321..a8854c7ec 100644 --- a/src/tallies/filter_meshsurface.cpp +++ b/src/tallies/filter_meshsurface.cpp @@ -11,14 +11,14 @@ void MeshSurfaceFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - meshes[mesh_]->surface_bins_crossed(p, match.bins_); + model::meshes[mesh_]->surface_bins_crossed(p, match.bins_); for (auto b : match.bins_) match.weights_.push_back(1.0); } std::string MeshSurfaceFilter::text_label(int bin) const { - auto& mesh = *meshes[mesh_]; + auto& mesh = *model::meshes[mesh_]; int n_dim = mesh.n_dimension_; // Get flattend mesh index and surface index. @@ -76,8 +76,8 @@ void MeshSurfaceFilter::set_mesh(int32_t mesh) { mesh_ = mesh; - n_bins_ = 4 * meshes[mesh_]->n_dimension_; - for (auto dim : meshes[mesh_]->shape_) n_bins_ *= dim; + n_bins_ = 4 * model::meshes[mesh_]->n_dimension_; + for (auto dim : model::meshes[mesh_]->shape_) n_bins_ *= dim; } //============================================================================== diff --git a/src/tallies/filter_surface.cpp b/src/tallies/filter_surface.cpp index 5458d2051..5b175ca70 100644 --- a/src/tallies/filter_surface.cpp +++ b/src/tallies/filter_surface.cpp @@ -20,8 +20,8 @@ SurfaceFilter::initialize() { // Convert surface IDs to indices of the global array. for (auto& s : surfaces_) { - auto search = surface_map.find(s); - if (search != surface_map.end()) { + auto search = model::surface_map.find(s); + if (search != model::surface_map.end()) { s = search->second; } else { std::stringstream err_msg; @@ -58,7 +58,7 @@ SurfaceFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); std::vector surface_ids; - for (auto c : surfaces_) surface_ids.push_back(surfaces[c]->id_); + for (auto c : surfaces_) surface_ids.push_back(model::surfaces[c]->id_); write_dataset(filter_group, "bins", surface_ids); } @@ -66,7 +66,7 @@ std::string SurfaceFilter::text_label(int bin) const { //TODO: off-by-one - return "Surface " + std::to_string(surfaces[surfaces_[bin-1]]->id_); + return "Surface " + std::to_string(model::surfaces[surfaces_[bin-1]]->id_); } } // namespace openmc diff --git a/src/tallies/filter_universe.cpp b/src/tallies/filter_universe.cpp index 05d8c3e0d..ecd12b237 100644 --- a/src/tallies/filter_universe.cpp +++ b/src/tallies/filter_universe.cpp @@ -20,8 +20,8 @@ UniverseFilter::initialize() { // Convert universe IDs to indices of the global array. for (auto& u : universes_) { - auto search = universe_map.find(u); - if (search != universe_map.end()) { + auto search = model::universe_map.find(u); + if (search != model::universe_map.end()) { u = search->second; } else { std::stringstream err_msg; @@ -56,7 +56,7 @@ UniverseFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); std::vector universe_ids; - for (auto u : universes_) universe_ids.push_back(universes[u]->id_); + for (auto u : universes_) universe_ids.push_back(model::universes[u]->id_); write_dataset(filter_group, "bins", universe_ids); } @@ -64,7 +64,7 @@ std::string UniverseFilter::text_label(int bin) const { //TODO: off-by-one - return "Universe " + std::to_string(universes[universes_[bin-1]]->id_); + return "Universe " + std::to_string(model::universes[universes_[bin-1]]->id_); } } // namespace openmc diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 475c6a60c..b82837ea1 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -3792,11 +3792,6 @@ contains ! Accumulate on master only unless run is not reduced then do it on all if (master .or. (.not. reduce_tallies)) then - ! Accumulate results for each tally - do i = 1, active_tallies % size() - call tallies(active_tallies % data(i)) % obj % accumulate() - end do - if (run_mode == MODE_EIGENVALUE) then if (current_batch > n_inactive) then ! Accumulate products of different estimators of k @@ -3820,6 +3815,11 @@ contains end do end if + ! Accumulate results for each tally + do i = 1, active_tallies % size() + call tallies(active_tallies % data(i)) % obj % accumulate() + end do + end subroutine accumulate_tallies !=============================================================================== diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 0b80255d5..f2339d85d 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -6,7 +6,7 @@ module tally_header use error use dict_header, only: DictIntInt use hdf5_interface, only: HID_T, HSIZE_T - use message_passing, only: n_procs + use message_passing, only: n_procs, master use nuclide_header, only: nuclide_dict use settings, only: reduce_tallies, run_mode use stl_vector, only: VectorInt @@ -17,7 +17,7 @@ module tally_header implicit none private - public :: configure_tallies + public :: allocate_tally_results public :: free_memory_tally public :: openmc_extend_tallies public :: openmc_get_tally_index @@ -182,25 +182,27 @@ contains this % n_realizations = this % n_realizations + n_procs end if - ! Calculate total source strength for normalization - if (run_mode == MODE_FIXEDSOURCE) then - total_source = total_source_strength() - else - total_source = ONE - end if + if (master .or. (.not. reduce_tallies)) then + ! Calculate total source strength for normalization + if (run_mode == MODE_FIXEDSOURCE) then + total_source = total_source_strength() + else + total_source = ONE + end if - ! Accumulate each result - do j = 1, size(this % results, 3) - do i = 1, size(this % results, 2) - val = this % results(RESULT_VALUE, i, j)/total_weight * total_source - this % results(RESULT_VALUE, i, j) = ZERO + ! Accumulate each result + do j = 1, size(this % results, 3) + do i = 1, size(this % results, 2) + val = this % results(RESULT_VALUE, i, j)/total_weight * total_source + this % results(RESULT_VALUE, i, j) = ZERO - this % results(RESULT_SUM, i, j) = & - this % results(RESULT_SUM, i, j) + val - this % results(RESULT_SUM_SQ, i, j) = & - this % results(RESULT_SUM_SQ, i, j) + val*val + this % results(RESULT_SUM, i, j) = & + this % results(RESULT_SUM, i, j) + val + this % results(RESULT_SUM_SQ, i, j) = & + this % results(RESULT_SUM_SQ, i, j) + val*val + end do end do - end do + end if end subroutine tally_accumulate subroutine tally_write_results_hdf5(this, group_id) @@ -275,9 +277,6 @@ contains allocate(this % results(3, this % total_score_bins, this % n_filter_bins)) end if - ! Initialize results array to zero - this % results(:,:,:) = ZERO - end subroutine tally_allocate_results @@ -392,21 +391,21 @@ contains ! tallies.xml file. !=============================================================================== - subroutine configure_tallies() bind(C) + subroutine allocate_tally_results() bind(C) integer :: i - ! Allocate and initialize global tallies + ! Allocate global tallies if (.not. allocated(global_tallies)) then allocate(global_tallies(3, N_GLOBAL_TALLIES)) end if - global_tallies(:,:) = ZERO + ! Allocate results arrays for tallies do i = 1, n_tallies call tallies(i) % obj % allocate_results() end do - end subroutine configure_tallies + end subroutine allocate_tally_results !=============================================================================== ! FREE_MEMORY_TALLY deallocates global arrays defined in this module diff --git a/src/timer.cpp b/src/timer.cpp index f1c211446..7dab802b0 100644 --- a/src/timer.cpp +++ b/src/timer.cpp @@ -2,6 +2,25 @@ namespace openmc { +//============================================================================== +// Global variables +//============================================================================== + +namespace simulation { + +Timer time_active; +Timer time_bank; +Timer time_bank_sample; +Timer time_bank_sendrecv; +Timer time_finalize; +Timer time_inactive; +Timer time_initialize; +Timer time_tallies; +Timer time_total; +Timer time_transport; + +} // namespace simulation + //============================================================================== // Timer implementation //============================================================================== @@ -34,35 +53,20 @@ double Timer::elapsed() } } -//============================================================================== -// Global variables -//============================================================================== - -Timer time_active; -Timer time_bank; -Timer time_bank_sample; -Timer time_bank_sendrecv; -Timer time_finalize; -Timer time_inactive; -Timer time_initialize; -Timer time_tallies; -Timer time_total; -Timer time_transport; - //============================================================================== // Fortran compatibility //============================================================================== -extern "C" double time_active_elapsed() { return time_active.elapsed(); } -extern "C" double time_bank_elapsed() { return time_bank.elapsed(); } -extern "C" double time_bank_sample_elapsed() { return time_bank_sample.elapsed(); } -extern "C" double time_bank_sendrecv_elapsed() { return time_bank_sendrecv.elapsed(); } -extern "C" double time_finalize_elapsed() { return time_finalize.elapsed(); } -extern "C" double time_inactive_elapsed() { return time_inactive.elapsed(); } -extern "C" double time_initialize_elapsed() { return time_initialize.elapsed(); } -extern "C" double time_tallies_elapsed() { return time_tallies.elapsed(); } -extern "C" double time_total_elapsed() { return time_total.elapsed(); } -extern "C" double time_transport_elapsed() { return time_transport.elapsed(); } +extern "C" double time_active_elapsed() { return simulation::time_active.elapsed(); } +extern "C" double time_bank_elapsed() { return simulation::time_bank.elapsed(); } +extern "C" double time_bank_sample_elapsed() { return simulation::time_bank_sample.elapsed(); } +extern "C" double time_bank_sendrecv_elapsed() { return simulation::time_bank_sendrecv.elapsed(); } +extern "C" double time_finalize_elapsed() { return simulation::time_finalize.elapsed(); } +extern "C" double time_inactive_elapsed() { return simulation::time_inactive.elapsed(); } +extern "C" double time_initialize_elapsed() { return simulation::time_initialize.elapsed(); } +extern "C" double time_tallies_elapsed() { return simulation::time_tallies.elapsed(); } +extern "C" double time_total_elapsed() { return simulation::time_total.elapsed(); } +extern "C" double time_transport_elapsed() { return simulation::time_transport.elapsed(); } //============================================================================== // Non-member functions @@ -70,16 +74,16 @@ extern "C" double time_transport_elapsed() { return time_transport.elapsed(); } void reset_timers() { - time_active.reset(); - time_bank.reset(); - time_bank_sample.reset(); - time_bank_sendrecv.reset(); - time_finalize.reset(); - time_inactive.reset(); - time_initialize.reset(); - time_tallies.reset(); - time_total.reset(); - time_transport.reset(); + simulation::time_active.reset(); + simulation::time_bank.reset(); + simulation::time_bank_sample.reset(); + simulation::time_bank_sendrecv.reset(); + simulation::time_finalize.reset(); + simulation::time_inactive.reset(); + simulation::time_initialize.reset(); + simulation::time_tallies.reset(); + simulation::time_total.reset(); + simulation::time_transport.reset(); } } // namespace openmc diff --git a/src/tracking.F90 b/src/tracking.F90 index 038c7fdb5..e322f31f2 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -33,10 +33,9 @@ module tracking implicit none interface - subroutine collision_mg(p, energy_bin_avg, material_xs) bind(C) + subroutine collision_mg(p, material_xs) bind(C) import Particle, C_DOUBLE, MaterialMacroXS type(Particle), intent(inout) :: p - real(C_DOUBLE), intent(in) :: energy_bin_avg(*) type(MaterialMacroXS), intent(in) :: material_xs end subroutine collision_mg @@ -235,7 +234,7 @@ contains if (run_CE) then call collision(p) else - call collision_mg(p, energy_bin_avg, material_xs) + call collision_mg(p, material_xs) end if ! Score collision estimator tallies -- this is done after a collision