Merge branch 'openmc-dev:develop' into universe_cloning

This commit is contained in:
Josh May 2023-01-16 18:43:34 -08:00 committed by GitHub
commit d9ce50dc9f
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
65 changed files with 1234 additions and 168 deletions

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

View file

@ -37,6 +37,7 @@ option(OPENMC_USE_DAGMC "Enable support for DAGMC (CAD) geometry"
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")
@ -99,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
#===============================================================================
@ -338,6 +348,7 @@ list(APPEND libopenmc_SOURCES
src/message_passing.cpp
src/mgxs.cpp
src/mgxs_interface.cpp
src/ncrystal_interface.cpp
src/nuclide.cpp
src/output.cpp
src/particle.cpp
@ -498,6 +509,11 @@ if (OPENMC_USE_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
#===============================================================================

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

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

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

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

View file

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

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)

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/

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

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.

View file

@ -14,7 +14,8 @@ namespace openmc {
inline bool dir_exists(const std::string& path)
{
struct stat s;
if (stat(path.c_str(), &s) != 0) return false;
if (stat(path.c_str(), &s) != 0)
return false;
return s.st_mode & S_IFDIR;
}

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]

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

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

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

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

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)

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

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

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

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

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

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)

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.

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)

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

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:
@ -254,16 +256,20 @@ class Geometry:
raise ValueError('Error determining root universe.')
@classmethod
def from_xml(cls, path='geometry.xml', materials=None):
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 : str, optional
path : PathLike, optional
Path to geometry XML file
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.
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
-------
@ -271,10 +277,13 @@ class Geometry:
Geometry object
"""
# Create dictionary to easily look up materials
if materials is None:
filename = Path(path).parent / 'materials.xml'
materials = openmc.Materials.from_xml(str(filename))
# 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()
@ -525,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.

View file

@ -42,6 +42,9 @@ 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

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))
@ -1566,4 +1650,4 @@ class Materials(cv.CheckedList):
tree = ET.parse(path)
root = tree.getroot()
return cls.from_xml_element(root)
return cls.from_xml_element(root)

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

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)

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:

View file

@ -225,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
@ -296,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
@ -459,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
@ -866,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)
@ -1207,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")
@ -1532,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:
@ -1630,6 +1655,7 @@ class Settings:
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)
@ -1722,6 +1748,7 @@ class Settings:
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)

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

View file

@ -91,7 +91,8 @@ Library::Library(pugi::xml_node node, const std::string& directory)
// Non-member functions
//==============================================================================
void read_cross_sections_xml() {
void read_cross_sections_xml()
{
pugi::xml_document doc;
std::string filename = settings::path_input + "materials.xml";
// Check if materials.xml exists

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

View file

@ -40,7 +40,8 @@ void update_universe_cell_count(int32_t a, int32_t b)
}
}
void read_geometry_xml() {
void read_geometry_xml()
{
// Display output message
write_message("Reading geometry XML file...", 5);

View file

@ -106,7 +106,8 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
openmc::openmc_set_seed(DEFAULT_SEED);
// Read XML input files
if (!read_model_xml()) read_separate_xml_files();
if (!read_model_xml())
read_separate_xml_files();
// Write some initial output under the header if needed
initial_output();
@ -299,7 +300,8 @@ int parse_command_line(int argc, char* argv[])
return 0;
}
bool read_model_xml() {
bool read_model_xml()
{
std::string model_filename =
settings::path_input.empty() ? "." : settings::path_input;
@ -314,7 +316,8 @@ bool read_model_xml() {
model_filename += "/model.xml";
// if this file doesn't exist, stop here
if (!file_exists(model_filename)) return false;
if (!file_exists(model_filename))
return false;
// try to process the path input as an XML file
pugi::xml_document doc;
@ -327,7 +330,7 @@ bool read_model_xml() {
// Read settings
if (!check_for_node(root, "settings")) {
fatal_error("No <settings> node present in the model.xml file.");
fatal_error("No <settings> node present in the model.xml file.");
}
auto settings_root = root.child("settings");
@ -343,13 +346,15 @@ bool read_model_xml() {
title();
}
write_message(fmt::format("Reading model XML file '{}' ...", model_filename), 5);
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"};
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 "
@ -384,7 +389,7 @@ bool read_model_xml() {
if (check_for_node(root, "tallies"))
read_tallies_xml(root.child("tallies"));
// Initialize distribcell_filters
// Initialize distribcell_filters
prepare_distribcell();
if (check_for_node(root, "plots"))
@ -416,7 +421,8 @@ void read_separate_xml_files()
read_plots_xml();
}
void initial_output() {
void initial_output()
{
// write initial output
if (settings::run_mode == RunMode::PLOTTING) {
// Read plots.xml if it exists

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);
}
}
// ======================================================================
@ -1264,7 +1283,8 @@ double density_effect(const vector<double>& f, const vector<double>& e_b_sq,
return delta - w_sq * (1.0 - beta_sq);
}
void read_materials_xml() {
void read_materials_xml()
{
write_message("Reading materials XML file...", 5);
pugi::xml_document doc;

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

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

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

View file

@ -124,7 +124,8 @@ extern "C" int openmc_plot_geometry()
return 0;
}
void read_plots_xml() {
void read_plots_xml()
{
// Check if plots.xml exists; this is only necessary when the plot runmode is
// initiated. Otherwise, we want to read plots.xml because it may be called
// later via the API. In that case, its ok for a plots.xml to not exist

View file

@ -44,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};
@ -214,7 +215,8 @@ void get_run_parameters(pugi::xml_node node_base)
}
}
void read_settings_xml() {
void read_settings_xml()
{
using namespace settings;
using namespace pugi;
// Check if settings.xml exists
@ -872,6 +874,12 @@ void read_settings_xml(pugi::xml_node root)
}
}
// 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")) {

View file

@ -705,7 +705,8 @@ std::string Tally::nuclide_name(int nuclide_idx) const
// Non-member functions
//==============================================================================
void read_tallies_xml() {
void read_tallies_xml()
{
// Check if tallies.xml exists. If not, just return since it is optional
std::string filename = settings::path_input + "tallies.xml";
if (!file_exists(filename))

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

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

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

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>

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

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

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

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

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

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

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

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)

View file

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

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

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

View file

@ -19,7 +19,7 @@ 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')
@ -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()

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

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