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

Implement recommended changes

Browse source
This commit is contained in:
Jose Ignacio Marquez Damian 2022-09-26 11:19:54 +02:00 • committed by Thomas Kittelmann
parent 29a10eafd3
commit 2e19260dea
6 changed files with 35 additions and 42 deletions
Download patch file Download diff file Expand all files Collapse all files

2
CMakeLists.txt
View file

@ -104,7 +104,7 @@ endmacro()
#===============================================================================
if(OPENMC_USE_NCRYSTAL)
find_package(NCrystal REQUIRED PATH_SUFFIXES ./)#fixme
find_package(NCrystal REQUIRED)
endif()
#===============================================================================

8
include/openmc/material.h
View file

@ -156,9 +156,9 @@ public:
double temperature() const;
#ifdef NCRYSTAL
//! Gwet pointer to NCrystal material object
//! Get pointer to NCrystal material object
//! \return Pointer to NCrystal material object
std::shared_ptr<const NCrystal::ProcImpl::Process> m_NCrystal_mat() const {return m_NCrystal_mat_; };
std::shared_ptr<const NCrystal::ProcImpl::Process> ncrystal_mat() const {return ncrystal_mat_; };
#endif
//----------------------------------------------------------------------------
@ -168,8 +168,8 @@ public:
vector<int> nuclide_; //!< Indices in nuclides vector
vector<int> element_; //!< Indices in elements vector
#ifdef NCRYSTAL
std::string cfg_; //!< NCrystal configuration string
std::shared_ptr<const NCrystal::ProcImpl::Process> m_NCrystal_mat_;
std::string ncrystal_cfg_; //!< NCrystal configuration string
std::shared_ptr<const NCrystal::ProcImpl::Process> ncrystal_mat_;
#endif
xt::xtensor<double, 1> atom_density_; //!< Nuclide atom density in [atom/b-cm]
double density_; //!< Total atom density in [atom/b-cm]

10
include/openmc/physics.h
View file

