Merge pull request #1150 from pshriwise/uwuw

UWUW Workflow for DAGMC

docs/source/io_formats/settings.rst

@@ -88,6 +88,14 @@ you care. This element has the following attributes/sub-elements:
 
     *Default*: 0.0
 
+--------------------------------
+``<dagmc>`` Element
+--------------------------------
+
+When the DAGMC mode is enabled, the OpenMC geometry will be read from the file
+``dagmc.h5m``. If a :ref:`geometry.xml <io_geometry>` file is present with
+``dagmc`` set to ``true``, it will be ignored.
+
 --------------------------------
 ``<electron_treatment>`` Element
 --------------------------------

docs/source/usersguide/geometry.rst

@@ -398,3 +398,32 @@ if needed, lattices, the last step is to create an instance of
 
 .. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry
 .. _quadratic surfaces: http://en.wikipedia.org/wiki/Quadric
+
+--------------------------
+Using CAD-based Geometry
+--------------------------
+
+OpenMC relies on the Direct Accelerated Geometry Monte Carlo toolkit (`DAGMC
+<https://svalinn.github.io/DAGMC/>`_) to represent CAD-based geometry in a
+surface mesh format. A DAGMC run can be enabled in OpenMC by setting the
+``dagmc`` property to ``True`` in the model Settings either via the Python
+:class:`openmc.settings` Python class::
+
+  settings = openmc.Settings()
+  settings.dagmc = True
+
+or in the :ref:`settings.xml <io_settings>` file::
+
+  <dagmc>true</dagmc>
+
+With ``dagmc`` set to true, OpenMC will load the DAGMC model (from a local file
+named ``dagmc.h5m``) when initializing a simulation. If a `geometry.xml
+<../io_formats/geometry.html>`_ is present as well, it will be ignored.
+
+  **Note:** DAGMC geometries used in OpenMC are currently required to be clean,
+  meaning that all surfaces have been `imprinted and merged
+  <https://svalinn.github.io/DAGMC/usersguide/trelis_workflow.html>`_
+  successfully and that the model is `watertight
+  <https://svalinn.github.io/DAGMC/usersguide/tools.html#make-watertight>`_. Future
+  implementations of DAGMC geometry will support small volume overlaps and
+  un-merged surfaces.

docs/source/usersguide/install.rst

@@ -159,14 +159,25 @@ Prerequisites
           sudo apt install mpich libmpich-dev
           sudo apt install openmpi-bin libopenmpi-dev
 
+    * DAGMC_ toolkit for simulation using CAD-based geometries
+
+      OpenMC supports particle tracking in CAD-based geometries via the Direct
+      Accelerated Geometry Monte Carlo (DAGMC) toolkit (`installation
+      instructions
+      <https://svalinn.github.io/DAGMC/install/dag_multiple.html>`_). For use in
+      OpenMC, only the ``MOAB_DIR`` and ``BUILD_TALLY`` variables need to be
+      specified in the CMake configuration step.
+    |
     * git_ version control software for obtaining source code
 
+
 .. _gfortran: http://gcc.gnu.org/wiki/GFortran
 .. _gcc: https://gcc.gnu.org/
 .. _CMake: http://www.cmake.org
 .. _OpenMPI: http://www.open-mpi.org
 .. _MPICH: http://www.mpich.org
 .. _HDF5: https://www.hdfgroup.org/solutions/hdf5/
+.. _DAGMC: https://svalinn.github.io/DAGMC/index.html
 
 Obtaining the Source
 --------------------
@@ -236,6 +247,12 @@ openmp
   Enables shared-memory parallelism using the OpenMP API. The Fortran compiler
   being used must support OpenMP. (Default: on)
 
+dagmc
+  Enables use of CAD-based DAGMC_ geometries. Please see the note about DAGMC in
+  the optional dependencies list for more information on this feature. The
+  installation directory for DAGMC should also be defined as `DAGMC_ROOT` in the
+  CMake configuration command. (Default: off)
+
 coverage
   Compile and link code instrumented for coverage analysis. This is typically
   used in conjunction with gcov_.

