Add a cpp driver test for comparing mesh tallies with internal and external moab ptrs: currently fails because feature not implemented

tests/regression_tests/external_moab/inputs_true.dat

@@ -0,0 +1,79 @@
+<?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 origin="0.0 0.0 0.0" type="spherical">
+      <r parameters="0.0 0.0" type="uniform" />
+      <theta type="discrete">
+        <parameters>0.0 1.0</parameters>
+      </theta>
+      <phi type="discrete">
+        <parameters>0.0 1.0</parameters>
+      </phi>
+    </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,50 @@
+#include "openmc/capi.h"
+#include "openmc/error.h"
+#include "openmc/mesh.h"
+#include <iostream>
+#include "moab/Core.hpp"
+
+int main(int argc, char * argv[]){
+
+  using namespace openmc;
+  int openmc_err;
+
+  // 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 mesh from file "<<filename<<std::endl;
+  }
+
+  // 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 a new unstructured mesh tally using new constructor
+
+  // Add unstructured mesh tally to openmc
+  
+  // 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,281 @@
+from pathlib import Path
+import os
+import shutil
+import subprocess
+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
+
+from subprocess import CalledProcessError
+
+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")
+
+            """.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])
+
+    ### Tallies ###
+
+    # 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)
+
+    # Create 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
+    r = openmc.stats.Uniform(a=0.0, b=0.0)
+    theta = openmc.stats.Discrete(x=[0.0], p=[1.0])
+    phi = openmc.stats.Discrete(x=[0.0], p=[1.0])
+
+    space = openmc.stats.SphericalIndependent(r, theta, phi)
+    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()