From ab85eb525b58726aa8c272a8559b2ee17699f2ff Mon Sep 17 00:00:00 2001 From: aprilnovak Date: Sun, 9 Jun 2019 15:59:28 -0700 Subject: [PATCH] Add current and previous Picard densities to CoupledDriver and update OpenmcNekDriver syntax accordingly. Refs #17 --- include/enrico/coupled_driver.h | 12 ++++++ include/enrico/openmc_nek_driver.h | 17 ++++----- src/openmc_nek_driver.cpp | 60 +++++++++++++----------------- 3 files changed, 45 insertions(+), 44 deletions(-) diff --git a/include/enrico/coupled_driver.h b/include/enrico/coupled_driver.h index 590f2ac1..c54d04a8 100644 --- a/include/enrico/coupled_driver.h +++ b/include/enrico/coupled_driver.h @@ -107,6 +107,9 @@ class CoupledDriver { //! Initialize current and previous Picard temperature fields virtual void init_temperatures() {} + //! Initialize current and previous Picard density fields + virtual void init_densities() {} + //! Initialize current and previous Picard heat source fields. Note that //! because the neutronics solver is assumed to run first, that no initial //! condition is required for the heat source. So, unlike init_temperatures(), @@ -122,6 +125,15 @@ class CoupledDriver { xt::xtensor temperatures_prev_; //!< Previous Picard iteration temperature + //! Current Picard iteration density; this density is the density + //! computed by the thermal-hydraulic solver, and data mappings may result in + //! a different density actually used in the neutronics solver. For example, + //! the entries in this xtensor may be averaged over neutronics cells to give + //! the density used by the neutronics solver. + xt::xtensor densities_; + + xt::xtensor densities_prev_; //!< Previous Picard iteration density + //! Current Picard iteration heat source; this heat source is the heat source //! computed by the neutronics solver, and data mappings may result in a different //! heat source actually used in the heat solver. For example, the entries in this diff --git a/include/enrico/openmc_nek_driver.h b/include/enrico/openmc_nek_driver.h index a1cd0951..ff255857 100644 --- a/include/enrico/openmc_nek_driver.h +++ b/include/enrico/openmc_nek_driver.h @@ -67,6 +67,13 @@ class OpenmcNekDriver : public CoupledDriver { //! to an MPI_Gatherv operation on Nek5000's local elements. void init_heat_source() override; + //! Initialize global density buffers on all OpenMC ranks. + //! + //! These arrays store the ensity of Nek's global elements. These are **not** + //! ordered by Nek's global element indices. Rather, these are ordered according + //! to an MPI_Gatherv operation on Nek5000's local elements. + void init_densities() override; + //! Initialize global fluid masks on all OpenMC ranks. //! //! These arrays store the dimensionless source of Nek's global elements. These are **not** @@ -90,19 +97,9 @@ class OpenmcNekDriver : public CoupledDriver { //! Initialize global volume buffers for OpenMC ranks void init_volumes(); - //! Allocate space for the global volume buffers in OpenMC ranks - //! Currently, the density values are uninitialized. - void init_elem_densities(); - //! Frees the MPI datatypes (currently, only position_mpi_datatype) void free_mpi_datatypes(); - //! Updates elem_densities_ from Nek5000's density data - void update_elem_densities(); - - //! Updates cell_densities_ from elem_densities_. - void update_cell_densities(); - std::unique_ptr openmc_driver_; //!< The OpenMC driver std::unique_ptr nek_driver_; //!< The Nek5000 driver diff --git a/src/openmc_nek_driver.cpp b/src/openmc_nek_driver.cpp index 2b885e26..bfdd24e6 100644 --- a/src/openmc_nek_driver.cpp +++ b/src/openmc_nek_driver.cpp @@ -52,7 +52,7 @@ OpenmcNekDriver::OpenmcNekDriver(MPI_Comm comm, pugi::xml_node node) init_cell_fluid_mask(); init_temperatures(); - init_elem_densities(); + init_densities(); init_heat_source(); } @@ -252,12 +252,13 @@ void OpenmcNekDriver::init_volumes() } } -void OpenmcNekDriver::init_elem_densities() +void OpenmcNekDriver::init_densities() { if (this->has_global_coupling_data()) { - elem_densities_.resize({gsl::narrow(n_global_elem_)}); + densities_.resize({gsl::narrow(n_global_elem_)}); + densities_prev_.resize({gsl::narrow(n_global_elem_)}); } - update_elem_densities(); + update_density(); } void OpenmcNekDriver::init_elem_fluid_mask() @@ -375,47 +376,38 @@ void OpenmcNekDriver::update_temperature() void OpenmcNekDriver::update_density() { - // Must be done in this order. Hence, update_density is public, whereas - // update_elem_densities and update_cell_densities are private - update_elem_densities(); - update_cell_densities(); -} + if (this->has_global_coupling_data()) { + std::copy(densities_.begin(), densities_.end(), densities_prev_.begin()); + } -void OpenmcNekDriver::update_elem_densities() -{ if (nek_driver_->active()) { // On Nek's master rank, d gets global data. On Nek's other ranks, d is empty. auto d = nek_driver_->density(); - // Update elem_densities_ on Nek's master rank only. if (nek_driver_->comm_.rank == 0) { - elem_densities_ = d; + densities_ = d; } + // Since OpenMC's and Nek's master ranks are the same, we know that elem_densities_ on // OpenMC's master rank were updated. Now we broadcast to the other OpenMC ranks. // TODO: This won't work if the Nek/OpenMC communicators are disjoint if (openmc_driver_->active()) { - openmc_driver_->comm_.Bcast(elem_densities_.data(), n_global_elem_, MPI_DOUBLE); - } - } -} - -void OpenmcNekDriver::update_cell_densities() -{ - if (nek_driver_->active() && openmc_driver_->active()) { - // Update cell densities for fluid cells. Use cell_fluid_mask_ to do this - // TODO: Might be able to use xtensor masking to do some of this - for (int i = 0; i < openmc_driver_->cells_.size(); ++i) { - if (cell_fluid_mask_[i] == 1) { - auto& c = openmc_driver_->cells_[i]; - double average_density = 0.0; - double total_vol = 0.0; - for (int e : mat_to_elems_.at(c.material_index_)) { - average_density += elem_densities_[e] * elem_volumes_[e]; - total_vol += elem_volumes_[e]; + openmc_driver_->comm_.Bcast(densities_.data(), n_global_elem_, MPI_DOUBLE); + + // For each OpenMC material in a fluid cell, volume average the densities and set + // TODO: Might be able to use xtensor masking to do some of this + for (int i = 0; i < openmc_driver_->cells_.size(); ++i) { + if (cell_fluid_mask_[i] == 1) { + auto& c = openmc_driver_->cells_[i]; + double average_density = 0.0; + double total_vol = 0.0; + for (int e : mat_to_elems_.at(c.material_index_)) { + average_density += densities_[e] * elem_volumes_[e]; + total_vol += elem_volumes_[e]; + } + double density = average_density / total_vol; + Ensures(density > 0.0); + c.set_density(average_density / total_vol); } - double density = average_density / total_vol; - Ensures(density > 0.0); - c.set_density(average_density / total_vol); } } }