@ -110,17 +110,11 @@ void split_particle(Particle& p);
// NCrystal wrapper class for the OpenMC random number generator
//==============================================================================
class NcrystalRNG_Wrapper : public NCrystal::RNGStream {
class NCrystalRNGWrapper : public NCrystal::RNGStream {
uint64_t* openmc_seed_;
public:
constexpr NcrystalRNG_Wrapper(uint64_t* s) noexcept : openmc_seed_(s) {}
//Can be cheaply created on the stack just before being used in calls to
//ProcImpl::Scatter objects, like:
//
// RNG_Wrapper rng(seed);
//
constexpr NCrystalRNGWrapper(uint64_t* seed) noexcept : openmc_seed_(seed) {}
protected:
//double actualGenerate() override { return prn(openmc_seed_); }
double actualGenerate() override {
return std::max<double>( std::numeric_limits<double>::min(), prn(openmc_seed_) );
}

31
openmc/material.py
View file

@ -118,7 +118,7 @@ class Material(IDManagerMixin):
self._volume = None
self._atoms = {}
self._isotropic = []
self._NCrystal_cfg = None
self._ncrystal_cfg = None
# A list of tuples (nuclide, percent, percent type)
self._nuclides = []
@ -141,7 +141,7 @@ class Material(IDManagerMixin):
string += '{: <16}\n'.format('\tS(a,b) Tables')
string += '{: <16}=\t{}\n'.format('\tNCrystal conf', 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)
@ -335,7 +335,7 @@ class Material(IDManagerMixin):
return material
@classmethod
def from_NCrystal(cls, cfg):
def from_ncrystal(cls, cfg):
"""Create material from NCrystal configuration string
Density is set from the NCrystal value,
and material temperature from the configuration string.
@ -352,26 +352,23 @@ class Material(IDManagerMixin):
"""
try:
import NCrystal
except ImportError:
raise SystemExit("ERROR: NCrystal Python module is required.")
import NCrystal
NC_mat = NCrystal.createInfo(cfg)
NC_comp = NC_mat.getComposition()
nc_mat = NCrystal.createInfo(cfg)
nc_comp = nc_mat.getComposition()
# Create the Material
material = cls()
for frac, atom in NC_comp:
for frac, atom in nc_comp:
if not atom.isNaturalElement():
raise ValueError('NCrystal-OpenMC interface only works with natural elements for now.')
material.add_element(atom.elementName(), frac, 'ao')
material._NCrystal_cfg = cfg
material._ncrystal_cfg = cfg
material._density_units = "g/cm3"
material._density = NC_mat.getDensity()
material.temperature = NC_mat.getTemperature()
material._density = nc_mat.getDensity()
material.temperature = nc_mat.getTemperature()
return material
@ -449,7 +446,7 @@ class Material(IDManagerMixin):
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if self._NCrystal_cfg is not None:
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
@ -656,7 +653,7 @@ 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:
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
@ -1185,13 +1182,13 @@ class Material(IDManagerMixin):
if self._volume:
element.set("volume", str(self._volume))
if self._NCrystal_cfg:
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 macroscopic cross sections.")
element.set("cfg", str(self._NCrystal_cfg))
element.set("cfg", str(self._ncrystal_cfg))
# Create temperature XML subelement
if self.temperature is not None:

16
src/material.cpp
View file

@ -63,8 +63,8 @@ Material::Material(pugi::xml_node node)
#ifdef NCRYSTAL
if (check_for_node(node, "cfg")) {
cfg_ = get_node_value(node, "cfg");
write_message(5, "NCrystal config string: >>{}<< ", cfg_);
m_NCrystal_mat_ = NCrystal::FactImpl::createScatter(cfg_);
write_message(5, "NCrystal config string: '{}'", ncrystal_cfg_);
ncrystal_mat_ = NCrystal::FactImpl::createScatter(ncrystal_cfg_);
}
#endif
@ -803,11 +803,11 @@ void Material::calculate_neutron_xs(Particle& p) const
#ifdef NCRYSTAL
double ncrystal_xs = -1;
if (m_NCrystal_mat_ != nullptr && p.E() < settings::ncrystal_max_energy){
if (ncrystal_mat_ && p.E() < settings::ncrystal_max_energy){
// Calculate scattering XS per atom with NCrystal, only once per material
NCrystal::CachePtr dummyCache;
NCrystal::CachePtr dummy_cache;
auto nc_energy = NCrystal::NeutronEnergy{p.E()};
ncrystal_xs = m_NCrystal_mat_->crossSection( dummyCache, nc_energy, {p.u().x, p.u().y, p.u().z}).get();
ncrystal_xs = ncrystal_mat_->crossSection( dummy_cache, nc_energy, {p.u().x, p.u().y, p.u().z}).get();
}
#endif
@ -849,8 +849,10 @@ void Material::calculate_neutron_xs(Particle& p) const
int i_nuclide = nuclide_[i];
// Calculate microscopic cross section for this nuclide
#ifndef NCRYSTAL
const
#ifdef NCRYSTAL
auto& micro {p.neutron_xs(i_nuclide)};
#else
const auto& micro {p.neutron_xs(i_nuclide)};
#endif
auto& micro {p.neutron_xs(i_nuclide)};
if (p.E() != micro.last_E || p.sqrtkT() != micro.last_sqrtkT ||

10
src/physics.cpp
View file

@ -142,18 +142,18 @@ void sample_neutron_reaction(Particle& p)
// exiting neutron
#ifdef NCRYSTAL
if (model::materials[p.material()]->m_NCrystal_mat() != nullptr && p.E() < settings::ncrystal_max_energy){
NcrystalRNG_Wrapper rng(p.current_seed()); // Initialize RNG
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 dummyCache;//fixme: avoid recreating here (triggers malloc)
NCrystal::CachePtr dummy_cache;
auto nc_energy = NCrystal::NeutronEnergy{p.E()};
auto outcome = model::materials[p.material()]->m_NCrystal_mat()->sampleScatter( dummyCache, rng , nc_energy, {p.u().x, p.u().y, p.u().z});
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;//fixme: define a new label for NCrystal?
p.event_mt() = ELASTIC;
} else {
scatter(p, i_nuclide);
}