include/openmc/capi.h

@@ -7,6 +7,7 @@
 
 #ifdef __cplusplus
 #include "openmc/bank.h"
+
 extern "C" {
   int openmc_fission_bank(openmc::Bank** ptr, int64_t* n);
   int openmc_source_bank(openmc::Bank** ptr, int64_t* n);

include/openmc/dagmc.h

@@ -2,6 +2,10 @@
 #ifndef OPENMC_DAGMC_H
 #define OPENMC_DAGMC_H
 
+namespace openmc {
+extern "C" const bool dagmc_enabled;
+}
+
 #ifdef DAGMC
 
 #include "DagMC.hpp"
@@ -10,10 +14,6 @@
 
 namespace openmc {
 
-//==============================================================================
-// Global variables
-//==============================================================================
-
 namespace model {
 
 extern moab::DagMC* DAG;
@@ -26,15 +26,11 @@ extern moab::DagMC* DAG;
 
 extern "C" void load_dagmc_geometry();
 extern "C" void free_memory_dagmc();
+extern "C" pugi::xml_document* read_uwuw_materials();
+bool get_uwuw_materials_xml(std::string& s);
 
 } // namespace openmc
 
 #endif // DAGMC
 
 #endif // OPENMC_DAGMC_H
-
-#ifdef DAGMC
-extern "C" constexpr bool dagmc_enabled = true;
-#else
-extern "C" constexpr bool dagmc_enabled = false;
-#endif

src/cross_sections.cpp

@@ -2,6 +2,9 @@
 
 #include "openmc/constants.h"
 #include "openmc/container_util.h"
+#ifdef DAGMC
+#include "openmc/dagmc.h"
+#endif
 #include "openmc/error.h"
 #include "openmc/file_utils.h"
 #include "openmc/settings.h"
@@ -81,15 +84,30 @@ extern "C" void read_mg_cross_sections_header();
 
 void read_cross_sections_xml()
 {
-  // Check if materials.xml exists
+  pugi::xml_document doc;
   std::string filename = settings::path_input + "materials.xml";
+#ifdef DAGMC
+  std::string s;
+  bool found_uwuw_mats = false;
+  if (settings::dagmc) {
+    found_uwuw_mats = get_uwuw_materials_xml(s);
+  }
+
+  if (found_uwuw_mats) {
+    // if we found uwuw materials, load those
+    doc.load_file(s.c_str());
+  } else {
+#endif
+  // Check if materials.xml exists
   if (!file_exists(filename)) {
     fatal_error("Material XML file '" + filename + "' does not exist.");
   }
-
   // Parse materials.xml file
-  pugi::xml_document doc;
   doc.load_file(filename.c_str());
+#ifdef DAGMC
+  }
+#endif
+
   auto root = doc.document_element();
 
   // Find cross_sections.xml file -- the first place to look is the

src/dagmc.cpp

@@ -1,16 +1,38 @@
 #include "openmc/dagmc.h"
 
 #include "openmc/cell.h"
+#include "openmc/constants.h"
 #include "openmc/error.h"
+#include "openmc/file_utils.h"
 #include "openmc/string_utils.h"
 #include "openmc/settings.h"
 #include "openmc/geometry.h"
 
+#ifdef DAGMC
+
+#include "uwuw.hpp"
+#include "dagmcmetadata.hpp"
+
+#endif
+
 #include <string>
 #include <sstream>
 #include <algorithm>
+#include <fstream>
+
+namespace openmc {
 
 #ifdef DAGMC
+const bool dagmc_enabled = true;
+#else
+const bool dagmc_enabled = false;
+#endif
+
+}
+
+#ifdef DAGMC
+
+const std::string DAGMC_FILENAME = "dagmc.h5m";
 
 namespace openmc {
 
@@ -20,28 +42,96 @@ moab::DagMC* DAG;
 
 } // namespace model
 
+
+bool get_uwuw_materials_xml(std::string& s) {
+  UWUW uwuw(DAGMC_FILENAME.c_str());
+
+  std::stringstream ss;
+  bool uwuw_mats_present = false;
+  if (uwuw.material_library.size() != 0) {
+    uwuw_mats_present = true;
+    // write header
+    ss << "<?xml version=\"1.0\"?>\n";
+    ss << "<materials>\n";
+    const auto& mat_lib = uwuw.material_library;
+    // write materials
+    for (auto mat : mat_lib) { ss << mat.second.openmc("atom"); }
+    // write footer
+    ss << "</materials>";
+    s = ss.str();
+  }
+
+  return uwuw_mats_present;
+}
+
+pugi::xml_document* read_uwuw_materials() {
+  pugi::xml_document* doc = nullptr;
+
+  std::string s;
+  bool found_uwuw_mats = get_uwuw_materials_xml(s);
+  if (found_uwuw_mats) {
+    doc = new pugi::xml_document();
+    pugi::xml_parse_result result = doc->load_string(s.c_str());
+  }
+  return doc;
+}
+
+bool write_uwuw_materials_xml() {
+  std::string s;
+  bool found_uwuw_mats = get_uwuw_materials_xml(s);
+    // if there is a material library in the file
+  if (found_uwuw_mats) {
+    // write a material.xml file
+    std::ofstream mats_xml("materials.xml");
+    mats_xml << s;
+    mats_xml.close();
+  }
+
+  return found_uwuw_mats;
+}
+
 void load_dagmc_geometry()
 {
   if (!model::DAG) {
     model::DAG = new moab::DagMC();
   }
 
-  int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC
+  /// Materials \\\
 
-  moab::ErrorCode rval = model::DAG->load_file("dagmc.h5m");
+  // create uwuw instance
+  UWUW uwuw(DAGMC_FILENAME.c_str());
+
+  // check for uwuw material definitions
+  bool using_uwuw = !uwuw.material_library.empty();
+
+  // notify user if UWUW materials are going to be used
+  if (using_uwuw) {
+    std::cout << "Found UWUW Materials in the DAGMC geometry file.\n";
+  }
+
+  int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC runs
+
+  // load the DAGMC geometry
+  moab::ErrorCode rval = model::DAG->load_file(DAGMC_FILENAME.c_str());
   MB_CHK_ERR_CONT(rval);
 
+  // initialize acceleration data structures
   rval = model::DAG->init_OBBTree();
   MB_CHK_ERR_CONT(rval);
 
-  std::vector<std::string> prop_keywords;
-  prop_keywords.push_back("mat");
-  prop_keywords.push_back("boundary");
-
-  std::map<std::string, std::string> ph;
-  model::DAG->parse_properties(prop_keywords, ph, ":");
+  // parse model metadata
+  dagmcMetaData DMD(model::DAG);
+  if (using_uwuw) {
+    DMD.load_property_data();
+  }
+  std::vector<std::string> keywords {"temp", "mat", "density", "boundary"};
+  std::map<std::string, std::string> dum;
+  std::string delimiters = ":/";
+  rval = model::DAG->parse_properties(keywords, dum, delimiters.c_str());
   MB_CHK_ERR_CONT(rval);
 
+  /// Cells (Volumes) \\\
+
   // initialize cell objects
   model::n_cells = model::DAG->num_entities(3);
 
@@ -69,35 +159,80 @@ void load_dagmc_geometry()
       model::universes[it->second]->cells_.push_back(i);
     }
 
+    // check for temperature assignment
+    std::string temp_value;
+    if (model::DAG->has_prop(vol_handle, "temp")) {
+      rval = model::DAG->prop_value(vol_handle, "temp", temp_value);
+      MB_CHK_ERR_CONT(rval);
+      double temp = std::stod(temp_value);
+      c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp));
+    } else {
+      c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * settings::temperature_default));
+    }
+
+    // MATERIALS
+
     if (model::DAG->is_implicit_complement(vol_handle)) {
-      // assuming implicit complement is void for now
-      c->material_.push_back(MATERIAL_VOID);
+      if (model::DAG->has_prop(vol_handle, "mat")) {
+        // if the implicit complement has been assigned a material, use it
+        std::string comp_mat = DMD.volume_material_property_data_eh[vol_handle];
+        // Note: material numbers are set by UWUW
+        int mat_number = uwuw.material_library[comp_mat].metadata["mat_number"].asInt();
+        c->material_.push_back(mat_number);
+      } else {
+        // if no material is found, the implicit complement is void
+        c->material_.push_back(MATERIAL_VOID);
+      }
       continue;
     }
 
