diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 260ee2ad78..5eb56a3795 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -366,6 +366,7 @@ public: MOABMesh() = default; MOABMesh(pugi::xml_node); MOABMesh(const std::string& filename); + MOABMesh(std::shared_ptr 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 mbi_; //!< MOAB instance + std::shared_ptr mbi_; //!< MOAB instance unique_ptr kdtree_; //!< MOAB KDTree instance vector baryc_data_; //!< Barycentric data for tetrahedra vector tag_names_; //!< Names of score tags added to the mesh diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index aaaff0e700..25208f6fb8 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -53,6 +53,9 @@ public: void set_filters(gsl::span 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_;} diff --git a/src/mesh.cpp b/src/mesh.cpp index 09ae82c707..4f85c58c80 100644 --- a/src/mesh.cpp +++ b/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 external_mbi) { + mbi_ = external_mbi; + filename_ = "unknown (external file)"; + this->initialize(); +} + void MOABMesh::initialize() { - // create MOAB instance - mbi_ = make_unique(); - // 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(); + + // 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) { diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index b159f43b03..153440c506 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -366,9 +366,18 @@ Tally::set_filters(gsl::span 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) { diff --git a/tests/regression_tests/external_moab/__init__.py b/tests/regression_tests/external_moab/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/external_moab/inputs_true.dat b/tests/regression_tests/external_moab/inputs_true.dat new file mode 100644 index 0000000000..af259eee47 --- /dev/null +++ b/tests/regression_tests/external_moab/inputs_true.dat @@ -0,0 +1,73 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 100 + 10 + + + 0.0 0.0 0.0 + + + + 15000000.0 1.0 + + + + + + + test_mesh_tets.h5m + + + 1 + + + 1 + flux + tracklength + + diff --git a/tests/regression_tests/external_moab/main.cpp b/tests/regression_tests/external_moab/main.cpp new file mode 100644 index 0000000000..6e3709b6a1 --- /dev/null +++ b/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 + +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 moabPtrLocal = + std::make_shared(); + + // 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(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(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(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 filters(1, filter_ptr); + model::tallies.back()->set_filters(filters); + + // Set tally estimator + model::tallies.back()->estimator_ = TallyEstimator::TRACKLENGTH; + + // Set tally score + std::vector 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; +} diff --git a/tests/regression_tests/external_moab/test.py b/tests/regression_tests/external_moab/test.py new file mode 100644 index 0000000000..caff6b5dbd --- /dev/null +++ b/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() diff --git a/tests/regression_tests/external_moab/test_mesh_tets.h5m b/tests/regression_tests/external_moab/test_mesh_tets.h5m new file mode 100644 index 0000000000..dc573142f6 Binary files /dev/null and b/tests/regression_tests/external_moab/test_mesh_tets.h5m differ