-
Notifications
You must be signed in to change notification settings - Fork 31
/
NeutronInelastic.test.cc
179 lines (157 loc) · 7.14 KB
/
NeutronInelastic.test.cc
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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
//----------------------------------*-C++-*----------------------------------//
// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: (Apache-2.0 OR MIT)
//---------------------------------------------------------------------------//
//! \file celeritas/neutron/NeutronInelastic.test.cc
//---------------------------------------------------------------------------//
#include <memory>
#include "corecel/cont/Range.hh"
#include "corecel/data/Ref.hh"
#include "corecel/math/ArrayUtils.hh"
#include "celeritas/Quantities.hh"
#include "celeritas/grid/GenericGridData.hh"
#include "celeritas/io/NeutronXsReader.hh"
#include "celeritas/mat/MaterialTrackView.hh"
#include "celeritas/neutron/NeutronTestBase.hh"
#include "celeritas/neutron/interactor/NeutronInelasticInteractor.hh"
#include "celeritas/neutron/model/NeutronInelasticModel.hh"
#include "celeritas/neutron/xs/NeutronInelasticMicroXsCalculator.hh"
#include "celeritas/neutron/xs/NucleonNucleonXsCalculator.hh"
#include "celeritas/phys/MacroXsCalculator.hh"
#include "celeritas_test.hh"
namespace celeritas
{
namespace test
{
//---------------------------------------------------------------------------//
class NeutronInelasticTest : public NeutronTestBase
{
protected:
using MevEnergy = units::MevEnergy;
using SPConstNInelasticModel = std::shared_ptr<NeutronInelasticModel const>;
void SetUp() override
{
using namespace units;
// Load neutron elastic cross section data
std::string data_path = this->test_data_path("celeritas", "");
NeutronXsReader read_el_data(NeutronXsType::inel, data_path.c_str());
// Set up the default particle: 100 MeV neutron along +z direction
auto const& particles = *this->particle_params();
this->set_inc_particle(pdg::neutron(), MevEnergy{100});
this->set_inc_direction({0, 0, 1});
// Set up the default material
this->set_material("HeCu");
model_ = std::make_shared<NeutronInelasticModel>(
ActionId{0}, particles, *this->material_params(), read_el_data);
}
protected:
SPConstNInelasticModel model_;
};
//---------------------------------------------------------------------------//
// TESTS
//---------------------------------------------------------------------------//
TEST_F(NeutronInelasticTest, micro_xs)
{
// Calculate the elastic neutron-nucleus microscopic cross section
using XsCalculator = NeutronInelasticMicroXsCalculator;
// Set the target element: Cu
ElementId el_id{1};
// Check the size of the element cross section data (G4PARTICLEXS4.0)
// The neutron/inelZ data are pruned by a factor of 5 for this test
NeutronInelasticRef shared = model_->host_ref();
GenericGridData grid = shared.micro_xs[el_id];
EXPECT_EQ(grid.grid.size(), 61);
// Microscopic cross section (units::BarnXs) in [1e-04:1e+4] (MeV)
std::vector<real_type> const expected_micro_xs = {2.499560005640001e-06,
2.499560005640001e-06,
2.499560005640001e-06,
2.499560005640001e-06,
0.2170446680979802,
1.3677671823188946,
0.81016638725225387,
0.84789596907525477};
real_type energy = 1e-4;
real_type const factor = 1e+1;
for (auto i : range(expected_micro_xs.size()))
{
XsCalculator calc_micro_xs(shared, MevEnergy{energy});
EXPECT_SOFT_EQ(calc_micro_xs(el_id).value(), expected_micro_xs[i]);
energy *= factor;
}
// Check the elastic cross section at the upper bound (20 GeV)
XsCalculator calc_upper_xs(shared, MevEnergy{2e+4});
EXPECT_SOFT_EQ(calc_upper_xs(el_id).value(), 0.80300000000000027);
}
TEST_F(NeutronInelasticTest, macro_xs)
{
// Calculate the inelastic neutron-nucleus macroscopic cross section
auto material = this->material_track().make_material_view();
auto calc_xs = MacroXsCalculator<NeutronInelasticMicroXsCalculator>(
model_->host_ref(), material);
// Macroscopic cross section (\f$ cm^{-1} \f$) in [1e-04:1e+4] (MeV)
std::vector<real_type> const expected_macro_xs = {1.0577605656430734e-06,
4.4447010621996484e-07,
2.5134945234021254e-07,
1.9270371228950039e-07,
0.015057496086707027,
0.094888935102106969,
0.056850427191922973,
0.059657345679963072};
real_type energy = 1e-4;
real_type const factor = 1e+1;
for (auto i : range(expected_macro_xs.size()))
{
EXPECT_SOFT_EQ(
native_value_to<units::InvCmXs>(calc_xs(MevEnergy{energy})).value(),
expected_macro_xs[i]);
energy *= factor;
}
// Check the neutron inelastic interaction cross section at the upper bound
// (20 GeV)
EXPECT_SOFT_EQ(
native_value_to<units::InvCmXs>(calc_xs(MevEnergy{2000})).value(),
0.061219850473480573);
}
TEST_F(NeutronInelasticTest, nucleon_xs)
{
// Calculate nucleon-nucleon cross sections (units::BarnXs)
NeutronInelasticRef shared = model_->host_ref();
NucleonNucleonXsCalculator calc_xs(shared);
size_type num_channels = shared.scalars.num_channels();
EXPECT_EQ(num_channels, 3);
std::vector<real_type> xs_zero;
std::vector<real_type> xs;
for (auto ch_id : range(ChannelId{num_channels}))
{
xs_zero.push_back(shared.xs_params[ch_id].xs_zero);
for (real_type inc_e : {0.01, 0.1, 1., 10., 100., 320.})
{
xs.push_back(calc_xs(ch_id, MevEnergy{inc_e}).value());
}
}
real_type const expected_xs_zero[] = {17.613, 20.36, 17.613};
real_type const expected_xs[] = {17.613,
17.613,
4.0,
0.8633,
0.0691,
0.0351,
20.36,
19.2,
1.92,
0.3024,
0.0316,
0.0233,
17.613,
17.613,
4.0,
0.8633,
0.0691,
0.0351};
EXPECT_VEC_SOFT_EQ(expected_xs_zero, xs_zero);
EXPECT_VEC_SOFT_EQ(expected_xs, xs);
}
//---------------------------------------------------------------------------//
} // namespace test
} // namespace celeritas