crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1884 from nelsonag/model_control

Browse source 
Expanded capabilities of openmc.model.Model container
This commit is contained in:
Paul Romano 2021-10-05 12:46:51 -05:00 • committed by GitHub
commit c3bb352f8f
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
24 changed files with 1387 additions and 135 deletions
Download patch file Download diff file Expand all files Collapse all files

7
include/openmc/plot.h
View file

@ -30,9 +30,7 @@ namespace model {
extern std::unordered_map<int, int> plot_map; //!< map of plot ids to index
extern vector<Plot> plots; //!< Plot instance container
extern uint64_t
plotter_prn_seeds[N_STREAMS]; // Random number seeds used for plotter
extern int plotter_stream; // Stream index used by the plotter
extern uint64_t plotter_seed; // Stream index used by the plotter
} // namespace model
@ -274,6 +272,9 @@ void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
//! Read plot specifications from a plots.xml file
void read_plots_xml();
//! Clear memory
void free_memory_plot();
//! Create a ppm image for a plot object
//! \param[in] plot object
void create_ppm(Plot const& pl);

2
include/openmc/settings.h
View file

@ -32,7 +32,7 @@ extern "C" bool cmfd_run; //!< is a CMFD run?
extern bool
delayed_photon_scaling; //!< Scale fission photon yield to include delayed
extern "C" bool entropy_on; //!< calculate Shannon entropy?
extern bool event_based; //!< use event-based mode (instead of history-based)
extern "C" bool event_based; //!< use event-based mode (instead of history-based)
extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular?
extern bool material_cell_offsets; //!< create material cells offsets?
extern "C" bool output_summary; //!< write summary.h5?

1
openmc/__init__.py
View file

@ -34,4 +34,5 @@ from . import examples
# Import a few names from the model module
from openmc.model import rectangular_prism, hexagonal_prism, Model
__version__ = '0.13.0-dev'

11
openmc/deplete/__init__.py
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 *

2
openmc/deplete/abc.py
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.mpi import comm
from .results import Results
from .chain import Chain
from .results_list import ResultsList

2
openmc/deplete/helpers.py
View file

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

11
openmc/deplete/integrators.py
View file

@ -581,3 +581,14 @@ class SILEQIIntegrator(SIIntegrator):
proc_time += time1 + time2
return proc_time, [eos_conc, inter_conc], [res_bar]
integrator_by_name = {
'cecm': CECMIntegrator,
'predictor': PredictorIntegrator,
'cf4': CF4Integrator,
'epc_rk4': EPCRK4Integrator,
'si_celi': SICELIIntegrator,
'si_leqi': SILEQIIntegrator,
'celi': CELIIntegrator,
'leqi': LEQIIntegrator}

