Moved the dummy communicator to being within openmc vice openmc.deplete, added test of model.calculate_volumes and model.deplete, fixed issues identified from those tests, and added a DLL quieting method (still needs to be propagated to integrator.integrate

This commit is contained in:
agnelson 2021-09-30 11:18:35 -05:00
parent dc80b799ac
commit ad483b3fea
14 changed files with 407 additions and 284 deletions

View file

@ -30,8 +30,10 @@ from openmc.plotter import *
from openmc.search import *
from openmc.polynomial import *
from . import examples
from openmc.mpi import DummyCommunicator, MPI, comm
# Import a few names from the model module
from openmc.model import rectangular_prism, hexagonal_prism, Model
__version__ = '0.13.0-dev'

View file

@ -4,17 +4,6 @@ openmc.deplete
A depletion front-end tool.
"""
import sys
from unittest.mock import Mock
from .dummy_comm import DummyCommunicator
try:
from mpi4py import MPI
comm = MPI.COMM_WORLD
except ImportError:
MPI = Mock()
comm = DummyCommunicator()
from .nuclide import *
from .chain import *

View file

@ -22,7 +22,7 @@ from uncertainties import ufloat
from openmc.data import DataLibrary
from openmc.lib import MaterialFilter, Tally
from openmc.checkvalue import check_type, check_greater_than
from . import comm
from openmc import comm
from .results import Results
from .chain import Chain
from .results_list import ResultsList

View file

@ -10,7 +10,7 @@ import sys
from numpy import dot, zeros, newaxis, asarray
from . import comm
from openmc import comm
from openmc.checkvalue import check_type, check_greater_than
from openmc.data import JOULE_PER_EV, REACTION_MT
from openmc.lib import (

View file

@ -20,7 +20,7 @@ from uncertainties import ufloat
import openmc
from openmc.checkvalue import check_value
import openmc.lib
from . import comm
from openmc import comm
from .abc import TransportOperator, OperatorResult
from .atom_number import AtomNumber
from .reaction_rates import ReactionRates
@ -176,6 +176,9 @@ class Operator(TransportOperator):
results are to be used.
diff_burnable_mats : bool
Whether to differentiate burnable materials with multiple instances
cleanup_when_done : bool
Whether to finalize and clear the shared library memory when the
depletion operation is complete. Defaults to clearing the library.
"""
_fission_helpers = {
"average": AveragedFissionYieldHelper,
@ -202,7 +205,7 @@ class Operator(TransportOperator):
self.settings = settings
self.geometry = geometry
self.diff_burnable_mats = diff_burnable_mats
self._cleanup_when_done = True
self.cleanup_when_done = True
# Reduce the chain before we create more materials
if reduce_chain:
@ -318,6 +321,10 @@ class Operator(TransportOperator):
# Reset results in OpenMC
openmc.lib.reset()
# Update the number densities regardless of the source rate
self.number.set_density(vec)
self._update_materials()
# If the source rate is zero, return zero reaction rates without running
# a transport solve
if source_rate == 0.0:
@ -328,11 +335,7 @@ class Operator(TransportOperator):
# Prevent OpenMC from complaining about re-creating tallies
openmc.reset_auto_ids()
# Update status
self.number.set_density(vec)
# Update material compositions and tally nuclides
self._update_materials()
# Update tally nuclides data in preparation for transport solve
nuclides = self._get_tally_nuclides()
self._rate_helper.nuclides = nuclides
self._normalization_helper.nuclides = nuclides
@ -556,7 +559,7 @@ class Operator(TransportOperator):
def finalize(self):
"""Finalize a depletion simulation and release resources."""
if self._cleanup_when_done:
if self.cleanup_when_done:
openmc.lib.finalize()
def _update_materials(self):

View file

@ -9,7 +9,7 @@ import copy
import h5py
import numpy as np
from . import comm, MPI
from openmc import comm, MPI
from .reaction_rates import ReactionRates
VERSION_RESULTS = (1, 1)

View file

@ -2,6 +2,7 @@ from contextlib import contextmanager
from ctypes import (c_bool, c_int, c_int32, c_int64, c_double, c_char_p,
c_char, POINTER, Structure, c_void_p, create_string_buffer)
import sys
import os
import numpy as np
from numpy.ctypeslib import as_array
@ -110,9 +111,14 @@ def global_bounding_box():
return llc, urc
def calculate_volumes():
def calculate_volumes(output=True):
"""Run stochastic volume calculation"""
_dll.openmc_calculate_volumes()
if output:
_dll.openmc_calculate_volumes()
else:
with quiet_dll():
_dll.openmc_calculate_volumes()
def current_batch():
@ -127,13 +133,17 @@ def current_batch():
return c_int.in_dll(_dll, 'current_batch').value
def export_properties(filename=None):
def export_properties(filename=None, output=True):
"""Export physical properties.
.. versionchanged:: 0.13.0
Parameters
----------
filename : str or None
Filename to export properties to (defaults to "properties.h5")
output: bool, optional
Whether or not to show output. Defaults to showing output
See Also
--------
@ -142,7 +152,11 @@ def export_properties(filename=None):
"""
if filename is not None:
filename = c_char_p(filename.encode())
_dll.openmc_properties_export(filename)
if output:
_dll.openmc_properties_export(filename)
else:
with quiet_dll():
_dll.openmc_properties_export(filename)
def finalize():
@ -220,15 +234,19 @@ def import_properties(filename):
_dll.openmc_properties_import(filename.encode())
def init(args=None, intracomm=None):
def init(args=None, intracomm=None, output=True):
"""Initialize OpenMC
.. versionchanged:: 0.13.0
Parameters
----------
args : list of str
args : list of str, optional
Command-line arguments
intracomm : mpi4py.MPI.Intracomm or None
intracomm : mpi4py.MPI.Intracomm or None, optional
MPI intracommunicator
output: bool, optional
Whether or not to show output. Defaults to showing output
"""
if args is not None:
@ -254,7 +272,11 @@ def init(args=None, intracomm=None):
address = MPI._addressof(intracomm)
intracomm = c_void_p(address)
_dll.openmc_init(argc, argv, intracomm)
if output:
_dll.openmc_init(argc, argv, intracomm)
else:
with quiet_dll():
_dll.openmc_init(argc, argv, intracomm)
openmc.lib.is_initialized = True
@ -345,9 +367,22 @@ def next_batch():
return status.value
def plot_geometry():
"""Plot geometry"""
_dll.openmc_plot_geometry()
def plot_geometry(output=True):
"""Plot geometry
.. versionchanged:: 0.13.0
Parameters
----------
output: bool, optional
Whether or not to show output. Defaults to showing output
"""
if output:
_dll.openmc_plot_geometry()
else:
with quiet_dll():
_dll.openmc_plot_geometry()
def reset():
@ -360,9 +395,22 @@ def reset_timers():
_dll.openmc_reset_timers()
def run():
"""Run simulation"""
_dll.openmc_run()
def run(output=True):
"""Run simulation
.. versionchanged:: 0.13.0
Parameters
----------
output: bool, optional
Whether or not to show output. Defaults to showing output
"""
if output:
_dll.openmc_run()
else:
with quiet_dll():
_dll.openmc_run()
def simulation_init():
@ -478,3 +526,30 @@ class _FortranObjectWithID(_FortranObject):
# assigned. If the array index of the object is out of bounds, an
# OutOfBoundsError will be raised here by virtue of referencing self.id
self.id
@contextmanager
def quiet_dll():
"""This context manager allows us to suppress standard output from DLLs"""
sys.stdout.flush()
# Save the initial file descriptor states
initial_stdout = sys.stdout
initial_stdout_fno = os.dup(sys.stdout.fileno())
# Get a garbage descriptor so we can throw away output
devnull = os.open(os.devnull, os.O_WRONLY)
# Get the current stdout stream and make a duplicate of it
new_stdout = os.dup(1)
# Copy the garbage output to the stdout stream
os.dup2(devnull, 1)
os.close(devnull)
# Now point stdout to the re-defined stdout
sys.stdout = os.fdopen(new_stdout, 'w')
try:
yield
finally:
# Now we just clean up after ourselves and reset the streams
sys.stdout = initial_stdout
sys.stdout.flush()
os.dup2(initial_stdout_fno, 1)

View file

@ -191,17 +191,11 @@ class Model:
@intracomm.setter
def intracomm(self, intracomm):
try:
from mpi4py import MPI
mpi_avail = True
except ImportError:
mpi_avail = False
if intracomm is None or not mpi_avail:
# TODO: move dummy_comm from openmc.deplete to openmc
from openmc.deplete.dummy_comm import DummyCommunicator
self._intracomm = DummyCommunicator()
if intracomm is None:
self._intracomm = openmc.comm
else:
check_type('intracomm', intracomm, MPI.Comm)
check_type('intracomm', intracomm,
(openmc.MPI.Comm, openmc.DummyCommunicator))
self._intracomm = intracomm
@classmethod
@ -263,25 +257,12 @@ class Model:
# where it exists (here in init), but in the future the functionality
# should be exposed so that it can be accessed via model.run(...)
# TODO: the output flag is not yet implemented for the openmc.lib.init
# command. This will be the subject of future work.
args = []
if isinstance(threads, Integral) and threads > 0:
args += ['-s', str(threads)]
if geometry_debug:
args.append('-g')
if event_based:
args.append('-e')
if isinstance(restart_file, str):
args += ['-r', restart_file]
if tracks:
args.append('-t')
args = openmc.process_CLI_arguments(
volume=False, geometry_debug=geometry_debug,
restart_file=restart_file, threads=threads, tracks=tracks,
event_based=event_based)
# Args adds the openmc_exec command in the first entry; remove it
args = args[1:]
self.finalize_lib()
@ -289,7 +270,13 @@ class Model:
self.export_to_xml()
self.intracomm.barrier()
openmc.lib.init(args=args, intracomm=self.intracomm)
if isinstance(self.intracomm, openmc.DummyCommunicator):
# openmc.lib.init does not accept DummyCommunicator, and importing
# the DummyCommunicator class is overkill. Filter it here
intracomm = None
else:
intracomm = self.intracomm
openmc.lib.init(args=args, intracomm=intracomm, output=output)
def finalize_lib(self):
"""Finalize simulation and free memory allocated for the C API
@ -301,7 +288,7 @@ class Model:
openmc.lib.finalize()
def deplete(self, timesteps, method='cecm', final_step=True,
operator_kwargs=None, integrator_kwargs=None, directory='.'):
operator_kwargs=None, directory='.', **integrator_kwargs):
"""Deplete model using specified timesteps/power
.. versionchanged:: 0.13.0
@ -320,25 +307,34 @@ class Model:
operator_kwargs : dict
Keyword arguments passed to the depletion Operator initializer
(e.g., :func:`openmc.deplete.Operator`)
integrator_kwargs : dict
Keyword arguments passed to the depletion Operator initializer
(e.g., :func:`openmc.deplete.integrator.cecm`)
directory : str, optional
Directory to write XML files to. If it doesn't exist already, it
will be created. Defaults to the current working directory
**kwargs
Keyword arguments passed to integration function (e.g.,
:func:`openmc.deplete.integrator.cecm`)
integrator_kwargs : dict
Remaining keyword arguments passed to the depletion Integrator
initializer (e.g., :func:`openmc.deplete.integrator.cecm`).
"""
if operator_kwargs is None:
op_kwargs = {}
elif isinstance(operator_kwargs, dict):
op_kwargs = operator_kwargs
else:
msg = "operator_kwargs must be a dict or None"
raise ValueError(msg)
# Import openmc.deplete here so the Model can be used even if the
# shared library is unavailable.
import openmc.deplete as dep
# Store whether or not the library was initialized when we started
started_initialized = self.is_initialized
with _change_directory(Path(directory)):
depletion_operator = \
dep.Operator(self.geometry, self.settings, **operator_kwargs)
dep.Operator(self.geometry, self.settings, **op_kwargs)
# Tell depletion_operator.finalize NOT to clear C API memory when
# it is done
depletion_operator.cleanup_when_done = False
@ -351,13 +347,9 @@ class Model:
**integrator_kwargs)
# Now perform the depletion
# TODO: add output parameter to integrate
integrator.integrate(final_step)
# If we did not perform a transport calculation on the final step,
# then make the code update the C API material inventory
if not final_step:
depletion_operator._update_materials()
# Now make the python Materials match the C API material data
for mat_id, mat in self._materials_by_id.items():
if mat.depletable:
@ -370,6 +362,11 @@ class Model:
mat.add_nuclide(nuc, density)
mat.set_density('atom/b-cm', sum(densities))
# If we didnt start intialized, we should cleanup after ourselves
if not started_initialized:
depletion_operator.cleanup_when_done = True
depletion_operator.finalize()
def export_to_xml(self, directory='.'):
"""Export model to XML files.
@ -531,25 +528,16 @@ class Model:
if isinstance(particles, Integral) and particles > 0:
openmc.lib.settings.particles = particles
# If we dont want output, make the verbosity quiet
if not output:
init_verbosity = openmc.lib.settings.verbosity
if not output:
openmc.lib.settings.verbosity = 1
# Event-based can be set at init-time or on a case-basis.
# Handle the argument here.
# TODO This will be dealt with in a future change to the C API
# Then run using the C API
openmc.lib.run()
openmc.lib.run(output)
# Reset changes for the openmc.run kwargs handling
if particles is not None:
openmc.lib.settings.particles = init_particles
if not output:
# Then re-set the initial verbosity
openmc.lib.settings.verbosity = init_verbosity
else:
# Then run via the command line
@ -603,7 +591,7 @@ class Model:
to the model. Defaults to applying the volumes.
"""
if len(self.settings.volume_calculation) == 0:
if len(self.settings.volume_calculations) == 0:
# Then there is no volume calculation specified
raise ValueError("The Settings.volume_calculation attribute must"
" be specified before executing this method!")
@ -618,18 +606,9 @@ class Model:
"with the call to mpi4py"
raise ValueError(msg)
# Apply the output settings
if not output:
init_verbosity = openmc.lib.settings.verbosity
if not output:
openmc.lib.settings.verbosity = 1
# Compute the volumes
openmc.lib.calculate_volumes()
openmc.lib.calculate_volumes(output)
# Reset the output verbosity
if not output:
openmc.lib.settings.verbosity = init_verbosity
else:
self.export_to_xml()
openmc.calculate_volumes(threads=threads, output=output,
@ -638,23 +617,19 @@ class Model:
# Now we apply the volumes
if apply_volumes:
# Load the results
f_names = \
[f"volume_{i + 1}.h5"
for i in range(len(self.settings.volume_calculations))]
vol_calcs = [openmc.VolumeCalculation.load_results(f_name)
for f_name in f_names]
# And now we can add them to the model
for vol_calc in vol_calcs:
# Load the results and add them to the model
for i, vol_calc in enumerate(self.settings.volume_calculations):
f_name = f"volume_{i + 1}.h5"
vol_calc.load_results(f_name)
# First add them to the Python side
self.geometry.add_volume_information(vol_calc)
# And now repeat for the C API
if vol_calc.domain == 'material':
if self.is_initialized and vol_calc.domain_type == 'material':
# Then we can do this in the C API
for domain_id in vol_calc.domains:
self.update_material_volumes(
[domain_id], vol_calc.volumes[domain_id])
for domain_id in vol_calc.ids:
openmc.lib.materials[domain_id].volume = \
vol_calc.volumes[domain_id].n
def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc',
convert=True, convert_exec='convert'):
@ -689,25 +664,15 @@ class Model:
"before executing this method!")
with _change_directory(Path(cwd)):
# TODO: openmis_initialized doesnt read plots.xml unless it is in plot mode
# so the following will not work. Commented out for now and
# replacing with non-C API code
# TODO: openmc.is_initialized doesnt read plots.xml unless it is
# in plot mode so the following will not work. Commented out for
# now and replacing with non-C API code
self.export_to_xml()
openmc.plot_geometry(output=output, openmc_exec=openmc_exec)
# if self.is_initialized:
# # Apply the output settings
# if not output:
# init_verbosity = openmc.lib.settings.verbosity
# if not output:
# openmc.lib.settings.verbosity = 1
# # Compute the volumes
# openmc.lib.plot_geometry()
# # Reset the output verbosity
# if not output:
# openmc.lib.settings.verbosity = init_verbosity
# openmc.lib.plot_geometry(output)
# else:
# self.export_to_xml()
# openmc.plot_geometry(output=output, openmc_exec=openmc_exec)
@ -721,8 +686,8 @@ class Model:
png_file = ppm_file.replace('.ppm', '.png')
subprocess.check_call([convert_exec, ppm_file, png_file])
def _change_py_C_attribs(self, names_or_ids, value, obj_type, attrib_name,
density_units='atom/b-cm'):
def _change_py_lib_attribs(self, names_or_ids, value, obj_type,
attrib_name, density_units='atom/b-cm'):
# Method to do the same work whether it is a cell or material and
# a temperature or volume
check_type('names_or_ids', names_or_ids, Iterable, (Integral, str))
@ -819,7 +784,7 @@ class Model:
"""
self._change_py_C_attribs(names_or_ids, vector, 'cell', 'rotation')
self._change_py_lib_attribs(names_or_ids, vector, 'cell', 'rotation')
def translate_cells(self, names_or_ids, vector):
"""Translate the identified cell(s) by the specified translation vector.
@ -841,7 +806,7 @@ class Model:
"""
self._change_py_C_attribs(names_or_ids, vector, 'cell', 'translation')
self._change_py_lib_attribs(names_or_ids, vector, 'cell', 'translation')
def update_densities(self, names_or_ids, density, density_units='atom/b-cm'):
"""Update the density of a given set of materials to a new value
@ -863,7 +828,7 @@ class Model:
"""
self._change_py_C_attribs(names_or_ids, density, 'material', 'density',
self._change_py_lib_attribs(names_or_ids, density, 'material', 'density',
density_units)
def update_cell_temperatures(self, names_or_ids, temperature):
@ -884,7 +849,7 @@ class Model:
"""
self._change_py_C_attribs(names_or_ids, temperature, 'cell',
self._change_py_lib_attribs(names_or_ids, temperature, 'cell',
'temperature')
def update_material_temperatures(self, names_or_ids, temperature):
@ -905,7 +870,7 @@ class Model:
"""
self._change_py_C_attribs(names_or_ids, temperature, 'material',
self._change_py_lib_attribs(names_or_ids, temperature, 'material',
'temperature')
def update_material_volumes(self, names_or_ids, volume):
@ -926,4 +891,4 @@ class Model:
"""
self._change_py_C_attribs(names_or_ids, volume, 'material', 'volume')
self._change_py_lib_attribs(names_or_ids, volume, 'material', 'volume')

View file

@ -1,4 +1,5 @@
import sys
from unittest.mock import Mock
class DummyCommunicator:
@ -31,3 +32,10 @@ class DummyCommunicator:
def Abort(self, exit_code_or_msg):
sys.exit(exit_code_or_msg)
try:
from mpi4py import MPI
comm = MPI.COMM_WORLD
except ImportError:
MPI = Mock()
comm = DummyCommunicator()

View file

@ -1,138 +1,9 @@
from math import pi
import openmc
import pytest
from tests.regression_tests import config
@pytest.fixture(scope='module')
def pin_model_attributes():
uo2 = openmc.Material(name='UO2')
uo2.set_density('g/cm3', 10.29769)
uo2.add_element('U', 1., enrichment=2.4)
uo2.add_element('O', 2.)
zirc = openmc.Material(name='Zirc')
zirc.set_density('g/cm3', 6.55)
zirc.add_element('Zr', 1.)
borated_water = openmc.Material(name='Borated water')
borated_water.set_density('g/cm3', 0.740582)
borated_water.add_element('B', 4.0e-5)
borated_water.add_element('H', 5.0e-2)
borated_water.add_element('O', 2.4e-2)
borated_water.add_s_alpha_beta('c_H_in_H2O')
mats = openmc.Materials([uo2, zirc, borated_water])
pitch = 1.25984
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
box = openmc.model.rectangular_prism(pitch, pitch,
boundary_type='reflective')
# Define cells
fuel_inf_cell = openmc.Cell(name='inf fuel', fill=uo2)
fuel_inf_univ = openmc.Universe(cells=[fuel_inf_cell])
fuel = openmc.Cell(name='fuel', fill=fuel_inf_univ, region=-fuel_or)
clad = openmc.Cell(fill=zirc, region=+fuel_or & -clad_or)
water = openmc.Cell(fill=borated_water, region=+clad_or & box)
# Define overall geometry
geom = openmc.Geometry([fuel, clad, water])
uo2.volume = pi * fuel_or.r**2
settings = openmc.Settings()
settings.batches = 100
settings.inactive = 10
settings.particles = 1000
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.source.Source(space=uniform_dist)
entropy_mesh = openmc.RegularMesh()
entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]
entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50]
entropy_mesh.dimension = [10, 10, 1]
settings.entropy_mesh = entropy_mesh
tals = openmc.Tallies()
tal = openmc.Tally(name='test')
tal.scores = ['flux']
tals.append(tal)
plot1 = openmc.Plot()
plot1.origin = (0., 0., 0.)
plot1.width = (pitch, pitch)
plot1.pixels = (300, 300)
plot1.color_by = 'material'
plot1.filename = 'test'
plot2 = openmc.Plot()
plot2.origin = (0., 0., 0.)
plot2.width = (pitch, pitch)
plot2.pixels = (300, 300)
plot2.color_by = 'cell'
plots = openmc.Plots((plot1, plot2))
chain = './chain_simple.xml'
fission_q = {'U235': 200e6}
chain_file_xml = """<?xml version="1.0"?>
<depletion_chain>
<nuclide name="I135" decay_modes="1" reactions="1" half_life="2.36520E+04">
<decay type="beta" target="Xe135" branching_ratio="1.0" />
<reaction type="(n,gamma)" Q="0.0" target="Xe136" /> <!-- Not precisely true, but whatever -->
</nuclide>
<nuclide name="Xe135" decay_modes="1" reactions="1" half_life="3.29040E+04">
<decay type=" beta" target="Cs135" branching_ratio="1.0" />
<reaction type="(n,gamma)" Q="0.0" target="Xe136" />
</nuclide>
<nuclide name="Xe136" decay_modes="0" reactions="0" />
<nuclide name="Cs135" decay_modes="0" reactions="0" />
<nuclide name="Gd157" decay_modes="0" reactions="1" >
<reaction type="(n,gamma)" Q="0.0" target="Nothing" />
</nuclide>
<nuclide name="Gd156" decay_modes="0" reactions="1">
<reaction type="(n,gamma)" Q="0.0" target="Gd157" />
</nuclide>
<nuclide name="U234" decay_modes="0" reactions="1">
<reaction type="fission" Q="191840000."/>
<neutron_fission_yields>
<energies>2.53000e-02</energies>
<fission_yields energy="2.53000e-02">
<products>Gd157 Gd156 I135 Xe135 Xe136 Cs135</products>
<data>1.093250e-04 2.087260e-04 2.780820e-02 6.759540e-03 2.392300e-02 4.356330e-05</data>
</fission_yields>
</neutron_fission_yields>
</nuclide>
<nuclide name="U235" decay_modes="0" reactions="1">
<reaction type="fission" Q="193410000."/>
<neutron_fission_yields>
<energies>2.53000e-02</energies>
<fission_yields energy="2.53000e-02">
<products>Gd157 Gd156 I135 Xe135 Xe136 Cs135</products>
<data>6.142710e-5 1.483250e-04 0.0292737 0.002566345 0.0219242 4.9097e-6</data>
</fission_yields>
</neutron_fission_yields>
</nuclide>
<nuclide name="U238" decay_modes="0" reactions="1">
<reaction type="fission" Q="197790000."/>
<neutron_fission_yields>
<energies>2.53000e-02</energies>
<fission_yields energy="2.53000e-02">
<products>Gd157 Gd156 I135 Xe135 Xe136 Cs135</products>
<data>4.141120e-04 7.605360e-04 0.0135457 0.00026864 0.0024432 3.7100E-07</data>
</fission_yields>
</neutron_fission_yields>
</nuclide>
</depletion_chain>
"""
operator_kwargs = {'chain_file': chain, 'fission_q': fission_q}
return (mats, geom, settings, tals, plots, operator_kwargs, chain_file_xml)
@pytest.fixture(scope='module')
def mpi_intracomm():
if config['mpi']:

View file

@ -7,7 +7,8 @@ import os
from pathlib import Path
import numpy as np
from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram, pool
from openmc import comm
from openmc.deplete import Chain, reaction_rates, nuclide, cram, pool
import pytest
from tests import cdtemp

View file

@ -14,11 +14,11 @@ import numpy as np
from uncertainties import ufloat
import pytest
from openmc import comm
from openmc.deplete import (
ReactionRates, Results, ResultsList, comm, OperatorResult,
PredictorIntegrator, CECMIntegrator, CF4Integrator, CELIIntegrator,
EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator,
cram)
ReactionRates, Results, ResultsList, OperatorResult, PredictorIntegrator,
CECMIntegrator, CF4Integrator, CELIIntegrator, EPCRK4Integrator,
LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator, cram)
from tests import dummy_operator

View file

@ -713,7 +713,6 @@ def test_trigger_set_n_batches(uo2_trigger_model, mpi_intracomm):
def test_cell_translation(pincell_model_w_univ, mpi_intracomm):
openmc.lib.finalize()
openmc.lib.init(intracomm=mpi_intracomm)
openmc.lib.simulation_init()
# Cell 1 is filled with a material so it has a translation, but we can't
# set it.
cell = openmc.lib.cells[1]
@ -727,9 +726,12 @@ def test_cell_translation(pincell_model_w_univ, mpi_intracomm):
# This time we *can* set it
cell.translation = (1., 0., -1.)
assert cell.translation == pytest.approx([1., 0., -1.])
openmc.lib.finalize()
def test_cell_rotation(pincell_model_w_univ):
def test_cell_rotation(pincell_model_w_univ, mpi_intracomm):
openmc.lib.finalize()
openmc.lib.init(intracomm=mpi_intracomm)
# Cell 1 is filled with a material so we cannot rotate it, but we can get
# its rotation matrix (which will be the identity matrix)
cell = openmc.lib.cells[1]
@ -742,3 +744,4 @@ def test_cell_rotation(pincell_model_w_univ):
assert cell.rotation == pytest.approx([0., 0., 0.])
cell.rotation = (180., 0., 0.)
assert cell.rotation == pytest.approx([180., 0., 0.])
openmc.lib.finalize()

View file

@ -1,10 +1,110 @@
from math import pi
import numpy as np
import pytest
from pathlib import Path
from shutil import which
import openmc
import openmc.lib
from openmc.deplete.dummy_comm import DummyCommunicator
from openmc.mpi import DummyCommunicator
@pytest.fixture(scope='function')
def pin_model_attributes():
uo2 = openmc.Material(material_id=1, name='UO2')
uo2.set_density('g/cm3', 10.29769)
uo2.add_element('U', 1., enrichment=2.4)
uo2.add_element('O', 2.)
uo2.depletable = True
zirc = openmc.Material(material_id=2, name='Zirc')
zirc.set_density('g/cm3', 6.55)
zirc.add_element('Zr', 1.)
zirc.depletable = False
borated_water = openmc.Material(material_id=3, name='Borated water')
borated_water.set_density('g/cm3', 0.740582)
borated_water.add_element('B', 4.0e-5)
borated_water.add_element('H', 5.0e-2)
borated_water.add_element('O', 2.4e-2)
borated_water.add_s_alpha_beta('c_H_in_H2O')
borated_water.depletable = False
mats = openmc.Materials([uo2, zirc, borated_water])
pitch = 1.25984
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
box = openmc.model.rectangular_prism(pitch, pitch,
boundary_type='reflective')
# Define cells
fuel_inf_cell = openmc.Cell(cell_id=1, name='inf fuel', fill=uo2)
fuel_inf_univ = openmc.Universe(universe_id=1, cells=[fuel_inf_cell])
fuel = openmc.Cell(cell_id=2, name='fuel',
fill=fuel_inf_univ, region=-fuel_or)
clad = openmc.Cell(cell_id=3, fill=zirc, region=+fuel_or & -clad_or)
water = openmc.Cell(cell_id=4, fill=borated_water, region=+clad_or & box)
# Define overall geometry
geom = openmc.Geometry([fuel, clad, water])
uo2.volume = pi * fuel_or.r**2
settings = openmc.Settings()
settings.batches = 100
settings.inactive = 10
settings.particles = 1000
# Create a uniform spatial source distribution over fissionable zones
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
uniform_dist = openmc.stats.Box(
bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.source.Source(space=uniform_dist)
entropy_mesh = openmc.RegularMesh()
entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]
entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50]
entropy_mesh.dimension = [10, 10, 1]
settings.entropy_mesh = entropy_mesh
tals = openmc.Tallies()
tal = openmc.Tally(tally_id=1, name='test')
tal.filters = [openmc.MaterialFilter(bins=[uo2])]
tal.scores = ['flux', 'fission']
tals.append(tal)
plot1 = openmc.Plot(plot_id=1)
plot1.origin = (0., 0., 0.)
plot1.width = (pitch, pitch)
plot1.pixels = (300, 300)
plot1.color_by = 'material'
plot1.filename = 'test'
plot2 = openmc.Plot(plot_id=2)
plot2.origin = (0., 0., 0.)
plot2.width = (pitch, pitch)
plot2.pixels = (300, 300)
plot2.color_by = 'cell'
plots = openmc.Plots((plot1, plot2))
chain = './test_chain.xml'
chain_file_xml = """<?xml version="1.0"?>
<depletion_chain>
<nuclide name="Xe136" decay_modes="0" reactions="0" />
<nuclide name="U235" decay_modes="0" reactions="1">
<reaction type="fission" Q="200000000."/>
<neutron_fission_yields>
<energies>2.53000e-02</energies>
<fission_yields energy="2.53000e-02">
<products>Xe136</products>
<data>1.0</data>
</fission_yields>
</neutron_fission_yields>
</nuclide>
</depletion_chain>
"""
operator_kwargs = {'chain_file': chain}
return (mats, geom, settings, tals, plots, operator_kwargs, chain_file_xml)
def test_init(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
@ -130,7 +230,7 @@ def test_init_finalize_lib(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
# We are going to init and then make sure data is loaded
mats, geom, settings, tals, plots, _, _ = pin_model_attributes
test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
test_model.init_lib()
test_model.init_lib(output=False)
# First check that the API is advertised as initialized
assert openmc.lib.is_initialized is True
@ -157,13 +257,13 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm):
# Create PWR pin cell model and write XML files
openmc.reset_auto_ids()
model = openmc.examples.pwr_pin_cell()
model.init_lib()
model.init_lib(output=False)
# Change fuel temperature and density and export properties
cell = openmc.lib.cells[1]
cell.set_temperature(600.0)
cell.fill.set_density(5.0, 'g/cm3')
openmc.lib.export_properties()
openmc.lib.export_properties(output=False)
# Import properties to existing model
model.import_properties("properties.h5")
@ -203,17 +303,20 @@ def test_run(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
sp_path = test_model.run(output=False)
with openmc.StatePoint(sp_path) as sp:
cli_keff = sp.k_combined
cli_flux = sp.get_tally(id=1).get_values()[0, 0, 0]
cli_flux = sp.get_tally(id=1).get_values(scores=['flux'])[0, 0, 0]
cli_fiss = sp.get_tally(id=1).get_values(scores=['fission'])[0, 0, 0]
test_model.init_lib()
test_model.init_lib(output=False)
sp_path = test_model.run(output=False)
with openmc.StatePoint(sp_path) as sp:
C_keff = sp.k_combined
C_flux = sp.get_tally(id=1).get_values()[0, 0, 0]
lib_keff = sp.k_combined
lib_flux = sp.get_tally(id=1).get_values(scores=['flux'])[0, 0, 0]
lib_fiss = sp.get_tally(id=1).get_values(scores=['fission'])[0, 0, 0]
# and lets compare results
assert abs(C_keff - cli_keff) < 1e-13
assert abs(C_flux - cli_flux) < 1e-13
assert abs(lib_keff - cli_keff) < 1e-13
assert abs(lib_flux - cli_flux) < 1e-13
assert abs(lib_fiss - cli_fiss) < 1e-13
# Now we should make sure that the flags for items which should be handled
# by init are properly set
@ -248,8 +351,8 @@ def test_plots(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
# We will run the test twice, the first time without C API, the second with
for i in range(2):
if i == 1:
test_model.init_lib()
test_model.plot_geometry(output=True, convert=convert)
test_model.init_lib(output=False)
test_model.plot_geometry(output=False, convert=convert)
# Now look for the files, expect to find test.ppm, plot_2.ppm, and if
# convert is True, test.png, plot_2.png
@ -262,11 +365,11 @@ def test_plots(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
test_model.finalize_lib()
def test_py_C_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
def test_py_lib_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
mats, geom, settings, tals, plots, _, _ = pin_model_attributes
test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
test_model.init_lib()
test_model.init_lib(output=False)
# Now we can call rotate_cells, translate_cells, update_densities,
# update_cell_temperatures, and update_material_temperatures and make sure
@ -358,3 +461,106 @@ def test_py_C_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
assert abs(test_model.materials[0].volume - 2.) < 1e-13
test_model.finalize_lib()
# def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
# mats, geom, settings, tals, plots, op_kwargs, chain_file_xml = \
# pin_model_attributes
# with open('test_chain.xml', 'w') as f:
# f.write(chain_file_xml)
# test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
# initial_mat = mats[0].clone()
# initial_u = initial_mat.get_nuclide_atom_densities()['U235'][1]
# # Note that the chain file includes only U-235 fission to a stable Xe136 w/
# # a yield of 100%. Thus all the U235 we lose becomes Xe136
# # In this test we first run without pre-initializing the shared library
# # data and then compare. Then we repeat with the C API already initialized
# # and make sure we get the same answer
# test_model.deplete([1e6], 'predictor', final_step=False,
# operator_kwargs=op_kwargs,
# power=1.)
# # Get the new Xe136 and U235 atom densities
# after_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1]
# after_u = mats[0].get_nuclide_atom_densities()['U235'][1]
# assert abs((after_xe + after_u) - initial_u) < 1e-15
# assert test_model.is_initialized is False
# # Reset the initial material densities
# mats[0].nuclides.clear()
# densities = initial_mat.get_nuclide_atom_densities()
# tot_density = 0.
# for nuc, density in densities.values():
# mats[0].add_nuclide(nuc, density)
# tot_density += density
# mats[0].set_density('atom/b-cm', tot_density)
# # Now we can re-run with the pre-initialized API
# test_model.init_lib(output=False)
# test_model.deplete([1e6], 'predictor', final_step=False,
# operator_kwargs=op_kwargs,
# power=1.)
# # Get the new Xe136 and U235 atom densities
# after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1]
# after_lib_u = mats[0].get_nuclide_atom_densities()['U235'][1]
# assert abs((after_lib_xe + after_lib_u) - initial_u) < 1e-15
# assert test_model.is_initialized is True
# # And end by comparing to the previous case
# assert abs(after_xe - after_lib_xe) < 1e-15
# assert abs(after_u - after_lib_u) < 1e-15
# test_model.finalize_lib()
def test_calc_volumes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
mats, geom, settings, tals, plots, _, _ = pin_model_attributes
test_model = openmc.Model(geom, mats, settings, tals, plots, mpi_intracomm)
# With no vol calcs, it should fail
with pytest.raises(ValueError):
test_model.calculate_volumes(output=False)
# Add a cell and mat volume calc
material_vol_calc = openmc.VolumeCalculation(
[mats[2]], samples=1000, lower_left=(-.63, -.63, -100.),
upper_right=(.63, .63, 100.))
cell_vol_calc = openmc.VolumeCalculation(
[geom.root_universe.cells[3]], samples=1000,
lower_left=(-.63, -.63, -100.), upper_right=(.63, .63, 100.))
test_model.settings.volume_calculations = \
[material_vol_calc, cell_vol_calc]
# Now lets compute the volumes and check to see if it was applied
# First lets do without using the C-API
# Make sure the volumes are unassigned first
assert mats[2].volume is None
assert geom.root_universe.cells[3].volume is None
test_model.calculate_volumes(output=False, apply_volumes=True)
# Now let's test that we have volumes assigned; we arent checking the
# value, just that the value was changed
assert mats[2].volume > 0.
assert geom.root_universe.cells[3].volume > 0.
# Now reset the values
mats[2].volume = None
geom.root_universe.cells[3].volume = None
# And do again with an initialized library
for file in ['volume_1.h5', 'volume_2.h5']:
file = Path(file)
file.unlink()
test_model.init_lib(output=False)
test_model.calculate_volumes(output=False, apply_volumes=True)
assert mats[2].volume > 0.
assert geom.root_universe.cells[3].volume > 0.
assert openmc.lib.materials[3].volume == mats[2].volume
test_model.finalize_lib()