Merge pull request #1569 from pshriwise/get_all_cells

Add option for setting temperature of contained cells

include/openmc/cell.h

@@ -43,6 +43,8 @@ constexpr int32_t OP_UNION        {std::numeric_limits<int32_t>::max() - 4};
 //==============================================================================
 
 class Cell;
+class ParentCell;
+class CellInstance;
 class Universe;
 class UniversePartitioner;
 
@@ -74,7 +76,6 @@ public:
 };
 
 //==============================================================================
-//! A geometry primitive that links surfaces, universes, and materials
 //==============================================================================
 
 class Cell {
@@ -135,7 +136,10 @@ public:
   //! \param[in] T Temperature in [K]
   //! \param[in] instance Instance index. If -1 is given, the temperature for
   //!   all instances is set.
-  void set_temperature(double T, int32_t instance = -1);
+  //! \param[in] set_contained If this cell is not filled with a material,
+  //!   collect all contained cells with material fills and set their
+  //!   temperatures.
+  void set_temperature(double T, int32_t instance = -1, bool set_contained = false);
 
   //! Get the name of a cell
   //! \return Cell name
@@ -145,6 +149,17 @@ public:
   //! \param[in] name Cell name
   void set_name(const std::string& name) { name_ = name; };
 
+  //! Get all cell instances contained by this cell
+  //! \return Map with cell indexes as keys and instances as values
+  std::unordered_map<int32_t, std::vector<int32_t>>
+  get_contained_cells() const;
+
+protected:
+  void
+  get_contained_cells_inner(std::unordered_map<int32_t, std::vector<int32_t>>& contained_cells,
+                            std::vector<ParentCell>& parent_cells) const;
+
+public:
   //----------------------------------------------------------------------------
   // Data members
 
@@ -191,6 +206,11 @@ public:
   std::vector<int32_t> offset_;  //!< Distribcell offset table
 };
 
+struct CellInstanceItem {
+  int32_t index {-1};        //! Index into global cells array
+  int     lattice_indx{-1};  //! Flat index value of the lattice cell
+};
+
 //==============================================================================
 
 class CSGCell : public Cell
@@ -289,6 +309,16 @@ private:
   std::vector<std::vector<int32_t>> partitions_;
 };
 
+
+//==============================================================================
+//! Define a containing (parent) cell
+//==============================================================================
+
+struct ParentCell {
+  gsl::index cell_index;
+  gsl::index lattice_index;
+};
+
 //==============================================================================
 //! Define an instance of a particular cell
 //==============================================================================

include/openmc/lattice.h

@@ -129,6 +129,13 @@ public:
   //!  cell found in the geometry tree under this lattice tile.
   virtual int32_t& offset(int map, const int i_xyz[3]) = 0;
 
+  //! \brief Get the distribcell offset for a lattice tile.
+  //! \param The map index for the target cell.
+  //! \param indx The index for a lattice tile.
+  //! \return Distribcell offset i.e. the largest instance number for the target
+  //!  cell found in the geometry tree for this lattice index.
+  virtual int32_t offset(int map, int indx) const = 0;
+
   //! \brief Convert an array index to a useful human-readable string.
   //! \param indx The index for a lattice tile.
   //! \return A string representing the lattice tile.
@@ -220,6 +227,8 @@ public:
 
   int32_t& offset(int map, const int i_xyz[3]);
 
+  int32_t offset(int map, int indx) const;
+
   std::string index_to_string(int indx) const;
 
   void to_hdf5_inner(hid_t group_id) const;
@@ -262,6 +271,8 @@ public:
 
   int32_t& offset(int map, const int i_xyz[3]);
 
+  int32_t offset(int map, int indx) const;
+
   std::string index_to_string(int indx) const;
 
   void to_hdf5_inner(hid_t group_id) const;

src/cell.cpp

@@ -5,8 +5,8 @@
 #include <cctype>
 #include <cmath>
 #include <iterator>
-#include <sstream>
 #include <set>
+#include <sstream>
 #include <string>
 
 #include <fmt/core.h>
@@ -240,7 +240,7 @@ Cell::temperature(int32_t instance) const
 }
 
 void
