diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 3f7029b72..44309faa5 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -51,9 +51,11 @@ extern std::vector> meshes; } // namespace model -//============================================================================== -//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes -//============================================================================== +#ifdef LIBMESH +namespace settings { +extern std::unique_ptr LMI; +} +#endif class Mesh { @@ -302,7 +304,6 @@ public: //! Write the unstructured mesh to file // - // //! \param[in] filename Name of the file to write virtual void write(std::string filename) const = 0; @@ -323,16 +324,10 @@ public: std::string bin_label(int bin) const override; - void to_hdf5(hid_t group) const override; + void surface_bins_crossed(const Particle& p, + std::vector& bins) const override; - //! Count number of bank sites in each mesh bin / energy bin - // - //! \param[in] bank Array of bank sites - //! \param[out] Whether any bank sites are outside the mesh - //! \return Array indicating number of sites in each mesh/energy bin - xt::xtensor count_sites(const Particle::Bank* bank, - int64_t length, - bool* outside) const; + void to_hdf5(hid_t group) const override; // Data members bool output_ {true}; @@ -341,25 +336,19 @@ public: #ifdef DAGMC -class MOABUnstructuredMesh : public UnstructuredMesh { +class MOABMesh : public UnstructuredMesh { public: - MOABUnstructuredMesh() = default; - MOABUnstructuredMesh(pugi::xml_node); - MOABUnstructuredMesh(const std::string& filename); - ~MOABUnstructuredMesh() = default; + MOABMesh() = default; + MOABMesh(pugi::xml_node); + MOABMesh(const std::string& filename); + ~MOABMesh() = default; void bins_crossed(const Particle& p, std::vector& bins, std::vector& lengths) const override; - //! Determine which surface bins were crossed by a particle. - // - //! \param[in] p Particle to check - //! \param[out] bins Surface bins that were crossed - void surface_bins_crossed(const Particle& p, std::vector& bins) const; - int get_bin(Position r) const; int n_bins() const override; @@ -502,8 +491,6 @@ private: #ifdef LIBMESH class LibMesh : public UnstructuredMesh { - typedef std::vector> UnstructuredMeshHits; - public: // Constructors LibMesh(pugi::xml_node node); @@ -514,16 +501,14 @@ public: std::vector& bins, std::vector& lengths) const override; - void surface_bins_crossed(const Particle& p, std::vector& bins) const override; - int get_bin(Position r) const override; int n_bins() const override; int n_surface_bins() const override; - std::pair, std::vector> plot(Position plot_ll, - Position plot_ur) const override; + std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const override; void add_score(const std::string& var_name) override; @@ -552,12 +537,12 @@ private: int get_bin_from_element(const libMesh::Elem* elem) const; // Data members - std::unique_ptr m_; //!< pointer to the libmesh mesh instance + std::unique_ptr m_; //!< pointer to the libMesh mesh instance std::vector> point_locators_; //!< pointers to locators for each thread - std::unique_ptr equation_systems_; //!< pointer to the equation systems of the mesh (for result output) - std::map variable_map_; //!< mapping of variable names (scores) to their numbers on the mesh - libMesh::BoundingBox bbox_; //!< bounding box of the mesh + std::unique_ptr equation_systems_; //!< pointer to the equation systems of the mesh std::string eq_system_name_; //!< name of the equation system holding OpenMC results + std::map variable_map_; //!< mapping of variable names (tally scores) to libMesh variable numbers + libMesh::BoundingBox bbox_; //!< bounding box of the mesh libMesh::Elem* first_element_; //!< pointer to the first element in the mesh (maybe should be a key?) }; #endif diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 69231321d..650940c77 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -27,10 +27,6 @@ namespace openmc { namespace settings { -#ifdef LIBMESH -extern std::unique_ptr LMI; -extern const libMesh::Parallel::Communicator* libmesh_comm; -#endif // Boolean flags extern bool assume_separate; //!< assume tallies are spatially separate? diff --git a/src/initialize.cpp b/src/initialize.cpp index ca392073d..515442f80 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -75,9 +75,7 @@ if (!settings::LMI && !libMesh::initialized()) #else settings::LMI = std::make_unique(argc, argv, n_threads); #endif - settings::libmesh_comm = &(settings::LMI->comm()); } - #endif // Start total and initialization timer diff --git a/src/mesh.cpp b/src/mesh.cpp index efd627b6d..37f6a3a7c 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -53,6 +53,12 @@ std::vector> meshes; } // namespace model +#ifdef LIBMESH +namespace settings { +std::unique_ptr LMI; +} +#endif + //============================================================================== // Helper functions //============================================================================== @@ -116,7 +122,7 @@ StructuredMesh::bin_label(int bin) const { } //============================================================================== -// Untructured Mesh implementation +// Unstructured Mesh implementation //============================================================================== UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) { @@ -147,6 +153,12 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) { } +void +UnstructuredMesh::surface_bins_crossed(const Particle& p, + std::vector& bins) const { + throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."}; +} + std::string UnstructuredMesh::bin_label(int bin) const { return fmt::format("Mesh Index ({})", bin); @@ -175,96 +187,6 @@ UnstructuredMesh::to_hdf5(hid_t group) const close_group(mesh_group); } - -//============================================================================== -// RegularMesh implementation -//============================================================================== - -// RegularMesh::RegularMesh(pugi::xml_node node) -// : StructuredMesh {node} -// { -// // Determine number of dimensions for mesh -// if (check_for_node(node, "dimension")) { -// shape_ = get_node_xarray(node, "dimension"); -// int n = n_dimension_ = shape_.size(); -// if (n != 1 && n != 2 && n != 3) { -// fatal_error("Mesh must be one, two, or three dimensions."); -// } - -// // Check that dimensions are all greater than zero -// if (xt::any(shape_ <= 0)) { -// fatal_error("All entries on the element for a tally " -// "mesh must be positive."); -// } -// } - -// // Check for lower-left coordinates -// if (check_for_node(node, "lower_left")) { -// // Read mesh lower-left corner location -// lower_left_ = get_node_xarray(node, "lower_left"); -// } else { -// fatal_error("Must specify on a mesh."); -// } - -// if (check_for_node(node, "width")) { -// // Make sure both upper-right or width were specified -// if (check_for_node(node, "upper_right")) { -// fatal_error("Cannot specify both and on a mesh."); -// } - -// width_ = get_node_xarray(node, "width"); - -// // Check to ensure width has same dimensions -// auto n = width_.size(); -// if (n != lower_left_.size()) { -// fatal_error("Number of entries on must be the same as " -// "the number of entries on ."); -// } - -// // Check for negative widths -// if (xt::any(width_ < 0.0)) { -// fatal_error("Cannot have a negative on a tally mesh."); -// } - -// // Set width and upper right coordinate -// upper_right_ = xt::eval(lower_left_ + shape_ * width_); - -// } else if (check_for_node(node, "upper_right")) { -// upper_right_ = get_node_xarray(node, "upper_right"); - -// // Check to ensure width has same dimensions -// auto n = upper_right_.size(); -// if (n != lower_left_.size()) { -// fatal_error("Number of entries on must be the " -// "same as the number of entries on ."); -// } - -// // Check that upper-right is above lower-left -// if (xt::any(upper_right_ < lower_left_)) { -// fatal_error("The coordinates must be greater than " -// "the coordinates on a tally mesh."); -// } - -// // Set width -// if (shape_.size() > 0) { -// width_ = xt::eval((upper_right_ - lower_left_) / shape_); -// } -// } else { -// fatal_error("Must specify either and on a mesh."); -// } - -// // Make sure lower_left and dimension match -// if (shape_.size() > 0) { -// if (shape_.size() != lower_left_.size()) { -// fatal_error("Number of entries on must be the same " -// "as the number of entries on ."); -// } - -// // Set volume fraction -// volume_frac_ = 1.0/xt::prod(shape_)(); -// } -// } - void StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const { *in_mesh = true; @@ -1441,32 +1363,35 @@ extern "C" int openmc_add_unstructured_mesh(const char filename[], const char library[], int* id) { + + std::string lib_name(library); + std::string mesh_file(filename); bool valid_lib = false; #ifdef DAGMC - if (library == "moab") { - model::meshes.push_back(std::move(std::make_unique(filename))); + if (lib_name == "moab") { + model::meshes.push_back(std::move(std::make_unique(mesh_file))); valid_lib = true; } #endif #ifdef LIBMESH - if (library == "libmesh") { - model::meshes.push_back(std::move(std::make_unique(filename))); + if (lib_name == "libmesh") { + model::meshes.push_back(std::move(std::make_unique(mesh_file))); valid_lib = true; } #endif if (!valid_lib) { - set_errmsg("Mesh library " + std::string(library) + \ - "is not supported by this build of OpenMC"); + set_errmsg(fmt::format("Mesh library {} is not supported " + "by this build of OpenMC", lib_name)); return OPENMC_E_INVALID_ARGUMENT; } int mesh_id = 0; for (const auto& m : model::meshes) { mesh_id = std::max(m->id_, mesh_id); } - mesh_id += 1; - model::meshes.back()->id_ = mesh_id; + + model::meshes.back()->id_ = mesh_id + 1; *id = mesh_id; return 0; @@ -1630,16 +1555,16 @@ openmc_rectilinear_mesh_set_grid(int32_t index, const double* grid_x, #ifdef DAGMC -MOABUnstructuredMesh::MOABUnstructuredMesh(pugi::xml_node node) : UnstructuredMesh(node) { +MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node) { initialize(); } -MOABUnstructuredMesh::MOABUnstructuredMesh(const std::string& filename) { +MOABMesh::MOABMesh(const std::string& filename) { filename_ = filename; initialize(); } -void MOABUnstructuredMesh::initialize() { +void MOABMesh::initialize() { // create MOAB instance mbi_ = std::make_unique(); // load unstructured mesh file @@ -1677,7 +1602,7 @@ void MOABUnstructuredMesh::initialize() { } void -MOABUnstructuredMesh::build_kdtree(const moab::Range& all_tets) +MOABMesh::build_kdtree(const moab::Range& all_tets) { moab::Range all_tris; int adj_dim = 2; @@ -1712,7 +1637,7 @@ MOABUnstructuredMesh::build_kdtree(const moab::Range& all_tets) } void -MOABUnstructuredMesh::intersect_track(const moab::CartVect& start, +MOABMesh::intersect_track(const moab::CartVect& start, const moab::CartVect& dir, double track_len, std::vector& hits) const { @@ -1743,7 +1668,7 @@ MOABUnstructuredMesh::intersect_track(const moab::CartVect& start, } void -MOABUnstructuredMesh::bins_crossed(const Particle& p, +MOABMesh::bins_crossed(const Particle& p, std::vector& bins, std::vector& lengths) const { @@ -1820,7 +1745,7 @@ MOABUnstructuredMesh::bins_crossed(const Particle& p, }; moab::EntityHandle -MOABUnstructuredMesh::get_tet(const Position& r) const +MOABMesh::get_tet(const Position& r) const { moab::CartVect pos(r.x, r.y, r.z); // find the leaf of the kd-tree for this position @@ -1847,13 +1772,13 @@ MOABUnstructuredMesh::get_tet(const Position& r) const return 0; } -double MOABUnstructuredMesh::volume(int bin) const { +double MOABMesh::volume(int bin) const { return tet_volume(get_ent_handle_from_bin(bin)); } -std::string MOABUnstructuredMesh::library() const { return "moab"; } +std::string MOABMesh::library() const { return "moab"; } -double MOABUnstructuredMesh::tet_volume(moab::EntityHandle tet) const { +double MOABMesh::tet_volume(moab::EntityHandle tet) const { std::vector conn; moab::ErrorCode rval = mbi_->get_connectivity(&tet, 1, conn); if (rval != moab::MB_SUCCESS) { @@ -1869,13 +1794,8 @@ double MOABUnstructuredMesh::tet_volume(moab::EntityHandle tet) const { return 1.0 / 6.0 * (((p[1] - p[0]) * (p[2] - p[0])) % (p[3] - p[0])); } -void MOABUnstructuredMesh::surface_bins_crossed(const Particle& p, std::vector& bins) const { - // TODO: Implement triangle crossings here - throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."}; -} - int -MOABUnstructuredMesh::get_bin(Position r) const { +MOABMesh::get_bin(Position r) const { moab::EntityHandle tet = get_tet(r); if (tet == 0) { return -1; @@ -1885,7 +1805,7 @@ MOABUnstructuredMesh::get_bin(Position r) const { } void -MOABUnstructuredMesh::compute_barycentric_data(const moab::Range& tets) { +MOABMesh::compute_barycentric_data(const moab::Range& tets) { moab::ErrorCode rval; baryc_data_.clear(); @@ -1915,7 +1835,7 @@ MOABUnstructuredMesh::compute_barycentric_data(const moab::Range& tets) { } bool -MOABUnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet) const { +MOABMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet) const { moab::ErrorCode rval; @@ -1950,7 +1870,7 @@ MOABUnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle t } int -MOABUnstructuredMesh::get_bin_from_index(int idx) const { +MOABMesh::get_bin_from_index(int idx) const { if (idx >= n_bins()) { fatal_error(fmt::format("Invalid bin index: {}", idx)); } @@ -1958,25 +1878,25 @@ MOABUnstructuredMesh::get_bin_from_index(int idx) const { } int -MOABUnstructuredMesh::get_index(const Position& r, +MOABMesh::get_index(const Position& r, bool* in_mesh) const { int bin = get_bin(r); *in_mesh = bin != -1; return bin; } -int MOABUnstructuredMesh::get_index_from_bin(int bin) const { +int MOABMesh::get_index_from_bin(int bin) const { return bin; } std::pair, std::vector> -MOABUnstructuredMesh::plot(Position plot_ll, Position plot_ur) const { +MOABMesh::plot(Position plot_ll, Position plot_ur) const { // TODO: Implement mesh lines return {}; } int -MOABUnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { +MOABMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { int bin = eh - ehs_[0]; if (bin >= n_bins()) { fatal_error(fmt::format("Invalid bin: {}", bin)); @@ -1985,18 +1905,18 @@ MOABUnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { } moab::EntityHandle -MOABUnstructuredMesh::get_ent_handle_from_bin(int bin) const { +MOABMesh::get_ent_handle_from_bin(int bin) const { if (bin >= n_bins()) { fatal_error(fmt::format("Invalid bin index: ", bin)); } return ehs_[0] + bin; } -int MOABUnstructuredMesh::n_bins() const { +int MOABMesh::n_bins() const { return ehs_.size(); } -int MOABUnstructuredMesh::n_surface_bins() const { +int MOABMesh::n_surface_bins() const { // collect all triangles in the set of tets for this mesh moab::Range tris; moab::ErrorCode rval; @@ -2009,7 +1929,7 @@ int MOABUnstructuredMesh::n_surface_bins() const { } Position -MOABUnstructuredMesh::centroid(int bin) const { +MOABMesh::centroid(int bin) const { moab::ErrorCode rval; auto tet = this->get_ent_handle_from_bin(bin); @@ -2041,7 +1961,7 @@ MOABUnstructuredMesh::centroid(int bin) const { } std::pair -MOABUnstructuredMesh::get_score_tags(std::string score) const { +MOABMesh::get_score_tags(std::string score) const { moab::ErrorCode rval; // add a tag to the mesh // all scores are treated as a single value @@ -2083,11 +2003,11 @@ MOABUnstructuredMesh::get_score_tags(std::string score) const { } void -MOABUnstructuredMesh::add_score(const std::string& score) { +MOABMesh::add_score(const std::string& score) { auto score_tags = get_score_tags(score); } -void MOABUnstructuredMesh::remove_score(const std::string& score) { +void MOABMesh::remove_score(const std::string& score) { auto value_name = score + "_mean"; moab::Tag tag; moab::ErrorCode rval = mbi_->tag_get_handle(value_name.c_str(), tag); @@ -2116,7 +2036,7 @@ void MOABUnstructuredMesh::remove_score(const std::string& score) { } void -MOABUnstructuredMesh::set_score_data(const std::string& score, +MOABMesh::set_score_data(const std::string& score, std::vector values, std::vector std_dev) { auto score_tags = this->get_score_tags(score); @@ -2148,7 +2068,7 @@ MOABUnstructuredMesh::set_score_data(const std::string& score, } void -MOABUnstructuredMesh::write(std::string base_filename) const { +MOABMesh::write(std::string base_filename) const { // add extension to the base name auto filename = base_filename + ".vtk"; write_message(5, "Writing unstructured mesh {}...", filename); @@ -2165,14 +2085,6 @@ MOABUnstructuredMesh::write(std::string base_filename) const { } } -xt::xtensor -UnstructuredMesh::count_sites(const Particle::Bank* bank, - int64_t length, - bool* outside) const { - xt::xtensor out; - return out; -} - #endif #ifdef LIBMESH @@ -2195,10 +2107,7 @@ LibMesh::centroid(int bin) const { std::string LibMesh::library() const { return "libmesh"; } void LibMesh::initialize() { - // always 3 for unstructured meshes - n_dimension_ = 3; - - m_ = std::make_unique(*settings::libmesh_comm, 3); + m_ = std::make_unique(settings::LMI->comm(), 3); m_->read(filename_); m_->prepare_for_use(); @@ -2227,7 +2136,6 @@ void LibMesh::initialize() { // bounding box for the mesh bbox_ = libMesh::MeshTools::create_bounding_box(*m_); - } int LibMesh::n_bins() const { @@ -2235,21 +2143,22 @@ int LibMesh::n_bins() const { } int LibMesh::n_surface_bins() const { - // TODO: Return number of faces in the mesh here - throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."}; -} - -void -LibMesh::surface_bins_crossed(const Particle& p, - std::vector& bins) const -{ - // TODO: Implement triangle crossings here - throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."}; + int n_bins = 0; + for (int i = 0; i < this->n_bins(); i++) { + const libMesh::Elem& e = m_->elem_ref(i); + n_bins += e.n_faces(); + // if this is a boundary element, sure to count boundary faces + // twice + for (int j = 0; j < e.n_sides(); j++) { + if (e.neighbor_ptr(j) == NULL) { n_bins++; } + } + } + return n_bins; } void LibMesh::add_score(const std::string& var_name) { - // check if this is a new varaible + // check if this is a new variable std::string value_name = var_name + "_mean"; if (!variable_map_.count(value_name)) { auto& eqn_sys = equation_systems_->get_system(eq_system_name_); @@ -2258,7 +2167,7 @@ LibMesh::add_score(const std::string& var_name) { } std::string std_dev_name = var_name + "_std_dev"; - // check if this is a new varaible + // check if this is a new variable if (!variable_map_.count(std_dev_name)) { auto& eqn_sys = equation_systems_->get_system(eq_system_name_); auto var_num = eqn_sys.add_variable(std_dev_name, libMesh::CONSTANT, libMesh::MONOMIAL); @@ -2309,8 +2218,7 @@ LibMesh::set_score_data(const std::string& var_name, } void LibMesh::write(std::string filename) const { - std::cout << "Writing file : " << filename + ".e" << std::endl; - + write_message(fmt::format("Writing file: {}.e for unstructured mesh {}", filename, this->id_)); libMesh::ExodusII_IO exo(*m_); std::set systems_out = {eq_system_name_}; exo.write_discontinuous_exodusII(filename + ".e", *equation_systems_, &systems_out); @@ -2322,7 +2230,7 @@ LibMesh::bins_crossed(const Particle& p, std::vector& lengths) const { // TODO: Implement triangle crossings here - throw std::runtime_error{"LibMesh tracklength tallies are not implemented."}; + throw std::runtime_error{"Tracklength tallies on libMesh instances are not implemented."}; } int @@ -2366,7 +2274,9 @@ LibMesh::get_element_from_bin(int bin) const { return m_->elem_ptr(bin); } -double LibMesh::volume(int bin) const { return m_->elem_ref(bin).volume(); } +double LibMesh::volume(int bin) const { + return m_->elem_ref(bin).volume(); +} #endif // LIBMESH @@ -2397,7 +2307,7 @@ void read_meshes(pugi::xml_node root) model::meshes.push_back(std::make_unique(node)); #ifdef DAGMC } else if (mesh_type == "unstructured" && mesh_lib == "moab") { - model::meshes.push_back(std::make_unique(node)); + model::meshes.push_back(std::make_unique(node)); #endif #ifdef LIBMESH } else if (mesh_type == "unstructured" && mesh_lib == "libmesh") { diff --git a/src/settings.cpp b/src/settings.cpp index d54543f72..41b18942c 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -37,11 +37,6 @@ namespace openmc { namespace settings { -#ifdef LIBMESH -std::unique_ptr LMI; -const libMesh::Parallel::Communicator* libmesh_comm; -#endif - // Default values for boolean flags bool assume_separate {false}; bool check_overlaps {false};