Merge pull request enrico-dev#29 from aprilnovak/move-convergence-check

Move convergence check

include/enrico/coupled_driver.h

@@ -30,6 +30,11 @@ class CoupledDriver {
   //! Execute the coupled driver
   virtual void execute();
 
+  //! Whether the calling rank has access to the coupled solution
+  //! fields for the heat source, temperature, density, and other protected member
+  //! variables of this class.
+  virtual bool has_global_coupling_data() const = 0;
+
   //! Update the heat source for the thermal-hydraulics solver
   virtual void update_heat_source() final;
 
@@ -43,8 +48,7 @@ class CoupledDriver {
   virtual void update_density() {}
 
   //! Check convergence of the coupled solve for the current Picard iteration.
-  //! By default, derived classes will run up to the maximum number of Picard iterations.
-  virtual bool is_converged() { return false; }
+  virtual bool is_converged();
 
   enum class Norm { L1, L2, LINF }; //! Types of norms
 

include/enrico/openmc_heat_driver.h

@@ -24,12 +24,16 @@ class OpenmcHeatDriver : public CoupledDriver {
   //! \param node  XML node containing settings
   explicit OpenmcHeatDriver(MPI_Comm comm, pugi::xml_node node);
 
+  //! Whether the calling rank has access to global coupling fields. Because OpenMC
+  //! and the surrogate driver share the same communicator, and the surrogate driver does not
+  //! split up its computation among multiple ranks, we only need to check that
+  //! both communicators are active (which they always should be).
+  bool has_global_coupling_data() const override;
+
   void set_heat_source() override;
 
   void update_temperature() override;
 
-  bool is_converged() override;
-
   NeutronicsDriver& get_neutronics_driver() const override;
 
   HeatFluidsDriver & get_heat_driver() const override;

include/enrico/openmc_nek_driver.h

@@ -29,6 +29,14 @@ class OpenmcNekDriver : public CoupledDriver {
   //! Frees any data structures that need manual freeing.
   ~OpenmcNekDriver();
 
+  //! Whether the calling rank has access to global coupling fields. Because the OpenMC
+  //! and Nek communicators are assumed to overlap (though they are not the same), and
+  //! Nek broadcasts its solution onto the OpenMC ranks, we need to check that both
+  //! communicators are active.
+  //!
+  //! TODO: This won't work if the OpenMC and Nek communicators are disjoint
+  bool has_global_coupling_data() const override;
+
   void set_heat_source() override;
 
   void update_temperature() override;
@@ -39,13 +47,6 @@ class OpenmcNekDriver : public CoupledDriver {
 
   HeatFluidsDriver & get_heat_driver() const override;
 
-  //! Check convergence based on temperature field and specified epsilon
-  //!
-  //! Currently compares L_inf norm of temperature data between iterations
-  //!
-  //! \return true if L_inf norm of temperature data is less than epsilon
-  bool is_converged() override;
-
   Comm intranode_comm_; //!< The communicator representing intranode ranks
   std::unique_ptr<OpenmcDriver> openmc_driver_; //!< The OpenMC driver
   std::unique_ptr<NekDriver> nek_driver_;       //!< The Nek5000 driver

src/coupled_driver.cpp

@@ -117,6 +117,24 @@ void CoupledDriver::compute_temperature_norm(const CoupledDriver::Norm& n,
   converged = norm < epsilon_;
 }
 
+bool CoupledDriver::is_converged()
+{
+  bool converged;
+
+  // assumes that the rank 0 process is in the communicator that has access
+  // to global coupling data
+  if (comm_.rank == 0) {
+    double norm;
+    this->compute_temperature_norm(Norm::LINF, norm, converged);
+
+    std::string msg = "temperature norm_linf: " + std::to_string(norm);
+    comm_.message(msg);
+  }
+
+  comm_.Bcast(&converged, 1, MPI_CXX_BOOL);
+  return converged;
+}
+
 void CoupledDriver::update_heat_source()
 {
   // Store previous heat source solution if more than one iteration has been performed

src/openmc_heat_driver.cpp

@@ -28,6 +28,11 @@ OpenmcHeatDriver::OpenmcHeatDriver(MPI_Comm comm, pugi::xml_node node)
   init_heat_source();
 }
 
+bool OpenmcHeatDriver::has_global_coupling_data() const
+{
+  return openmc_driver_->active() && heat_driver_->active();
+}
+
 NeutronicsDriver& OpenmcHeatDriver::get_neutronics_driver() const
 {
   return *openmc_driver_;
@@ -110,59 +115,60 @@ void OpenmcHeatDriver::init_tallies()
 
 void OpenmcHeatDriver::init_temperatures()
 {
-  std::size_t n =
-    heat_driver_->n_pins_ * heat_driver_->n_axial_ * heat_driver_->n_rings();
-  temperatures_.resize({n});
-  temperatures_prev_.resize({n});
-
-  if (temperature_ic_ == Initial::neutronics) {
-    // Loop over all of the rings in the heat transfer model and set the temperature IC
-    // based on temperatures used in the OpenMC input file. More than one OpenMC cell may
-    // correspond to a particular ring, so the initial temperature set for that ring
-    // should be a volume average of the OpenMC cell temperatures.
-
-    // TODO: This initial condition used in the coupled driver does not truly represent
-    // the actual initial condition used in the OpenMC input file, since the surrogate
-    // heat solver only includes a radial dependence, though the various azimuthal
-    // segments corresponding to a single heat transfer ring may run with different
-    // initial temperatures. For this reading of the initial condition to be completely
-    // accurate, the heat solver must either be modified to include an azimuthal
-    // dependence, or a check inserted here to ensure that the initial temperatures in the
-    // azimuthal segments in the OpenMC model are all the same (i.e. no initial azimuthal
-    // dependence).
-
-    int ring_index = 0;
-    for (int i = 0; i < heat_driver_->n_pins_; ++i) {
-      for (int j = 0; j < heat_driver_->n_axial_; ++j) {
-        for (int k = 0; k < heat_driver_->n_rings(); ++k) {
-          const auto& cell_instances = ring_to_cell_inst_[ring_index];
-
-          double T_avg = 0.0;
-          double total_vol = 0.0;
-
-          for (auto idx : cell_instances) {
-            const auto& c = openmc_driver_->cells_[idx];
-            double vol = c.volume_;
-
-            total_vol += vol;
-            T_avg += c.get_temperature() * vol;
+  if (this->has_global_coupling_data()) {
+    std::size_t n =
+      heat_driver_->n_pins_ * heat_driver_->n_axial_ * heat_driver_->n_rings();
+    temperatures_.resize({n});
+    temperatures_prev_.resize({n});
+
+    if (temperature_ic_ == Initial::neutronics) {
+      // Loop over all of the rings in the heat transfer model and set the temperature IC
+      // based on temperatures used in the OpenMC input file. More than one OpenMC cell may
+      // correspond to a particular ring, so the initial temperature set for that ring should
+      // be a volume average of the OpenMC cell temperatures.
+
+      // TODO: This initial condition used in the coupled driver does not truly represent
+      // the actual initial condition used in the OpenMC input file, since the surrogate heat
+      // solver only includes a radial dependence, though the various azimuthal segments
+      // corresponding to a single heat transfer ring may run with different initial
+      // temperatures. For this reading of the initial condition to be completely accurate,
+      // the heat solver must either be modified to include an azimuthal dependence, or
+      // a check inserted here to ensure that the initial temperatures in the azimuthal
+      // segments in the OpenMC model are all the same (i.e. no initial azimuthal dependence).
+
+      int ring_index = 0;
+      for (int i = 0; i < heat_driver_->n_pins_; ++i) {
+        for (int j = 0; j < heat_driver_->n_axial_; ++j) {
+          for (int k = 0; k < heat_driver_->n_rings(); ++k) {
+            const auto& cell_instances = ring_to_cell_inst_[ring_index];
+
+            double T_avg = 0.0;
+            double total_vol = 0.0;
+
+            for (auto idx : cell_instances) {
+              const auto& c = openmc_driver_->cells_[idx];
+              double vol = c.volume_;
+
+              total_vol += vol;
+              T_avg += c.get_temperature() * vol;
+            }
+
+            T_avg /= total_vol;
+
+            int t_index = (i * heat_driver_->n_axial_ + j) * heat_driver_->n_rings() + k;
+            temperatures_prev_[t_index] = T_avg;
+            temperatures_[t_index] = T_avg;
+
+            ++ring_index;
           }
-
-          T_avg /= total_vol;
-
-          int t_index = (i * heat_driver_->n_axial_ + j) * heat_driver_->n_rings() + k;
-          temperatures_prev_[t_index] = T_avg;
-          temperatures_[t_index] = T_avg;
-
-          ++ring_index;
         }
       }
     }
-  }
 
-  if (temperature_ic_ == Initial::heat) {
-    throw std::runtime_error{"Temperature initial conditions from surrogate heat-fluids "
-                             "solver not supported."};
+    if (temperature_ic_ == Initial::heat) {
+      throw std::runtime_error{"Temperature initial conditions from surrogate heat-fluids "
+                               "solver not supported."};
+    }
   }
 }
 
@@ -244,18 +250,4 @@ void OpenmcHeatDriver::update_temperature()
   }
 }
 
-bool OpenmcHeatDriver::is_converged()
-{
-  bool converged;
-  if (comm_.rank == 0) {
-    double norm;
-    compute_temperature_norm(Norm::LINF, norm, converged);
-
-    std::string msg = "temperature norm_linf: " + std::to_string(norm);
-    comm_.message(msg);
-  }
-  err_chk(comm_.Bcast(&converged, 1, MPI_CXX_BOOL));
-  return converged;
-}
-
 } // namespace enrico

src/openmc_nek_driver.cpp

@@ -61,6 +61,11 @@ OpenmcNekDriver::~OpenmcNekDriver()
   free_mpi_datatypes();
 }
 
+bool OpenmcNekDriver::has_global_coupling_data() const
+{
+  return openmc_driver_->active() && nek_driver_->active();
+}
+
 NeutronicsDriver& OpenmcNekDriver::get_neutronics_driver() const
 {
   return *openmc_driver_;
@@ -96,14 +101,12 @@ void OpenmcNekDriver::init_mpi_datatypes()
 
 void OpenmcNekDriver::init_mappings()
 {
-  // TODO: This won't work if the Nek/OpenMC communicators are disjoint
+  if (this->has_global_coupling_data()) {
+    elem_centroids_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
+    elem_fluid_mask_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
+  }
 
   if (nek_driver_->active()) {
-    // Only the OpenMC procs get the global element centroids/fluid-identities
-    if (openmc_driver_->active()) {
-      elem_centroids_.resize(n_global_elem_);
-      elem_fluid_mask_.resize({n_global_elem_});
-    }
     // Step 1: Get global element centroids/fluid-identities on all OpenMC ranks
     // Each Nek proc finds the centroids/fluid-identities of its local elements
     Position local_element_centroids[n_local_elem_];
@@ -207,14 +210,11 @@ void OpenmcNekDriver::init_tallies()
 
 void OpenmcNekDriver::init_temperatures()
 {
-  if (nek_driver_->active() && openmc_driver_->active()) {
+  if (this->has_global_coupling_data()) {
     temperatures_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
     temperatures_prev_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
-  }
-  // TODO: This won't work if the Nek/OpenMC communicators are disjoint
-  // Only the OpenMC procs get the global temperatures
-  if (temperature_ic_ == Initial::neutronics) {
-    if (nek_driver_->active() && openmc_driver_->active()) {
+
+    if (temperature_ic_ == Initial::neutronics) {
       // Loop over the OpenMC cells, then loop over the global Nek elements
       // corresponding to that cell and assign the OpenMC cell temperature to
       // the correct index in the temperatures_ array. This mapping assumes that
@@ -245,13 +245,11 @@ void OpenmcNekDriver::init_temperatures()
 
 void OpenmcNekDriver::init_volumes()
 {
-  // TODO: This won't work if the Nek/OpenMC communicators are disjoint
+  if (this->has_global_coupling_data()) {
+    elem_volumes_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
+  }
 
-  // Only the OpenMC procs get the global element volumes (gev)
   if (nek_driver_->active()) {
-    if (openmc_driver_->active()) {
-      elem_volumes_.resize({n_global_elem_});
-    }
     // Every Nek proc gets its local element volumes (lev)
     double local_elem_volumes[n_local_elem_];
     for (int i = 0; i < n_local_elem_; ++i) {
@@ -274,16 +272,15 @@ void OpenmcNekDriver::init_volumes()
 
 void OpenmcNekDriver::init_elem_densities()
 {
-  // TODO: This won't work if the Nek/OpenMC communicators are disjoint
-  if (nek_driver_->active() && openmc_driver_->active()) {
+  if (this->has_global_coupling_data()) {
     elem_densities_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
   }
   update_elem_densities();
 }
 
 void OpenmcNekDriver::init_elem_fluid_mask()
 {
-  if (nek_driver_->active() && openmc_driver_->active()) {
+  if (this->has_global_coupling_data()) {
     elem_fluid_mask_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
   }
   if (nek_driver_->active()) {
@@ -355,12 +352,11 @@ void OpenmcNekDriver::set_heat_source()
 
 void OpenmcNekDriver::update_temperature()
 {
-  if (nek_driver_->active()) {
-    // Copy previous
-    if (openmc_driver_->active()) {
-      std::copy(temperatures_.begin(), temperatures_.end(), temperatures_prev_.begin());
-    }
+  if (this->has_global_coupling_data()) {
+    std::copy(temperatures_.begin(), temperatures_.end(), temperatures_prev_.begin());
+  }
 
+  if (nek_driver_->active()) {
     auto t = nek_driver_->temperature();
     if (nek_driver_->comm_.rank == 0) {
       temperatures_ = t;
@@ -443,21 +439,6 @@ void OpenmcNekDriver::update_cell_densities()
   }
 }
 
-bool OpenmcNekDriver::is_converged()
-{
-  bool converged;
-  // WARNING: Assumes that OpenmcNekDriver rank 0 is in openmc_driver_->comm
-  if (comm_.rank == 0) {
-    double norm;
-    compute_temperature_norm(Norm::LINF, norm, converged);
-
-    std::string msg = "temperature norm_linf: " + std::to_string(norm);
-    comm_.message(msg);
-  }
-  err_chk(comm_.Bcast(&converged, 1, MPI_CXX_BOOL));
-  return converged;
-}
-
 void OpenmcNekDriver::free_mpi_datatypes()
 {
   MPI_Type_free(&position_mpi_datatype);