-Cell::set_temperature(double T, int32_t instance)
+Cell::set_temperature(double T, int32_t instance, bool set_contained)
 {
   if (settings::temperature_method == TemperatureMethod::INTERPOLATION) {
     if (T < data::temperature_min) {
@@ -252,16 +252,33 @@ Cell::set_temperature(double T, int32_t instance)
     }
   }
 
-  if (instance >= 0) {
-    // If temperature vector is not big enough, resize it first
-    if (sqrtkT_.size() != n_instances_) sqrtkT_.resize(n_instances_, sqrtkT_[0]);
+  if (type_ == Fill::MATERIAL) {
+    if (instance >= 0) {
+      // If temperature vector is not big enough, resize it first
+      if (sqrtkT_.size() != n_instances_) sqrtkT_.resize(n_instances_, sqrtkT_[0]);
 
-    // Set temperature for the corresponding instance
-    sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T);
+      // Set temperature for the corresponding instance
+      sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T);
+    } else {
+      // Set temperature for all instances
+      for (auto& T_ : sqrtkT_) {
+        T_ = std::sqrt(K_BOLTZMANN * T);
+      }
+    }
   } else {
-    // Set temperature for all instances
-    for (auto& T_ : sqrtkT_) {
-      T_ = std::sqrt(K_BOLTZMANN * T);
+    if (!set_contained) {
+      throw std::runtime_error{fmt::format("Attempted to set the temperature of cell {} "
+                                           "which is not filled by a material.", id_)};
+    }
+
+    auto contained_cells = this->get_contained_cells();
+    for (const auto& entry : contained_cells) {
+      auto& cell = model::cells[entry.first];
+      Expects(cell->type_ == Fill::MATERIAL);
+      auto& instances =  entry.second;
+      for (auto instance : instances) {
+        cell->set_temperature(T, instance);
+      }
     }
   }
 }
@@ -1159,6 +1176,67 @@ openmc_cell_set_name(int32_t index, const char* name) {
 }
 
 
+std::unordered_map<int32_t, std::vector<int32_t>>
+Cell::get_contained_cells() const {
+  std::unordered_map<int32_t, std::vector<int32_t>> contained_cells;
+  std::vector<ParentCell> parent_cells;
+
+  // if this cell is filled w/ a material, it contains no other cells
+  if (type_ != Fill::MATERIAL) {
+    this->get_contained_cells_inner(contained_cells, parent_cells);
+  }
+
+  return contained_cells;
+}
+
+//! Get all cells within this cell
+void
+Cell::get_contained_cells_inner(std::unordered_map<int32_t, std::vector<int32_t>>& contained_cells,
+                                std::vector<ParentCell>& parent_cells) const
+{
+
+  // filled by material, determine instance based on parent cells
+  if (type_ == Fill::MATERIAL) {
+    int instance = 0;
+    if (this->distribcell_index_ >= 0) {
+      for (auto& parent_cell : parent_cells) {
+        auto& cell = openmc::model::cells[parent_cell.cell_index];
+        if (cell->type_ == Fill::UNIVERSE) {
+          instance += cell->offset_[distribcell_index_];
+        } else if (cell->type_ == Fill::LATTICE) {
+          auto& lattice = model::lattices[cell->fill_];
+          instance += lattice->offset(this->distribcell_index_, parent_cell.lattice_index);
+        }
+      }
+    }
+    // add entry to contained cells
+    contained_cells[model::cell_map[id_]].push_back(instance);
+  // filled with universe, add the containing cell to the parent cells
+  // and recurse
+  } else if (type_ == Fill::UNIVERSE) {
+    parent_cells.push_back({model::cell_map[id_], -1});
+    auto& univ = model::universes[fill_];
+    for(auto cell_index : univ->cells_) {
+      auto& cell = model::cells[cell_index];
+      cell->get_contained_cells_inner(contained_cells, parent_cells);
+    }
+    parent_cells.pop_back();
+  // filled with a lattice, visit each universe in the lattice
+  // with a recursive call to collect the cell instances
+  } else if (type_ == Fill::LATTICE) {
+    auto& lattice = model::lattices[fill_];
+    for (auto i = lattice->begin(); i != lattice->end(); ++i) {
+      auto& univ = model::universes[*i];
+      parent_cells.push_back({model::cell_map[id_], i.indx_});
+      for (auto cell_index : univ->cells_) {
+        auto& cell = model::cells[cell_index];
+        cell->get_contained_cells_inner(contained_cells, parent_cells);
+      }
+      parent_cells.pop_back();
+    }
+  }
+}
+
 //! Return the index in the cells array of a cell with a given ID
 extern "C" int
 openmc_get_cell_index(int32_t id, int32_t* index)

