Use subchannel density and fluid temperature solutions in OpenMC from…

… surrogate. Refs #44

include/enrico/openmc_heat_driver.h

@@ -34,6 +34,8 @@ class OpenmcHeatDriver : public CoupledDriver {
 
   void set_temperature() override;
 
+  void update_density() override;
+
   NeutronicsDriver& get_neutronics_driver() const override;
 
   HeatFluidsDriver & get_heat_driver() const override;
@@ -42,9 +44,15 @@ class OpenmcHeatDriver : public CoupledDriver {
   std::unordered_map<int, std::vector<int>> ring_to_cell_inst_;
   std::unordered_map<int, std::vector<int>> cell_inst_to_ring_;
 
+  // Mapping of surrogate fluid elements to OpenMC cell instances and vice versa
+  std::unordered_map<int, std::vector<int>> elem_to_cell_inst_;
+  std::unordered_map<int, std::vector<int>> cell_inst_to_elem_;
+
 protected:
   void init_temperatures() override;
 
+  void init_densities() override;
+
   void init_heat_source() override;
 
 private:
@@ -61,7 +69,11 @@ class OpenmcHeatDriver : public CoupledDriver {
 
   std::unique_ptr<SurrogateHeatDriver> heat_driver_; //!< The heat surrogate driver
 
-  int32_t n_materials_; //! Number of materials in OpenMC model
+  //! Number of OpenMC cells that are solid
+  int32_t n_solid_cells_;
+
+  //! Number of OpenMC cells that are fluid
+  int32_t n_fluid_cells_;
 };
 
 } // namespace enrico

src/openmc_heat_driver.cpp

@@ -25,6 +25,7 @@ OpenmcHeatDriver::OpenmcHeatDriver(MPI_Comm comm, pugi::xml_node node)
   init_tallies();
 
   init_temperatures();
+  init_densities();
   init_heat_source();
 }
 
@@ -52,22 +53,28 @@ void OpenmcHeatDriver::init_mappings()
   const auto& z = heat_driver_->z_;
 
   std::unordered_map<CellInstance, int> tracked;
+
+  // index for rings in solid
   int ring_index = 0;
+
+  // index for elements in fluid
+  int elem_index = 0;
+
   // TODO: Don't hardcode number of azimuthal segments
   int n_azimuthal = 4;
 
   int n_fissionable_mapped = 0;
 
