Merge pull request #1851 from helen-brooks/external-moab-interface

Support for external MOAB instance via new constructor

include/openmc/mesh.h

@@ -366,6 +366,7 @@ public:
   MOABMesh() = default;
   MOABMesh(pugi::xml_node);
   MOABMesh(const std::string& filename);
+  MOABMesh(std::shared_ptr<moab::Interface> external_mbi);
 
   // Overridden Methods
 
@@ -412,6 +413,9 @@ private:
 
   // Methods
 
+  //! Create the MOAB interface pointer
+  void create_interface();
+
   //! Find all intersections with faces of the mesh.
   //
   //! \param[in] start Staring location
@@ -503,7 +507,7 @@ private:
   moab::Range ehs_; //!< Range of tetrahedra EntityHandle's in the mesh
   moab::EntityHandle tetset_; //!< EntitySet containing all tetrahedra
   moab::EntityHandle kdtree_root_; //!< Root of the MOAB KDTree
-  unique_ptr<moab::Interface> mbi_;         //!< MOAB instance
+  std::shared_ptr<moab::Interface> mbi_;    //!< MOAB instance
   unique_ptr<moab::AdaptiveKDTree> kdtree_; //!< MOAB KDTree instance
   vector<moab::Matrix3> baryc_data_;        //!< Barycentric data for tetrahedra
   vector<std::string> tag_names_; //!< Names of score tags added to the mesh

include/openmc/tallies/tally.h

@@ -53,6 +53,9 @@ public:
 
   void set_filters(gsl::span<Filter*> filters);
 
+  //! Given already-set filters, set the stride lengths
+  void set_strides();
+
   int32_t strides(int i) const {return strides_[i];}
 
   int32_t n_filter_bins() const {return n_filter_bins_;}

src/mesh.cpp

