OpenMC/tests/regression_tests/external_moab/test.py
Patrick Shriwise eda39ad9ca
Python source class refactor (#2524)
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
2023-06-20 21:27:55 -05:00

274 lines
8.7 KiB
Python

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)
target_compile_features(main PUBLIC cxx_std_14)
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']:
mpi_arg = "On"
else:
mpi_arg = "Off"
try:
print("Building driver")
# Run cmake/make to build the shared libary
subprocess.run(['cmake', os.path.pardir, f'-DOPENMC_USE_MPI={mpi_arg}'], 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)
# 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.IndependentSource(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()