+  // Establish mappings between solid regions and OpenMC cells. The center
+  // coordinate for each region in the T/H model is obtained and used to
+  // determine the OpenMC cell at that position.
   for (int i = 0; i < heat_driver_->n_pins_; ++i) {
-    // Get coordinate of pin center
     double x_center = heat_driver_->pin_centers_(i, 0);
     double y_center = heat_driver_->pin_centers_(i, 1);
 
     for (int j = 0; j < heat_driver_->n_axial_; ++j) {
-      // Get average z value
       double zavg = 0.5 * (z(j) + z(j + 1));
 
-      // Loop over radial rings
       for (int k = 0; k < heat_driver_->n_rings(); ++k) {
         double ravg;
         if (k < heat_driver_->n_fuel_rings_) {
@@ -124,36 +131,137 @@ void OpenmcHeatDriver::init_mappings()
   // capture all fissionable cells in OpenMC.
   Ensures(openmc_driver_->n_fissionable_cells_ == n_fissionable_mapped);
 
+  // set the number of solid cells before defining mappings for fluid cells
+  if (openmc_driver_->active()) {
+    n_solid_cells_ = openmc_driver_->cells_.size();
+  }
+
+  // Establish mappings between fluid regions and OpenMC cells; it is assumed that there
+  // are no azimuthal divisions in the fluid phase in the OpenMC model so that we
+  // can take a point on a 45 degree ray from the pin center. TODO: add a check to make
+  // sure that the T/H model is finer than the OpenMC model.
+  for (int i = 0; i < heat_driver_->n_pins_; ++i) {
+    double x_center = heat_driver_->pin_centers_(i, 0);
+    double y_center = heat_driver_->pin_centers_(i, 1);
+
+    for (int j = 0; j < heat_driver_->n_axial_; ++j) {
+      double zavg = 0.5 * (z(j) + z(j + 1));
+      double l = heat_driver_->pin_pitch() / std::sqrt(2.0);
+      double d = (l - heat_driver_->clad_outer_radius_) / 2.0;
+      double x = x_center + (heat_driver_->clad_outer_radius_ + d) * std::sqrt(2.0) / 2.0;
+      double y = y_center + (heat_driver_->clad_outer_radius_ + d) * std::sqrt(2.0) / 2.0;
+
+      // Determine cell instance corresponding to given fluid location
+      Position r{x, y, zavg};
+      CellInstance c{r};
+      if (tracked.find(c) == tracked.end()) {
+        openmc_driver_->cells_.push_back(c);
+        tracked[c] = openmc_driver_->cells_.size() - 1;
+
+        Ensures(!c.is_fissionable());
+      }
+
+      // Map OpenMC material to fluid element and vice versa
+      int32_t array_index = tracked[c];
+      cell_inst_to_elem_[array_index].push_back(elem_index);
+      elem_to_cell_inst_[elem_index].push_back(array_index);
+
+      ++elem_index;
+    }
+  }
+
   if (openmc_driver_->active()) {
-    n_materials_ = openmc_driver_->cells_.size();
+    n_fluid_cells_ = openmc_driver_->cells_.size() - n_solid_cells_;
   }
 }
 
 void OpenmcHeatDriver::init_tallies()
 {
   if (openmc_driver_->active()) {
-    // Build vector of material indices
+    // Build vector of material indices to construct tallies; tallies are only
+    // used in the solid regions
     std::vector<int32_t> mats;
-    for (const auto& c : openmc_driver_->cells_) {
-      mats.push_back(c.material_index_);
+    for (int c = 0; c < n_solid_cells_; ++c) {
+      const auto& cell = openmc_driver_->cells_[c];
+      mats.push_back(cell.material_index_);
     }
+
     openmc_driver_->create_tallies(mats);
   }
 }
 
+void OpenmcHeatDriver::init_densities()
+{
+  if (this->has_global_coupling_data()) {
+    std::size_t n_solid =
+      heat_driver_->n_pins_ * heat_driver_->n_axial_ * heat_driver_->n_rings();
+    std::size_t n_fluid = heat_driver_->n_pins_ * heat_driver_->n_axial_;
+    densities_.resize({n_solid + n_fluid});
+    densities_prev_.resize({n_solid + n_fluid});
+
+    if (density_ic_ == Initial::neutronics) {
+      // Loop over all of the fluid regions in the heat transfer model and set
+      // the density IC based on the densities used in the OpenMC input file. Set
+      // the initial solid density to 0.0 because it is not used in the simulation
+      // at all anyways.
+      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, ++ring_index) {
+            densities_(ring_index) = 0.0;
+          }
+        }
+      }
+
+      int elem_index = 0;
+      int solid_offset = ring_index;
+      for (int i = 0; i < heat_driver_->n_pins_; ++i) {
+        for (int j = 0; j < heat_driver_->n_axial_; ++j) {
+
+          double rho_avg = 0.0;
+          double total_vol = 0.0;
+
+          for (const auto& idx : elem_to_cell_inst_[elem_index++]) {
+            const auto& c = openmc_driver_->cells_[idx];
+            double vol = c.volume_;
+
+            total_vol += vol;
+            rho_avg += c.get_density() * vol;
+          }
+
+          densities_(elem_index + solid_offset) = rho_avg / total_vol;
+        }
+      }
+
+      std::copy(densities_.begin(), densities_.end(), densities_prev_.begin());
+    }
+
+    if (density_ic_ == Initial::heat) {
+      throw std::runtime_error{
+        "Density initial conditions from surrogate heat-fluids "
+        "solver not supported."};
+    }
+  }
+}
+
 void OpenmcHeatDriver::init_temperatures()
 {
   if (this->has_global_coupling_data()) {
-    std::size_t n =
+    std::size_t n_solid =
       heat_driver_->n_pins_ * heat_driver_->n_axial_ * heat_driver_->n_rings();
-    temperatures_.resize({n});
-    temperatures_prev_.resize({n});
+    std::size_t n_fluid = heat_driver_->n_pins_ * heat_driver_->n_axial_;
+    temperatures_.resize({n_solid + n_fluid});
+    temperatures_prev_.resize({n_solid + n_fluid});
 
     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.
+      // should be a volume average of the OpenMC cell temperatures. Then, loop over all
+      // axial fluid cells for each pin and set the fluid temperature IC based on the
+      // fluid temperature in the OpenMC model. No more than one OpenMC cell should correspond
+      // to a given axial fluid cell, so no averaging is needed (though we retain the averaging
+      // syntax below: TODO - update elem_to_cell_inst_ to permit multiple OpenMC cells per axial.
 
       // 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
@@ -165,11 +273,14 @@ void OpenmcHeatDriver::init_temperatures()
       // the azimuthal segments in the OpenMC model are all the same (i.e. no initial
       // azimuthal dependence).
 
+      int index = 0;
+
+      // set temperatures for solid cells
       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];
+            const auto& cell_instances = ring_to_cell_inst_[ring_index++];
 
             double T_avg = 0.0;
             double total_vol = 0.0;
@@ -182,16 +293,35 @@ void OpenmcHeatDriver::init_temperatures()
               T_avg += c.get_temperature() * vol;
             }
 
-            T_avg /= total_vol;
+            temperatures_[index] = T_avg / total_vol;
+            ++index;
+          }
+        }
+      }
 