@@ -160,7 +160,6 @@ StructuredMesh::bin_label(int bin) const {
 //==============================================================================
 
 UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) {
-  n_dimension_ = 3;
 
   // check the mesh type
   if (check_for_node(node, "type")) {
@@ -1582,14 +1581,22 @@ MOABMesh::MOABMesh(const std::string& filename) {
   initialize();
 }
 
+MOABMesh::MOABMesh(std::shared_ptr<moab::Interface> external_mbi) {
+  mbi_ = external_mbi;
+  filename_ = "unknown (external file)";
+  this->initialize();
+}
+
 void MOABMesh::initialize() {
-  // create MOAB instance
-  mbi_ = make_unique<moab::Core>();
-  // load unstructured mesh file
-  moab::ErrorCode rval = mbi_->load_file(filename_.c_str());
-  if (rval != moab::MB_SUCCESS) {
-    fatal_error("Failed to load the unstructured mesh file: " + filename_);
-  }
+
+  // Create the MOAB interface and load data from file
+  this->create_interface();
+
+  // Initialise MOAB error code
+  moab::ErrorCode rval = moab::MB_SUCCESS;
+
+  // Set the dimension
+  n_dimension_ = 3;
 
   // set member range of tetrahedral entities
   rval = mbi_->get_entities_by_dimension(0, n_dimension_, ehs_);
@@ -1619,6 +1626,22 @@ void MOABMesh::initialize() {
   build_kdtree(ehs_);
 }
 
+void
+MOABMesh::create_interface()
+{
+  // Do not create a MOAB instance if one is already in memory
+  if (mbi_) return;
+
+  // create MOAB instance
+  mbi_ = std::make_shared<moab::Core>();
+
+  // load unstructured mesh file
+  moab::ErrorCode rval = mbi_->load_file(filename_.c_str());
+  if (rval != moab::MB_SUCCESS) {
+    fatal_error("Failed to load the unstructured mesh file: " + filename_);
+  }
+}
+
 void
 MOABMesh::build_kdtree(const moab::Range& all_tets)
 {

src/tallies/tally.cpp

@@ -366,9 +366,18 @@ Tally::set_filters(gsl::span<Filter*> filters)
     }
   }
 
+  // Set the strides.
+  set_strides();
+
+}
+
+void
+Tally::set_strides()
+{
   // Set the strides.  Filters are traversed in reverse so that the last filter
   // has the shortest stride in memory and the first filter has the longest
   // stride.
+  auto n = filters_.size();
   strides_.resize(n, 0);
   int stride = 1;
   for (int i = n-1; i >= 0; --i) {

tests/regression_tests/external_moab/inputs_true.dat

@@ -0,0 +1,73 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="1" name="fuel" region="1 -2 3 -4 5 -6" universe="1" />
+  <cell id="2" material="2" name="clad" region="(-1 | 2 | -3 | 4 | -5 | 6) (7 -8 9 -10 11 -12)" universe="1" />
+  <cell id="3" material="3" name="water" region="(-7 | 8 | -9 | 10 | -11 | 12) (13 -14 15 -16 17 -18)" universe="1" />
+  <surface coeffs="-5.0" id="1" name="minimum x" type="x-plane" />
+  <surface coeffs="5.0" id="2" name="maximum x" type="x-plane" />
+  <surface coeffs="-5.0" id="3" name="minimum y" type="y-plane" />
+  <surface coeffs="5.0" id="4" name="maximum y" type="y-plane" />
+  <surface coeffs="-5.0" id="5" name="minimum z" type="z-plane" />
+  <surface coeffs="5.0" id="6" name="maximum z" type="z-plane" />
+  <surface coeffs="-6.0" id="7" name="minimum x" type="x-plane" />
+  <surface coeffs="6.0" id="8" name="maximum x" type="x-plane" />
+  <surface coeffs="-6.0" id="9" name="minimum y" type="y-plane" />
+  <surface coeffs="6.0" id="10" name="maximum y" type="y-plane" />
+  <surface coeffs="-6.0" id="11" name="minimum z" type="z-plane" />
+  <surface coeffs="6.0" id="12" name="maximum z" type="z-plane" />
+  <surface boundary="vacuum" coeffs="-10" id="13" name="minimum x" type="x-plane" />
+  <surface boundary="vacuum" coeffs="10" id="14" name="maximum x" type="x-plane" />
+  <surface boundary="vacuum" coeffs="-10" id="15" name="minimum y" type="y-plane" />
+  <surface boundary="vacuum" coeffs="10" id="16" name="maximum y" type="y-plane" />
+  <surface boundary="vacuum" coeffs="-10" id="17" name="minimum z" type="z-plane" />
+  <surface boundary="vacuum" coeffs="10" id="18" name="maximum z" type="z-plane" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material depletable="true" id="1" name="fuel">
+    <density units="g/cc" value="4.5" />
+    <nuclide ao="1.0" name="U235" />
+  </material>
+  <material id="2" name="zircaloy">
+    <density units="g/cc" value="5.77" />
+    <nuclide ao="0.5145" name="Zr90" />
+    <nuclide ao="0.1122" name="Zr91" />
+    <nuclide ao="0.1715" name="Zr92" />
+    <nuclide ao="0.1738" name="Zr94" />
+    <nuclide ao="0.028" name="Zr96" />
+  </material>
+  <material id="3" name="water">
+    <density units="atom/b-cm" value="0.07416" />
+    <nuclide ao="2.0" name="H1" />
+    <nuclide ao="1.0" name="O16" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>fixed source</run_mode>
+  <particles>100</particles>
+  <batches>10</batches>
+  <source strength="1.0">
+    <space type="point">
+      <parameters>0.0 0.0 0.0</parameters>
+    </space>
+    <angle reference_uvw="-1.0 0.0 0.0" type="monodirectional" />
+    <energy type="discrete">
+      <parameters>15000000.0 1.0</parameters>
+    </energy>
+  </source>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <mesh id="1" library="moab" type="unstructured">
+    <filename>test_mesh_tets.h5m</filename>
+  </mesh>
+  <filter id="1" type="mesh">
+    <bins>1</bins>
+  </filter>
+  <tally id="1" name="unstructured mesh tally">
+    <filters>1</filters>
+    <scores>flux</scores>
+    <estimator>tracklength</estimator>
+  </tally>
+</tallies>

tests/regression_tests/external_moab/main.cpp

@@ -0,0 +1,98 @@
+#include "moab/Core.hpp"
+#include "openmc/capi.h"
+#include "openmc/error.h"
+#include "openmc/mesh.h"
+#include "openmc/tallies/filter_mesh.h"
+#include "openmc/tallies/tally.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 MOAB interface
+  std::shared_ptr<moab::Interface> moabPtrLocal =
+    std::make_shared<moab::Core>();
+
+  // Load unstructured mesh file
+  std::string filename = "test_mesh_tets.h5m";
+  moab::ErrorCode rval = moabPtrLocal->load_file(filename.c_str());
+  if (rval != moab::MB_SUCCESS) {
+    fatal_error("Failed to load the unstructured mesh file: " + filename);
+  } else {
+    std::cout << "Loaded external MOAB mesh from file " << filename
+              << std::endl;
+  }
+
+  // Add a new unstructured mesh to openmc using new constructor
+  model::meshes.push_back(std::make_unique<MOABMesh>(moabPtrLocal));
+
+  // Check we now have 2 copies of the shared ptr
+  if (moabPtrLocal.use_count() != 2) {
+    fatal_error("Incorrect number of MOAB shared pointers");
+  }
+
+  // Auto-assign mesh ID
+  model::meshes.back()->set_id(C_NONE);
+  int mesh_id = model::meshes.back()->id_;
+
+  // Check we now have 2 meshes and id was correctly set
+  if (model::meshes.size() != 2)
+    fatal_error("Wrong number of meshes.");
+  else if (mesh_id != 2)
+    fatal_error("Mesh ID is incorrect");
+
+  // Add a new mesh filter with auto-assigned ID
+  Filter* filter_ptr = Filter::create("mesh", C_NONE);
+
+  // Upcast pointer type
+  MeshFilter* mesh_filter = dynamic_cast<MeshFilter*>(filter_ptr);
+
+  if (!mesh_filter) {
+    fatal_error("Failed to create mesh filter");
+  }
+
+  // Pass in the index of our mesh to the filter
+  int32_t mesh_idx = model::meshes.size() - 1;
+  mesh_filter->set_mesh(mesh_idx);
+
+  // Create a tally with auto-assigned ID
+  model::tallies.push_back(make_unique<Tally>(C_NONE));
+
+  // Set tally name - matches that in test.py
+  model::tallies.back()->name_ = "external mesh tally";
+
+  // Set tally filter to our mesh filter
+  std::vector<Filter*> filters(1, filter_ptr);
+  model::tallies.back()->set_filters(filters);
+
+  // Set tally estimator
+  model::tallies.back()->estimator_ = TallyEstimator::TRACKLENGTH;
+
+  // Set tally score
+  std::vector<std::string> score_names(1, "flux");
+  model::tallies.back()->set_scores(score_names);
+
+  // 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/external_moab/test.py

@@ -0,0 +1,275 @@
+from pathlib import Path
+import os
+import shutil
+import subprocess
+from subprocess import CalledProcessError
+import textwrap
+import glob
+from itertools import product
+
+import openmc
+import openmc.lib
+import numpy as np
+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.")
+
+TETS_PER_VOXEL = 12
+
+# Test that an external moab 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)
+            set_target_properties(main PROPERTIES CXX_STANDARD
+            14 CXX_STANDARD_REQUIRED YES CXX_EXTENSIONS NO)
+            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(str(local_builddir))
+
+    if config['mpi']:
+        os.environ['CXX'] = 'mpicxx'
+
+    try:
+        print("Building driver")
+        # Run cmake/make to build the shared libary
+        subprocess.run(['cmake', os.path.pardir], 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')
+
+
+class ExternalMoabTest(PyAPITestHarness):
+    def __init__(self, executable, statepoint_name, model):
+        super().__init__(statepoint_name, model)
+        self.executable = executable
+
+    def _run_openmc(self):
+        if config['mpi']:
+            mpi_args = [config['mpiexec'], '-n', config['mpi_np']]
+            openmc.run(openmc_exec=self.executable,
+                       mpi_args=mpi_args,
+                       event_based=config['event'])
+        else:
+            openmc.run(openmc_exec=self.executable,
+                       event_based=config['event'])
+
+    # Override some methods to do nothing
+    def _get_results(self):
+        pass
+
+    def _write_results(self, results_string):
+        pass
+
+    def _overwrite_results(self):
+        pass
+
+    def _test_output_created(self):
+        pass
+
+    # Directly compare results of unstructured  mesh with internal and
+    # external moab
+    def _compare_results(self):
+
+        with openmc.StatePoint(self._sp_name) as sp:
+            # loop over the tallies and get data
+
+            ext_data = []
+            unstr_data = []
+
+            for tally in sp.tallies.values():
+
+                # Safety check that mesh filter is correct
+                if tally.contains_filter(openmc.MeshFilter):
+                    flt = tally.find_filter(openmc.MeshFilter)
+
+                    if isinstance(flt.mesh, openmc.UnstructuredMesh):
+
+                        if tally.name == "external mesh tally":
+                            ext_data = tally.get_reshaped_data(value='mean')
+
+                        elif tally.name == "unstructured mesh tally":
+                            unstr_data = tally.get_reshaped_data(value='mean')
+
+            # we expect these results to be the same to within at 8
+            # decimal places
+            decimals = 8
+            np.testing.assert_array_almost_equal(unstr_data,
+                                                 ext_data, decimals)
+
+    @staticmethod
+    def get_mesh_tally_data(tally):
+        data = tally.get_reshaped_data(value='mean')
+        std_dev = tally.get_reshaped_data(value='std_dev')
+        data.shape = (data.size, 1)
+        std_dev.shape = (std_dev.size, 1)
+        return np.sum(data, axis=1), np.sum(std_dev, axis=1)
+
+    def _cleanup(self):
+        super()._cleanup()
+        output = glob.glob('tally*.vtk')
+        for f in output:
+            if os.path.exists(f):
+                os.remove(f)
+
+
+def test_external_mesh(cpp_driver):
+
+    # Materials
+    materials = openmc.Materials()
+
+    fuel_mat = openmc.Material(name="fuel")
+    fuel_mat.add_nuclide("U235", 1.0)
+    fuel_mat.set_density('g/cc', 4.5)
+    materials.append(fuel_mat)
+
+    zirc_mat = openmc.Material(name="zircaloy")
+    zirc_mat.add_element("Zr", 1.0)
+    zirc_mat.set_density("g/cc", 5.77)
+    materials.append(zirc_mat)
+
+    water_mat = openmc.Material(name="water")
+    water_mat.add_nuclide("H1", 2.0)
+    water_mat.add_nuclide("O16", 1.0)
+    water_mat.set_density("atom/b-cm", 0.07416)
+    materials.append(water_mat)
+
+    materials.export_to_xml()
+
+    # Geometry
+    fuel_min_x = openmc.XPlane(-5.0, name="minimum x")
+    fuel_max_x = openmc.XPlane(5.0, name="maximum x")
+
+    fuel_min_y = openmc.YPlane(-5.0, name="minimum y")
+    fuel_max_y = openmc.YPlane(5.0, name="maximum y")
+
+    fuel_min_z = openmc.ZPlane(-5.0, name="minimum z")
+    fuel_max_z = openmc.ZPlane(5.0, name="maximum z")
+
+    fuel_cell = openmc.Cell(name="fuel")
+    fuel_cell.region = +fuel_min_x & -fuel_max_x & \
+                       +fuel_min_y & -fuel_max_y & \
+                       +fuel_min_z & -fuel_max_z
+    fuel_cell.fill = fuel_mat
+
+    clad_min_x = openmc.XPlane(-6.0, name="minimum x")
+    clad_max_x = openmc.XPlane(6.0, name="maximum x")
+
+    clad_min_y = openmc.YPlane(-6.0, name="minimum y")
+    clad_max_y = openmc.YPlane(6.0, name="maximum y")
+
+    clad_min_z = openmc.ZPlane(-6.0, name="minimum z")
+    clad_max_z = openmc.ZPlane(6.0, name="maximum z")
+
+    clad_cell = openmc.Cell(name="clad")
+    clad_cell.region = (-fuel_min_x | +fuel_max_x |
+                        -fuel_min_y | +fuel_max_y |
+                        -fuel_min_z | +fuel_max_z) & \
+        (+clad_min_x & -clad_max_x &
+         +clad_min_y & -clad_max_y &
+         +clad_min_z & -clad_max_z)
+    clad_cell.fill = zirc_mat
+
+    bounds = (10, 10, 10)
+
+    water_min_x = openmc.XPlane(x0=-bounds[0],
+                                name="minimum x",
+                                boundary_type='vacuum')
+    water_max_x = openmc.XPlane(x0=bounds[0],
+                                name="maximum x",
+                                boundary_type='vacuum')
+
+    water_min_y = openmc.YPlane(y0=-bounds[1],
+                                name="minimum y",
+                                boundary_type='vacuum')
+    water_max_y = openmc.YPlane(y0=bounds[1],
+                                name="maximum y",
+                                boundary_type='vacuum')
+
+    water_min_z = openmc.ZPlane(z0=-bounds[2],
+                                name="minimum z",
+                                boundary_type='vacuum')
+    water_max_z = openmc.ZPlane(z0=bounds[2],
+                                name="maximum z",
+                                boundary_type='vacuum')
+
+    water_cell = openmc.Cell(name="water")
+    water_cell.region = (-clad_min_x | +clad_max_x |
+                         -clad_min_y | +clad_max_y |
+                         -clad_min_z | +clad_max_z) & \
+        (+water_min_x & -water_max_x &
+         +water_min_y & -water_max_y &
+         +water_min_z & -water_max_z)
+    water_cell.fill = water_mat
+
+    # create a containing universe
+    geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell])
+
+    # Meshes
+    mesh_filename = "test_mesh_tets.h5m"
+
+    # Create a normal unstructured mesh to compare to
+    uscd_mesh = openmc.UnstructuredMesh(mesh_filename, 'moab')
+
+    # Create filters
+    uscd_filter = openmc.MeshFilter(mesh=uscd_mesh)
+
+    # Tallies
+    tallies = openmc.Tallies()
+    uscd_tally = openmc.Tally(name="unstructured mesh tally")
+    uscd_tally.filters = [uscd_filter]
+    uscd_tally.scores = ['flux']
+    uscd_tally.estimator = 'tracklength'
+    tallies.append(uscd_tally)
+
+    # Settings
+    settings = openmc.Settings()
+    settings.run_mode = 'fixed source'
+    settings.particles = 100
+    settings.batches = 10
+
+    # Source setup
+    space = openmc.stats.Point()
+    angle = openmc.stats.Monodirectional((-1.0, 0.0, 0.0))
+    energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0])
+    source = openmc.Source(space=space, energy=energy, angle=angle)
+    settings.source = source
+
+    model = openmc.model.Model(geometry=geometry,
+                               materials=materials,
+                               tallies=tallies,
+                               settings=settings)
+
+    harness = ExternalMoabTest(cpp_driver,
+                               'statepoint.10.h5',
+                               model)
+
+    # Run open MC and check results
+    harness.main()