30839048241

Committed 03 Aug 2026 05:57PM UTC coverage: 81.424% (-0.001%) from 81.425%

Build # 30839048241

Build Type

Pull #3682

github

Committed by

web-flow

Commit Message

Merge dafb00d5f into 8202ef6fb

Pull Request Pull Request #3682: Extend Kalbach Mann systematics to support incident photons

Coverage Stats

18530 of 26853 branches covered (69.01%)

Branch coverage included in aggregate %.

33 of 45 new or added lines in 4 files covered. (73.33%)

60341 of 70012 relevant lines covered (86.19%)

50042705.93 hits per line

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

79.91

/src/secondary_kalbach.cpp

#include "openmc/secondary_kalbach.h"

#include <algorithm> // for copy, move
#include <cmath>     // for log, sqrt, sinh
#include <cstddef>   // for size_t
#include <iterator>  // for back_inserter

#include "openmc/tensor.h"

#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/particle_type.h"
#include "openmc/random_dist.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/vector.h"

namespace openmc {

//==============================================================================
//! KalbachMann implementation
//==============================================================================

KalbachMann::KalbachMann(hid_t group)
{
  is_photon_ = false;
  // Check if projectile is a photon
  if (attribute_exists(group, "particle")) {
    std::string temp;
    read_attribute(group, "particle", temp);
    auto type = ParticleType(temp);
    if (type.is_photon())
      is_photon_ = true;
  }
  // Open incoming energy dataset
  hid_t dset = open_dataset(group, "energy");

  // Get interpolation parameters
  tensor::Tensor<int> temp;
  read_attribute(dset, "interpolation", temp);

  tensor::View<int> temp_b = temp.slice(0); // breakpoints
  tensor::View<int> temp_i = temp.slice(1); // interpolation parameters

  std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_));
  for (const auto i : temp_i)
    interpolation_.push_back(int2interp(i));
  n_region_ = breakpoints_.size();

  // Get incoming energies
  read_dataset(dset, energy_);
  std::size_t n_energy = energy_.size();
  close_dataset(dset);

  // Get outgoing energy distribution data
  dset = open_dataset(group, "distribution");
  vector<int> offsets;
  vector<int> interp;
  vector<int> n_discrete;
  read_attribute(dset, "offsets", offsets);
  read_attribute(dset, "interpolation", interp);
  read_attribute(dset, "n_discrete_lines", n_discrete);

  tensor::Tensor<double> eout;
  read_dataset(dset, eout);
  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 = eout.shape(1) - j;
    }

    // Assign interpolation scheme and number of discrete lines
    KMTable d;
    d.interpolation = int2interp(interp[i]);
    d.n_discrete = n_discrete[i];

    // Copy data
    d.e_out = eout.slice(0, tensor::range(j, j + n));
    d.p = eout.slice(1, tensor::range(j, j + n));
    d.c = eout.slice(2, tensor::range(j, j + n));
    d.r = eout.slice(3, tensor::range(j, j + n));
    d.a = eout.slice(4, tensor::range(j, j + n));

    // 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
    if (false) {
      // Calculate cumulative distribution function -- discrete portion
      for (int k = 0; k < d.n_discrete; ++k) {
        if (k == 0) {
          d.c[k] = d.p[k];
        } else {
          d.c[k] = d.c[k - 1] + d.p[k];
        }
      }

      // Continuous portion
      for (int k = d.n_discrete; k < n; ++k) {
        if (k == d.n_discrete) {
          d.c[k] = d.c[k - 1] + d.p[k];
        } else {
          if (d.interpolation == Interpolation::histogram) {
            d.c[k] = d.c[k - 1] + d.p[k - 1] * (d.e_out[k] - d.e_out[k - 1]);
          } else if (d.interpolation == Interpolation::lin_lin) {
            d.c[k] = d.c[k - 1] + 0.5 * (d.p[k - 1] + d.p[k]) *
                                    (d.e_out[k] - d.e_out[k - 1]);
          }
        }
      }

      // Normalize density and distribution functions
      d.p /= d.c[n - 1];
      d.c /= d.c[n - 1];
    }

    distribution_.push_back(std::move(d));
  } // incoming energies
}