-            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;
+      // set temperatures for fluid cells
+      int axial_index = 0;
+      for (int i = 0; i < heat_driver_->n_pins_; ++i) {
+        for (int j = 0; j < heat_driver_->n_axial_; ++j) {
+          const auto& cell_instances = elem_to_cell_inst_[axial_index++];
 
-            ++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;
           }
+
+          temperatures_[index] = T_avg / total_vol;
+          ++index;
         }
       }
+
+      std::copy(temperatures_.begin(), temperatures_.end(), temperatures_prev_.begin());
     }
 
     if (temperature_ic_ == Initial::heat) {
@@ -204,8 +334,8 @@ void OpenmcHeatDriver::init_temperatures()
 
 void OpenmcHeatDriver::init_heat_source()
 {
-  heat_source_ = xt::empty<double>({n_materials_});
-  heat_source_prev_ = xt::empty<double>({n_materials_});
+  heat_source_ = xt::empty<double>({n_solid_cells_});
+  heat_source_prev_ = xt::empty<double>({n_solid_cells_});
 }
 
 void OpenmcHeatDriver::set_heat_source()
@@ -238,42 +368,73 @@ void OpenmcHeatDriver::set_heat_source()
   }
 }
 
+void OpenmcHeatDriver::update_density()
+{
+  if (this->has_global_coupling_data()) {
+    std::copy(densities_.begin(), densities_.end(), densities_prev_.begin());
+  }
+
+  densities_ = heat_driver_->density();
+  int solid_offset = heat_driver_->n_pins_ * heat_driver_->n_rings() * heat_driver_->n_axial_;
+
+  for (int i = 0; i < n_fluid_cells_; ++i) {
+    int index = i + n_solid_cells_;
+    const auto& c = openmc_driver_->cells_[index];
+    const auto& elems = cell_inst_to_elem_[index];
+
+    double average_rho = 0.0;
+    double total_vol = 0.0;
+    for (int elem_index : elems) {
+      double vol = heat_driver_->z_(elem_index + 1) - heat_driver_->z_(elem_index);
+      average_rho += densities_(elem_index + solid_offset);
+      total_vol += vol;
+    }
+
+    c.set_density(average_rho / total_vol);
+  }
+}
+
 void OpenmcHeatDriver::set_temperature()
 {
   const auto& r_fuel = heat_driver_->r_grid_fuel_;
   const auto& r_clad = heat_driver_->r_grid_clad_;
 
-  // For each OpenMC material, volume average temperatures and set
-  for (int i = 0; i < openmc_driver_->cells_.size(); ++i) {
-    // Get cell instance
+  // For each OpenMC solid cell, volume average solid temperatures and set
+  // by grabbing the T/H regions corresponding to the cell
+  for (int i = 0; i < n_solid_cells_; ++i) {
     const auto& c = openmc_driver_->cells_[i];
-
-    // Get rings corresponding to this cell instance
     const auto& rings = cell_inst_to_ring_[i];
 
-    // Get volume-average temperature for this material
     double average_temp = 0.0;
     double total_vol = 0.0;
     for (int ring_index : rings) {
-      // Use difference in r**2 as a proxy for volume. This is only used for
-      // averaging, so the absolute value doesn't matter
-      int j = ring_index % heat_driver_->n_rings();
-      double vol;
-      if (j < heat_driver_->n_fuel_rings_) {
-        vol = r_fuel(j + 1) * r_fuel(j + 1) - r_fuel(j) * r_fuel(j);
-      } else {
-        j -= heat_driver_->n_fuel_rings_;
-        vol = r_clad(j + 1) * r_clad(j + 1) - r_clad(j) * r_clad(j);
-      }
-
+      double vol = heat_driver_->solid_areas_(ring_index);
       average_temp += temperatures_(ring_index) * vol;
       total_vol += vol;
     }
-    average_temp /= total_vol;
 
-    // Set temperature for cell instance
+    average_temp /= total_vol;
     c.set_temperature(average_temp);
   }
+
+  // For each OpenMC fluid cell, set the fluid temperature by grabbing
+  // the T/H region corresponding to the cell
+  int solid_offset = heat_driver_->n_pins_ * heat_driver_->n_rings() * heat_driver_->n_axial_;
+  for (int i = 0; i < n_fluid_cells_; ++i) {
+    int index = n_solid_cells_ + i;
+    const auto& c = openmc_driver_->cells_[index];
+    const auto& elems = cell_inst_to_elem_[index];
+
+    double average_temp = 0.0;
+    double total_vol = 0.0;
+    for (int elem_index : elems) {
+      double vol = heat_driver_->z_(elem_index + 1) - heat_driver_->z_(elem_index);
+      average_temp += temperatures_(elem_index + solid_offset) * vol;
+      total_vol += vol;
+    }
+
+    c.set_temperature(average_temp / total_vol);
+  }
 }
 
 } // namespace enrico