• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 18538152141

15 Oct 2025 06:02PM UTC coverage: 81.97% (-3.2%) from 85.194%
18538152141

Pull #3417

github

web-flow
Merge 4604e1321 into e9077b137
Pull Request #3417: Addition of a collision tracking feature

16794 of 23357 branches covered (71.9%)

Branch coverage included in aggregate %.

480 of 522 new or added lines in 13 files covered. (91.95%)

457 existing lines in 53 files now uncovered.

54128 of 63165 relevant lines covered (85.69%)

42776927.64 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

55.49
/src/thermal.cpp
1
#include "openmc/thermal.h"
2

3
#include <algorithm> // for sort, move, min, max, find
4
#include <cmath>     // for round, sqrt, abs
5

6
#include "xtensor/xarray.hpp"
7
#include "xtensor/xbuilder.hpp"
8
#include "xtensor/xmath.hpp"
9
#include "xtensor/xsort.hpp"
10
#include "xtensor/xtensor.hpp"
11
#include "xtensor/xview.hpp"
12
#include <fmt/core.h>
13

14
#include "openmc/constants.h"
15
#include "openmc/endf.h"
16
#include "openmc/error.h"
17
#include "openmc/random_lcg.h"
18
#include "openmc/search.h"
19
#include "openmc/secondary_correlated.h"
20
#include "openmc/secondary_thermal.h"
21
#include "openmc/settings.h"
22
#include "openmc/string_utils.h"
23

24
namespace openmc {
25

26
//==============================================================================
27
// Global variables
28
//==============================================================================
29

30
namespace data {
31
std::unordered_map<std::string, int> thermal_scatt_map;
32
vector<unique_ptr<ThermalScattering>> thermal_scatt;
33
} // namespace data
34

35
//==============================================================================
36
// ThermalScattering implementation
37
//==============================================================================
38

39
ThermalScattering::ThermalScattering(
1,161✔
40
  hid_t group, const vector<double>& temperature)
1,161✔
41
{
42
  // Get name of table from group
43
  name_ = object_name(group);
1,161✔
44

45
  // Get rid of leading '/'
46
  name_ = name_.substr(1);
1,161✔
47

48
  read_attribute(group, "atomic_weight_ratio", awr_);
1,161✔
49
  read_attribute(group, "energy_max", energy_max_);
1,161✔
50
  read_attribute(group, "nuclides", nuclides_);
1,161✔
51

52
  // Read temperatures
53
  hid_t kT_group = open_group(group, "kTs");
1,161✔
54

55
  // Determine temperatures available
56
  auto dset_names = dataset_names(kT_group);
1,161✔
57
  auto n = dset_names.size();
1,161✔
58
  auto temps_available = xt::empty<double>({n});
1,161✔
59
  for (int i = 0; i < dset_names.size(); ++i) {
2,322✔
60
    // Read temperature value
61
    double T;
62
    read_dataset(kT_group, dset_names[i].data(), T);
1,161✔
63
    temps_available[i] = std::round(T / K_BOLTZMANN);
1,161✔
64
  }
65
  std::sort(temps_available.begin(), temps_available.end());
1,161✔
66

67
  // Determine actual temperatures to read -- start by checking whether a
68
  // temperature range was given, in which case all temperatures in the range
69
  // are loaded irrespective of what temperatures actually appear in the model
70
  vector<int> temps_to_read;
1,161✔
71
  if (settings::temperature_range[1] > 0.0) {
1,161!
72
    for (const auto& T : temps_available) {
×
73
      if (settings::temperature_range[0] <= T &&
×
74
          T <= settings::temperature_range[1]) {
×
75
        temps_to_read.push_back(std::round(T));
×
76
      }
77
    }
78
  }
79

80
  switch (settings::temperature_method) {
1,161!
81
  case TemperatureMethod::NEAREST:
1,139✔
82
    // Determine actual temperatures to read
83
    for (const auto& T : temperature) {
2,278✔
84

85
      auto i_closest = xt::argmin(xt::abs(temps_available - T))[0];
1,139✔
86
      auto temp_actual = temps_available[i_closest];
1,139✔
87
      if (std::abs(temp_actual - T) < settings::temperature_tolerance) {
1,139!
88
        if (std::find(temps_to_read.begin(), temps_to_read.end(),
1,139✔
89
              std::round(temp_actual)) == temps_to_read.end()) {
2,278!
90
          temps_to_read.push_back(std::round(temp_actual));
1,139✔
91
        }
92
      } else {
93
        fatal_error(fmt::format(
×
94
          "Nuclear data library does not contain cross sections "
95
          "for {}  at or near {} K. Available temperatures "
96
          "are {} K. Consider making use of openmc.Settings.temperature "
97
          "to specify how intermediate temperatures are treated.",
98
          name_, std::round(T), concatenate(temps_available)));
×
99
      }
100
    }
101
    break;
1,139✔
102

103
  case TemperatureMethod::INTERPOLATION:
22✔
104
    // If temperature interpolation or multipole is selected, get a list of
105
    // bounding temperatures for each actual temperature present in the model
106
    for (const auto& T : temperature) {
44✔
107
      bool found = false;
22✔
108
      for (int j = 0; j < temps_available.size() - 1; ++j) {
22!
109
        if (temps_available[j] <= T && T < temps_available[j + 1]) {
×
110
          int T_j = temps_available[j];
×
111
          int T_j1 = temps_available[j + 1];
×
112
          if (std::find(temps_to_read.begin(), temps_to_read.end(), T_j) ==
×
113
              temps_to_read.end()) {
×
114
            temps_to_read.push_back(T_j);
×
115
          }
116
          if (std::find(temps_to_read.begin(), temps_to_read.end(), T_j1) ==
×
117
              temps_to_read.end()) {
×
118
            temps_to_read.push_back(T_j1);
×
119
          }
120
          found = true;
×
121
        }
122
      }
123
      if (!found) {
22!
124
        // If no pairs found, check if the desired temperature falls within
125
        // bounds' tolerance
126
        if (std::abs(T - temps_available[0]) <=
22!
127
            settings::temperature_tolerance) {
128
          if (std::find(temps_to_read.begin(), temps_to_read.end(),
22✔
129
                temps_available[0]) == temps_to_read.end()) {
44!
130
            temps_to_read.push_back(temps_available[0]);
22✔
131
          }
132
        } else if (std::abs(T - temps_available[n - 1]) <=
×
133
                   settings::temperature_tolerance) {
134
          if (std::find(temps_to_read.begin(), temps_to_read.end(),
×
135
                temps_available[n - 1]) == temps_to_read.end()) {
×
136
            temps_to_read.push_back(temps_available[n - 1]);
×
137
          }
138
        } else {
139
          fatal_error(
×
140
            fmt::format("Nuclear data library does not contain cross "
×
141
                        "sections for {} at temperatures that bound {} K.",
142
              name_, std::round(T)));
×
143
        }
144
      }
145
    }
146
  }
147

148
  // Sort temperatures to read
149
  std::sort(temps_to_read.begin(), temps_to_read.end());
1,161✔
150

151
  auto n_temperature = temps_to_read.size();
1,161✔
152
  kTs_.reserve(n_temperature);
1,161✔
153
  data_.reserve(n_temperature);
1,161✔
154

155
  for (auto T : temps_to_read) {
2,322✔
156
    // Get temperature as a string
157
    std::string temp_str = fmt::format("{}K", T);
963✔
158

159
    // Read exact temperature value
160
    double kT;
161
    read_dataset(kT_group, temp_str.data(), kT);
1,161✔
162
    kTs_.push_back(kT);
1,161✔
163

164
    // Open group for this temperature
165
    hid_t T_group = open_group(group, temp_str.data());
1,161✔
166
    data_.emplace_back(T_group);
1,161✔
167
    close_group(T_group);
1,161✔
168
  }
1,161✔
169

170
  close_group(kT_group);
1,161✔
171
}
1,161✔
172

173
void ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
87,481,507✔
174
  double* elastic, double* inelastic, uint64_t* seed) const
