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

Merge branch 'develop' into distribution-alias-sampling

Browse source
This commit is contained in:
Patrick Myers 2023-01-13 10:33:19 -05:00 • committed by GitHub
commit 2eae8b8f8b
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
94 changed files with 2977 additions and 470 deletions
Download patch file Download diff file Expand all files Collapse all files

8
.github/workflows/ci.yml vendored
View file

@ -29,6 +29,7 @@ jobs:
mpi: [n, y]
omp: [n, y]
dagmc: [n]
ncrystal: [n]
libmesh: [n]
event: [n]
vectfit: [n]
@ -47,6 +48,10 @@ jobs:
python-version: '3.10'
mpi: y
omp: y
- ncrystal: y
python-version: '3.10'
mpi: n
omp: n
- libmesh: y
python-version: '3.10'
mpi: y
@ -64,7 +69,7 @@ jobs:
omp: n
mpi: y
name: "Python ${{ matrix.python-version }} (omp=${{ matrix.omp }},
mpi=${{ matrix.mpi }}, dagmc=${{ matrix.dagmc }},
mpi=${{ matrix.mpi }}, dagmc=${{ matrix.dagmc }}, ncrystal=${{ matrix.ncrystal }},
libmesh=${{ matrix.libmesh }}, event=${{ matrix.event }}
vectfit=${{ matrix.vectfit }})"
@ -73,6 +78,7 @@ jobs:
PHDF5: ${{ matrix.mpi }}
OMP: ${{ matrix.omp }}
DAGMC: ${{ matrix.dagmc }}
NCRYSTAL: ${{ matrix.ncrystal }}
EVENT: ${{ matrix.event }}
VECTFIT: ${{ matrix.vectfit }}
LIBMESH: ${{ matrix.libmesh }}

32
CMakeLists.txt
View file

@ -36,6 +36,8 @@ 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
foreach(OLD_OPT IN ITEMS "openmp" "profile" "coverage" "dagmc" "libmesh")
@ -98,6 +100,15 @@ macro(find_package_write_status pkg)
endif()
endmacro()
#===============================================================================
# NCrystal Scattering Support
#===============================================================================
if(OPENMC_USE_NCRYSTAL)
find_package(NCrystal REQUIRED)
message(STATUS "Found NCrystal: ${NCrystal_DIR} (version ${NCrystal_VERSION})")
endif()
#===============================================================================
# DAGMC Geometry Support - need DAGMC/MOAB
#===============================================================================
@ -158,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
#===============================================================================
@ -331,10 +351,12 @@ 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
src/mgxs_interface.cpp
src/ncrystal_interface.cpp
src/nuclide.cpp
src/output.cpp
src/particle.cpp
@ -494,6 +516,16 @@ endif()
include(CTest)
add_subdirectory(tests/cpp_unit_tests)
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
#===============================================================================

9
cmake/OpenMCConfig.cmake.in
View file

@ -9,6 +9,11 @@ if(@OPENMC_USE_DAGMC@)
find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@)
endif()
if(@OPENMC_USE_NCRYSTAL@)
find_package(NCrystal REQUIRED)
message(STATUS "Found NCrystal: ${NCrystal_DIR} (version ${NCrystal_VERSION})")
endif()
if(@OPENMC_USE_LIBMESH@)
include(FindPkgConfig)
list(APPEND CMAKE_PREFIX_PATH @LIBMESH_PREFIX@)
@ -25,3 +30,7 @@ endif()
if(@OPENMC_USE_MPI@)
find_package(MPI REQUIRED)
endif()
if(@OPENMC_USE_MCPL@)
find_package(MCPL REQUIRED)
endif()

4
docs/source/conf.py
View file

@ -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)
}

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

@ -32,6 +32,17 @@ standard deviation.
*Default*: false
-------------------------------------
``<create_delayed_neutrons>`` Element
-------------------------------------
The ``<create_delayed_neutrons>`` element indicates whether delayed neutrons
are created in fission. If this element is set to "true", delayed neutrons
will be created in fission events; otherwise only prompt neutrons will be
created.
*Default*: true
-------------------------------------
``<create_fission_neutrons>`` Element
-------------------------------------
@ -732,6 +743,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 +786,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
------------------------------

28
docs/source/methods/cross_sections.rst
View file

@ -178,6 +178,27 @@ been selected. There are three methods available:
section data is loaded for a single temperature and is used in the
unresolved resonance and fast energy ranges.
------------------
NCrystal materials
------------------
As an alternative of the standard thermal scattering treatment using
:math:`S(\alpha,\beta)` tables, OpenMC allows to create materials using
NCrystal_. In addition to the regular thermal elastic, and thermal inelastic
processes, NCrystal allows the generation of models for materials that cannot
currently included in ACE files such as oriented single crystals (see the
`NCrystal paper`_), and further extend the physics `using plugins`_. Thermal
scattering kernels are generated on the fly from dynamic and structural data, or
loaded from :math:`S(\alpha,\beta)` tables converted from ENDF6 evaluations.
These kernels are sampled in a direct way using a fast `rejection algorithm`_
that does not require previous processing. A `large library`_ of materials is
already included in the NCrystal distribution, and new materials can be easily
defined from scratch in the `NCMAT format`_ or `combining existing files`_.
The compositions of the materials defined in NCrystal are passed on to OpenMC
all other reactions except for thermal neutron scattering are handled by
continuous energy ACE libraries.
----------------
Multi-Group Data
----------------
@ -279,3 +300,10 @@ or even isotropic scattering.
.. _ENDF/B data: https://www.nndc.bnl.gov/endf-b8.0/
.. _Leppanen: https://doi.org/10.1016/j.anucene.2009.03.019
.. _algorithms: http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package
.. _NCrystal: https://github.com/mctools/ncrystal
.. _NCrystal paper: https://doi.org/10.1016/j.cpc.2019.07.015
.. _using plugins: https://doi.org/10.1016/j.cpc.2021.108082
.. _rejection algorithm: https://doi.org/10.1016/j.jcp.2018.11.043
.. _large library: https://github.com/mctools/ncrystal/wiki/Data-library
.. _NCMAT format: https://github.com/mctools/ncrystal/wiki/NCMAT-format
.. _combining existing files: https://github.com/mctools/ncrystal/wiki/Announcement-Release3.0.0#2-multiphase-materials

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

@ -66,7 +66,7 @@ Core Functions
decay_energy
decay_photon_energy
dose_coefficients
gnd_name
gnds_name
half_life
isotopes
kalbach_slope

18
docs/source/usersguide/install.rst
View file

@ -271,6 +271,17 @@ Prerequisites
cmake -DOPENMC_USE_DAGMC=on -DCMAKE_PREFIX_PATH=/path/to/dagmc/installation ..
* NCrystal_ library for defining materials with enhanced thermal neutron transport
Adding this option allows the creation of materials from NCrystal, which
replaces the scattering kernel treatment of ACE files with a modular,
on-the-fly approach. To use it `install
<https://github.com/mctools/ncrystal/wiki/Get-NCrystal>`_ and `initialize
<https://github.com/mctools/ncrystal/wiki/Using-NCrystal#setting-up>`_
NCrystal and turn on the option in the CMake configuration step:
cmake -DOPENMC_USE_NCRYSTAL=on ..
* libMesh_ mesh library framework for numerical simulations of partial differential equations
This optional dependency enables support for unstructured mesh tally
@ -294,6 +305,7 @@ Prerequisites
.. _MOAB: https://bitbucket.org/fathomteam/moab
.. _libMesh: https://libmesh.github.io/
.. _libpng: http://www.libpng.org/pub/png/libpng.html
.. _NCrystal: https://github.com/mctools/ncrystal
Obtaining the Source
--------------------
@ -362,6 +374,12 @@ OPENMC_USE_DAGMC
should also be defined as `DAGMC_ROOT` in the CMake configuration command.
(Default: off)
OPENMC_USE_NCRYSTAL
Turns on support for NCrystal materials. NCrystal must be
`installed <https://github.com/mctools/ncrystal/wiki/Get-NCrystal>`_ and
`initialized <https://github.com/mctools/ncrystal/wiki/Using-NCrystal#setting-up>`_.
(Default: off)
OPENMC_USE_LIBMESH
Enables the use of unstructured mesh tallies with libMesh_. (Default: off)

37
docs/source/usersguide/materials.rst
View file

@ -101,6 +101,42 @@ you would need to add hydrogen and oxygen to a material and then assign the
.. _usersguide_naming:
-------------------------
Adding NCrystal materials
-------------------------
Additional support for thermal scattering can be added by using NCrystal_. The
:meth:`Material.from_ncrystal` class method generates a :class:`openmc.Material`
object from an `NCrystal configuration string
<https://github.com/mctools/ncrystal/wiki/Using-NCrystal#uniform-material-configuration-syntax>`_.
Temperature, material composition, and density are passed from the configuration
string and the `NCMAT file
<https://github.com/mctools/ncrystal/wiki/NCMAT-format>`_ that define the
material, e.g.::
mat = openmc.Material.from_ncrystal('Al_sg225.ncmat;temp=300K')
defines a material containing polycrystalline alumnium,
::
mat = openmc.Material.from_ncrystal("""Ge_sg227.ncmat;dcutoff=0.5;mos=40arcsec;
dir1=@crys_hkl:5,1,1@lab:0,0,1;
dir2=@crys_hkl:0,-1,1@lab:0,1,0""")
defines an oriented germanium single crystal with 40 arcsec mosaicity.
NCrystal only handles low energy neutron interactions. Other interactions are
provided by standard ACE files. NCrystal_ comes with a `predefined library
<https://github.com/mctools/ncrystal/wiki/Data-library>`_ but more materials can
be added by creating NCMAT files or on-the-fly in the configuration string.
.. warning:: Currently, NCrystal_ materials cannot be modified after they are created.
Density, temperature and composition should be defined in the
configuration string or the NCMAT file.
.. _NCrystal: https://github.com/mctools/ncrystal
------------------
Naming Conventions
------------------
@ -222,3 +258,4 @@ been generated, you can tell OpenMC to use this file either by setting
materials.cross_sections = '/path/to/cross_sections.xml'
.. _MCNP: https://mcnp.lanl.gov/

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
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

10
include/openmc/material.h
View file

@ -12,6 +12,7 @@
#include "openmc/bremsstrahlung.h"
#include "openmc/constants.h"
#include "openmc/memory.h" // for unique_ptr
#include "openmc/ncrystal_interface.h"
#include "openmc/particle.h"
#include "openmc/vector.h"
@ -151,12 +152,17 @@ public:
//! \return Temperature in [K]
double temperature() const;
//! Get pointer to NCrystal material object
//! \return Pointer to NCrystal material object
const NCrystalMat& ncrystal_mat() const { return ncrystal_mat_; };
//----------------------------------------------------------------------------
// Data
int32_t id_ {C_NONE}; //!< Unique ID
std::string name_; //!< Name of material
vector<int> nuclide_; //!< Indices in nuclides vector
vector<int> element_; //!< Indices in elements vector
NCrystalMat ncrystal_mat_; //!< NCrystal material object
xt::xtensor<double, 1> atom_density_; //!< Nuclide atom density in [atom/b-cm]
double density_; //!< Total atom density in [atom/b-cm]
double density_gpcc_; //!< Total atom density in [g/cm^3]
@ -221,6 +227,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

94
include/openmc/ncrystal_interface.h Normal file
View file

@ -0,0 +1,94 @@
#ifndef OPENMC_NCRYSTAL_INTERFACE_H
#define OPENMC_NCRYSTAL_INTERFACE_H
#ifdef NCRYSTAL
#include "NCrystal/NCRNG.hh"
#include "NCrystal/NCrystal.hh"
#endif
#include "openmc/particle.h"
#include <cstdint> // for uint64_t
#include <limits> // for numeric_limits
#include <string>
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
extern "C" const bool NCRYSTAL_ENABLED;
//! Energy in [eV] to switch between NCrystal and ENDF
constexpr double NCRYSTAL_MAX_ENERGY {5.0};
//==============================================================================
// Wrapper class an NCrystal material
//==============================================================================
class NCrystalMat {
public:
//----------------------------------------------------------------------------
// Constructors
NCrystalMat() = default;
explicit NCrystalMat(const std::string& cfg);
//----------------------------------------------------------------------------
// Methods
#ifdef NCRYSTAL
//! Return configuration string
std::string cfg() const;
//! Get cross section from NCrystal material
//
//! \param[in] p Particle object
//! \return Cross section in [b]
double xs(const Particle& p) const;
// Process scattering event
//
//! \param[in] p Particle object
void scatter(Particle& p) const;
//! Whether the object holds a valid NCrystal material
operator bool() const;
#else
//----------------------------------------------------------------------------
// Trivial methods when compiling without NCRYSTAL
std::string cfg() const
{
return "";
}
double xs(const Particle& p) const
{
return -1.0;
}
void scatter(Particle& p) const {}
operator bool() const
{
return false;
}
#endif
private:
//----------------------------------------------------------------------------
// Data members (only present when compiling with NCrystal support)
#ifdef NCRYSTAL
std::string cfg_; //!< NCrystal configuration string
std::shared_ptr<const NCrystal::ProcImpl::Process>
ptr_; //!< Pointer to NCrystal material object
#endif
};
//==============================================================================
// Functions
//==============================================================================
void ncrystal_update_micro(double xs, NuclideMicroXS& micro);
} // namespace openmc
#endif // OPENMC_NCRYSTAL_INTERFACE_H

5
include/openmc/physics.h
View file

@ -96,8 +96,9 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p);
void sample_secondary_photons(Particle& p, int i_nuclide);
//！Split or Roulette particles based their weight and the lower weight window
// bound.
//! Split or Roulette particles based their weight and the lower weight window
//! bound.
//
//! \param[in] p, particle to be split or rouletted with the weight window.
void split_particle(Particle& p);

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();

8
include/openmc/settings.h
View file

@ -28,6 +28,7 @@ extern bool check_overlaps; //!< check overlaps in geometry?
extern bool confidence_intervals; //!< use confidence intervals for results?
extern bool
create_fission_neutrons; //!< create fission neutrons (fixed source)?
extern bool create_delayed_neutrons; //!< create delayed fission neutrons?
extern "C" bool cmfd_run; //!< is a CMFD run?
extern bool
delayed_photon_scaling; //!< Scale fission photon yield to include delayed
@ -47,7 +48,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?
@ -127,9 +130,12 @@ extern double weight_survive; //!< Survival weight after Russian roulette
//==============================================================================
//! 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;

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();

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

