Address Patrick's second round of review comments.

Author: nuclearkevin

include/openmc/cell.h

@@ -223,6 +223,12 @@ class Cell {
   //! \return Density multiplier [-]
   double density_mult(int32_t instance = -1) const;
 
+  //! Get the density of a cell instance in g/cm3
+  //! \param[in] instance Instance index. If -1 is given, the density
+  //! for the first instance is returned.
+  //! \return Density multiplier g/cm3
+  double density(int32_t instance = -1) const;
+
   //! Set the temperature of a cell instance
   //! \param[in] T Temperature in [K]
   //! \param[in] instance Instance index. If -1 is given, the temperature for
@@ -233,15 +239,14 @@ class Cell {
   void set_temperature(
     double T, int32_t instance = -1, bool set_contained = false);
 
-  //! Set the density multiplier of a cell instance
-  //! \param[in] rho Density multiplier [-]
-  //! \param[in] instance Instance index. If -1 is given, the density multiplier
-  //! for
-  //!   all instances is set.
+  //! Set the density of a cell instance
+  //! \param[in] rho Density [g/cm3]
+  //! \param[in] instance Instance index. If -1 is given, the density
+  //! for all instances is set.
   //! \param[in] set_contained If this cell is not filled with a material,
   //!   collect all contained cells with material fills and set their
-  //!   density multipliers.
-  void set_density_mult(
+  //!   densities.
+  void set_density(
     double rho, int32_t instance = -1, bool set_contained = false);
 
   int32_t n_instances() const;

openmc/cell.py

@@ -74,7 +74,8 @@ class Cell(IDManagerMixin):
         to give each distributed cell instance a unique temperature.
     density : float or iterable of float
         Density of the cell in g/cm3. Multiple densities can be given to give
-        each distributed cell instance a unique density.
+        each distributed cell instance a unique density. Densities set here will
+        override the density set on materials used to fill the cell.
     translation : Iterable of float
         If the cell is filled with a universe, this array specifies a vector
         that is used to translate (shift) the universe.

src/cell.cpp

@@ -112,6 +112,21 @@ double Cell::density_mult(int32_t instance) const
   }
 }
 
+double Cell::density(int32_t instance) const
+{
+  const int32_t mat_index = material(instance);
+  if (mat_index == MATERIAL_VOID)
+    return 0.0;
+
+  if (instance >= 0) {
+    double rho =
+      rho_mult_.size() == 1 ? rho_mult_.at(0) : rho_mult_.at(instance);
+    return rho * model::materials[mat_index]->density_gpcc();
+  } else {
+    return rho_mult_[0];
+  }
+}
+
 void Cell::set_temperature(double T, int32_t instance, bool set_contained)
 {
   if (settings::temperature_method == TemperatureMethod::INTERPOLATION) {
@@ -162,7 +177,7 @@ void Cell::set_temperature(double T, int32_t instance, bool set_contained)
   }
 }
 
-void Cell::set_density_mult(double rho, int32_t instance, bool set_contained)
+void Cell::set_density(double rho, int32_t instance, bool set_contained)
 {
   if (type_ != Fill::MATERIAL && !set_contained) {
     fatal_error(
@@ -172,17 +187,21 @@ void Cell::set_density_mult(double rho, int32_t instance, bool set_contained)
   }
 
   if (type_ == Fill::MATERIAL) {
+    const int32_t mat_index = material(instance);
+    if (mat_index == MATERIAL_VOID)
+      return;
+
     if (instance >= 0) {
       // If density multiplier vector is not big enough, resize it first
       if (rho_mult_.size() != n_instances())
         rho_mult_.resize(n_instances(), rho_mult_[0]);
 
       // Set density multiplier for the corresponding instance
-      rho_mult_.at(instance) = rho;
+      rho_mult_.at(instance) = rho / model::materials[mat_index]->density_gpcc();
     } else {
       // Set density multiplier for all instances
       for (auto& Rho_ : rho_mult_) {
-        Rho_ = rho;
+        Rho_ = rho / model::materials[mat_index]->density_gpcc();
       }
     }
   } else {
@@ -192,7 +211,7 @@ void Cell::set_density_mult(double rho, int32_t instance, bool set_contained)
       assert(cell->type_ == Fill::MATERIAL);
       auto& instances = entry.second;
       for (auto instance : instances) {
-        cell->set_density_mult(rho, instance);
+        cell->set_density(rho, instance);
       }
     }
   }
@@ -1243,18 +1262,8 @@ extern "C" int openmc_cell_set_density(
 
   int32_t instance_index = instance ? *instance : -1;
   try {
-    if (model::cells[index]->type_ != Fill::MATERIAL) {
-      fatal_error(
-        fmt::format("Cell {}, instance {} is not filled with a material.",
-          model::cells[index]->id_, instance_index));
-    }
-
-    int32_t mat_index = model::cells[index]->material(instance_index);
-    if (mat_index != MATERIAL_VOID) {
-      model::cells[index]->set_density_mult(
-        rho / model::materials[mat_index]->density_gpcc(), instance_index,
-        set_contained);
-    }
+    model::cells[index]->set_density(
+      rho, instance_index, set_contained);
   } catch (const std::exception& e) {
     set_errmsg(e.what());
     return OPENMC_E_UNASSIGNED;

tests/regression_tests/cpp_driver/driver.cpp

@@ -65,8 +65,8 @@ int main(int argc, char** argv)
   // model)
   model::cells[model::cell_map[4]]->set_temperature(400.0, 3, true);
 
-  // set a larger density multiplier for another lattice cell
-  model::cells[model::cell_map[4]]->set_density_mult(2.0, 2, true);
+  // set the density of another lattice cell to 2
+  model::cells[model::cell_map[4]]->set_density(2.0, 2, true);
 
   // the summary file will be used to check that
   // temperatures were set correctly so clear

tests/regression_tests/cpp_driver/results_true.dat

@@ -1,13 +1,13 @@
 k-combined:
-1.964706E+00 2.868566E-02
+1.953962E+00 1.828426E-02
 tally 1:
-8.017467E+01
-7.150161E+02
-2.300257E+01
-5.904224E+01
-7.992670E+01
-7.135847E+02
-1.831039E+02
-3.734941E+03
-1.831039E+02
-3.734941E+03
+9.607953E+01
+1.031898E+03
+2.853683E+01
+9.085469E+01
+9.745011E+01
+1.058928E+03
+2.220665E+02
+5.496813E+03
+2.220665E+02
+5.496813E+03