Add current and previous Picard densities to CoupledDriver and update…

… OpenmcNekDriver syntax accordingly. Refs enrico-dev#17

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<double, 1> 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<double, 1> densities_;
+
+  xt::xtensor<double, 1> 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

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();
-
   //! MPI datatype for sending/receiving Position objects.
   MPI_Datatype position_mpi_datatype;
 

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();
 }
 
@@ -270,12 +270,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<std::size_t>(n_global_elem_)});
+    densities_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
+    densities_prev_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
   }
-  update_elem_densities();
+  update_density();
 }
 
 void OpenmcNekDriver::init_elem_fluid_mask()
@@ -393,42 +394,33 @@ 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);
+      openmc_driver_->comm_.Bcast(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
+    // 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 += elem_densities_[e] * elem_volumes_[e];
+          average_density += densities_[e] * elem_volumes_[e];
           total_vol += elem_volumes_[e];
         }
         double density = average_density / total_vol;