13
openmc/cell.py
View file

@ -1,6 +1,5 @@
from collections import OrderedDict
from collections.abc import Iterable
from copy import deepcopy
from math import cos, sin, pi
from numbers import Real
from xml.etree import ElementTree as ET
@ -519,8 +518,14 @@ class Cell(IDManagerMixin):
paths = self._paths
self._paths = None
clone = deepcopy(self)
clone.id = None
clone = openmc.Cell(name=self.name)
clone.volume = self.volume
if self.temperature is not None:
clone.temperature = self.temperature
if self.translation is not None:
clone.translation = self.translation
if self.rotation is not None:
clone.rotation = self.rotation
clone._num_instances = None
# Restore paths on original instance
@ -650,7 +655,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/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)

19
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 = {}
@ -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")

2
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']

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

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

4
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

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)

10
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.

16
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:
@ -904,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)

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
----------

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)

157
openmc/geometry.py
View file

@ -1,3 +1,5 @@
import os
import typing
from collections import OrderedDict, defaultdict
from collections.abc import Iterable
from copy import deepcopy
@ -7,7 +9,7 @@ import warnings
import openmc
import openmc._xml as xml
from .checkvalue import check_type, check_less_than, check_greater_than
from .checkvalue import check_type, check_less_than, check_greater_than, PathLike
class Geometry:
@ -103,6 +105,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.
@ -117,25 +152,7 @@ class Geometry:
.. versionadded:: 0.12
"""
# 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
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)
@ -143,18 +160,17 @@ 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(cls, path='geometry.xml', materials=None):
"""Generate geometry from XML file
def from_xml_element(cls, elem, materials=None):
"""Generate geometry from an XML element
Parameters
----------
path : str, optional
Path to geometry XML file
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.
@ -165,6 +181,11 @@ class 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):
@ -173,13 +194,10 @@ class Geometry:
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'):
for surface in elem.findall('surface'):
s = openmc.Surface.from_xml_element(surface)
surfaces[s.id] = s
@ -193,24 +211,24 @@ class Geometry:
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)
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 root)
# contain it (needed to determine which universe is the elem)
child_of = defaultdict(list)
for elem in root.findall('lattice'):
lat = openmc.RectLattice.from_xml_element(elem, get_universe)
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 elem in root.findall('hex_lattice'):
lat = openmc.HexLattice.from_xml_element(elem, get_universe)
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)
@ -224,15 +242,8 @@ class Geometry:
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)
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)
@ -244,6 +255,41 @@ class Geometry:
else:
raise ValueError('Error determining root universe.')
@classmethod
def from_xml(
cls,
path: PathLike = 'geometry.xml',
materials: typing.Optional[typing.Union[PathLike, 'openmc.Materials']] = 'materials.xml'
):
"""Generate geometry from XML file
Parameters
----------
path : PathLike, optional
Path to geometry XML file
materials : openmc.Materials or PathLike
Materials used to assign to cells. If PathLike, an attempt is made
to generate materials from the provided xml file.
Returns
-------
openmc.Geometry
Geometry object
"""
# Using str and os.Pathlike here to avoid error when using just the imported PathLike
# TypeError: Subscripted generics cannot be used with class and instance checks
check_type('materials', materials, (str, os.PathLike, openmc.Materials))
if isinstance(materials, (str, os.PathLike)):
materials = openmc.Materials.from_xml(materials)
tree = ET.parse(path)
root = tree.getroot()
return cls.from_xml_element(root, materials)
def find(self, point):
"""Find cells/universes/lattices which contain a given point
@ -488,6 +534,27 @@ class Geometry:
"""
return self._get_domains_by_name(name, case_sensitive, matching, 'cell')
def get_surfaces_by_name(self, name, case_sensitive=False, matching=False):
"""Return a list of surfaces with matching names.
Parameters
----------
name : str
The name to search match
case_sensitive : bool
Whether to distinguish upper and lower case letters in each
surface's name (default is False)
matching : bool
Whether the names must match completely (default is False)
Returns
-------
list of openmc.Surface
Surfaces matching the queried name
"""
return self._get_domains_by_name(name, case_sensitive, matching, 'surface')
def get_cells_by_fill_name(self, name, case_sensitive=False, matching=False):
"""Return a list of cells with fills with matching names.

6
openmc/lib/__init__.py
View file

@ -42,12 +42,18 @@ else:
def _dagmc_enabled():
return c_bool.in_dll(_dll, "DAGMC_ENABLED").value
def _ncrystal_enabled():
return c_bool.in_dll(_dll, "NCRYSTAL_ENABLED").value
def _coord_levels():
return c_int.in_dll(_dll, "n_coord_levels").value
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 *

203
openmc/material.py
View file

@ -99,6 +99,10 @@ class Material(IDManagerMixin):
[decay/sec].
.. versionadded:: 0.13.2
ncrystal_cfg : str
NCrystal configuration string
.. versionadded:: 0.13.3
"""
@ -118,6 +122,7 @@ class Material(IDManagerMixin):
self._volume = None
self._atoms = {}
self._isotropic = []
self._ncrystal_cfg = None
# A list of tuples (nuclide, percent, percent type)
self._nuclides = []
@ -140,6 +145,9 @@ class Material(IDManagerMixin):
string += '{: <16}\n'.format('\tS(a,b) Tables')
if self._ncrystal_cfg:
string += '{: <16}=\t{}\n'.format('\tNCrystal conf', self._ncrystal_cfg)
for sab in self._sab:
string += '{: <16}=\t{}\n'.format('\tS(a,b)', sab)
@ -219,6 +227,10 @@ class Material(IDManagerMixin):
def volume(self):
return self._volume
@property
def ncrystal_cfg(self):
return self._ncrystal_cfg
@name.setter
def name(self, name: Optional[str]):
if name is not None:
@ -331,6 +343,64 @@ class Material(IDManagerMixin):
return material
@classmethod
def from_ncrystal(cls, cfg, **kwargs):
"""Create material from NCrystal configuration string.
Density, temperature, and material composition, and (ultimately) thermal
neutron scattering will be automatically be provided by NCrystal based
on this string. The name and material_id parameters are simply passed on
to the Material constructor.
Parameters
----------
cfg : str
NCrystal configuration string
**kwargs
Keyword arguments passed to :class:`openmc.Material`
Returns
-------
openmc.Material
Material instance
"""
import NCrystal
nc_mat = NCrystal.createInfo(cfg)
def openmc_natabund(Z):
#nc_mat.getFlattenedComposition might need natural abundancies.
#This call-back function is used so NCrystal can flatten composition
#using OpenMC's natural abundancies. In practice this function will
#only get invoked in the unlikely case where a material is specified
#by referring both to natural elements and specific isotopes of the
#same element.
elem_name = openmc.data.ATOMIC_SYMBOL[Z]
return [
(int(iso_name[len(elem_name):]), abund)
for iso_name, abund in openmc.data.isotopes(elem_name)
]
flat_compos = nc_mat.getFlattenedComposition(
preferNaturalElements=True, naturalAbundProvider=openmc_natabund)
# Create the Material
material = cls(temperature=nc_mat.getTemperature(), **kwargs)
for Z, A_vals in flat_compos:
elemname = openmc.data.ATOMIC_SYMBOL[Z]
for A, frac in A_vals:
if A:
material.add_nuclide(f'{elemname}{A}', frac)
else:
material.add_element(elemname, frac)
material.set_density('g/cm3', nc_mat.getDensity())
material._ncrystal_cfg = NCrystal.normaliseCfg(cfg)
return material
def add_volume_information(self, volume_calc):
"""Add volume information to a material.
@ -405,6 +475,9 @@ class Material(IDManagerMixin):
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if self._ncrystal_cfg is not None:
raise ValueError("Cannot add nuclides to NCrystal material")
# If nuclide name doesn't look valid, give a warning
try:
Z, _, _ = openmc.data.zam(nuclide)
@ -609,6 +682,9 @@ class Material(IDManagerMixin):
raise ValueError("Element name should be given by the "
"element's symbol or name, e.g., 'Zr', 'zirconium'")
if self._ncrystal_cfg is not None:
raise ValueError("Cannot add elements to NCrystal material")
# Allow for element identifier to be given as a symbol or name
if len(element) > 2:
el = element.lower()
@ -990,7 +1066,7 @@ 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].
@ -1024,16 +1100,16 @@ class Material(IDManagerMixin):
multiplier = 1
elif units == 'W/g':
multiplier = 1.0 / self.get_mass_density()
decayheat = {}
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())
return decayheat if by_nuclide else sum(decayheat.values())
def get_nuclide_atoms(self):
"""Return number of atoms of each nuclide in the material
@ -1182,6 +1258,14 @@ class Material(IDManagerMixin):
if self._volume:
element.set("volume", str(self._volume))
if self._ncrystal_cfg:
if self._sab:
raise ValueError("NCrystal materials are not compatible with S(a,b).")
if self._macroscopic is not None:
raise ValueError("NCrystal materials are not compatible with macroscopic cross sections.")
element.set("cfg", str(self._ncrystal_cfg))
# Create temperature XML subelement
if self.temperature is not None:
element.set("temperature", str(self.temperature))
@ -1449,6 +1533,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.
@ -1468,32 +1603,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'):
@ -1513,14 +1650,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)

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)']

114
openmc/model/model.py
View file

@ -6,10 +6,12 @@ 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
@ -238,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
@ -399,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
----------
@ -418,14 +449,7 @@ class Model:
d.mkdir(parents=True)
self.settings.export_to_xml(d)
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()
self.geometry.export_to_xml(d)
self.geometry.export_to_xml(d, remove_surfs=remove_surfs)
# If a materials collection was specified, export it. Otherwise, look
# for all materials in the geometry and use that to automatically build
@ -442,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

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

2
openmc/plotter.py
View file

@ -563,7 +563,7 @@ def _calculate_cexs_elem_mat(this, types, temperature=294.,
for nuclide in nuclides.items():
sabs[nuclide[0]] = None
if isinstance(this, openmc.Material):
for sab_name in this._sab:
for sab_name, _ in this._sab:
sab = openmc.data.ThermalScattering.from_hdf5(
library.get_by_material(sab_name, data_type='thermal')['path'])
for nuc in sab.nuclides:

374
openmc/settings.py
View file

@ -157,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:
@ -172,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
@ -223,6 +225,10 @@ class Settings:
Weight windows to use for variance reduction
.. versionadded:: 0.13
create_delayed_neutrons : bool
Whether delayed neutrons are created in fission.
.. versionadded:: 0.13.3
weight_windows_on : bool
Whether weight windows are enabled
@ -294,6 +300,7 @@ class Settings:
VolumeCalculation, 'volume calculations')
self._create_fission_neutrons = None
self._create_delayed_neutrons = None
self._delayed_photon_scaling = None
self._material_cell_offsets = None
self._log_grid_bins = None
@ -457,6 +464,10 @@ class Settings:
def create_fission_neutrons(self) -> bool:
return self._create_fission_neutrons
@property
def create_delayed_neutrons(self) -> bool:
return self._create_delayed_neutrons
@property
def delayed_photon_scaling(self) -> bool:
return self._delayed_photon_scaling
@ -622,6 +633,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.")
@ -655,7 +668,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)
@ -667,6 +680,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
@ -859,6 +875,12 @@ class Settings:
create_fission_neutrons, bool)
self._create_fission_neutrons = create_fission_neutrons
@create_delayed_neutrons.setter
def create_delayed_neutrons(self, create_delayed_neutrons: bool):
cv.check_type('Whether create only prompt neutrons',
create_delayed_neutrons, bool)
self._create_delayed_neutrons = create_delayed_neutrons
@delayed_photon_scaling.setter
def delayed_photon_scaling(self, value: bool):
cv.check_type('delayed photon scaling', value, bool)
@ -1018,6 +1040,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")
@ -1035,6 +1061,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:
@ -1073,25 +1102,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:
@ -1142,15 +1181,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
@ -1177,6 +1222,11 @@ class Settings:
elem = ET.SubElement(root, "create_fission_neutrons")
elem.text = str(self._create_fission_neutrons).lower()
def _create_create_delayed_neutrons_subelement(self, root):
if self._create_delayed_neutrons is not None:
elem = ET.SubElement(root, "create_delayed_neutrons")
elem.text = str(self._create_delayed_neutrons).lower()
def _create_delayed_photon_scaling_subelement(self, root):
if self._delayed_photon_scaling is not None:
elem = ET.SubElement(root, "delayed_photon_scaling")
@ -1207,15 +1257,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")
@ -1315,10 +1371,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')
@ -1337,13 +1393,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):
@ -1396,13 +1454,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')
@ -1462,13 +1522,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')
@ -1490,6 +1552,11 @@ class Settings:
if text is not None:
self.create_fission_neutrons = text in ('true', '1')
def _create_delayed_neutrons_from_xml_element(self, root):
text = get_text(root, 'create_delayed_neutrons')
if text is not None:
self.create_delayed_neutrons = text in ('true', '1')
def _delayed_photon_scaling_from_xml_element(self, root):
text = get_text(root, 'delayed_photon_scaling')
if text is not None:
@ -1520,7 +1587,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)
@ -1529,6 +1596,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:
@ -1539,6 +1609,69 @@ 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_create_delayed_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.
@ -1548,57 +1681,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)
@ -1606,10 +1689,79 @@ 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._create_delayed_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
@ -1629,54 +1781,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)

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)

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).
"""

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

1
src/finalize.cpp
View file

@ -74,6 +74,7 @@ int openmc_finalize()
settings::check_overlaps = false;
settings::confidence_intervals = false;
settings::create_fission_neutrons = true;
settings::create_delayed_neutrons = true;
settings::electron_treatment = ElectronTreatment::LED;
settings::delayed_photon_scaling = true;
settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0};

5
src/geometry_aux.cpp
View file

