Read/write densities instead of density multipliers in properties.h5

Author: nuclearkevin

docs/source/io_formats/properties.rst

@@ -25,9 +25,7 @@ The current version of the properties file format is 1.0.
 **/geometry/cells/cell <uid>/**
 
 :Datasets: - **temperature** (*double[]*) -- Temperature of the cell in [K].
-           - **density_mult** (*double[]*) -- Unitless density multipliers for
-             the cell. The cell density is equal to the density multiplier
-             times the density of the material filling the cell.
+           - **density** (*double[]*) -- Density of the cell in [g/cm3].
 
 **/materials/**
 

docs/source/io_formats/summary.rst

@@ -38,9 +38,7 @@ The current version of the summary file format is 6.0.
              is an array if the cell uses distributed materials, otherwise it is
              a scalar.
            - **temperature** (*double[]*) -- Temperature of the cell in Kelvin.
-           - **density_mult** (*double[]*) -- Unitless density multipliers for the cell.
-             The cell density is equal to the density multiplier times the density of the
-             material filling the cell.
+           - **density** (*double[]*) -- Density of the cell in [g/cm3].
            - **translation** (*double[3]*) -- Translation applied to the fill
              universe. This dataset is present only if fill_type is set to
              'universe'.

openmc/model/model.py

@@ -648,19 +648,19 @@ def import_properties(self, filename: PathLike):
                         else:
                             lib_cell.set_temperature(temperature[0])
 
-                    density_mult = group['density_mult'][()]
-                    mat_density = cell.fill.get_mass_density()
-                    if density_mult.size > 1:
-                        cell.density = [mat_density * m for m in density_mult]
-                    else:
-                        cell.density = density_mult * mat_density
-                    if self.is_initialized:
-                        lib_cell = openmc.lib.cells[cell_id]
-                        if density_mult.size > 1:
-                            for i, rho_mult in enumerate(density_mult):
-                                lib_cell.set_density(rho_mult * mat_density, i)
-                        else:
-                            lib_cell.set_density(density_mult[0] * mat_density)
+                    if group['density']:
+                      density = group['density'][()]
+                      if density.size > 1:
+                          cell.density = [rho for rho in density]
+                      else:
+                          cell.density = density
+                      if self.is_initialized:
+                          lib_cell = openmc.lib.cells[cell_id]
+                          if density.size > 1:
+                              for i, rho in enumerate(density):
+                                  lib_cell.set_density(rho, i)
+                          else:
+                              lib_cell.set_density(density[0])
 
             # Make sure number of materials matches
             mats_group = fh['materials']

src/cell.cpp

@@ -114,6 +114,9 @@ double Cell::density_mult(int32_t instance) const
 
 double Cell::density(int32_t instance) const
 {
+  if (instance >= material_.size() && material_.size() != 1)
+    fatal_error(fmt::format("Instances and materials don't match! {} {}", instance, material_.size()));
+
   const int32_t mat_index = material(instance);
   if (mat_index == MATERIAL_VOID)
     return 0.0;
@@ -123,7 +126,7 @@ double Cell::density(int32_t instance) const
       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];
+    return rho_mult_[0] * model::materials[mat_index]->density_gpcc();
   }
 }
 
@@ -229,8 +232,14 @@ void Cell::export_properties_hdf5(hid_t group) const
     temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);
   write_dataset(cell_group, "temperature", temps);
 
-  // Write density multipliers for one or more cell instances
-  write_dataset(cell_group, "density_mult", rho_mult_);
+  // Write density for one or more cell instances
+  if (type_ == Fill::MATERIAL && material_.size() > 0) {
+    vector<double> rho;
+    for (int32_t i = 0; i < rho_mult_.size(); ++i)
+      rho.push_back(density(i));
+
+    write_dataset(cell_group, "density", rho);
+  }
 
   close_group(cell_group);
 }
@@ -258,17 +267,22 @@ void Cell::import_properties_hdf5(hid_t group)
     this->set_temperature(temps[i], i);
   }
 
-  // Read density multipliers
-  if (object_exists(cell_group, "density_mult")) {
-    read_dataset(cell_group, "density_mult", rho_mult_);
+  // Read densities
+  if (object_exists(cell_group, "density")) {
+    vector<double> rho;
+    read_dataset(cell_group, "density", rho);
 
-    // Ensure number of density multipliers makes sense
-    auto n_rho = rho_mult_.size();
+    // Ensure number of densities makes sense
+    auto n_rho = rho.size();
     if (n_rho > 1 && n_rho != n_instances()) {
-      fatal_error(fmt::format("Number of density multipliers for cell {} "
+      fatal_error(fmt::format("Number of densities for cell {} "
                               "doesn't match number of instances",
         id_));
     }
+
+    // Set densities.
+    for (int32_t i = 0; i < n_rho; ++i)
+      set_density(rho[i], i);
   }
 
   close_group(cell_group);

tests/unit_tests/test_model.py

@@ -254,8 +254,8 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm):
     # Change cell fuel temperature, density, material density and export properties
     cell = openmc.lib.cells[1]
     cell.set_temperature(600.0)
-    cell.set_density(10.0)
     cell.fill.set_density(5.0, 'g/cm3')
+    cell.set_density(10.0)
     openmc.lib.export_properties(output=False)
 
     # Import properties to existing model