17475377863

Committed 04 Sep 2025 08:05PM UTC coverage: 84.707% (-0.5%) from 85.209%

Build # 17475377863

Build Type

Pull #3550

github

Committed by

web-flow

Commit Message

Merge ff83fc4b9 into 591856472

Pull Request Pull Request #3550: [Point Detector] Add distribution get_pdf functionality

Run Details

0 of 370 new or added lines in 10 files covered. (0.0%)

3 existing lines in 1 file now uncovered.

52947 of 62506 relevant lines covered (84.71%)

38165206.92 hits per line

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

67.24

/src/distribution_angle.cpp

#include "openmc/distribution_angle.h"

#include <cmath> // for abs, copysign

#include "xtensor/xarray.hpp"
#include "xtensor/xview.hpp"

#include "openmc/endf.h"
#include "openmc/hdf5_interface.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/vector.h" // for vector

namespace openmc {

//==============================================================================
// AngleDistribution implementation
//==============================================================================

AngleDistribution::AngleDistribution(hid_t group)
{
  // Get incoming energies
  read_dataset(group, "energy", energy_);
  int n_energy = energy_.size();

  // Get outgoing energy distribution data
  vector<int> offsets;
  vector<int> interp;
  hid_t dset = open_dataset(group, "mu");
  read_attribute(dset, "offsets", offsets);
  read_attribute(dset, "interpolation", interp);
  xt::xarray<double> temp;
  read_dataset(dset, temp);
  close_dataset(dset);

  for (int i = 0; i < n_energy; ++i) {
    // Determine number of outgoing energies
    int j = offsets[i];
    int n;
    if (i < n_energy - 1) {
      n = offsets[i + 1] - j;
    } else {
      n = temp.shape()[1] - j;
    }

    // Create and initialize tabular distribution
    auto xs = xt::view(temp, 0, xt::range(j, j + n));
    auto ps = xt::view(temp, 1, xt::range(j, j + n));
    auto cs = xt::view(temp, 2, xt::range(j, j + n));
    vector<double> x {xs.begin(), xs.end()};
    vector<double> p {ps.begin(), ps.end()};
    vector<double> c {cs.begin(), cs.end()};

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

    distribution_.emplace_back(mudist);
  }
}

double AngleDistribution::sample(double E, uint64_t* seed) const
{
  // Determine number of incoming energies
  auto n = energy_.size();

  // Find energy bin and calculate interpolation factor -- if the energy is
  // outside the range of the tabulated energies, choose the first or last bins
  int i;
  double r;
  if (E < energy_[0]) {
    i = 0;
    r = 0.0;
  } else if (E > energy_[n - 1]) {
    i = n - 2;
    r = 1.0;
  } else {
    i = lower_bound_index(energy_.begin(), energy_.end(), E);
    r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);
  }

  // Sample between the ith and (i+1)th bin
  if (r > prn(seed))
    ++i;

  // Sample i-th distribution
  double mu = distribution_[i]->sample(seed);

  // Make sure mu is in range [-1,1] and return
  if (std::abs(mu) > 1.0)
    mu = std::copysign(1.0, mu);
  return mu;
}

double AngleDistribution::pdf(double E, double mu) const
{
  // Determine number of incoming energies
  auto n = energy_.size();

  // Find energy bin and calculate interpolation factor -- if the energy is
  // outside the range of the tabulated energies, choose the first or last bins
  int i;
  double r;
  if (E < energy_[0]) {
    i = 0;
    r = 0.0;
  } else if (E > energy_[n - 1]) {
    i = n - 2;
    r = 1.0;
  } else {
    i = lower_bound_index(energy_.begin(), energy_.end(), E);
    r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);
  }

  double pdf = 0.0;
  if (r > 0.0)
    pdf += r * distribution_[i + 1]->pdf(mu);
  if (r < 1.0)
    pdf += (1.0 - r) * distribution_[i]->pdf(mu);
  return pdf;
}

} // namespace openmc

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/random_lcg.h"
11	#include "openmc/search.h"
12	#include "openmc/vector.h" // for vector
13
14	namespace openmc {
15
16	//==============================================================================
17	// AngleDistribution implementation
18	//==============================================================================
19
20	AngleDistribution::AngleDistribution(hid_t group)	3,028,703✔
21	{
22	// Get incoming energies
23	read_dataset(group, "energy", energy_);	3,028,703✔
24	int n_energy = energy_.size();	3,028,703✔
25
26	// Get outgoing energy distribution data
27	vector<int> offsets;	3,028,703✔
28	vector<int> interp;	3,028,703✔
29	hid_t dset = open_dataset(group, "mu");	3,028,703✔
30	read_attribute(dset, "offsets", offsets);	3,028,703✔
31	read_attribute(dset, "interpolation", interp);	3,028,703✔
32	xt::xarray<double> temp;	3,028,703✔
33	read_dataset(dset, temp);	3,028,703✔
34	close_dataset(dset);	3,028,703✔
35
36	for (int i = 0; i < n_energy; ++i) {	28,709,620✔
37	// Determine number of outgoing energies
38	int j = offsets[i];	25,680,917✔
39	int n;
40	if (i < n_energy - 1) {	25,680,917✔
41	n = offsets[i + 1] - j;	22,652,214✔
42	} else {
43	n = temp.shape()[1] - j;	3,028,703✔
44	}
45
46	// Create and initialize tabular distribution
47	auto xs = xt::view(temp, 0, xt::range(j, j + n));	25,680,917✔
48	auto ps = xt::view(temp, 1, xt::range(j, j + n));	25,680,917✔
49	auto cs = xt::view(temp, 2, xt::range(j, j + n));	25,680,917✔
50	vector<double> x {xs.begin(), xs.end()};	25,680,917✔
51	vector<double> p {ps.begin(), ps.end()};	25,680,917✔
52	vector<double> c {cs.begin(), cs.end()};	25,680,917✔
53
54	// To get answers that match ACE data, for now we still use the tabulated
55	// CDF values that were passed through to the HDF5 library. At a later
56	// time, we can remove the CDF values from the HDF5 library and
57	// reconstruct them using the PDF
58	Tabular* mudist =
59	new Tabular {x.data(), p.data(), n, int2interp(interp[i]), c.data()};	25,680,917✔
60
61	distribution_.emplace_back(mudist);	25,680,917✔
62	}	25,680,917✔
63	}	3,028,703✔
64
65	double AngleDistribution::sample(double E, uint64_t* seed) const	600,692,901✔
66	{
67	// Determine number of incoming energies
68	auto n = energy_.size();	600,692,901✔
69
70	// Find energy bin and calculate interpolation factor -- if the energy is
71	// outside the range of the tabulated energies, choose the first or last bins
72	int i;
73	double r;
74	if (E < energy_[0]) {	600,692,901✔
75	i = 0;	336✔
76	r = 0.0;	336✔
77	} else if (E > energy_[n - 1]) {	600,692,565✔
78	i = n - 2;	×
79	r = 1.0;	×
80	} else {
81	i = lower_bound_index(energy_.begin(), energy_.end(), E);	600,692,565✔
82	r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);	600,692,565✔
83	}
84
85	// Sample between the ith and (i+1)th bin
86	if (r > prn(seed))	600,692,901✔
87	++i;	240,049,284✔
88
89	// Sample i-th distribution
90	double mu = distribution_[i]->sample(seed);	600,692,901✔
91
92	// Make sure mu is in range [-1,1] and return
93	if (std::abs(mu) > 1.0)	600,692,901✔
94	mu = std::copysign(1.0, mu);	×
95	return mu;	600,692,901✔
96	}
97
NEW 98	double AngleDistribution::pdf(double E, double mu) const	×
99	{
100	// Determine number of incoming energies
NEW 101	auto n = energy_.size();	×
102
103	// Find energy bin and calculate interpolation factor -- if the energy is
104	// outside the range of the tabulated energies, choose the first or last bins
105	int i;
106	double r;
NEW 107	if (E < energy_[0]) {	×
NEW 108	i = 0;	×
NEW 109	r = 0.0;	×
NEW 110	} else if (E > energy_[n - 1]) {	×
NEW 111	i = n - 2;	×
NEW 112	r = 1.0;	×
113	} else {
NEW 114	i = lower_bound_index(energy_.begin(), energy_.end(), E);	×
NEW 115	r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);	×
116	}
117
NEW 118	double pdf = 0.0;	×
NEW 119	if (r > 0.0)	×
NEW 120	pdf += r * distribution_[i + 1]->pdf(mu);	×
NEW 121	if (r < 1.0)	×
NEW 122	pdf += (1.0 - r) * distribution_[i]->pdf(mu);	×
NEW 123	return pdf;	×
124	}
125
126	} // namespace openmc

openmc-dev / openmc / 17475377863

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

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