Merge pull request openmc-dev#1851 from helen-brooks/external-moab-in…

…terface

Support for external MOAB instance via new constructor

include/openmc/mesh.h

@@ -366,6 +366,7 @@ class MOABMesh : public UnstructuredMesh {
   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 @@ class MOABMesh : public UnstructuredMesh {
 
   // 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 @@ class MOABMesh : public UnstructuredMesh {
   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 @@ class Tally {
 
   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;
+}