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

openmc-dev / openmc / 23460042813

23 Mar 2026 09:04PM UTC coverage: 81.334% (-0.2%) from 81.556%
23460042813

Pull #3832

github

web-flow
Merge 1144e8851 into 6cd39073b
Pull Request #3832: Allowing multiple cells in surface source banking

17652 of 25502 branches covered (69.22%)

Branch coverage included in aggregate %.

94 of 97 new or added lines in 4 files covered. (96.91%)

868 existing lines in 32 files now uncovered.

58133 of 67676 relevant lines covered (85.9%)

44449215.59 hits per line

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

75.66
/src/secondary_thermal.cpp
1
#include "openmc/secondary_thermal.h"
2

3
#include "openmc/hdf5_interface.h"
4
#include "openmc/math_functions.h"
5
#include "openmc/random_lcg.h"
6
#include "openmc/search.h"
7
#include "openmc/vector.h"
8

9
#include "openmc/tensor.h"
10

11
#include <cassert>
12
#include <cmath> // for log, exp
13

14
namespace openmc {
15

16
//==============================================================================
17
// CoherentElasticAE implementation
18
//==============================================================================
19

20
CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs) : xs_ {xs}
134✔
21
{
22
  const auto& bragg = xs_.bragg_edges();
134✔
23
  auto n = bragg.size();
134✔
24
  bragg_edges_ = tensor::Tensor<double>(bragg.data(), n);
134✔
25

26
  const auto& factors = xs_.factors();
134✔
27
  factors_diff_ = tensor::zeros<double>({n});
134✔
28
  factors_diff_.slice(0) = factors[0];
134✔
29
  for (int i = 1; i < n; ++i) {
19,878✔
30
    factors_diff_.slice(i) = factors[i] - factors[i - 1];
39,488✔
31
  }
32
}
134✔
33

34
void CoherentElasticAE::sample(
178,013✔
35
  double E_in, double& E_out, double& mu, uint64_t* seed) const
36
{
37
  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
38
  E_out = E_in;
178,013✔
39
  const auto& energies {xs_.bragg_edges()};
178,013!
40
  assert(E_in >= energies.front());
178,013!
41

42
  const int i = lower_bound_index(energies.begin(), energies.end(), E_in);
178,013✔
43

44
  // Sample a Bragg edge between 1 and i
45
  // E[0] < E_in < E[i+1] -> can scatter in bragg edges 0..i
46
  const auto& factors = xs_.factors();
178,013✔
47
  const double prob = prn(seed) * factors[i];
178,013✔
48

49
  const int k = std::lower_bound(factors.begin(), factors.begin() + i, prob) -
178,013✔
50
                factors.begin();
178,013✔
51

52
  // Characteristic scattering cosine for this Bragg edge (ENDF-102, Eq. 7-2)
53
  mu = 1.0 - 2.0 * energies[k] / E_in;
178,013✔
54
}
178,013✔
55

UNCOV
56
double CoherentElasticAE::sample_energy_and_pdf(
×
57
  double E_in, double mu, double& E_out, uint64_t* seed) const
58
{
59
  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
UNCOV
60
  E_out = E_in;
×
UNCOV
61
  const auto& factors = xs_.factors();
×
62

63
  if (E_in < bragg_edges_.front())
×
64
    return 0.0;
65

UNCOV
66
  const int i =
×
UNCOV
67
    lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E_in);
×
UNCOV
68
  double E = 0.5 * (1 - mu) * E_in;
×
UNCOV
69
  double C = 0.5 * E_in / factors[i];
×
70

UNCOV
71
  return C * get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)),
×
UNCOV
72
               factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in);
×
73
}
74

75
//==============================================================================
76
// IncoherentElasticAE implementation
77
//==============================================================================
78

UNCOV
79
IncoherentElasticAE::IncoherentElasticAE(hid_t group)
×
80
{
UNCOV
81
  read_dataset(group, "debye_waller", debye_waller_);
×
UNCOV
82
}
×
83

UNCOV
84
void IncoherentElasticAE::sample(
×
85
  double E_in, double& E_out, double& mu, uint64_t* seed) const