@ -61,6 +61,11 @@ void read_geometry_xml()
// Get root element
pugi::xml_node root = doc.document_element();
read_geometry_xml(root);
}
void read_geometry_xml(pugi::xml_node root)
{
// Read surfaces, cells, lattice
read_surfaces(root);
read_cells(root);

124
src/initialize.cpp
View file

@ -14,6 +14,7 @@
#include "openmc/constants.h"
#include "openmc/cross_sections.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/geometry_aux.h"
#include "openmc/hdf5_interface.h"
#include "openmc/material.h"
@ -105,7 +106,11 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
openmc::openmc_set_seed(DEFAULT_SEED);
// Read XML input files
read_input_xml();
if (!read_model_xml())
read_separate_xml_files();
// Write some initial output under the header if needed
initial_output();
// Check for particle restart run
if (settings::particle_restart_run)
@ -177,10 +182,8 @@ int parse_command_line(int argc, char* argv[])
} else if (arg == "-e" || arg == "--event") {
settings::event_based = true;
} else if (arg == "-r" || arg == "--restart") {
i += 1;
// Check what type of file this is
hid_t file_id = file_open(argv[i], 'r', true);
std::string filetype;
@ -280,7 +283,15 @@ int parse_command_line(int argc, char* argv[])
if (argc > 1 && last_flag < argc - 1) {
settings::path_input = std::string(argv[last_flag + 1]);
// Add slash at end of directory if it isn't there
// check that the path is either a valid directory or file
if (!dir_exists(settings::path_input) &&
!file_exists(settings::path_input)) {
fatal_error(fmt::format(
"The path specified to the OpenMC executable '{}' does not exist.",
settings::path_input));
}
// Add slash at end of directory if it isn't the
if (!ends_with(settings::path_input, "/")) {
settings::path_input += "/";
}
@ -289,7 +300,105 @@ int parse_command_line(int argc, char* argv[])
return 0;
}
void read_input_xml()
bool read_model_xml()
{
std::string model_filename =
settings::path_input.empty() ? "." : settings::path_input;
// some string cleanup
// a trailing "/" is applied to path_input if it's specified,
// remove it for the first attempt at reading the input file
if (ends_with(model_filename, "/"))
model_filename.pop_back();
// if the current filename is a directory, append the default model filename
if (dir_exists(model_filename))
model_filename += "/model.xml";
// if this file doesn't exist, stop here
if (!file_exists(model_filename))
return false;
// try to process the path input as an XML file
pugi::xml_document doc;
if (!doc.load_file(model_filename.c_str())) {
fatal_error(fmt::format(
"Error reading from single XML input file '{}'", model_filename));
}
pugi::xml_node root = doc.document_element();
// Read settings
if (!check_for_node(root, "settings")) {
fatal_error("No <settings> node present in the model.xml file.");
}
auto settings_root = root.child("settings");
// Verbosity
if (check_for_node(settings_root, "verbosity")) {
settings::verbosity = std::stoi(get_node_value(settings_root, "verbosity"));
}
// To this point, we haven't displayed any output since we didn't know what
// the verbosity is. Now that we checked for it, show the title if necessary
if (mpi::master) {
if (settings::verbosity >= 2)
title();
}
write_message(
fmt::format("Reading model XML file '{}' ...", model_filename), 5);
read_settings_xml(settings_root);
// If other XML files are present, display warning
// that they will be ignored
auto other_inputs = {"materials.xml", "geometry.xml", "settings.xml",
"tallies.xml", "plots.xml"};
for (const auto& input : other_inputs) {
if (file_exists(settings::path_input + input)) {
warning((fmt::format("Other XML file input(s) are present. These files "
"will be ignored in favor of the {} file.",
model_filename)));
break;
}
}
// Read materials and cross sections
if (!check_for_node(root, "materials")) {
fatal_error(fmt::format(
"No <materials> node present in the {} file.", model_filename));
}
read_cross_sections_xml(root.child("materials"));
read_materials_xml(root.child("materials"));
// Read geometry
if (!check_for_node(root, "geometry")) {
fatal_error(fmt::format(
"No <geometry> node present in the {} file.", model_filename));
}
read_geometry_xml(root.child("geometry"));
// Final geometry setup and assign temperatures
finalize_geometry();
// Finalize cross sections having assigned temperatures
finalize_cross_sections();
if (check_for_node(root, "tallies"))
read_tallies_xml(root.child("tallies"));
// Initialize distribcell_filters
prepare_distribcell();
if (check_for_node(root, "plots"))
read_plots_xml(root.child("plots"));
return true;
}
void read_separate_xml_files()
{
read_settings_xml();
read_cross_sections_xml();
@ -310,6 +419,11 @@ void read_input_xml()
// Read the plots.xml regardless of plot mode in case plots are requested
// via the API
read_plots_xml();
}
void initial_output()
{
// write initial output
if (settings::run_mode == RunMode::PLOTTING) {
// Read plots.xml if it exists
if (mpi::master && settings::verbosity >= 5)

31
src/material.cpp
View file

@ -60,6 +60,13 @@ Material::Material(pugi::xml_node node)
name_ = get_node_value(node, "name");
}
if (check_for_node(node, "cfg")) {
auto cfg = get_node_value(node, "cfg");
write_message(
5, "NCrystal config string for material #{}: '{}'", this->id(), cfg);
ncrystal_mat_ = NCrystalMat(cfg);
}
if (check_for_node(node, "depletable")) {
depletable_ = get_node_value_bool(node, "depletable");
}
@ -367,8 +374,8 @@ void Material::finalize()
this->init_thermal();
}
// Normalize density
this->normalize_density();
// Normalize density
this->normalize_density();
}
void Material::normalize_density()
@ -792,6 +799,12 @@ void Material::calculate_neutron_xs(Particle& p) const
// Initialize position in i_sab_nuclides
int j = 0;
// Calculate NCrystal cross section
double ncrystal_xs = -1.0;
if (ncrystal_mat_ && p.E() < NCRYSTAL_MAX_ENERGY) {
ncrystal_xs = ncrystal_mat_.xs(p);
}
// Add contribution from each nuclide in material
for (int i = 0; i < nuclide_.size(); ++i) {
// ======================================================================
@ -830,10 +843,16 @@ void Material::calculate_neutron_xs(Particle& p) const
int i_nuclide = nuclide_[i];
// Calculate microscopic cross section for this nuclide
const auto& micro {p.neutron_xs(i_nuclide)};
auto& micro {p.neutron_xs(i_nuclide)};
if (p.E() != micro.last_E || p.sqrtkT() != micro.last_sqrtkT ||
i_sab != micro.index_sab || sab_frac != micro.sab_frac) {
data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p);
// If NCrystal is being used, update micro cross section cache
if (ncrystal_xs >= 0.0) {
data::nuclides[i_nuclide]->calculate_elastic_xs(p);
ncrystal_update_micro(ncrystal_xs, micro);
}
}
// ======================================================================
@ -1281,6 +1300,12 @@ void read_materials_xml()
// Loop over XML material elements and populate the array.
pugi::xml_node root = doc.document_element();
read_materials_xml(root);
}
void read_materials_xml(pugi::xml_node root)
{
for (pugi::xml_node material_node : root.children("material")) {
model::materials.push_back(make_unique<Material>(material_node));
}

247
src/mcpl_interface.cpp Normal file
View file

@ -0,0 +1,247 @@
#include "openmc/mcpl_interface.h"
#include "openmc/bank.h"
#include "openmc/error.h"
#include "openmc/message_passing.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/state_point.h"
#include "openmc/vector.h"
#include <fmt/core.h>
#ifdef OPENMC_MCPL
#include <mcpl.h>
#endif
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
#ifdef OPENMC_MCPL
const bool MCPL_ENABLED = true;
#else
const bool MCPL_ENABLED = false;
#endif
//==============================================================================
// Functions
//==============================================================================
#ifdef OPENMC_MCPL
SourceSite mcpl_particle_to_site(const mcpl_particle_t* particle)
{
SourceSite site;
switch (particle->pdgcode) {
case 2112:
site.particle = ParticleType::neutron;
break;
case 22:
site.particle = ParticleType::photon;
break;
case 11:
site.particle = ParticleType::electron;
break;
case -11:
site.particle = ParticleType::positron;
break;
}
// Copy position and direction
site.r.x = particle->position[0];
site.r.y = particle->position[1];
site.r.z = particle->position[2];
site.u.x = particle->direction[0];
site.u.y = particle->direction[1];
site.u.z = particle->direction[2];
// MCPL stores kinetic energy in [MeV], time in [ms]
site.E = particle->ekin * 1e6;
site.time = particle->time * 1e-3;
site.wgt = particle->weight;
return site;
}
#endif
//==============================================================================
vector<SourceSite> mcpl_source_sites(std::string path)
{
vector<SourceSite> sites;
#ifdef OPENMC_MCPL
// Open MCPL file and determine number of particles
auto mcpl_file = mcpl_open_file(path.c_str());
size_t n_sites = mcpl_hdr_nparticles(mcpl_file);
for (int i = 0; i < n_sites; i++) {
// Extract particle from mcpl-file, checking if it is a neutron, photon,
// electron, or positron. Otherwise skip.
const mcpl_particle_t* particle;
int pdg = 0;
while (pdg != 2112 && pdg != 22 && pdg != 11 && pdg != -11) {
particle = mcpl_read(mcpl_file);
pdg = particle->pdgcode;
}
// Convert to source site and add to vector
sites.push_back(mcpl_particle_to_site(particle));
}
// Check that some sites were read
if (sites.empty()) {
fatal_error("MCPL file contained no neutron, photon, electron, or positron "
"source particles.");
}
mcpl_close_file(mcpl_file);
#else
fatal_error(
"Your build of OpenMC does not support reading MCPL source files.");
#endif
return sites;
}
//==============================================================================
#ifdef OPENMC_MCPL
void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank)
{
int64_t dims_size = settings::n_particles;
int64_t count_size = simulation::work_per_rank;
// Set vectors for source bank and starting bank index of each process
vector<int64_t>* bank_index = &simulation::work_index;
vector<SourceSite>* source_bank = &simulation::source_bank;
vector<int64_t> surf_source_index_vector;
vector<SourceSite> surf_source_bank_vector;
if (surf_source_bank) {
surf_source_index_vector = calculate_surf_source_size();
dims_size = surf_source_index_vector[mpi::n_procs];
count_size = simulation::surf_source_bank.size();
bank_index = &surf_source_index_vector;
// Copy data in a SharedArray into a vector.
surf_source_bank_vector.resize(count_size);
surf_source_bank_vector.assign(simulation::surf_source_bank.data(),
simulation::surf_source_bank.data() + count_size);
source_bank = &surf_source_bank_vector;
}
if (mpi::master) {
// Particles are writeen to disk from the master node only
// Save source bank sites since the array is overwritten below
#ifdef OPENMC_MPI
vector<SourceSite> temp_source {source_bank->begin(), source_bank->end()};
#endif
// loop over the other nodes and receive data - then write those.
for (int i = 0; i < mpi::n_procs; ++i) {
// number of particles for node node i
size_t count[] {
static_cast<size_t>((*bank_index)[i + 1] - (*bank_index)[i])};
#ifdef OPENMC_MPI
if (i > 0)
MPI_Recv(source_bank->data(), count[0], mpi::source_site, i, i,
mpi::intracomm, MPI_STATUS_IGNORE);
#endif
// now write the source_bank data again.
for (const auto& site : *source_bank) {
// particle is now at the iterator
// write it to the mcpl-file
mcpl_particle_t p;
p.position[0] = site.r.x;
p.position[1] = site.r.y;
p.position[2] = site.r.z;
// mcpl requires that the direction vector is unit length
// which is also the case in openmc
p.direction[0] = site.u.x;
p.direction[1] = site.u.y;
p.direction[2] = site.u.z;
// MCPL stores kinetic energy in [MeV], time in [ms]
p.ekin = site.E * 1e-6;
p.time = site.time * 1e3;
p.weight = site.wgt;
switch (site.particle) {
case ParticleType::neutron:
p.pdgcode = 2112;
break;
case ParticleType::photon:
p.pdgcode = 22;
break;
case ParticleType::electron:
p.pdgcode = 11;
break;
case ParticleType::positron:
p.pdgcode = -11;
break;
}
mcpl_add_particle(file_id, &p);
}
}
#ifdef OPENMC_MPI
// Restore state of source bank
std::copy(temp_source.begin(), temp_source.end(), source_bank->begin());
#endif
} else {
#ifdef OPENMC_MPI
MPI_Send(source_bank->data(), count_size, mpi::source_site, 0, mpi::rank,
mpi::intracomm);
#endif
}
}
#endif
//==============================================================================
void write_mcpl_source_point(const char* filename, bool surf_source_bank)
{
std::string filename_;
if (filename) {
filename_ = filename;
} else {
// Determine width for zero padding
int w = std::to_string(settings::n_max_batches).size();
filename_ = fmt::format("{0}source.{1:0{2}}.mcpl", settings::path_output,
simulation::current_batch, w);
}
#ifdef OPENMC_MCPL
mcpl_outfile_t file_id;
std::string line;
if (mpi::master) {
file_id = mcpl_create_outfile(filename_.c_str());
if (VERSION_DEV) {
line = fmt::format("OpenMC {0}.{1}.{2}-development", VERSION_MAJOR,
VERSION_MINOR, VERSION_RELEASE);
} else {
line = fmt::format(
"OpenMC {0}.{1}.{2}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE);
}
mcpl_hdr_set_srcname(file_id, line.c_str());
}
write_mcpl_source_bank(file_id, surf_source_bank);
if (mpi::master) {
mcpl_close_outfile(file_id);
}
#endif
}
} // namespace openmc

108
src/ncrystal_interface.cpp Normal file
View file

