Add density initialization to OpenmcNekDriver. Refs enrico-dev#17

include/enrico/coupled_driver.h

@@ -103,6 +103,10 @@ class CoupledDriver {
   //! temperatures in the neutronics input file.
   Initial temperature_ic_{Initial::neutronics};
 
+  //! Where to obtain the density initial condition from. Defaults to the densities
+  //! in the neutronics input file.
+  Initial density_ic_{Initial::neutronics};
+
 protected:
   //! Initialize current and previous Picard temperature fields
   virtual void init_temperatures() {}

src/cell_instance.cpp

@@ -56,7 +56,7 @@ void CellInstance::set_density(double rho) const
 double CellInstance::get_density() const
 {
   double rho;
-  err_chk(openmc_material_get_density(material_index));
+  err_chk(openmc_material_get_density(material_index_));
   return rho;
 }
 

src/coupled_driver.cpp

@@ -34,6 +34,18 @@ CoupledDriver::CoupledDriver(MPI_Comm comm, pugi::xml_node node)
     }
   }
 
+  if (node.child("density_ic")) {
+    auto s = std::string{node.child_value("density_ic")};
+
+    if (s == "neutronics") {
+      density_ic_ = Initial::neutronics;
+    } else if (s == "heat") {
+      density_ic_ = Initial::heat;
+    } else {
+      throw std::runtime_error{"Invalid value for <density_ic>"};
+    }
+  }
+
   Expects(power_ > 0);
   Expects(max_timesteps_ >= 0);
   Expects(max_picard_iter_ >= 0);

src/openmc_nek_driver.cpp

@@ -275,8 +275,35 @@ void OpenmcNekDriver::init_densities()
   if (this->has_global_coupling_data()) {
     densities_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
     densities_prev_.resize({gsl::narrow<std::size_t>(n_global_elem_)});
+
+    if (density_ic_ == Initial::neutronics) {
+      // Loop over the OpenMC cells, then loop over the global Nek elements
+      // corresponding to that cell and assign the OpenMC cell density to
+      // the correct index in the densities_ array. This mapping assumes that
+      // each Nek element is fully contained within an OpenMC cell, i.e. Nek
+      // elements are not split between multiple OpenMC cells.
+      for (const auto& c: openmc_driver_->cells_) {
+        const auto& global_elems = mat_to_elems_.at(c.material_index_);
+
+        for (int elem: global_elems) {
+          double rho = c.get_density();
+          densities_[elem] = rho;
+          densities_prev_[elem] = rho;
+        }
+      }
+    }
+    else if (density_ic_ == Initial::heat) {
+      // Use whatever density is in Nek's internal arrays, either from a restart
+      // file or from a useric fortran routine
+      update_density();
+
+      // update_density() begins by saving densities_ to densities_prev_, and
+      // then changes densities_. We need to save densities_ here to densities_prev_
+      // manually because init_densities() initializes both densities_ and
+      // densities_prev_
+      std::copy(densities_.begin(), densities_.end(), densities_prev_.begin());
+    }
   }
-  update_density();
 }
 
 void OpenmcNekDriver::init_elem_fluid_mask()