diff --git a/LICENSE b/LICENSE index 848a23845..3523dd370 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,5 @@ -Copyright (c) 2011-2021 Massachusetts Institute of Technology and OpenMC contributors +Copyright (c) 2011-2021 Massachusetts Institute of Technology, UChicago Argonne +LLC, and OpenMC contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/docs/source/conf.py b/docs/source/conf.py index 52d872abf..a10494260 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -62,7 +62,7 @@ master_doc = 'index' # General information about the project. project = 'OpenMC' -copyright = '2011-2021, Massachusetts Institute of Technology and OpenMC contributors' +copyright = '2011-2021, Massachusetts Institute of Technology, UChicago Argonne LLC, and OpenMC contributors' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the diff --git a/docs/source/license.rst b/docs/source/license.rst index 40ab106ad..d8ea319aa 100644 --- a/docs/source/license.rst +++ b/docs/source/license.rst @@ -4,7 +4,8 @@ License Agreement ================= -Copyright © 2011-2021 Massachusetts Institute of Technology and OpenMC contributors +Copyright © 2011-2021 Massachusetts Institute of Technology, UChicago Argonne +LLC, and OpenMC contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 93a5f42ad..650b7e16a 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -30,9 +30,7 @@ namespace model { extern std::unordered_map plot_map; //!< map of plot ids to index extern vector 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 @@ -281,6 +279,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 an image for a plot object //! \param[in] plot object void create_image(Plot const& pl); diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 2d0261960..f14f6c15f 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -32,7 +32,8 @@ 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? diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index 8a93af514..f99389069 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -57,8 +57,8 @@ Indicates the default path to an HDF5 file that contains multi-group cross section libraries if the user has not specified the tag in .I materials.xml\fP. .SH LICENSE -Copyright \(co 2011-2021 Massachusetts Institute of Technology and OpenMC -contributors. +Copyright \(co 2011-2021 Massachusetts Institute of Technology, UChicago +Argonne LLC, and OpenMC contributors. .PP Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/openmc/__init__.py b/openmc/__init__.py index 8985099ba..d118e6a05 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -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' diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 00b29e319..b8c1cdfff 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -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 * diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index ca16926ea..b9198696d 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -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 diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 6c3badc09..f22382923 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -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 ( diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index fbf43a808..585deabba 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -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} diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 0e5cc9efa..29370eb8c 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -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.""" diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index ce8331e82..111c8e638 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -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) diff --git a/openmc/deplete/dummy_comm.py b/openmc/dummy_comm.py similarity index 100% rename from openmc/deplete/dummy_comm.py rename to openmc/dummy_comm.py diff --git a/openmc/executor.py b/openmc/executor.py index 882764730..b64d50125 100644 --- a/openmc/executor.py +++ b/openmc/executor.py @@ -6,6 +6,83 @@ import openmc from .plots import _get_plot_image +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, @@ -116,8 +193,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 @@ -140,12 +217,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) @@ -161,8 +235,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. @@ -191,27 +265,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) diff --git a/openmc/lib/__init__.py b/openmc/lib/__init__.py index 82dc92ba4..c14b0d9c2 100644 --- a/openmc/lib/__init__.py +++ b/openmc/lib/__init__.py @@ -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 diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 0da24aab8..e277d527a 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -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) diff --git a/openmc/lib/material.py b/openmc/lib/material.py index 6770721ab..fde197d3d 100644 --- a/openmc/lib/material.py +++ b/openmc/lib/material.py @@ -172,7 +172,6 @@ class Material(_FortranObjectWithID): @property def nuclides(self): return self._get_densities()[0] - return nuclides @property def densities(self): diff --git a/openmc/lib/settings.py b/openmc/lib/settings.py index 1eab4aa8a..2e9dd18df 100644 --- a/openmc/lib/settings.py +++ b/openmc/lib/settings.py @@ -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): diff --git a/openmc/model/model.py b/openmc/model/model.py index d6c91e73e..ad898fd0a 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -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 `_ 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') diff --git a/openmc/mpi.py b/openmc/mpi.py new file mode 100644 index 000000000..cfc10c0d1 --- /dev/null +++ b/openmc/mpi.py @@ -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() diff --git a/src/finalize.cpp b/src/finalize.cpp index b1e09a37e..2f6c65aaf 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -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 diff --git a/src/initialize.cpp b/src/initialize.cpp index e53767691..75ec20518 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -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(); diff --git a/src/output.cpp b/src/output.cpp index 73f3a0e9a..7bc012296 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -72,26 +72,26 @@ void title() // Write version information fmt::print( - " | The OpenMC Monte Carlo Code\n" - " Copyright | 2011-2021 MIT and OpenMC contributors\n" - " License | https://docs.openmc.org/en/latest/license.html\n" - " Version | {}.{}.{}{}\n", + " | The OpenMC Monte Carlo Code\n" + " Copyright | 2011-2021 MIT, UChicago Argonne LLC, and contributors\n" + " License | https://docs.openmc.org/en/latest/license.html\n" + " Version | {}.{}.{}{}\n", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE, VERSION_DEV ? "-dev" : ""); #ifdef GIT_SHA1 - fmt::print(" Git SHA1 | {}\n", GIT_SHA1); + fmt::print(" Git SHA1 | {}\n", GIT_SHA1); #endif // Write the date and time - fmt::print(" Date/Time | {}\n", time_stamp()); + fmt::print(" Date/Time | {}\n", time_stamp()); #ifdef OPENMC_MPI // Write number of processors - fmt::print(" MPI Processes | {}\n", mpi::n_procs); + fmt::print(" MPI Processes | {}\n", mpi::n_procs); #endif #ifdef _OPENMP // Write number of OpenMP threads - fmt::print(" OpenMP Threads | {}\n", omp_get_max_threads()); + fmt::print(" OpenMP Threads | {}\n", omp_get_max_threads()); #endif std::cout << std::endl; } @@ -335,8 +335,8 @@ void print_version() #ifdef GIT_SHA1 fmt::print("Git SHA1: {}\n", GIT_SHA1); #endif - fmt::print("Copyright (c) 2011-2021 Massachusetts Institute of " - "Technology and OpenMC contributors\nMIT/X license at " + fmt::print("Copyright (c) 2011-2021 MIT, UChicago Argonne LLC, and " + "contributors\nMIT/X license at " "\n"); } } diff --git a/src/plot.cpp b/src/plot.cpp index 27a01bc07..918300e09 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -125,9 +125,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)); } @@ -144,6 +146,12 @@ void read_plots_xml() } } +void free_memory_plot() +{ + model::plots.clear(); + model::plot_map.clear(); +} + //============================================================================== // CREATE_IMAGE creates an image based on user input from a plots.xml // specification in the PNG/PPM format diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index fa90b35f2..9769a49e3 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -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 diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 3195c5976..a08d8738c 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -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 diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 41613832f..ccff33ac3 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -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() diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 4b658d5ab..b1021a269 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -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 = """ + + + + + + 2.53000e-02 + + Xe136 + 1.0 + + + + +""" + 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()