@ -0,0 +1,108 @@
#include "openmc/ncrystal_interface.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/random_lcg.h"
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
#ifdef NCRYSTAL
const bool NCRYSTAL_ENABLED = true;
#else
const bool NCRYSTAL_ENABLED = false;
#endif
//==============================================================================
// NCrystal wrapper class for the OpenMC random number generator
//==============================================================================
#ifdef NCRYSTAL
class NCrystalRNGWrapper : public NCrystal::RNGStream {
public:
constexpr NCrystalRNGWrapper(uint64_t* seed) noexcept : openmc_seed_(seed) {}
protected:
double actualGenerate() override
{
return std::max<double>(
std::numeric_limits<double>::min(), prn(openmc_seed_));
}
private:
uint64_t* openmc_seed_;
};
#endif
//==============================================================================
// NCrystal implementation
//==============================================================================
NCrystalMat::NCrystalMat(const std::string& cfg)
{
#ifdef NCRYSTAL
cfg_ = cfg;
ptr_ = NCrystal::FactImpl::createScatter(cfg);
#else
fatal_error("Your build of OpenMC does not support NCrystal materials.");
#endif
}
#ifdef NCRYSTAL
std::string NCrystalMat::cfg() const
{
return cfg_;
}
double NCrystalMat::xs(const Particle& p) const
{
// Calculate scattering XS per atom with NCrystal, only once per material
NCrystal::CachePtr dummy_cache;
auto nc_energy = NCrystal::NeutronEnergy {p.E()};
return ptr_->crossSection(dummy_cache, nc_energy, {p.u().x, p.u().y, p.u().z})
.get();
}
void NCrystalMat::scatter(Particle& p) const
{
NCrystalRNGWrapper rng(p.current_seed()); // Initialize RNG
// create a cache pointer for multi thread physics
NCrystal::CachePtr dummy_cache;
auto nc_energy = NCrystal::NeutronEnergy {p.E()};
auto outcome = ptr_->sampleScatter(
dummy_cache, rng, nc_energy, {p.u().x, p.u().y, p.u().z});
// Modify attributes of particle
p.E() = outcome.ekin.get();
Direction u_old {p.u()};
p.u() =
Direction(outcome.direction[0], outcome.direction[1], outcome.direction[2]);
p.mu() = u_old.dot(p.u());
p.event_mt() = ELASTIC;
}
NCrystalMat::operator bool() const
{
return ptr_.get();
}
#endif
//==============================================================================
// Functions
//==============================================================================
void ncrystal_update_micro(double xs, NuclideMicroXS& micro)
{
if (micro.thermal > 0 || micro.thermal_elastic > 0) {
fatal_error("S(a,b) treatment and NCrystal are not compatible.");
}
// remove free atom cross section
// and replace it by scattering cross section per atom from NCrystal
micro.total = micro.total - micro.elastic + xs;
micro.elastic = xs;
}
} // namespace openmc

4
src/nuclide.cpp
View file

@ -519,7 +519,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const
case EmissionMode::prompt:
return (*fission_rx_[0]->products_[0].yield_)(E);
case EmissionMode::delayed:
if (n_precursor_ > 0) {
if (n_precursor_ > 0 && settings::create_delayed_neutrons) {
auto rx = fission_rx_[0];
if (group >= 1 && group < rx->products_.size()) {
// If delayed group specified, determine yield immediately
@ -544,7 +544,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const
return 0.0;
}
case EmissionMode::total:
if (total_nu_) {
if (total_nu_ && settings::create_delayed_neutrons) {
return (*total_nu_)(E);
} else {
return (*fission_rx_[0]->products_[0].yield_)(E);

8
src/physics.cpp
View file

@ -10,6 +10,7 @@
#include "openmc/material.h"
#include "openmc/math_functions.h"
#include "openmc/message_passing.h"
#include "openmc/ncrystal_interface.h"
#include "openmc/nuclide.h"
#include "openmc/photon.h"
#include "openmc/physics_common.h"
@ -140,7 +141,12 @@ void sample_neutron_reaction(Particle& p)
// Sample a scattering reaction and determine the secondary energy of the
// exiting neutron
scatter(p, i_nuclide);
const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat();
if (ncrystal_mat && p.E() < NCRYSTAL_MAX_ENERGY) {
ncrystal_mat.scatter(p);
} else {
scatter(p, i_nuclide);
}
// Advance URR seed stream 'N' times after energy changes
if (p.E() != p.E_last()) {

6
src/plot.cpp
View file

@ -141,6 +141,12 @@ void read_plots_xml()
doc.load_file(filename.c_str());
pugi::xml_node root = doc.document_element();
read_plots_xml(root);
}
void read_plots_xml(pugi::xml_node root)
{
for (auto node : root.children("plot")) {
model::plots.emplace_back(node);
model::plot_map[model::plots.back().id_] = model::plots.size() - 1;

52
src/settings.cpp
View file

@ -18,6 +18,7 @@
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/mcpl_interface.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/output.h"
@ -43,6 +44,7 @@ bool assume_separate {false};
bool check_overlaps {false};
bool cmfd_run {false};
bool confidence_intervals {false};
bool create_delayed_neutrons {true};
bool create_fission_neutrons {true};
bool delayed_photon_scaling {true};
bool entropy_on {false};
@ -60,7 +62,9 @@ bool run_CE {true};
bool source_latest {false};
bool source_separate {false};
bool source_write {true};
bool source_mcpl_write {false};
bool surf_source_write {false};
bool surf_mcpl_write {false};
bool surf_source_read {false};
bool survival_biasing {false};
bool temperature_multipole {false};
@ -215,16 +219,16 @@ void read_settings_xml()
{
using namespace settings;
using namespace pugi;
// Check if settings.xml exists
std::string filename = path_input + "settings.xml";
std::string filename = settings::path_input + "settings.xml";
if (!file_exists(filename)) {
if (run_mode != RunMode::PLOTTING) {
fatal_error(
fmt::format("Settings XML file '{}' does not exist! In order "
"to run OpenMC, you first need a set of input files; at a "
"minimum, this "
"includes settings.xml, geometry.xml, and materials.xml. "
"includes settings.xml, geometry.xml, and materials.xml "
"or a single XML file containing all of these files. "
"Please consult "
"the user's guide at https://docs.openmc.org for further "
"information.",
@ -256,8 +260,17 @@ void read_settings_xml()
if (verbosity >= 2)
title();
}
write_message("Reading settings XML file...", 5);
read_settings_xml(root);
}
void read_settings_xml(pugi::xml_node root)
{
using namespace settings;
using namespace pugi;
// Find if a multi-group or continuous-energy simulation is desired
if (check_for_node(root, "energy_mode")) {
std::string temp_str = get_node_value(root, "energy_mode", true, true);
@ -429,7 +442,12 @@ void read_settings_xml()
for (pugi::xml_node node : root.children("source")) {
if (check_for_node(node, "file")) {
auto path = get_node_value(node, "file", false, true);
model::external_sources.push_back(make_unique<FileSource>(path));
if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) {
auto sites = mcpl_source_sites(path);
model::external_sources.push_back(make_unique<FileSource>(sites));
} else {
model::external_sources.push_back(make_unique<FileSource>(path));
}
} else if (check_for_node(node, "library")) {
// Get shared library path and parameters
auto path = get_node_value(node, "library", false, true);
@ -647,6 +665,15 @@ void read_settings_xml()
if (check_for_node(node_sp, "write")) {
source_write = get_node_value_bool(node_sp, "write");
}
if (check_for_node(node_sp, "mcpl")) {
source_mcpl_write = get_node_value_bool(node_sp, "mcpl");
// Make sure MCPL support is enabled
if (source_mcpl_write && !MCPL_ENABLED) {
fatal_error(
"Your build of OpenMC does not support writing MCPL source files.");
}
}
if (check_for_node(node_sp, "overwrite_latest")) {
source_latest = get_node_value_bool(node_sp, "overwrite_latest");
source_separate = source_latest;
@ -677,9 +704,18 @@ void read_settings_xml()
max_surface_particles =
std::stoll(get_node_value(node_ssw, "max_particles"));
}
if (check_for_node(node_ssw, "mcpl")) {
surf_mcpl_write = get_node_value_bool(node_ssw, "mcpl");
// Make sure MCPL support is enabled
if (surf_mcpl_write && !MCPL_ENABLED) {
fatal_error("Your build of OpenMC does not support writing MCPL "
"surface source files.");
}
}
}
// If source is not seperate and is to be written out in the statepoint file,
// If source is not separate and is to be written out in the statepoint file,
// make sure that the sourcepoint batch numbers are contained in the
// statepoint list
if (!source_separate) {
@ -838,6 +874,12 @@ void read_settings_xml()
}
}
// Check whether create delayed neutrons in fission
if (check_for_node(root, "create_delayed_neutrons")) {
create_delayed_neutrons =
get_node_value_bool(root, "create_delayed_neutrons");
}
// Check whether create fission sites
if (run_mode == RunMode::FIXED_SOURCE) {
if (check_for_node(root, "create_fission_neutrons")) {

25
src/simulation.cpp
View file

@ -8,6 +8,7 @@
#include "openmc/event.h"
#include "openmc/geometry_aux.h"
#include "openmc/material.h"
#include "openmc/mcpl_interface.h"
#include "openmc/message_passing.h"
#include "openmc/nuclide.h"
#include "openmc/output.h"
@ -391,21 +392,35 @@ void finalize_batch()
// Write out a separate source point if it's been specified for this batch
if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
settings::source_write && settings::source_separate) {
write_source_point(nullptr);
if (settings::source_mcpl_write) {
write_mcpl_source_point(nullptr);
} else {
write_source_point(nullptr);
}
}
// Write a continously-overwritten source point if requested.
if (settings::source_latest) {
auto filename = settings::path_output + "source.h5";
write_source_point(filename.c_str());
if (settings::source_mcpl_write) {
auto filename = settings::path_output + "source.mcpl";
write_mcpl_source_point(filename.c_str());
} else {
auto filename = settings::path_output + "source.h5";
write_source_point(filename.c_str());
}
}
}
// Write out surface source if requested.
if (settings::surf_source_write &&
simulation::current_batch == settings::n_batches) {
auto filename = settings::path_output + "surface_source.h5";
write_source_point(filename.c_str(), true);
if (settings::surf_mcpl_write) {
auto filename = settings::path_output + "surface_source.mcpl";
write_mcpl_source_point(filename.c_str(), true);
} else {
auto filename = settings::path_output + "surface_source.h5";
write_source_point(filename.c_str(), true);
}
}
}

5
src/tallies/tally.cpp
View file

@ -719,6 +719,11 @@ void read_tallies_xml()
doc.load_file(filename.c_str());
pugi::xml_node root = doc.document_element();
read_tallies_xml(root);
}
void read_tallies_xml(pugi::xml_node root)
{
// Check for <assume_separate> setting
if (check_for_node(root, "assume_separate")) {
settings::assume_separate = get_node_value_bool(root, "assume_separate");

87
src/tallies/trigger.cpp
View file

@ -37,6 +37,10 @@ std::pair<double, double> get_tally_uncertainty(
int n = tally->n_realizations_;
auto mean = sum / n;
// if the result has no contributions, return an invalid pair
if (mean == 0) return {-1 , -1};
double std_dev = std::sqrt((sum_sq / n - mean * mean) / (n - 1));
double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.;
@ -68,42 +72,49 @@ void check_tally_triggers(double& ratio, int& tally_id, int& score)
const auto& results = t.results_;
for (auto filter_index = 0; filter_index < results.shape()[0];
++filter_index) {
for (auto score_index = 0; score_index < results.shape()[1];
++score_index) {
// Compute the tally uncertainty metrics.
auto uncert_pair =
get_tally_uncertainty(i_tally, score_index, filter_index);
double std_dev = uncert_pair.first;
double rel_err = uncert_pair.second;
// Compute the tally uncertainty metrics.
auto uncert_pair =
get_tally_uncertainty(i_tally, trigger.score_index, filter_index);
// Pick out the relevant uncertainty metric for this trigger.
double uncertainty;
switch (trigger.metric) {
case TriggerMetric::variance:
uncertainty = std_dev * std_dev;
break;
case TriggerMetric::standard_deviation:
uncertainty = std_dev;
break;
case TriggerMetric::relative_error:
uncertainty = rel_err;
break;
case TriggerMetric::not_active:
UNREACHABLE();
}
// if there is a score without contributions, set ratio to inf and
// exit early
if (uncert_pair.first == -1) {
ratio = INFINITY;
score = t.scores_[trigger.score_index];
tally_id = t.id_;
return;
}
// Compute the uncertainty / threshold ratio.
double this_ratio = uncertainty / trigger.threshold;
if (trigger.metric == TriggerMetric::variance) {
this_ratio = std::sqrt(ratio);
}
double std_dev = uncert_pair.first;
double rel_err = uncert_pair.second;
// If this is the most uncertain value, set the output variables.
if (this_ratio > ratio) {
ratio = this_ratio;
score = t.scores_[trigger.score_index];
tally_id = t.id_;
}
// Pick out the relevant uncertainty metric for this trigger.
double uncertainty;
switch (trigger.metric) {
case TriggerMetric::variance:
uncertainty = std_dev * std_dev;
break;
case TriggerMetric::standard_deviation:
uncertainty = std_dev;
break;
case TriggerMetric::relative_error:
uncertainty = rel_err;
break;
case TriggerMetric::not_active:
UNREACHABLE();
}
// Compute the uncertainty / threshold ratio.
double this_ratio = uncertainty / trigger.threshold;
if (trigger.metric == TriggerMetric::variance) {
this_ratio = std::sqrt(ratio);
}
// If this is the most uncertain value, set the output variables.
if (this_ratio > ratio) {
ratio = this_ratio;
score = t.scores_[trigger.score_index];
tally_id = t.id_;
}
}
}
@ -181,9 +192,13 @@ void check_triggers()
"eigenvalue",
keff_ratio);
} else {
msg = fmt::format(
"Triggers unsatisfied, max unc./thresh. is {} for {} in tally {}",
tally_ratio, reaction_name(score), tally_id);
if (tally_ratio == INFINITY) {
msg = fmt::format("Triggers unsatisfied, no result tallied for score {} in tally {}", reaction_name(score), tally_id);
} else{
msg = fmt::format(
"Triggers unsatisfied, max unc./thresh. is {} for {} in tally {}",
tally_ratio, reaction_name(score), tally_id);
}
}
write_message(msg, 7);

32
src/thermal.cpp
View file

@ -119,22 +119,24 @@ ThermalScattering::ThermalScattering(
if (!found) {
// If no pairs found, check if the desired temperature falls within
// bounds' tolerance
if (std::abs(T - temps_available[0]) <= settings::temperature_tolerance){
if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temps_available[0])) ==
temps_to_read.end()) {
temps_to_read.push_back(std::round(temps_available[0]));
}}
else if (std::abs(T - temps_available[n - 1]) <= settings::temperature_tolerance){
if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temps_available[n - 1])) ==
temps_to_read.end()){
temps_to_read.push_back(std::round(temps_available[n - 1]));
}}
else {
if (std::abs(T - temps_available[0]) <=
settings::temperature_tolerance) {
if (std::find(temps_to_read.begin(), temps_to_read.end(),
std::round(temps_available[0])) == temps_to_read.end()) {
temps_to_read.push_back(std::round(temps_available[0]));
}
} else if (std::abs(T - temps_available[n - 1]) <=
settings::temperature_tolerance) {
if (std::find(temps_to_read.begin(), temps_to_read.end(),
std::round(temps_available[n - 1])) == temps_to_read.end()) {
temps_to_read.push_back(std::round(temps_available[n - 1]));
}
} else {
fatal_error(
fmt::format("Nuclear data library does not contain cross "
"sections for {} at temperatures that bound {} K.",
name_, std::round(T)));
}
fmt::format("Nuclear data library does not contain cross "
"sections for {} at temperatures that bound {} K.",
name_, std::round(T)));
}
}
}
}

5
tests/regression_tests/diff_tally/test.py
View file

@ -103,9 +103,8 @@ class DiffTallyTestHarness(PyAPITestHarness):
sp = openmc.StatePoint(statepoint)
# Extract the tally data as a Pandas DataFrame.
df = pd.DataFrame()
for t in sp.tallies.values():
df = df.append(t.get_pandas_dataframe(), ignore_index=True)
tally_dfs = [t.get_pandas_dataframe() for t in sp.tallies.values()]
df = pd.concat(tally_dfs, ignore_index=True)
# Extract the relevant data as a CSV string.
cols = ('d_material', 'd_nuclide', 'd_variable', 'score', 'mean',

0
tests/regression_tests/model_xml/__init__.py Normal file
View file

38
tests/regression_tests/model_xml/adj_cell_rotation_inputs_true.dat Normal file
View file

@ -0,0 +1,38 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material depletable="true" id="1">
<density units="g/cc" value="10.0" />
<nuclide ao="1.0" name="U235" />
</material>
<material id="2">
<density units="g/cc" value="0.1" />
<nuclide ao="0.1" name="H1" />
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1" />
<cell id="2" material="2" region="1" universe="1" />
<cell fill="1" id="3" region="2 -3 4 -5 6 -8" rotation="10 20 30" universe="2" />
<cell fill="1" id="4" region="2 -3 4 -5 8 -7" translation="0 0 15" universe="2" />
<surface coeffs="1.0 0.0 0.0 5.0" id="1" type="sphere" />
<surface boundary="vacuum" coeffs="-7.5" id="2" name="minimum x" type="x-plane" />
<surface boundary="vacuum" coeffs="7.5" id="3" name="maximum x" type="x-plane" />
<surface boundary="vacuum" coeffs="-7.5" id="4" name="minimum y" type="y-plane" />
<surface boundary="vacuum" coeffs="7.5" id="5" name="maximum y" type="y-plane" />
<surface boundary="vacuum" coeffs="-7.5" id="6" type="z-plane" />
<surface boundary="vacuum" coeffs="22.5" id="7" type="z-plane" />
<surface coeffs="7.5" id="8" type="z-plane" />
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>10000</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>-4.0 -4.0 -4.0 4.0 4.0 4.0</parameters>
</space>
</source>
</settings>
</model>

23
tests/regression_tests/model_xml/energy_laws_inputs_true.dat Normal file
View file

@ -0,0 +1,23 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material depletable="true" id="1">
<density units="g/cm3" value="20.0" />
<nuclide ao="1.0" name="U233" />
<nuclide ao="1.0" name="Am244" />
<nuclide ao="1.0" name="H2" />
<nuclide ao="1.0" name="Na23" />
<nuclide ao="1.0" name="Ta181" />
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1" />
<surface boundary="reflective" coeffs="0.0 0.0 0.0 100.0" id="1" type="sphere" />
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
</settings>
</model>

67
tests/regression_tests/model_xml/inputs_true.dat Normal file
View file

@ -0,0 +1,67 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="6">
<density units="g/cm3" value="2.6989" />
<nuclide ao="1.0" name="Al27" />
</material>
</materials>
<geometry>
<cell id="12" material="void" region="-16 17 -18" universe="10" />
<cell id="13" material="6" region="-16 18 -19" universe="10" />
<cell id="14" material="void" region="~(-16 17 -19)" universe="10" />
<surface id="16" type="x-cylinder" boundary="vacuum" coeffs="0.0 0.0 1.0" />
<surface id="17" type="x-plane" boundary="vacuum" coeffs="-1.0" />
<surface id="18" type="x-plane" coeffs="1.0" />
<surface id="19" type="x-plane" boundary="vacuum" coeffs="1000000000.0" />
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>10000</particles>
<batches>1</batches>
<source strength="1.0">
<space type="point">
<parameters>0 0 0</parameters>
</space>
<angle type="monodirectional" reference_uvw="1.0 0.0 0.0" />
<energy type="discrete">
<parameters>14000000.0 1.0</parameters>
</energy>
</source>
<electron_treatment>ttb</electron_treatment>
<photon_transport>true</photon_transport>
<cutoff>
<energy_photon>1000.0</energy_photon>
</cutoff>
</settings>
<tallies>
<filter id="1" type="surface">
<bins>16</bins>
</filter>
<filter id="2" type="particle">
<bins>neutron photon electron positron</bins>
</filter>
<tally id="1">
<filters>1 2</filters>
<scores>current</scores>
</tally>
<tally id="2">
<filters>2</filters>
<nuclides>Al27 total</nuclides>
<scores>total (n,gamma)</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="3">
<filters>2</filters>
<nuclides>Al27 total</nuclides>
<scores>total heating (n,gamma)</scores>
<estimator>collision</estimator>
</tally>
<tally id="4">
<filters>2</filters>
<nuclides>Al27 total</nuclides>
<scores>total heating (n,gamma)</scores>
<estimator>analog</estimator>
</tally>
</tallies>
</model>

53
tests/regression_tests/model_xml/lattice_multiple_inputs_true.dat Normal file
View file

@ -0,0 +1,53 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material depletable="true" id="1" name="UO2">
<density units="g/cm3" value="10.0" />
<nuclide ao="1.0" name="U235" />
<nuclide ao="2.0" name="O16" />
</material>
<material id="2" name="light water">
<density units="g/cm3" value="1.0" />
<nuclide ao="2.0" name="H1" />
<nuclide ao="1.0" name="O16" />
<sab name="c_H_in_H2O" />
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1" />
<cell id="2" material="2" region="1" universe="1" />
<cell id="3" material="1" region="-2" universe="2" />
<cell id="4" material="2" region="2" universe="2" />
<cell fill="3" id="5" universe="4" />
<cell fill="5" id="6" region="3 -4 5 -6" universe="6" />
<lattice id="3">
<pitch>1.2 1.2</pitch>
<outer>1</outer>
<dimension>2 2</dimension>
<lower_left>-1.2 -1.2</lower_left>
<universes>
2 1
1 1 </universes>
</lattice>
<lattice id="5">
<pitch>2.4 2.4</pitch>
<dimension>2 2</dimension>
<lower_left>-2.4 -2.4</lower_left>
<universes>
4 4
4 4 </universes>
</lattice>
<surface coeffs="0.0 0.0 0.4" id="1" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.5" id="2" type="z-cylinder" />
<surface boundary="reflective" coeffs="-2.4" id="3" name="minimum x" type="x-plane" />
<surface boundary="reflective" coeffs="2.4" id="4" name="maximum x" type="x-plane" />
<surface boundary="reflective" coeffs="-2.4" id="5" name="minimum y" type="y-plane" />
<surface boundary="reflective" coeffs="2.4" id="6" name="maximum y" type="y-plane" />
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
</settings>
</model>

67
tests/regression_tests/model_xml/photon_production_inputs_true.dat Normal file
View file

@ -0,0 +1,67 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1">
<density units="g/cm3" value="2.6989" />
<nuclide ao="1.0" name="Al27" />
</material>
</materials>
<geometry>
<cell id="1" material="void" region="-1 2 -3" universe="1" />
<cell id="2" material="1" region="-1 3 -4" universe="1" />
<cell id="3" material="void" region="~(-1 2 -4)" universe="1" />
<surface boundary="vacuum" coeffs="0.0 0.0 1.0" id="1" type="x-cylinder" />
<surface boundary="vacuum" coeffs="-1.0" id="2" type="x-plane" />
<surface coeffs="1.0" id="3" type="x-plane" />
<surface boundary="vacuum" coeffs="1000000000.0" id="4" type="x-plane" />
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>10000</particles>
<batches>1</batches>
<source strength="1.0">
<space type="point">
<parameters>0 0 0</parameters>
</space>
<angle reference_uvw="1.0 0.0 0.0" type="monodirectional" />
<energy type="discrete">
<parameters>14000000.0 1.0</parameters>
</energy>
</source>
<electron_treatment>ttb</electron_treatment>
<photon_transport>true</photon_transport>
<cutoff>
<energy_photon>1000.0</energy_photon>
</cutoff>
</settings>
<tallies>
<filter id="1" type="surface">
<bins>1</bins>
</filter>
<filter id="2" type="particle">
<bins>neutron photon electron positron</bins>
</filter>
<tally id="1">
<filters>1 2</filters>
<scores>current</scores>
</tally>
<tally id="2">
<filters>2</filters>
<nuclides>Al27 total</nuclides>
<scores>total (n,gamma)</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="3">
<filters>2</filters>
<nuclides>Al27 total</nuclides>
<scores>total heating (n,gamma)</scores>
<estimator>collision</estimator>
</tally>
<tally id="4">
<filters>2</filters>
<nuclides>Al27 total</nuclides>
<scores>total heating (n,gamma)</scores>
<estimator>analog</estimator>
</tally>
</tallies>
</model>

100
tests/regression_tests/model_xml/test.py Normal file
View file

@ -0,0 +1,100 @@
from difflib import unified_diff
import glob
import filecmp
import os
from pathlib import Path
import openmc
import pytest
from tests.testing_harness import PyAPITestHarness, colorize
# use a few models from other tests to make sure the same results are
# produced when using a single model.xml file as input
from ..adj_cell_rotation.test import model as adj_cell_rotation_model
from ..lattice_multiple.test import model as lattice_multiple_model
from ..energy_laws.test import model as energy_laws_model
from ..photon_production.test import model as photon_production_model
class ModelXMLTestHarness(PyAPITestHarness):
"""Accept a results file to check against and assume inputs_true is the contents of a model.xml file.
"""
def __init__(self, model=None, inputs_true=None, results_true=None):
statepoint_name = f'statepoint.{model.settings.batches}.h5'
super().__init__(statepoint_name, model, inputs_true)
self.results_true = 'results_true.dat' if results_true is None else results_true
def _build_inputs(self):
self._model.export_to_model_xml()
def _get_inputs(self):
return open('model.xml').read()
def _compare_results(self):
"""Make sure the current results agree with the reference."""
compare = filecmp.cmp('results_test.dat', self.results_true)
if not compare:
expected = open(self.results_true).readlines()
actual = open('results_test.dat').readlines()
diff = unified_diff(expected, actual, self.results_true,
'results_test.dat')
print('Result differences:')
print(''.join(colorize(diff)))
os.rename('results_test.dat', 'results_error.dat')
assert compare, 'Results do not agree'
def _cleanup(self):
super()._cleanup()
if os.path.exists('model.xml'):
os.remove('model.xml')
test_names = [
'adj_cell_rotation',
'lattice_multiple',
'energy_laws',
'photon_production'
]
@pytest.mark.parametrize("test_name", test_names, ids=lambda test: test)
def test_model_xml(test_name, request):
openmc.reset_auto_ids()
test_path = '../' + test_name
results = test_path + "/results_true.dat"
inputs = test_name + "_inputs_true.dat"
model_name = test_name + "_model"
harness = ModelXMLTestHarness(request.getfixturevalue(model_name), inputs, results)
harness.main()
def test_input_arg(run_in_tmpdir):
pincell = openmc.examples.pwr_pin_cell()
pincell.settings.particles = 100
# export to separate XML files and run
pincell.export_to_xml()
openmc.run()
# make sure the executable isn't falling back on the separate XMLs
for f in glob.glob('*.xml'):
os.remove(f)
# now export to a single XML file with a custom name
pincell.export_to_model_xml('pincell.xml')
assert Path('pincell.xml').exists()
# run by specifying that single file
openmc.run(path_input='pincell.xml')
# check that this works for plotting too
openmc.plot_geometry(path_input='pincell.xml')
# now ensure we get an error for an incorrect filename,
# even in the presence of other, valid XML files
pincell.export_to_model_xml()
with pytest.raises(RuntimeError, match='ex-em-ell.xml'):
openmc.run(path_input='ex-em-ell.xml')

0
tests/regression_tests/ncrystal/__init__.py Normal file
View file

45
tests/regression_tests/ncrystal/inputs_true.dat Normal file
View file

@ -0,0 +1,45 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="-1" universe="1" />
<cell id="2" material="void" region="1 -2" universe="1" />
<surface coeffs="0.0 0.0 0.0 0.1" id="1" type="sphere" />
<surface boundary="vacuum" coeffs="0.0 0.0 0.0 100" id="2" type="sphere" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material cfg="Al_sg225.ncmat" id="1" temperature="293.15">
<density units="g/cm3" value="2.6986455176922477" />
<nuclide ao="1.0" name="Al27" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>100000</particles>
<batches>10</batches>
<source strength="1.0">
<space type="point">
<parameters>0 0 -20</parameters>
</space>
<angle reference_uvw="0.0 0.0 1.0" type="monodirectional" />
<energy type="discrete">
<parameters>0.012 1.0</parameters>
</energy>
</source>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<filter id="1" type="surface">
<bins>1</bins>
</filter>
<filter id="2" type="polar">
<bins>0.0 0.017453292519943295 0.03490658503988659 0.05235987755982989 0.06981317007977318 0.08726646259971647 0.10471975511965978 0.12217304763960307 0.13962634015954636 0.15707963267948966 0.17453292519943295 0.19198621771937624 0.20943951023931956 0.22689280275926285 0.24434609527920614 0.2617993877991494 0.2792526803190927 0.29670597283903605 0.3141592653589793 0.33161255787892263 0.3490658503988659 0.3665191429188092 0.3839724354387525 0.4014257279586958 0.4188790204786391 0.4363323129985824 0.4537856055185257 0.47123889803846897 0.4886921905584123 0.5061454830783556 0.5235987755982988 0.5410520681182421 0.5585053606381855 0.5759586531581288 0.5934119456780721 0.6108652381980153 0.6283185307179586 0.6457718232379019 0.6632251157578453 0.6806784082777885 0.6981317007977318 0.7155849933176751 0.7330382858376184 0.7504915783575618 0.767944870877505 0.7853981633974483 0.8028514559173916 0.8203047484373349 0.8377580409572782 0.8552113334772214 0.8726646259971648 0.8901179185171081 0.9075712110370514 0.9250245035569946 0.9424777960769379 0.9599310885968813 0.9773843811168246 0.9948376736367679 1.0122909661567112 1.0297442586766545 1.0471975511965976 1.064650843716541 1.0821041362364843 1.0995574287564276 1.117010721276371 1.1344640137963142 1.1519173063162575 1.1693705988362009 1.1868238913561442 1.2042771838760873 1.2217304763960306 1.239183768915974 1.2566370614359172 1.2740903539558606 1.2915436464758039 1.3089969389957472 1.3264502315156905 1.3439035240356338 1.361356816555577 1.3788101090755203 1.3962634015954636 1.413716694115407 1.4311699866353502 1.4486232791552935 1.4660765716752369 1.4835298641951802 1.5009831567151235 1.5184364492350666 1.53588974175501 1.5533430342749532 1.5707963267948966 1.5882496193148399 1.6057029118347832 1.6231562043547265 1.6406094968746698 1.6580627893946132 1.6755160819145565 1.6929693744344996 1.710422666954443 1.7278759594743862 1.7453292519943295 1.7627825445142729 1.7802358370342162 1.7976891295541595 1.8151424220741028 1.8325957145940461 1.8500490071139892 1.8675022996339325 1.8849555921538759 1.9024088846738192 1.9198621771937625 1.9373154697137058 1.9547687622336491 1.9722220547535925 1.9896753472735358 2.007128639793479 2.0245819323134224 2.0420352248333655 2.059488517353309 2.076941809873252 2.0943951023931953 2.111848394913139 2.129301687433082 2.1467549799530254 2.1642082724729685 2.181661564992912 2.199114857512855 2.2165681500327987 2.234021442552742 2.251474735072685 2.2689280275926285 2.2863813201125716 2.303834612632515 2.321287905152458 2.3387411976724017 2.356194490192345 2.3736477827122884 2.3911010752322315 2.4085543677521746 2.426007660272118 2.443460952792061 2.4609142453120048 2.478367537831948 2.4958208303518914 2.5132741228718345 2.530727415391778 2.548180707911721 2.5656340004316642 2.5830872929516078 2.600540585471551 2.6179938779914944 2.6354471705114375 2.652900463031381 2.670353755551324 2.6878070480712677 2.705260340591211 2.722713633111154 2.7401669256310974 2.7576202181510405 2.775073510670984 2.792526803190927 2.8099800957108707 2.827433388230814 2.8448866807507573 2.8623399732707004 2.8797932657906435 2.897246558310587 2.91469985083053 2.9321531433504737 2.949606435870417 2.9670597283903604 2.9845130209103035 3.001966313430247 3.01941960595019 3.036872898470133 3.0543261909900767 3.07177948351002 3.0892327760299634 3.1066860685499065 3.12413936106985 3.141592653589793</bins>
</filter>
<filter id="3" type="cellfrom">
<bins>1</bins>
</filter>
<tally id="1" name="angular distribution">
<filters>1 2 3</filters>
<scores>current</scores>
</tally>
</tallies>

181
tests/regression_tests/ncrystal/results_true.dat Normal file
View file

@ -0,0 +1,181 @@
surface polar low [rad] polar high [rad] cellfrom nuclide score mean std. dev.
0 1 0.00e+00 1.75e-02 1 total current 9.82e-01 1.43e-04
1 1 1.75e-02 3.49e-02 1 total current 0.00e+00 0.00e+00
2 1 3.49e-02 5.24e-02 1 total current 1.00e-06 1.00e-06
3 1 5.24e-02 6.98e-02 1 total current 0.00e+00 0.00e+00
4 1 6.98e-02 8.73e-02 1 total current 0.00e+00 0.00e+00
5 1 8.73e-02 1.05e-01 1 total current 1.00e-06 1.00e-06
6 1 1.05e-01 1.22e-01 1 total current 0.00e+00 0.00e+00
7 1 1.22e-01 1.40e-01 1 total current 1.00e-06 1.00e-06
8 1 1.40e-01 1.57e-01 1 total current 1.00e-06 1.00e-06
9 1 1.57e-01 1.75e-01 1 total current 1.00e-06 1.00e-06
10 1 1.75e-01 1.92e-01 1 total current 0.00e+00 0.00e+00
11 1 1.92e-01 2.09e-01 1 total current 2.00e-06 1.33e-06
12 1 2.09e-01 2.27e-01 1 total current 0.00e+00 0.00e+00
13 1 2.27e-01 2.44e-01 1 total current 1.00e-06 1.00e-06
14 1 2.44e-01 2.62e-01 1 total current 1.00e-06 1.00e-06
15 1 2.62e-01 2.79e-01 1 total current 2.00e-06 1.33e-06
16 1 2.79e-01 2.97e-01 1 total current 0.00e+00 0.00e+00
17 1 2.97e-01 3.14e-01 1 total current 2.00e-06 2.00e-06
18 1 3.14e-01 3.32e-01 1 total current 1.00e-06 1.00e-06
19 1 3.32e-01 3.49e-01 1 total current 1.00e-06 1.00e-06
20 1 3.49e-01 3.67e-01 1 total current 2.00e-06 1.33e-06
21 1 3.67e-01 3.84e-01 1 total current 0.00e+00 0.00e+00
22 1 3.84e-01 4.01e-01 1 total current 2.00e-06 1.33e-06
23 1 4.01e-01 4.19e-01 1 total current 2.00e-06 1.33e-06
24 1 4.19e-01 4.36e-01 1 total current 1.00e-06 1.00e-06
25 1 4.36e-01 4.54e-01 1 total current 2.00e-06 2.00e-06
26 1 4.54e-01 4.71e-01 1 total current 5.00e-06 2.24e-06
27 1 4.71e-01 4.89e-01 1 total current 4.00e-06 1.63e-06
28 1 4.89e-01 5.06e-01 1 total current 3.00e-06 1.53e-06
29 1 5.06e-01 5.24e-01 1 total current 3.00e-06 1.53e-06
30 1 5.24e-01 5.41e-01 1 total current 3.00e-06 1.53e-06
31 1 5.41e-01 5.59e-01 1 total current 7.00e-06 2.13e-06
32 1 5.59e-01 5.76e-01 1 total current 3.00e-06 1.53e-06
33 1 5.76e-01 5.93e-01 1 total current 3.00e-06 1.53e-06
34 1 5.93e-01 6.11e-01 1 total current 2.00e-06 1.33e-06
35 1 6.11e-01 6.28e-01 1 total current 2.00e-06 1.33e-06
36 1 6.28e-01 6.46e-01 1 total current 3.00e-06 2.13e-06
37 1 6.46e-01 6.63e-01 1 total current 3.00e-06 1.53e-06
38 1 6.63e-01 6.81e-01 1 total current 2.00e-06 1.33e-06
39 1 6.81e-01 6.98e-01 1 total current 3.00e-06 1.53e-06
40 1 6.98e-01 7.16e-01 1 total current 1.00e-06 1.00e-06
41 1 7.16e-01 7.33e-01 1 total current 6.00e-06 2.67e-06
42 1 7.33e-01 7.50e-01 1 total current 6.00e-06 2.21e-06
43 1 7.50e-01 7.68e-01 1 total current 7.00e-06 3.35e-06
44 1 7.68e-01 7.85e-01 1 total current 6.00e-06 2.67e-06
45 1 7.85e-01 8.03e-01 1 total current 6.00e-06 2.21e-06
46 1 8.03e-01 8.20e-01 1 total current 7.00e-06 2.13e-06
47 1 8.20e-01 8.38e-01 1 total current 5.00e-06 2.24e-06
48 1 8.38e-01 8.55e-01 1 total current 3.00e-06 1.53e-06
49 1 8.55e-01 8.73e-01 1 total current 8.00e-06 3.27e-06
50 1 8.73e-01 8.90e-01 1 total current 7.00e-06 2.13e-06
51 1 8.90e-01 9.08e-01 1 total current 6.00e-06 2.21e-06
52 1 9.08e-01 9.25e-01 1 total current 4.00e-06 2.21e-06
53 1 9.25e-01 9.42e-01 1 total current 1.20e-05 3.27e-06
54 1 9.42e-01 9.60e-01 1 total current 8.00e-06 3.89e-06
55 1 9.60e-01 9.77e-01 1 total current 1.20e-05 4.42e-06
56 1 9.77e-01 9.95e-01 1 total current 7.00e-06 3.67e-06
57 1 9.95e-01 1.01e+00 1 total current 1.20e-05 3.89e-06
58 1 1.01e+00 1.03e+00 1 total current 9.00e-06 2.33e-06
59 1 1.03e+00 1.05e+00 1 total current 6.00e-06 2.67e-06
60 1 1.05e+00 1.06e+00 1 total current 7.00e-06 2.13e-06
61 1 1.06e+00 1.08e+00 1 total current 5.00e-06 2.24e-06
62 1 1.08e+00 1.10e+00 1 total current 1.20e-05 2.91e-06
63 1 1.10e+00 1.12e+00 1 total current 1.30e-05 3.35e-06
64 1 1.12e+00 1.13e+00 1 total current 9.00e-06 3.14e-06
65 1 1.13e+00 1.15e+00 1 total current 1.20e-05 3.89e-06
66 1 1.15e+00 1.17e+00 1 total current 6.00e-06 2.21e-06
67 1 1.17e+00 1.19e+00 1 total current 5.11e-03 4.20e-05
68 1 1.19e+00 1.20e+00 1 total current 8.00e-06 3.27e-06
69 1 1.20e+00 1.22e+00 1 total current 1.30e-05 4.48e-06
70 1 1.22e+00 1.24e+00 1 total current 1.20e-05 3.89e-06
71 1 1.24e+00 1.26e+00 1 total current 1.50e-05 4.28e-06
72 1 1.26e+00 1.27e+00 1 total current 7.00e-06 3.00e-06
73 1 1.27e+00 1.29e+00 1 total current 1.60e-05 3.71e-06
74 1 1.29e+00 1.31e+00 1 total current 9.00e-06 3.79e-06
75 1 1.31e+00 1.33e+00 1 total current 1.30e-05 2.60e-06
76 1 1.33e+00 1.34e+00 1 total current 1.60e-05 3.06e-06
77 1 1.34e+00 1.36e+00 1 total current 1.40e-05 2.67e-06
78 1 1.36e+00 1.38e+00 1 total current 1.40e-05 6.86e-06
79 1 1.38e+00 1.40e+00 1 total current 1.50e-05 4.28e-06
80 1 1.40e+00 1.41e+00 1 total current 3.26e-03 6.93e-05
81 1 1.41e+00 1.43e+00 1 total current 1.40e-05 3.71e-06
82 1 1.43e+00 1.45e+00 1 total current 1.30e-05 3.00e-06
83 1 1.45e+00 1.47e+00 1 total current 1.30e-05 3.67e-06
84 1 1.47e+00 1.48e+00 1 total current 1.70e-05 5.59e-06
85 1 1.48e+00 1.50e+00 1 total current 1.70e-05 3.67e-06
86 1 1.50e+00 1.52e+00 1 total current 1.40e-05 3.40e-06
87 1 1.52e+00 1.54e+00 1 total current 2.50e-05 4.53e-06
88 1 1.54e+00 1.55e+00 1 total current 1.30e-05 5.39e-06
89 1 1.55e+00 1.57e+00 1 total current 1.80e-05 4.67e-06
90 1 1.57e+00 1.59e+00 1 total current 1.40e-05 4.52e-06
91 1 1.59e+00 1.61e+00 1 total current 1.70e-05 4.23e-06
92 1 1.61e+00 1.62e+00 1 total current 1.40e-05 3.71e-06
93 1 1.62e+00 1.64e+00 1 total current 1.00e-05 2.11e-06
94 1 1.64e+00 1.66e+00 1 total current 2.00e-05 4.22e-06
95 1 1.66e+00 1.68e+00 1 total current 2.30e-05 5.59e-06
96 1 1.68e+00 1.69e+00 1 total current 1.70e-05 6.51e-06
97 1 1.69e+00 1.71e+00 1 total current 1.30e-05 3.00e-06
98 1 1.71e+00 1.73e+00 1 total current 1.50e-05 4.01e-06
99 1 1.73e+00 1.75e+00 1 total current 1.70e-05 3.96e-06
100 1 1.75e+00 1.76e+00 1 total current 1.80e-05 5.12e-06
101 1 1.76e+00 1.78e+00 1 total current 2.50e-05 6.54e-06
102 1 1.78e+00 1.80e+00 1 total current 1.80e-05 3.59e-06
103 1 1.80e+00 1.82e+00 1 total current 1.50e-05 4.01e-06
104 1 1.82e+00 1.83e+00 1 total current 1.10e-05 4.07e-06
105 1 1.83e+00 1.85e+00 1 total current 1.50e-05 4.01e-06
106 1 1.85e+00 1.87e+00 1 total current 1.90e-05 4.82e-06
107 1 1.87e+00 1.88e+00 1 total current 2.20e-05 3.89e-06
108 1 1.88e+00 1.90e+00 1 total current 2.10e-05 3.79e-06
109 1 1.90e+00 1.92e+00 1 total current 1.50e-05 3.42e-06
110 1 1.92e+00 1.94e+00 1 total current 2.20e-05 5.12e-06
111 1 1.94e+00 1.95e+00 1 total current 2.10e-05 5.86e-06
112 1 1.95e+00 1.97e+00 1 total current 2.60e-05 4.27e-06
113 1 1.97e+00 1.99e+00 1 total current 2.20e-05 4.16e-06
114 1 1.99e+00 2.01e+00 1 total current 2.40e-05 5.42e-06
115 1 2.01e+00 2.02e+00 1 total current 1.60e-05 4.52e-06
116 1 2.02e+00 2.04e+00 1 total current 1.30e-05 3.35e-06
117 1 2.04e+00 2.06e+00 1 total current 1.90e-05 4.07e-06
118 1 2.06e+00 2.08e+00 1 total current 1.30e-05 3.00e-06
119 1 2.08e+00 2.09e+00 1 total current 1.40e-05 4.00e-06
120 1 2.09e+00 2.11e+00 1 total current 3.10e-05 5.47e-06
121 1 2.11e+00 2.13e+00 1 total current 2.30e-05 5.39e-06
122 1 2.13e+00 2.15e+00 1 total current 2.20e-05 4.67e-06
123 1 2.15e+00 2.16e+00 1 total current 1.70e-05 4.48e-06
124 1 2.16e+00 2.18e+00 1 total current 1.60e-05 4.00e-06
125 1 2.18e+00 2.20e+00 1 total current 1.80e-05 4.16e-06
126 1 2.20e+00 2.22e+00 1 total current 1.80e-05 4.42e-06
127 1 2.22e+00 2.23e+00 1 total current 1.80e-05 5.54e-06
128 1 2.23e+00 2.25e+00 1 total current 1.90e-05 4.33e-06
129 1 2.25e+00 2.27e+00 1 total current 1.00e-05 3.33e-06
130 1 2.27e+00 2.29e+00 1 total current 1.80e-05 4.42e-06
131 1 2.29e+00 2.30e+00 1 total current 4.15e-03 5.40e-05
132 1 2.30e+00 2.32e+00 1 total current 1.90e-05 2.33e-06
133 1 2.32e+00 2.34e+00 1 total current 2.30e-05 3.96e-06
134 1 2.34e+00 2.36e+00 1 total current 1.90e-05 3.48e-06
135 1 2.36e+00 2.37e+00 1 total current 1.60e-05 3.71e-06
136 1 2.37e+00 2.39e+00 1 total current 1.60e-05 4.27e-06
137 1 2.39e+00 2.41e+00 1 total current 2.30e-05 4.48e-06
138 1 2.41e+00 2.43e+00 1 total current 2.10e-05 3.48e-06
139 1 2.43e+00 2.44e+00 1 total current 1.60e-05 2.21e-06
140 1 2.44e+00 2.46e+00 1 total current 9.00e-06 2.77e-06
141 1 2.46e+00 2.48e+00 1 total current 1.30e-05 3.00e-06
142 1 2.48e+00 2.50e+00 1 total current 2.70e-05 3.67e-06
143 1 2.50e+00 2.51e+00 1 total current 1.90e-05 4.07e-06
144 1 2.51e+00 2.53e+00 1 total current 1.20e-05 4.16e-06
145 1 2.53e+00 2.55e+00 1 total current 1.30e-05 2.13e-06
146 1 2.55e+00 2.57e+00 1 total current 1.10e-05 2.33e-06
147 1 2.57e+00 2.58e+00 1 total current 1.50e-05 3.07e-06
148 1 2.58e+00 2.60e+00 1 total current 1.20e-05 2.49e-06
149 1 2.60e+00 2.62e+00 1 total current 1.80e-05 5.54e-06
150 1 2.62e+00 2.64e+00 1 total current 1.30e-05 3.67e-06
151 1 2.64e+00 2.65e+00 1 total current 1.60e-05 3.40e-06
152 1 2.65e+00 2.67e+00 1 total current 7.00e-06 3.35e-06
153 1 2.67e+00 2.69e+00 1 total current 1.00e-05 2.98e-06
154 1 2.69e+00 2.71e+00 1 total current 7.00e-06 3.35e-06
155 1 2.71e+00 2.72e+00 1 total current 1.20e-05 2.91e-06
156 1 2.72e+00 2.74e+00 1 total current 9.00e-06 2.33e-06
157 1 2.74e+00 2.76e+00 1 total current 1.00e-05 3.33e-06
158 1 2.76e+00 2.78e+00 1 total current 1.10e-05 3.14e-06
159 1 2.78e+00 2.79e+00 1 total current 1.00e-05 3.33e-06
160 1 2.79e+00 2.81e+00 1 total current 1.40e-05 4.76e-06
161 1 2.81e+00 2.83e+00 1 total current 8.00e-06 2.91e-06
162 1 2.83e+00 2.84e+00 1 total current 5.00e-06 2.69e-06
163 1 2.84e+00 2.86e+00 1 total current 6.00e-06 2.21e-06
164 1 2.86e+00 2.88e+00 1 total current 5.00e-06 1.67e-06
165 1 2.88e+00 2.90e+00 1 total current 4.00e-06 2.21e-06
166 1 2.90e+00 2.91e+00 1 total current 7.00e-06 2.13e-06
167 1 2.91e+00 2.93e+00 1 total current 6.00e-06 2.67e-06
168 1 2.93e+00 2.95e+00 1 total current 7.00e-06 2.13e-06
169 1 2.95e+00 2.97e+00 1 total current 5.00e-06 1.67e-06
170 1 2.97e+00 2.98e+00 1 total current 3.00e-06 1.53e-06
171 1 2.98e+00 3.00e+00 1 total current 6.00e-06 2.21e-06
172 1 3.00e+00 3.02e+00 1 total current 3.00e-06 1.53e-06
173 1 3.02e+00 3.04e+00 1 total current 1.00e-05 2.98e-06
174 1 3.04e+00 3.05e+00 1 total current 2.00e-06 1.33e-06
175 1 3.05e+00 3.07e+00 1 total current 1.00e-06 1.00e-06
176 1 3.07e+00 3.09e+00 1 total current 2.00e-06 1.33e-06
177 1 3.09e+00 3.11e+00 1 total current 0.00e+00 0.00e+00
178 1 3.11e+00 3.12e+00 1 total current 1.00e-06 1.00e-06
179 1 3.12e+00 3.14e+00 1 total current 0.00e+00 0.00e+00

80
tests/regression_tests/ncrystal/test.py Normal file
View file

@ -0,0 +1,80 @@
from math import pi
import numpy as np
import openmc
import openmc.lib
import pytest
from tests.testing_harness import PyAPITestHarness
pytestmark = pytest.mark.skipif(
not openmc.lib._ncrystal_enabled(),
reason="NCrystal materials are not enabled.")
def pencil_beam_model(cfg, E0, N):
"""Return an openmc.Model() object for a monoenergetic pencil
beam hitting a 1 mm sphere filled with the material defined by
the cfg string, and compute the angular distribution"""
# Material definition
m1 = openmc.Material.from_ncrystal(cfg)
materials = openmc.Materials([m1])
# Geometry definition
sample_sphere = openmc.Sphere(r=0.1)
outer_sphere = openmc.Sphere(r=100, boundary_type="vacuum")
cell1 = openmc.Cell(region=-sample_sphere, fill=m1)
cell2_region = +sample_sphere & -outer_sphere
cell2 = openmc.Cell(region=cell2_region, fill=None)
geometry = openmc.Geometry([cell1, cell2])
# Source definition
source = openmc.Source()
source.space = openmc.stats.Point((0, 0, -20))
source.angle = openmc.stats.Monodirectional(reference_uvw=(0, 0, 1))
source.energy = openmc.stats.Discrete([E0], [1.0])
# Execution settings
settings = openmc.Settings()
settings.source = source
settings.run_mode = "fixed source"
settings.batches = 10
settings.particles = N
# Tally definition
tally1 = openmc.Tally(name="angular distribution")
tally1.scores = ["current"]
filter1 = openmc.SurfaceFilter(sample_sphere)
filter2 = openmc.PolarFilter(np.linspace(0, pi, 180+1))
filter3 = openmc.CellFromFilter(cell1)
tally1.filters = [filter1, filter2, filter3]
tallies = openmc.Tallies([tally1])
return openmc.Model(geometry, materials, settings, tallies)
class NCrystalTest(PyAPITestHarness):
def _get_results(self):
"""Digest info in the statepoint and return as a string."""
# Read the statepoint file.
with openmc.StatePoint(self._sp_name) as sp:
tal = sp.get_tally(name='angular distribution')
df = tal.get_pandas_dataframe()
return df.to_string()
def test_ncrystal():
n_particles = 100000
T = 293.6 # K
E0 = 0.012 # eV
cfg = 'Al_sg225.ncmat'
test = pencil_beam_model(cfg, E0, n_particles)
harness = NCrystalTest('statepoint.10.h5', model=test)
harness.main()

0
tests/regression_tests/source_mcpl_file/__init__.py Normal file
View file

8
tests/regression_tests/source_mcpl_file/geometry.xml Normal file
View file

@ -0,0 +1,8 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

9
tests/regression_tests/source_mcpl_file/materials.xml Normal file
View file

@ -0,0 +1,9 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
</materials>

2
tests/regression_tests/source_mcpl_file/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
3.009416E-01 3.229999E-03

15
tests/regression_tests/source_mcpl_file/settings.xml Normal file
View file

@ -0,0 +1,15 @@
<?xml version="1.0"?>
<settings>
<state_point batches="10" />
<source_point mcpl="true" separate="true" />
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

101
tests/regression_tests/source_mcpl_file/test.py Normal file
View file

@ -0,0 +1,101 @@
#!/usr/bin/env python
import openmc.lib
import pytest
import glob
import os
from tests.testing_harness import *
pytestmark = pytest.mark.skipif(
not openmc.lib._mcpl_enabled(),
reason="MCPL is not enabled.")
settings1="""<?xml version="1.0"?>
<settings>
<state_point batches="10" />
<source_point mcpl="true" separate="true" />
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>
"""
settings2 = """<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<file>source.10.{}</file>
</source>
</settings>
"""
class SourceFileTestHarness(TestHarness):
def execute_test(self):
"""Run OpenMC with the appropriate arguments and check the outputs."""
try:
self._run_openmc()
self._test_output_created()
self._run_openmc_restart()
results = self._get_results()
self._write_results(results)
self._compare_results()
finally:
self._cleanup()
def update_results(self):
"""Update the results_true using the current version of OpenMC."""
try:
self._run_openmc()
self._test_output_created()
self._run_openmc_restart()
results = self._get_results()
self._write_results(results)
self._overwrite_results()
finally:
self._cleanup()
def _test_output_created(self):
"""Make sure statepoint and source files have been created."""
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))
assert len(statepoint) == 1, 'Either multiple or no statepoint files ' \
'exist.'
assert statepoint[0].endswith('h5'), \
'Statepoint file is not a HDF5 file.'
source = glob.glob(os.path.join(os.getcwd(), 'source.10.mcpl*'))
assert len(source) == 1, 'Either multiple or no source files exist.'
assert source[0].endswith('mcpl') or source[0].endswith('mcpl.gz'), \
'Source file is not a MCPL file.'
def _run_openmc_restart(self):
# Get the name of the source file.
source = glob.glob(os.path.join(os.getcwd(), 'source.10.*'))
# Write the new settings.xml file.
with open('settings.xml','w') as fh:
fh.write(settings2.format(source[0].split('.')[-1]))
# Run OpenMC.
self._run_openmc()
def _cleanup(self):
TestHarness._cleanup(self)
output = glob.glob(os.path.join(os.getcwd(), 'source.*'))
with open('settings.xml','w') as fh:
fh.write(settings1)
def test_source_file():
harness = SourceFileTestHarness('statepoint.10.h5')
harness.main()