22
openmc/deplete/operator.py
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.mpi 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,6 +205,7 @@ class Operator(TransportOperator):
self.settings = settings
self.geometry = geometry
self.diff_burnable_mats = diff_burnable_mats
self.cleanup_when_done = True
# Reduce the chain before we create more materials
if reduce_chain:
@ -317,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:
@ -327,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
@ -536,7 +540,8 @@ class Operator(TransportOperator):
# Initialize OpenMC library
comm.barrier()
openmc.lib.init(intracomm=comm)
if not openmc.lib.is_initialized:
openmc.lib.init(intracomm=comm)
# Generate tallies in memory
materials = [openmc.lib.materials[int(i)]
@ -554,7 +559,8 @@ class Operator(TransportOperator):
def finalize(self):
"""Finalize a depletion simulation and release resources."""
openmc.lib.finalize()
if self.cleanup_when_done:
openmc.lib.finalize()
def _update_materials(self):
"""Updates material compositions in OpenMC on all processes."""

2
openmc/deplete/results.py
View file

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

0
openmc/deplete/dummy_comm.py → openmc/dummy_comm.py
View file

117
openmc/executor.py
View file

@ -5,6 +5,83 @@ import subprocess
import openmc
def _process_CLI_arguments(volume=False, geometry_debug=False, particles=None,
plot=False, restart_file=None, threads=None,
tracks=False, event_based=None,
openmc_exec='openmc', mpi_args=None):
"""Converts user-readable flags in to command-line arguments to be run with
the OpenMC executable via subprocess.
Parameters
----------
volume : bool, optional
Run in stochastic volume calculation mode. Defaults to False.
geometry_debug : bool, optional
Turn on geometry debugging during simulation. Defaults to False.
particles : int, optional
Number of particles to simulate per generation.
plot : bool, optional
Run in plotting mode. Defaults to False.
restart_file : str, optional
Path to restart file to use
threads : int, optional
Number of OpenMP threads. If OpenMC is compiled with OpenMP threading
enabled, the default is implementation-dependent but is usually equal
to the number of hardware threads available (or a value set by the
:envvar:`OMP_NUM_THREADS` environment variable).
tracks : bool, optional
Write tracks for all particles. Defaults to False.
event_based : None or bool, optional
Turns on event-based parallelism if True. If None, the value in
the Settings will be used.
openmc_exec : str, optional
Path to OpenMC executable. Defaults to 'openmc'.
mpi_args : list of str, optional
MPI execute command and any additional MPI arguments to pass,
e.g. ['mpiexec', '-n', '8'].
.. versionadded:: 0.13.0
Returns
-------
args : Iterable of str
The runtime flags converted to CLI arguments of the OpenMC executable
"""
args = [openmc_exec]
if volume:
args.append('--volume')
if isinstance(particles, Integral) and particles > 0:
args += ['-n', str(particles)]
if isinstance(threads, Integral) and threads > 0:
args += ['-s', str(threads)]
if geometry_debug:
args.append('-g')
if event_based is not None:
if event_based:
args.append('-e')
if isinstance(restart_file, str):
args += ['-r', restart_file]
if tracks:
args.append('-t')
if plot:
args.append('-p')
if mpi_args is not None:
args = mpi_args + args
return args
def _run(args, output, cwd):
# Launch a subprocess
p = subprocess.Popen(args, cwd=cwd, stdout=subprocess.PIPE,
@ -126,8 +203,8 @@ def calculate_volumes(threads=None, output=True, cwd='.',
----------
threads : int, optional
Number of OpenMP threads. If OpenMC is compiled with OpenMP threading
enabled, the default is implementation-dependent but is usually equal to
the number of hardware threads available (or a value set by the
enabled, the default is implementation-dependent but is usually equal
to the number of hardware threads available (or a value set by the
:envvar:`OMP_NUM_THREADS` environment variable).
output : bool, optional
Capture OpenMC output from standard out
@ -150,12 +227,9 @@ def calculate_volumes(threads=None, output=True, cwd='.',
openmc.VolumeCalculation
"""
args = [openmc_exec, '--volume']
if isinstance(threads, Integral) and threads > 0:
args += ['-s', str(threads)]
if mpi_args is not None:
args = mpi_args + args
args = _process_CLI_arguments(volume=True, threads=threads,
openmc_exec=openmc_exec, mpi_args=mpi_args)
_run(args, output, cwd)
@ -171,8 +245,8 @@ def run(particles=None, threads=None, geometry_debug=False,
Number of particles to simulate per generation.
threads : int, optional
Number of OpenMP threads. If OpenMC is compiled with OpenMP threading
enabled, the default is implementation-dependent but is usually equal to
the number of hardware threads available (or a value set by the
enabled, the default is implementation-dependent but is usually equal
to the number of hardware threads available (or a value set by the
:envvar:`OMP_NUM_THREADS` environment variable).
geometry_debug : bool, optional
Turn on geometry debugging during simulation. Defaults to False.
@ -201,27 +275,10 @@ def run(particles=None, threads=None, geometry_debug=False,
If the `openmc` executable returns a non-zero status
"""
args = [openmc_exec]
if isinstance(particles, Integral) and particles > 0:
args += ['-n', str(particles)]
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')
if mpi_args is not None:
args = mpi_args + args
args = _process_CLI_arguments(
volume=False, geometry_debug=geometry_debug, particles=particles,
restart_file=restart_file, threads=threads, tracks=tracks,
event_based=event_based, openmc_exec=openmc_exec, mpi_args=mpi_args)
_run(args, output, cwd)

5
openmc/lib/__init__.py
View file

@ -59,3 +59,8 @@ from .tally import *
from .settings import settings
from .math import *
from .plot import *
# Flag to denote whether or not openmc.lib.init has been called
# TODO: Establish and use a flag in the C++ code to represent the status of the
# openmc_init and openmc_finalize methods
is_initialized = False

124
openmc/lib/core.py
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
@ -109,9 +110,22 @@ def global_bounding_box():
return llc, urc
def calculate_volumes():
"""Run stochastic volume calculation"""
_dll.openmc_calculate_volumes()
def calculate_volumes(output=True):
"""Run stochastic volume calculation
.. versionchanged:: 0.13.0
The *output* argument was added.
Parameters
----------
output : bool, optional
Whether or not to show output. Defaults to showing output
"""
with quiet_dll(output):
_dll.openmc_calculate_volumes()
def current_batch():
@ -126,13 +140,18 @@ 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
The *output* argument was added.
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
--------
@ -141,12 +160,15 @@ def export_properties(filename=None):
"""
if filename is not None:
filename = c_char_p(filename.encode())
_dll.openmc_properties_export(filename)
with quiet_dll(output):
_dll.openmc_properties_export(filename)
def finalize():
"""Finalize simulation and free memory"""
_dll.openmc_finalize()
openmc.lib.is_initialized = False
def find_cell(xyz):
@ -218,15 +240,20 @@ 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
The *output* argument was added.
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:
@ -252,7 +279,9 @@ def init(args=None, intracomm=None):
address = MPI._addressof(intracomm)
intracomm = c_void_p(address)
_dll.openmc_init(argc, argv, intracomm)
with quiet_dll(output):
_dll.openmc_init(argc, argv, intracomm)
openmc.lib.is_initialized = True
def is_statepoint_batch():
@ -342,9 +371,20 @@ 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
The *output* argument was added.
Parameters
----------
output : bool, optional
Whether or not to show output. Defaults to showing output
"""
with quiet_dll(output):
_dll.openmc_plot_geometry()
def reset():
@ -357,9 +397,20 @@ def reset_timers():
_dll.openmc_reset_timers()
def run():
"""Run simulation"""
_dll.openmc_run()
def run(output=True):
"""Run simulation
.. versionchanged:: 0.13.0
The *output* argument was added.
Parameters
----------
output : bool, optional
Whether or not to show output. Defaults to showing output
"""
with quiet_dll(output):
_dll.openmc_run()
def simulation_init():
@ -475,3 +526,46 @@ 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(output=True):
"""This context manager allows us to suppress standard output from DLLs
Parameters
----------
output : bool
Denotes whether the output should be displayed (True) or not (False)
.. versionadded:: 0.13.0
"""
# This contextmanager is modified from that provided here:
# https://stackoverflow.com/a/14797594
if output:
yield
else:
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)

1
openmc/lib/material.py
View file

@ -172,7 +172,6 @@ class Material(_FortranObjectWithID):
@property
def nuclides(self):
return self._get_densities()[0]
return nuclides
@property
def densities(self):

1
openmc/lib/settings.py
View file

@ -34,6 +34,7 @@ class _Settings:
restart_run = _DLLGlobal(c_bool, 'restart_run')
run_CE = _DLLGlobal(c_bool, 'run_CE')
verbosity = _DLLGlobal(c_int, 'verbosity')
event_based = _DLLGlobal(c_bool, 'event_based')
@property
def run_mode(self):

645
openmc/model/model.py
View file

@ -1,11 +1,32 @@
from collections.abc import Iterable
import operator
import os
from pathlib import Path
from numbers import Integral
import time
import warnings
import subprocess
from contextlib import contextmanager
import h5py
import openmc
from openmc.dummy_comm import DummyCommunicator
from openmc.executor import _process_CLI_arguments
from openmc.checkvalue import check_type, check_value
from openmc.exceptions import InvalidIDError
@contextmanager
def _change_directory(working_dir):
"""A context manager for executing in a provided working directory"""
start_dir = Path.cwd()
Path.mkdir(working_dir, exist_ok=True)
os.chdir(working_dir)
try:
yield
finally:
os.chdir(start_dir)
class Model:
@ -19,6 +40,10 @@ class Model:
not set, it will attempt to create a ``materials.xml`` file based on all
materials appearing in the geometry.
.. versionchanged:: 0.13.0
The model information can now be loaded in to OpenMC directly via
openmc.lib
Parameters
----------
geometry : openmc.Geometry, optional
@ -66,6 +91,30 @@ class Model:
if plots is not None:
self.plots = plots
# Store dictionaries to the materials and cells by ID and names
if materials is None:
mats = self.geometry.get_all_materials().values()
else:
mats = self.materials
cells = self.geometry.get_all_cells()
# Get the ID maps
self._materials_by_id = {mat.id: mat for mat in mats}
self._cells_by_id = {cell.id: cell for cell in cells.values()}
# Get the names maps, but since names are not unique, store a list for
# each name key. In this way when the user requests a change by a name,
# the change will be applied to all of the same name.
self._cells_by_name = {}
for cell in cells.values():
if cell.name not in self._cells_by_name:
self._cells_by_name[cell.name] = set()
self._cells_by_name[cell.name].add(cell)
self._materials_by_name = {}
for mat in mats:
if mat.name not in self._materials_by_name:
self._materials_by_name[mat.name] = set()
self._materials_by_name[mat.name].add(mat)
@property
def geometry(self):
return self._geometry
@ -86,6 +135,10 @@ class Model:
def plots(self):
return self._plots
@property
def is_initialized(self):
return openmc.lib.is_initialized
@geometry.setter
def geometry(self, geometry):
check_type('geometry', geometry, openmc.Geometry)
@ -130,6 +183,8 @@ class Model:
def from_xml(cls, geometry='geometry.xml', materials='materials.xml',
settings='settings.xml'):
"""Create model from existing XML files
When initializing this way, the user must manually load plots and
tallies.
Parameters
----------
@ -151,44 +206,159 @@ class Model:
settings = openmc.Settings.from_xml(settings)
return cls(geometry, materials, settings)
def deplete(self, timesteps, chain_file=None, method='cecm',
fission_q=None, **kwargs):
def init_lib(self, threads=None, geometry_debug=False, restart_file=None,
tracks=False, output=True, event_based=None, intracomm=None):
"""Initializes the model in memory via the C API
.. versionadded:: 0.13.0
Parameters
----------
threads : int, optional
Number of OpenMP threads. If OpenMC is compiled with OpenMP
threading enabled, the default is implementation-dependent but is
usually equal to the number of hardware threads available
(or a value set by the :envvar:`OMP_NUM_THREADS` environment
variable).
geometry_debug : bool, optional
Turn on geometry debugging during simulation. Defaults to False.
restart_file : str, optional
Path to restart file to use
tracks : bool, optional
Write tracks for all particles. Defaults to False.
output : bool
Capture OpenMC output from standard out
event_based : None or bool, optional
Turns on event-based parallelism if True. If None, the value in
the Settings will be used.
intracomm : mpi4py.MPI.Intracomm or None, optional
MPI intracommunicator
"""
# TODO: right now the only way to set most of the above parameters via
# the C API are at initialization time despite use-cases existing to
# set them for individual runs. For now this functionality is exposed
# where it exists (here in init), but in the future the functionality
# should be exposed so that it can be accessed via model.run(...)
args = _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()
# The Model object needs to be aware of the communicator so it can
# use it in certain cases, therefore lets store the communicator
if intracomm is not None:
self._intracomm = intracomm
else:
self._intracomm = DummyCommunicator()
if self._intracomm.rank == 0:
self.export_to_xml()
self._intracomm.barrier()
# We cannot pass DummyCommunicator to openmc.lib.init so pass instead
# the user-provided intracomm which will either be None or an mpi4py
# communicator
openmc.lib.init(args=args, intracomm=intracomm, output=output)
def finalize_lib(self):
"""Finalize simulation and free memory allocated for the C API
.. versionadded:: 0.13.0
"""
openmc.lib.finalize()
def deplete(self, timesteps, method='cecm', final_step=True,
operator_kwargs=None, directory='.', output=True,
**integrator_kwargs):
"""Deplete model using specified timesteps/power
.. versionchanged:: 0.13.0
The *final_step*, *operator_kwargs*, *directory*, and *output*
arguments were added.
Parameters
----------
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not
cumulative.
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the chain
found under the ``depletion_chain`` in the
:envvar:`OPENMC_CROSS_SECTIONS` environment variable if it exists.
method : str
Integration method used for depletion (e.g., 'cecm', 'predictor')
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV].
If not given, values will be pulled from the ``chain_file``.
**kwargs
Keyword arguments passed to integration function (e.g.,
:func:`openmc.deplete.integrator.cecm`)
method : str, optional
Integration method used for depletion (e.g., 'cecm', 'predictor').
Defaults to 'cecm'.
final_step : bool, optional
Indicate whether or not a transport solve should be run at the end
of the last timestep. Defaults to running this transport solve.
operator_kwargs : dict
Keyword arguments passed to the depletion Operator initializer
(e.g., :func:`openmc.deplete.Operator`)
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
output : bool
Capture OpenMC output from standard out
integrator_kwargs : dict
Remaining keyword arguments passed to the depletion Integrator
initializer (e.g., :func:`openmc.deplete.integrator.cecm`).
"""
# Import the depletion module. This is done here rather than the module
# header to delay importing openmc.lib (through openmc.deplete) which
# can be tough to install properly.
if operator_kwargs is None:
op_kwargs = {}
elif isinstance(operator_kwargs, dict):
op_kwargs = operator_kwargs
else:
raise ValueError("operator_kwargs must be a dict or None")
# Import openmc.deplete here so the Model can be used even if the
# shared library is unavailable.
import openmc.deplete as dep
# Create OpenMC transport operator
op = dep.Operator(
self.geometry, self.settings, chain_file,
fission_q=fission_q,
)
# Store whether or not the library was initialized when we started
started_initialized = self.is_initialized
# Perform depletion
check_value('method', method, ('cecm', 'predictor', 'cf4', 'epc_rk4',
'si_celi', 'si_leqi', 'celi', 'leqi'))
getattr(dep.integrator, method)(op, timesteps, **kwargs)
with _change_directory(Path(directory)):
with openmc.lib.quiet_dll(output):
depletion_operator = \
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
# Set up the integrator
check_value('method', method,
dep.integrators.integrator_by_name.keys())
integrator_class = dep.integrators.integrator_by_name[method]
integrator = integrator_class(depletion_operator, timesteps,
**integrator_kwargs)
# Now perform the depletion
with openmc.lib.quiet_dll(output):
integrator.integrate(final_step)
# Now make the python Materials match the C API material data
for mat_id, mat in self._materials_by_id.items():
if mat.depletable:
# Get the C data
c_mat = openmc.lib.materials[mat_id]
nuclides, densities = c_mat._get_densities()
# And now we can remove isotopes and add these ones in
mat.nuclides.clear()
for nuc, density in zip(nuclides, densities):
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.
@ -226,6 +396,9 @@ class Model:
def import_properties(self, filename):
"""Import physical properties
.. versionchanged:: 0.13.0
This method now updates values as loaded in memory with the C API
Parameters
----------
filename : str
@ -254,32 +427,69 @@ class Model:
cell = cells[cell_id]
if cell.fill_type in ('material', 'distribmat'):
cell.temperature = group['temperature'][()]
if self.is_initialized:
lib_cell = openmc.lib.cells[cell_id]
lib_cell.set_temperature(group['temperature'][()])
# Make sure number of materials matches
mats_group = fh['materials']
n_cells = mats_group.attrs['n_materials']
if n_cells != len(materials):
raise ValueError("Number of materials in properties file doesn't "
"match current model.")
raise ValueError("Number of materials in properties file "
"doesn't match current model.")
# Update material densities
for name, group in mats_group.items():
mat_id = int(name.split()[1])
atom_density = group.attrs['atom_density']
materials[mat_id].set_density('atom/b-cm', atom_density)
if self.is_initialized:
C_mat = openmc.lib.materials[mat_id]
C_mat.set_density(atom_density, 'atom/b-cm')
def run(self, **kwargs):
"""Creates the XML files, runs OpenMC, and returns the path to the last
def run(self, particles=None, threads=None, geometry_debug=False,
restart_file=None, tracks=False, output=True, cwd='.',
openmc_exec='openmc', mpi_args=None, event_based=None):
"""Runs OpenMC. If the C API has been initialized, then the C API is
used, otherwise, this method creates the XML files and runs OpenMC via
a system call. In both cases this method returns the path to the last
statepoint file generated.
.. versionchanged:: 0.12
Instead of returning the final k-effective value, this function now
returns the path to the final statepoint written.
.. versionchanged:: 0.13.0
This method can utilize the C API for execution
Parameters
----------
**kwargs
Keyword arguments passed to :func:`openmc.run`
particles : int, optional
Number of particles to simulate per generation.
threads : int, optional
Number of OpenMP threads. If OpenMC is compiled with OpenMP
threading enabled, the default is implementation-dependent but is
usually equal to the number of hardware threads available (or a
value set by the :envvar:`OMP_NUM_THREADS` environment variable).
geometry_debug : bool, optional
Turn on geometry debugging during simulation. Defaults to False.
restart_file : str, optional
Path to restart file to use
tracks : bool, optional
Write tracks for all particles. Defaults to False.
output : bool, optional
Capture OpenMC output from standard out
cwd : str, optional
Path to working directory to run in. Defaults to the current
working directory.
openmc_exec : str, optional
Path to OpenMC executable. Defaults to 'openmc'.
mpi_args : list of str, optional
MPI execute command and any additional MPI arguments to pass,
e.g. ['mpiexec', '-n', '8'].
event_based : None or bool, optional
Turns on event-based parallelism if True. If None, the value in
the Settings will be used.
Returns
-------
@ -289,23 +499,364 @@ class Model:
"""
self.export_to_xml()
# Setting tstart here ensures we don't pick up any pre-existing statepoint
# files in the output directory
# Setting tstart here ensures we don't pick up any pre-existing
# statepoint files in the output directory
tstart = time.time()
last_statepoint = None
openmc.run(**kwargs)
# Operate in the provided working directory
with _change_directory(Path(cwd)):
if self.is_initialized:
# Handle the run options as applicable
# First dont allow ones that must be set via init
for arg_name, arg, default in zip(
['threads', 'geometry_debug', 'restart_file', 'tracks'],
[threads, geometry_debug, restart_file, tracks],
[None, False, None, False]):
if arg != default:
msg = f"{arg_name} must be set via Model.is_initialized(...)"
raise ValueError(msg)
# Get output directory and return the last statepoint written by this run
if self.settings.output and 'path' in self.settings.output:
output_dir = Path(self.settings.output['path'])
else:
output_dir = Path.cwd()
for sp in output_dir.glob('statepoint.*.h5'):
mtime = sp.stat().st_mtime
if mtime >= tstart: # >= allows for poor clock resolution
tstart = mtime
last_statepoint = sp
init_particles = openmc.lib.settings.particles
if particles is not None:
if isinstance(particles, Integral) and particles > 0:
openmc.lib.settings.particles = particles
init_event_based = openmc.lib.settings.event_based
if event_based is not None:
openmc.lib.settings.event_based = event_based
# Then run using the C API
openmc.lib.run(output)
# Reset changes for the openmc.run kwargs handling
openmc.lib.settings.particles = init_particles
openmc.lib.settings.event_based = init_event_based
else:
# Then run via the command line
self.export_to_xml()
openmc.run(particles, threads, geometry_debug, restart_file,
tracks, output, Path('.'), openmc_exec, mpi_args,
event_based)
# Get output directory and return the last statepoint written
if self.settings.output and 'path' in self.settings.output:
output_dir = Path(self.settings.output['path'])
else:
output_dir = Path.cwd()
for sp in output_dir.glob('statepoint.*.h5'):
mtime = sp.stat().st_mtime
if mtime >= tstart: # >= allows for poor clock resolution
tstart = mtime
last_statepoint = sp
return last_statepoint
def calculate_volumes(self, threads=None, output=True, cwd='.',
openmc_exec='openmc', mpi_args=None,
apply_volumes=True):
"""Runs an OpenMC stochastic volume calculation and, if requested,
applies volumes to the model
.. versionadded:: 0.13.0
Parameters
----------
threads : int, optional
Number of OpenMP threads. If OpenMC is compiled with OpenMP
threading enabled, the default is implementation-dependent but is
usually equal to the number of hardware threads available (or a
value set by the :envvar:`OMP_NUM_THREADS` environment variable).
This currenty only applies to the case when not using the C API.
output : bool, optional
Capture OpenMC output from standard out
openmc_exec : str, optional
Path to OpenMC executable. Defaults to 'openmc'.
This only applies to the case when not using the C API.
mpi_args : list of str, optional
MPI execute command and any additional MPI arguments to pass,
e.g. ['mpiexec', '-n', '8'].
This only applies to the case when not using the C API.
cwd : str, optional
Path to working directory to run in. Defaults to the current
working directory.
apply_volumes : bool, optional
Whether apply the volume calculation results from this calculation
to the model. Defaults to applying the volumes.
"""
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!")
with _change_directory(Path(cwd)):
if self.is_initialized:
if threads is not None:
msg = "Threads must be set via Model.is_initialized(...)"
raise ValueError(msg)
if mpi_args is not None:
msg = "The MPI environment must be set otherwise such as" \
"with the call to mpi4py"
raise ValueError(msg)
# Compute the volumes
openmc.lib.calculate_volumes(output)
else:
self.export_to_xml()
openmc.calculate_volumes(threads=threads, output=output,
openmc_exec=openmc_exec,
mpi_args=mpi_args)
# Now we apply the volumes
if apply_volumes:
# Load the results and add them to the model
for i, vol_calc in enumerate(self.settings.volume_calculations):
vol_calc.load_results(f"volume_{i + 1}.h5")
# First add them to the Python side
self.geometry.add_volume_information(vol_calc)
# And now repeat for the C API
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.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'):
"""Creates plot images as specified by the Model.plots attribute
If convert is True, this function requires that a program is installed
to convert PPM files to PNG files. Typically, that would be
`ImageMagick <https://www.imagemagick.org>`_ which includes a
`convert` command.
.. versionadded:: 0.13.0
Parameters
----------
output : bool, optional
Capture OpenMC output from standard out
cwd : str, optional
Path to working directory to run in. Defaults to the current
working directory.
openmc_exec : str, optional
Path to OpenMC executable. Defaults to 'openmc'.
This only applies to the case when not using the C API.
convert : bool, optional
Whether or not to attempt to convert from PPM to PNG
convert_exec : str, optional
Command that can convert PPM files into PNG files
"""
if len(self.plots) == 0:
# Then there is no volume calculation specified
raise ValueError("The Model.plots attribute must be specified "
"before executing this method!")
with _change_directory(Path(cwd)):
if self.is_initialized:
# Compute the volumes
openmc.lib.plot_geometry(output)
else:
self.export_to_xml()
openmc.plot_geometry(output=output, openmc_exec=openmc_exec)
if convert:
for p in self.plots:
if p.filename is not None:
ppm_file = f'{p.filename}.ppm'
else:
ppm_file = f'plot_{p.id}.ppm'
png_file = ppm_file.replace('.ppm', '.png')
subprocess.check_call([convert_exec, ppm_file, png_file])
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))
check_type('obj_type', obj_type, str)
obj_type = obj_type.lower()
check_value('obj_type', obj_type, ('material', 'cell'))
check_value('attrib_name', attrib_name,
('temperature', 'volume', 'density', 'rotation',
'translation'))
# The C API only allows setting density units of atom/b-cm and g/cm3
check_value('density_units', density_units, ('atom/b-cm', 'g/cm3'))
# The C API has no way to set cell volume or material temperature
# so lets raise exceptions as needed
if obj_type == 'cell' and attrib_name == 'volume':
raise NotImplementedError(
'Setting a Cell volume is not supported!')
if obj_type == 'material' and attrib_name == 'temperature':
raise NotImplementedError(
'Setting a material temperature is not supported!')
# And some items just dont make sense
if obj_type == 'cell' and attrib_name == 'density':
raise ValueError('Cannot set a Cell density!')
if obj_type == 'material' and attrib_name in ('rotation',
'translation'):
raise ValueError('Cannot set a material rotation/translation!')
# Set the
if obj_type == 'cell':
by_name = self._cells_by_name
by_id = self._cells_by_id
obj_by_id = openmc.lib.cells
else:
by_name = self._materials_by_name
by_id = self._materials_by_id
obj_by_id = openmc.lib.materials
# Get the list of ids to use if converting from names and accepting
# only values that have actual ids
ids = []
for name_or_id in names_or_ids:
if isinstance(name_or_id, Integral):
if name_or_id in by_id:
ids.append(int(name_or_id))
else:
cap_obj = obj_type.capitalize()
msg = f'{cap_obj} ID {name_or_id} " \
"is not present in the model!'
raise InvalidIDError(msg)
elif isinstance(name_or_id, str):
if name_or_id in by_name:
# Then by_name[name_or_id] is a list so we need to add all
# entries
ids.extend([obj.id for obj in by_name[name_or_id]])
else:
cap_obj = obj_type.capitalize()
msg = f'{cap_obj} {name_or_id} " \
"is not present in the model!'
raise InvalidIDError(msg)
# Now perform the change to both python and C API
for id_ in ids:
obj = by_id[id_]
if attrib_name == 'density':
obj.set_density(density_units, value)
else:
setattr(obj, attrib_name, value)
# Next lets keep what is in C API memory up to date as well
if self.is_initialized:
lib_obj = obj_by_id[id_]
if attrib_name == 'density':
lib_obj.set_density(value, density_units)
elif attrib_name == 'temperature':
lib_obj.set_temperature(value)
else:
setattr(lib_obj, attrib_name, value)
def rotate_cells(self, names_or_ids, vector):
"""Rotate the identified cell(s) by the specified rotation vector.
The rotation is only applied to cells filled with a universe.
.. note:: If applying this change to a name that is not unique, then
the change will be applied to all objects of that name.
.. versionadded:: 0.13.0
Parameters
----------
names_or_ids : Iterable of str or int
The cell names (if str) or id (if int) that are to be translated
or rotated. This parameter can include a mix of names and ids.
vector : Iterable of float
The rotation vector of length 3 to apply. This array specifies the
angles in degrees about the x, y, and z axes, respectively.
"""
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.
The translation is only applied to cells filled with a universe.
.. note:: If applying this change to a name that is not unique, then
the change will be applied to all objects of that name.
.. versionadded:: 0.13.0
Parameters
----------
names_or_ids : Iterable of str or int
The cell names (if str) or id (if int) that are to be translated
or rotated. This parameter can include a mix of names and ids.
vector : Iterable of float
The translation vector of length 3 to apply. This array specifies
the x, y, and z dimensions of the 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
.. note:: If applying this change to a name that is not unique, then
the change will be applied to all objects of that name.
.. versionadded:: 0.13.0
Parameters
----------
names_or_ids : Iterable of str or int
The material names (if str) or id (if int) that are to be updated.
This parameter can include a mix of names and ids.
density : float
The density to apply in the units specified by `density_units`
density_units : {'atom/b-cm', 'g/cm3'}, optional
Units for `density`. Defaults to 'atom/b-cm'
"""
self._change_py_lib_attribs(names_or_ids, density, 'material',
'density', density_units)
def update_cell_temperatures(self, names_or_ids, temperature):
"""Update the temperature of a set of cells to the given value
.. note:: If applying this change to a name that is not unique, then
the change will be applied to all objects of that name.
.. versionadded:: 0.13.0
Parameters
----------
names_or_ids : Iterable of str or int
The cell names (if str) or id (if int) that are to be updated.
This parameter can include a mix of names and ids.
temperature : float
The temperature to apply in units of Kelvin
"""
self._change_py_lib_attribs(names_or_ids, temperature, 'cell',
'temperature')
def update_material_volumes(self, names_or_ids, volume):
"""Update the volume of a set of materials to the given value
.. note:: If applying this change to a name that is not unique, then
the change will be applied to all objects of that name.
.. versionadded:: 0.13.0
Parameters
----------
names_or_ids : Iterable of str or int
The material names (if str) or id (if int) that are to be updated.
This parameter can include a mix of names and ids.
volume : float
The volume to apply in units of cm^3
"""
self._change_py_lib_attribs(names_or_ids, volume, 'material', 'volume')

8
openmc/mpi.py Normal file
View file

@ -0,0 +1,8 @@
try:
from mpi4py import MPI
comm = MPI.COMM_WORLD
except ImportError:
from unittest.mock import Mock
MPI = Mock()
from openmc.dummy_comm import DummyCommunicator
comm = DummyCommunicator()

3
src/finalize.cpp
View file

@ -15,6 +15,7 @@
#include "openmc/message_passing.h"
#include "openmc/nuclide.h"
#include "openmc/photon.h"
#include "openmc/plot.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
@ -45,6 +46,7 @@ void free_memory()
free_memory_mesh();
free_memory_tally();
free_memory_bank();
free_memory_plot();
if (mpi::master) {
free_memory_cmfd();
}
@ -128,6 +130,7 @@ int openmc_finalize()
data::temperature_min = 0.0;
data::temperature_max = INFTY;
model::root_universe = -1;
model::plotter_seed = 1;
openmc::openmc_set_seed(DEFAULT_SEED);
// Deallocate arrays

4
src/initialize.cpp
View file

@ -305,9 +305,11 @@ void read_input_xml()
// Initialize distribcell_filters
prepare_distribcell();
// Read the plots.xml regardless of plot mode in case plots are requested
// via the API
read_plots_xml();
if (settings::run_mode == RunMode::PLOTTING) {
// Read plots.xml if it exists
read_plots_xml();
if (mpi::master && settings::verbosity >= 5)
print_plot();

12
src/plot.cpp
View file

@ -119,9 +119,11 @@ extern "C" int openmc_plot_geometry()
void read_plots_xml()
{
// Check if plots.xml exists
// Check if plots.xml exists; this is only necessary when the plot runmode is
// initiated. Otherwise, we want to read plots.xml because it may be called
// later via the API. In that case, its ok for a plots.xml to not exist
std::string filename = settings::path_input + "plots.xml";
if (!file_exists(filename)) {
if (!file_exists(filename) && settings::run_mode == RunMode::PLOTTING) {
fatal_error(fmt::format("Plots XML file '{}' does not exist!", filename));
}
@ -138,6 +140,12 @@ void read_plots_xml()
}
}
void free_memory_plot()
{
model::plots.clear();
model::plot_map.clear();
}
//==============================================================================
// CREATE_PPM creates an image based on user input from a plots.xml <plot>
// specification in the portable pixmap format (PPM)

3
tests/unit_tests/test_deplete_chain.py
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.mpi import comm
from openmc.deplete import Chain, reaction_rates, nuclide, cram, pool
import pytest
from tests import cdtemp

8
tests/unit_tests/test_deplete_integrator.py
View file

@ -14,11 +14,11 @@ import numpy as np
from uncertainties import ufloat
import pytest
from openmc.mpi 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

8
tests/unit_tests/test_lib.py
View file

@ -263,6 +263,7 @@ def test_settings(lib_init):
assert settings.generations_per_batch == 1
assert settings.particles == 100
assert settings.seed == 1
assert settings.event_based is False
settings.seed = 11
@ -713,7 +714,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 +727,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 +745,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()

523
tests/unit_tests/test_model.py
View file

@ -1,29 +1,281 @@
from math import pi
from pathlib import Path
from shutil import which
import numpy as np
import pytest
import openmc
import openmc.lib
@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):
mats, geom, settings, tals, plots, _, _ = \
pin_model_attributes
openmc.reset_auto_ids()
# Check blank initialization of a model
test_model = openmc.Model()
assert test_model.geometry.root_universe is None
assert len(test_model.materials) == 0
ref_settings = openmc.Settings()
assert sorted(test_model.settings.__dict__.keys()) == \
sorted(ref_settings.__dict__.keys())
for ref_k, ref_v in ref_settings.__dict__.items():
assert test_model.settings.__dict__[ref_k] == ref_v
assert len(test_model.tallies) == 0
assert len(test_model.plots) == 0
assert test_model._materials_by_id == {}
assert test_model._materials_by_name == {}
assert test_model._cells_by_id == {}
assert test_model._cells_by_name == {}
assert test_model.is_initialized is False
# Now check proper init of an actual model. Assume no interference between
# parameters and so we can apply them all at once instead of testing one
# parameter initialization at a time
test_model = openmc.Model(geom, mats, settings, tals, plots)
assert test_model.geometry is geom
assert test_model.materials is mats
assert test_model.settings is settings
assert test_model.tallies is tals
assert test_model.plots is plots
assert test_model._materials_by_id == {1: mats[0], 2: mats[1], 3: mats[2]}
assert test_model._materials_by_name == {
'UO2': {mats[0]}, 'Zirc': {mats[1]}, 'Borated water': {mats[2]}}
# The last cell is the one that contains the infinite fuel
assert test_model._cells_by_id == \
{2: geom.root_universe.cells[2], 3: geom.root_universe.cells[3],
4: geom.root_universe.cells[4],
1: geom.root_universe.cells[2].fill.cells[1]}
# No cell name for 2 and 3, so we expect a blank name to be assigned to
# cell 3 due to overwriting
assert test_model._cells_by_name == {
'fuel': {geom.root_universe.cells[2]},
'': {geom.root_universe.cells[3], geom.root_universe.cells[4]},
'inf fuel': {geom.root_universe.cells[2].fill.cells[1]}}
assert test_model.is_initialized is False
# Finally test the parameter type checking by passing bad types and
# obtaining the right exception types
def_params = [geom, mats, settings, tals, plots]
for i in range(len(def_params)):
args = def_params.copy()
# Try an integer, as that is a bad type for all arguments
args[i] = i
with pytest.raises(TypeError):
test_model = openmc.Model(*args)
def test_from_xml(run_in_tmpdir, pin_model_attributes):
mats, geom, settings, tals, plots, _, _ = pin_model_attributes
# This test will write the individual files to xml and then init that way
# and run the same sort of test as in test_init
mats.export_to_xml()
geom.export_to_xml()
settings.export_to_xml()
tals.export_to_xml()
plots.export_to_xml()
# This from_xml method cannot load chain and fission_q
test_model = openmc.Model.from_xml()
assert test_model.geometry.root_universe.cells.keys() == \
geom.root_universe.cells.keys()
assert [c.fill.name for c in
test_model.geometry.root_universe.cells.values()] == \
[c.fill.name for c in geom.root_universe.cells.values()]
assert [mat.name for mat in test_model.materials] == \
[mat.name for mat in mats]
# We will assume the attributes of settings that are custom objects are
# OK if the others are so we dotn need to implement explicit comparisons
no_test = ['_source', '_entropy_mesh']
assert sorted(k for k in test_model.settings.__dict__.keys()
if k not in no_test) == \
sorted(k for k in settings.__dict__.keys() if k not in no_test)
keys = sorted(k for k in settings.__dict__.keys() if k not in no_test)
for ref_k in keys:
assert test_model.settings.__dict__[ref_k] == settings.__dict__[ref_k]
assert len(test_model.tallies) == 0
assert len(test_model.plots) == 0
assert test_model._materials_by_id == \
{1: test_model.materials[0], 2: test_model.materials[1],
3: test_model.materials[2]}
assert test_model._materials_by_name == {
'UO2': {test_model.materials[0]}, 'Zirc': {test_model.materials[1]},
'Borated water': {test_model.materials[2]}}
assert test_model._cells_by_id == {
2: test_model.geometry.root_universe.cells[2],
3: test_model.geometry.root_universe.cells[3],
4: test_model.geometry.root_universe.cells[4],
1: test_model.geometry.root_universe.cells[2].fill.cells[1]}
# No cell name for 2 and 3, so we expect a blank name to be assigned to
# cell 3 due to overwriting
assert test_model._cells_by_name == {
'fuel': {test_model.geometry.root_universe.cells[2]},
'': {test_model.geometry.root_universe.cells[3],
test_model.geometry.root_universe.cells[4]},
'inf fuel': {test_model.geometry.root_universe.cells[2].fill.cells[1]}}
assert test_model.is_initialized is False
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)
test_model.init_lib(output=False, intracomm=mpi_intracomm)
# First check that the API is advertised as initialized
assert openmc.lib.is_initialized is True
assert test_model.is_initialized is True
# Now make sure it actually is initialized by making a call to the lib
c_mat = openmc.lib.find_material((0.6, 0., 0.))
# This should be Borated water
assert c_mat.name == 'Borated water'
assert c_mat.id == 3
# Ok, now lets test that we can clear the data and check that it is cleared
test_model.finalize_lib()
# First check that the API is advertised as initialized
assert openmc.lib.is_initialized is False
assert test_model.is_initialized is False
# Note we cant actually test that a sys call fails because we should get a
# seg fault
def test_import_properties(run_in_tmpdir, mpi_intracomm):
"""Test importing properties on the Model class """
# Create PWR pin cell model and write XML files
openmc.reset_auto_ids()
model = openmc.examples.pwr_pin_cell()
model.export_to_xml()
model.init_lib(output=False, intracomm=mpi_intracomm)
# Change fuel temperature and density and export properties
openmc.lib.init(intracomm=mpi_intracomm)
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")
# Check to see that values are assigned to the C and python representations
# First python
cell = model.geometry.get_all_cells()[1]
assert cell.temperature == [600.0]
assert cell.fill.get_mass_density() == pytest.approx(5.0)
# Now C
assert openmc.lib.cells[1].get_temperature() == 600.
assert openmc.lib.materials[1].get_density('g/cm3') == pytest.approx(5.0)
# Clear the C API
openmc.lib.finalize()
# Import properties to existing model and re-export to new directory
model.import_properties("properties.h5")
# Verify the attributes survived by exporting to XML and re-creating
model.export_to_xml("with_properties")
# Load model with properties and confirm temperature/density has been changed
# Load model with properties and confirm temperature/density changed
model_with_properties = openmc.Model.from_xml(
'with_properties/geometry.xml',
'with_properties/materials.xml',
@ -32,3 +284,262 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm):
cell = model_with_properties.geometry.get_all_cells()[1]
assert cell.temperature == [600.0]
assert cell.fill.get_mass_density() == pytest.approx(5.0)
def test_run(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)
# This case will run by getting the k-eff and tallies for command-line and
# C API execution modes and ensuring they give the same result.
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(scores=['flux'])[0, 0, 0]
cli_fiss = sp.get_tally(id=1).get_values(scores=['fission'])[0, 0, 0]
test_model.init_lib(output=False, intracomm=mpi_intracomm)
sp_path = test_model.run(output=False)
with openmc.StatePoint(sp_path) as sp:
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 lib_keff.n == pytest.approx(cli_keff.n, abs=1e-13)
assert lib_flux == pytest.approx(cli_flux, abs=1e-13)
assert lib_fiss == pytest.approx(cli_fiss, abs=1e-13)
# Now we should make sure that the flags for items which should be handled
# by init are properly set
with pytest.raises(ValueError):
test_model.run(threads=1)
with pytest.raises(ValueError):
test_model.run(geometry_debug=True)
with pytest.raises(ValueError):
test_model.run(restart_file='1.h5')
with pytest.raises(ValueError):
test_model.run(tracks=True)
test_model.finalize_lib()
def test_plots(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)
# This test cannot check the correctness of the plot, but it can
# check that a plot was made and that the expected ppm and png files are
# there
# We will only test convert if it is on the system, so as not to add an
# extra dependency just for tests
convert = which('convert') is not None
if convert:
exts = ['ppm', 'png']
else:
exts = ['ppm']
# 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(output=False, intracomm=mpi_intracomm)
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
for fname in ['test.', 'plot_2.']:
for ext in exts:
test_file = Path(f'./{fname}{ext}')
assert test_file.exists()
test_file.unlink()
test_model.finalize_lib()
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)
test_model.init_lib(output=False, intracomm=mpi_intracomm)
# Now we can call rotate_cells, translate_cells, update_densities,
# and update_cell_temperatures and make sure the changes have taken hold.
# For each we will first try bad inputs to make sure we get the right
# errors and then we do a good one which calls the material by name and
# then id to make sure it worked
# The rotate_cells and translate_cells will work on the cell named fill, as
# it is filled with a universe and thus the operation will be valid
# First rotate_cells
with pytest.raises(TypeError):
# Make sure it tells us we have a bad names_or_ids type
test_model.rotate_cells(None, (0, 0, 90))
with pytest.raises(TypeError):
test_model.rotate_cells([None], (0, 0, 90))
with pytest.raises(openmc.exceptions.InvalidIDError):
# Make sure it tells us we had a bad id
test_model.rotate_cells([7200], (0, 0, 90))
with pytest.raises(openmc.exceptions.InvalidIDError):
# Make sure it tells us we had a bad id
test_model.rotate_cells(['bad_name'], (0, 0, 90))
# Now a good one
assert np.all(openmc.lib.cells[2].rotation == (0., 0., 0.))
test_model.rotate_cells([2], (0, 0, 90))
assert np.all(openmc.lib.cells[2].rotation == (0., 0., 90.))
# And same thing by name
test_model.rotate_cells(['fuel'], (0, 0, 180))
# Now translate_cells. We dont need to re-check the TypeErrors/bad ids,
# because the other functions use the same hidden method as rotate_cells
assert np.all(openmc.lib.cells[2].translation == (0., 0., 0.))
test_model.translate_cells([2], (0, 0, 10))
assert np.all(openmc.lib.cells[2].translation == (0., 0., 10.))
# Now lets do the density updates.
# Check initial conditions
assert openmc.lib.materials[1].get_density('atom/b-cm') == \
pytest.approx(0.06891296988603757, abs=1e-13)
mat_a_dens = np.sum(
[v[1] for v in test_model.materials[0].
get_nuclide_atom_densities().values()])
assert mat_a_dens == pytest.approx(0.06891296988603757, abs=1e-8)
# Change the density
test_model.update_densities(['UO2'], 2.)
assert openmc.lib.materials[1].get_density('atom/b-cm') == \
pytest.approx(2., abs=1e-13)
mat_a_dens = np.sum(
[v[1] for v in test_model.materials[0].
get_nuclide_atom_densities().values()])
assert mat_a_dens == pytest.approx(2., abs=1e-8)
# Now lets do the cell temperature updates.
# Check initial conditions
assert test_model._cells_by_id == \
{2: geom.root_universe.cells[2], 3: geom.root_universe.cells[3],
4: geom.root_universe.cells[4],
1: geom.root_universe.cells[2].fill.cells[1]}
assert openmc.lib.cells[3].get_temperature() == \
pytest.approx(293.6, abs=1e-13)
assert test_model.geometry.root_universe.cells[3].temperature is None
# Change the temperature
test_model.update_cell_temperatures([3], 600.)
assert openmc.lib.cells[3].get_temperature() == \
pytest.approx(600., abs=1e-13)
assert test_model.geometry.root_universe.cells[3].temperature == \
pytest.approx(600., abs=1e-13)
# And finally material volume
assert openmc.lib.materials[1].volume == \
pytest.approx(0.4831931368640985, abs=1e-13)
# The temperature on the material will be None because its just the default
assert test_model.materials[0].volume == \
pytest.approx(0.4831931368640985, abs=1e-13)
# Change the temperature
test_model.update_material_volumes(['UO2'], 2.)
assert openmc.lib.materials[1].volume == pytest.approx(2., abs=1e-13)
assert test_model.materials[0].volume == pytest.approx(2., abs=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)
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., output=False)
# 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 after_xe + after_u == pytest.approx(initial_u, abs=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, intracomm=mpi_intracomm)
test_model.deplete([1e6], 'predictor', final_step=False,
operator_kwargs=op_kwargs,
power=1., output=False)
# 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 after_lib_xe + after_lib_u == pytest.approx(initial_u, abs=1e-15)
assert test_model.is_initialized is True
# And end by comparing to the previous case
assert after_xe == pytest.approx(after_lib_xe, abs=1e-15)
assert after_u == pytest.approx(after_lib_u, abs=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)
# 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, intracomm=mpi_intracomm)
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()