86
{
UNCOV
87
  E_out = E_in;
×
88

89
  // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
UNCOV
90
  double c = 2 * E_in * debye_waller_;
×
UNCOV
91
  mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0;
×
UNCOV
92
}
×
UNCOV
93
double IncoherentElasticAE::sample_energy_and_pdf(
×
94
  double E_in, double mu, double& E_out, uint64_t* seed) const
95
{
UNCOV
96
  E_out = E_in;
×
97

98
  // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
UNCOV
99
  double c = 2 * E_in * debye_waller_;
×
UNCOV
100
  double A = c / (1 - std::exp(-2.0 * c)); // normalization factor
×
UNCOV
101
  return A * std::exp(-c * (1 - mu));
×
102
}
103

104
//==============================================================================
105
// IncoherentElasticAEDiscrete implementation
106
//==============================================================================
107

108
IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(
44✔
109
  hid_t group, const vector<double>& energy)
44✔
110
  : energy_ {energy}
44✔
111
{
112
  read_dataset(group, "mu_out", mu_out_);
44✔
113
}
44✔
114

115
void IncoherentElasticAEDiscrete::sample(
1,342✔
116
  double E_in, double& E_out, double& mu, uint64_t* seed) const
117
{
118
  // Get index and interpolation factor for elastic grid
119
  int i;
1,342✔
120
  double f;
1,342✔
121
  get_energy_index(energy_, E_in, i, f);
1,342✔
122

123
  // Interpolate between two discrete cosines corresponding to neighboring
124
  // incoming energies.
125

126
  // Sample outgoing cosine bin
127
  int n_mu = mu_out_.shape(1);
1,342!
128
  int k = prn(seed) * n_mu;
1,342✔
129

130
  // Rather than use the sampled discrete mu directly, it is smeared over
131
  // a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
132
  // discrete mu value itself.
133

134
  // Interpolate kth mu value between distributions at energies i and i+1
135
  mu = mu_out_(i, k) + f * (mu_out_(i + 1, k) - mu_out_(i, k));
1,342✔
136

137
  // Inteprolate (k-1)th mu value between distributions at energies i and i+1.
138
  // When k==0, pick a value that will smear the cosine out to a minimum of -1.
139
  double mu_left = (k == 0) ? -1.0 - (mu + 1.0)
1,342✔
140
                            : mu_out_(i, k - 1) +
990✔
141
                                f * (mu_out_(i + 1, k - 1) - mu_out_(i, k - 1));
990✔
142

143
  // Inteprolate (k+1)th mu value between distributions at energies i and i+1.
144
  // When k is the last discrete value, pick a value that will smear the cosine
145
  // out to a maximum of 1.
146
  double mu_right =
1,342✔
147
    (k == n_mu - 1)
1,342✔
148
      ? 1.0 + (1.0 - mu)
1,342✔
149
      : mu_out_(i, k + 1) + f * (mu_out_(i + 1, k + 1) - mu_out_(i, k + 1));
990✔
150

151
  // Smear cosine
152
  mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);
1,947✔
153

154
  // Energy doesn't change in elastic scattering
155
  E_out = E_in;
1,342✔
156
}
1,342✔
157

UNCOV
158
double IncoherentElasticAEDiscrete::sample_energy_and_pdf(
×
159
  double E_in, double mu, double& E_out, uint64_t* seed) const
160
{
161
  // Get index and interpolation factor for elastic grid
UNCOV
162
  int i;
×
UNCOV
163
  double f;
×
UNCOV
164
  get_energy_index(energy_, E_in, i, f);
×
165
  // Energy doesn't change in elastic scattering
UNCOV
166
  E_out = E_in;
×
167

UNCOV
168
  return get_pdf_discrete_interpolated(
×
UNCOV
169
    mu_out_.slice(i, tensor::all), mu_out_.slice(i + 1, tensor::all), f, mu);
×
170
}
171

172
//==============================================================================
173
// IncoherentInelasticAEDiscrete implementation
174
//==============================================================================
175

176
IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(
1,262✔
177
  hid_t group, const vector<double>& energy)
1,262✔
178
  : energy_ {energy}
1,262✔
179
{
180
  read_dataset(group, "energy_out", energy_out_);
1,262✔
181
  read_dataset(group, "mu_out", mu_out_);
1,262✔
182
  read_dataset(group, "skewed", skewed_);
1,262✔
183
}
1,262✔
184

185
void IncoherentInelasticAEDiscrete::sample_params(
123,683,823✔
186
  double E_in, double& E_out, int& j, uint64_t* seed) const
187
{
188
  // Get index and interpolation factor for inelastic grid
189
  int i;
123,683,823✔
190
  double f;
123,683,823✔
191
  get_energy_index(energy_, E_in, i, f);
123,683,823✔
192

193
  // Now that we have an incoming energy bin, we need to determine the outgoing
194
  // energy bin. This will depend on whether the outgoing energy distribution is
195
  // skewed. If it is skewed, then the first two and last two bins have lower
196
  // probabilities than the other bins (0.1 for the first and last bins and 0.4
197
  // for the second and second to last bins, relative to a normal bin
198
  // probability of 1). Otherwise, each bin is equally probable.
199

200
  int n = energy_out_.shape(1);
123,683,823!
201
  if (!skewed_) {
123,683,823!
202
    // All bins equally likely
UNCOV
203
    j = prn(seed) * n;
×
204
  } else {
205
    // Distribution skewed away from edge points
206
    double r = prn(seed) * (n - 3);
123,683,823✔
207
    if (r > 1.0) {
123,683,823✔
208
      // equally likely N-4 middle bins
209
      j = r + 1;
121,650,936✔
210
    } else if (r > 0.6) {
2,032,887✔
211
      // second to last bin has relative probability of 0.4
212
      j = n - 2;
806,605✔
213
    } else if (r > 0.5) {
1,226,282✔
214
      // last bin has relative probability of 0.1
215
      j = n - 1;
203,402✔
216
    } else if (r > 0.1) {
1,022,880✔
217
      // second bin has relative probability of 0.4
218
      j = 1;
819,261✔
219
    } else {
220
      // first bin has relative probability of 0.1
221
      j = 0;
203,619✔
222
    }
223
  }
224

225
  // Determine outgoing energy corresponding to E_in[i] and E_in[i+1]
226
  double E_ij = energy_out_(i, j);
123,683,823✔
227
  double E_i1j = energy_out_(i + 1, j);
123,683,823✔
228

229
  // Outgoing energy
230
  E_out = (1 - f) * E_ij + f * E_i1j;
123,683,823✔
231
}
123,683,823✔
232

233
void IncoherentInelasticAEDiscrete::sample(
123,683,823✔
234
  double E_in, double& E_out, double& mu, uint64_t* seed) const
235
{
236
  // Get index and interpolation factor for inelastic grid
237
  int i;
123,683,823✔
238
  double f;
123,683,823✔
239
  get_energy_index(energy_, E_in, i, f);
123,683,823✔
240

241
  int j;
123,683,823✔
242
  sample_params(E_in, E_out, j, seed);
123,683,823✔
243

244
  // Sample outgoing cosine bin
245
  int m = mu_out_.shape(2);
123,683,823!
246
  int k = prn(seed) * m;
123,683,823✔
247

248
  // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
249
  double mu_ijk = mu_out_(i, j, k);
123,683,823✔
250
  double mu_i1jk = mu_out_(i + 1, j, k);
123,683,823✔
251

252
  // Cosine of angle between incoming and outgoing neutron
253
  mu = (1 - f) * mu_ijk + f * mu_i1jk;
123,683,823✔
254
}
123,683,823✔
255

UNCOV
256
double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(
×
257
  double E_in, double mu, double& E_out, uint64_t* seed) const
258
{
259
  // Get index and interpolation factor for inelastic grid
UNCOV
260
  int i;
×
UNCOV
261
  double f;
×
UNCOV
262
  get_energy_index(energy_, E_in, i, f);
×
UNCOV
263
  int j;
×
UNCOV
264
  sample_params(E_in, E_out, j, seed);
×
265

UNCOV
266
  return get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all),
×
UNCOV
267
    mu_out_.slice(i + 1, j, tensor::all), f, mu);
×
268
}
269

270
//==============================================================================
271
// IncoherentInelasticAE implementation
272
//==============================================================================
273

274
IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)
44✔
275
{
276
  // Read correlated angle-energy distribution
277
  CorrelatedAngleEnergy dist {group};
44✔
278

279
  // Copy incident energies
280
  energy_ = dist.energy();
44✔
281

282
  // Convert to S(a,b) native format
283
  for (const auto& edist : dist.distribution()) {
176✔
284
    // Create temporary distribution
285
    DistEnergySab d;
132✔
286

287
    // Copy outgoing energy distribution
288
    d.n_e_out = edist.e_out.size();
132✔
289
    d.e_out = edist.e_out;
132✔
290
    d.e_out_pdf = edist.p;
132✔
291
    d.e_out_cdf = edist.c;
132✔
292

293
    for (int j = 0; j < d.n_e_out; ++j) {
660✔
294
      auto adist = dynamic_cast<Tabular*>(edist.angle[j].get());
528!
295
      if (adist) {
528!
296
        // On first pass, allocate space for angles
297
        if (j == 0) {
528✔
298
          auto n_mu = adist->x().size();
132✔
299
          d.mu = tensor::Tensor<double>({d.n_e_out, n_mu});
264✔
300
        }
301

302
        // Copy outgoing angles
303
        tensor::View<double> mu_j = d.mu.slice(j);
528✔
304
        std::copy(adist->x().begin(), adist->x().end(), mu_j.begin());
528✔
305
      }
528✔
306
    }
307

308
    distribution_.emplace_back(std::move(d));
132✔
309
  }
132✔
310
}
44✔
311

312
void IncoherentInelasticAE::sample_params(
3,614,897✔
313
  double E_in, double& E_out, double& f, int& l, int& j, uint64_t* seed) const
314
{
315
  // Get index and interpolation factor for inelastic grid
316
  int i;
3,614,897✔
317
  double f0;
3,614,897✔
318
  get_energy_index(energy_, E_in, i, f0);
3,614,897✔
319

320
  // Pick closer energy based on interpolation factor
321
  l = f0 > 0.5 ? i + 1 : i;
3,614,897✔
322

323
  // Determine outgoing energy bin
324
  // (First reset n_energy_out to the right value)
325
  int n = distribution_[l].n_e_out;
3,614,897✔
326
  double r1 = prn(seed);
3,614,897✔
327
  double c_j = distribution_[l].e_out_cdf[0];
3,614,897✔
328
  double c_j1;
3,614,897✔
329
  for (j = 0; j < n - 1; ++j) {
8,721,537!
330
    c_j1 = distribution_[l].e_out_cdf[j + 1];
8,721,537✔
331
    if (r1 < c_j1)
8,721,537✔
332
      break;
333
    c_j = c_j1;
5,106,640✔
334
  }
335

336
  // check to make sure j is <= n_energy_out - 2
337
  j = std::min(j, n - 2);
3,614,897!
338

339
  // Get the data to interpolate between
340
  double E_l_j = distribution_[l].e_out[j];
3,614,897!
341
  double p_l_j = distribution_[l].e_out_pdf[j];
3,614,897✔
342

343
  // Next part assumes linear-linear interpolation in standard
344
  double E_l_j1 = distribution_[l].e_out[j + 1];
3,614,897✔
345
  double p_l_j1 = distribution_[l].e_out_pdf[j + 1];
3,614,897✔
346

347
  // Find secondary energy (variable E)
348
  double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j);
3,614,897✔
349
  if (frac == 0.0) {
3,614,897!
350
    E_out = E_l_j + (r1 - c_j) / p_l_j;
3,614,897✔
351
  } else {
UNCOV
352
    E_out = E_l_j +
×
UNCOV
353
            (std::sqrt(std::max(0.0, p_l_j * p_l_j + 2.0 * frac * (r1 - c_j))) -
×
UNCOV
354
              p_l_j) /
×
355
              frac;
356
  }
357

358
  // Adjustment of outgoing energy
359
  double E_l = energy_[l];
3,614,897✔
360
  if (E_out < 0.5 * E_l) {
3,614,897✔
361
    E_out *= 2.0 * E_in / E_l - 1.0;
344,476✔
362
  } else {
363
    E_out += E_in - E_l;
3,270,421✔
364
  }
365

366
  f = (r1 - c_j) / (c_j1 - c_j);
3,614,897✔
367
}
3,614,897✔
368
void IncoherentInelasticAE::sample(
3,614,897✔
369
  double E_in, double& E_out, double& mu, uint64_t* seed) const
370
{
371
  double f;
3,614,897✔
372
  int l, j;
3,614,897✔
373
  sample_params(E_in, E_out, f, l, j, seed);
3,614,897✔
374

375
  // Sample outgoing cosine bin
376
  int n_mu = distribution_[l].mu.shape(1);
3,614,897!
377
  std::size_t k = prn(seed) * n_mu;
3,614,897✔
378

379
  // Rather than use the sampled discrete mu directly, it is smeared over
380
  // a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
381
  // discrete mu value itself.
382
  const auto& mu_l = distribution_[l].mu;
3,614,897✔
383

384
  // Interpolate kth mu value between distributions at energies j and j+1
385
  mu = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k));
3,614,897✔
386

387
  // Inteprolate (k-1)th mu value between distributions at energies j and j+1.
388
  // When k==0, pick a value that will smear the cosine out to a minimum of -1.
389
  double mu_left =
3,614,897✔
390
    (k == 0)
391
      ? mu_left = -1.0 - (mu + 1.0)
3,614,897✔
392
      : mu_left = mu_l(j, k - 1) + f * (mu_l(j + 1, k - 1) - mu_l(j, k - 1));
3,159,464✔
393

394
  // Inteprolate (k+1)th mu value between distributions at energies j and j+1.
395
  // When k is the last discrete value, pick a value that will smear the cosine
396
  // out to a maximum of 1.
397
  double mu_right =
3,614,897✔
398
    (k == n_mu - 1)
3,614,897✔
399
      ? mu_right = 1.0 + (1.0 - mu)
3,614,897✔
400
      : mu_right = mu_l(j, k + 1) + f * (mu_l(j + 1, k + 1) - mu_l(j, k + 1));
3,160,762✔
401

402
  // Smear cosine
403
  mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);
5,427,510✔
404
}
3,614,897✔
405

UNCOV
406
double IncoherentInelasticAE::sample_energy_and_pdf(
×
407
  double E_in, double mu, double& E_out, uint64_t* seed) const
408
{
UNCOV
409
  double f;
×
UNCOV
410
  int l, j;
×
UNCOV
411
  sample_params(E_in, E_out, f, l, j, seed);
×
412

UNCOV
413
  const auto& mu_l = distribution_[l].mu;
×
414

UNCOV
415
  return get_pdf_discrete_interpolated(
×
UNCOV
416
    mu_l.slice(j, tensor::all), mu_l.slice(j + 1, tensor::all), f, mu);
×
417
}
418

419
//==============================================================================
420
// MixedElasticAE implementation
421
//==============================================================================
422

423
MixedElasticAE::MixedElasticAE(
44✔
424
  hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs)
44✔
425
  : coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs)
44✔
426
{
427
  // Read incoherent elastic distribution
428
  hid_t incoherent_group = open_group(group, "incoherent");
44✔
429
  std::string temp;
44✔
430
  read_attribute(incoherent_group, "type", temp);
44✔
431
  if (temp == "incoherent_elastic") {
44!
UNCOV
432
    incoherent_dist_ = make_unique<IncoherentElasticAE>(incoherent_group);
×
433
  } else if (temp == "incoherent_elastic_discrete") {
44!
434
    auto xs = dynamic_cast<const Tabulated1D*>(&incoh_xs);
44✔
435
    incoherent_dist_ =
44✔
436
      make_unique<IncoherentElasticAEDiscrete>(incoherent_group, xs->x());
44✔
437
  }
438
  close_group(incoherent_group);
44✔
439
}
44✔
440

441
const AngleEnergy& MixedElasticAE::sample_dist(
46,486✔
442
  double E_in, uint64_t* seed) const
443
{
444
  // Evaluate coherent and incoherent elastic cross sections
445
  double xs_coh = coherent_xs_(E_in);
46,486✔
446
  double xs_incoh = incoherent_xs_(E_in);
46,486✔
447

448
  if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) {
46,486✔
449
    return coherent_dist_;
45,144✔
450
  } else {
451
    return *incoherent_dist_;
1,342✔
452
  }
453
}
454

455
void MixedElasticAE::sample(
46,486✔
456
  double E_in, double& E_out, double& mu, uint64_t* seed) const
457
{
458
  sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);
46,486✔
459
}
46,486✔
460

UNCOV
461
double MixedElasticAE::sample_energy_and_pdf(
×
462
  double E_in, double mu, double& E_out, uint64_t* seed) const
463
{
UNCOV
464
  return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
×
465
}
466

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