5
tests/regression_tests/surface_tally/test.py
View file

@ -164,9 +164,8 @@ class SurfaceTallyTestHarness(PyAPITestHarness):
sp = openmc.StatePoint(self._sp_name)
# Extract the tally data as a Pandas DataFrame.
df = pd.DataFrame()
for t in sp.tallies.values():
df = df.append(t.get_pandas_dataframe(), ignore_index=True)
tally_dfs = [t.get_pandas_dataframe() for t in sp.tallies.values()]
df = pd.concat(tally_dfs, ignore_index=True)
# Extract the relevant data as a CSV string.
cols = ('mean', 'std. dev.')

34
tests/unit_tests/test_cell.py
View file

@ -57,25 +57,51 @@ def test_clone():
m = openmc.Material()
cyl = openmc.ZCylinder()
c = openmc.Cell(fill=m, region=-cyl)
c.temperature = 650.
# Check cloning with all optional params as the defaults
c2 = c.clone()
assert c2.id != c.id
assert c2.fill != c.fill
assert c2.region != c.region
assert c2.temperature == c.temperature
c3 = c.clone(clone_materials=False)
assert c3.id != c.id
assert c3.fill == c.fill
assert c3.region != c.region
assert c3.temperature == c.temperature
c4 = c.clone(clone_regions=False)
assert c4.id != c.id
assert c4.fill != c.fill
assert c4.region == c.region
assert c4.temperature == c.temperature
# Add optional properties to the original cell to ensure they're cloned successfully
c.temperature = 650.
c.translation = (1., 2., 3.)
c.rotation = (4., 5., 6.)
c.volume = 100
c5 = c.clone(clone_materials=False, clone_regions=False)
assert c5.id != c.id
assert c5.fill == c.fill
assert c5.region == c.region
assert c5.temperature == c.temperature
assert c5.volume == c.volume
assert all(c5.translation == c.translation)
assert all(c5.rotation == c.rotation)
# Mutate the original to ensure the changes are not seen in the clones
c.fill = openmc.Material()
c.region = +openmc.ZCylinder()
c.translation = (-1., -2., -3.)
c.rotation = (-4., -5., -6.)
c.temperature = 1
c.volume = 1
assert c5.fill != c.fill
assert c5.region != c.region
assert c5.temperature != c.temperature
assert c5.volume != c.volume
assert all(c5.translation != c.translation)
assert all(c5.rotation != c.rotation)
def test_temperature(cell_with_lattice):

