mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Merge pull request #2 from paulromano/mixed_ncrystal_pr
Refactor interface between NCrystal and OpenMC
This commit is contained in:
Thomas Kittelmann
GitHub
commit
fe4a1de57d
14 changed files with 297 additions and 171 deletions
|
|
@ -348,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
|
||||
|
|
|
|||
|
|
@ -182,19 +182,22 @@ been selected. There are three methods available:
|
|||
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`_.
|
||||
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.
|
||||
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
|
||||
|
|
|
|||
|
|
@ -274,9 +274,10 @@ Prerequisites
|
|||
* 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>`_
|
||||
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 ..
|
||||
|
|
@ -374,8 +375,8 @@ OPENMC_USE_DAGMC
|
|||
(Default: off)
|
||||
|
||||
OPENMC_USE_NCRYSTAL
|
||||
Turns on support for NCrystal materials. NCrystal must be
|
||||
`installed <https://github.com/mctools/ncrystal/wiki/Get-NCrystal>`_ and
|
||||
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)
|
||||
|
||||
|
|
|
|||
|
|
@ -105,12 +105,14 @@ you would need to add hydrogen and oxygen to a material and then assign the
|
|||
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.::
|
||||
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')
|
||||
|
||||
|
|
@ -125,12 +127,12 @@ defines a material containing polycrystalline alumnium,
|
|||
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
|
||||
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
|
||||
.. 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
|
||||
|
|
|
|||
|
|
@ -12,13 +12,10 @@
|
|||
#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"
|
||||
|
||||
#ifdef NCRYSTAL
|
||||
#include "NCrystal/NCrystal.hh"
|
||||
#endif
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -155,25 +152,17 @@ public:
|
|||
//! \return Temperature in [K]
|
||||
double temperature() const;
|
||||
|
||||
#ifdef NCRYSTAL
|
||||
//! Get pointer to NCrystal material object
|
||||
//! \return Pointer to NCrystal material object
|
||||
std::shared_ptr<const NCrystal::ProcImpl::Process> ncrystal_mat() const
|
||||
{
|
||||
return ncrystal_mat_;
|
||||
};
|
||||
#endif
|
||||
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
|
||||
#ifdef NCRYSTAL
|
||||
std::string ncrystal_cfg_; //!< NCrystal configuration string
|
||||
std::shared_ptr<const NCrystal::ProcImpl::Process> ncrystal_mat_;
|
||||
#endif
|
||||
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]
|
||||
|
|
|
|||
91
include/openmc/ncrystal_interface.h
Normal file
91
include/openmc/ncrystal_interface.h
Normal file
|
|
@ -0,0 +1,91 @@
|
|||
#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;
|
||||
|
||||
//==============================================================================
|
||||
// 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
|
||||
|
|
@ -8,10 +8,6 @@
|
|||
#include "openmc/reaction.h"
|
||||
#include "openmc/vector.h"
|
||||
|
||||
#ifdef NCRYSTAL
|
||||
#include "NCrystal/NCRNG.hh"
|
||||
#endif
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -100,31 +96,12 @@ 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);
|
||||
|
||||
#ifdef NCRYSTAL
|
||||
//==============================================================================
|
||||
// NCrystal wrapper class for the OpenMC random number generator
|
||||
//==============================================================================
|
||||
|
||||
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
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_PHYSICS_H
|
||||
|
|
|
|||
|
|
@ -20,8 +20,6 @@ namespace openmc {
|
|||
// Global variable declarations
|
||||
//==============================================================================
|
||||
|
||||
extern "C" const bool NCRYSTAL_ENABLED;
|
||||
|
||||
namespace settings {
|
||||
|
||||
// Boolean flags
|
||||
|
|
@ -93,6 +91,8 @@ extern int n_log_bins; //!< number of bins for logarithmic energy grid
|
|||
extern int n_batches; //!< number of (inactive+active) batches
|
||||
extern int n_max_batches; //!< Maximum number of batches
|
||||
extern int max_tracks; //!< Maximum number of particle tracks written to file
|
||||
extern double
|
||||
ncrystal_max_energy; //!< Energy in eV to switch between NCrystal and ENDF
|
||||
extern ResScatMethod res_scat_method; //!< resonance upscattering method
|
||||
extern double res_scat_energy_min; //!< Min energy in [eV] for res. upscattering
|
||||
extern double res_scat_energy_max; //!< Max energy in [eV] for res. upscattering
|
||||
|
|
@ -124,10 +124,6 @@ extern int trigger_batch_interval; //!< Batch interval for triggers
|
|||
extern "C" int verbosity; //!< How verbose to make output
|
||||
extern double weight_cutoff; //!< Weight cutoff for Russian roulette
|
||||
extern double weight_survive; //!< Survival weight after Russian roulette
|
||||
#ifdef NCRYSTAL
|
||||
extern double
|
||||
ncrystal_max_energy; //!< Energy in eV to switch between NCrystal and ENDF
|
||||
#endif
|
||||
} // namespace settings
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -99,6 +99,10 @@ class Material(IDManagerMixin):
|
|||
[decay/sec].
|
||||
|
||||
.. versionadded:: 0.13.2
|
||||
ncrystal_cfg : str
|
||||
NCrystal configuration string
|
||||
|
||||
.. versionadded:: 0.13.3
|
||||
|
||||
"""
|
||||
|
||||
|
|
@ -141,7 +145,8 @@ class Material(IDManagerMixin):
|
|||
|
||||
string += '{: <16}\n'.format('\tS(a,b) Tables')
|
||||
|
||||
string += '{: <16}=\t{}\n'.format('\tNCrystal conf', self._ncrystal_cfg)
|
||||
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)
|
||||
|
|
@ -222,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:
|
||||
|
|
@ -335,9 +344,9 @@ class Material(IDManagerMixin):
|
|||
return material
|
||||
|
||||
@classmethod
|
||||
def from_ncrystal(cls, cfg, material_id=None, name=''):
|
||||
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
|
||||
|
|
@ -347,12 +356,8 @@ class Material(IDManagerMixin):
|
|||
----------
|
||||
cfg : str
|
||||
NCrystal configuration string
|
||||
material_id : int, optional
|
||||
Unique identifier for the material. If not specified, an identifier will
|
||||
automatically be assigned.
|
||||
name : str, optional
|
||||
Name of the material. If not specified, the name will be the empty
|
||||
string.
|
||||
**kwargs
|
||||
Keyword arguments passed to :class:`openmc.Material`
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -371,26 +376,20 @@ class Material(IDManagerMixin):
|
|||
#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.get(Z)
|
||||
if not elem_name:
|
||||
raise ValueError( f'Element with Z={Z} is not known' )
|
||||
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)
|
||||
flat_compos = nc_mat.getFlattenedComposition(
|
||||
preferNaturalElements=True, naturalAbundProvider=openmc_natabund)
|
||||
|
||||
# Create the Material
|
||||
material = cls(material_id=material_id,
|
||||
name=name,
|
||||
temperature=nc_mat.getTemperature())
|
||||
material = cls(temperature=nc_mat.getTemperature(), **kwargs)
|
||||
|
||||
for Z, A_vals in flat_compos:
|
||||
elemname = openmc.data.ATOMIC_SYMBOL.get(Z)
|
||||
if not elemname:
|
||||
raise ValueError(f'Element with Z={Z} is not known')
|
||||
elemname = openmc.data.ATOMIC_SYMBOL[Z]
|
||||
for A, frac in A_vals:
|
||||
if A:
|
||||
material.add_nuclide(f'{elemname}{A}', frac)
|
||||
|
|
@ -1067,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].
|
||||
|
|
@ -1101,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
|
||||
|
||||
|
|
@ -1651,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)
|
||||
|
|
|
|||
|
|
@ -60,13 +60,12 @@ Material::Material(pugi::xml_node node)
|
|||
name_ = get_node_value(node, "name");
|
||||
}
|
||||
|
||||
#ifdef NCRYSTAL
|
||||
if (check_for_node(node, "cfg")) {
|
||||
ncrystal_cfg_ = get_node_value(node, "cfg");
|
||||
write_message(5, "NCrystal config string for material #{}: '{}'", this->id(), ncrystal_cfg_);
|
||||
ncrystal_mat_ = NCrystal::FactImpl::createScatter(ncrystal_cfg_);
|
||||
auto cfg = get_node_value(node, "cfg");
|
||||
write_message(
|
||||
5, "NCrystal config string for material #{}: '{}'", this->id(), cfg);
|
||||
ncrystal_mat_ = NCrystalMat(cfg);
|
||||
}
|
||||
#endif
|
||||
|
||||
if (check_for_node(node, "depletable")) {
|
||||
depletable_ = get_node_value_bool(node, "depletable");
|
||||
|
|
@ -800,19 +799,11 @@ void Material::calculate_neutron_xs(Particle& p) const
|
|||
// Initialize position in i_sab_nuclides
|
||||
int j = 0;
|
||||
|
||||
#ifdef NCRYSTAL
|
||||
double ncrystal_xs = -1;
|
||||
|
||||
// Calculate NCrystal cross section
|
||||
double ncrystal_xs = -1.0;
|
||||
if (ncrystal_mat_ && p.E() < settings::ncrystal_max_energy) {
|
||||
// Calculate scattering XS per atom with NCrystal, only once per material
|
||||
NCrystal::CachePtr dummy_cache;
|
||||
auto nc_energy = NCrystal::NeutronEnergy {p.E()};
|
||||
ncrystal_xs =
|
||||
ncrystal_mat_
|
||||
->crossSection(dummy_cache, nc_energy, {p.u().x, p.u().y, p.u().z})
|
||||
.get();
|
||||
ncrystal_xs = ncrystal_mat_.xs(p);
|
||||
}
|
||||
#endif
|
||||
|
||||
// Add contribution from each nuclide in material
|
||||
for (int i = 0; i < nuclide_.size(); ++i) {
|
||||
|
|
@ -852,26 +843,16 @@ void Material::calculate_neutron_xs(Particle& p) const
|
|||
int i_nuclide = nuclide_[i];
|
||||
|
||||
// Calculate microscopic cross section for this nuclide
|
||||
#ifdef NCRYSTAL
|
||||
auto& micro {p.neutron_xs(i_nuclide)};
|
||||
#else
|
||||
const auto& micro {p.neutron_xs(i_nuclide)};
|
||||
#endif
|
||||
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);
|
||||
#ifdef NCRYSTAL
|
||||
|
||||
// If NCrystal is being used, update micro cross section cache
|
||||
if (ncrystal_xs >= 0.0) {
|
||||
if (micro.thermal > 0 || micro.thermal_elastic > 0) {
|
||||
fatal_error("S(a,b) treatment and NCrystal are not compatible.");
|
||||
}
|
||||
data::nuclides[i_nuclide]->calculate_elastic_xs(p);
|
||||
// remove free atom cross section
|
||||
// and replace it by scattering cross section per atom from NCrystal
|
||||
micro.total = micro.total - micro.elastic + ncrystal_xs;
|
||||
micro.elastic = ncrystal_xs;
|
||||
ncrystal_update_micro(ncrystal_xs, micro);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
// ======================================================================
|
||||
|
|
|
|||
108
src/ncrystal_interface.cpp
Normal file
108
src/ncrystal_interface.cpp
Normal 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
|
||||
|
|
@ -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"
|
||||
|
|
@ -138,32 +139,15 @@ void sample_neutron_reaction(Particle& p)
|
|||
if (!p.alive())
|
||||
return;
|
||||
|
||||
// Sample a scattering reaction and determine the secondary energy of the
|
||||
// exiting neutron
|
||||
|
||||
#ifdef NCRYSTAL
|
||||
if (model::materials[p.material()]->ncrystal_mat() &&
|
||||
p.E() < settings::ncrystal_max_energy) {
|
||||
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 =
|
||||
model::materials[p.material()]->ncrystal_mat()->sampleScatter(
|
||||
dummy_cache, rng, nc_energy, {p.u().x, p.u().y, p.u().z});
|
||||
p.E_last() = p.E();
|
||||
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;
|
||||
// Sample a scattering reaction and determine the secondary energy of the
|
||||
// exiting neutron
|
||||
const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat();
|
||||
if (ncrystal_mat && p.E() < settings::ncrystal_max_energy) {
|
||||
ncrystal_mat.scatter(p);
|
||||
} else {
|
||||
scatter(p, i_nuclide);
|
||||
}
|
||||
#else
|
||||
scatter(p, i_nuclide);
|
||||
#endif
|
||||
|
||||
// Advance URR seed stream 'N' times after energy changes
|
||||
if (p.E() != p.E_last()) {
|
||||
p.stream() = STREAM_URR_PTABLE;
|
||||
|
|
|
|||
|
|
@ -37,13 +37,6 @@ namespace openmc {
|
|||
// Global variables
|
||||
//==============================================================================
|
||||
|
||||
#ifdef NCRYSTAL
|
||||
const bool NCRYSTAL_ENABLED = true;
|
||||
#else
|
||||
const bool NCRYSTAL_ENABLED = false;
|
||||
#endif
|
||||
|
||||
|
||||
namespace settings {
|
||||
|
||||
// Default values for boolean flags
|
||||
|
|
@ -106,6 +99,7 @@ int n_batches;
|
|||
int n_max_batches;
|
||||
int max_splits {1000};
|
||||
int max_tracks {1000};
|
||||
double ncrystal_max_energy {5.0};
|
||||
ResScatMethod res_scat_method {ResScatMethod::rvs};
|
||||
double res_scat_energy_min {0.01};
|
||||
double res_scat_energy_max {1000.0};
|
||||
|
|
@ -128,10 +122,6 @@ int verbosity {7};
|
|||
double weight_cutoff {0.25};
|
||||
double weight_survive {1.0};
|
||||
|
||||
#ifdef NCRYSTAL
|
||||
double ncrystal_max_energy {5.0};
|
||||
#endif
|
||||
|
||||
} // namespace settings
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -1,3 +1,5 @@
|
|||
from math import pi
|
||||
|
||||
import numpy as np
|
||||
import openmc
|
||||
import openmc.lib
|
||||
|
|
@ -9,6 +11,7 @@ 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
|
||||
|
|
@ -24,7 +27,7 @@ def pencil_beam_model(cfg, E0, N):
|
|||
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_region = +sample_sphere & -outer_sphere
|
||||
cell2 = openmc.Cell(region=cell2_region, fill=None)
|
||||
geometry = openmc.Geometry([cell1, cell2])
|
||||
|
||||
|
|
@ -48,7 +51,7 @@ def pencil_beam_model(cfg, E0, N):
|
|||
tally1 = openmc.Tally(name="angular distribution")
|
||||
tally1.scores = ["current"]
|
||||
filter1 = openmc.SurfaceFilter(sample_sphere)
|
||||
filter2 = openmc.PolarFilter(np.linspace(0, np.pi, 180+1))
|
||||
filter2 = openmc.PolarFilter(np.linspace(0, pi, 180+1))
|
||||
filter3 = openmc.CellFromFilter(cell1)
|
||||
tally1.filters = [filter1, filter2, filter3]
|
||||
tallies = openmc.Tallies([tally1])
|
||||
|
|
@ -66,11 +69,12 @@ class NCrystalTest(PyAPITestHarness):
|
|||
df = tal.get_pandas_dataframe()
|
||||
return df.to_string()
|
||||
|
||||
|
||||
def test_ncrystal():
|
||||
NParticles = 100000
|
||||
T = 293.6 # K
|
||||
E0 = 0.012 # eV
|
||||
n_particles = 100000
|
||||
T = 293.6 # K
|
||||
E0 = 0.012 # eV
|
||||
cfg = 'Al_sg225.ncmat'
|
||||
test = pencil_beam_model(cfg, E0, NParticles)
|
||||
test = pencil_beam_model(cfg, E0, n_particles)
|
||||
harness = NCrystalTest('statepoint.10.h5', model=test)
|
||||
harness.main()
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue