Merge pull request #2036 from helen-brooks/refactor-dagmc-universes

Support external DAGMC instance

include/openmc/cell.h

@@ -300,5 +300,8 @@ struct CellInstanceHash {
 
 void read_cells(pugi::xml_node node);
 
+//!  Add cells to universes
+void populate_universes();
+
 } // namespace openmc
 #endif // OPENMC_CELL_H

include/openmc/dagmc.h

@@ -22,6 +22,7 @@ void check_dagmc_root_univ();
 #ifdef DAGMC
 
 #include "DagMC.hpp"
+#include "dagmcmetadata.hpp"
 
 #include "openmc/cell.h"
 #include "openmc/particle.h"
@@ -87,6 +88,11 @@ public:
   explicit DAGUniverse(const std::string& filename, bool auto_geom_ids = false,
     bool auto_mat_ids = false);
 
+  //! Alternative DAGMC universe constructor for external DAGMC instance
+  explicit DAGUniverse(std::shared_ptr<moab::DagMC> external_dagmc_ptr,
+    const std::string& filename = "", bool auto_geom_ids = false,
+    bool auto_mat_ids = false);
+
   //! Initialize the DAGMC accel. data structures, indices, material
   //! assignments, etc.
   void initialize();
@@ -139,10 +145,16 @@ public:
   bool has_graveyard() const { return has_graveyard_; }
 
 private:
+  void set_id();        //!< Deduce the universe id from model::universes
+  void init_dagmc();    //!< Create and initialise DAGMC pointer
+  void init_metadata(); //!< Create and initialise dagmcMetaData pointer
+  void init_geometry(); //!< Create cells and surfaces from DAGMC entities
+
   std::string
     filename_; //!< Name of the DAGMC file used to create this universe
   std::shared_ptr<UWUW>
-    uwuw_;                   //!< Pointer to the UWUW instance for this universe
+    uwuw_; //!< Pointer to the UWUW instance for this universe
+  std::unique_ptr<dagmcMetaData> dmd_ptr; //! Pointer to DAGMC metadata object
   bool adjust_geometry_ids_; //!< Indicates whether or not to automatically
                              //!< generate new cell and surface IDs for the
                              //!< universe

include/openmc/material.h

@@ -113,6 +113,9 @@ public:
   //! \param[in] units Units of density
   void set_density(double density, gsl::cstring_span units);
 
+  //! Set temperature of the material
+  void set_temperature(double temperature) { temperature_ = temperature; };
+
   //! Get nuclides in material
   //! \return Indices into the global nuclides vector
   gsl::span<const int> nuclides() const

src/cell.cpp

@@ -846,6 +846,20 @@ void read_cells(pugi::xml_node node)
 
   read_dagmc_universes(node);
 
+  populate_universes();
+
+  // Allocate the cell overlap count if necessary.
+  if (settings::check_overlaps) {
+    model::overlap_check_count.resize(model::cells.size(), 0);
+  }
+
+  if (model::cells.size() == 0) {
+    fatal_error("No cells were found in the geometry.xml file");
+  }
+}
+
+void populate_universes()
+{
   // Populate the Universe vector and map.
   for (int i = 0; i < model::cells.size(); i++) {
     int32_t uid = model::cells[i]->universe_;
@@ -860,15 +874,6 @@ void read_cells(pugi::xml_node node)
     }
   }
   model::universes.shrink_to_fit();
-
-  // Allocate the cell overlap count if necessary.
-  if (settings::check_overlaps) {
-    model::overlap_check_count.resize(model::cells.size(), 0);
-  }
-
-  if (model::cells.size() == 0) {
-    fatal_error("No cells were found in the geometry.xml file");
-  }
 }
 
 //==============================================================================

src/dagmc.cpp

@@ -12,7 +12,6 @@
 #include "openmc/string_utils.h"
 
 #ifdef DAGMC
-#include "dagmcmetadata.hpp"
 #include "uwuw.hpp"
 #endif
 #include <fmt/core.h>
@@ -74,6 +73,23 @@ DAGUniverse::DAGUniverse(
   const std::string& filename, bool auto_geom_ids, bool auto_mat_ids)
   : filename_(filename), adjust_geometry_ids_(auto_geom_ids),
     adjust_material_ids_(auto_mat_ids)
+{
+  set_id();
+  initialize();
+}
+
+DAGUniverse::DAGUniverse(std::shared_ptr<moab::DagMC> dagmc_ptr,
+  const std::string& filename, bool auto_geom_ids, bool auto_mat_ids)
+  : dagmc_instance_(dagmc_ptr), filename_(filename),
+    adjust_geometry_ids_(auto_geom_ids), adjust_material_ids_(auto_mat_ids)
+{
+  set_id();
+  init_metadata();
+  read_uwuw_materials();
+  init_geometry();
+}
+
+void DAGUniverse::set_id()
 {
   // determine the next universe id
   int32_t next_univ_id = 0;
@@ -85,48 +101,27 @@ DAGUniverse::DAGUniverse(
 
   // set the universe id
   id_ = next_univ_id;
-
-  initialize();
 }
 
 void DAGUniverse::initialize()
 {
   geom_type() = GeometryType::DAG;
 
-  // determine the next cell id
-  int32_t next_cell_id = 0;
-  for (const auto& c : model::cells) {
-    if (c->id_ > next_cell_id)
-      next_cell_id = c->id_;
-  }
-  cell_idx_offset_ = model::cells.size();
-  next_cell_id++;
+  init_dagmc();
 
-  // determine the next surface id
-  int32_t next_surf_id = 0;
-  for (const auto& s : model::surfaces) {
-    if (s->id_ > next_surf_id)
-      next_surf_id = s->id_;
-  }
-  surf_idx_offset_ = model::surfaces.size();
-  next_surf_id++;
+  init_metadata();
+
+  read_uwuw_materials();
+
+  init_geometry();
+}
+
+void DAGUniverse::init_dagmc()
+{
 
   // create a new DAGMC instance
   dagmc_instance_ = std::make_shared<moab::DagMC>();
 
-  // --- Materials ---
-
-  // read any UWUW materials from the file
-  read_uwuw_materials();
-
-  // check for uwuw material definitions
-  bool using_uwuw = uses_uwuw();
-
-  // notify user if UWUW materials are going to be used
-  if (using_uwuw) {
-    write_message("Found UWUW Materials in the DAGMC geometry file.", 6);
-  }
-
   // load the DAGMC geometry
   filename_ = settings::path_input + filename_;
   if (!file_exists(filename_)) {
@@ -138,18 +133,35 @@ void DAGUniverse::initialize()
   // initialize acceleration data structures
   rval = dagmc_instance_->init_OBBTree();
   MB_CHK_ERR_CONT(rval);
+}
 
+void DAGUniverse::init_metadata()
+{
   // parse model metadata
-  dagmcMetaData DMD(dagmc_instance_.get(), false, false);
-  DMD.load_property_data();
+  dmd_ptr =
+    std::make_unique<dagmcMetaData>(dagmc_instance_.get(), false, false);
+  dmd_ptr->load_property_data();
 
   std::vector<std::string> keywords {"temp"};
   std::map<std::string, std::string> dum;
   std::string delimiters = ":/";
+  moab::ErrorCode rval;
   rval = dagmc_instance_->parse_properties(keywords, dum, delimiters.c_str());
   MB_CHK_ERR_CONT(rval);
+}
 
-  // --- Cells (Volumes) ---
+void DAGUniverse::init_geometry()
+{
+  moab::ErrorCode rval;
+
+  // determine the next cell id
+  int32_t next_cell_id = 0;
+  for (const auto& c : model::cells) {
+    if (c->id_ > next_cell_id)
+      next_cell_id = c->id_;
+  }
+  cell_idx_offset_ = model::cells.size();
+  next_cell_id++;
 
   // initialize cell objects
   int n_cells = dagmc_instance_->num_entities(3);
@@ -178,7 +190,7 @@ void DAGUniverse::initialize()
     // --- Materials ---
 
     // determine volume material assignment
-    std::string mat_str = DMD.get_volume_property("material", vol_handle);
+    std::string mat_str = dmd_ptr->get_volume_property("material", vol_handle);
 
     if (mat_str.empty()) {
       fatal_error(fmt::format("Volume {} has no material assignment.", c->id_));
@@ -194,9 +206,10 @@ void DAGUniverse::initialize()
     if (mat_str == "void" || mat_str == "vacuum" || mat_str == "graveyard") {
       c->material_.push_back(MATERIAL_VOID);
     } else {
-      if (using_uwuw) {
+      if (uses_uwuw()) {
         // lookup material in uwuw if present
-        std::string uwuw_mat = DMD.volume_material_property_data_eh[vol_handle];
+        std::string uwuw_mat =
+          dmd_ptr->volume_material_property_data_eh[vol_handle];
         if (uwuw_->material_library.count(uwuw_mat) != 0) {
           // Note: material numbers are set by UWUW
           int mat_number = uwuw_->material_library.get_material(uwuw_mat)
@@ -246,7 +259,14 @@ void DAGUniverse::initialize()
 
   has_graveyard_ = graveyard;
 
-  // --- Surfaces ---
+  // determine the next surface id
+  int32_t next_surf_id = 0;
+  for (const auto& s : model::surfaces) {
+    if (s->id_ > next_surf_id)
+      next_surf_id = s->id_;
+  }
+  surf_idx_offset_ = model::surfaces.size();
+  next_surf_id++;
 
   // initialize surface objects
   int n_surfaces = dagmc_instance_->num_entities(2);
@@ -259,7 +279,8 @@ void DAGUniverse::initialize()
                                   : dagmc_instance_->id_by_index(2, i + 1);
 
     // set BCs
-    std::string bc_value = DMD.get_surface_property("boundary", surf_handle);
+    std::string bc_value =
+      dmd_ptr->get_surface_property("boundary", surf_handle);
     to_lower(bc_value);
     if (bc_value.empty() || bc_value == "transmit" ||
         bc_value == "transmission") {
@@ -398,7 +419,7 @@ void DAGUniverse::to_hdf5(hid_t universes_group) const
 
 bool DAGUniverse::uses_uwuw() const
 {
-  return !uwuw_->material_library.empty();
+  return uwuw_ && !uwuw_->material_library.empty();
 }
 
 std::string DAGUniverse::get_uwuw_materials_xml() const
@@ -487,33 +508,32 @@ void DAGUniverse::legacy_assign_material(
 
 void DAGUniverse::read_uwuw_materials()
 {
-
-  int32_t next_material_id = 0;
-  for (const auto& m : model::materials) {
-    next_material_id = std::max(m->id_, next_material_id);
-  }
-  next_material_id++;
+  // If no filename was provided, don't read UWUW materials
+  if (filename_ == "")
+    return;
 
   uwuw_ = std::make_shared<UWUW>(filename_.c_str());
-  const auto& mat_lib = uwuw_->material_library;
-  if (mat_lib.size() == 0)
+
+  if (!uses_uwuw())
     return;
 
+  // Notify user if UWUW materials are going to be used
+  write_message("Found UWUW Materials in the DAGMC geometry file.", 6);
+
   // if we're using automatic IDs, update the UWUW material metadata
   if (adjust_material_ids_) {
+    int32_t next_material_id = 0;
+    for (const auto& m : model::materials) {
+      next_material_id = std::max(m->id_, next_material_id);
+    }
+    next_material_id++;
+
     for (auto& mat : uwuw_->material_library) {
       mat.second->metadata["mat_number"] = next_material_id++;
     }
   }
 
-  std::stringstream ss;
-  ss << "<?xml version=\"1.0\"?>\n";
-  ss << "<materials>\n";
-  for (auto mat : mat_lib) {
-    ss << mat.second->openmc("atom");
-  }
-  ss << "</materials>";
-  std::string mat_xml_string = ss.str();
+  std::string mat_xml_string = get_uwuw_materials_xml();
 
   // create a pugi XML document from this string
   pugi::xml_document doc;

tests/regression_tests/dagmc/external/dagmc.h5m

@@ -0,0 +1 @@
+../legacy/dagmc.h5m

tests/regression_tests/dagmc/external/inputs_true.dat

@@ -0,0 +1,40 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <dagmc_universe filename="dagmc.h5m" id="1" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material depletable="true" id="40" name="no-void fuel">
+    <density units="g/cc" value="11" />
+    <nuclide ao="1.0" name="U235" />
+  </material>
+  <material id="41" name="water">
+    <density units="g/cc" value="1.0" />
+    <nuclide ao="2.0" name="H1" />
+    <nuclide ao="1.0" name="O16" />
+    <sab name="c_H_in_H2O" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>eigenvalue</run_mode>
+  <particles>100</particles>
+  <batches>5</batches>
+  <inactive>0</inactive>
+  <source strength="1.0">
+    <space type="box">
+      <parameters>-4 -4 -4 4 4 4</parameters>
+    </space>
+  </source>
+  <temperature_default>293</temperature_default>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <filter id="1" type="cell">
+    <bins>1</bins>
+  </filter>
+  <tally id="1">
+    <filters>1</filters>
+    <scores>total</scores>
+  </tally>
+</tallies>

tests/regression_tests/dagmc/external/main.cpp

@@ -0,0 +1,94 @@
+#include "openmc/capi.h"
+#include "openmc/cross_sections.h"
+#include "openmc/dagmc.h"
+#include "openmc/error.h"
+#include "openmc/geometry.h"
+#include "openmc/geometry_aux.h"
+#include "openmc/material.h"
+#include "openmc/nuclide.h"
+#include <iostream>
+
+int main(int argc, char* argv[])
+{
+  using namespace openmc;
+  int openmc_err;
+
+  // Initialise OpenMC
+  openmc_err = openmc_init(argc, argv, nullptr);
+  if (openmc_err == -1) {
+    // This happens for the -h and -v flags
+    return EXIT_SUCCESS;
+  } else if (openmc_err) {
+    fatal_error(openmc_err_msg);
+  }
+
+  // Create DAGMC ptr
+  std::string filename = "dagmc.h5m";
+  std::shared_ptr<moab::DagMC> dag_ptr = std::make_shared<moab::DagMC>();
+  moab::ErrorCode rval = dag_ptr->load_file(filename.c_str());
+  if (rval != moab::MB_SUCCESS) {
+    fatal_error("Failed to load file");
+  }
+
+  // Initialize acceleration data structures
+  rval = dag_ptr->init_OBBTree();
+  if (rval != moab::MB_SUCCESS) {
+    fatal_error("Failed to initialise OBB tree");
+  }
+
+  // Get rid of existing geometry
+  std::unordered_map<std::string, int> nuclide_map_copy =
+    openmc::data::nuclide_map;
+  openmc::data::nuclides.clear();
+  openmc::data::nuclide_map = nuclide_map_copy;
+  openmc::model::surfaces.clear();
+  openmc::model::surface_map.clear();
+  openmc::model::cells.clear();
+  openmc::model::cell_map.clear();
+  openmc::model::universes.clear();
+  openmc::model::universe_map.clear();
+
+  // Update materials (emulate an external program)
+  for (auto& mat_ptr : openmc::model::materials) {
+    mat_ptr->set_temperature(300);
+  }
+
+  // Create new DAGMC universe
+  openmc::model::universes.push_back(
+    std::make_unique<openmc::DAGUniverse>(dag_ptr));
+  model::universe_map[model::universes.back()->id_] =
+    model::universes.size() - 1;
+
+  // Add cells to universes
+  openmc::populate_universes();
+
+  // Set root universe
+  openmc::model::root_universe = openmc::find_root_universe();
+  openmc::check_dagmc_root_univ();
+
+  // Final geometry setup and assign temperatures
+  openmc::finalize_geometry();
+
+  // Finalize cross sections having assigned temperatures
+  openmc::finalize_cross_sections();
+
+  // Check that we correctly assigned cell temperatures with non-void fill
+  for (auto& cell_ptr : openmc::model::cells) {
+    if (cell_ptr->material_.front() != openmc::C_NONE &&
+        cell_ptr->temperature() != 300) {
+      fatal_error("Failed to set cell temperature");
+    }
+  }
+
+  // Run OpenMC
+  openmc_err = openmc_run();
+  if (openmc_err)
+    fatal_error(openmc_err_msg);
+
+  // Deallocate memory
+  openmc_err = openmc_finalize();
+  if (openmc_err)
+    fatal_error(openmc_err_msg);
+
+  return EXIT_SUCCESS;
+}

tests/regression_tests/dagmc/external/results_true.dat

@@ -0,0 +1,5 @@
+k-combined:
+8.426936E-01 5.715847E-02
+tally 1:
+8.093843E+00
+1.328829E+01

tests/regression_tests/dagmc/external/test.py

@@ -0,0 +1,128 @@
+from pathlib import Path
+import os
+import shutil
+import subprocess
+import textwrap
+
+import openmc
+import openmc.lib
+import pytest
+
+from tests.regression_tests import config
+from tests.testing_harness import PyAPITestHarness
+
+pytestmark = pytest.mark.skipif(
+    not openmc.lib._dagmc_enabled(),
+    reason="DAGMC is not enabled.")
+
+# Test that an external DAGMC instance can be passed in through the C API
+
+@pytest.fixture
+def cpp_driver(request):
+    """Compile the external source"""
+
+    # Get build directory and write CMakeLists.txt file
+    openmc_dir = Path(str(request.config.rootdir)) / 'build'
+    with open('CMakeLists.txt', 'w') as f:
+        f.write(textwrap.dedent("""
+            cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
+            project(openmc_cpp_driver CXX)
+            add_executable(main main.cpp)
+            find_package(OpenMC REQUIRED HINTS {})
+            target_link_libraries(main OpenMC::libopenmc)
+            target_compile_features(main PUBLIC cxx_std_14)
+            set(CMAKE_CXX_FLAGS "-pedantic-errors")
+            add_compile_definitions(DAGMC=1)
+            """.format(openmc_dir)))
+
+    # Create temporary build directory and change to there
+    local_builddir = Path('build')
+    local_builddir.mkdir(exist_ok=True)
+    os.chdir(local_builddir)
+
+    mpi_arg = "On" if config['mpi'] else "Off"
+
+    try:
+        # Run cmake/make to build the shared libary
+        subprocess.run(['cmake', os.path.pardir, f'-DOPENMC_USE_MPI={mpi_arg}'], check=True)
+        subprocess.run(['make'], check=True)
+        os.chdir(os.path.pardir)
+
+        yield "./build/main"
+
+    finally:
+        # Remove local build directory when test is complete
+        shutil.rmtree('build')
+        os.remove('CMakeLists.txt')
+
+@pytest.fixture
+def model():
+    model = openmc.model.Model()
+
+    # Settings
+    model.settings.batches = 5
+    model.settings.inactive = 0
+    model.settings.particles = 100
+    source_box = openmc.stats.Box([-4, -4, -4],
+                                  [ 4,  4,  4])
+    source = openmc.Source(space=source_box)
+    model.settings.source = source
+    model.settings.temperature['default'] = 293
+
+    # Geometry
+    dag_univ = openmc.DAGMCUniverse("dagmc.h5m")
+    model.geometry = openmc.Geometry(dag_univ)
+
+    # Tallies
+    tally = openmc.Tally()
+    tally.scores = ['total']
+    tally.filters = [openmc.CellFilter(1)]
+    model.tallies = [tally]
+
+    # Materials
+    u235 = openmc.Material(name="no-void fuel")
+    u235.add_nuclide('U235', 1.0, 'ao')
+    u235.set_density('g/cc', 11)
+    u235.id = 40
+    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)
+    water.add_s_alpha_beta('c_H_in_H2O')
+    water.id = 41
+    mats = openmc.Materials([u235, water])
+    model.materials = mats
+
+    return model
+
+class ExternalDAGMCTest(PyAPITestHarness):
+    def __init__(self, executable, statepoint_name, model):
+        super().__init__(statepoint_name, model)
+        self.executable = executable
+
+    def _run_openmc(self):
+        """
+        Just test if results generated with the external C++ API are
+        self-consistent with the internal python API.
+        We generate the "truth" results with the python API but
+        the main test produces the results file by running the
+        executable compiled from main.cpp. This future-proofs
+        the test - we only care that the two routes are equivalent.
+        """
+        if config['update']:
+            # Generate the results file with internal python API
+            openmc.run(openmc_exec=config['exe'], event_based=config['event'])
+        elif config['mpi']:
+            mpi_args = [config['mpiexec'], '-n', config['mpi_np']]
+            # Run main cpp executable with MPI
+            openmc.run(openmc_exec=self.executable,
+                       mpi_args=mpi_args,
+                       event_based=config['event'])
+        else:
+            # Run main cpp executable
+            openmc.run(openmc_exec=self.executable,
+                       event_based=config['event'])
+
+def test_external_dagmc(cpp_driver, model):
+    harness = ExternalDAGMCTest(cpp_driver, 'statepoint.5.h5', model)
+    harness.main()