src/lattice.cpp

@@ -342,6 +342,14 @@ RectLattice::offset(int map, const int i_xyz[3])
 
 //==============================================================================
 
+int32_t
+RectLattice::offset(int map, int indx) const
+{
+  return offsets_[nx*ny*nz*map + indx];
+}
+
+//==============================================================================
+
 std::string
 RectLattice::index_to_string(int indx) const
 {
@@ -973,6 +981,15 @@ HexLattice::offset(int map, const int i_xyz[3])
   return offsets_[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]];
 }
 
+int32_t
+HexLattice::offset(int map, int indx) const
+{
+  int nx {2*n_rings_ - 1};
+  int ny {2*n_rings_ - 1};
+  int nz {n_axial_};
+  return offsets_[nx*ny*nz*map + indx];
+}
+
 //==============================================================================
 
 std::string

tests/regression_tests/cpp_driver/driver.cpp

@@ -1,6 +1,8 @@
 
 #include "openmc/capi.h"
 #include "openmc/cell.h"
+#include "openmc/geometry.h"
+#include "openmc/summary.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/filter_cell.h"
 #include "openmc/tallies/tally.h"
@@ -29,6 +31,25 @@ int main(int argc, char** argv) {
   tally->set_filters(filters);
   tally->set_scores({"flux"});
 
+  // set the temperature of the cell containing
+  // the lattice
+  auto& root_univ = openmc::model::universes[openmc::model::root_universe];
+  auto& lattice_cell = openmc::model::cells[root_univ->cells_[0]];
+  lattice_cell->set_temperature(300, 1, true);
+
+  // check that material-filled cells return no contained cells
+  for (auto& cell : openmc::model::cells) {
+    if (cell->type_ == Fill::MATERIAL) {
+      auto contained_cells = cell->get_contained_cells();
+      assert(contained_cells.empty());
+    }
+  }
+
+  // the summary file will be used to check that
+  // temperatures were set correctly so clear
+  // error output can be provided
+  openmc::write_summary();
+
   openmc_run();
   openmc_finalize();
   return 0;

tests/regression_tests/cpp_driver/inputs_true.dat

@@ -2,8 +2,22 @@
 <geometry>
   <cell id="1" material="1" region="-1" universe="1" />
   <cell id="2" material="2" region="1 -2" universe="1" />
+  <cell id="3" material="3" region="2" universe="1" />
+  <cell fill="2" id="4" region="3 -4 5 -6" universe="3" />
+  <lattice id="2">
+    <pitch>4.0 4.0</pitch>
+    <dimension>2 2</dimension>
+    <lower_left>-4.0 -4.0</lower_left>
+    <universes>
+1 1 
+1 1 </universes>
+  </lattice>
   <surface coeffs="0.0 0.0 1.5" id="1" type="z-cylinder" />
-  <surface boundary="reflective" coeffs="0.0 0.0 3.0" id="2" type="z-cylinder" />
+  <surface coeffs="0.0 0.0 1.7" id="2" type="z-cylinder" />
+  <surface boundary="reflective" coeffs="-4.0" id="3" name="minimum x" type="x-plane" />
+  <surface boundary="reflective" coeffs="4.0" id="4" name="maximum x" type="x-plane" />
+  <surface boundary="reflective" coeffs="-4.0" id="5" name="minimum y" type="y-plane" />
+  <surface boundary="reflective" coeffs="4.0" id="6" name="maximum y" type="y-plane" />
 </geometry>
 <?xml version='1.0' encoding='utf-8'?>
 <materials>
@@ -11,7 +25,11 @@
     <density units="g/cc" value="11" />
     <nuclide ao="1.0" name="U235" />
   </material>
-  <material id="2" name="water">
+  <material id="2" name="cladding">
+    <density units="g/cc" value="6.44" />
+    <nuclide ao="1.0" name="Zr90" />
+  </material>
+  <material id="3" name="water">
     <density units="g/cc" value="1.0" />
     <nuclide ao="2.0" name="H1" />
     <nuclide ao="1.0" name="O16" />

tests/regression_tests/cpp_driver/results_true.dat

@@ -1,7 +1,11 @@
 k-combined:
-1.857752E+00 2.922425E-02
+1.934690E+00 2.593158E-02
 tally 1:
-5.337194E+01
-3.209877E+02
-1.621671E+02
-2.939588E+03
+8.771138E+01
+8.618680E+02
+2.643695E+01
+7.811923E+01
+8.917373E+01
+8.880318E+02
+2.033221E+02
+4.619301E+03

tests/regression_tests/cpp_driver/test.py

@@ -4,6 +4,7 @@ import shutil
 import subprocess
 import textwrap
 
+from numpy.testing import assert_allclose
 import openmc
 import pytest
 
@@ -46,6 +47,7 @@ def cpp_driver(request):
     finally:
         # Remove local build directory when test is complete
         shutil.rmtree('build')
+        os.remove('CMakeLists.txt')
 
 
 @pytest.fixture
@@ -57,22 +59,39 @@ def model():
     u235.add_nuclide('U235', 1.0, 'ao')
     u235.set_density('g/cc', 11)
 
+    zirc = openmc.Material(name='cladding')
+    zirc.add_nuclide('Zr90', 1.0)
+    zirc.set_density('g/cc', 6.44)
+
     water = openmc.Material(name="water")
     water.add_nuclide('H1', 2.0, 'ao')
     water.add_nuclide('O16', 1.0, 'ao')
     water.set_density('g/cc', 1.0)
 
-    mats = openmc.Materials([u235, water])
+    mats = openmc.Materials([u235, zirc, water])
     model.materials = mats
 
     # geometry
     fuel_or = openmc.ZCylinder(r=1.5)
-    coolant_or = openmc.ZCylinder(r=3.0, boundary_type='reflective')
+    cladding_or = openmc.ZCylinder(r=1.7)
 
     fuel = openmc.Cell(fill=u235, region=-fuel_or)
-    coolant = openmc.Cell(fill=water, region=+fuel_or & -coolant_or)
+    cladding = openmc.Cell(fill=zirc, region=+fuel_or & -cladding_or)
+    moderator = openmc.Cell(fill=water, region=+cladding_or)
 
-    model.geometry = openmc.Geometry([fuel, coolant])
+    pincell_univ = openmc.Universe(cells=[fuel, cladding, moderator])
+
+    # lattice
+    lattice = openmc.RectLattice()
+    lattice.pitch = (4.0, 4.0)
+    lattice.lower_left = (-4.0, -4.0)
+    lattice.universes = [[pincell_univ, pincell_univ], [pincell_univ, pincell_univ]]
+    lattice_region = openmc.model.rectangular_prism(8.0,
+                                                    8.0,
+                                                    boundary_type='reflective')
+    lattice_cell = openmc.Cell(fill=lattice, region=lattice_region)
+
+    model.geometry = openmc.Geometry([lattice_cell])
 
     model.settings.particles = 100
     model.settings.batches = 10
@@ -97,6 +116,22 @@ class ExternalDriverTestHarness(PyAPITestHarness):
             openmc.run(openmc_exec=self.executable,
                        event_based=config['event'])
 
+    def _compare_results(self):
+        super()._compare_results()
+
+        # load the summary file
+        summary = openmc.Summary('summary.h5')
+
+        # get the summary cells
+        cells = summary.geometry.get_all_cells()
+
+        # for the 2 by 2 lattice, each cell should have 4
+        # temperature values set to 300 K
+        for cell in cells.values():
+            if isinstance(cell.fill, openmc.Material):
+                assert len(cell.temperature) == 4
+                assert_allclose(cell.temperature, 300.0)
+
 
 def test_cpp_driver(cpp_driver, model):
     harness = ExternalDriverTestHarness(cpp_driver, 'statepoint.10.h5', model)

tests/regression_tests/source_dlopen/test.py

@@ -39,6 +39,7 @@ def compile_source(request):
 
     # Remove local build directory when test is complete
     shutil.rmtree('build')
+    os.remove('CMakeLists.txt')
 
 
 @pytest.fixture