175
{
176
  // Determine temperature for S(a,b) table
177
  double kT = sqrtkT * sqrtkT;
87,481,507✔
178
  int i = 0;
87,481,507✔
179

180
  auto n = kTs_.size();
87,481,507✔
181
  if (n > 1) {
87,481,507!
182
    if (settings::temperature_method == TemperatureMethod::NEAREST) {
×
183
      while (kTs_[i + 1] < kT && i + 1 < n - 1)
×
184
        ++i;
×
185
      // Pick closer of two bounding temperatures
186
      if (kT - kTs_[i] > kTs_[i + 1] - kT)
×
187
        ++i;
×
188
    } else {
189
      // If current kT outside of the bounds of available, snap to the bound
190
      if (kT < kTs_.front()) {
×
191
        i = 0;
×
192
      } else if (kT > kTs_.back()) {
×
193
        i = kTs_.size() - 1;
×
194
      } else {
195
        // Find temperatures that bound the actual temperature
196
        while (kTs_[i + 1] < kT && i + 1 < n - 1)
×
197
          ++i;
×
198
        // Randomly sample between temperature i and i+1
199
        double f = (kT - kTs_[i]) / (kTs_[i + 1] - kTs_[i]);
×
200
        if (f > prn(seed))
×
201
          ++i;
×
202
      }
203
    }
204
  }
205

206
  // Set temperature index
207
  *i_temp = i;
87,481,507✔
208

209
  // Calculate cross sections for ith temperature
210
  data_[i].calculate_xs(E, elastic, inelastic);
87,481,507✔
211
}
87,481,507✔
212

213
bool ThermalScattering::has_nuclide(const char* name) const
×
214
{
215
  std::string nuc {name};
×
216
  return std::find(nuclides_.begin(), nuclides_.end(), nuc) != nuclides_.end();
×
UNCOV
217
}
×
218

