crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Enable adaptive mesh support on libMesh tallies (#3185)

Browse source 
Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>
This commit is contained in:
Kevin Sawatzky 2024-11-23 08:09:47 -06:00 • committed by GitHub
parent ae37d6c0da
commit dd01c40ae1
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
2 changed files with 82 additions and 15 deletions
Download patch file Download diff file Expand all files Collapse all files

11
include/openmc/mesh.h
View file

@ -948,6 +948,7 @@ public:
private:
void initialize() override;
void set_mesh_pointer_from_filename(const std::string& filename);
void build_eqn_sys();
// Methods
@ -966,7 +967,8 @@ private:
vector<unique_ptr<libMesh::PointLocatorBase>>
pl_; //!< per-thread point locators
unique_ptr<libMesh::EquationSystems>
equation_systems_; //!< pointer to the equation systems of the mesh
equation_systems_; //!< 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>
@ -975,6 +977,13 @@ private:
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

86
src/mesh.cpp
View file

@ -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
@ -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);
@ -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);
@ -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()
{
@ -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);
@ -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 ID space, 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(), -1);
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();
@ -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
@ -3079,6 +3098,18 @@ int LibMesh::n_surface_bins() const
void LibMesh::add_score(const std::string& var_name)
{
if (adaptive_) {
warning(fmt::format(
"Exodus output cannot be provided as unstructured mesh {} is adaptive.",
this->id_));
return;
}
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)) {
@ -3100,14 +3131,28 @@ 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_) {
warning(fmt::format(
"Exodus output cannot be provided as unstructured mesh {} is adaptive.",
this->id_));
return;
}
if (!equation_systems_) {
build_eqn_sys();
}
auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
if (!eqn_sys.is_initialized()) {
@ -3125,6 +3170,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
@ -3143,6 +3192,14 @@ void LibMesh::set_score_data(const std::string& var_name,
void LibMesh::write(const std::string& filename) const
{
if (adaptive_) {
warning(fmt::format(
"Exodus output cannot be provided as unstructured mesh {} is adaptive.",
this->id_));
return;
}
write_message(fmt::format(
"Writing file: {}.e for unstructured mesh {}", filename, this->id_));
libMesh::ExodusII_IO exo(*m_);
@ -3176,7 +3233,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));
}
@ -3191,7 +3249,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