Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Enable adaptive mesh support on libMesh tallies #3185

Open
wants to merge 6 commits into
base: develop
Choose a base branch
from
Open

Enable adaptive mesh support on libMesh tallies #3185

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
34cc1c1
Allow libMesh meshes with an AMR hierarchy.
nuclearkevin Nov 3, 2024
7d569de
Swap the unordered map for a vector
nuclearkevin Nov 3, 2024
046aff5
Slight performance improvements and improved code clarity
nuclearkevin Nov 3, 2024
06e87a8
Improve error messages and run clang-format.
nuclearkevin Nov 4, 2024
df50cd2
Remove constructor flag
nuclearkevin Nov 7, 2024
7fa37b1
Add missing build_eqn_sys call.
nuclearkevin Nov 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 11 additions & 2 deletions include/openmc/mesh.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -948,6 +948,7 @@ class LibMesh : public UnstructuredMesh {
private:
void initialize() override;
void set_mesh_pointer_from_filename(const std::string& filename);
void build_eqn_sys();

// Methods

Expand All @@ -965,8 +966,9 @@ class LibMesh : public UnstructuredMesh {
//!< during intialization
vector<unique_ptr<libMesh::PointLocatorBase>>
pl_; //!< per-thread point locators
unique_ptr<libMesh::EquationSystems>
equation_systems_; //!< pointer to the equation systems of the mesh
unique_ptr<libMesh::EquationSystems> equation_systems_ =
nullptr; //!< pointer to the libMesh EquationSystems
//!< instance
std::string
eq_system_name_; //!< name of the equation system holding OpenMC results
std::unordered_map<std::string, unsigned int>
Expand All @@ -975,6 +977,13 @@ class LibMesh : public UnstructuredMesh {
libMesh::BoundingBox bbox_; //!< bounding box of the mesh
libMesh::dof_id_type
first_element_id_; //!< id of the first element in the mesh

const bool adaptive_; //!< whether this mesh has adaptivity enabled or not
std::vector<libMesh::dof_id_type>
bin_to_elem_map_; //!< mapping bin indices to dof indices for active
//!< elements
std::vector<int> elem_to_bin_map_; //!< mapping dof indices to bin indices for
//!< active elements
};

#endif
Expand Down
80 changes: 66 additions & 14 deletions src/mesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2922,7 +2922,7 @@ void MOABMesh::write(const std::string& base_filename) const

const std::string LibMesh::mesh_lib_type = "libmesh";

LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)
LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node), adaptive_(false)
{
// filename_ and length_multiplier_ will already be set by the
// UnstructuredMesh constructor
Expand All @@ -2933,6 +2933,7 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)

// create the mesh from a pointer to a libMesh Mesh
LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier)
: adaptive_(input_mesh.n_active_elem() != input_mesh.n_elem())
{
m_ = &input_mesh;
set_length_multiplier(length_multiplier);
Expand All @@ -2941,6 +2942,7 @@ LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier)

// create the mesh from an input file
LibMesh::LibMesh(const std::string& filename, double length_multiplier)
: adaptive_(false)
{
set_mesh_pointer_from_filename(filename);
set_length_multiplier(length_multiplier);
Expand All @@ -2955,6 +2957,15 @@ void LibMesh::set_mesh_pointer_from_filename(const std::string& filename)
m_->read(filename_);
}

// build a libMesh equation system for storing values
void LibMesh::build_eqn_sys()
{
eq_system_name_ = fmt::format("mesh_{}_system", id_);
equation_systems_ = make_unique<libMesh::EquationSystems>(*m_);
libMesh::ExplicitSystem& eq_sys =
equation_systems_->add_system<libMesh::ExplicitSystem>(eq_system_name_);
}

// intialize from mesh file
void LibMesh::initialize()
{
Expand Down Expand Up @@ -2982,13 +2993,6 @@ void LibMesh::initialize()
filename_));
}

// create an equation system for storing values
eq_system_name_ = fmt::format("mesh_{}_system", id_);

equation_systems_ = make_unique<libMesh::EquationSystems>(*m_);
libMesh::ExplicitSystem& eq_sys =
equation_systems_->add_system<libMesh::ExplicitSystem>(eq_system_name_);

