From 6973ce28d0be007ec4c39f03e33ede1be1e29957 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 4 Apr 2020 13:06:31 -0500 Subject: [PATCH] Updates to the LibMesh methods for adding variables and writing. --- CMakeLists.txt | 2 +- include/openmc/mesh.h | 18 +++++++++-- include/openmc/state_point.h | 2 +- src/mesh.cpp | 60 +++++++++++++++++++++++++++--------- src/state_point.cpp | 34 +++++++++----------- 5 files changed, 77 insertions(+), 39 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index fa7988304..606a8518e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -395,7 +395,7 @@ if (MPI_ENABLED) target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI) endif() if (libmesh) - target_compile_definitions(libopenmc PUBLIC -DLIBMESH) + target_compile_definitions(libopenmc PRIVATE LIBMESH) endif() # Set git SHA1 hash as a compile definition diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 318fa3fe7..e8ed7d320 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -282,6 +282,16 @@ public: std::string bin_label(int bin) const override; + //! Add a variable to the libmesh mesh instance + virtual void add_score(const std::string& var_name) = 0; + + //! Set the value of a bin for a variable on the libmesh mesh instance + virtual void set_score_data(const std::string& var_name, + std::vector values, + std::vector std_dev) const = 0; + + virtual void write(std::string filename) const = 0; + std::string filename_; }; @@ -495,15 +505,17 @@ public: bool intersects(Position& r0, Position r1, int* ijk) const; - void write(const std::string& filename) const; + void write(std::string filename) const override; void to_hdf5(hid_t group) const; //! Add a variable to the libmesh mesh instance - void add_variable(const std::string& var_name); + void add_score(const std::string& var_name) override; //! Set the value of a bin for a variable on the libmesh mesh instance - void set_variable(const std::string& var_name, int bin, double value); + void set_score_data(const std::string& var_name, + std::vector values, + std::vector std_dev) const override; private: diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index 53f072b85..e3cd790db 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -19,7 +19,7 @@ void read_source_bank(hid_t group_id, std::vector& sites, bool d void write_tally_results_nr(hid_t file_id); void restart_set_keff(); -#ifdef DAGMC +#ifdef LIBMESH void write_unstructured_mesh_results(); #endif diff --git a/src/mesh.cpp b/src/mesh.cpp index 40bd5858c..3050bdeae 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2277,30 +2277,60 @@ LibMesh::get_element_from_bin(int bin) const { } void -LibMesh::add_variable(const std::string& var_name) { +LibMesh::add_score(const std::string& var_name) { // check if this is a new varaible - if (!variable_map_.count(var_name)) { + std::string value_name = var_name + "_value"; + if (!variable_map_.count(value_name)) { auto& eqn_sys = equation_systems_->get_system(eq_system_name_); - auto var_num = eqn_sys.add_variable(var_name, libMesh::CONSTANT, libMesh::MONOMIAL); - variable_map_[var_name] = var_num; + auto var_num = eqn_sys.add_variable(value_name, libMesh::CONSTANT, libMesh::MONOMIAL); + variable_map_[value_name] = var_num; + equation_systems_->init(); + } + + std::string std_dev_name = var_name + "_std_dev"; + // check if this is a new varaible + 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); + variable_map_[std_dev_name] = var_num; equation_systems_->init(); } } void -LibMesh::set_variable(const std::string& var_name, int bin, double value) { - // look up the variable - unsigned int var_num = variable_map_.at(var_name); - auto& eqn_sys = equation_systems_->get_system(eq_system_name_); - const libMesh::DofMap& dof_map = eqn_sys.get_dof_map(); - auto e = m_->elem_ptr(bin); - std::vector dof_indices; - dof_map.dof_indices(e, dof_indices, var_num); - Ensures(dof_indices.size() == 1); - eqn_sys.solution->set(dof_indices[0], value); +LibMesh::set_score_data(const std::string& var_name, + std::vector values, + std::vector std_dev) const +{ + auto& eqn_sys = equation_systems_->get_system(eq_system_name_); + const libMesh::DofMap& dof_map = eqn_sys.get_dof_map(); + + // look up the value variable + std::string value_name = var_name + "_value"; + unsigned int value_num = variable_map_.at(value_name); + // look up the std dev variable + std::string std_dev_name = var_name + "_std_dev"; + unsigned int std_dev_num = variable_map_.at(std_dev_name); + + for (int bin = 0; bin < this->n_bins(); bin++) { + auto e = m_->elem_ptr(bin); + + // set value + std::vector value_dof_indices; + dof_map.dof_indices(e, value_dof_indices, value_num); + Ensures(value_dof_indices.size() == 1); + eqn_sys.solution->set(value_dof_indices[0], values[bin]); + + // set std dev + std::vector std_dev_dof_indices; + dof_map.dof_indices(e, std_dev_dof_indices, std_dev_num); + Ensures(std_dev_dof_indices.size() == 1); + eqn_sys.solution->set(std_dev_dof_indices[0], std_dev[bin]); + + } } -void LibMesh::write(const std::string& filename) const { +void LibMesh::write(std::string filename) const { libMesh::ExodusII_IO exo(*m_); std::set systems_out = {eq_system_name_}; exo.write_discontinuous_exodusII(filename, *equation_systems_, &systems_out); diff --git a/src/state_point.cpp b/src/state_point.cpp index 956518c1c..b50b752fa 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -770,31 +770,27 @@ void read_source_bank(hid_t group_id, std::vector& sites, bool d H5Tclose(banktype); } -#ifdef DAGMC +#ifdef LIBMESH void write_unstructured_mesh_results() { - + std::cout << "Checking for unstructured meshes to write..." << std::endl; for (auto& tally : model::tallies) { std::vector tally_scores; for (auto filter_idx : tally->filters()) { auto& filter = model::tally_filters[filter_idx]; - if (filter->type() != "mesh") continue; - - // check if the filter uses an unstructured mesh - auto mesh_filter = dynamic_cast(filter.get()); - auto mesh_idx = mesh_filter->mesh(); - auto umesh = dynamic_cast(model::meshes[mesh_idx].get()); - - if (!umesh) continue; - - // if this tally has more than one filter, print - // warning and skip writing the mesh - if (tally->filters().size() > 1) { - warning(fmt::format("Skipping unstructured mesh writing for tally " - "{}. More than one filter is present on the tally.", - tally->id_)); - break; - } + if (filter->type() == "mesh") { + auto mesh_filter = dynamic_cast(filter.get()); + auto& mesh = model::meshes[mesh_filter->mesh()]; + auto umesh = dynamic_cast(mesh.get()); + if (umesh) { + for (int i = 0; i < tally->scores_.size(); i++) { + // get values for this score and create vectors for marking + // up the mesh + std::vector vals_vec, sum_sq_vec; + for (int j = 0; j < tally->results_.shape()[0]; j++) { + vals_vec.push_back(tally->results_(j, i, TallyResult::SUM) / tally->n_realizations_); + sum_sq_vec.push_back(tally->results_(j , i, TallyResult::SUM_SQ) - std::pow(vals_vec[i], 2)/ tally->n_realizations_); + } int n_realizations = tally->n_realizations_;