Move update of density to base coupled class. Refs enrico-dev#15

Author: aprilnovak

include/enrico/coupled_driver.h

@@ -48,7 +48,10 @@ class CoupledDriver {
   virtual void set_temperature() {};
 
   //! Update the density for the neutronics solver
-  virtual void update_density() {}
+  void update_density();
+
+  //! Update the density for the neutronics solver
+  virtual void set_density() {}
 
   //! Check convergence of the coupled solve for the current Picard iteration.
   virtual bool is_converged();

include/enrico/openmc_heat_driver.h

@@ -34,7 +34,7 @@ class OpenmcHeatDriver : public CoupledDriver {
 
   void set_temperature() override;
 
-  void update_density() override;
+  void set_density() override;
 
   NeutronicsDriver& get_neutronics_driver() const override;
 

include/enrico/openmc_nek_driver.h

@@ -42,7 +42,7 @@ class OpenmcNekDriver : public CoupledDriver {
 
   void set_temperature() override;
 
-  void update_density() override;
+  void set_density() override;
 
   NeutronicsDriver& get_neutronics_driver() const override;
 

src/coupled_driver.cpp

@@ -194,4 +194,24 @@ void CoupledDriver::update_temperature()
   set_temperature();
 }
 
+void CoupledDriver::update_density()
+{
+  // Store previous density solution; a previous solution will always be present
+  // because a density IC is set and the neutronics solver runs first
+  if (has_global_coupling_data()) {
+    std::copy(densities_.begin(), densities_.end(), densities_prev_.begin());
+  }
+
+  auto& heat = get_heat_driver();
+  if (heat.active()) {
+    auto d = heat.density();
+
+    if (heat.has_coupling_data())
+      densities_ = d;
+  }
+
+  // Set density in the neutronics solver
+  set_density();
+}
+
 } // namespace enrico

src/openmc_heat_driver.cpp

@@ -339,14 +339,8 @@ void OpenmcHeatDriver::set_heat_source()
   }
 }
 
-void OpenmcHeatDriver::update_density()
+void OpenmcHeatDriver::set_density()
 {
-  if (this->has_global_coupling_data()) {
-    std::copy(densities_.begin(), densities_.end(), densities_prev_.begin());
-  }
-
-  densities_ = heat_driver_->density();
-
   for (gsl::index i = n_solid_cells_; i < n_solid_cells_ + n_fluid_cells_; ++i) {
     const auto& c = openmc_driver_->cells_[i];
     const auto& elems = cell_inst_to_elem_[i];

src/openmc_nek_driver.cpp

@@ -416,21 +416,9 @@ void OpenmcNekDriver::set_temperature()
   }
 }
 
-void OpenmcNekDriver::update_density()
+void OpenmcNekDriver::set_density()
 {
-  if (this->has_global_coupling_data()) {
-    std::copy(densities_.begin(), densities_.end(), densities_prev_.begin());
-  }
-
   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_->has_coupling_data()) {
-      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