Added underrelaxation on the temperature. Refs enrico-dev#9

include/enrico/coupled_driver.h

@@ -42,7 +42,10 @@ class CoupledDriver {
   virtual void set_heat_source() {}
 
   //! Update the temperature for the neutronics solver
-  virtual void update_temperature() {}
+  virtual void update_temperature() final;
+
+  //! Set the temperature in the neutronics solver
+  virtual void set_temperature() {};
 
   //! Update the density for the neutronics solver
   virtual void update_density() {}
@@ -91,9 +94,14 @@ class CoupledDriver {
   //! Picard iteration convergence tolerance, defaults to 1e-3 if not set
   double epsilon_{1e-3};
 
-  //! Constant relaxation factor, defaults to 1.0 (standard Picard) if not set
+  //! Constant relaxation factor for the heat source,
+  //! defaults to 1.0 (standard Picard) if not set
   double alpha_{1.0};
 
+  //! Constant relaxation factor for the temperature, defaults to the
+  //! relaxation aplied to the heat source if not set
+  double alpha_T_{alpha_};
+
   //! Enumeration of available temperature initial condition specifications.
   //! 'neutronics' sets temperature condition from the neutronics input files,
   //! while 'heat' sets temperature based on a thermal-fluids input (or restart) file.

include/enrico/openmc_heat_driver.h

@@ -32,7 +32,7 @@ class OpenmcHeatDriver : public CoupledDriver {
 
   void set_heat_source() override;
 
-  void update_temperature() override;
+  void set_temperature() override;
 
   NeutronicsDriver& get_neutronics_driver() const override;
 

include/enrico/openmc_nek_driver.h

@@ -40,7 +40,7 @@ class OpenmcNekDriver : public CoupledDriver {
 
   void set_heat_source() override;
 
-  void update_temperature() override;
+  void set_temperature() override;
 
   void update_density() override;
 

src/coupled_driver.cpp

@@ -22,6 +22,9 @@ CoupledDriver::CoupledDriver(MPI_Comm comm, pugi::xml_node node)
   if (node.child("alpha"))
     alpha_ = node.child("alpha").text().as_double();
 
+  if (node.child("alpha_T"))
+    alpha_T_ = node.child("alpha_T").text().as_double();
+
   if (node.child("temperature_ic")) {
     auto s = std::string{node.child_value("temperature_ic")};
 
@@ -158,4 +161,25 @@ void CoupledDriver::update_heat_source()
   set_heat_source();
 }
 
+void CoupledDriver::update_temperature()
+{
+  // Store previous temperature solution; a previous solution will always be present
+  // because a temperature IC is set and the neutronics solver runs first
+  if (has_global_coupling_data()) {
+    std::copy(temperatures_.begin(), temperatures_.end(), temperatures_prev_.begin());
+  }
+
+  // Compute the next iterate of the temperature
+  auto& heat = get_heat_driver();
+  if (heat.active()) {
+    auto t = heat.temperature();
+
+    if (heat.has_coupling_data())
+      temperatures_ = alpha_T_ * t + (1.0 - alpha_T_) * temperatures_prev_;
+  }
+
+  // Set temperature in the neutronics solver
+  set_temperature();
+}
+
 } // namespace enrico

src/openmc_heat_driver.cpp

@@ -208,15 +208,11 @@ void OpenmcHeatDriver::set_heat_source()
   }
 }
 
-void OpenmcHeatDriver::update_temperature()
+void OpenmcHeatDriver::set_temperature()
 {
-  std::copy(temperatures_.begin(), temperatures_.end(), temperatures_prev_.begin());
-
   const auto& r_fuel = heat_driver_->r_grid_fuel_;
   const auto& r_clad = heat_driver_->r_grid_clad_;
 
-  temperatures_ = heat_driver_->temperature();
-
   // For each OpenMC material, volume average temperatures and set
   for (int i = 0; i < openmc_driver_->cells_.size(); ++i) {
     // Get cell instance

src/openmc_nek_driver.cpp

@@ -332,18 +332,9 @@ void OpenmcNekDriver::set_heat_source()
   }
 }
 
-void OpenmcNekDriver::update_temperature()
+void OpenmcNekDriver::set_temperature()
 {
-  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;
-    }
-
     if (openmc_driver_->active()) {
       // Broadcast global_element_temperatures onto all the OpenMC procs
       openmc_driver_->comm_.Bcast(temperatures_.data(), n_global_elem_, MPI_DOUBLE);