Repository navigation
Expand file tree
/
Copy pathncrystal_interface.cpp
More file actions
108 lines (90 loc) · 2.96 KB
/
Copy pathncrystal_interface.cpp
File metadata and controls
108 lines (90 loc) · 2.96 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
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