20501343020

Committed 25 Dec 2025 07:38AM UTC coverage: 82.185% (+0.05%) from 82.139%

Build # 20501343020

Build Type

Pull #3692

github

Committed by

web-flow

Commit Message

Merge fc35e8018 into 3f06a42ab

Pull Request Pull Request #3692: Fix a bug in rotational periodic boundary conditions

Run Details

17096 of 23665 branches covered (72.24%)

Branch coverage included in aggregate %.

36 of 37 new or added lines in 3 files covered. (97.3%)

198 existing lines in 11 files now uncovered.

55209 of 64313 relevant lines covered (85.84%)

43488756.38 hits per line

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

91.84

/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/math_functions.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)
{
  // 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
  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];
  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 (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/math_functions.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)	67,490✔
25	{
26	// Open incoming energy dataset
27	hid_t dset = open_dataset(group, "energy");	67,490✔
28
29	// Get interpolation parameters
30	xt::xarray<int> temp;	67,490✔
31	read_attribute(dset, "interpolation", temp);	67,490✔
32
33	auto temp_b = xt::view(temp, 0); // view of breakpoints	67,490✔
34	auto temp_i = xt::view(temp, 1); // view of interpolation parameters	67,490✔
35
36	std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_));	67,490✔
37	for (const auto i : temp_i)	134,980✔
38	interpolation_.push_back(int2interp(i));	67,490✔
39	n_region_ = breakpoints_.size();	67,490✔
40
41	// Get incoming energies
42	read_dataset(dset, energy_);	67,490✔
43	std::size_t n_energy = energy_.size();	67,490✔
44	close_dataset(dset);	67,490✔
45
46	// Get outgoing energy distribution data
47	dset = open_dataset(group, "distribution");	67,490✔
48	vector<int> offsets;	67,490✔
49	vector<int> interp;	67,490✔
50	vector<int> n_discrete;	67,490✔
51	read_attribute(dset, "offsets", offsets);	67,490✔
52	read_attribute(dset, "interpolation", interp);	67,490✔
53	read_attribute(dset, "n_discrete_lines", n_discrete);	67,490✔
54
55	xt::xarray<double> eout;	67,490✔
56	read_dataset(dset, eout);	67,490✔
57	close_dataset(dset);	67,490✔
58
59	for (int i = 0; i < n_energy; ++i) {	1,374,973✔
60	// Determine number of outgoing energies
61	int j = offsets[i];	1,307,483✔
62	int n;
63	if (i < n_energy - 1) {	1,307,483✔
64	n = offsets[i + 1] - j;	1,239,993✔
65	} else {
66	n = eout.shape()[1] - j;	67,490✔
67	}
68
69	// Assign interpolation scheme and number of discrete lines
70	KMTable d;	1,307,483✔
71	d.interpolation = int2interp(interp[i]);	1,307,483✔
72	d.n_discrete = n_discrete[i];	1,307,483✔
73
74	// Copy data
75	d.e_out = xt::view(eout, 0, xt::range(j, j + n));	1,307,483✔
76	d.p = xt::view(eout, 1, xt::range(j, j + n));	1,307,483✔
77	d.c = xt::view(eout, 2, xt::range(j, j + n));	1,307,483✔
78	d.r = xt::view(eout, 3, xt::range(j, j + n));	1,307,483✔
79	d.a = xt::view(eout, 4, xt::range(j, j + n));	1,307,483✔
80
81	// To get answers that match ACE data, for now we still use the tabulated
82	// CDF values that were passed through to the HDF5 library. At a later
83	// time, we can remove the CDF values from the HDF5 library and
84	// reconstruct them using the PDF
85	if (false) {
86	// Calculate cumulative distribution function -- discrete portion
87	for (int k = 0; k < d.n_discrete; ++k) {
88	if (k == 0) {
89	d.c[k] = d.p[k];
90	} else {
91	d.c[k] = d.c[k - 1] + d.p[k];
92	}
93	}
94
95	// Continuous portion
96	for (int k = d.n_discrete; k < n; ++k) {
97	if (k == d.n_discrete) {
98	d.c[k] = d.c[k - 1] + d.p[k];
99	} else {
100	if (d.interpolation == Interpolation::histogram) {
101	d.c[k] = d.c[k - 1] + d.p[k - 1] * (d.e_out[k] - d.e_out[k - 1]);
102	} else if (d.interpolation == Interpolation::lin_lin) {
103	d.c[k] = d.c[k - 1] + 0.5 * (d.p[k - 1] + d.p[k]) *
104	(d.e_out[k] - d.e_out[k - 1]);
105	}
106	}
107	}
108
109	// Normalize density and distribution functions
110	d.p /= d.c[n - 1];
111	d.c /= d.c[n - 1];
112	}
113
114	distribution_.push_back(std::move(d));	1,307,483✔
115	} // incoming energies	1,307,483✔
116	}	67,490✔
117
118	void KalbachMann::sample(	4,985,796✔
119	double E_in, double& E_out, double& mu, uint64_t* seed) const
120	{
121	// Find energy bin and calculate interpolation factor
122	int i;
123	double r;
124	get_energy_index(energy_, E_in, i, r);	4,985,796✔
125
126	// Sample between the ith and [i+1]th bin
127	int l = r > prn(seed) ? i + 1 : i;	4,985,796✔
128
129	// Interpolation for energy E1 and EK
130	int n_energy_out = distribution_[i].e_out.size();	4,985,796✔
131	int n_discrete = distribution_[i].n_discrete;	4,985,796✔
132	double E_i_1 = distribution_[i].e_out[n_discrete];	4,985,796✔
133	double E_i_K = distribution_[i].e_out[n_energy_out - 1];	4,985,796✔
134
135	n_energy_out = distribution_[i + 1].e_out.size();	4,985,796✔
136	n_discrete = distribution_[i + 1].n_discrete;	4,985,796✔
137	double E_i1_1 = distribution_[i + 1].e_out[n_discrete];	4,985,796✔
138	double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];	4,985,796✔
139
140	double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);	4,985,796✔
141	double E_K = E_i_K + r * (E_i1_K - E_i_K);	4,985,796✔
142
143	// Determine outgoing energy bin
144	n_energy_out = distribution_[l].e_out.size();	4,985,796✔
145	n_discrete = distribution_[l].n_discrete;	4,985,796✔
146	double r1 = prn(seed);	4,985,796✔
147	double c_k = distribution_[l].c[0];	4,985,796✔
148	int k = 0;	4,985,796✔
149	int end = n_energy_out - 2;	4,985,796✔
150
151	// Discrete portion
152	for (int j = 0; j < n_discrete; ++j) {	4,985,796!
UNCOV 153	k = j;	×
UNCOV 154	c_k = distribution_[l].c[k];	×
UNCOV 155	if (r1 < c_k) {	×
UNCOV 156	end = j;	×
UNCOV 157	break;	×
158	}
159	}
160
161	// Continuous portion
162	double c_k1;
163	for (int j = n_discrete; j < end; ++j) {	125,233,746✔
164	k = j;	125,040,331✔
165	c_k1 = distribution_[l].c[k + 1];	125,040,331✔
166	if (r1 < c_k1)	125,040,331✔
167	break;	4,792,381✔
168	k = j + 1;	120,247,950✔
169	c_k = c_k1;	120,247,950✔
170	}
171
172	double E_l_k = distribution_[l].e_out[k];	4,985,796✔
173	double p_l_k = distribution_[l].p[k];	4,985,796✔
174	double km_r, km_a;
175	if (distribution_[l].interpolation == Interpolation::histogram) {	4,985,796✔
176	// Histogram interpolation
177	if (p_l_k > 0.0 && k >= n_discrete) {	4,966,040!
178	E_out = E_l_k + (r1 - c_k) / p_l_k;	4,966,040✔
179	} else {
UNCOV 180	E_out = E_l_k;	×
181	}
182
183	// Determine Kalbach-Mann parameters
184	km_r = distribution_[l].r[k];	4,966,040✔
185	km_a = distribution_[l].a[k];	4,966,040✔
186
187	} else {
188	// Linear-linear interpolation
189	double E_l_k1 = distribution_[l].e_out[k + 1];	19,756✔
190	double p_l_k1 = distribution_[l].p[k + 1];	19,756✔
191
192	double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);	19,756✔
193	if (frac == 0.0) {	19,756!
UNCOV 194	E_out = E_l_k + (r1 - c_k) / p_l_k;	×
195	} else {
196	E_out =	19,756✔
197	E_l_k +	19,756✔
198	(std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -	19,756✔
199	p_l_k) /	19,756✔
200	frac;
201	}
202
203	// Determine Kalbach-Mann parameters
204	km_r = distribution_[l].r[k] +	19,756✔
205	(E_out - E_l_k) / (E_l_k1 - E_l_k) *	39,512✔
206	(distribution_[l].r[k + 1] - distribution_[l].r[k]);	19,756✔
207	km_a = distribution_[l].a[k] +	19,756✔
208	(E_out - E_l_k) / (E_l_k1 - E_l_k) *	39,512✔
209	(distribution_[l].a[k + 1] - distribution_[l].a[k]);	19,756✔
210	}
211
212	// Now interpolate between incident energy bins i and i + 1
213	if (k >= n_discrete) {	4,985,796!
214	if (l == i) {	4,985,796✔
215	E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);	2,480,119✔
216	} else {
217	E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);	2,505,677✔
218	}
219	}
220
221	// Sampled correlated angle from Kalbach-Mann parameters
222	if (prn(seed) > km_r) {	4,985,796✔
223	double T = uniform_distribution(-1., 1., seed) * std::sinh(km_a);	4,906,208✔
224	mu = std::log(T + std::sqrt(T * T + 1.0)) / km_a;	4,906,208✔
225	} else {
226	double r1 = prn(seed);	79,588✔
227	mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;	79,588✔
228	}
229	}	4,985,796✔
230
231	} // namespace openmc

openmc-dev / openmc / 20501343020

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