20197994533

Committed 13 Dec 2025 09:17PM UTC coverage: 82.153% (+0.004%) from 82.149%

Build # 20197994533

Build Type

Pull #3682

github

Committed by

web-flow

Commit Message

Merge 238b58cd8 into bbfa18d72

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

Coverage Stats

17023 of 23581 branches covered (72.19%)

Branch coverage included in aggregate %.

33 of 44 new or added lines in 4 files covered. (75.0%)

7 existing lines in 1 file now uncovered.

55123 of 64238 relevant lines covered (85.81%)

43459371.47 hits per line

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

82.76

/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 "xtensor/xarray.hpp"
#include "xtensor/xview.hpp"

#include "openmc/hdf5_interface.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);
    if (temp == "photon")
      is_photon_ = true;
  }
  // Open incoming energy dataset
  hid_t dset = open_dataset(group, "energy");

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

  auto temp_b = xt::view(temp, 0); // view of breakpoints
  auto temp_i = xt::view(temp, 1); // view of 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);

  xt::xarray<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 = xt::view(eout, 0, xt::range(j, j + n));
    d.p = xt::view(eout, 1, xt::range(j, j + n));
    d.c = xt::view(eout, 2, xt::range(j, j + n));
    d.r = xt::view(eout, 3, xt::range(j, j + n));
    d.a = xt::view(eout, 4, xt::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(
  double E_in, double& E_out, double& mu, uint64_t* seed) const
{
  // Find energy bin and calculate interpolation factor -- if the energy is
  // outside the range of the tabulated energies, choose the first or last bins
  auto n_energy_in = energy_.size();
  int i;
  double r;
  if (E_in < energy_[0]) {
    i = 0;
    r = 0.0;
  } else if (E_in > energy_[n_energy_in - 1]) {
    i = n_energy_in - 2;
    r = 1.0;
  } else {
    i = lower_bound_index(energy_.begin(), energy_.end(), E_in);
    r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);
  }

  // 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];
  double km_r, km_a;
  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);
    }
  }

  // Sampled correlated angle from Kalbach-Mann parameters
  if (is_photon_) {
    double T = uniform_distribution(0., 1., seed);
    double sinha = std::sinh(km_a);
    mu = sinha * ((1 + T) - 2 * km_r) /
         (sinha * T + std::cosh(km_a) - std::exp(km_a * T));
  } 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;
    }
  }
}

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

openmc-dev / openmc / 20197994533

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