-    if (model::DAG->has_prop(vol_handle, "mat")){
-      std::string mat_value;
+    // determine volume material assignment
+    std::string mat_value;
+    if (model::DAG->has_prop(vol_handle, "mat")) {
       rval = model::DAG->prop_value(vol_handle, "mat", mat_value);
       MB_CHK_ERR_CONT(rval);
-      to_lower(mat_value);
-
-      if (mat_value == "void" || mat_value == "vacuum") {
-        c->material_.push_back(MATERIAL_VOID);
-      } else {
-        c->material_.push_back(std::stoi(mat_value));
-      }
     } else {
       std::stringstream err_msg;
       err_msg << "Volume " << c->id_ << " has no material assignment.";
       fatal_error(err_msg.str());
     }
+
+    std::string cmp_str = mat_value;
+    to_lower(cmp_str);
+    // material void checks
+    if (cmp_str.find("void") != std::string::npos   ||
+        cmp_str.find("vacuum") != std::string::npos ||
+        cmp_str.find("graveyard") != std::string::npos) {
+      c->material_.push_back(MATERIAL_VOID);
+    } else {
+      if (using_uwuw) {
+        // lookup material in uwuw if the were present
+        std::string uwuw_mat = DMD.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[uwuw_mat].metadata["mat_number"].asInt();
+          c->material_.push_back(mat_number);
+        } else {
+          std::stringstream err_msg;
+          err_msg << "Material with value " << mat_value << " not found ";
+          err_msg << "in the UWUW material library";
+          fatal_error(err_msg);
+        }
+      } else {
+        // if not using UWUW materials, we'll find this material
+        // later in the materials.xml
+        c->material_.push_back(std::stoi(mat_value));
+      }
+    }
   }
 