void KalbachMann::sample_params(
  double E_in, double& E_out, double& km_a, double& km_r, uint64_t* seed) const
{
  // Find energy bin and calculate interpolation factor
  int i;
  double r;
  get_energy_index(energy_, E_in, i, r);

  // Sample between the ith and [i+1]th bin
  int l = r > prn(seed) ? i + 1 : i;

  // Interpolation for energy E1 and EK
  int n_energy_out = distribution_[i].e_out.size();
  int n_discrete = distribution_[i].n_discrete;
  double E_i_1 = distribution_[i].e_out[n_discrete];
  double E_i_K = distribution_[i].e_out[n_energy_out - 1];

  n_energy_out = distribution_[i + 1].e_out.size();
  n_discrete = distribution_[i + 1].n_discrete;
  double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
  double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];

  double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
  double E_K = E_i_K + r * (E_i1_K - E_i_K);

  // Determine outgoing energy bin
  n_energy_out = distribution_[l].e_out.size();
  n_discrete = distribution_[l].n_discrete;
  double r1 = prn(seed);
  double c_k = distribution_[l].c[0];
  int k = 0;
  int end = n_energy_out - 2;

  // Discrete portion
  for (int j = 0; j < n_discrete; ++j) {
    k = j;
    c_k = distribution_[l].c[k];
    if (r1 < c_k) {
      end = j;
      break;
    }
  }

  // Continuous portion
  double c_k1;
  for (int j = n_discrete; j < end; ++j) {
    k = j;
    c_k1 = distribution_[l].c[k + 1];
    if (r1 < c_k1)
      break;
    k = j + 1;
    c_k = c_k1;
  }

  double E_l_k = distribution_[l].e_out[k];
  double p_l_k = distribution_[l].p[k];
  if (distribution_[l].interpolation == Interpolation::histogram) {
    // Histogram interpolation
    if (p_l_k > 0.0 && k >= n_discrete) {
      E_out = E_l_k + (r1 - c_k) / p_l_k;
    } else {
      E_out = E_l_k;
    }

    // Determine Kalbach-Mann parameters
    km_r = distribution_[l].r[k];
    km_a = distribution_[l].a[k];

  } else {
    // Linear-linear interpolation
    double E_l_k1 = distribution_[l].e_out[k + 1];
    double p_l_k1 = distribution_[l].p[k + 1];

    double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
    if (frac == 0.0) {
      E_out = E_l_k + (r1 - c_k) / p_l_k;
    } else {
      E_out =
        E_l_k +
        (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
          p_l_k) /
          frac;
    }

    // Determine Kalbach-Mann parameters
    km_r = distribution_[l].r[k] +
           (E_out - E_l_k) / (E_l_k1 - E_l_k) *
             (distribution_[l].r[k + 1] - distribution_[l].r[k]);
    km_a = distribution_[l].a[k] +
           (E_out - E_l_k) / (E_l_k1 - E_l_k) *
             (distribution_[l].a[k + 1] - distribution_[l].a[k]);
  }

  // Now interpolate between incident energy bins i and i + 1
  if (k >= n_discrete) {
    if (l == i) {
      E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
    } else {
      E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
    }
  }
}

void KalbachMann::sample(
  double E_in, double& E_out, double& mu, uint64_t* seed) const
{
  double km_r, km_a;
  sample_params(E_in, E_out, km_a, km_r, seed);

  // Sampled correlated angle from Kalbach-Mann parameters
  if (is_photon_) {
    if (prn(seed) > km_r) {
      double T = uniform_distribution(-1., 1., seed);
      mu = std::log(std::cosh(km_a) + T * std::sinh(km_a)) / km_a;
    } else {
      mu = uniform_distribution(-1., 1., seed);
    }
  } else {
    if (prn(seed) > km_r) {
      double T = uniform_distribution(-1., 1., seed) * std::sinh(km_a);
      mu = std::log(T + std::sqrt(T * T + 1.0)) / km_a;
    } else {
      double r1 = prn(seed);
      mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;
    }
  }
}
double KalbachMann::sample_energy_and_pdf(
  double E_in, double mu, double& E_out, uint64_t* seed) const
{
  double km_r, km_a;
  sample_params(E_in, E_out, km_a, km_r, seed);

  // https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle
  return km_a / (2 * std::sinh(km_a)) *
         (std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu));
}

} // namespace openmc