219
//==============================================================================
220
// ThermalData implementation
221
//==============================================================================
222

223
ThermalData::ThermalData(hid_t group)
1,161✔
224
{
225
  // Coherent/incoherent elastic data
226
  if (object_exists(group, "elastic")) {
1,161✔
227
    // Read cross section data
228
    hid_t elastic_group = open_group(group, "elastic");
140✔
229

230
    // Read elastic cross section
231
    elastic_.xs = read_function(elastic_group, "xs");
140✔
232

233
    // Read angle-energy distribution
234
    hid_t dgroup = open_group(elastic_group, "distribution");
140✔
235
    std::string temp;
140✔
236
    read_attribute(dgroup, "type", temp);
140✔
237
    if (temp == "coherent_elastic") {
140✔
238
      auto xs = dynamic_cast<CoherentElasticXS*>(elastic_.xs.get());
96!
239
      elastic_.distribution = make_unique<CoherentElasticAE>(*xs);
96✔
240
    } else if (temp == "incoherent_elastic") {
44!
241
      elastic_.distribution = make_unique<IncoherentElasticAE>(dgroup);
×
242
    } else if (temp == "incoherent_elastic_discrete") {
44!
243
      auto xs = dynamic_cast<Tabulated1D*>(elastic_.xs.get());
×
244
      elastic_.distribution =
245
        make_unique<IncoherentElasticAEDiscrete>(dgroup, xs->x());
×
246
    } else if (temp == "mixed_elastic") {
44!
247
      // Get coherent/incoherent cross sections
248
      auto mixed_xs = dynamic_cast<Sum1D*>(elastic_.xs.get());
44!
249
      const auto& coh_xs =
250
        dynamic_cast<const CoherentElasticXS*>(mixed_xs->functions(0).get());
44!
251
      const auto& incoh_xs = mixed_xs->functions(1).get();
44✔
252

253
      // Create mixed elastic distribution
254
      elastic_.distribution =
255
        make_unique<MixedElasticAE>(dgroup, *coh_xs, *incoh_xs);
44✔
256
    }
257

258
    close_group(elastic_group);
140✔
259
  }
140✔
260

261
  // Inelastic data
262
  if (object_exists(group, "inelastic")) {
1,161!
263
    // Read type of inelastic data
264
    hid_t inelastic_group = open_group(group, "inelastic");
1,161✔
265

266
    // Read inelastic cross section
267
    inelastic_.xs = read_function(inelastic_group, "xs");
1,161✔
268

269
    // Read angle-energy distribution
270
    hid_t dgroup = open_group(inelastic_group, "distribution");
1,161✔
271
    std::string temp;
1,161✔
272
    read_attribute(dgroup, "type", temp);
1,161✔
273
    if (temp == "incoherent_inelastic") {
1,161✔
274
      inelastic_.distribution = make_unique<IncoherentInelasticAE>(dgroup);
44✔
275
    } else if (temp == "incoherent_inelastic_discrete") {
1,117!
276
      auto xs = dynamic_cast<Tabulated1D*>(inelastic_.xs.get());
1,117!
277
      inelastic_.distribution =
278
        make_unique<IncoherentInelasticAEDiscrete>(dgroup, xs->x());
1,117✔
279
    }
280

281
    close_group(inelastic_group);
1,161✔
282
  }
1,161✔
283
}
1,161✔
284

285
void ThermalData::calculate_xs(
87,481,507✔
286
  double E, double* elastic, double* inelastic) const
287
{
288
  // Calculate thermal elastic scattering cross section
289
  if (elastic_.xs) {
87,481,507✔
290
    *elastic = (*elastic_.xs)(E);
4,893,471✔
291
  } else {
292
    *elastic = 0.0;
82,588,036✔
293
  }
294

295
  // Calculate thermal inelastic scattering cross section
296
  *inelastic = (*inelastic_.xs)(E);
87,481,507✔
297
}
87,481,507✔
298

299
void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
72,177,629✔
300
  double* E_out, double* mu, uint64_t* seed)
301
{
302
  // Determine whether inelastic or elastic scattering will occur
303
  if (prn(seed) < micro_xs.thermal_elastic / micro_xs.thermal) {
72,177,629✔
304
    elastic_.distribution->sample(E, *E_out, *mu, seed);
179,355✔
305
  } else {
306
    inelastic_.distribution->sample(E, *E_out, *mu, seed);
71,998,274✔
307
  }
308

309
  // Because of floating-point roundoff, it may be possible for mu to be
310
  // outside of the range [-1,1). In these cases, we just set mu to exactly
311
  // -1 or 1
312
  if (std::abs(*mu) > 1.0)
72,177,629!
313
    *mu = std::copysign(1.0, *mu);
×
314
}
72,177,629✔
315

316
void free_memory_thermal()
7,877✔
317
{
318
  data::thermal_scatt.clear();
7,877✔
319
  data::thermal_scatt_map.clear();
7,877✔
320
}
7,877✔
321

322
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc