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

openmc-dev / openmc / 20470828925

23 Dec 2025 08:24PM UTC coverage: 81.898% (-0.04%) from 81.94%
20470828925

Pull #3550

github

web-flow
Merge b067e7f24 into 3f06a42ab
Pull Request #3550: [Point Detector] Add distribution get_pdf functionality

17078 of 23776 branches covered (71.83%)

Branch coverage included in aggregate %.

48 of 212 new or added lines in 13 files covered. (22.64%)

4 existing lines in 2 files now uncovered.

55208 of 64487 relevant lines covered (85.61%)

43370386.44 hits per line

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

75.9
/src/distribution_angle.cpp
1
#include "openmc/distribution_angle.h"
2

3
#include <cmath> // for abs, copysign
4

5
#include "xtensor/xarray.hpp"
6
#include "xtensor/xview.hpp"
7

8
#include "openmc/endf.h"
9
#include "openmc/hdf5_interface.h"
10
#include "openmc/math_functions.h"
11
#include "openmc/random_lcg.h"
12
#include "openmc/search.h"
13
#include "openmc/vector.h" // for vector
14

15
namespace openmc {
16

17
//==============================================================================
18
// AngleDistribution implementation
19
//==============================================================================
20

21
AngleDistribution::AngleDistribution(hid_t group)
2,963,266✔
22
{
23
  // Get incoming energies
24
  read_dataset(group, "energy", energy_);
2,963,266✔
25
  int n_energy = energy_.size();
2,963,266✔
26

27
  // Get outgoing energy distribution data
28
  vector<int> offsets;
2,963,266✔
29
  vector<int> interp;
2,963,266✔
30
  hid_t dset = open_dataset(group, "mu");
2,963,266✔
31
  read_attribute(dset, "offsets", offsets);
2,963,266✔
32
  read_attribute(dset, "interpolation", interp);
2,963,266✔
33
  xt::xarray<double> temp;
2,963,266✔
34
  read_dataset(dset, temp);
2,963,266✔
35
  close_dataset(dset);
2,963,266✔
36

37
  for (int i = 0; i < n_energy; ++i) {
28,475,883✔
38
    // Determine number of outgoing energies
39
    int j = offsets[i];
25,512,617✔
40
    int n;
41
    if (i < n_energy - 1) {
25,512,617✔
42
      n = offsets[i + 1] - j;
22,549,351✔
43
    } else {
44
      n = temp.shape()[1] - j;
2,963,266✔
45
    }
46

47
    // Create and initialize tabular distribution
48
    auto xs = xt::view(temp, 0, xt::range(j, j + n));
25,512,617✔
49
    auto ps = xt::view(temp, 1, xt::range(j, j + n));
25,512,617✔
50
    auto cs = xt::view(temp, 2, xt::range(j, j + n));
25,512,617✔
51
    vector<double> x {xs.begin(), xs.end()};
25,512,617✔
52
    vector<double> p {ps.begin(), ps.end()};
25,512,617✔
53
    vector<double> c {cs.begin(), cs.end()};
25,512,617✔
54

55
    // To get answers that match ACE data, for now we still use the tabulated
56
    // CDF values that were passed through to the HDF5 library. At a later
57
    // time, we can remove the CDF values from the HDF5 library and
58
    // reconstruct them using the PDF
59
    Tabular* mudist =
60
      new Tabular {x.data(), p.data(), n, int2interp(interp[i]), c.data()};
25,512,617✔
61

62
    distribution_.emplace_back(mudist);
25,512,617✔
63
  }
25,512,617✔
64
}
2,963,266✔
65

66
double AngleDistribution::sample(double E, uint64_t* seed) const
691,691,651✔
67
{
68
  // Find energy bin and calculate interpolation factor
69
  int i;
70
  double r;
71
  get_energy_index(energy_, E, i, r);
691,691,651✔
72

73
  // Sample between the ith and (i+1)th bin
74
  if (r > prn(seed))
691,691,651✔
75
    ++i;
280,574,831✔
76

77
  // Sample i-th distribution
78
  double mu = distribution_[i]->sample(seed);
691,691,651✔
79

80
  // Make sure mu is in range [-1,1] and return
81
  if (std::abs(mu) > 1.0)
691,691,651!
82
    mu = std::copysign(1.0, mu);
×
83
  return mu;
691,691,651✔
84
}
85

NEW
86
double AngleDistribution::evaluate(double E, double mu) const
×
87
{
88
  // Find energy bin and calculate interpolation factor
89
  int i;
90
  double r;
NEW
91
  get_energy_index(energy_, E, i, r);
×
92

NEW
93
  double pdf = 0.0;
×
NEW
94
  if (r > 0.0)
×
NEW
95
    pdf += r * distribution_[i + 1]->evaluate(mu);
×
NEW
96
  if (r < 1.0)
×
NEW
97
    pdf += (1.0 - r) * distribution_[i]->evaluate(mu);
×
NEW
98
  return pdf;
×
99
}
100

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