for (int i = 0; i < num_threads(); i++) {
pl_.emplace_back(m_->sub_point_locator());
pl_.back()->set_contains_point_tol(FP_COINCIDENT);
Expand All @@ -2999,6 +3003,21 @@ void LibMesh::initialize()
auto first_elem = *m_->elements_begin();
first_element_id_ = first_elem->id();

// if the mesh is adaptive elements aren't guaranteed by libMesh to be
// contiguous in memory, so we need to map from bin indices (defined over
// active elements) to global dof ids
if (adaptive_) {
bin_to_elem_map_.reserve(m_->n_active_local_elem());
elem_to_bin_map_.resize(m_->n_local_elem(), 0);
for (auto it = m_->active_local_elements_begin();
it != m_->active_local_elements_end(); it++) {
auto elem = *it;

bin_to_elem_map_.push_back(elem->id());
elem_to_bin_map_[elem->id()] = bin_to_elem_map_.size() - 1;
}
}

// bounding box for the mesh for quick rejection checks
bbox_ = libMesh::MeshTools::create_bounding_box(*m_);
libMesh::Point ll = bbox_.min();
Expand Down Expand Up @@ -3056,7 +3075,7 @@ std::string LibMesh::library() const

int LibMesh::n_bins() const
{
return m_->n_elem();
return m_->n_active_elem();
}

int LibMesh::n_surface_bins() const
Expand All @@ -3079,6 +3098,16 @@ int LibMesh::n_surface_bins() const

void LibMesh::add_score(const std::string& var_name)
{
if (adaptive_) {
fatal_error(fmt::format(
"Exodus output cannot be provided as unstructured mesh {} is adaptive.",
this->id_));
}

if (!equation_systems_) {
build_eqn_sys();
}

// check if this is a new variable
std::string value_name = var_name + "_mean";
if (!variable_map_.count(value_name)) {
Expand All @@ -3100,14 +3129,26 @@ void LibMesh::add_score(const std::string& var_name)

void LibMesh::remove_scores()
{
auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
eqn_sys.clear();
variable_map_.clear();
if (equation_systems_) {
auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
eqn_sys.clear();
variable_map_.clear();
}
}

void LibMesh::set_score_data(const std::string& var_name,
const vector<double>& values, const vector<double>& std_dev)
{
if (adaptive_) {
fatal_error(fmt::format(
"Exodus output cannot be provided as unstructured mesh {} is adaptive.",
this->id_));
}

if (!equation_systems_) {
build_eqn_sys();
}

auto& eqn_sys = equation_systems_->get_system(eq_system_name_);

if (!eqn_sys.is_initialized()) {
Expand All @@ -3125,6 +3166,10 @@ void LibMesh::set_score_data(const std::string& var_name,

for (auto it = m_->local_elements_begin(); it != m_->local_elements_end();
it++) {
if (!(*it)->active()) {
continue;
}

auto bin = get_bin_from_element(*it);

// set value
Expand All @@ -3143,6 +3188,12 @@ void LibMesh::set_score_data(const std::string& var_name,

void LibMesh::write(const std::string& filename) const
{
if (adaptive_) {
fatal_error(fmt::format(
"Exodus output cannot be provided as unstructured mesh {} is adaptive.",
this->id_));
}

write_message(fmt::format(
"Writing file: {}.e for unstructured mesh {}", filename, this->id_));
libMesh::ExodusII_IO exo(*m_);
Expand Down Expand Up @@ -3176,7 +3227,8 @@ int LibMesh::get_bin(Position r) const

int LibMesh::get_bin_from_element(const libMesh::Elem* elem) const
{
int bin = elem->id() - first_element_id_;
int bin =
adaptive_ ? elem_to_bin_map_[elem->id()] : elem->id() - first_element_id_;
if (bin >= n_bins() || bin < 0) {
fatal_error(fmt::format("Invalid bin: {}", bin));
}
Expand All @@ -3191,7 +3243,7 @@ std::pair<vector<double>, vector<double>> LibMesh::plot(

const libMesh::Elem& LibMesh::get_element_from_bin(int bin) const
{
return m_->elem_ref(bin);
return adaptive_ ? m_->elem_ref(bin_to_elem_map_.at(bin)) : m_->elem_ref(bin);
}

double LibMesh::volume(int bin) const
Expand Down