12
tests/unit_tests/test_data_misc.py
View file

@ -101,12 +101,12 @@ def test_water_density():
assert dens(500.0, 3.0) == pytest.approx(1e-3/0.120241800e-2, 1e-6)
def test_gnd_name():
assert openmc.data.gnd_name(1, 1) == 'H1'
assert openmc.data.gnd_name(40, 90) == ('Zr90')
assert openmc.data.gnd_name(95, 242, 0) == ('Am242')
assert openmc.data.gnd_name(95, 242, 1) == ('Am242_m1')
assert openmc.data.gnd_name(95, 242, 10) == ('Am242_m10')
def test_gnds_name():
assert openmc.data.gnds_name(1, 1) == 'H1'
assert openmc.data.gnds_name(40, 90) == ('Zr90')
assert openmc.data.gnds_name(95, 242, 0) == ('Am242')
assert openmc.data.gnds_name(95, 242, 1) == ('Am242_m1')
assert openmc.data.gnds_name(95, 242, 10) == ('Am242_m10')
def test_isotopes():

34
tests/unit_tests/test_geometry.py
View file

@ -1,4 +1,5 @@
import xml.etree.ElementTree as ET
from pathlib import Path
import numpy as np
import openmc
@ -159,12 +160,13 @@ def test_get_by_name():
m2 = openmc.Material(name='Zirconium')
m2.add_element('Zr', 1.0)
c1 = openmc.Cell(fill=m1, name='cell1')
s1 = openmc.Sphere(name='surface1')
c1 = openmc.Cell(fill=m1, region=-s1, name='cell1')
u1 = openmc.Universe(name='Zircaloy universe', cells=[c1])
cyl = openmc.ZCylinder()
c2 = openmc.Cell(fill=u1, region=-cyl, name='cell2')
c3 = openmc.Cell(fill=m2, region=+cyl, name='Cell3')
s2 = openmc.ZCylinder(name='surface2')
c2 = openmc.Cell(fill=u1, region=-s2, name='cell2')
c3 = openmc.Cell(fill=m2, region=+s2, name='Cell3')
root = openmc.Universe(name='root Universe', cells=[c2, c3])
geom = openmc.Geometry(root)
@ -177,6 +179,13 @@ def test_get_by_name():
mats = geom.get_materials_by_name('zirconium', True, True)
assert not mats
surfaces = set(geom.get_surfaces_by_name('surface'))
assert not surfaces ^ {s1, s2}
surfaces = set(geom.get_surfaces_by_name('Surface2', False, True))
assert not surfaces ^ {s2}
surfaces = geom.get_surfaces_by_name('Surface2', True, True)
assert not surfaces
cells = set(geom.get_cells_by_name('cell'))
assert not cells ^ {c1, c2, c3}
cells = set(geom.get_cells_by_name('cell', True))
@ -274,7 +283,22 @@ def test_from_xml(run_in_tmpdir, mixed_lattice_model):
# Export model
mixed_lattice_model.export_to_xml()
# Import geometry
mats_from_xml = openmc.Materials.from_xml('materials.xml')
# checking string a Path are both acceptable
for path in ['geometry.xml', Path('geometry.xml')]:
for materials in [mats_from_xml, 'materials.xml']:
# Import geometry from file
geom = openmc.Geometry.from_xml(path=path, materials=materials)
assert isinstance(geom, openmc.Geometry)
ll, ur = geom.bounding_box
assert ll == pytest.approx((-6.0, -6.0, -np.inf))
assert ur == pytest.approx((6.0, 6.0, np.inf))
with pytest.raises(TypeError) as excinfo:
geom = openmc.Geometry.from_xml(path='geometry.xml', materials=None)
assert 'Unable to set "materials" to "None"' in str(excinfo.value)
# checking that the default args also work
geom = openmc.Geometry.from_xml()
assert isinstance(geom, openmc.Geometry)
ll, ur = geom.bounding_box