1	#include "openmc/secondary_kalbach.h"
2
3	#include <algorithm> // for copy, move
4	#include <cmath> // for log, sqrt, sinh
5	#include <cstddef> // for size_t
6	#include <iterator> // for back_inserter
7
8	#include "openmc/tensor.h"
9
10	#include "openmc/hdf5_interface.h"
11	#include "openmc/math_functions.h"
12	#include "openmc/particle_type.h"
13	#include "openmc/random_dist.h"
14	#include "openmc/random_lcg.h"
15	#include "openmc/search.h"
16	#include "openmc/vector.h"
17
18	namespace openmc {
19
20	//==============================================================================
21	//! KalbachMann implementation
22	//==============================================================================
23
24	KalbachMann::KalbachMann(hid_t group)	78,047✔
25	{
26	is_photon_ = false;	78,047✔
27	// Check if projectile is a photon
28	if (attribute_exists(group, "particle")) {	78,047!
NEW 29	std::string temp;	×
NEW 30	read_attribute(group, "particle", temp);	×
NEW 31	auto type = ParticleType(temp);	×
NEW 32	if (type.is_photon())	×
NEW 33	is_photon_ = true;	×
NEW 34	}	×
35	// Open incoming energy dataset
36	hid_t dset = open_dataset(group, "energy");	78,047✔
37
38	// Get interpolation parameters
39	tensor::Tensor<int> temp;	78,047✔
40	read_attribute(dset, "interpolation", temp);	78,047✔
41
42	tensor::View<int> temp_b = temp.slice(0); // breakpoints	78,047✔
43	tensor::View<int> temp_i = temp.slice(1); // interpolation parameters	78,047✔
44
45	std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_));	78,047✔
46	for (const auto i : temp_i)	156,094✔
47	interpolation_.push_back(int2interp(i));	78,047✔
48	n_region_ = breakpoints_.size();	78,047✔
49
50	// Get incoming energies
51	read_dataset(dset, energy_);	78,047✔
52	std::size_t n_energy = energy_.size();	78,047✔
53	close_dataset(dset);	78,047✔
54
55	// Get outgoing energy distribution data
56	dset = open_dataset(group, "distribution");	78,047✔
57	vector<int> offsets;	78,047✔
58	vector<int> interp;	78,047✔
59	vector<int> n_discrete;	78,047✔
60	read_attribute(dset, "offsets", offsets);	78,047✔
61	read_attribute(dset, "interpolation", interp);	78,047✔
62	read_attribute(dset, "n_discrete_lines", n_discrete);	78,047✔
63
64	tensor::Tensor<double> eout;	78,047✔
65	read_dataset(dset, eout);	78,047✔
66	close_dataset(dset);	78,047✔
67
68	for (int i = 0; i < n_energy; ++i) {	1,564,018✔
69	// Determine number of outgoing energies
70	int j = offsets[i];	1,485,971✔
71	int n;	1,485,971✔
72	if (i < n_energy - 1) {	1,485,971✔
73	n = offsets[i + 1] - j;	1,407,924✔
74	} else {
75	n = eout.shape(1) - j;	156,094!
76	}
77
78	// Assign interpolation scheme and number of discrete lines
79	KMTable d;	1,485,971✔
80	d.interpolation = int2interp(interp[i]);	1,485,971✔
81	d.n_discrete = n_discrete[i];	1,485,971✔
82
83	// Copy data
84	d.e_out = eout.slice(0, tensor::range(j, j + n));	1,485,971✔
85	d.p = eout.slice(1, tensor::range(j, j + n));	1,485,971✔
86	d.c = eout.slice(2, tensor::range(j, j + n));	1,485,971✔
87	d.r = eout.slice(3, tensor::range(j, j + n));	1,485,971✔
88	d.a = eout.slice(4, tensor::range(j, j + n));	1,485,971✔
89
90	// To get answers that match ACE data, for now we still use the tabulated
91	// CDF values that were passed through to the HDF5 library. At a later
92	// time, we can remove the CDF values from the HDF5 library and
93	// reconstruct them using the PDF
94	if (false) {	1,485,971✔
95	// Calculate cumulative distribution function -- discrete portion
96	for (int k = 0; k < d.n_discrete; ++k) {
97	if (k == 0) {
98	d.c[k] = d.p[k];
99	} else {
100	d.c[k] = d.c[k - 1] + d.p[k];
101	}
102	}
103
104	// Continuous portion
105	for (int k = d.n_discrete; k < n; ++k) {
106	if (k == d.n_discrete) {
107	d.c[k] = d.c[k - 1] + d.p[k];
108	} else {
109	if (d.interpolation == Interpolation::histogram) {
110	d.c[k] = d.c[k - 1] + d.p[k - 1] * (d.e_out[k] - d.e_out[k - 1]);
111	} else if (d.interpolation == Interpolation::lin_lin) {
112	d.c[k] = d.c[k - 1] + 0.5 * (d.p[k - 1] + d.p[k]) *
113	(d.e_out[k] - d.e_out[k - 1]);
114	}
115	}
116	}
117
118	// Normalize density and distribution functions
119	d.p /= d.c[n - 1];
120	d.c /= d.c[n - 1];
121	}
122
123	distribution_.push_back(std::move(d));	1,485,971✔
124	} // incoming energies	1,485,971✔
125	}	312,188✔
126
127	void KalbachMann::sample_params(	6,863,296✔
128	double E_in, double& E_out, double& km_a, double& km_r, uint64_t* seed) const
129	{
130	// Find energy bin and calculate interpolation factor
131	int i;	6,863,296✔
132	double r;	6,863,296✔
133	get_energy_index(energy_, E_in, i, r);	6,863,296✔
134
135	// Sample between the ith and [i+1]th bin
136	int l = r > prn(seed) ? i + 1 : i;	6,863,296✔
137
138	// Interpolation for energy E1 and EK
139	int n_energy_out = distribution_[i].e_out.size();	6,863,296✔
140	int n_discrete = distribution_[i].n_discrete;	6,863,296✔
141	double E_i_1 = distribution_[i].e_out[n_discrete];	6,863,296✔
142	double E_i_K = distribution_[i].e_out[n_energy_out - 1];	6,863,296✔
143
144	n_energy_out = distribution_[i + 1].e_out.size();	6,863,296✔
145	n_discrete = distribution_[i + 1].n_discrete;	6,863,296✔
146	double E_i1_1 = distribution_[i + 1].e_out[n_discrete];	6,863,296✔
147	double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];	6,863,296✔
148
149	double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);	6,863,296✔
150	double E_K = E_i_K + r * (E_i1_K - E_i_K);	6,863,296✔
151
152	// Determine outgoing energy bin
153	n_energy_out = distribution_[l].e_out.size();	6,863,296✔
154	n_discrete = distribution_[l].n_discrete;	6,863,296✔
155	double r1 = prn(seed);	6,863,296✔
156	double c_k = distribution_[l].c[0];	6,863,296✔
157	int k = 0;	6,863,296✔
158	int end = n_energy_out - 2;	6,863,296✔
159
160	// Discrete portion
161	for (int j = 0; j < n_discrete; ++j) {	6,863,296!
162	k = j;	×
163	c_k = distribution_[l].c[k];	×
164	if (r1 < c_k) {	×
165	end = j;
166	break;
167	}
168	}
169
170	// Continuous portion
171	double c_k1;	6,863,296✔
172	for (int j = n_discrete; j < end; ++j) {	193,089,975✔
173	k = j;	192,857,084✔
174	c_k1 = distribution_[l].c[k + 1];	192,857,084✔
175	if (r1 < c_k1)	192,857,084✔
176	break;
177	k = j + 1;
178	c_k = c_k1;
179	}
180
181	double E_l_k = distribution_[l].e_out[k];	6,863,296✔
182	double p_l_k = distribution_[l].p[k];	6,863,296✔
183	if (distribution_[l].interpolation == Interpolation::histogram) {	6,863,296✔
184	// Histogram interpolation
185	if (p_l_k > 0.0 && k >= n_discrete) {	6,843,540!
186	E_out = E_l_k + (r1 - c_k) / p_l_k;	6,843,540✔
187	} else {
188	E_out = E_l_k;	×
189	}
190
191	// Determine Kalbach-Mann parameters
192	km_r = distribution_[l].r[k];	6,843,540✔
193	km_a = distribution_[l].a[k];	6,843,540✔
194
195	} else {
196	// Linear-linear interpolation
197	double E_l_k1 = distribution_[l].e_out[k + 1];	19,756!
198	double p_l_k1 = distribution_[l].p[k + 1];	19,756✔
199
200	double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);	19,756✔
201	if (frac == 0.0) {	19,756!
202	E_out = E_l_k + (r1 - c_k) / p_l_k;	×
203	} else {
204	E_out =	39,512✔
205	E_l_k +	19,756✔
206	(std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -	39,512!
207	p_l_k) /	19,756✔
208	frac;
209	}
210
211	// Determine Kalbach-Mann parameters
212	km_r = distribution_[l].r[k] +	19,756✔
213	(E_out - E_l_k) / (E_l_k1 - E_l_k) *	19,756✔
214	(distribution_[l].r[k + 1] - distribution_[l].r[k]);	19,756✔
215	km_a = distribution_[l].a[k] +	19,756✔
216	(E_out - E_l_k) / (E_l_k1 - E_l_k) *	19,756✔
217	(distribution_[l].a[k + 1] - distribution_[l].a[k]);	19,756✔
218	}
219
220	// Now interpolate between incident energy bins i and i + 1
221	if (k >= n_discrete) {	6,863,296!
222	if (l == i) {	6,863,296✔
223	E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);	3,430,491✔
224	} else {
225	E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);	3,432,805✔
226	}
227	}
228	}	6,863,296✔
229
230	void KalbachMann::sample(	6,863,296✔
231	double E_in, double& E_out, double& mu, uint64_t* seed) const
232	{
233	double km_r, km_a;	6,863,296✔
234	sample_params(E_in, E_out, km_a, km_r, seed);	6,863,296✔
235
236	// Sampled correlated angle from Kalbach-Mann parameters
237	if (is_photon_) {	6,863,296!
NEW 238	if (prn(seed) > km_r) {	×
NEW 239	double T = uniform_distribution(-1., 1., seed);	×
NEW 240	mu = std::log(std::cosh(km_a) + T * std::sinh(km_a)) / km_a;	×
241	} else {
NEW 242	mu = uniform_distribution(-1., 1., seed);	×
243	}
244	} else {
245	if (prn(seed) > km_r) {	6,863,296✔
246	double T = uniform_distribution(-1., 1., seed) * std::sinh(km_a);	6,717,006✔
247	mu = std::log(T + std::sqrt(T * T + 1.0)) / km_a;	6,717,006✔
248	} else {
249	double r1 = prn(seed);	146,290✔
250	mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;	146,290✔
251	}
252	}
253	}	6,863,296✔
254	double KalbachMann::sample_energy_and_pdf(	×
255	double E_in, double mu, double& E_out, uint64_t* seed) const
256	{
257	double km_r, km_a;	×
258	sample_params(E_in, E_out, km_a, km_r, seed);	×
259
260	// https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle
261	return km_a / (2 * std::sinh(km_a)) *	×
262	(std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu));	×
263	}
264
265	} // namespace openmc

openmc-dev / openmc / 30839048241

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