-  // Allocate the cell overlap count if necessary.
+  // allocate the cell overlap count if necessary
   if (settings::check_overlaps) {
     model::overlap_check_count.resize(model::cells.size(), 0);
   }
 
+  /// Surfaces \\\
+
   // initialize surface objects
   int n_surfaces = model::DAG->num_entities(2);
   model::surfaces.resize(n_surfaces);
@@ -110,8 +245,9 @@ void load_dagmc_geometry()
     s->id_ = model::DAG->id_by_index(2, i+1);
     s->dagmc_ptr_ = model::DAG;
 
+    // set BCs
+    std::string bc_value;
     if (model::DAG->has_prop(surf_handle, "boundary")) {
-      std::string bc_value;
       rval = model::DAG->prop_value(surf_handle, "boundary", bc_value);
       MB_CHK_ERR_CONT(rval);
       to_lower(bc_value);
@@ -130,7 +266,8 @@ void load_dagmc_geometry()
                 << "\" specified on surface " << s->id_;
         fatal_error(err_msg);
       }
-    } else {   // if no BC property is found, set to transmit
+    } else {
+      // if no condition is found, set to transmit
       s->bc_ = BC_TRANSMIT;
     }
 
@@ -147,5 +284,6 @@ void free_memory_dagmc()
   delete model::DAG;
 }
 
+
 }
 #endif

src/dagmc_header.F90

@@ -13,6 +13,11 @@ module dagmc_header
     subroutine free_memory_dagmc() bind(C)
     end subroutine free_memory_dagmc
 
+    function read_uwuw_materials() result(doc) bind(C)
+      import C_PTR
+      type(C_PTR) :: doc
+    end function read_uwuw_materials
+
   end interface
 
 end module dagmc_header

src/input_xml.F90

@@ -520,18 +520,29 @@ contains
     ! Display output message
     call write_message("Reading materials XML file...", 5)
 