4
tests/unit_tests/test_mesh_to_vtk.py
View file

@ -77,8 +77,8 @@ def test_write_data_to_vtk_size_mismatch(mesh):
# by regex. These are needed to make the test string match the error message
# string when using the match argument as that uses regular expression
expected_error_msg = (
f"The size of the dataset 'label' \({len(data)}\) should be equal to "
f"the number of mesh cells \({mesh.num_mesh_cells}\)"
fr"The size of the dataset 'label' \({len(data)}\) should be equal to "
fr"the number of mesh cells \({mesh.num_mesh_cells}\)"
)
with pytest.raises(ValueError, match=expected_error_msg):
mesh.write_data_to_vtk(filename="out.vtk", datasets={"label": data})

38
tests/unit_tests/test_model.py
View file

@ -1,5 +1,6 @@
from math import pi
from pathlib import Path
import os
import numpy as np
import pytest
@ -529,3 +530,40 @@ def test_calc_volumes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
assert openmc.lib.materials[3].volume == mats[2].volume
test_model.finalize_lib()
def test_model_xml(run_in_tmpdir):
# load a model from examples
pwr_model = openmc.examples.pwr_core()
# export to separate XMLs manually
pwr_model.settings.export_to_xml('settings_ref.xml')
pwr_model.materials.export_to_xml('materials_ref.xml')
pwr_model.geometry.export_to_xml('geometry_ref.xml')
# now write and read a model.xml file
pwr_model.export_to_model_xml()
new_model = openmc.Model.from_model_xml()
# make sure we can also export this again to separate
# XML files
new_model.export_to_xml()
def test_single_xml_exec(run_in_tmpdir):
pincell_model = openmc.examples.pwr_pin_cell()
pincell_model.export_to_model_xml('pwr_pincell.xml')
openmc.run(path_input='pwr_pincell.xml')
with pytest.raises(RuntimeError, match='ex-em-ell.xml'):
openmc.run(path_input='ex-em-ell.xml')
# test that a file in a different directory can be used
os.mkdir('inputs')
pincell_model.export_to_model_xml('./inputs/pincell.xml')
openmc.run(path_input='./inputs/pincell.xml')
with pytest.raises(RuntimeError, match='input_dir'):
openmc.run(path_input='input_dir/pincell.xml')

27
tests/unit_tests/test_plotter.py Normal file
View file

@ -0,0 +1,27 @@
import openmc
import numpy as np
def test_calculate_cexs_elem_mat_sab():
"""Checks that sab cross sections are included in the
_calculate_cexs_elem_mat method and have the correct shape"""
mat_1 = openmc.Material()
mat_1.add_element("H", 4.0, "ao")
mat_1.add_element("O", 4.0, "ao")
mat_1.add_element("C", 4.0, "ao")
mat_1.add_s_alpha_beta("c_C6H6")
mat_1.set_density("g/cm3", 0.865)
energy_grid, data = openmc.plotter._calculate_cexs_elem_mat(
mat_1,
["inelastic"],
sab_name="c_C6H6",
)
assert isinstance(energy_grid, np.ndarray)
assert isinstance(data, np.ndarray)
assert len(energy_grid) > 1
assert len(data) == 1
assert len(data[0]) == len(energy_grid)

6
tests/unit_tests/test_settings.py
View file

@ -17,7 +17,7 @@ def test_export_to_xml(run_in_tmpdir):
s.output = {'summary': True, 'tallies': False, 'path': 'here'}
s.verbosity = 7
s.sourcepoint = {'batches': [50, 150, 500, 1000], 'separate': True,
'write': True, 'overwrite': True}
'write': True, 'overwrite': True, 'mcpl': True}
s.statepoint = {'batches': [50, 150, 500, 1000]}
s.surf_source_read = {'path': 'surface_source_1.h5'}
s.surf_source_write = {'surface_ids': [2], 'max_particles': 200}
@ -49,6 +49,7 @@ def test_export_to_xml(run_in_tmpdir):
domains=[openmc.Cell()], samples=1000, lower_left=(-10., -10., -10.),
upper_right = (10., 10., 10.))
s.create_fission_neutrons = True
s.create_delayed_neutrons = False
s.log_grid_bins = 2000
s.photon_transport = False
s.electron_treatment = 'led'
@ -75,7 +76,7 @@ def test_export_to_xml(run_in_tmpdir):
assert s.output == {'summary': True, 'tallies': False, 'path': 'here'}
assert s.verbosity == 7
assert s.sourcepoint == {'batches': [50, 150, 500, 1000], 'separate': True,
'write': True, 'overwrite': True}
'write': True, 'overwrite': True, 'mcpl': True}
assert s.statepoint == {'batches': [50, 150, 500, 1000]}
assert s.surf_source_read == {'path': 'surface_source_1.h5'}
assert s.surf_source_write == {'surface_ids': [2], 'max_particles': 200}
@ -107,6 +108,7 @@ def test_export_to_xml(run_in_tmpdir):
'energy_min': 1.0, 'energy_max': 1000.0,
'nuclides': ['U235', 'U238', 'Pu239']}
assert s.create_fission_neutrons
assert not s.create_delayed_neutrons
assert s.log_grid_bins == 2000
assert not s.photon_transport
assert s.electron_treatment == 'led'

75
tests/unit_tests/test_triggers.py Normal file
View file

@ -0,0 +1,75 @@
import openmc
def test_tally_trigger(run_in_tmpdir):
pincell = openmc.examples.pwr_pin_cell()
# create a tally filter on the materials
mat_filter = openmc.MaterialFilter(pincell.materials)
# create a tally with triggers applied
tally = openmc.Tally()
tally.filters = [mat_filter]
tally.scores = ['scatter']
trigger = openmc.Trigger('rel_err', 0.05)
trigger.scores = ['scatter']
tally.triggers = [trigger]
pincell.tallies = [tally]
pincell.settings.trigger_active = True
pincell.settings.trigger_max_batches = 100
pincell.settings.trigger_batch_interval = 5
sp_file = pincell.run()
with openmc.StatePoint(sp_file) as sp:
expected_realizations = sp.n_realizations
# adding other scores to the tally should not change the
# number of batches required to satisfy the trigger
tally.scores = ['total', 'absorption', 'scatter']
sp_file = pincell.run()
with openmc.StatePoint(sp_file) as sp:
realizations = sp.n_realizations
assert realizations == expected_realizations
def test_tally_trigger_null_score(run_in_tmpdir):
pincell = openmc.examples.pwr_pin_cell()
# create a tally filter on the materials
mat_filter = openmc.MaterialFilter(pincell.materials)
# apply a tally with a score that be tallied in this model
tally = openmc.Tally()
tally.filters = [mat_filter]
tally.scores = ['pair-production']
trigger = openmc.Trigger('rel_err', 0.05)
trigger.scores = ['pair-production']
tally.triggers = [trigger]
pincell.tallies = [tally]
pincell.settings.trigger_active = True
pincell.settings.trigger_max_batches = 50
pincell.settings.trigger_batch_interval = 5
sp_file = pincell.run()
with openmc.StatePoint(sp_file) as sp:
# verify that the tally mean is zero
tally_out = sp.get_tally(id=tally.id)
assert all(tally_out.mean == 0.0)
# we expect that this simulation will run
# up to the max allowed batches
total_batches = sp.n_realizations + sp.n_inactive
assert total_batches == pincell.settings.trigger_max_batches

7
tools/ci/gha-install-mcpl.sh Executable file
View file

@ -0,0 +1,7 @@
#!/bin/bash
set -ex
cd $HOME
git clone https://github.com/mctools/mcpl
cd mcpl
mkdir build && cd build
cmake .. && make 2>/dev/null && sudo make install

46
tools/ci/gha-install-ncrystal.sh Executable file
View file

@ -0,0 +1,46 @@
#!/bin/bash
set -ex
cd $HOME
#Use the NCrystal develop branch (in the near future we can move this to master):
git clone https://github.com/mctools/ncrystal --branch develop --single-branch --depth 1 ncrystal_src
SRC_DIR="$PWD/ncrystal_src"
BLD_DIR="$PWD/ncrystal_bld"
INST_DIR="$PWD/ncrystal_inst"
PYTHON=$(which python3)
CPU_COUNT=1
mkdir "$BLD_DIR"
cd ncrystal_bld
cmake \
"${SRC_DIR}" \
-DBUILD_SHARED_LIBS=ON \
-DNCRYSTAL_NOTOUCH_CMAKE_BUILD_TYPE=ON \
-DNCRYSTAL_MODIFY_RPATH=OFF \
-DCMAKE_BUILD_TYPE=Release \
-DNCRYSTAL_ENABLE_EXAMPLES=OFF \
-DNCRYSTAL_ENABLE_SETUPSH=OFF \
-DNCRYSTAL_ENABLE_DATA=EMBED \
-DCMAKE_INSTALL_PREFIX="${INST_DIR}" \
-DPython3_EXECUTABLE="$PYTHON"
make -j${CPU_COUNT:-1}
make install
#Note: There is no "make test" or "make ctest" functionality for NCrystal
# yet. If it appears in the future, we should add it here.
# Output the configuration to the log
"${INST_DIR}/bin/ncrystal-config" --setup
# Change environmental variables
eval $( "${INST_DIR}/bin/ncrystal-config" --setup )
# Check installation worked
nctool --test

14
tools/ci/gha-install.py
View file

@ -19,14 +19,14 @@ def which(program):
return None
def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False):
def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False, ncrystal=False):
# Create build directory and change to it
shutil.rmtree('build', ignore_errors=True)
os.mkdir('build')
os.chdir('build')
# Build in debug mode by default
cmake_cmd = ['cmake', '-DCMAKE_BUILD_TYPE=Debug']
# Build in debug mode by default with support for MCPL
cmake_cmd = ['cmake', '-DCMAKE_BUILD_TYPE=Debug', '-DOPENMC_USE_MCPL=on']
# Turn off OpenMP if specified
if not omp:
@ -54,6 +54,11 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False):
libmesh_path = os.environ.get('HOME') + '/LIBMESH'
cmake_cmd.append('-DCMAKE_PREFIX_PATH=' + libmesh_path)
if ncrystal:
cmake_cmd.append('-DOPENMC_USE_NCRYSTAL=ON')
ncrystal_cmake_path = os.environ.get('HOME') + '/ncrystal_inst/lib/cmake'
cmake_cmd.append(f'-DCMAKE_PREFIX_PATH={ncrystal_cmake_path}')
# Build in coverage mode for coverage testing
cmake_cmd.append('-DOPENMC_ENABLE_COVERAGE=on')
@ -70,10 +75,11 @@ def main():
mpi = (os.environ.get('MPI') == 'y')
phdf5 = (os.environ.get('PHDF5') == 'y')
dagmc = (os.environ.get('DAGMC') == 'y')
ncrystal = (os.environ.get('NCRYSTAL') == 'y')
libmesh = (os.environ.get('LIBMESH') == 'y')
# Build and install
install(omp, mpi, phdf5, dagmc, libmesh)
install(omp, mpi, phdf5, dagmc, libmesh, ncrystal)
if __name__ == '__main__':
main()

8
tools/ci/gha-install.sh
View file

@ -17,6 +17,11 @@ if [[ $DAGMC = 'y' ]]; then
./tools/ci/gha-install-dagmc.sh
fi
# Install NCrystal if needed
if [[ $NCRYSTAL = 'y' ]]; then
./tools/ci/gha-install-ncrystal.sh
fi
# Install vectfit for WMP generation if needed
if [[ $VECTFIT = 'y' ]]; then
./tools/ci/gha-install-vectfit.sh
@ -27,6 +32,9 @@ if [[ $LIBMESH = 'y' ]]; then
./tools/ci/gha-install-libmesh.sh
fi
# Install MCPL
./tools/ci/gha-install-mcpl.sh
# For MPI configurations, make sure mpi4py and h5py are built against the
# correct version of MPI
if [[ $MPI == 'y' ]]; then

7
tools/ci/gha-script.sh
View file

@ -14,5 +14,12 @@ if [[ $EVENT == 'y' ]]; then
args="${args} --event "
fi
# Check NCrystal installation
if [[ $NCRYSTAL = 'y' ]]; then
# Change environmental variables
eval $( "${HOME}/ncrystal_inst/bin/ncrystal-config" --setup )
nctool --test
fi
# Run regression and unit tests
pytest --cov=openmc -v $args tests