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

openmc-dev / openmc / 25631610076

10 May 2026 02:44PM UTC coverage: 81.026% (-0.4%) from 81.388%
25631610076

Pull #3757

github

web-flow
Merge debc5a921 into d56cda254
Pull Request #3757: Testing point detectors

17769 of 25812 branches covered (68.84%)

Branch coverage included in aggregate %.

51 of 373 new or added lines in 25 files covered. (13.67%)

3 existing lines in 2 files now uncovered.

58805 of 68694 relevant lines covered (85.6%)

40257430.53 hits per line

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

62.14
/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 "openmc/tensor.h"
7
#include <fmt/core.h>
8

9
#include "openmc/constants.h"
10
#include "openmc/endf.h"
11
#include "openmc/error.h"
12
#include "openmc/random_lcg.h"
13
#include "openmc/search.h"
14
#include "openmc/secondary_correlated.h"
15
#include "openmc/secondary_thermal.h"
16
#include "openmc/settings.h"
17
#include "openmc/string_utils.h"
18

19
namespace openmc {
20

21
//==============================================================================
22
// Global variables
23
//==============================================================================
24

25
namespace data {
26
std::unordered_map<std::string, int> thermal_scatt_map;
27
vector<unique_ptr<ThermalScattering>> thermal_scatt;
28
} // namespace data
29

30
//==============================================================================
31
// ThermalScattering implementation
32
//==============================================================================
33

34
ThermalScattering::ThermalScattering(
963✔
35
  hid_t group, const vector<double>& temperature)
963✔
36
{
37
  // Get name of table from group
38
  name_ = object_name(group);
963✔
39

40
  // Get rid of leading '/'
41
  name_ = name_.substr(1);
963✔
42

43
  read_attribute(group, "atomic_weight_ratio", awr_);
963✔
44
  read_attribute(group, "energy_max", energy_max_);
963✔
45
  read_attribute(group, "nuclides", nuclides_);
963✔
46

47
  // Read temperatures
48
  hid_t kT_group = open_group(group, "kTs");
963✔
49

50
  // Determine temperatures available
51
  auto dset_names = dataset_names(kT_group);
963✔
52
  auto n = dset_names.size();
963✔
53
  auto temps_available = tensor::Tensor<double>({n});
963✔
54
  for (int i = 0; i < dset_names.size(); ++i) {
1,926✔
55
    // Read temperature value
56
    double T;
963✔
57
    read_dataset(kT_group, dset_names[i].data(), T);
963✔
58
    temps_available[i] = std::round(T / K_BOLTZMANN);
963✔
59
  }
60
  std::sort(temps_available.begin(), temps_available.end());
963✔
61

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

75
  switch (settings::temperature_method) {
963!
76
  case TemperatureMethod::NEAREST:
947✔
77
    // Determine actual temperatures to read
78
    for (const auto& T : temperature) {
1,894✔
79

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

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

143
  // Sort temperatures to read
144
  std::sort(temps_to_read.begin(), temps_to_read.end());
963✔
145

146
  auto n_temperature = temps_to_read.size();
963✔
147
  kTs_.reserve(n_temperature);
963✔
148
  data_.reserve(n_temperature);
963✔
149

150
  for (auto T : temps_to_read) {
1,926✔
151
    // Get temperature as a string
152
    std::string temp_str = fmt::format("{}K", T);
963✔
153

154
    // Read exact temperature value
155
    double kT;
963✔
156
    read_dataset(kT_group, temp_str.data(), kT);
963✔
157
    kTs_.push_back(kT);
963✔
158

159
    // Open group for this temperature
160
    hid_t T_group = open_group(group, temp_str.data());
963✔
161
    data_.emplace_back(T_group);
963✔
162
    close_group(T_group);
963✔
163
  }
963✔
164

165
  close_group(kT_group);
963✔
166
}
1,926✔
167

168
void ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
108,121,558✔
169
  double* elastic, double* inelastic, uint64_t* seed) const
170
{
171
  // Determine temperature for S(a,b) table
172
  double kT = sqrtkT * sqrtkT;
108,121,558✔
173
  int i = 0;
108,121,558✔
174

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

201
  // Set temperature index
202
  *i_temp = i;
108,121,558✔
203

204
  // Calculate cross sections for ith temperature
205
  data_[i].calculate_xs(E, elastic, inelastic);
108,121,558✔
206
}
108,121,558✔
207

208
bool ThermalScattering::has_nuclide(const char* name) const
×
209
{
210
  std::string nuc {name};
×
211
  return std::find(nuclides_.begin(), nuclides_.end(), nuc) != nuclides_.end();
×
212
}
×
213

214
//==============================================================================
215
// ThermalData implementation
216
//==============================================================================
217

218
ThermalData::ThermalData(hid_t group)
963✔
219
{
220
  // Coherent/incoherent elastic data
221
  if (object_exists(group, "elastic")) {
963✔
222
    // Read cross section data
223
    hid_t elastic_group = open_group(group, "elastic");
98✔
224

225
    // Read elastic cross section
226
    elastic_.xs = read_function(elastic_group, "xs");
98✔
227

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

248
      // Create mixed elastic distribution
249
      elastic_.distribution =
32✔
250
        make_unique<MixedElasticAE>(dgroup, *coh_xs, *incoh_xs);
32✔
251
    }
252

253
    close_group(elastic_group);
98✔
254
  }
98✔
255

256
  // Inelastic data
257
  if (object_exists(group, "inelastic")) {
963!
258
    // Read type of inelastic data
259
    hid_t inelastic_group = open_group(group, "inelastic");
963✔
260

261
    // Read inelastic cross section
262
    inelastic_.xs = read_function(inelastic_group, "xs");
963✔
263

264
    // Read angle-energy distribution
265
    hid_t dgroup = open_group(inelastic_group, "distribution");
963✔
266
    std::string temp;
963✔
267
    read_attribute(dgroup, "type", temp);
963✔
268
    if (temp == "incoherent_inelastic") {
963✔
269
      inelastic_.distribution = make_unique<IncoherentInelasticAE>(dgroup);
32✔
270
    } else if (temp == "incoherent_inelastic_discrete") {
931!
271
      auto xs = dynamic_cast<Tabulated1D*>(inelastic_.xs.get());
931!
272
      inelastic_.distribution =
931✔
273
        make_unique<IncoherentInelasticAEDiscrete>(dgroup, xs->x());
931✔
274
    }
275

276
    close_group(inelastic_group);
963✔
277
  }
963✔
278
}
963✔
279

280
void ThermalData::calculate_xs(
108,121,558✔
281
  double E, double* elastic, double* inelastic) const
282
{
283
  // Calculate thermal elastic scattering cross section
284
  if (elastic_.xs) {
108,121,558✔
285
    *elastic = (*elastic_.xs)(E);
3,558,888✔
286
  } else {
287
    *elastic = 0.0;
104,562,670✔
288
  }
289

290
  // Calculate thermal inelastic scattering cross section
291
  *inelastic = (*inelastic_.xs)(E);
108,121,558✔
292
}
108,121,558✔
293

294
AngleEnergy& ThermalData::sample_dist(
92,714,744✔
295
  const NuclideMicroXS& micro_xs, double E, uint64_t* seed) const
296
{
297
  // Determine whether inelastic or elastic scattering will occur
298
  if (prn(seed) < micro_xs.thermal_elastic / micro_xs.thermal) {
92,714,744✔
299
    return *elastic_.distribution;
130,440✔
300
  } else {
301
    return *inelastic_.distribution;
92,584,304✔
302
  }
303
}
304

305
void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
92,714,744✔
306
  double* E_out, double* mu, uint64_t* seed) const
307
{
308
  sample_dist(micro_xs, E, seed).sample(E, *E_out, *mu, seed);
92,714,744✔
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)
92,714,744!
313
    *mu = std::copysign(1.0, *mu);
×
314
}
92,714,744✔
315

316
double ThermalData::sample_energy_and_pdf(const NuclideMicroXS& micro_xs,
×
317
  double E_in, double mu, double& E_out, uint64_t* seed, bool is_com,
318
  double awr) const
319
{
320
  return sample_dist(micro_xs, E_in, seed)
×
NEW
321
    .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
×
322
}
323

324
void free_memory_thermal()
6,204✔
325
{
326
  data::thermal_scatt.clear();
6,204✔
327
  data::thermal_scatt_map.clear();
6,204✔
328
}
6,204✔
329

330
} // 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

© 2026 Coveralls, Inc