-    ! Check if materials.xml exists
-    filename = trim(path_input) // "materials.xml"
-    inquire(FILE=filename, EXIST=file_exists)
-    if (.not. file_exists) then
-      call fatal_error("Material XML file '" // trim(filename) // "' does not &
-           &exist!")
+    doc % ptr = C_NULL_PTR
+
+#ifdef DAGMC
+    if (dagmc) then
+      doc % ptr = read_uwuw_materials()
+    end if
+#endif
+
+    if (.not. c_associated(doc % ptr)) then
+      ! Check if materials.xml exists
+      filename = trim(path_input) // "materials.xml"
+      inquire(FILE=filename, EXIST=file_exists)
+      if (.not. file_exists) then
+        call fatal_error("Material XML file '" // trim(filename) // "' does not &
+             &exist!")
     end if
 
     ! Parse materials.xml file
     call doc % load_file(filename)
-    root = doc % document_element()
 
+    end if
+
+    root = doc % document_element()
     call read_materials(root % ptr)
 
     ! Get pointer to list of XML <material>

src/settings.cpp

@@ -405,13 +405,6 @@ void read_settings_xml()
 #endif
   }
 
-#ifdef _OPENMP
-  if (dagmc && omp_get_max_threads() > 1) {
-    warning("Forcing number of threads to 1 for DAGMC simulation.");
-    omp_set_num_threads(1);
-  }
-#endif
-
   // ==========================================================================
   // EXTERNAL SOURCE
 

src/tracking.F90

@@ -499,7 +499,7 @@ contains
       p % coord(1) % cell = i_cell-1 ! decrement for C++ indexing
       p % cell_instance = 1
       p % material = cells(i_cell) % material(1)
-      p % sqrtKT = cells(i_cell) % sqrtKT(1)
+      p % sqrtKT = cells(i_cell) % sqrtKT(0)
       return
     end if
 #endif

tests/regression_tests/dagmc/results_true.dat

@@ -1,5 +1,5 @@
 k-combined:
-1.115067E+00 5.423808E-02
+1.028803E+00 3.340602E-02
 tally 1:
-8.543144E+00
-1.530584E+01
+8.346847E+00
+1.453407E+01

tests/regression_tests/dagmc/test.py

@@ -38,8 +38,8 @@ def test_dagmc():
     water = openmc.Material()
     water.add_nuclide('H1', 2.0, 'ao')
     water.add_nuclide('O16', 1.0, 'ao')
-    water.add_s_alpha_beta('c_H_in_H2O')
     water.set_density('g/cc', 1.0)
+    water.add_s_alpha_beta('c_H_in_H2O')
     water.id = 41
 
     mats = openmc.Materials([u235, water])

tests/regression_tests/uwuw/inputs_true.dat

@@ -0,0 +1,23 @@
+<?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>
+  <dagmc>true</dagmc>
+</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/uwuw/results_true.dat

@@ -0,0 +1 @@
+../dagmc/results_true.dat

tests/regression_tests/uwuw/test.py

@@ -0,0 +1,41 @@
+import openmc
+import openmc.capi
+from openmc.stats import Box
+from openmc.material import Materials
+
+import pytest
+from tests.testing_harness import PyAPITestHarness
+
+pytestmark = pytest.mark.skipif(
+    not openmc.capi._dagmc_enabled(),
+    reason="DAGMC CAD geometry is not enabled.")
+
+class UWUWTest(PyAPITestHarness):
+
+    def _build_inputs(self):
+        model = openmc.model.Model()
+
+        # settings
+        model.settings.batches = 5
+        model.settings.inactive = 0
+        model.settings.particles = 100
+
+        source = openmc.Source(space=Box([-4, -4, -4],
+                                         [ 4,  4,  4]))
+        model.settings.source = source
+
+        model.settings.dagmc = True
+
+        model.settings.export_to_xml()
+
+        # tally
+        tally = openmc.Tally()
+        tally.scores = ['total']
+        tally.filters = [openmc.CellFilter(1)]
+        model.tallies = [tally]
+
+        model.tallies.export_to_xml()
+
+def test_uwuw():
+    harness = UWUWTest('statepoint.5.h5')
+    harness.main()