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

Merge branch 'develop' into mixed_ncrystal_pr

Browse source
This commit is contained in:
Paul Romano 2023-01-09 13:13:11 +07:00
commit b6d6bdf3f5
166 changed files with 7699 additions and 4693 deletions
Download patch file Download diff file Expand all files Collapse all files

33
CMakeLists.txt
View file

@ -3,8 +3,8 @@ project(openmc C CXX)
# Set version numbers
set(OPENMC_VERSION_MAJOR 0)
set(OPENMC_VERSION_MINOR 14)
set(OPENMC_VERSION_RELEASE 0)
set(OPENMC_VERSION_MINOR 13)
set(OPENMC_VERSION_RELEASE 3)
set(OPENMC_VERSION ${OPENMC_VERSION_MAJOR}.${OPENMC_VERSION_MINOR}.${OPENMC_VERSION_RELEASE})
configure_file(include/openmc/version.h.in "${CMAKE_BINARY_DIR}/include/openmc/version.h" @ONLY)
@ -36,6 +36,7 @@ option(OPENMC_ENABLE_COVERAGE "Compile with coverage analysis flags"
option(OPENMC_USE_DAGMC "Enable support for DAGMC (CAD) geometry" OFF)
option(OPENMC_USE_LIBMESH "Enable support for libMesh unstructured mesh tallies" OFF)
option(OPENMC_USE_MPI "Enable MPI" OFF)
option(OPENMC_USE_MCPL "Enable MCPL" OFF)
option(OPENMC_USE_NCRYSTAL "Enable support for NCrystal scattering" OFF)
# Warnings for deprecated options
@ -168,6 +169,15 @@ if(${HDF5_VERSION} VERSION_GREATER_EQUAL 1.12.0)
list(APPEND cxxflags -DH5Oget_info_by_idx_vers=1 -DH5O_info_t_vers=1)
endif()
#===============================================================================
# MCPL
#===============================================================================
if (OPENMC_USE_MCPL)
find_package(MCPL REQUIRED)
message(STATUS "Found MCPL: ${MCPL_DIR} (found version \"${MCPL_VERSION}\")")
endif()
#===============================================================================
# Set compile/link flags based on which compiler is being used
#===============================================================================
@ -333,6 +343,7 @@ list(APPEND libopenmc_SOURCES
src/lattice.cpp
src/material.cpp
src/math_functions.cpp
src/mcpl_interface.cpp
src/mesh.cpp
src/message_passing.cpp
src/mgxs.cpp
@ -492,11 +503,29 @@ if (OPENMC_USE_MPI)
target_link_libraries(libopenmc MPI::MPI_CXX)
endif()
if (OPENMC_USE_MCPL)
target_compile_definitions(libopenmc PUBLIC OPENMC_MCPL)
target_link_libraries(libopenmc MCPL::mcpl)
endif()
if(OPENMC_USE_NCRYSTAL)
target_compile_definitions(libopenmc PRIVATE NCRYSTAL)
target_link_libraries(libopenmc NCrystal::NCrystal)
endif()
#===============================================================================
# Log build info that this executable can report later
#===============================================================================
target_compile_definitions(libopenmc PRIVATE BUILD_TYPE=${CMAKE_BUILD_TYPE})
target_compile_definitions(libopenmc PRIVATE COMPILER_ID=${CMAKE_CXX_COMPILER_ID})
target_compile_definitions(libopenmc PRIVATE COMPILER_VERSION=${CMAKE_CXX_COMPILER_VERSION})
if (OPENMC_ENABLE_PROFILE)
target_compile_definitions(libopenmc PRIVATE PROFILINGBUILD)
endif()
if (OPENMC_ENABLE_COVERAGE)
target_compile_definitions(libopenmc PRIVATE COVERAGEBUILD)
endif()
#===============================================================================
# openmc executable
#===============================================================================

4
cmake/OpenMCConfig.cmake.in
View file

@ -30,3 +30,7 @@ endif()
if(@OPENMC_USE_MPI@)
find_package(MPI REQUIRED)
endif()
if(@OPENMC_USE_MCPL@)
find_package(MCPL REQUIRED)
endif()

8
docs/source/conf.py
View file

@ -69,9 +69,9 @@ copyright = '2011-2022, Massachusetts Institute of Technology, UChicago Argonne
# built documents.
#
# The short X.Y version.
version = "0.14"
version = "0.13"
# The full version, including alpha/beta/rc tags.
release = "0.14.0"
release = "0.13.3"
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
@ -247,7 +247,7 @@ napoleon_use_ivar = True
intersphinx_mapping = {
'python': ('https://docs.python.org/3', None),
'numpy': ('https://numpy.org/doc/stable/', None),
'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
'scipy': ('https://docs.scipy.org/doc/scipy/', None),
'pandas': ('https://pandas.pydata.org/pandas-docs/stable/', None),
'matplotlib': ('https://matplotlib.org/', None)
'matplotlib': ('https://matplotlib.org/stable/', None)
}

8
docs/source/devguide/tests.rst
View file

@ -31,9 +31,13 @@ Prerequisites
- The test suite requires a specific set of cross section data in order for
tests to pass. A download URL for the data that OpenMC expects can be found
within ``tools/ci/download-xs.sh``.
within ``tools/ci/download-xs.sh``. Once the tarball is downloaded and
unpacked, set the :envvar:`OPENMC_CROSS_SECTIONS` environment variable to the
path of the ``cross_sections.xml`` file within the unpacked data.
- In addition to the HDF5 data, some tests rely on ENDF files. A download URL
for those can also be found in ``tools/ci/download-xs.sh``.
for those can also be found in ``tools/ci/download-xs.sh``. Once the tarball
is downloaded and unpacked, set the :envvar:`OPENMC_ENDF_DATA` environment
variable to the top-level directory of the unpacked tarball.
- Some tests require `NJOY <https://www.njoy21.io/NJOY2016>`_ to preprocess
cross section data. The test suite assumes that you have an ``njoy``
executable available on your :envvar:`PATH`.

29
docs/source/io_formats/settings.rst
View file

@ -732,6 +732,14 @@ attributes/sub-elements:
*Default*: false
:mcpl:
If this element is set to "true", the source point file containing the
source bank will be written as an MCPL_ file name ``source.mcpl`` instead of
an HDF5 file. This option is only applicable if the ``<separate>`` element
is set to true.
*Default*: false
------------------------------
``<surf_source_read>`` Element
------------------------------
@ -767,6 +775,16 @@ certain surfaces and write out the source bank in a separate file called
*Default*: None
:mcpl:
An optional boolean which indicates if the banked particles should be
written to a file in the MCPL_-format instead of the native HDF5-based
format. If activated the output file name is changed to
``surface_source.mcpl``.
*Default*: false
.. _MCPL: https://mctools.github.io/mcpl/mcpl.pdf
------------------------------
``<survival_biasing>`` Element
------------------------------
@ -832,7 +850,9 @@ cell, the nearest temperature at which cross sections are given is to be
applied, within a given tolerance (see :ref:`temperature_tolerance`). A value of
"interpolation" indicates that cross sections are to be linear-linear
interpolated between temperatures at which nuclear data are present (see
:ref:`temperature_treatment`).
:ref:`temperature_treatment`). With the "interpolation" method, temperatures
outside of the bounds of the nuclear data may be accepted, provided they still
fall within the tolerance (see :ref:`temperature_tolerance`).
*Default*: "nearest"
@ -871,7 +891,12 @@ The ``<temperature_tolerance>`` element specifies a tolerance in Kelvin that is
to be applied when the "nearest" temperature method is used. For example, if a
cell temperature is 340 K and the tolerance is 15 K, then the closest
temperature in the range of 325 K to 355 K will be used to evaluate cross
sections.
sections. If the ``<temperature_method>`` is "interpolation", the tolerance
specified applies to cell temperatures outside of the data bounds. For example,
if a cell is specified at 695K, a tolerance of 15K and data is only available
at 700K and 1000K, the cell's cross sections will be evaluated at 700K, since
the desired temperature of 695K is within the tolerance of the actual data
despite not being bounded on both sides.
*Default*: 10 K

3
docs/source/io_formats/statepoint.rst
View file

@ -109,7 +109,8 @@ The current version of the statepoint file format is 17.0.
- **y** (*double[]*) -- Interpolant values for energyfunction
interpolation. Only used for 'energyfunction' filters.
:Attributes: - **interpolation** (*int*) -- Interpolation type. Only used for
:Attributes:
- **interpolation** (*int*) -- Interpolation type. Only used for
'energyfunction' filters.
**/tallies/derivatives/derivative <id>/**

2
docs/source/pythonapi/base.rst
View file

@ -120,7 +120,7 @@ Constructing Tallies
openmc.MaterialFilter
openmc.CellFilter
openmc.CellFromFilter
openmc.CellbornFilter
openmc.CellBornFilter
openmc.CellInstanceFilter
openmc.CollisionFilter
openmc.SurfaceFilter

3
docs/source/pythonapi/data.rst
View file

@ -63,9 +63,10 @@ Core Functions
atomic_weight
combine_distributions
decay_constant
decay_energy
decay_photon_energy
dose_coefficients
gnd_name
gnds_name
half_life
isotopes
kalbach_slope

1
docs/source/pythonapi/model.rst
View file

@ -31,6 +31,7 @@ Composite Surfaces
openmc.model.XConeOneSided
openmc.model.YConeOneSided
openmc.model.ZConeOneSided
openmc.model.Polygon
TRISO Fuel Modeling
-------------------

6
docs/source/pythonapi/stats.rst
View file

@ -22,6 +22,12 @@ Univariate Probability Distributions
openmc.stats.Legendre
openmc.stats.Mixture
openmc.stats.Normal
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.stats.muir
Angular Distributions

99
docs/source/releasenotes/0.13.2.rst Normal file
View file

@ -0,0 +1,99 @@
====================
What's New in 0.13.2
====================
.. currentmodule:: openmc
-------
Summary
-------
This release of OpenMC includes several bug fixes, performance improvements for
complex geometries and depletion simulations, and other general enhancements.
Notably, a capability has been added to compute the photon spectra from decay of
unstable nuclides. Alongside that, a new :data:`openmc.config` configuration
variable has been introduced that allows easier configuration of data sources.
Additionally, users can now perform cell or material rejection when sampling
external source distributions. Finally,
------------------------------------
Compatibility Notes and Deprecations
------------------------------------
- If you are building against libMesh for unstructured mesh tally support,
version 1.6 or higher is now required.
- The ``openmc.stats.Muir`` class has been replaced by a
:func:`openmc.stats.muir` function that returns an instance of
:class:`openmc.stats.Normal`.
------------
New Features
------------
- The :meth:`openmc.Material.get_nuclide_atom_densities` method now takes an
optional ``nuclide`` argument.
- Functions/methods in the :mod:`openmc.deplete` module now accept times in
Julian years (``'a'``).
- The :meth:`openmc.Universe.plot` method now allows a pre-existing axes object
to be passed in.
- Performance optimization for geometries with many complex regions.
- Performance optimization for depletion by avoiding deepcopies and caching
reaction rates.
- The :class:`openmc.RegularMesh` class now has a
:meth:`~openmc.RegularMesh.from_domain` classmethod.
- The :class:`openmc.CylindricalMesh` class now has a
:meth:`~openmc.CylindricalMesh.from_domain` classmethod.
- Improved method to condense diffusion coefficients from the :mod:`openmc.mgxs`
module.
- A new :data:`openmc.config` configuration variable has been introduced that
allows data sources to be specified at runtime or via environment variables.
- The :class:`openmc.EnergyFunctionFilter` class now supports multiple
interpolation schemes, not just linear-linear interpolation.
- The :class:`openmc.DAGMCUniverse` class now has ``material_names``,
``n_cells``, and ``n_surfaces`` attributes.
- A new :func:`openmc.data.decay_photon_energy` function has been added that
returns the energy spectrum of photons emitted from the decay of an unstable
nuclide.
- The :class:`openmc.Material` class also has a new
:attr:`~openmc.Material.decay_photon_energy` attribute that gives the decay
photon energy spectrum from the material based on its constituent nuclides.
- The :class:`openmc.deplete.StepResult` now has a
:meth:`~openmc.deplete.StepResult.get_material` method.
- The :class:`openmc.Source` class now takes a ``domains`` argument that
specifies a list of cells, materials, or universes that is used to reject
source sites (i.e., if the sampled sites are not within the specified domain,
they are rejected).
---------
Bug Fixes
---------
- `Delay call to Tally::set_strides <https://github.com/openmc-dev/openmc/pull/2183>`_
- `Fix reading reference direction from XML for angular distributions <https://github.com/openmc-dev/openmc/pull/2204>`_
- `Fix erroneous behavior in Material.add_components <https://github.com/openmc-dev/openmc/pull/2205>`_
- `Fix reading thermal elastic data from ACE <https://github.com/openmc-dev/openmc/pull/2226>`_
- `Fix reading source file with time attribute <https://github.com/openmc-dev/openmc/pull/2228>`_
- `Fix conversion of multiple thermal scattering data files from ACE <https://github.com/openmc-dev/openmc/pull/2232>`_
- `Fix reading values from wwinp file <https://github.com/openmc-dev/openmc/pull/2242>`_
- `Handle possibility of .ppm file in Universe.plot <https://github.com/openmc-dev/openmc/pull/2251>`_
- `Update volume calc types to mitigate overflow issues <https://github.com/openmc-dev/openmc/pull/2270>`_
------------
Contributors
------------
- `Lewis Gross <https://github.com/lewisgross1296>`_
- `Andrew Johnson <https://github.com/drewejohnson>`_
- `Miriam Kreher <https://github.com/mkreher13>`_
- `James Logan <https://github.com/jlogan03>`_
- `Jose Ignacio Marquez Damien <https://github.com/marquezj>`_
- `Josh May <https://github.com/joshmay1>`_
- `Patrick Myers <https://github.com/myerspat>`_
- `Adam Nelson <https://github.com/nelsonag>`_
- `April Novak <https://github.com/aprilnovak>`_
- `Ethan Peterson <https://github.com/eepeterson>`_
- `Gavin Ridley <https://github.com/gridley>`_
- `Paul Romano <https://github.com/paulromano>`_
- `Patrick Shriwise <https://github.com/pshriwise>`_
- `Jonathan Shimwell <https://github.com/Shimwell>`_
- `Olek Yardas <https://github.com/yardasol>`_

1
docs/source/releasenotes/index.rst
View file

@ -7,6 +7,7 @@ Release Notes
.. toctree::
:maxdepth: 1
0.13.2
0.13.1
0.13.0
0.12.2

4
docs/source/usersguide/scripts.rst
View file

@ -157,9 +157,9 @@ geometry.xml
added. Any 'surfaces' attributes/elements on a cell will be renamed 'region'.
materials.xml
Nuclide names will be changed from ACE aliases (e.g., Am-242m) to HDF5/GND
Nuclide names will be changed from ACE aliases (e.g., Am-242m) to HDF5/GNDS
names (e.g., Am242_m1). Thermal scattering table names will be changed from
ACE aliases (e.g., HH2O) to HDF5/GND names (e.g., c_H_in_H2O).
ACE aliases (e.g., HH2O) to HDF5/GNDS names (e.g., c_H_in_H2O).
----------------------
``openmc-update-mgxs``

1
docs/source/usersguide/settings.rst
View file

@ -474,7 +474,6 @@ selected::
Some features related to photon transport are not currently implemented,
including:
* Tallying photon energy deposition.
* Generating a photon source from a neutron calculation that can be used
for a later fixed source photon calculation.
* Photoneutron reactions.

5
docs/source/usersguide/tallies.rst
View file

@ -269,9 +269,8 @@ The following tables show all valid scores:
|heating |Total nuclear heating in units of eV per source |
| |particle. For neutrons, this corresponds to MT=301 |
| |produced by NJOY's HEATR module while for photons, |
| |this is tallied from either direct photon energy |
| |deposition (analog estimator) or pre-generated |
| |photon heating number. See :ref:`methods_heating` |
| |this is tallied from direct photon energy |
| |deposition. See :ref:`methods_heating`. |
+----------------------+---------------------------------------------------+
|heating-local |Total nuclear heating in units of eV per source |
| |particle assuming energy from secondary photons is |

146
include/openmc/cell.h
View file

@ -52,6 +52,100 @@ extern vector<unique_ptr<Cell>> cells;
} // namespace model
//==============================================================================
class Region {
public:
//----------------------------------------------------------------------------
// Constructors
Region() {}
explicit Region(std::string region_spec, int32_t cell_id);
//----------------------------------------------------------------------------
// Methods
//! \brief Determine if a cell contains the particle at a given location.
//!
//! The bounds of the cell are determined by a logical expression involving
//! surface half-spaces. The expression used is given in infix notation
//!
//! The function is split into two cases, one for simple cells (those
//! involving only the intersection of half-spaces) and one for complex cells.
//! Both cases use short circuiting; however, in the case fo complex cells,
//! the complexity increases with the binary operators involved.
//! \param r The 3D Cartesian coordinate to check.
//! \param u A direction used to "break ties" the coordinates are very
//! close to a surface.
//! \param on_surface The signed index of a surface that the coordinate is
//! known to be on. This index takes precedence over surface sense
//! calculations.
bool contains(Position r, Direction u, int32_t on_surface) const;
//! Find the oncoming boundary of this cell.
std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const;
//! Get the BoundingBox for this cell.
BoundingBox bounding_box(int32_t cell_id) const;
//! Get the CSG expression as a string
std::string str() const;
//! Get a vector containing all the surfaces in the region expression
vector<int32_t> surfaces() const;
//----------------------------------------------------------------------------
// Accessors
//! Get Boolean of if the cell is simple or not
bool is_simple() const { return simple_; }
private:
//----------------------------------------------------------------------------
// Private Methods
//! Get a vector of the region expression in postfix notation
vector<int32_t> generate_postfix(int32_t cell_id) const;
//! Determine if a particle is inside the cell for a simple cell (only
//! intersection operators)
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
//! Determine if a particle is inside the cell for a complex cell.
//!
//! Uses the comobination of half-spaces and binary operators to determine
//! if short circuiting can be used. Short cicuiting uses the relative and
//! absolute depth of parenthases in the expression.
bool contains_complex(Position r, Direction u, int32_t on_surface) const;
//! BoundingBox if the paritcle is in a simple cell.
BoundingBox bounding_box_simple() const;
//! BoundingBox if the particle is in a complex cell.
BoundingBox bounding_box_complex(vector<int32_t> postfix) const;
//! Enfource precedence: Parenthases, Complement, Intersection, Union
void add_precedence();
//! Add parenthesis to enforce precedence
std::vector<int32_t>::iterator add_parentheses(
std::vector<int32_t>::iterator start);
//! Remove complement operators from the expression
void remove_complement_ops();
//! Remove complement operators by using DeMorgan's laws
void apply_demorgan(
vector<int32_t>::iterator start, vector<int32_t>::iterator stop);
//----------------------------------------------------------------------------
// Private Data
//! Definition of spatial region as Boolean expression of half-spaces
// TODO: Should this be a vector of some other type
vector<int32_t> expression_;
bool simple_; //!< Does the region contain only intersections?
};
//==============================================================================
class Cell {
@ -108,6 +202,12 @@ public:
//! Get the BoundingBox for this cell.
virtual BoundingBox bounding_box() const = 0;
//! Get a vector of surfaces in the cell
virtual vector<int32_t> surfaces() const { return vector<int32_t>(); }
//! Check if the cell region expression is simple
virtual bool is_simple() const { return true; }
//----------------------------------------------------------------------------
// Accessors
@ -198,10 +298,6 @@ public:
//! T. The units are sqrt(eV).
vector<double> sqrtkT_;
//! Definition of spatial region as Boolean expression of half-spaces
vector<int32_t> region_;
bool simple_; //!< Does the region contain only intersections?
//! \brief Neighboring cells in the same universe.
NeighborList neighbors_;
@ -227,35 +323,36 @@ struct CellInstanceItem {
class CSGCell : public Cell {
public:
//----------------------------------------------------------------------------
// Constructors
CSGCell();
explicit CSGCell(pugi::xml_node cell_node);
bool contains(Position r, Direction u, int32_t on_surface) const override;
//----------------------------------------------------------------------------
// Methods
vector<int32_t> surfaces() const override { return region_.surfaces(); }
std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const override;
Position r, Direction u, int32_t on_surface, Particle* p) const override
{
return region_.distance(r, u, on_surface, p);
}
bool contains(Position r, Direction u, int32_t on_surface) const override
{
return region_.contains(r, u, on_surface);
}
BoundingBox bounding_box() const override
{
return region_.bounding_box(id_);
}
void to_hdf5_inner(hid_t group_id) const override;
BoundingBox bounding_box() const override;
bool is_simple() const override { return region_.is_simple(); }
protected:
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
bool contains_complex(Position r, Direction u, int32_t on_surface) const;
BoundingBox bounding_box_simple() const;
static BoundingBox bounding_box_complex(vector<int32_t> postfix);
//! Applies DeMorgan's laws to a section of the RPN
//! \param start Starting point for token modification
//! \param stop Stopping point for token modification
static void apply_demorgan(
vector<int32_t>::iterator start, vector<int32_t>::iterator stop);
//! Removes complement operators from the RPN
//! \param rpn The rpn to remove complement operators from.
static void remove_complement_ops(vector<int32_t>& rpn);
//! Returns the beginning position of a parenthesis block (immediately before
//! two surface tokens) in the RPN given a starting position at the end of
//! that block (immediately after two surface tokens)
@ -263,6 +360,9 @@ protected:
//! \param rpn The rpn being searched
static vector<int32_t>::iterator find_left_parenthesis(
vector<int32_t>::iterator start, const vector<int32_t>& rpn);
private:
Region region_;
};
//==============================================================================

1
include/openmc/constants.h
View file

@ -5,6 +5,7 @@
#define OPENMC_CONSTANTS_H
#include <cmath>
#include <cstdint>
#include <limits>
#include "openmc/array.h"

5
include/openmc/cross_sections.h
View file

@ -62,6 +62,11 @@ extern vector<Library> libraries;
//! libraries
void read_cross_sections_xml();
//! Read cross sections file (either XML or multigroup H5) and populate data
//! libraries
//! \param[in] root node of the cross_sections.xml
void read_cross_sections_xml(pugi::xml_node root);
//! Load nuclide and thermal scattering data from HDF5 files
//
//! \param[in] nuc_temps Temperatures for each nuclide in [K]

18
include/openmc/file_utils.h
View file

@ -3,14 +3,32 @@
#include <fstream> // for ifstream
#include <string>
#include <sys/stat.h>
namespace openmc {
// TODO: replace with std::filesystem when switch to C++17 is made
//! Determine if a path is a directory
//! \param[in] path Path to check
//! \return Whether the path is a directory
inline bool dir_exists(const std::string& path)
{
struct stat s;
if (stat(path.c_str(), &s) != 0)
return false;
return s.st_mode & S_IFDIR;
}
//! Determine if a file exists
//! \param[in] filename Path to file
//! \return Whether file exists
inline bool file_exists(const std::string& filename)
{
// rule out file being a path to a directory
if (dir_exists(filename))
return false;
std::ifstream s {filename};
return s.good();
}

6
include/openmc/geometry_aux.h
View file

@ -10,6 +10,7 @@
#include <vector>
#include "openmc/vector.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -19,8 +20,13 @@ extern std::unordered_map<int32_t, std::unordered_map<int32_t, int32_t>>
extern std::unordered_map<int32_t, int32_t> universe_level_counts;
} // namespace model
//! Read geometry from XML file
void read_geometry_xml();
//! Read geometry from XML node
//! \param[in] root node of geometry XML element
void read_geometry_xml(pugi::xml_node root);
//==============================================================================
//! Replace Universe, Lattice, and Material IDs with indices.
//==============================================================================

10
include/openmc/initialize.h
View file

@ -1,6 +1,8 @@
#ifndef OPENMC_INITIALIZE_H
#define OPENMC_INITIALIZE_H
#include <string>
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
@ -11,7 +13,13 @@ int parse_command_line(int argc, char* argv[]);
#ifdef OPENMC_MPI
void initialize_mpi(MPI_Comm intracomm);
#endif
void read_input_xml();
//! Read material, geometry, settings, and tallies from a single XML file
bool read_model_xml();
//! Read inputs from separate XML files
void read_separate_xml_files();
//! Write some output that occurs right after initialization
void initial_output();
} // namespace openmc

4
include/openmc/material.h
View file

@ -238,6 +238,10 @@ double density_effect(const vector<double>& f, const vector<double>& e_b_sq,
//! Read material data from materials.xml
void read_materials_xml();
//! Read material data XML node
//! \param[in] root node of materials XML element
void read_materials_xml(pugi::xml_node root);
void free_memory_material();
} // namespace openmc

36
include/openmc/mcpl_interface.h Normal file
View file

@ -0,0 +1,36 @@
#ifndef OPENMC_MCPL_INTERFACE_H
#define OPENMC_MCPL_INTERFACE_H
#include "openmc/particle_data.h"
#include "openmc/vector.h"
#include <string>
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
extern "C" const bool MCPL_ENABLED;
//==============================================================================
// Functions
//==============================================================================
//! Get a vector of source sites from an MCPL file
//
//! \param[in] path Path to MCPL file
//! \return Vector of source sites
vector<SourceSite> mcpl_source_sites(std::string path);
//! Write an MCPL source file
//
//! \param[in] filename Path to MCPL file
//! \param[in] surf_source_bank Whether to use the surface source bank
void write_mcpl_source_point(
const char* filename, bool surf_source_bank = false);
} // namespace openmc
#endif // OPENMC_MCPL_INTERFACE_H

3
include/openmc/output.h
View file

@ -40,6 +40,9 @@ void print_usage();
//! Display current version and copright/license information
void print_version();
//! Display compile flags employed, etc
void print_build_info();
//! Display header listing what physical values will displayed
void print_columns();

3
include/openmc/photon.h
View file

@ -100,6 +100,9 @@ public:
// Bremsstrahlung scaled DCS
xt::xtensor<double, 2> dcs_;
// Whether atomic relaxation data is present
bool has_atomic_relaxation_ {false};
// Constant data
static constexpr int MAX_STACK_SIZE =
7; //!< maximum possible size of atomic relaxation stack

4
include/openmc/plot.h
View file

@ -279,6 +279,10 @@ void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
//! Read plot specifications from a plots.xml file
void read_plots_xml();
//! Read plot specifications from an XML Node
//! \param[in] XML node containing plot info
void read_plots_xml(pugi::xml_node root);
//! Clear memory
void free_memory_plot();

7
include/openmc/settings.h
View file

@ -49,7 +49,9 @@ extern "C" bool run_CE; //!< run with continuous-energy data?
extern bool source_latest; //!< write latest source at each batch?
extern bool source_separate; //!< write source to separate file?
extern bool source_write; //!< write source in HDF5 files?
extern bool source_mcpl_write; //!< write source in mcpl files?
extern bool surf_source_write; //!< write surface source file?
extern bool surf_mcpl_write; //!< write surface mcpl file?
extern bool surf_source_read; //!< read surface source file?
extern bool survival_biasing; //!< use survival biasing?
extern bool temperature_multipole; //!< use multipole data?
@ -133,9 +135,12 @@ extern double
//==============================================================================
//! Read settings from XML file
//! \param[in] root XML node for <settings>
void read_settings_xml();
//! Read settings from XML node
//! \param[in] root XML node for <settings>
void read_settings_xml(pugi::xml_node root);
void free_memory_settings();
} // namespace openmc

1
include/openmc/source.h
View file

@ -102,6 +102,7 @@ class FileSource : public Source {
public:
// Constructors
explicit FileSource(std::string path);
explicit FileSource(const vector<SourceSite>& sites) : sites_ {sites} {}
// Methods
SourceSite sample(uint64_t* seed) const override;

33
include/openmc/tallies/filter.h
View file

@ -17,6 +17,34 @@
namespace openmc {
enum class FilterType {
AZIMUTHAL,
CELLBORN,
CELLFROM,
CELL,
CELL_INSTANCE,
COLLISION,
DELAYED_GROUP,
DISTRIBCELL,
ENERGY_FUNCTION,
ENERGY,
ENERGY_OUT,
LEGENDRE,
MATERIAL,
MESH,
MESH_SURFACE,
MU,
PARTICLE,
POLAR,
SPHERICAL_HARMONICS,
SPATIAL_LEGENDRE,
SURFACE,
TIME,
UNIVERSE,
ZERNIKE,
ZERNIKE_RADIAL
};
//==============================================================================
//! Modifies tally score events.
//==============================================================================
@ -58,7 +86,8 @@ public:
//----------------------------------------------------------------------------
// Methods
virtual std::string type() const = 0;
virtual std::string type_str() const = 0;
virtual FilterType type() const = 0;
//! Matches a tally event to a set of filter bins and weights.
//!
@ -72,7 +101,7 @@ public:
//! Writes data describing this filter to an HDF5 statepoint group.
virtual void to_statepoint(hid_t filter_group) const
{
write_dataset(filter_group, "type", type());
write_dataset(filter_group, "type", type_str());
write_dataset(filter_group, "n_bins", n_bins_);
}

3
include/openmc/tallies/filter_azimuthal.h
View file

@ -24,7 +24,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "azimuthal"; }
std::string type_str() const override { return "azimuthal"; }
FilterType type() const override { return FilterType::AZIMUTHAL; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_cell.h
View file

@ -25,7 +25,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "cell"; }
std::string type_str() const override { return "cell"; }
FilterType type() const override { return FilterType::CELL; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_cell_instance.h
View file

@ -28,7 +28,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "cellinstance"; }
std::string type_str() const override { return "cellinstance"; }
FilterType type() const override { return FilterType::CELL_INSTANCE; }
void from_xml(pugi::xml_node node) override;

5
include/openmc/tallies/filter_cellborn.h
View file

@ -11,12 +11,13 @@ namespace openmc {
//! Specifies which cell the particle was born in.
//==============================================================================
class CellbornFilter : public CellFilter {
class CellBornFilter : public CellFilter {
public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "cellborn"; }
std::string type_str() const override { return "cellborn"; }
FilterType type() const override { return FilterType::CELLBORN; }
void get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const override;

3
include/openmc/tallies/filter_cellfrom.h
View file

@ -16,7 +16,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "cellfrom"; }
std::string type_str() const override { return "cellfrom"; }
FilterType type() const override { return FilterType::CELLFROM; }
void get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const override;

3
include/openmc/tallies/filter_collision.h
View file

@ -23,7 +23,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "collision"; }
std::string type_str() const override { return "collision"; }
FilterType type() const override { return FilterType::COLLISION; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_delayedgroup.h
View file

@ -25,7 +25,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "delayedgroup"; }
std::string type_str() const override { return "delayedgroup"; }
FilterType type() const override { return FilterType::DELAYED_GROUP; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_distribcell.h
View file

@ -21,7 +21,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "distribcell"; }
std::string type_str() const override { return "distribcell"; }
FilterType type() const override { return FilterType::DISTRIBCELL; }
void from_xml(pugi::xml_node node) override;

6
include/openmc/tallies/filter_energy.h
View file

@ -22,7 +22,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "energy"; }
std::string type_str() const override { return "energy"; }
FilterType type() const override { return FilterType::ENERGY; }
void from_xml(pugi::xml_node node) override;
@ -63,7 +64,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "energyout"; }
std::string type_str() const override { return "energyout"; }
FilterType type() const override { return FilterType::ENERGY_OUT; }
void get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const override;

3
include/openmc/tallies/filter_energyfunc.h
View file

@ -24,7 +24,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "energyfunction"; }
std::string type_str() const override { return "energyfunction"; }
FilterType type() const override { return FilterType::ENERGY_FUNCTION; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_legendre.h
View file

@ -21,7 +21,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "legendre"; }
std::string type_str() const override { return "legendre"; }
FilterType type() const override { return FilterType::LEGENDRE; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_material.h
View file

@ -25,7 +25,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "material"; }
std::string type_str() const override { return "material"; }
FilterType type() const override { return FilterType::MATERIAL; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_mesh.h
View file

@ -24,7 +24,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "mesh"; }
std::string type_str() const override { return "mesh"; }
FilterType type() const override { return FilterType::MESH; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_meshsurface.h
View file

@ -10,7 +10,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "meshsurface"; }
std::string type_str() const override { return "meshsurface"; }
FilterType type() const override { return FilterType::MESH_SURFACE; }
void get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const override;

3
include/openmc/tallies/filter_mu.h
View file

@ -23,7 +23,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "mu"; }
std::string type_str() const override { return "mu"; }
FilterType type() const override { return FilterType::MU; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_particle.h
View file

@ -21,7 +21,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "particle"; }
std::string type_str() const override { return "particle"; }
FilterType type() const override { return FilterType::PARTICLE; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_polar.h
View file

@ -24,7 +24,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "polar"; }
std::string type_str() const override { return "polar"; }
FilterType type() const override { return FilterType::POLAR; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_sph_harm.h
View file

@ -25,7 +25,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "sphericalharmonics"; }
std::string type_str() const override { return "sphericalharmonics"; }
FilterType type() const override { return FilterType::SPHERICAL_HARMONICS; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_sptl_legendre.h
View file

@ -23,7 +23,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "spatiallegendre"; }
std::string type_str() const override { return "spatiallegendre"; }
FilterType type() const override { return FilterType::SPATIAL_LEGENDRE; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_surface.h
View file

@ -25,7 +25,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "surface"; }
std::string type_str() const override { return "surface"; }
FilterType type() const override { return FilterType::SURFACE; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_time.h
View file

@ -22,7 +22,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "time"; }
std::string type_str() const override { return "time"; }
FilterType type() const override { return FilterType::TIME; }
void from_xml(pugi::xml_node node) override;

3
include/openmc/tallies/filter_universe.h
View file

@ -25,7 +25,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "universe"; }
std::string type_str() const override { return "universe"; }
FilterType type() const override { return FilterType::UNIVERSE; }
void from_xml(pugi::xml_node node) override;

6
include/openmc/tallies/filter_zernike.h
View file

@ -21,7 +21,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "zernike"; }
std::string type_str() const override { return "zernike"; }
FilterType type() const override { return FilterType::ZERNIKE; }
void from_xml(pugi::xml_node node) override;
@ -72,7 +73,8 @@ public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override { return "zernikeradial"; }
std::string type_str() const override { return "zernikeradial"; }
FilterType type() const override { return FilterType::ZERNIKE_RADIAL; }
void get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const override;

8
include/openmc/tallies/tally.h
View file

@ -48,10 +48,10 @@ public:
void set_nuclides(const vector<std::string>& nuclides);
//! returns vector of indices corresponding to the tally this is called on
const vector<int32_t>& filters() const { return filters_; }
const vector<int32_t>& filters() const { return filters_; }
//! \brief Returns the tally filter at index i
int32_t filters(int i) const { return filters_[i]; }
int32_t filters(int i) const { return filters_[i]; }
void set_filters(gsl::span<Filter*> filters);
@ -178,6 +178,10 @@ extern double global_tally_leakage;
//! Read tally specification from tallies.xml
void read_tallies_xml();
//! Read tally specification from an XML node
//! \param[in] root node of tallies XML element
void read_tallies_xml(pugi::xml_node root);
//! \brief Accumulate the sum of the contributions from each history within the
//! batch to a new random variable
void accumulate_tallies();

7
include/openmc/volume_calc.h
View file

@ -1,6 +1,9 @@
#ifndef OPENMC_VOLUME_CALC_H
#define OPENMC_VOLUME_CALC_H
#include <cstdint>
#include <string>
#include "openmc/array.h"
#include "openmc/position.h"
#include "openmc/tallies/trigger.h"
@ -10,7 +13,6 @@
#include "xtensor/xtensor.hpp"
#include <gsl/gsl-lite.hpp>
#include <string>
namespace openmc {
@ -69,7 +71,8 @@ private:
//! \param[in] i_material Index in global materials vector
//! \param[in,out] indices Vector of material indices
//! \param[in,out] hits Number of hits corresponding to each material
void check_hit(int i_material, vector<int>& indices, vector<int>& hits) const;
void check_hit(
int i_material, vector<uint64_t>& indices, vector<uint64_t>& hits) const;
};
//==============================================================================

6
openmc/__init__.py
View file

@ -10,11 +10,11 @@ from openmc.material import *
from openmc.plots import *
from openmc.region import *
from openmc.volume import *
from openmc.source import *
from openmc.weight_windows import *
from openmc.settings import *
from openmc.surface import *
from openmc.universe import *
from openmc.source import *
from openmc.settings import *
from openmc.lattice import *
from openmc.filter import *
from openmc.filter_expansion import *
@ -38,4 +38,4 @@ from .config import *
from openmc.model import rectangular_prism, hexagonal_prism, Model
__version__ = '0.14.0-dev'
__version__ = '0.13.3-dev'

32
openmc/_xml.py
View file

@ -1,22 +1,40 @@
def clean_indentation(element, level=0, spaces_per_level=2):
"""
copy and paste from https://effbot.org/zone/element-lib.htm#prettyprint
it basically walks your tree and adds spaces and newlines so the tree is
printed in a nice way
def clean_indentation(element, level=0, spaces_per_level=2, trailing_indent=True):
"""Set indentation of XML element and its sub-elements.
Copied and pasted from https://effbot.org/zone/element-lib.htm#prettyprint.
It walks your tree and adds spaces and newlines so the tree is
printed in a nice way.
Parameters
----------
level : int
Indentation level for the element passed in (default 0)
spaces_per_level : int
Number of spaces per indentation level (default 2)
trailing_indent : bool
Whether or not to add indentation after closing the element
"""
i = "\n" + level*spaces_per_level*" "
# ensure there's always some tail for the element passed in
if not element.tail:
element.tail = ""
if len(element):
if not element.text or not element.text.strip():
element.text = i + spaces_per_level*" "
if not element.tail or not element.tail.strip():
if trailing_indent and (not element.tail or not element.tail.strip()):
element.tail = i
for sub_element in element:
# `trailing_indent` is intentionally not forwarded to the recursive
# call. Any child element of the topmost element should add
# indentation at the end to ensure its parent's indentation is
# correct.
clean_indentation(sub_element, level+1, spaces_per_level)
if not sub_element.tail or not sub_element.tail.strip():
sub_element.tail = i
else:
if level and (not element.tail or not element.tail.strip()):
if trailing_indent and level and (not element.tail or not element.tail.strip()):
element.tail = i

6
openmc/cell.py
View file

@ -318,12 +318,12 @@ class Cell(IDManagerMixin):
@temperature.setter
def temperature(self, temperature):
# Make sure temperatures are positive
cv.check_type('cell temperature', temperature, (Iterable, Real))
cv.check_type('cell temperature', temperature, (Iterable, Real), none_ok=True)
if isinstance(temperature, Iterable):
cv.check_type('cell temperature', temperature, Iterable, Real)
for T in temperature:
cv.check_greater_than('cell temperature', T, 0.0, True)
else:
elif isinstance(temperature, Real):
cv.check_greater_than('cell temperature', temperature, 0.0, True)
# If this cell is filled with a universe or lattice, propagate
@ -650,7 +650,7 @@ class Cell(IDManagerMixin):
surfaces : dict
Dictionary mapping surface IDs to :class:`openmc.Surface` instances
materials : dict
Dictionary mapping material IDs to :class:`openmc.Material`
Dictionary mapping material ID strings to :class:`openmc.Material`
instances (defined in :math:`openmc.Geometry.from_xml`)
get_universe : function
Function returning universe (defined in

4
openmc/checkvalue.py
View file

@ -1,12 +1,12 @@
import copy
import os
from typing import Union
import typing # required to prevent typing.Union namespace overwriting Union
from collections.abc import Iterable
import numpy as np
# Type for arguments that accept file paths
PathLike = Union[str, os.PathLike]
PathLike = typing.Union[str, os.PathLike]
def check_type(name, value, expected_type, expected_iter_type=None, *, none_ok=False):

4
openmc/config.py
View file

@ -42,7 +42,9 @@ class _Config(MutableMapping):
# Reset photon source data since it relies on chain file
_DECAY_PHOTON_ENERGY.clear()
else:
raise KeyError(f'Unrecognized config key: {key}')
raise KeyError(f'Unrecognized config key: {key}. Acceptable keys '
'are "cross_sections", "mg_cross_sections" and '
'"chain_file"')
def __iter__(self):
return iter(self._mapping)

4
openmc/data/ace.py
View file

@ -24,7 +24,7 @@ import numpy as np
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from .data import ATOMIC_SYMBOL, gnd_name, EV_PER_MEV, K_BOLTZMANN
from .data import ATOMIC_SYMBOL, gnds_name, EV_PER_MEV, K_BOLTZMANN
from .endf import ENDF_FLOAT_RE
@ -88,7 +88,7 @@ def get_metadata(zaid, metastable_scheme='nndc'):
# Determine name
element = ATOMIC_SYMBOL[Z]
name = gnd_name(Z, mass_number, metastable)
name = gnds_name(Z, mass_number, metastable)
return (name, element, Z, mass_number, metastable)

33
openmc/data/data.py
View file

@ -197,8 +197,8 @@ NEUTRON_MASS = 1.00866491595
# Used in atomic_mass function as a cache
_ATOMIC_MASS = {}
# Regex for GND nuclide names (used in zam function)
_GND_NAME_RE = re.compile(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)')
# Regex for GNDS nuclide names (used in zam function)
_GNDS_NAME_RE = re.compile(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)')
# Used in half_life function as a cache
_HALF_LIFE = {}
@ -223,19 +223,19 @@ def atomic_mass(isotope):
"""
if not _ATOMIC_MASS:
# Load data from AME2016 file
mass_file = os.path.join(os.path.dirname(__file__), 'mass16.txt')
# Load data from AME2020 file
mass_file = os.path.join(os.path.dirname(__file__), 'mass_1.mas20.txt')
with open(mass_file, 'r') as ame:
# Read lines in file starting at line 40
for line in itertools.islice(ame, 39, None):
# Read lines in file starting at line 37
for line in itertools.islice(ame, 36, None):
name = f'{line[20:22].strip()}{int(line[16:19])}'
mass = float(line[96:99]) + 1e-6*float(
line[100:106] + '.' + line[107:112])
mass = float(line[106:109]) + 1e-6*float(
line[110:116] + '.' + line[117:123])
_ATOMIC_MASS[name.lower()] = mass
# For isotopes found in some libraries that represent all natural
# isotopes of their element (e.g. C0), calculate the atomic mass as
# the sum of the atomic mass times the natural abudance of the isotopes
# the sum of the atomic mass times the natural abundance of the isotopes
# that make up the element.
for element in ['C', 'Zn', 'Pt', 'Os', 'Tl']:
isotope_zero = element.lower() + '0'
@ -436,8 +436,11 @@ def water_density(temperature, pressure=0.1013):
return coeff / pi / gamma1_pi
def gnd_name(Z, A, m=0):
"""Return nuclide name using GND convention
def gnds_name(Z, A, m=0):
"""Return nuclide name using GNDS convention
.. versionchanged:: 0.14.0
Function name changed from ``gnd_name`` to ``gnds_name``
Parameters
----------
@ -451,7 +454,7 @@ def gnd_name(Z, A, m=0):
Returns
-------
str
Nuclide name in GND convention, e.g., 'Am242_m1'
Nuclide name in GNDS convention, e.g., 'Am242_m1'
"""
if m > 0:
@ -502,7 +505,7 @@ def zam(name):
Parameters
----------
name : str
Name of nuclide using GND convention, e.g., 'Am242_m1'
Name of nuclide using GNDS convention, e.g., 'Am242_m1'
Returns
-------
@ -511,10 +514,10 @@ def zam(name):
"""
try:
symbol, A, state = _GND_NAME_RE.match(name).groups()
symbol, A, state = _GNDS_NAME_RE.match(name).groups()
except AttributeError:
raise ValueError(f"'{name}' does not appear to be a nuclide name in "
"GND format")
"GNDS format")
if symbol not in ATOMIC_NUMBER:
raise ValueError(f"'{symbol}' is not a recognized element symbol")

50
openmc/data/decay.py
View file

@ -144,7 +144,7 @@ class FissionProductYields(EqualityMixin):
# Assign basic nuclide properties
self.nuclide = {
'name': ev.gnd_name,
'name': ev.gnds_name,
'atomic_number': ev.target['atomic_number'],
'mass_number': ev.target['mass_number'],
'isomeric_state': ev.target['isomeric_state']
@ -587,7 +587,7 @@ def decay_photon_energy(nuclide: str) -> Optional[Univariate]:
for the first time, you need to ensure that a depletion chain has been
specified in openmc.config['chain_file'].
.. versionadded:: 0.14.0
.. versionadded:: 0.13.2
Parameters
----------
@ -621,3 +621,49 @@ def decay_photon_energy(nuclide: str) -> Optional[Univariate]:
"sources listed.")
return _DECAY_PHOTON_ENERGY.get(nuclide)
_DECAY_ENERGY = {}
def decay_energy(nuclide: str):
"""Get decay energy value resulting from the decay of a nuclide
This function relies on data stored in a depletion chain. Before calling it
for the first time, you need to ensure that a depletion chain has been
specified in openmc.config['chain_file'].
.. versionadded:: 0.13.3
Parameters
----------
nuclide : str
Name of nuclide, e.g., 'H3'
Returns
-------
float
Decay energy of nuclide in [eV]. If the nuclide is stable, a value of
0.0 is returned.
"""
if not _DECAY_ENERGY:
chain_file = openmc.config.get('chain_file')
if chain_file is None:
raise DataError(
"A depletion chain file must be specified with "
"openmc.config['chain_file'] in order to load decay data."
)
from openmc.deplete import Chain
chain = Chain.from_xml(chain_file)
for nuc in chain.nuclides:
if nuc.decay_energy:
_DECAY_ENERGY[nuc.name] = nuc.decay_energy
# If the chain file contained no decay energy, warn the user
if not _DECAY_ENERGY:
warn(f"Chain file '{chain_file}' does not have any decay energy.")
return _DECAY_ENERGY.get(nuclide, 0.0)

12
openmc/data/endf.py
View file

@ -12,7 +12,7 @@ import re
import numpy as np
from .data import gnd_name
from .data import gnds_name
from .function import Tabulated1D
try:
from ._endf import float_endf
@ -520,10 +520,10 @@ class Evaluation:
self.reaction_list.append((mf, mt, nc, mod))
@property
def gnd_name(self):
return gnd_name(self.target['atomic_number'],
self.target['mass_number'],
self.target['isomeric_state'])
def gnds_name(self):
return gnds_name(self.target['atomic_number'],
self.target['mass_number'],
self.target['isomeric_state'])
class Tabulated2D:
@ -531,7 +531,7 @@ class Tabulated2D:
This is a dummy class that is not really used other than to store the
interpolation information for a two-dimensional function. Once we refactor
to adopt GND-like data containers, this will probably be removed or
to adopt GNDS-like data containers, this will probably be removed or
extended.
Parameters

3475
openmc/data/mass16.txt
View file

File diff suppressed because it is too large Load diff

3594
openmc/data/mass_1.mas20.txt Normal file
View file

File diff suppressed because it is too large Load diff

4
openmc/data/multipole.py
View file

@ -748,12 +748,12 @@ class WindowedMultipole(EqualityMixin):
Parameters
----------
name : str
Name of the nuclide using the GND naming convention
Name of the nuclide using the GNDS naming convention
Attributes
----------
name : str
Name of the nuclide using the GND naming convention
Name of the nuclide using the GNDS naming convention
spacing : float
The width of each window in sqrt(E)-space. For example, the frst window
will end at (sqrt(E_min) + spacing)**2 and the second window at

7
openmc/data/neutron.py
View file

@ -44,7 +44,7 @@ class IncidentNeutron(EqualityMixin):
Parameters
----------
name : str
Name of the nuclide using the GND naming convention
Name of the nuclide using the GNDS naming convention
atomic_number : int
Number of protons in the target nucleus
mass_number : int
@ -75,7 +75,7 @@ class IncidentNeutron(EqualityMixin):
Metastable state of the target nucleus. A value of zero indicates ground
state.
name : str
Name of the nuclide using the GND naming convention
Name of the nuclide using the GNDS naming convention
reactions : collections.OrderedDict
Contains the cross sections, secondary angle and energy distributions,
and other associated data for each reaction. The keys are the MT values
@ -868,9 +868,8 @@ class IncidentNeutron(EqualityMixin):
heatr_evals = get_evaluations(kwargs["heatr"])
heatr_local_evals = get_evaluations(kwargs["heatr"] + "_local")
for ev, ev_local in zip(heatr_evals, heatr_local_evals):
temp = "{}K".format(round(ev.target["temperature"]))
for ev, ev_local, temp in zip(heatr_evals, heatr_local_evals, data.temperatures):
# Get total KERMA (originally from ACE file) and energy grid
kerma = data.reactions[301].xs[temp]
E = kerma.x

4
openmc/data/photon.py
View file

@ -13,7 +13,7 @@ from scipy.interpolate import CubicSpline
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from . import HDF5_VERSION
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Table, get_metadata, get_table
from .data import ATOMIC_SYMBOL, EV_PER_MEV
from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record
@ -143,6 +143,8 @@ class AtomicRelaxation(EqualityMixin):
Dictionary indicating the number of electrons in a subshell when neutral
(values) for given subshells (keys). The subshells should be given as
strings, e.g., 'K', 'L1', 'L2', etc.
subshells : list
List of subshells as strings, e.g. ``['K', 'L1', ...]``
transitions : pandas.DataFrame
Dictionary indicating allowed transitions and their probabilities
(values) for given subshells (keys). The subshells should be given as

2
openmc/data/reaction.py
View file

@ -555,7 +555,7 @@ def _get_activation_products(ev, rx):
Z, A = divmod(items[2], 1000)
excited_state = items[3]
# Get GND name for product
# Get GNDS name for product
symbol = ATOMIC_SYMBOL[Z]
if excited_state > 0:
name = '{}{}_e{}'.format(symbol, A, excited_state)

39
openmc/data/thermal.py
View file

@ -104,7 +104,7 @@ def get_thermal_name(name):
Returns
-------
str
GND-format thermal scattering name
GNDS-format thermal scattering name
"""
if name in _THERMAL_NAMES:
@ -396,7 +396,7 @@ class ThermalScattering(EqualityMixin):
Parameters
----------
name : str
Name of the material using GND convention, e.g. c_H_in_H2O
Name of the material using GNDS convention, e.g. c_H_in_H2O
atomic_weight_ratio : float
Atomic mass ratio of the target nuclide.
kTs : Iterable of float
@ -415,7 +415,7 @@ class ThermalScattering(EqualityMixin):
Inelastic scattering cross section derived in the incoherent
approximation
name : str
Name of the material using GND convention, e.g. c_H_in_H2O
Name of the material using GNDS convention, e.g. c_H_in_H2O
temperatures : Iterable of str
List of string representations the temperatures of the target nuclide
in the data set. The temperatures are strings of the temperature,
@ -491,7 +491,7 @@ class ThermalScattering(EqualityMixin):
ACE table to read from. If given as a string, it is assumed to be
the filename for the ACE file.
name : str
GND-conforming name of the material, e.g. c_H_in_H2O. If none is
GNDS-conforming name of the material, e.g. c_H_in_H2O. If none is
passed, the appropriate name is guessed based on the name of the ACE
table.
@ -596,7 +596,7 @@ class ThermalScattering(EqualityMixin):
ACE table to read from. If given as a string, it is assumed to be
the filename for the ACE file.
name : str
GND-conforming name of the material, e.g. c_H_in_H2O. If none is
GNDS-conforming name of the material, e.g. c_H_in_H2O. If none is
passed, the appropriate name is guessed based on the name of the ACE
table.
@ -670,12 +670,41 @@ class ThermalScattering(EqualityMixin):
mu_i = []
for j in range(n_energy_out[i]):
mu = ace.xss[idx + 4:idx + 4 + n_mu]
# The equiprobable angles produced by NJOY are not always
# sorted. This is problematic when the smearing algorithm
# is applied when sampling the angles. We sort the angles
# here, because they are equiprobable, so the order
# doesn't matter.
mu.sort()
p_mu = 1. / n_mu * np.ones(n_mu)
mu_ij = Discrete(mu, p_mu)
mu_ij.c = np.cumsum(p_mu)
mu_i.append(mu_ij)
idx += 3 + n_mu
# Check if the CDF for the outgoing energy distribution starts
# at 0. For NJOY and FRENDY evaluations, this is never the case,
# and can very rarely lead to negative energies when sampling
# the outgoing energy. From Eq. 7.6 of the ENDF manual, we can
# add an outgoing energy 0 eV that has a PDF of 0 (and of
# course, a CDF of 0 as well).
if eout_i.c[0] > 0.:
eout_i.x = np.insert(eout_i.x, 0, 0.)
eout_i.p = np.insert(eout_i.p, 0, 0.)
eout_i.c = np.insert(eout_i.c, 0, 0.)
# For this added outgoing energy (of 0 eV) we add a set of
# isotropic discrete angles.
dmu = 2. / n_mu
mu = np.linspace(-1. + 0.5*dmu, 1. - 0.5*dmu, n_mu)
p_mu = 1. / n_mu * np.ones(n_mu)
mu_0 = Discrete(mu, p_mu)
mu_0.c = np.cumsum(p_mu)
mu_i.insert(0, mu_0)
# We don't worry about renormalizing the outgoing energy PDF/CDF
# after this manipulation, because it never seems to be
# normalized to begin with (at least with NJOY).
energy_out.append(eout_i)
mu_out.append(mu_i)

10
openmc/deplete/abc.py
View file

@ -265,8 +265,8 @@ class ReactionRateHelper(ABC):
Ordering of reactions
"""
def divide_by_adens(self, number):
"""Normalize reaction rates by number of nuclides
def divide_by_atoms(self, number):
"""Normalize reaction rates by number of atoms
Acts on the current material examined by :meth:`get_material_rates`
@ -798,7 +798,7 @@ class Integrator(ABC):
t, self._i_res = self._get_start_data()
for i, (dt, source_rate) in enumerate(self):
if output:
if output and comm.rank == 0:
print(f"[openmc.deplete] t={t} s, dt={dt} s, source={source_rate}")
# Solve transport equation (or obtain result from restart)
@ -818,7 +818,7 @@ class Integrator(ABC):
conc = conc_list.pop()
StepResult.save(self.operator, conc_list, res_list, [t, t + dt],
source_rate, self._i_res + i, proc_time)
source_rate, self._i_res + i, proc_time)
t += dt
@ -826,7 +826,7 @@ class Integrator(ABC):
# source rate is passed to the transport operator (which knows to
# just return zero reaction rates without actually doing a transport
# solve)
if output and final_step:
if output and final_step and comm.rank == 0:
print(f"[openmc.deplete] t={t} (final operator evaluation)")
res_list = [self.operator(conc, source_rate if final_step else 0.0)]
StepResult.save(self.operator, [conc], res_list, [t, t],

39
openmc/deplete/chain.py
View file

@ -15,7 +15,7 @@ from numbers import Real, Integral
from warnings import warn
from openmc.checkvalue import check_type, check_greater_than
from openmc.data import gnd_name, zam, DataLibrary
from openmc.data import gnds_name, zam, DataLibrary
from openmc.exceptions import DataError
from .nuclide import FissionYieldDistribution
@ -135,14 +135,14 @@ def replace_missing(product, decay_data):
Parameters
----------
product : str
Name of product in GND format, e.g. 'Y86_m1'.
Name of product in GNDS format, e.g. 'Y86_m1'.
decay_data : dict
Dictionary of decay data
Returns
-------
product : str
Replacement for missing product in GND format.
Replacement for missing product in GNDS format.
"""
# Determine atomic number, mass number, and metastable state
@ -213,7 +213,7 @@ def replace_missing_fpy(actinide, fpy_data, decay_data):
# Check if metastable state has data (e.g., Am242m)
Z, A, m = zam(actinide)
if m == 0:
metastable = gnd_name(Z, A, 1)
metastable = gnds_name(Z, A, 1)
if metastable in fpy_data:
return metastable
@ -222,7 +222,7 @@ def replace_missing_fpy(actinide, fpy_data, decay_data):
while isotone in decay_data:
Z += 1
A += 1
isotone = gnd_name(Z, A, 0)
isotone = gnds_name(Z, A, 0)
if isotone in fpy_data:
return isotone
@ -231,7 +231,7 @@ def replace_missing_fpy(actinide, fpy_data, decay_data):
while isotone in decay_data:
Z -= 1
A -= 1
isotone = gnd_name(Z, A, 0)
isotone = gnds_name(Z, A, 0)
if isotone in fpy_data:
return isotone
@ -357,7 +357,7 @@ class Chain:
reactions = {}
for f in neutron_files:
evaluation = openmc.data.endf.Evaluation(f)
name = evaluation.gnd_name
name = evaluation.gnds_name
reactions[name] = {}
for mf, mt, nc, mod in evaluation.reaction_list:
if mf == 3:
@ -630,15 +630,27 @@ class Chain:
# Gain from radioactive decay
if nuc.n_decay_modes != 0:
for _, target, branching_ratio in nuc.decay_modes:
# Allow for total annihilation for debug purposes
if target is not None:
branch_val = branching_ratio * decay_constant
for decay_type, target, branching_ratio in nuc.decay_modes:
branch_val = branching_ratio * decay_constant
if branch_val != 0.0:
# Allow for total annihilation for debug purposes
if branch_val != 0.0:
if target is not None:
k = self.nuclide_dict[target]
matrix[k, i] += branch_val
# Produce alphas and protons from decay
if 'alpha' in decay_type:
k = self.nuclide_dict.get('He4')
if k is not None:
count = decay_type.count('alpha')
matrix[k, i] += count * branch_val
elif 'p' in decay_type:
k = self.nuclide_dict.get('H1')
if k is not None:
count = decay_type.count('p')
matrix[k, i] += count * branch_val
if nuc.name in rates.index_nuc:
# Extract all reactions for this nuclide in this cell
nuc_ind = rates.index_nuc[nuc.name]
@ -892,7 +904,7 @@ class Chain:
ground_target = grounds.get(parent_name)
if ground_target is None:
pz, pa, pm = zam(parent_name)
ground_target = gnd_name(pz, pa + 1, 0)
ground_target = gnds_name(pz, pa + 1, 0)
new_ratios[ground_target] = ground_br
parent.add_reaction(reaction, ground_target, rxn_Q, ground_br)
@ -1025,6 +1037,7 @@ class Chain:
new_nuclide = Nuclide(previous.name)
new_nuclide.half_life = previous.half_life
new_nuclide.decay_energy = previous.decay_energy
new_nuclide.sources = previous.sources.copy()
if hasattr(previous, '_fpy'):
new_nuclide._fpy = previous._fpy

3
openmc/deplete/coupled_operator.py
View file

@ -460,7 +460,6 @@ class CoupledOperator(OpenMCOperator):
# Run OpenMC
openmc.lib.run()
openmc.lib.reset_timers()
# Extract results
rates = self._calculate_reaction_rates(source_rate)
@ -529,6 +528,8 @@ class CoupledOperator(OpenMCOperator):
"openmc_simulation_n{}.h5".format(step),
write_source=False)
openmc.lib.reset_timers()
def finalize(self):
"""Finalize a depletion simulation and release resources."""
if self.cleanup_when_done:

22
openmc/deplete/independent_operator.py
View file

@ -31,6 +31,9 @@ class IndependentOperator(OpenMCOperator):
passed to an integrator class, such as
:class:`openmc.deplete.CECMIntegrator`.
Note that passing an empty :class:`~openmc.deplete.MicroXS` instance to the
``micro_xs`` argument allows a decay-only calculation to be run.
.. versionadded:: 0.13.1
Parameters
@ -38,9 +41,11 @@ class IndependentOperator(OpenMCOperator):
materials : openmc.Materials
Materials to deplete.
micro_xs : MicroXS
One-group microscopic cross sections in [b] .
One-group microscopic cross sections in [b]. If the
:class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation will
be run.
chain_file : str
Path to the depletion chain XML file. Defaults to
Path to the depletion chain XML file. Defaults to
``openmc.config['chain_file']``.
keff : 2-tuple of float, optional
keff eigenvalue and uncertainty from transport calculation.
@ -150,7 +155,7 @@ class IndependentOperator(OpenMCOperator):
@classmethod
def from_nuclides(cls, volume, nuclides,
micro_xs,
chain_file,
chain_file=None,
nuc_units='atom/b-cm',
keff=None,
normalization_mode='fission-q',
@ -169,10 +174,13 @@ class IndependentOperator(OpenMCOperator):
Dictionary with nuclide names as keys and nuclide concentrations as
values.
micro_xs : MicroXS
One-group microscopic cross sections.
chain_file : str
Path to the depletion chain XML file.
nuc_units : {'atom/cm3', 'atom/b-cm'}
One-group microscopic cross sections in [b]. If the
:class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation
will be run.
chain_file : str, optional
Path to the depletion chain XML file. Defaults to
``openmc.config['chain_file']``.
nuc_units : {'atom/cm3', 'atom/b-cm'}, optional
Units for nuclide concentration.
keff : 2-tuple of float, optional
keff eigenvalue and uncertainty from transport calculation.

2
openmc/deplete/nuclide.py
View file

@ -82,7 +82,7 @@ class Nuclide:
Parameters
----------
name : str, optional
GND name of this nuclide, e.g. ``"He4"``, ``"Am242_m1"``
GNDS name of this nuclide, e.g. ``"He4"``, ``"Am242_m1"``
Attributes
----------

6
openmc/deplete/openmc_operator.py
View file

@ -541,10 +541,10 @@ class OpenMCOperator(TransportOperator):
self._normalization_helper.update(
tally_rates[:, fission_ind])
# Divide by total number and store
rates[i] = self._rate_helper.divide_by_adens(number)
# Divide by total number of atoms and store
rates[i] = self._rate_helper.divide_by_atoms(number)
# Scale reaction rates to obtain units of reactions/sec
# Scale reaction rates to obtain units of (reactions/sec)/atom
rates *= self._normalization_helper.factor(source_rate)
# Store new fission yields on the chain

7
openmc/deplete/results.py
View file

@ -1,7 +1,8 @@
import numbers
import bisect
import math
from typing import Iterable, Optional, Tuple, Union
import typing # required to prevent typing.Union namespace overwriting Union
from typing import Iterable, Optional, Tuple
from warnings import warn
import h5py
@ -95,7 +96,7 @@ class Results(list):
def get_atoms(
self,
mat: Union[Material, str],
mat: typing.Union[Material, str],
nuc: str,
nuc_units: str = "atoms",
time_units: str = "s"
@ -165,7 +166,7 @@ class Results(list):
def get_reaction_rate(
self,
mat: Union[Material, str],
mat: typing.Union[Material, str],
nuc: str,
rx: str
) -> Tuple[np.ndarray, np.ndarray]:

28
openmc/deplete/stepresult.py
View file

@ -202,7 +202,7 @@ class StepResult:
def get_material(self, mat_id):
"""Return material object for given depleted composition
.. versionadded:: 0.14.0
.. versionadded:: 0.13.2
Parameters
----------
@ -318,10 +318,11 @@ class StepResult:
chunks=(1, 1, n_mats, n_nuc_number),
dtype='float64')
handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn),
maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn),
chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn),
dtype='float64')
if n_nuc_rxn > 0 and n_rxn > 0:
handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn),
maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn),
chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn),
dtype='float64')
handle.create_dataset("eigenvalues", (1, n_stages, 2),
maxshape=(None, n_stages, 2), dtype='float64')
@ -358,7 +359,9 @@ class StepResult:
# Grab handles
number_dset = handle["/number"]
rxn_dset = handle["/reaction rates"]
has_reactions = ("reaction rates" in handle)
if has_reactions:
rxn_dset = handle["/reaction rates"]
eigenvalues_dset = handle["/eigenvalues"]
time_dset = handle["/time"]
source_rate_dset = handle["/source_rate"]
@ -375,9 +378,10 @@ class StepResult:
number_shape[0] = new_shape
number_dset.resize(number_shape)
rxn_shape = list(rxn_dset.shape)
rxn_shape[0] = new_shape
rxn_dset.resize(rxn_shape)
if has_reactions:
rxn_shape = list(rxn_dset.shape)
rxn_shape[0] = new_shape
rxn_dset.resize(rxn_shape)
eigenvalues_shape = list(eigenvalues_dset.shape)
eigenvalues_shape[0] = new_shape
@ -407,7 +411,8 @@ class StepResult:
high = max(inds)
for i in range(n_stages):
number_dset[index, i, low:high+1] = self.data[i]
rxn_dset[index, i, low:high+1] = self.rates[i]
if has_reactions:
rxn_dset[index, i, low:high+1] = self.rates[i]
if comm.rank == 0:
eigenvalues_dset[index, i] = self.k[i]
if comm.rank == 0:
@ -483,7 +488,8 @@ class StepResult:
for i in range(results.n_stages):
rate = ReactionRates(results.index_mat, rxn_nuc_to_ind, rxn_to_ind, True)
rate[:] = handle["/reaction rates"][step, i, :, :, :]
if "reaction rates" in handle:
rate[:] = handle["/reaction rates"][step, i, :, :, :]
results.rates.append(rate)
return results

59
openmc/executor.py
View file

@ -9,7 +9,7 @@ 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):
openmc_exec='openmc', mpi_args=None, path_input=None):
"""Converts user-readable flags in to command-line arguments to be run with
the OpenMC executable via subprocess.
@ -42,6 +42,9 @@ def _process_CLI_arguments(volume=False, geometry_debug=False, particles=None,
mpi_args : list of str, optional
MPI execute command and any additional MPI arguments to pass,
e.g. ['mpiexec', '-n', '8'].
path_input : str or Pathlike
Path to a single XML file or a directory containing XML files for the
OpenMC executable to read.
.. versionadded:: 0.13.0
@ -82,6 +85,9 @@ def _process_CLI_arguments(volume=False, geometry_debug=False, particles=None,
if mpi_args is not None:
args = mpi_args + args
if path_input is not None:
args += [path_input]
return args
@ -118,7 +124,7 @@ def _run(args, output, cwd):
raise RuntimeError(error_msg)
def plot_geometry(output=True, openmc_exec='openmc', cwd='.'):
def plot_geometry(output=True, openmc_exec='openmc', cwd='.', path_input=None):
"""Run OpenMC in plotting mode
Parameters
@ -129,6 +135,9 @@ def plot_geometry(output=True, openmc_exec='openmc', cwd='.'):
Path to OpenMC executable
cwd : str, optional
Path to working directory to run in
path_input : str
Path to a single XML file or a directory containing XML files for the
OpenMC executable to read.
Raises
------
@ -136,10 +145,13 @@ def plot_geometry(output=True, openmc_exec='openmc', cwd='.'):
If the `openmc` executable returns a non-zero status
"""
_run([openmc_exec, '-p'], output, cwd)
args = [openmc_exec, '-p']
if path_input is not None:
args += [path_input]
_run(args, output, cwd)
def plot_inline(plots, openmc_exec='openmc', cwd='.'):
def plot_inline(plots, openmc_exec='openmc', cwd='.', path_input=None):
"""Display plots inline in a Jupyter notebook.
.. versionchanged:: 0.13.0
@ -155,6 +167,9 @@ def plot_inline(plots, openmc_exec='openmc', cwd='.'):
Path to OpenMC executable
cwd : str, optional
Path to working directory to run in
path_input : str
Path to a single XML file or a directory containing XML files for the
OpenMC executable to read.
Raises
------
@ -171,7 +186,7 @@ def plot_inline(plots, openmc_exec='openmc', cwd='.'):
openmc.Plots(plots).export_to_xml(cwd)
# Run OpenMC in geometry plotting mode
plot_geometry(False, openmc_exec, cwd)
plot_geometry(False, openmc_exec, cwd, path_input)
if plots is not None:
images = [_get_plot_image(p, cwd) for p in plots]
@ -179,7 +194,8 @@ def plot_inline(plots, openmc_exec='openmc', cwd='.'):
def calculate_volumes(threads=None, output=True, cwd='.',
openmc_exec='openmc', mpi_args=None):
openmc_exec='openmc', mpi_args=None,
path_input=None):
"""Run stochastic volume calculations in OpenMC.
This function runs OpenMC in stochastic volume calculation mode. To specify
@ -210,6 +226,10 @@ def calculate_volumes(threads=None, output=True, cwd='.',
cwd : str, optional
Path to working directory to run in. Defaults to the current working
directory.
path_input : str or Pathlike
Path to a single XML file or a directory containing XML files for the
OpenMC executable to read.
Raises
------
@ -223,14 +243,16 @@ def calculate_volumes(threads=None, output=True, cwd='.',
"""
args = _process_CLI_arguments(volume=True, threads=threads,
openmc_exec=openmc_exec, mpi_args=mpi_args)
openmc_exec=openmc_exec, mpi_args=mpi_args,
path_input=path_input)
_run(args, output, cwd)
def run(particles=None, threads=None, geometry_debug=False,
restart_file=None, tracks=False, output=True, cwd='.',
openmc_exec='openmc', mpi_args=None, event_based=False):
openmc_exec='openmc', mpi_args=None, event_based=False,
path_input=None):
"""Run an OpenMC simulation.
Parameters
@ -239,17 +261,17 @@ 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.
restart_file : str, optional
Path to restart file to use
tracks : bool, optional
Enables the writing of particles tracks. The number of particle
tracks written to tracks.h5 is limited to 1000 unless
Settings.max_tracks is set. Defaults to False.
Enables the writing of particles tracks. The number of particle tracks
written to tracks.h5 is limited to 1000 unless Settings.max_tracks is
set. Defaults to False.
output : bool
Capture OpenMC output from standard out
cwd : str, optional
@ -258,13 +280,17 @@ def run(particles=None, threads=None, geometry_debug=False,
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'].
MPI execute command and any additional MPI arguments to pass, e.g.
['mpiexec', '-n', '8'].
event_based : bool, optional
Turns on event-based parallelism, instead of default history-based
.. versionadded:: 0.12
path_input : str or Pathlike
Path to a single XML file or a directory containing XML files for the
OpenMC executable to read.
Raises
------
RuntimeError
@ -275,6 +301,7 @@ def run(particles=None, threads=None, geometry_debug=False,
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)
event_based=event_based, openmc_exec=openmc_exec, mpi_args=mpi_args,
path_input=path_input)
_run(args, output, cwd)

10
openmc/filter.py
View file

@ -533,7 +533,7 @@ class CellFromFilter(WithIDFilter):
expected_type = Cell
class CellbornFilter(WithIDFilter):
class CellBornFilter(WithIDFilter):
"""Bins tally events based on which cell the particle was born in.
Parameters
@ -557,6 +557,14 @@ class CellbornFilter(WithIDFilter):
expected_type = Cell
# Temporary alias for CellbornFilter
def CellbornFilter(*args, **kwargs):
warnings.warn('The name of "CellbornFilter" has changed to '
'"CellBornFilter". "CellbornFilter" will be '
'removed in the future.', FutureWarning)
return CellBornFilter(*args, **kwargs)
class CellInstanceFilter(Filter):
"""Bins tally events based on which cell instance a particle is in.

308
openmc/geometry.py
View file

@ -3,10 +3,11 @@ from collections.abc import Iterable
from copy import deepcopy
from pathlib import Path
from xml.etree import ElementTree as ET
import warnings
import openmc
import openmc._xml as xml
from .checkvalue import check_type
from .checkvalue import check_type, check_less_than, check_greater_than
class Geometry:
@ -25,12 +26,19 @@ class Geometry:
bounding_box : 2-tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
of the universe.
merge_surfaces : bool
Whether to remove redundant surfaces when the geometry is exported.
surface_precision : int
Number of decimal places to round to for comparing the coefficients of
surfaces for considering them topologically equivalent.
"""
def __init__(self, root=None):
self._root_universe = None
self._offsets = {}
self.merge_surfaces = False
self.surface_precision = 10
if root is not None:
if isinstance(root, openmc.UniverseBase):
self.root_universe = root
@ -48,11 +56,31 @@ class Geometry:
def bounding_box(self):
return self.root_universe.bounding_box
@property
def merge_surfaces(self):
return self._merge_surfaces
@property
def surface_precision(self):
return self._surface_precision
@root_universe.setter
def root_universe(self, root_universe):
check_type('root universe', root_universe, openmc.UniverseBase)
self._root_universe = root_universe
@merge_surfaces.setter
def merge_surfaces(self, merge_surfaces):
check_type('merge surfaces', merge_surfaces, bool)
self._merge_surfaces = merge_surfaces
@surface_precision.setter
def surface_precision(self, surface_precision):
check_type('surface precision', surface_precision, int)
check_less_than('surface_precision', surface_precision, 16)
check_greater_than('surface_precision', surface_precision, 0)
self._surface_precision = surface_precision
def add_volume_information(self, volume_calc):
"""Add volume information from a stochastic volume calculation.
@ -75,6 +103,39 @@ class Geometry:
if universe.id in volume_calc.volumes:
universe.add_volume_information(volume_calc)
def to_xml_element(self, remove_surfs=False):
"""Creates a 'geometry' element to be written to an XML file.
Parameters
----------
remove_surfs : bool
Whether or not to remove redundant surfaces from the geometry when
exporting
"""
# Find and remove redundant surfaces from the geometry
if remove_surfs:
warnings.warn("remove_surfs kwarg will be deprecated soon, please "
"set the Geometry.merge_surfaces attribute instead.")
self.merge_surfaces = True
if self.merge_surfaces:
self.remove_redundant_surfaces()
# Create XML representation
element = ET.Element("geometry")
self.root_universe.create_xml_subelement(element, memo=set())
# Sort the elements in the file
element[:] = sorted(element, key=lambda x: (
x.tag, int(x.get('id'))))
# Clean the indentation in the file to be user-readable
xml.clean_indentation(element)
xml.reorder_attributes(element) # TODO: Remove when support is Python 3.8+
return element
def export_to_xml(self, path='geometry.xml', remove_surfs=False):
"""Export geometry to an XML file.
@ -89,20 +150,7 @@ class Geometry:
.. versionadded:: 0.12
"""
# Find and remove redundant surfaces from the geometry
if remove_surfs:
self.remove_redundant_surfaces()
# Create XML representation
root_element = ET.Element("geometry")
self.root_universe.create_xml_subelement(root_element, memo=set())
# Sort the elements in the file
root_element[:] = sorted(root_element, key=lambda x: (
x.tag, int(x.get('id'))))
# Clean the indentation in the file to be user-readable
xml.clean_indentation(root_element)
root_element = self.to_xml_element(remove_surfs)
# Check if path is a directory
p = Path(path)
@ -110,10 +158,101 @@ class Geometry:
p /= 'geometry.xml'
# Write the XML Tree to the geometry.xml file
xml.reorder_attributes(root_element) # TODO: Remove when support is Python 3.8+
tree = ET.ElementTree(root_element)
tree.write(str(p), xml_declaration=True, encoding='utf-8')
@classmethod
def from_xml_element(cls, elem, materials=None):
"""Generate geometry from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
materials : openmc.Materials or None
Materials used to assign to cells. If None, an attempt is made to
generate it from the materials.xml file.
Returns
-------
openmc.Geometry
Geometry object
"""
mats = dict()
if materials is not None:
mats.update({str(m.id): m for m in materials})
mats['void'] = None
# Helper function for keeping a cache of Universe instances
universes = {}
def get_universe(univ_id):
if univ_id not in universes:
univ = openmc.Universe(univ_id)
universes[univ_id] = univ
return universes[univ_id]
# Get surfaces
surfaces = {}
periodic = {}
for surface in elem.findall('surface'):
s = openmc.Surface.from_xml_element(surface)
surfaces[s.id] = s
# Check for periodic surface
other_id = xml.get_text(surface, 'periodic_surface_id')
if other_id is not None:
periodic[s.id] = int(other_id)
# Apply periodic surfaces
for s1, s2 in periodic.items():
surfaces[s1].periodic_surface = surfaces[s2]
# Add any DAGMC universes
for e in elem.findall('dagmc_universe'):
dag_univ = openmc.DAGMCUniverse.from_xml_element(e)
universes[dag_univ.id] = dag_univ
# Dictionary that maps each universe to a list of cells/lattices that
# contain it (needed to determine which universe is the elem)
child_of = defaultdict(list)
for e in elem.findall('lattice'):
lat = openmc.RectLattice.from_xml_element(e, get_universe)
universes[lat.id] = lat
if lat.outer is not None:
child_of[lat.outer].append(lat)
for u in lat.universes.ravel():
child_of[u].append(lat)
for e in elem.findall('hex_lattice'):
lat = openmc.HexLattice.from_xml_element(e, get_universe)
universes[lat.id] = lat
if lat.outer is not None:
child_of[lat.outer].append(lat)
if lat.ndim == 2:
for ring in lat.universes:
for u in ring:
child_of[u].append(lat)
else:
for axial_slice in lat.universes:
for ring in axial_slice:
for u in ring:
child_of[u].append(lat)
for e in elem.findall('cell'):
c = openmc.Cell.from_xml_element(e, surfaces, mats, get_universe)
if c.fill_type in ('universe', 'lattice'):
child_of[c.fill].append(c)
# Determine which universe is the root by finding one which is not a
# child of any other object
for u in universes.values():
if not child_of[u]:
return cls(u)
else:
raise ValueError('Error determining root universe.')
@classmethod
def from_xml(cls, path='geometry.xml', materials=None):
"""Generate geometry from XML file
@ -132,84 +271,15 @@ class Geometry:
Geometry object
"""
# Helper function for keeping a cache of Universe instances
universes = {}
def get_universe(univ_id):
if univ_id not in universes:
univ = openmc.Universe(univ_id)
universes[univ_id] = univ
return universes[univ_id]
tree = ET.parse(path)
root = tree.getroot()
# Get surfaces
surfaces = {}
periodic = {}
for surface in root.findall('surface'):
s = openmc.Surface.from_xml_element(surface)
surfaces[s.id] = s
# Check for periodic surface
other_id = xml.get_text(surface, 'periodic_surface_id')
if other_id is not None:
periodic[s.id] = int(other_id)
# Apply periodic surfaces
for s1, s2 in periodic.items():
surfaces[s1].periodic_surface = surfaces[s2]
# Add any DAGMC universes
for elem in root.findall('dagmc_universe'):
dag_univ = openmc.DAGMCUniverse.from_xml_element(elem)
universes[dag_univ.id] = dag_univ
# Dictionary that maps each universe to a list of cells/lattices that
# contain it (needed to determine which universe is the root)
child_of = defaultdict(list)
for elem in root.findall('lattice'):
lat = openmc.RectLattice.from_xml_element(elem, get_universe)
universes[lat.id] = lat
if lat.outer is not None:
child_of[lat.outer].append(lat)
for u in lat.universes.ravel():
child_of[u].append(lat)
for elem in root.findall('hex_lattice'):
lat = openmc.HexLattice.from_xml_element(elem, get_universe)
universes[lat.id] = lat
if lat.outer is not None:
child_of[lat.outer].append(lat)
if lat.ndim == 2:
for ring in lat.universes:
for u in ring:
child_of[u].append(lat)
else:
for axial_slice in lat.universes:
for ring in axial_slice:
for u in ring:
child_of[u].append(lat)
# Create dictionary to easily look up materials
if materials is None:
filename = Path(path).parent / 'materials.xml'
materials = openmc.Materials.from_xml(str(filename))
mats = {str(m.id): m for m in materials}
mats['void'] = None
for elem in root.findall('cell'):
c = openmc.Cell.from_xml_element(elem, surfaces, mats, get_universe)
if c.fill_type in ('universe', 'lattice'):
child_of[c.fill].append(c)
tree = ET.parse(path)
root = tree.getroot()
# Determine which universe is the root by finding one which is not a
# child of any other object
for u in universes.values():
if not child_of[u]:
return cls(u)
else:
raise ValueError('Error determining root universe.')
return cls.from_xml_element(root, materials)
def find(self, point):
"""Find cells/universes/lattices which contain a given point
@ -396,28 +466,6 @@ class Geometry:
surfaces = cell.region.get_surfaces(surfaces)
return surfaces
def get_redundant_surfaces(self):
"""Return all of the topologically redundant surface IDs
.. versionadded:: 0.12
Returns
-------
dict
Dictionary whose keys are the ID of a redundant surface and whose
values are the topologically equivalent :class:`openmc.Surface`
that should replace it.
"""
tally = defaultdict(list)
for surf in self.get_all_surfaces().values():
coeffs = tuple(surf._coefficients[k] for k in surf._coeff_keys)
key = (surf._type,) + coeffs
tally[key].append(surf)
return {replace.id: keep
for keep, *redundant in tally.values()
for replace in redundant}
def _get_domains_by_name(self, name, case_sensitive, matching, domain_type):
if not case_sensitive:
name = name.lower()
@ -565,16 +613,42 @@ class Geometry:
return self._get_domains_by_name(name, case_sensitive, matching, 'lattice')
def remove_redundant_surfaces(self):
"""Remove redundant surfaces from the geometry"""
"""Remove and return all of the redundant surfaces.
Uses surface_precision attribute of Geometry instance for rounding and
comparing surface coefficients.
.. versionadded:: 0.12
Returns
-------
redundant_surfaces
Dictionary whose keys are the ID of a redundant surface and whose
values are the topologically equivalent :class:`openmc.Surface`
that should replace it.
"""
# Get redundant surfaces
redundant_surfaces = self.get_redundant_surfaces()
redundancies = defaultdict(list)
for surf in self.get_all_surfaces().values():
coeffs = tuple(round(surf._coefficients[k],
self.surface_precision)
for k in surf._coeff_keys)
key = (surf._type,) + coeffs
redundancies[key].append(surf)
# Iterate through all cells contained in the geometry
for cell in self.get_all_cells().values():
# Recursively remove redundant surfaces from regions
if cell.region:
cell.region.remove_redundant_surfaces(redundant_surfaces)
redundant_surfaces = {replace.id: keep
for keep, *redundant in redundancies.values()
for replace in redundant}
if redundant_surfaces:
# Iterate through all cells contained in the geometry
for cell in self.get_all_cells().values():
# Recursively remove redundant surfaces from regions
if cell.region:
cell.region.remove_redundant_surfaces(redundant_surfaces)
return redundant_surfaces
def determine_paths(self, instances_only=False):
"""Determine paths through CSG tree for cells and materials.

3
openmc/lib/__init__.py
View file

@ -51,6 +51,9 @@ def _coord_levels():
def _libmesh_enabled():
return c_bool.in_dll(_dll, "LIBMESH_ENABLED").value
def _mcpl_enabled():
return c_bool.in_dll(_dll, "MCPL_ENABLED").value
from .error import *
from .core import *
from .nuclide import *

162
openmc/material.py
View file

@ -6,7 +6,7 @@ from pathlib import Path
import re
import typing # imported separately as py3.8 requires typing.Iterable
import warnings
from typing import Optional, Union
from typing import Optional
from xml.etree import ElementTree as ET
import h5py
@ -98,7 +98,7 @@ class Material(IDManagerMixin):
this distribution is the total intensity of the photon source in
[decay/sec].
.. versionadded:: 0.14.0
.. versionadded:: 0.13.2
"""
@ -902,6 +902,8 @@ class Material(IDManagerMixin):
element : str
Specifies the element to match when searching through the nuclides
.. versionadded:: 0.13.2
Returns
-------
nuclides : list of str
@ -954,6 +956,8 @@ class Material(IDManagerMixin):
Nuclide for which atom density is desired. If not specified, the
atom density for each nuclide in the material is given.
.. versionadded:: 0.13.2
Returns
-------
nuclides : dict
@ -1041,7 +1045,7 @@ class Material(IDManagerMixin):
Returns
-------
Union[dict, float]
typing.Union[dict, float]
If by_nuclide is True then a dictionary whose keys are nuclide
names and values are activity is returned. Otherwise the activity
of the material is returned as a float.
@ -1063,7 +1067,50 @@ class Material(IDManagerMixin):
activity[nuclide] = inv_seconds * 1e24 * atoms_per_bcm * multiplier
return activity if by_nuclide else sum(activity.values())
def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False):
"""Returns the decay heat of the material or for each nuclide in the
material in units of [W], [W/g] or [W/cm3].
.. versionadded:: 0.13.3
Parameters
----------
units : {'W', 'W/g', 'W/cm3'}
Specifies the units of decay heat to return. Options include total
heat [W], specific [W/g] or volumetric heat [W/cm3].
Default is total heat [W].
by_nuclide : bool
Specifies if the decay heat should be returned for the material as a
whole or per nuclide. Default is False.
Returns
-------
Union[dict, float]
If `by_nuclide` is True then a dictionary whose keys are nuclide
names and values are decay heat is returned. Otherwise the decay heat
of the material is returned as a float.
"""
cv.check_value('units', units, {'W', 'W/g', 'W/cm3'})
cv.check_type('by_nuclide', by_nuclide, bool)
if units == 'W':
multiplier = self.volume
elif units == 'W/cm3':
multiplier = 1
elif units == 'W/g':
multiplier = 1.0 / self.get_mass_density()
decayheat = {}
for nuclide, atoms_per_bcm in self.get_nuclide_atom_densities().items():
decay_erg = openmc.data.decay_energy(nuclide)
inv_seconds = openmc.data.decay_constant(nuclide)
decay_erg *= openmc.data.JOULE_PER_EV
decayheat[nuclide] = inv_seconds * decay_erg * 1e24 * atoms_per_bcm * multiplier
return decayheat if by_nuclide else sum(decayheat.values())
def get_nuclide_atoms(self):
"""Return number of atoms of each nuclide in the material
@ -1487,6 +1534,57 @@ class Materials(cv.CheckedList):
for material in self:
material.make_isotropic_in_lab()
def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_indent=True):
"""Writes XML content of the materials to an open file handle.
Parameters
----------
file : IOTextWrapper
Open file handle to write content into.
header : bool
Whether or not to write the XML header
level : int
Indentation level of materials element
spaces_per_level : int
Number of spaces per indentation
trailing_indentation : bool
Whether or not to write a trailing indentation for the materials element
"""
indentation = level*spaces_per_level*' '
# Write the header and the opening tag for the root element.
if header:
file.write("<?xml version='1.0' encoding='utf-8'?>\n")
file.write(indentation+'<materials>\n')
# Write the <cross_sections> element.
if self.cross_sections is not None:
element = ET.Element('cross_sections')
element.text = str(self.cross_sections)
clean_indentation(element, level=level+1)
element.tail = element.tail.strip(' ')
file.write((level+1)*spaces_per_level*' ')
reorder_attributes(element) # TODO: Remove when support is Python 3.8+
ET.ElementTree(element).write(file, encoding='unicode')
# Write the <material> elements.
for material in sorted(self, key=lambda x: x.id):
element = material.to_xml_element()
clean_indentation(element, level=level+1)
element.tail = element.tail.strip(' ')
file.write((level+1)*spaces_per_level*' ')
reorder_attributes(element) # TODO: Remove when support is Python 3.8+
ET.ElementTree(element).write(file, encoding='unicode')
# Write the closing tag for the root element.
file.write(indentation+'</materials>\n')
# Write a trailing indentation for the next element
# at this level if needed
if trailing_indent:
file.write(indentation)
def export_to_xml(self, path: PathLike = 'materials.xml'):
"""Export material collection to an XML file.
@ -1506,32 +1604,34 @@ class Materials(cv.CheckedList):
# one go.
with open(str(p), 'w', encoding='utf-8',
errors='xmlcharrefreplace') as fh:
self._write_xml(fh)
# Write the header and the opening tag for the root element.
fh.write("<?xml version='1.0' encoding='utf-8'?>\n")
fh.write('<materials>\n')
@classmethod
def from_xml_element(cls, elem):
"""Generate materials collection from XML file
# Write the <cross_sections> element.
if self.cross_sections is not None:
element = ET.Element('cross_sections')
element.text = str(self.cross_sections)
clean_indentation(element, level=1)
element.tail = element.tail.strip(' ')
fh.write(' ')
reorder_attributes(element) # TODO: Remove when support is Python 3.8+
ET.ElementTree(element).write(fh, encoding='unicode')
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
# Write the <material> elements.
for material in sorted(self, key=lambda x: x.id):
element = material.to_xml_element()
clean_indentation(element, level=1)
element.tail = element.tail.strip(' ')
fh.write(' ')
reorder_attributes(element) # TODO: Remove when support is Python 3.8+
ET.ElementTree(element).write(fh, encoding='unicode')
Returns
-------
openmc.Materials
Materials collection
# Write the closing tag for the root element.
fh.write('</materials>\n')
"""
# Generate each material
materials = cls()
for material in elem.findall('material'):
materials.append(Material.from_xml_element(material))
# Check for cross sections settings
xs = elem.find('cross_sections')
if xs is not None:
materials.cross_sections = xs.text
return materials
@classmethod
def from_xml(cls, path: PathLike = 'materials.xml'):
@ -1551,14 +1651,4 @@ class Materials(cv.CheckedList):
tree = ET.parse(path)
root = tree.getroot()
# Generate each material
materials = cls()
for material in root.findall('material'):
materials.append(Material.from_xml_element(material))
# Check for cross sections settings
xs = tree.find('cross_sections')
if xs is not None:
materials.cross_sections = xs.text
return materials
return cls.from_xml_element(root)

15
openmc/mesh.py
View file

@ -172,10 +172,11 @@ class StructuredMesh(MeshBase):
-------
vertices : numpy.ndarray
Returns a numpy.ndarray representing the coordinates of the mesh
vertices with a shape equal to (dim1 + 1, ..., dimn + 1, ndim).
vertices with a shape equal to (ndim, dim1 + 1, ..., dimn + 1). Can be
unpacked along the first dimension with xx, yy, zz = mesh.vertices.
"""
return np.stack(np.meshgrid(*self._grids, indexing='ij'), axis=-1)
return np.stack(np.meshgrid(*self._grids, indexing='ij'), axis=0)
@property
def centroids(self):
@ -185,13 +186,15 @@ class StructuredMesh(MeshBase):
-------
centroids : numpy.ndarray
Returns a numpy.ndarray representing the mesh element centroid
coordinates with a shape equal to (dim1, ..., dimn, ndim).
coordinates with a shape equal to (ndim, dim1, ..., dimn). Can be
unpacked along the first dimension with xx, yy, zz = mesh.centroids.
"""
ndim = self.n_dimension
vertices = self.vertices
s0 = (slice(0, -1),)*ndim + (slice(None),)
s1 = (slice(1, None),)*ndim + (slice(None),)
s0 = (slice(None),) + (slice(0, -1),)*ndim
s1 = (slice(None),) + (slice(1, None),)*ndim
return (vertices[s0] + vertices[s1]) / 2
@property
@ -324,7 +327,7 @@ class RegularMesh(StructuredMesh):
@property
def dimension(self):
return self._dimension
return tuple(self._dimension)
@property
def n_dimension(self):

4
openmc/mgxs/mgxs.py
View file

@ -925,7 +925,7 @@ class MGXS:
# Tabulate the atomic number densities for all nuclides
elif nuclides == 'all':
nuclides = self.get_nuclides()
densities = np.zeros(self.num_nuclides, dtype=np.float)
densities = np.zeros(self.num_nuclides, dtype=float)
for i, nuclide in enumerate(nuclides):
densities[i] += self.get_nuclide_density(nuclide)
@ -3019,7 +3019,7 @@ class DiffusionCoefficient(TransportXS):
new_filt = openmc.EnergyFilter(old_filt.values)
p1_tally.filters[-2] = new_filt
p1_tally = p1_tally.get_slice(filters=[openmc.LegendreFilter],
p1_tally = p1_tally.get_slice(filters=[openmc.LegendreFilter],
filter_bins=[('P1',)],squeeze=True)
p1_tally._scores = ['scatter-1']
total_xs = self.tallies['total'] / self.tallies['flux (tracklength)']

106
openmc/model/model.py
View file

@ -5,10 +5,13 @@ import os
from pathlib import Path
from numbers import Integral
from tempfile import NamedTemporaryFile
import warnings
from xml.etree import ElementTree as ET
import h5py
import openmc
import openmc._xml as xml
from openmc.dummy_comm import DummyCommunicator
from openmc.executor import _process_CLI_arguments
from openmc.checkvalue import check_type, check_value
@ -237,6 +240,35 @@ class Model:
plots = openmc.Plots.from_xml(plots) if Path(plots).exists() else None
return cls(geometry, materials, settings, tallies, plots)
@classmethod
def from_model_xml(cls, path='model.xml'):
"""Create model from single XML file
.. vesionadded:: 0.13.3
Parameters
----------
path : str or Pathlike
Path to model.xml file
"""
tree = ET.parse(path)
root = tree.getroot()
model = cls()
meshes = {}
model.settings = openmc.Settings.from_xml_element(root.find('settings'), meshes)
model.materials = openmc.Materials.from_xml_element(root.find('materials'))
model.geometry = openmc.Geometry.from_xml_element(root.find('geometry'), model.materials)
if root.find('tallies'):
model.tallies = openmc.Tallies.from_xml_element(root.find('tallies'), meshes)
if root.find('plots'):
model.plots = openmc.Plots.from_xml_element(root.find('plots'))
return model
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
@ -398,7 +430,7 @@ class Model:
depletion_operator.finalize()
def export_to_xml(self, directory='.', remove_surfs=False):
"""Export model to XML files.
"""Export model to separate XML files.
Parameters
----------
@ -434,6 +466,78 @@ class Model:
if self.plots:
self.plots.export_to_xml(d)
def export_to_model_xml(self, path='model.xml', remove_surfs=False):
"""Export model to a single XML file.
.. versionadded:: 0.13.3
Parameters
----------
path : str or Pathlike
Location of the XML file to write (default is 'model.xml'). Can be a
directory or file path.
remove_surfs : bool
Whether or not to remove redundant surfaces from the geometry when
exporting.
"""
xml_path = Path(path)
# if the provided path doesn't end with the XML extension, assume the
# input path is meant to be a directory. If the directory does not
# exist, create it and place a 'model.xml' file there.
if not str(xml_path).endswith('.xml') and not xml_path.exists():
os.mkdir(xml_path)
xml_path /= 'model.xml'
# if this is an XML file location and the file's parent directory does
# not exist, create it before continuing
elif not xml_path.parent.exists():
os.mkdir(xml_path.parent)
if remove_surfs:
warnings.warn("remove_surfs kwarg will be deprecated soon, please "
"set the Geometry.merge_surfaces attribute instead.")
self.geometry.merge_surfaces = True
# Can be used to modify tallies in case any surfaces are redundant
redundant_surfaces = self.geometry.remove_redundant_surfaces()
# provide a memo to track which meshes have been written
mesh_memo = set()
settings_element = self.settings.to_xml_element(mesh_memo)
geometry_element = self.geometry.to_xml_element()
xml.clean_indentation(geometry_element, level=1)
xml.clean_indentation(settings_element, level=1)
# If a materials collection was specified, export it. Otherwise, look
# for all materials in the geometry and use that to automatically build
# a collection.
if self.materials:
materials = self.materials
else:
materials = openmc.Materials(self.geometry.get_all_materials()
.values())
with open(xml_path, 'w', encoding='utf-8', errors='xmlcharrefreplace') as fh:
# write the XML header
fh.write("<?xml version='1.0' encoding='utf-8'?>\n")
fh.write("<model>\n")
# Write the materials collection to the open XML file first.
# This will write the XML header also
materials._write_xml(fh, False, level=1)
# Write remaining elements as a tree
ET.ElementTree(geometry_element).write(fh, encoding='unicode')
ET.ElementTree(settings_element).write(fh, encoding='unicode')
if self.tallies:
tallies_element = self.tallies.to_xml_element(mesh_memo)
xml.clean_indentation(tallies_element, level=1, trailing_indent=self.plots)
ET.ElementTree(tallies_element).write(fh, encoding='unicode')
if self.plots:
plots_element = self.plots.to_xml_element()
xml.clean_indentation(plots_element, level=1, trailing_indent=False)
ET.ElementTree(plots_element).write(fh, encoding='unicode')
fh.write("</model>\n")
def import_properties(self, filename):
"""Import physical properties

340
openmc/model/surface_composite.py
View file

@ -1,10 +1,16 @@
from abc import ABC, abstractmethod
from copy import copy
from math import sqrt, pi, sin, cos
from math import sqrt, pi, sin, cos, isclose
import warnings
import operator
import numpy as np
from scipy.spatial import ConvexHull, Delaunay
from matplotlib.path import Path
import openmc
from openmc.checkvalue import check_greater_than, check_value
from openmc.checkvalue import (check_greater_than, check_value,
check_iterable_type, check_length)
class CompositeSurface(ABC):
@ -621,3 +627,333 @@ class ZConeOneSided(CompositeSurface):
__neg__ = XConeOneSided.__neg__
__pos__ = XConeOneSided.__pos__
class Polygon(CompositeSurface):
"""Create a polygon composite surface from a path of closed points.
Parameters
----------
points : np.ndarray
An Nx2 array of points defining the vertices of the polygon.
basis : {'rz', 'xy', 'yz', 'xz'}, optional
2D basis set for the polygon. The polygon is two dimensional and has
infinite extent in the third (unspecified) dimension. For example, the
'xy' basis produces a polygon with infinite extent in the +/- z
direction. For the 'rz' basis the phi extent is infinite, thus forming
an axisymmetric surface.
Attributes
----------
points : np.ndarray
An Nx2 array of points defining the vertices of the polygon.
basis : {'rz', 'xy', 'yz', 'xz'}
2D basis set for the polygon.
regions : list of openmc.Region
A list of :class:`openmc.Region` objects, one for each of the convex polygons
formed during the decomposition of the input polygon.
region : openmc.Union
The union of all the regions comprising the polygon.
"""
def __init__(self, points, basis='rz'):
check_value('basis', basis, ('xy', 'yz', 'xz', 'rz'))
self._basis = basis
points = np.asarray(points, dtype=float)
check_iterable_type('points', points, float, min_depth=2, max_depth=2)
check_length('points', points[0, :], 2, 2)
# If the last point is the same as the first, remove it and make sure
# there are still at least 3 points for a valid polygon.
if np.allclose(points[0, :], points[-1, :]):
points = points[:-1, :]
check_length('points', points, 3)
# Order the points counter-clockwise (necessary for offset method)
self._points = self._make_ccw(points)
# Create a triangulation of the points.
self._tri = Delaunay(self._points, qhull_options='QJ')
# Decompose the polygon into groups of simplices forming convex subsets
# and get the sets of (surface, operator) pairs defining the polygon
self._surfsets = self._decompose_polygon_into_convex_sets()
# Set surface names as required by CompositeSurface protocol
surfnames = []
i = 0
for surfset in self._surfsets:
for surf, op, on_boundary in surfset:
if on_boundary:
setattr(self, f'surface_{i}', surf)
surfnames.append(f'surface_{i}')
i += 1
self._surfnames = tuple(surfnames)
# Generate a list of regions whose union represents the polygon.
regions = []
for surfs_ops in self._surfsets:
regions.append([op(surf) for surf, op, _ in surfs_ops])
self._regions = [openmc.Intersection(regs) for regs in regions]
# Create the union of all the convex subsets
self._region = openmc.Union(self._regions)
def __neg__(self):
return self._region
def __pos__(self):
return ~self._region
@property
def _surface_names(self):
return self._surfnames
@CompositeSurface.boundary_type.setter
def boundary_type(self, boundary_type):
if boundary_type != 'transmission':
warnings.warn("Setting boundary_type to a value other than "
"'transmission' on Polygon composite surfaces can "
"result in unintended behavior. Please use the "
"regions property of the Polygon to generate "
"individual openmc.Cell objects to avoid unwanted "
"behavior.")
for name in self._surface_names:
getattr(self, name).boundary_type = boundary_type
@property
def points(self):
return self._points
@property
def basis(self):
return self._basis
@property
def _normals(self):
"""Generate the outward normal unit vectors for the polygon."""
# Rotation matrix for 90 degree clockwise rotation (-90 degrees about z
# axis for an 'xy' basis).
rotation = np.array([[0., 1.], [-1., 0.]])
# Get the unit vectors that point from one point in the polygon to the
# next given that they are ordered counterclockwise and that the final
# point is connected to the first point
tangents = np.diff(self._points, axis=0, append=[self._points[0, :]])
tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True)
# Rotate the tangent vectors clockwise by 90 degrees, which for a
# counter-clockwise ordered polygon will produce the outward normal
# vectors.
return rotation.dot(tangents.T).T
@property
def _equations(self):
normals = self._normals
equations = np.empty((normals.shape[0], 3))
equations[:, :2] = normals
equations[:, 2] = -np.sum(normals*self.points, axis=-1)
return equations
@property
def regions(self):
return self._regions
@property
def region(self):
return self._region
def _make_ccw(self, points):
"""Order a set of points counter-clockwise.
Parameters
----------
points : np.ndarray (Nx2)
An Nx2 array of coordinate pairs describing the vertices.
Returns
-------
ordered_points : the input points ordered counter-clockwise
"""
# Calculates twice the signed area of the polygon using the "Shoelace
# Formula" https://en.wikipedia.org/wiki/Shoelace_formula
xpts, ypts = points.T
# If signed area is positive the curve is oriented counter-clockwise
if np.sum(ypts*(np.roll(xpts, 1) - np.roll(xpts, -1))) > 0:
return points
return points[::-1, :]
def _group_simplices(self, neighbor_map, group=None):
"""Generate a convex grouping of simplices.
Parameters
----------
neighbor_map : dict
A map whose keys are simplex indices for simplices inside the polygon
and whose values are a list of simplex indices that neighbor this
simplex and are also inside the polygon.
group : list
A list of simplex indices that comprise the current convex group.
Returns
-------
group : list
The list of simplex indices that comprise the complete convex group.
"""
# If neighbor_map is empty there's nothing left to do
if not neighbor_map:
return group
# If group is empty, grab the next simplex in the dictionary and recurse
if group is None:
sidx = next(iter(neighbor_map))
return self._group_simplices(neighbor_map, group=[sidx])
# Otherwise use the last simplex in the group
else:
sidx = group[-1]
# Remove current simplex from dictionary since it is in a group
neighbors = neighbor_map.pop(sidx, [])
# For each neighbor check if it is part of the same convex
# hull as the rest of the group. If yes, recurse. If no, continue on.
for n in neighbors:
if n in group or neighbor_map.get(n, None) is None:
continue
test_group = group + [n]
test_point_idx = np.unique(self._tri.simplices[test_group, :])
test_points = self._tri.points[test_point_idx]
# If test_points are convex keep adding to this group
if len(test_points) == len(ConvexHull(test_points).vertices):
group = self._group_simplices(neighbor_map, group=test_group)
return group
def _get_convex_hull_surfs(self, qhull):
"""Generate a list of surfaces given by a set of linear equations
Parameters
----------
qhull : scipy.spatial.ConvexHull
A ConvexHull object representing the sub-region of the polygon.
Returns
-------
surfs_ops : list of (surface, operator) tuples
"""
basis = self.basis
boundary_eqns = self._equations
# Collect surface/operator pairs such that the intersection of the
# regions defined by these pairs is the inside of the polygon.
surfs_ops = []
# hull facet equation: dx*x + dy*y + c = 0
for dx, dy, c in qhull.equations:
# check if this facet is on the boundary of the polygon
facet_eq = np.array([dx, dy, c])
on_boundary = any([np.allclose(facet_eq, eq) for eq in boundary_eqns])
# Check if the facet is horizontal
if isclose(dx, 0, abs_tol=1e-8):
if basis in ('xz', 'yz', 'rz'):
surf = openmc.ZPlane(z0=-c/dy)
else:
surf = openmc.YPlane(y0=-c/dy)
# if (0, 1).(dx, dy) < 0 we want positive halfspace instead
op = operator.pos if dy < 0 else operator.neg
# Check if the facet is vertical
elif isclose(dy, 0, abs_tol=1e-8):
if basis in ('xy', 'xz'):
surf = openmc.XPlane(x0=-c/dx)
elif basis == 'yz':
surf = openmc.YPlane(y0=-c/dx)
else:
surf = openmc.ZCylinder(r=-c/dx)
# if (1, 0).(dx, dy) < 0 we want positive halfspace instead
op = operator.pos if dx < 0 else operator.neg
# Otherwise the facet is at an angle
else:
op = operator.neg
if basis == 'xy':
surf = openmc.Plane(a=dx, b=dy, d=-c)
elif basis == 'yz':
surf = openmc.Plane(b=dx, c=dy, d=-c)
elif basis == 'xz':
surf = openmc.Plane(a=dx, c=dy, d=-c)
else:
y0 = -c/dy
r2 = dy**2 / dx**2
# Check if the *slope* of the facet is positive. If dy/dx < 0
# then we want up to be True for the one-sided cones.
up = dy / dx < 0
surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=up)
# if (1, -1).(dx, dy) < 0 for up cones we want positive halfspace
# if (1, 1).(dx, dy) < 0 for down cones we want positive halfspace
# otherwise we keep the negative halfspace operator
if (up and dx - dy < 0) or (not up and dx + dy < 0):
op = operator.pos
surfs_ops.append((surf, op, on_boundary))
return surfs_ops
def _decompose_polygon_into_convex_sets(self):
"""Decompose the Polygon into a set of convex polygons.
Returns
-------
surfsets : a list of lists of surface, operator pairs
"""
# Get centroids of all the simplices and determine if they are inside
# the polygon defined by input vertices or not.
centroids = np.mean(self._points[self._tri.simplices], axis=1)
in_polygon = Path(self._points).contains_points(centroids)
self._in_polygon = in_polygon
# Build a map with keys of simplex indices inside the polygon whose
# values are lists of that simplex's neighbors also inside the
# polygon
neighbor_map = {}
for i, nlist in enumerate(self._tri.neighbors):
if not in_polygon[i]:
continue
neighbor_map[i] = [n for n in nlist if in_polygon[n] and n >=0]
# Get the groups of simplices forming convex polygons whose union
# comprises the full input polygon. While there are still simplices
# left in the neighbor map, group them together into convex sets.
groups = []
while neighbor_map:
groups.append(self._group_simplices(neighbor_map))
self._groups = groups
# Generate lists of (surface, operator) pairs for each convex
# sub-region.
surfsets = []
for group in groups:
# Find all the unique points in the convex group of simplices,
# generate the convex hull and find the resulting surfaces and
# unary operators that represent this convex subset of the polygon.
idx = np.unique(self._tri.simplices[group, :])
qhull = ConvexHull(self._tri.points[idx, :])
surf_ops = self._get_convex_hull_surfs(qhull)
surfsets.append(surf_ops)
return surfsets
def offset(self, distance):
"""Offset this polygon by a set distance
Parameters
----------
distance : float
The distance to offset the polygon by. Positive is outward
(expanding) and negative is inward (shrinking).
Returns
-------
offset_polygon : openmc.model.Polygon
"""
normals = np.insert(self._normals, 0, self._normals[-1, :], axis=0)
cos2theta = np.sum(normals[1:, :]*normals[:-1, :], axis=-1, keepdims=True)
costheta = np.cos(np.arccos(cos2theta) / 2)
nvec = (normals[1:, :] + normals[:-1, :])
unit_nvec = nvec / np.linalg.norm(nvec, axis=-1, keepdims=True)
disp_vec = distance / costheta * unit_nvec
return type(self)(self.points + disp_vec, basis=self.basis)

2
openmc/nuclide.py
View file

@ -26,7 +26,7 @@ class Nuclide(str):
if name.endswith('m'):
name = name[:-1] + '_m1'
msg = ('OpenMC nuclides follow the GND naming convention. '
msg = ('OpenMC nuclides follow the GNDS naming convention. '
f'Nuclide "{orig_name}" is being renamed as "{name}".')
warnings.warn(msg)

54
openmc/plots.py
View file

@ -909,13 +909,13 @@ class Plots(cv.CheckedList):
self._plots_file.append(xml_element)
def export_to_xml(self, path='plots.xml'):
"""Export plot specifications to an XML file.
def to_xml_element(self):
"""Create a 'plots' element to be written to an XML file.
Parameters
----------
path : str
Path to file to write. Defaults to 'plots.xml'.
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing all plot elements
"""
# Reset xml element tree
@ -925,17 +925,50 @@ class Plots(cv.CheckedList):
# Clean the indentation in the file to be user-readable
clean_indentation(self._plots_file)
reorder_attributes(self._plots_file) # TODO: Remove when support is Python 3.8+
return self._plots_file
def export_to_xml(self, path='plots.xml'):
"""Export plot specifications to an XML file.
Parameters
----------
path : str
Path to file to write. Defaults to 'plots.xml'.
"""
# Check if path is a directory
p = Path(path)
if p.is_dir():
p /= 'plots.xml'
self.to_xml_element()
# Write the XML Tree to the plots.xml file
reorder_attributes(self._plots_file) # TODO: Remove when support is Python 3.8+
tree = ET.ElementTree(self._plots_file)
tree.write(str(p), xml_declaration=True, encoding='utf-8')
@classmethod
def from_xml_element(cls, elem):
"""Generate plots collection from XML file
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.Plots
Plots collection
"""
# Generate each plot
plots = cls()
for e in elem.findall('plot'):
plots.append(Plot.from_xml_element(e))
return plots
@classmethod
def from_xml(cls, path='plots.xml'):
"""Generate plots collection from XML file
@ -953,9 +986,6 @@ class Plots(cv.CheckedList):
"""
tree = ET.parse(path)
root = tree.getroot()
return cls.from_xml_element(root)
# Generate each plot
plots = cls()
for elem in root.findall('plot'):
plots.append(Plot.from_xml_element(elem))
return plots

372
openmc/settings.py
View file

@ -6,7 +6,8 @@ import itertools
from math import ceil
from numbers import Integral, Real
from pathlib import Path
from typing import Optional, Union
import typing # required to prevent typing.Union namespace overwriting Union
from typing import Optional
from xml.etree import ElementTree as ET
import openmc.checkvalue as cv
@ -156,6 +157,7 @@ class Settings:
:separate: bool indicating whether the source should be written as a
separate file
:write: bool indicating whether or not to write the source
:mcpl: bool indicating whether to write the source as an MCPL file
statepoint : dict
Options for writing state points. Acceptable keys are:
@ -171,6 +173,7 @@ class Settings:
banked (int)
:max_particles: Maximum number of particles to be banked on
surfaces per process (int)
:mcpl: Output in the form of an MCPL-file (bool)
survival_biasing : bool
Indicate whether survival biasing is to be used
tabular_legendre : dict
@ -186,13 +189,15 @@ class Settings:
'default', 'method', 'range', 'tolerance', and 'multipole'. The value
for 'default' should be a float representing the default temperature in
Kelvin. The value for 'method' should be 'nearest' or 'interpolation'.
If the method is 'nearest', 'tolerance' indicates a range of temperature
within which cross sections may be used. The value for 'range' should be
a pair a minimum and maximum temperatures which are used to indicate
that cross sections be loaded at all temperatures within the
range. 'multipole' is a boolean indicating whether or not the windowed
multipole method should be used to evaluate resolved resonance cross
sections.
If the method is 'nearest', 'tolerance' indicates a range of
temperature within which cross sections may be used. If the method is
'interpolation', 'tolerance' indicates the range of temperatures outside
of the available cross section temperatures where cross sections will
evaluate to the nearer bound. The value for 'range' should be a pair of
minimum and maximum temperatures which are used to indicate that cross
sections be loaded at all temperatures within the range. 'multipole' is
a boolean indicating whether or not the windowed multipole method should
be used to evaluate resolved resonance cross sections.
trace : tuple or list
Show detailed information about a single particle, indicated by three
integers: the batch number, generation number, and particle number
@ -582,7 +587,7 @@ class Settings:
self._max_order = max_order
@source.setter
def source(self, source: Union[Source, typing.Iterable[Source]]):
def source(self, source: typing.Union[Source, typing.Iterable[Source]]):
if not isinstance(source, MutableSequence):
source = [source]
self._source = cv.CheckedList(Source, 'source distributions', source)
@ -619,6 +624,8 @@ class Settings:
cv.check_type('sourcepoint write', value, bool)
elif key == 'overwrite':
cv.check_type('sourcepoint overwrite', value, bool)
elif key == 'mcpl':
cv.check_type('sourcepoint mcpl', value, bool)
else:
raise ValueError(f"Unknown key '{key}' encountered when "
"setting sourcepoint options.")
@ -652,7 +659,7 @@ class Settings:
cv.check_type('surface source writing options', surf_source_write, Mapping)
for key, value in surf_source_write.items():
cv.check_value('surface source writing key', key,
('surface_ids', 'max_particles'))
('surface_ids', 'max_particles', 'mcpl'))
if key == 'surface_ids':
cv.check_type('surface ids for source banking', value,
Iterable, Integral)
@ -664,6 +671,9 @@ class Settings:
value, Integral)
cv.check_greater_than('maximum particle banks on surfaces per process',
value, 0)
elif key == 'mcpl':
cv.check_type('write to an MCPL-format file', value, bool)
self._surf_source_write = surf_source_write
@confidence_intervals.setter
@ -843,7 +853,7 @@ class Settings:
@volume_calculations.setter
def volume_calculations(
self, vol_calcs: Union[VolumeCalculation, typing.Iterable[VolumeCalculation]]
self, vol_calcs: typing.Union[VolumeCalculation, typing.Iterable[VolumeCalculation]]
):
if not isinstance(vol_calcs, MutableSequence):
vol_calcs = [vol_calcs]
@ -889,7 +899,7 @@ class Settings:
self._write_initial_source = value
@weight_windows.setter
def weight_windows(self, value: Union[WeightWindows, typing.Iterable[WeightWindows]]):
def weight_windows(self, value: typing.Union[WeightWindows, typing.Iterable[WeightWindows]]):
if not isinstance(value, MutableSequence):
value = [value]
self._weight_windows = cv.CheckedList(WeightWindows, 'weight windows', value)
@ -1015,6 +1025,10 @@ class Settings:
subelement = ET.SubElement(element, "overwrite_latest")
subelement.text = str(self._sourcepoint['overwrite']).lower()
if 'mcpl' in self._sourcepoint:
subelement = ET.SubElement(element, "mcpl")
subelement.text = str(self._sourcepoint['mcpl']).lower()
def _create_surf_source_read_subelement(self, root):
if self._surf_source_read:
element = ET.SubElement(root, "surf_source_read")
@ -1032,6 +1046,9 @@ class Settings:
if 'max_particles' in self._surf_source_write:
subelement = ET.SubElement(element, "max_particles")
subelement.text = str(self._surf_source_write['max_particles'])
if 'mcpl' in self._surf_source_write:
subelement = ET.SubElement(element, "mcpl")
subelement.text = str(self._surf_source_write['mcpl']).lower()
def _create_confidence_intervals(self, root):
if self._confidence_intervals is not None:
@ -1070,25 +1087,35 @@ class Settings:
subelement = ET.SubElement(element, key)
subelement.text = str(value)
def _create_entropy_mesh_subelement(self, root):
if self.entropy_mesh is not None:
# use default heuristic for entropy mesh if not set by user
if self.entropy_mesh.dimension is None:
if self.particles is None:
raise RuntimeError("Number of particles must be set in order to " \
"use entropy mesh dimension heuristic")
else:
n = ceil((self.particles / 20.0)**(1.0 / 3.0))
d = len(self.entropy_mesh.lower_left)
self.entropy_mesh.dimension = (n,)*d
def _create_entropy_mesh_subelement(self, root, mesh_memo=None):
if self.entropy_mesh is None:
return
# See if a <mesh> element already exists -- if not, add it
path = f"./mesh[@id='{self.entropy_mesh.id}']"
if root.find(path) is None:
root.append(self.entropy_mesh.to_xml_element())
# use default heuristic for entropy mesh if not set by user
if self.entropy_mesh.dimension is None:
if self.particles is None:
raise RuntimeError("Number of particles must be set in order to " \
"use entropy mesh dimension heuristic")
else:
n = ceil((self.particles / 20.0)**(1.0 / 3.0))
d = len(self.entropy_mesh.lower_left)
self.entropy_mesh.dimension = (n,)*d
subelement = ET.SubElement(root, "entropy_mesh")
subelement.text = str(self.entropy_mesh.id)
# add mesh ID to this element
subelement = ET.SubElement(root, "entropy_mesh")
subelement.text = str(self.entropy_mesh.id)
# If this mesh has already been written outside the
# settings element, skip writing it again
if mesh_memo and self.entropy_mesh.id in mesh_memo:
return
# See if a <mesh> element already exists -- if not, add it
path = f"./mesh[@id='{self.entropy_mesh.id}']"
if root.find(path) is None:
root.append(self.entropy_mesh.to_xml_element())
if mesh_memo is not None:
mesh_memo.add(self.entropy_mesh.id)
def _create_trigger_subelement(self, root):
if self._trigger_active is not None:
@ -1139,15 +1166,21 @@ class Settings:
element = ET.SubElement(root, "track")
element.text = ' '.join(map(str, itertools.chain(*self._track)))
def _create_ufs_mesh_subelement(self, root):
if self.ufs_mesh is not None:
# See if a <mesh> element already exists -- if not, add it
path = f"./mesh[@id='{self.ufs_mesh.id}']"
if root.find(path) is None:
root.append(self.ufs_mesh.to_xml_element())
def _create_ufs_mesh_subelement(self, root, mesh_memo=None):
if self.ufs_mesh is None:
return
subelement = ET.SubElement(root, "ufs_mesh")
subelement.text = str(self.ufs_mesh.id)
subelement = ET.SubElement(root, "ufs_mesh")
subelement.text = str(self.ufs_mesh.id)
if mesh_memo and self.ufs_mesh.id in mesh_memo:
return
# See if a <mesh> element already exists -- if not, add it
path = f"./mesh[@id='{self.ufs_mesh.id}']"
if root.find(path) is None:
root.append(self.ufs_mesh.to_xml_element())
if mesh_memo is not None: mesh_memo.add(self.ufs_mesh.id)
def _create_resonance_scattering_subelement(self, root):
res = self.resonance_scattering
@ -1204,15 +1237,21 @@ class Settings:
elem = ET.SubElement(root, "write_initial_source")
elem.text = str(self._write_initial_source).lower()
def _create_weight_windows_subelement(self, root):
def _create_weight_windows_subelement(self, root, mesh_memo=None):
for ww in self._weight_windows:
# Add weight window information
root.append(ww.to_xml_element())
# if this mesh has already been written,
# skip writing the mesh element
if mesh_memo and ww.mesh.id in mesh_memo:
continue
# See if a <mesh> element already exists -- if not, add it
path = f"./mesh[@id='{ww.mesh.id}']"
if root.find(path) is None:
root.append(ww.mesh.to_xml_element())
if mesh_memo is not None: mesh_memo.add(ww.mesh.id)
if self._weight_windows_on is not None:
elem = ET.SubElement(root, "weight_windows_on")
@ -1312,10 +1351,10 @@ class Settings:
def _sourcepoint_from_xml_element(self, root):
elem = root.find('source_point')
if elem is not None:
for key in ('separate', 'write', 'overwrite_latest', 'batches'):
for key in ('separate', 'write', 'overwrite_latest', 'batches', 'mcpl'):
value = get_text(elem, key)
if value is not None:
if key in ('separate', 'write'):
if key in ('separate', 'write', 'mcpl'):
value = value in ('true', '1')
elif key == 'overwrite_latest':
value = value in ('true', '1')
@ -1334,13 +1373,15 @@ class Settings:
def _surf_source_write_from_xml_element(self, root):
elem = root.find('surf_source_write')
if elem is not None:
for key in ('surface_ids', 'max_particles'):
for key in ('surface_ids', 'max_particles','mcpl'):
value = get_text(elem, key)
if value is not None:
if key == 'surface_ids':
value = [int(x) for x in value.split()]
elif key in ('max_particles'):
value = int(value)
elif key == 'mcpl':
value = value in ('true', '1')
self.surf_source_write[key] = value
def _confidence_intervals_from_xml_element(self, root):
@ -1393,13 +1434,15 @@ class Settings:
if value is not None:
self.cutoff[key] = float(value)
def _entropy_mesh_from_xml_element(self, root):
def _entropy_mesh_from_xml_element(self, root, meshes=None):
text = get_text(root, 'entropy_mesh')
if text is not None:
path = f"./mesh[@id='{int(text)}']"
elem = root.find(path)
if elem is not None:
self.entropy_mesh = RegularMesh.from_xml_element(elem)
if meshes is not None and self.entropy_mesh is not None:
meshes[self.entropy_mesh.id] = self.entropy_mesh
def _trigger_from_xml_element(self, root):
elem = root.find('trigger')
@ -1459,13 +1502,15 @@ class Settings:
values = [int(x) for x in text.split()]
self.track = list(zip(values[::3], values[1::3], values[2::3]))
def _ufs_mesh_from_xml_element(self, root):
def _ufs_mesh_from_xml_element(self, root, meshes=None):
text = get_text(root, 'ufs_mesh')
if text is not None:
path = f"./mesh[@id='{int(text)}']"
elem = root.find(path)
if elem is not None:
self.ufs_mesh = RegularMesh.from_xml_element(elem)
if meshes is not None and self.ufs_mesh is not None:
meshes[self.ufs_mesh.id] = self.ufs_mesh
def _resonance_scattering_from_xml_element(self, root):
elem = root.find('resonance_scattering')
@ -1517,7 +1562,7 @@ class Settings:
if text is not None:
self.write_initial_source = text in ('true', '1')
def _weight_windows_from_xml_element(self, root):
def _weight_windows_from_xml_element(self, root, meshes=None):
for elem in root.findall('weight_windows'):
ww = WeightWindows.from_xml_element(elem, root)
self.weight_windows.append(ww)
@ -1526,6 +1571,9 @@ class Settings:
if text is not None:
self.weight_windows_on = text in ('true', '1')
if meshes is not None and self.weight_windows:
meshes.update({ww.mesh.id: ww.mesh for ww in self.weight_windows})
def _max_splits_from_xml_element(self, root):
text = get_text(root, 'max_splits')
if text is not None:
@ -1536,6 +1584,68 @@ class Settings:
if text is not None:
self.max_tracks = int(text)
def to_xml_element(self, mesh_memo=None):
"""Create a 'settings' element to be written to an XML file.
Parameters
----------
mesh_memo : set of ints
A set of mesh IDs to keep track of whether a mesh has already been written.
"""
# Reset xml element tree
element = ET.Element("settings")
self._create_run_mode_subelement(element)
self._create_particles_subelement(element)
self._create_batches_subelement(element)
self._create_inactive_subelement(element)
self._create_max_lost_particles_subelement(element)
self._create_rel_max_lost_particles_subelement(element)
self._create_generations_per_batch_subelement(element)
self._create_keff_trigger_subelement(element)
self._create_source_subelement(element)
self._create_output_subelement(element)
self._create_statepoint_subelement(element)
self._create_sourcepoint_subelement(element)
self._create_surf_source_read_subelement(element)
self._create_surf_source_write_subelement(element)
self._create_confidence_intervals(element)
self._create_electron_treatment_subelement(element)
self._create_energy_mode_subelement(element)
self._create_max_order_subelement(element)
self._create_photon_transport_subelement(element)
self._create_ptables_subelement(element)
self._create_seed_subelement(element)
self._create_survival_biasing_subelement(element)
self._create_cutoff_subelement(element)
self._create_entropy_mesh_subelement(element, mesh_memo)
self._create_trigger_subelement(element)
self._create_no_reduce_subelement(element)
self._create_verbosity_subelement(element)
self._create_tabular_legendre_subelements(element)
self._create_temperature_subelements(element)
self._create_trace_subelement(element)
self._create_track_subelement(element)
self._create_ufs_mesh_subelement(element, mesh_memo)
self._create_resonance_scattering_subelement(element)
self._create_volume_calcs_subelement(element)
self._create_create_fission_neutrons_subelement(element)
self._create_delayed_photon_scaling_subelement(element)
self._create_event_based_subelement(element)
self._create_max_particles_in_flight_subelement(element)
self._create_material_cell_offsets_subelement(element)
self._create_log_grid_bins_subelement(element)
self._create_write_initial_source_subelement(element)
self._create_weight_windows_subelement(element, mesh_memo)
self._create_max_splits_subelement(element)
self._create_max_tracks_subelement(element)
# Clean the indentation in the file to be user-readable
clean_indentation(element)
reorder_attributes(element) # TODO: Remove when support is Python 3.8+
return element
def export_to_xml(self, path: PathLike = 'settings.xml'):
"""Export simulation settings to an XML file.
@ -1545,57 +1655,7 @@ class Settings:
Path to file to write. Defaults to 'settings.xml'.
"""
# Reset xml element tree
root_element = ET.Element("settings")
self._create_run_mode_subelement(root_element)
self._create_particles_subelement(root_element)
self._create_batches_subelement(root_element)
self._create_inactive_subelement(root_element)
self._create_max_lost_particles_subelement(root_element)
self._create_rel_max_lost_particles_subelement(root_element)
self._create_generations_per_batch_subelement(root_element)
self._create_keff_trigger_subelement(root_element)
self._create_source_subelement(root_element)
self._create_output_subelement(root_element)
self._create_statepoint_subelement(root_element)
self._create_sourcepoint_subelement(root_element)
self._create_surf_source_read_subelement(root_element)
self._create_surf_source_write_subelement(root_element)
self._create_confidence_intervals(root_element)
self._create_electron_treatment_subelement(root_element)
self._create_energy_mode_subelement(root_element)
self._create_max_order_subelement(root_element)
self._create_photon_transport_subelement(root_element)
self._create_ptables_subelement(root_element)
self._create_seed_subelement(root_element)
self._create_survival_biasing_subelement(root_element)
self._create_cutoff_subelement(root_element)
self._create_entropy_mesh_subelement(root_element)
self._create_trigger_subelement(root_element)
self._create_no_reduce_subelement(root_element)
self._create_verbosity_subelement(root_element)
self._create_tabular_legendre_subelements(root_element)
self._create_temperature_subelements(root_element)
self._create_trace_subelement(root_element)
self._create_track_subelement(root_element)
self._create_ufs_mesh_subelement(root_element)
self._create_resonance_scattering_subelement(root_element)
self._create_volume_calcs_subelement(root_element)
self._create_create_fission_neutrons_subelement(root_element)
self._create_delayed_photon_scaling_subelement(root_element)
self._create_event_based_subelement(root_element)
self._create_max_particles_in_flight_subelement(root_element)
self._create_material_cell_offsets_subelement(root_element)
self._create_log_grid_bins_subelement(root_element)
self._create_write_initial_source_subelement(root_element)
self._create_weight_windows_subelement(root_element)
self._create_max_splits_subelement(root_element)
self._create_max_tracks_subelement(root_element)
# Clean the indentation in the file to be user-readable
clean_indentation(root_element)
root_element = self.to_xml_element()
# Check if path is a directory
p = Path(path)
@ -1603,10 +1663,78 @@ class Settings:
p /= 'settings.xml'
# Write the XML Tree to the settings.xml file
reorder_attributes(root_element) # TODO: Remove when support is Python 3.8+
tree = ET.ElementTree(root_element)
tree.write(str(p), xml_declaration=True, encoding='utf-8')
@classmethod
def from_xml_element(cls, elem, meshes=None):
"""Generate settings from XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
meshes : dict or None
A dictionary with mesh IDs as keys and mesh instances as values that
have already been read from XML. Pre-existing meshes are used
and new meshes are added to when creating tally objects.
Returns
-------
openmc.Settings
Settings object
"""
settings = cls()
settings._eigenvalue_from_xml_element(elem)
settings._run_mode_from_xml_element(elem)
settings._particles_from_xml_element(elem)
settings._batches_from_xml_element(elem)
settings._inactive_from_xml_element(elem)
settings._max_lost_particles_from_xml_element(elem)
settings._rel_max_lost_particles_from_xml_element(elem)
settings._generations_per_batch_from_xml_element(elem)
settings._keff_trigger_from_xml_element(elem)
settings._source_from_xml_element(elem)
settings._volume_calcs_from_xml_element(elem)
settings._output_from_xml_element(elem)
settings._statepoint_from_xml_element(elem)
settings._sourcepoint_from_xml_element(elem)
settings._surf_source_read_from_xml_element(elem)
settings._surf_source_write_from_xml_element(elem)
settings._confidence_intervals_from_xml_element(elem)
settings._electron_treatment_from_xml_element(elem)
settings._energy_mode_from_xml_element(elem)
settings._max_order_from_xml_element(elem)
settings._photon_transport_from_xml_element(elem)
settings._ptables_from_xml_element(elem)
settings._seed_from_xml_element(elem)
settings._survival_biasing_from_xml_element(elem)
settings._cutoff_from_xml_element(elem)
settings._entropy_mesh_from_xml_element(elem, meshes)
settings._trigger_from_xml_element(elem)
settings._no_reduce_from_xml_element(elem)
settings._verbosity_from_xml_element(elem)
settings._tabular_legendre_from_xml_element(elem)
settings._temperature_from_xml_element(elem)
settings._trace_from_xml_element(elem)
settings._track_from_xml_element(elem)
settings._ufs_mesh_from_xml_element(elem, meshes)
settings._resonance_scattering_from_xml_element(elem)
settings._create_fission_neutrons_from_xml_element(elem)
settings._delayed_photon_scaling_from_xml_element(elem)
settings._event_based_from_xml_element(elem)
settings._max_particles_in_flight_from_xml_element(elem)
settings._material_cell_offsets_from_xml_element(elem)
settings._log_grid_bins_from_xml_element(elem)
settings._write_initial_source_from_xml_element(elem)
settings._weight_windows_from_xml_element(elem, meshes)
settings._max_splits_from_xml_element(elem)
settings._max_tracks_from_xml_element(elem)
# TODO: Get volume calculations
return settings
@classmethod
def from_xml(cls, path: PathLike = 'settings.xml'):
"""Generate settings from XML file
@ -1626,54 +1754,4 @@ class Settings:
"""
tree = ET.parse(path)
root = tree.getroot()
settings = cls()
settings._eigenvalue_from_xml_element(root)
settings._run_mode_from_xml_element(root)
settings._particles_from_xml_element(root)
settings._batches_from_xml_element(root)
settings._inactive_from_xml_element(root)
settings._max_lost_particles_from_xml_element(root)
settings._rel_max_lost_particles_from_xml_element(root)
settings._generations_per_batch_from_xml_element(root)
settings._keff_trigger_from_xml_element(root)
settings._source_from_xml_element(root)
settings._volume_calcs_from_xml_element(root)
settings._output_from_xml_element(root)
settings._statepoint_from_xml_element(root)
settings._sourcepoint_from_xml_element(root)
settings._surf_source_read_from_xml_element(root)
settings._surf_source_write_from_xml_element(root)
settings._confidence_intervals_from_xml_element(root)
settings._electron_treatment_from_xml_element(root)
settings._energy_mode_from_xml_element(root)
settings._max_order_from_xml_element(root)
settings._photon_transport_from_xml_element(root)
settings._ptables_from_xml_element(root)
settings._seed_from_xml_element(root)
settings._survival_biasing_from_xml_element(root)
settings._cutoff_from_xml_element(root)
settings._entropy_mesh_from_xml_element(root)
settings._trigger_from_xml_element(root)
settings._no_reduce_from_xml_element(root)
settings._verbosity_from_xml_element(root)
settings._tabular_legendre_from_xml_element(root)
settings._temperature_from_xml_element(root)
settings._trace_from_xml_element(root)
settings._track_from_xml_element(root)
settings._ufs_mesh_from_xml_element(root)
settings._resonance_scattering_from_xml_element(root)
settings._create_fission_neutrons_from_xml_element(root)
settings._delayed_photon_scaling_from_xml_element(root)
settings._event_based_from_xml_element(root)
settings._max_particles_in_flight_from_xml_element(root)
settings._material_cell_offsets_from_xml_element(root)
settings._log_grid_bins_from_xml_element(root)
settings._write_initial_source_from_xml_element(root)
settings._weight_windows_from_xml_element(root)
settings._max_splits_from_xml_element(root)
settings._max_tracks_from_xml_element(root)
# TODO: Get volume calculations
return settings
return cls.from_xml_element(root)

45
openmc/source.py
View file

@ -2,6 +2,9 @@ from collections.abc import Iterable
from enum import Enum
from numbers import Real
import warnings
import typing # imported separately as py3.8 requires typing.Iterable
# also required to prevent typing.Union namespace overwriting Union
from typing import Optional, Sequence
from xml.etree import ElementTree as ET
import numpy as np
@ -9,6 +12,7 @@ import h5py
import openmc
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from openmc.stats.multivariate import UnitSphere, Spatial
from openmc.stats.univariate import Univariate
from ._xml import get_text
@ -70,9 +74,19 @@ class Source:
"""
def __init__(self, space=None, angle=None, energy=None, time=None, filename=None,
library=None, parameters=None, strength=1.0, particle='neutron',
domains=None):
def __init__(
self,
space: Optional[openmc.stats.Spatial] = None,
angle: Optional[openmc.stats.UnitSphere] = None,
energy: Optional[openmc.stats.Univariate] = None,
time: Optional[openmc.stats.Univariate] = None,
filename: Optional[str] = None,
library: Optional[str] = None,
parameters: Optional[str] = None,
strength: float = 1.0,
particle: str = 'neutron',
domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None
):
self._space = None
self._angle = None
self._energy = None
@ -209,7 +223,7 @@ class Source:
cv.check_value('source particle', particle, ['neutron', 'photon'])
self._particle = particle
def to_xml_element(self):
def to_xml_element(self) -> ET.Element:
"""Return XML representation of the source
Returns
@ -244,7 +258,7 @@ class Source:
return element
@classmethod
def from_xml_element(cls, elem):
def from_xml_element(cls, elem: ET.Element) -> 'openmc.Source':
"""Generate source from an XML element
Parameters
@ -349,8 +363,18 @@ class SourceParticle:
Type of the particle
"""
def __init__(self, r=(0., 0., 0.), u=(0., 0., 1.), E=1.0e6, time=0.0, wgt=1.0,
delayed_group=0, surf_id=0, particle=ParticleType.NEUTRON):
def __init__(
self,
r: typing.Iterable[float] = (0., 0., 0.),
u: typing.Iterable[float] = (0., 0., 1.),
E: float = 1.0e6,
time: float = 0.0,
wgt: float = 1.0,
delayed_group: int = 0,
surf_id: int = 0,
particle: ParticleType = ParticleType.NEUTRON
):
self.r = tuple(r)
self.u = tuple(u)
self.E = float(E)
@ -364,7 +388,7 @@ class SourceParticle:
name = self.particle.name.lower()
return f'<SourceParticle: {name} at E={self.E:.6e} eV>'
def to_tuple(self):
def to_tuple(self) -> tuple:
"""Return source particle attributes as a tuple
Returns
@ -377,7 +401,10 @@ class SourceParticle:
self.delayed_group, self.surf_id, self.particle.value)
def write_source_file(source_particles, filename, **kwargs):
def write_source_file(
source_particles: typing.Iterable[SourceParticle],
filename: PathLike, **kwargs
):
"""Write a source file using a collection of source particles
Parameters

2
openmc/stats/univariate.py
View file

@ -729,7 +729,7 @@ def muir(e0, m_rat, kt):
distribution: the mean energy of particles ``e0``, the mass of reactants
``m_rat``, and the ion temperature ``kt``.
.. versionadded:: 0.14.0
.. versionadded:: 0.13.2
Parameters
----------

119
openmc/tallies.py
View file

@ -3119,17 +3119,17 @@ class Tallies(cv.CheckedList):
for tally in self:
root_element.append(tally.to_xml_element())
def _create_mesh_subelements(self, root_element):
already_written = set()
def _create_mesh_subelements(self, root_element, memo=None):
already_written = memo if memo else set()
for tally in self:
for f in tally.filters:
if isinstance(f, openmc.MeshFilter):
if f.mesh.id not in already_written:
if len(f.mesh.name) > 0:
root_element.append(ET.Comment(f.mesh.name))
root_element.append(f.mesh.to_xml_element())
already_written.add(f.mesh.id)
if f.mesh.id in already_written:
continue
if len(f.mesh.name) > 0:
root_element.append(ET.Comment(f.mesh.name))
root_element.append(f.mesh.to_xml_element())
already_written.add(f.mesh.id)
def _create_filter_subelements(self, root_element):
already_written = dict()
@ -3155,6 +3155,22 @@ class Tallies(cv.CheckedList):
for d in derivs:
root_element.append(d.to_xml_element())
def to_xml_element(self, memo=None):
"""Creates a 'tallies' element to be written to an XML file.
"""
element = ET.Element("tallies")
self._create_mesh_subelements(element, memo)
self._create_filter_subelements(element)
self._create_tally_subelements(element)
self._create_derivative_subelements(element)
# Clean the indentation in the file to be user-readable
clean_indentation(element)
reorder_attributes(element) # TODO: Remove when support is Python 3.8+
return element
def export_to_xml(self, path='tallies.xml'):
"""Create a tallies.xml file that can be used for a simulation.
@ -3164,15 +3180,7 @@ class Tallies(cv.CheckedList):
Path to file to write. Defaults to 'tallies.xml'.
"""
root_element = ET.Element("tallies")
self._create_mesh_subelements(root_element)
self._create_filter_subelements(root_element)
self._create_tally_subelements(root_element)
self._create_derivative_subelements(root_element)
# Clean the indentation in the file to be user-readable
clean_indentation(root_element)
root_element = self.to_xml_element()
# Check if path is a directory
p = Path(path)
@ -3180,10 +3188,56 @@ class Tallies(cv.CheckedList):
p /= 'tallies.xml'
# Write the XML Tree to the tallies.xml file
reorder_attributes(root_element) # TODO: Remove when support is Python 3.8+
tree = ET.ElementTree(root_element)
tree.write(str(p), xml_declaration=True, encoding='utf-8')
@classmethod
def from_xml_element(cls, elem, meshes=None):
"""Generate tallies from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
meshes : dict or None
A dictionary with mesh IDs as keys and mesh instances as values that
have already been read from XML. Pre-existing meshes are used
and new meshes are added to when creating tally objects.
Returns
-------
openmc.Tallies
Tallies object
"""
# Read mesh elements
meshes = {} if meshes is None else meshes
for e in elem.findall('mesh'):
mesh = MeshBase.from_xml_element(e)
meshes[mesh.id] = mesh
# Read filter elements
filters = {}
for e in elem.findall('filter'):
filter = openmc.Filter.from_xml_element(e, meshes=meshes)
filters[filter.id] = filter
# Read derivative elements
derivatives = {}
for e in elem.findall('derivative'):
deriv = openmc.TallyDerivative.from_xml_element(e)
derivatives[deriv.id] = deriv
# Read tally elements
tallies = []
for e in elem.findall('tally'):
tally = openmc.Tally.from_xml_element(
e, filters=filters, derivatives=derivatives
)
tallies.append(tally)
return cls(tallies)
@classmethod
def from_xml(cls, path='tallies.xml'):
"""Generate tallies from XML file
@ -3201,31 +3255,4 @@ class Tallies(cv.CheckedList):
"""
tree = ET.parse(path)
root = tree.getroot()
# Read mesh elements
meshes = {}
for elem in root.findall('mesh'):
mesh = MeshBase.from_xml_element(elem)
meshes[mesh.id] = mesh
# Read filter elements
filters = {}
for elem in root.findall('filter'):
filter = openmc.Filter.from_xml_element(elem, meshes=meshes)
filters[filter.id] = filter
# Read derivative elements
derivatives = {}
for elem in root.findall('derivative'):
deriv = openmc.TallyDerivative.from_xml_element(elem)
derivatives[deriv.id] = deriv
# Read tally elements
tallies = []
for elem in root.findall('tally'):
tally = openmc.Tally.from_xml_element(
elem, filters=filters, derivatives=derivatives
)
tallies.append(tally)
return cls(tallies)
return cls.from_xml_element(root)

13
openmc/universe.py
View file

@ -361,7 +361,7 @@ class Universe(UniverseBase):
img_path = Path(tmpdir) / f'plot_{plot.id}.png'
if not img_path.is_file():
img_path = img_path.with_suffix('.ppm')
img = mpimg.imread(img_path)
img = mpimg.imread(str(img_path))
# Create a figure sized such that the size of the axes within
# exactly matches the number of pixels specified
@ -469,7 +469,7 @@ class Universe(UniverseBase):
"""
nuclides = OrderedDict()
if self._atoms is not None:
if self._atoms:
volume = self.volume
for name, atoms in self._atoms.items():
nuclide = openmc.Nuclide(name)
@ -650,19 +650,26 @@ class DAGMCUniverse(UniverseBase):
bounding_box : 2-tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
of the universe.
.. versionadded:: 0.13.1
material_names : list of str
Return a sorted list of materials names that are contained within the
DAGMC h5m file. This is useful when naming openmc.Material() objects
as each material name present in the DAGMC h5m file must have a
matching openmc.Material() with the same name.
.. versionadded:: 0.13.2
n_cells : int
The number of cells in the DAGMC model. This is the number of cells at
runtime and accounts for the implicit complement whether or not is it
present in the DAGMC file.
.. versionadded:: 0.13.2
n_surfaces : int
The number of surfaces in the model.
.. versionadded:: 0.13.1
.. versionadded:: 0.13.2
"""
def __init__(self,

41
openmc/weight_windows.py
View file

@ -159,6 +159,35 @@ class WeightWindows(IDManagerMixin):
string += '{: <16}=\t{}\n'.format('\tWeight Cutoff', self._weight_cutoff)
return string
def __eq__(self, other):
# ensure that `other` is a WeightWindows object
if not isinstance(other, WeightWindows):
return False
# TODO: add ability to check mesh equality
# check several attributes directly
attrs = ('particle_type',
'survival_ratio',
'max_lower_bound_ratio',
'max_split',
'weight_cutoff')
for attr in attrs:
if getattr(self, attr) != getattr(other, attr):
return False
# save most expensive checks for last
if not np.array_equal(self.energy_bounds, other.energy_bounds):
return False
if not np.array_equal(self.lower_ww_bounds, other.lower_ww_bounds):
return False
if not np.array_equal(self.upper_ww_bounds, other.upper_ww_bounds):
return False
return True
@property
def mesh(self):
return self._mesh
@ -206,9 +235,9 @@ class WeightWindows(IDManagerMixin):
# reshape data according to mesh and energy bins
bounds = np.asarray(bounds)
if isinstance(self.mesh, UnstructuredMesh):
bounds.reshape(-1, self.num_energy_bins)
bounds = bounds.reshape(-1, self.num_energy_bins)
else:
bounds.reshape(*self.mesh.dimension, self.num_energy_bins)
bounds = bounds.reshape(*self.mesh.dimension, self.num_energy_bins)
self._lower_ww_bounds = bounds
@property
@ -225,9 +254,9 @@ class WeightWindows(IDManagerMixin):
# reshape data according to mesh and energy bins
bounds = np.asarray(bounds)
if isinstance(self.mesh, UnstructuredMesh):
bounds.reshape(-1, self.num_energy_bins)
bounds = bounds.reshape(-1, self.num_energy_bins)
else:
bounds.reshape(*self.mesh.dimension, self.num_energy_bins)
bounds = bounds.reshape(*self.mesh.dimension, self.num_energy_bins)
self._upper_ww_bounds = bounds
@property
@ -341,6 +370,10 @@ class WeightWindows(IDManagerMixin):
particle_type = get_text(elem, 'particle_type')
survival_ratio = float(get_text(elem, 'survival_ratio'))
ww_shape = (len(e_bounds) - 1,) + mesh.dimension[::-1]
lower_ww_bounds = np.array(lower_ww_bounds).reshape(ww_shape).T
upper_ww_bounds = np.array(upper_ww_bounds).reshape(ww_shape).T
max_lower_bound_ratio = None
if get_text(elem, 'max_lower_bound_ratio'):
max_lower_bound_ratio = float(get_text(elem, 'max_lower_bound_ratio'))

5
scripts/openmc-update-inputs
View file

@ -4,7 +4,6 @@
"""
import argparse
from difflib import get_close_matches
from itertools import chain
from random import randint
from shutil import move
@ -27,8 +26,8 @@ geometry.xml: Lattices containing 'outside' attributes/tags will be replaced
will be renamed 'region'.
materials.xml: Nuclide names will be changed from ACE aliases (e.g., Am-242m) to
HDF5/GND names (e.g., Am242_m1). Thermal scattering table names will be
changed from ACE aliases (e.g., HH2O) to HDF5/GND names (e.g., c_H_in_H2O).
HDF5/GNDS names (e.g., Am242_m1). Thermal scattering table names will be
changed from ACE aliases (e.g., HH2O) to HDF5/GNDS names (e.g., c_H_in_H2O).
"""

2
setup.py
View file

@ -28,7 +28,7 @@ kwargs = {
# Data files and libraries
'package_data': {
'openmc.lib': ['libopenmc.{}'.format(suffix)],
'openmc.data': ['mass16.txt', 'BREMX.DAT', 'half_life.json', '*.h5'],
'openmc.data': ['mass_1.mas20.txt', 'BREMX.DAT', 'half_life.json', '*.h5'],
'openmc.data.effective_dose': ['*.txt']
},

899
src/cell.cpp
View file

@ -36,241 +36,6 @@ vector<unique_ptr<Cell>> cells;
} // namespace model
//==============================================================================
//! Convert region specification string to integer tokens.
//!
//! The characters (, ), |, and ~ count as separate tokens since they represent
//! operators.
//==============================================================================
vector<int32_t> tokenize(const std::string region_spec)
{
// Check for an empty region_spec first.
vector<int32_t> tokens;
if (region_spec.empty()) {
return tokens;
}
// Parse all halfspaces and operators except for intersection (whitespace).
for (int i = 0; i < region_spec.size();) {
if (region_spec[i] == '(') {
tokens.push_back(OP_LEFT_PAREN);
i++;
} else if (region_spec[i] == ')') {
tokens.push_back(OP_RIGHT_PAREN);
i++;
} else if (region_spec[i] == '|') {
tokens.push_back(OP_UNION);
i++;
} else if (region_spec[i] == '~') {
tokens.push_back(OP_COMPLEMENT);
i++;
} else if (region_spec[i] == '-' || region_spec[i] == '+' ||
std::isdigit(region_spec[i])) {
// This is the start of a halfspace specification. Iterate j until we
// find the end, then push-back everything between i and j.
int j = i + 1;
while (j < region_spec.size() && std::isdigit(region_spec[j])) {
j++;
}
tokens.push_back(std::stoi(region_spec.substr(i, j - i)));
i = j;
} else if (std::isspace(region_spec[i])) {
i++;
} else {
auto err_msg =
fmt::format("Region specification contains invalid character, \"{}\"",
region_spec[i]);
fatal_error(err_msg);
}
}
// Add in intersection operators where a missing operator is needed.
int i = 0;
while (i < tokens.size() - 1) {
bool left_compat {(tokens[i] < OP_UNION) || (tokens[i] == OP_RIGHT_PAREN)};
bool right_compat {(tokens[i + 1] < OP_UNION) ||
(tokens[i + 1] == OP_LEFT_PAREN) ||
(tokens[i + 1] == OP_COMPLEMENT)};
if (left_compat && right_compat) {
tokens.insert(tokens.begin() + i + 1, OP_INTERSECTION);
}
i++;
}
return tokens;
}
//==============================================================================
//! Add precedence for infix regions so intersections have higher
//! precedence than unions using parentheses.
//==============================================================================
std::vector<int32_t>::iterator add_parentheses(
std::vector<int32_t>::iterator start, std::vector<int32_t>& infix)
{
int32_t start_token = *start;
// Add left parenthesis
if (start_token == OP_INTERSECTION) {
start = infix.insert(start - 1, OP_LEFT_PAREN);
} else {
start = infix.insert(start + 1, OP_LEFT_PAREN);
}
start++;
// Initialize return iterator
auto return_iterator = infix.begin();
// Add right parenthesis
// While the start iterator is within the bounds of infix
while (start < infix.end()) {
start++;
// If the current token is an operator and is different than the start token
if (*start >= OP_UNION && *start != start_token) {
// Skip wrapped regions but save iterator position to check precedence and
// add right parenthesis, right parenthesis position depends on the
// operator, when the operator is a union then do not include the operator
// in the region, when the operator is an intersection then include the
// operato and next surface
if (*start == OP_LEFT_PAREN) {
return_iterator = start;
int depth = 1;
do {
start++;
if (*start > OP_COMPLEMENT) {
if (*start == OP_RIGHT_PAREN) {
depth--;
} else {
depth++;
}
}
} while (depth > 0);
} else {
start = infix.insert(
start_token == OP_UNION ? start - 1 : start, OP_RIGHT_PAREN);
if (return_iterator == infix.begin()) {
return_iterator = start - 1;
}
return return_iterator;
}
}
}
// If we get here a right parenthesis hasn't been placed,
// return iterator
infix.push_back(OP_RIGHT_PAREN);
if (return_iterator == infix.begin()) {
return_iterator = start - 1;
}
return return_iterator;
}
void add_precedence(std::vector<int32_t>& infix)
{
int32_t current_op = 0;
for (auto it = infix.begin(); it != infix.end(); it++) {
int32_t token = *it;
if (token == OP_UNION || token == OP_INTERSECTION) {
if (current_op == 0) {
// Set the current operator if is hasn't been set
current_op = token;
} else if (token != current_op) {
// If the current operator doesn't match the token, add parenthesis to assert precedence
it = add_parentheses(it, infix);
current_op = 0;
}
} else if (token > OP_COMPLEMENT) {
// If the token is a parenthesis reset the current operator
current_op = 0;
}
}
}
//==============================================================================
//! Convert infix region specification to Reverse Polish Notation (RPN)
//!
//! This function uses the shunting-yard algorithm.
//==============================================================================
vector<int32_t> generate_postfix(int32_t cell_id, vector<int32_t> infix)
{
vector<int32_t> rpn;
vector<int32_t> stack;
for (int32_t token : infix) {
if (token < OP_UNION) {
// If token is not an operator, add it to output
rpn.push_back(token);
} else if (token < OP_RIGHT_PAREN) {
// Regular operators union, intersection, complement
while (stack.size() > 0) {
int32_t op = stack.back();
if (op < OP_RIGHT_PAREN && ((token == OP_COMPLEMENT && token < op) ||
(token != OP_COMPLEMENT && token <= op))) {
// While there is an operator, op, on top of the stack, if the token
// is left-associative and its precedence is less than or equal to
// that of op or if the token is right-associative and its precedence
// is less than that of op, move op to the output queue and push the
// token on to the stack. Note that only complement is
// right-associative.
rpn.push_back(op);
stack.pop_back();
} else {
break;
}
}
stack.push_back(token);
} else if (token == OP_LEFT_PAREN) {
// If the token is a left parenthesis, push it onto the stack
stack.push_back(token);
} else {
// If the token is a right parenthesis, move operators from the stack to
// the output queue until reaching the left parenthesis.
for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) {
// If we run out of operators without finding a left parenthesis, it
// means there are mismatched parentheses.
if (it == stack.rend()) {
fatal_error(fmt::format(
"Mismatched parentheses in region specification for cell {}",
cell_id));
}
rpn.push_back(stack.back());
stack.pop_back();
}
// Pop the left parenthesis.
stack.pop_back();
}
}
while (stack.size() > 0) {
int32_t op = stack.back();
// If the operator is a parenthesis it is mismatched.
if (op >= OP_RIGHT_PAREN) {
fatal_error(fmt::format(
"Mismatched parentheses in region specification for cell {}", cell_id));
}
rpn.push_back(stack.back());
stack.pop_back();
}
return rpn;
}
//==============================================================================
// Cell implementation
//==============================================================================
@ -334,12 +99,16 @@ double Cell::temperature(int32_t instance) const
void Cell::set_temperature(double T, int32_t instance, bool set_contained)
{
if (settings::temperature_method == TemperatureMethod::INTERPOLATION) {
if (T < data::temperature_min) {
throw std::runtime_error {"Temperature is below minimum temperature at "
"which data is available."};
} else if (T > data::temperature_max) {
throw std::runtime_error {"Temperature is above maximum temperature at "
"which data is available."};
if (T < (data::temperature_min - settings::temperature_tolerance)) {
throw std::runtime_error {
fmt::format("Temperature of {} K is below minimum temperature at "
"which data is available of {} K.",
T, data::temperature_min)};
} else if (T > (data::temperature_max + settings::temperature_tolerance)) {
throw std::runtime_error {
fmt::format("Temperature of {} K is above maximum temperature at "
"which data is available of {} K.",
T, data::temperature_max)};
}
}
@ -586,44 +355,8 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
// Get a tokenized representation of the region specification and apply De
// Morgans law
region_ = tokenize(region_spec);
remove_complement_ops(region_);
// Convert user IDs to surface indices.
for (auto& r : region_) {
if (r < OP_UNION) {
const auto& it {model::surface_map.find(abs(r))};
if (it == model::surface_map.end()) {
throw std::runtime_error {
"Invalid surface ID " + std::to_string(abs(r)) +
" specified in region for cell " + std::to_string(id_) + "."};
}
r = (r > 0) ? it->second + 1 : -(it->second + 1);
}
}
// Check if this is a simple cell.
simple_ = true;
for (int32_t token : region_) {
if (token == OP_UNION) {
simple_ = false;
// Ensure intersections have precedence over unions
add_precedence(region_);
break;
}
}
region_.shrink_to_fit();
// If this cell is simple, remove all the superfluous operator tokens.
if (simple_) {
for (auto it = region_.begin(); it != region_.end(); it++) {
if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) {
region_.erase(it);
it--;
}
}
region_.shrink_to_fit();
}
Region region(region_spec, id_);
region_ = region;
// Read the translation vector.
if (check_for_node(cell_node, "translation")) {
@ -650,98 +383,13 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
//==============================================================================
bool CSGCell::contains(Position r, Direction u, int32_t on_surface) const
{
if (simple_) {
return contains_simple(r, u, on_surface);
} else {
return contains_complex(r, u, on_surface);
}
}
//==============================================================================
std::pair<double, int32_t> CSGCell::distance(
Position r, Direction u, int32_t on_surface, Particle* p) const
{
double min_dist {INFTY};
int32_t i_surf {std::numeric_limits<int32_t>::max()};
for (int32_t token : region_) {
// Ignore this token if it corresponds to an operator rather than a region.
if (token >= OP_UNION)
continue;
// Calculate the distance to this surface.
// Note the off-by-one indexing
bool coincident {std::abs(token) == std::abs(on_surface)};
double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
// Check if this distance is the new minimum.
if (d < min_dist) {
if (min_dist - d >= FP_PRECISION * min_dist) {
min_dist = d;
i_surf = -token;
}
}
}
return {min_dist, i_surf};
}
//==============================================================================
void CSGCell::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "geom_type", "csg", false);
// Write the region specification.
if (!region_.empty()) {
std::stringstream region_spec {};
for (int32_t token : region_) {
if (token == OP_LEFT_PAREN) {
region_spec << " (";
} else if (token == OP_RIGHT_PAREN) {
region_spec << " )";
} else if (token == OP_COMPLEMENT) {
region_spec << " ~";
} else if (token == OP_INTERSECTION) {
} else if (token == OP_UNION) {
region_spec << " |";
} else {
// Note the off-by-one indexing
auto surf_id = model::surfaces[abs(token) - 1]->id_;
region_spec << " " << ((token > 0) ? surf_id : -surf_id);
}
}
write_string(group_id, "region", region_spec.str(), false);
}
write_string(group_id, "region", region_.str(), false);
}
BoundingBox CSGCell::bounding_box_simple() const
{
BoundingBox bbox;
for (int32_t token : region_) {
bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0);
}
return bbox;
}
void CSGCell::apply_demorgan(
vector<int32_t>::iterator start, vector<int32_t>::iterator stop)
{
do {
if (*start < OP_UNION) {
*start *= -1;
} else if (*start == OP_UNION) {
*start = OP_INTERSECTION;
} else if (*start == OP_INTERSECTION) {
*start = OP_UNION;
}
start++;
} while (start < stop);
}
//==============================================================================
vector<int32_t>::iterator CSGCell::find_left_parenthesis(
vector<int32_t>::iterator start, const vector<int32_t>& infix)
@ -776,75 +424,409 @@ vector<int32_t>::iterator CSGCell::find_left_parenthesis(
return it;
}
void CSGCell::remove_complement_ops(vector<int32_t>& infix)
{
auto it = std::find(infix.begin(), infix.end(), OP_COMPLEMENT);
while (it != infix.end()) {
// Erase complement
infix.erase(it);
//==============================================================================
// Region implementation
//==============================================================================
// Define stop given left parenthesis or not
auto stop = it;
if (*it == OP_LEFT_PAREN) {
int depth = 1;
do {
stop++;
if (*stop > OP_COMPLEMENT) {
if (*stop == OP_RIGHT_PAREN) {
depth--;
} else {
depth++;
}
Region::Region(std::string region_spec, int32_t cell_id)
{
// Check if region_spec is not empty.
if (!region_spec.empty()) {
// Parse all halfspaces and operators except for intersection (whitespace).
for (int i = 0; i < region_spec.size();) {
if (region_spec[i] == '(') {
expression_.push_back(OP_LEFT_PAREN);
i++;
} else if (region_spec[i] == ')') {
expression_.push_back(OP_RIGHT_PAREN);
i++;
} else if (region_spec[i] == '|') {
expression_.push_back(OP_UNION);
i++;
} else if (region_spec[i] == '~') {
expression_.push_back(OP_COMPLEMENT);
i++;
} else if (region_spec[i] == '-' || region_spec[i] == '+' ||
std::isdigit(region_spec[i])) {
// This is the start of a halfspace specification. Iterate j until we
// find the end, then push-back everything between i and j.
int j = i + 1;
while (j < region_spec.size() && std::isdigit(region_spec[j])) {
j++;
}
} while (depth > 0);
it++;
expression_.push_back(std::stoi(region_spec.substr(i, j - i)));
i = j;
} else if (std::isspace(region_spec[i])) {
i++;
} else {
auto err_msg =
fmt::format("Region specification contains invalid character, \"{}\"",
region_spec[i]);
fatal_error(err_msg);
}
}
// apply DeMorgan's law to any surfaces/operators between these
// positions in the RPN
apply_demorgan(it, stop);
// update iterator position
it = std::find(infix.begin(), infix.end(), OP_COMPLEMENT);
}
}
BoundingBox CSGCell::bounding_box_complex(vector<int32_t> postfix)
{
vector<BoundingBox> stack(postfix.size());
int i_stack = -1;
for (auto& token : postfix) {
if (token == OP_UNION) {
stack[i_stack - 1] = stack[i_stack - 1] | stack[i_stack];
i_stack--;
} else if (token == OP_INTERSECTION) {
stack[i_stack - 1] = stack[i_stack - 1] & stack[i_stack];
i_stack--;
} else {
i_stack++;
stack[i_stack] = model::surfaces[abs(token) - 1]->bounding_box(token > 0);
// Add in intersection operators where a missing operator is needed.
int i = 0;
while (i < expression_.size() - 1) {
bool left_compat {
(expression_[i] < OP_UNION) || (expression_[i] == OP_RIGHT_PAREN)};
bool right_compat {(expression_[i + 1] < OP_UNION) ||
(expression_[i + 1] == OP_LEFT_PAREN) ||
(expression_[i + 1] == OP_COMPLEMENT)};
if (left_compat && right_compat) {
expression_.insert(expression_.begin() + i + 1, OP_INTERSECTION);
}
i++;
}
}
Ensures(i_stack == 0);
return stack.front();
}
// Remove complement operators using DeMorgan's laws
auto it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT);
while (it != expression_.end()) {
// Erase complement
expression_.erase(it);
// Define stop given left parenthesis or not
auto stop = it;
if (*it == OP_LEFT_PAREN) {
int depth = 1;
do {
stop++;
if (*stop > OP_COMPLEMENT) {
if (*stop == OP_RIGHT_PAREN) {
depth--;
} else {
depth++;
}
}
} while (depth > 0);
it++;
}
// apply DeMorgan's law to any surfaces/operators between these
// positions in the RPN
apply_demorgan(it, stop);
// update iterator position
it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT);
}
// Convert user IDs to surface indices.
for (auto& r : expression_) {
if (r < OP_UNION) {
const auto& it {model::surface_map.find(abs(r))};
if (it == model::surface_map.end()) {
throw std::runtime_error {
"Invalid surface ID " + std::to_string(abs(r)) +
" specified in region for cell " + std::to_string(cell_id) + "."};
}
r = (r > 0) ? it->second + 1 : -(it->second + 1);
}
}
// Check if this is a simple cell.
simple_ = true;
for (int32_t token : expression_) {
if (token == OP_UNION) {
simple_ = false;
// Ensure intersections have precedence over unions
add_precedence();
break;
}
}
// If this cell is simple, remove all the superfluous operator tokens.
if (simple_) {
for (auto it = expression_.begin(); it != expression_.end(); it++) {
if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) {
expression_.erase(it);
it--;
}
}
}
expression_.shrink_to_fit();
BoundingBox CSGCell::bounding_box() const
{
if (simple_) {
return bounding_box_simple();
} else {
auto postfix = generate_postfix(this->id_, this->region_);
return bounding_box_complex(postfix);
simple_ = true;
}
}
//==============================================================================
bool CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const
void Region::apply_demorgan(
vector<int32_t>::iterator start, vector<int32_t>::iterator stop)
{
for (int32_t token : region_) {
do {
if (*start < OP_UNION) {
*start *= -1;
} else if (*start == OP_UNION) {
*start = OP_INTERSECTION;
} else if (*start == OP_INTERSECTION) {
*start = OP_UNION;
}
start++;
} while (start < stop);
}
//==============================================================================
//! Add precedence for infix regions so intersections have higher
//! precedence than unions using parentheses.
//==============================================================================
std::vector<int32_t>::iterator Region::add_parentheses(
std::vector<int32_t>::iterator start)
{
int32_t start_token = *start;
// Add left parenthesis
if (start_token == OP_INTERSECTION) {
start = expression_.insert(start - 1, OP_LEFT_PAREN);
} else {
start = expression_.insert(start + 1, OP_LEFT_PAREN);
}
start++;
// Keep track of return iterator distance. If we don't encounter a left
// parenthesis, we return an iterator corresponding to wherever the right
// parenthesis is inserted. If a left parenthesis is encountered, an iterator
// corresponding to the left parenthesis is returned. Also note that we keep
// track of a *distance* instead of an iterator because the underlying memory
// allocation may change.
std::size_t return_it_dist = 0;
// Add right parenthesis
// While the start iterator is within the bounds of infix
while (start < expression_.end()) {
start++;
// If the current token is an operator and is different than the start token
if (*start >= OP_UNION && *start != start_token) {
// Skip wrapped regions but save iterator position to check precedence and
// add right parenthesis, right parenthesis position depends on the
// operator, when the operator is a union then do not include the operator
// in the region, when the operator is an intersection then include the
// operator and next surface
if (*start == OP_LEFT_PAREN) {
return_it_dist = std::distance(expression_.begin(), start);
int depth = 1;
do {
start++;
if (*start > OP_COMPLEMENT) {
if (*start == OP_RIGHT_PAREN) {
depth--;
} else {
depth++;
}
}
} while (depth > 0);
} else {
start = expression_.insert(
start_token == OP_UNION ? start - 1 : start, OP_RIGHT_PAREN);
if (return_it_dist > 0) {
return expression_.begin() + return_it_dist;
} else {
return start - 1;
}
}
}
}
// If we get here a right parenthesis hasn't been placed,
// return iterator
expression_.push_back(OP_RIGHT_PAREN);
if (return_it_dist > 0) {
return expression_.begin() + return_it_dist;
} else {
return start - 1;
}
}
//==============================================================================
void Region::add_precedence()
{
int32_t current_op = 0;
std::size_t current_dist = 0;
for (auto it = expression_.begin(); it != expression_.end(); it++) {
int32_t token = *it;
if (token == OP_UNION || token == OP_INTERSECTION) {
if (current_op == 0) {
// Set the current operator if is hasn't been set
current_op = token;
current_dist = std::distance(expression_.begin(), it);
} else if (token != current_op) {
// If the current operator doesn't match the token, add parenthesis to
// assert precedence
if (current_op == OP_INTERSECTION) {
it = add_parentheses(expression_.begin() + current_dist);
} else {
it = add_parentheses(it);
}
current_op = 0;
current_dist = 0;
}
} else if (token > OP_COMPLEMENT) {
// If the token is a parenthesis reset the current operator
current_op = 0;
current_dist = 0;
}
}
}
//==============================================================================
//! Convert infix region specification to Reverse Polish Notation (RPN)
//!
//! This function uses the shunting-yard algorithm.
//==============================================================================
vector<int32_t> Region::generate_postfix(int32_t cell_id) const
{
vector<int32_t> rpn;
vector<int32_t> stack;
for (int32_t token : expression_) {
if (token < OP_UNION) {
// If token is not an operator, add it to output
rpn.push_back(token);
} else if (token < OP_RIGHT_PAREN) {
// Regular operators union, intersection, complement
while (stack.size() > 0) {
int32_t op = stack.back();
if (op < OP_RIGHT_PAREN && ((token == OP_COMPLEMENT && token < op) ||
(token != OP_COMPLEMENT && token <= op))) {
// While there is an operator, op, on top of the stack, if the token
// is left-associative and its precedence is less than or equal to
// that of op or if the token is right-associative and its precedence
// is less than that of op, move op to the output queue and push the
// token on to the stack. Note that only complement is
// right-associative.
rpn.push_back(op);
stack.pop_back();
} else {
break;
}
}
stack.push_back(token);
} else if (token == OP_LEFT_PAREN) {
// If the token is a left parenthesis, push it onto the stack
stack.push_back(token);
} else {
// If the token is a right parenthesis, move operators from the stack to
// the output queue until reaching the left parenthesis.
for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) {
// If we run out of operators without finding a left parenthesis, it
// means there are mismatched parentheses.
if (it == stack.rend()) {
fatal_error(fmt::format(
"Mismatched parentheses in region specification for cell {}",
cell_id));
}
rpn.push_back(stack.back());
stack.pop_back();
}
// Pop the left parenthesis.
stack.pop_back();
}
}
while (stack.size() > 0) {
int32_t op = stack.back();
// If the operator is a parenthesis it is mismatched.
if (op >= OP_RIGHT_PAREN) {
fatal_error(fmt::format(
"Mismatched parentheses in region specification for cell {}", cell_id));
}
rpn.push_back(stack.back());
stack.pop_back();
}
return rpn;
}
//==============================================================================
std::string Region::str() const
{
std::stringstream region_spec {};
if (!expression_.empty()) {
for (int32_t token : expression_) {
if (token == OP_LEFT_PAREN) {
region_spec << " (";
} else if (token == OP_RIGHT_PAREN) {
region_spec << " )";
} else if (token == OP_COMPLEMENT) {
region_spec << " ~";
} else if (token == OP_INTERSECTION) {
} else if (token == OP_UNION) {
region_spec << " |";
} else {
// Note the off-by-one indexing
auto surf_id = model::surfaces[abs(token) - 1]->id_;
region_spec << " " << ((token > 0) ? surf_id : -surf_id);
}
}
}
return region_spec.str();
}
//==============================================================================
std::pair<double, int32_t> Region::distance(
Position r, Direction u, int32_t on_surface, Particle* p) const
{
double min_dist {INFTY};
int32_t i_surf {std::numeric_limits<int32_t>::max()};
for (int32_t token : expression_) {
// Ignore this token if it corresponds to an operator rather than a region.
if (token >= OP_UNION)
continue;
// Calculate the distance to this surface.
// Note the off-by-one indexing
bool coincident {std::abs(token) == std::abs(on_surface)};
double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
// Check if this distance is the new minimum.
if (d < min_dist) {
if (min_dist - d >= FP_PRECISION * min_dist) {
min_dist = d;
i_surf = -token;
}
}
}
return {min_dist, i_surf};
}
//==============================================================================
bool Region::contains(Position r, Direction u, int32_t on_surface) const
{
if (simple_) {
return contains_simple(r, u, on_surface);
} else {
return contains_complex(r, u, on_surface);
}
}
//==============================================================================
bool Region::contains_simple(Position r, Direction u, int32_t on_surface) const
{
for (int32_t token : expression_) {
// Assume that no tokens are operators. Evaluate the sense of particle with
// respect to the surface and see if the token matches the sense. If the
// particle's surface attribute is set and matches the token, that
@ -865,14 +847,13 @@ bool CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const
//==============================================================================
bool CSGCell::contains_complex(
Position r, Direction u, int32_t on_surface) const
bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const
{
bool in_cell = true;
int total_depth = 0;
// For each token
for (auto it = region_.begin(); it != region_.end(); it++) {
for (auto it = expression_.begin(); it != expression_.end(); it++) {
int32_t token = *it;
// If the token is a surface evaluate the sense
@ -923,6 +904,76 @@ bool CSGCell::contains_complex(
return in_cell;
}
//==============================================================================
BoundingBox Region::bounding_box(int32_t cell_id) const
{
if (simple_) {
return bounding_box_simple();
} else {
auto postfix = generate_postfix(cell_id);
return bounding_box_complex(postfix);
}
}
//==============================================================================
BoundingBox Region::bounding_box_simple() const
{
BoundingBox bbox;
for (int32_t token : expression_) {
bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0);
}
return bbox;
}
//==============================================================================
BoundingBox Region::bounding_box_complex(vector<int32_t> postfix) const
{
vector<BoundingBox> stack(postfix.size());
int i_stack = -1;
for (auto& token : postfix) {
if (token == OP_UNION) {
stack[i_stack - 1] = stack[i_stack - 1] | stack[i_stack];
i_stack--;
} else if (token == OP_INTERSECTION) {
stack[i_stack - 1] = stack[i_stack - 1] & stack[i_stack];
i_stack--;
} else {
i_stack++;
stack[i_stack] = model::surfaces[abs(token) - 1]->bounding_box(token > 0);
}
}
Ensures(i_stack == 0);
return stack.front();
}
//==============================================================================
vector<int32_t> Region::surfaces() const
{
if (simple_) {
return expression_;
}
vector<int32_t> surfaces = expression_;
auto it = std::find_if(surfaces.begin(), surfaces.end(),
[&](const auto& value) { return value >= OP_UNION; });
while (it != surfaces.end()) {
surfaces.erase(it);
it = std::find_if(surfaces.begin(), surfaces.end(),
[&](const auto& value) { return value >= OP_UNION; });
}
return surfaces;
}
//==============================================================================
// Non-method functions
//==============================================================================

5
src/cross_sections.cpp
View file

@ -104,6 +104,11 @@ void read_cross_sections_xml()
auto root = doc.document_element();
read_cross_sections_xml(root);
}
void read_cross_sections_xml(pugi::xml_node root)
{
// Find cross_sections.xml file -- the first place to look is the
// materials.xml file. If no file is found there, then we check the
// OPENMC_CROSS_SECTIONS environment variable

13
src/dagmc.cpp
View file

@ -486,13 +486,21 @@ void DAGUniverse::legacy_assign_material(
// if no material was set using a name, assign by id
if (!mat_found_by_name) {
bool found_by_id = true;
try {
auto id = std::stoi(mat_string);
if (model::material_map.find(id) == model::material_map.end())
found_by_id = false;
c->material_.emplace_back(id);
} catch (const std::invalid_argument&) {
fatal_error(fmt::format(
"No material '{}' found for volume (cell) {}", mat_string, c->id_));
found_by_id = false;
}
// report failure for failed int conversion or missing material
if (!found_by_id)
fatal_error(
fmt::format("Material with name/ID '{}' not found for volume (cell) {}",
mat_string, c->id_));
}
if (settings::verbosity >= 10) {
@ -557,7 +565,6 @@ DAGCell::DAGCell(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx)
: Cell {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx)
{
geom_type_ = GeometryType::DAG;
simple_ = true;
};
std::pair<double, int32_t> DAGCell::distance(

Some files were not shown because too many files have changed in this diff Show more