17271733860

Committed 27 Aug 2025 03:46PM UTC coverage: 84.486% (-0.7%) from 85.204%

Build # 17271733860

Build Type

Pull #3550

github

Committed by

web-flow

Commit Message

Merge 06ae298be into d1df80a21

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

Run Details

4 of 551 new or added lines in 11 files covered. (0.73%)

4 existing lines in 2 files now uncovered.

52940 of 62661 relevant lines covered (84.49%)

37443356.37 hits per line

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

55.19

/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/nuclide.h"
#include "openmc/particle.h"
#include "openmc/random_dist.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/tallies/tally_scoring.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 -- 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 (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;
  }
}
void KalbachMann::get_pdf(double det_pos[4], double E_in, double& E_out,
  uint64_t* seed, Particle& p, std::vector<double>& mu_cm,
  std::vector<double>& Js, std::vector<Particle>& ghost_particles,
  std::vector<double>& pdfs_lab) 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
  const auto& nuc {data::nuclides[p.event_nuclide()]};

  // double E_out;
  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;
    }
    //  Linear-linear interpolation
    // 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);
    }
  }
  /* function to transform pdf CM to lab frame - Future implementation
    get_pdf_to_point_elastic(det_pos, p, mu_cm, Js, ghost_particles, E_out /
    1e6); for (std::size_t i = 0; i < mu_cm.size(); ++i) {
      // Assuming Js.size() is the same as mu_cm.size()
      double mu_c = mu_cm[i];
      double derivative = Js[i];
      double pdf_cm = km_a / (2 * std::sinh(km_a)) *
                      (std::cosh(km_a * mu_c) +
                        km_r * std::sinh(km_a * mu_c)); // center of mass
      pdfs_lab.push_back(pdf_cm / std::abs(derivative));
    }
   */

  // https://docs.openmc.org/en/v0.8.0/methods/physics.html#equation-KM-pdf-angle
  // double pdf_mu = km_a / (2 * std::sinh(km_a)) * (std::cosh(km_a * mymu) +
  // km_r * std::sinh(km_a * mymu)); // center of mass return pdf_mu;

  const auto& rx {nuc->reactions_[p.event_index_mt()]};
  if (!rx->scatter_in_cm_) {
    fatal_error("didn't implement lab");
  }
}

} // 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/nuclide.h"
13	#include "openmc/particle.h"
14	#include "openmc/random_dist.h"
15	#include "openmc/random_lcg.h"
16	#include "openmc/search.h"
17	#include "openmc/tallies/tally_scoring.h"
18	#include "openmc/vector.h"
19
20	namespace openmc {
21
22	//==============================================================================
23	//! KalbachMann implementation
24	//==============================================================================
25
26	KalbachMann::KalbachMann(hid_t group)	59,993✔
27	{
28	// Open incoming energy dataset
29	hid_t dset = open_dataset(group, "energy");	59,993✔
30
31	// Get interpolation parameters
32	xt::xarray<int> temp;	59,993✔
33	read_attribute(dset, "interpolation", temp);	59,993✔
34
35	auto temp_b = xt::view(temp, 0); // view of breakpoints	59,993✔
36	auto temp_i = xt::view(temp, 1); // view of interpolation parameters	59,993✔
37
38	std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_));	59,993✔
39	for (const auto i : temp_i)	119,986✔
40	interpolation_.push_back(int2interp(i));	59,993✔
41	n_region_ = breakpoints_.size();	59,993✔
42
43	// Get incoming energies
44	read_dataset(dset, energy_);	59,993✔
45	std::size_t n_energy = energy_.size();	59,993✔
46	close_dataset(dset);	59,993✔
47
48	// Get outgoing energy distribution data
49	dset = open_dataset(group, "distribution");	59,993✔
50	vector<int> offsets;	59,993✔
51	vector<int> interp;	59,993✔
52	vector<int> n_discrete;	59,993✔
53	read_attribute(dset, "offsets", offsets);	59,993✔
54	read_attribute(dset, "interpolation", interp);	59,993✔
55	read_attribute(dset, "n_discrete_lines", n_discrete);	59,993✔
56
57	xt::xarray<double> eout;	59,993✔
58	read_dataset(dset, eout);	59,993✔
59	close_dataset(dset);	59,993✔
60
61	for (int i = 0; i < n_energy; ++i) {	1,219,057✔
62	// Determine number of outgoing energies
63	int j = offsets[i];	1,159,064✔
64	int n;
65	if (i < n_energy - 1) {	1,159,064✔
66	n = offsets[i + 1] - j;	1,099,071✔
67	} else {
68	n = eout.shape()[1] - j;	59,993✔
69	}
70
71	// Assign interpolation scheme and number of discrete lines
72	KMTable d;	1,159,064✔
73	d.interpolation = int2interp(interp[i]);	1,159,064✔
74	d.n_discrete = n_discrete[i];	1,159,064✔
75
76	// Copy data
77	d.e_out = xt::view(eout, 0, xt::range(j, j + n));	1,159,064✔
78	d.p = xt::view(eout, 1, xt::range(j, j + n));	1,159,064✔
79	d.c = xt::view(eout, 2, xt::range(j, j + n));	1,159,064✔
80	d.r = xt::view(eout, 3, xt::range(j, j + n));	1,159,064✔
81	d.a = xt::view(eout, 4, xt::range(j, j + n));	1,159,064✔
82
83	// To get answers that match ACE data, for now we still use the tabulated
84	// CDF values that were passed through to the HDF5 library. At a later
85	// time, we can remove the CDF values from the HDF5 library and
86	// reconstruct them using the PDF
87	if (false) {
88	// Calculate cumulative distribution function -- discrete portion
89	for (int k = 0; k < d.n_discrete; ++k) {
90	if (k == 0) {
91	d.c[k] = d.p[k];
92	} else {
93	d.c[k] = d.c[k - 1] + d.p[k];
94	}
95	}
96
97	// Continuous portion
98	for (int k = d.n_discrete; k < n; ++k) {
99	if (k == d.n_discrete) {
100	d.c[k] = d.c[k - 1] + d.p[k];
101	} else {
102	if (d.interpolation == Interpolation::histogram) {
103	d.c[k] = d.c[k - 1] + d.p[k - 1] * (d.e_out[k] - d.e_out[k - 1]);
104	} else if (d.interpolation == Interpolation::lin_lin) {
105	d.c[k] = d.c[k - 1] + 0.5 * (d.p[k - 1] + d.p[k]) *
106	(d.e_out[k] - d.e_out[k - 1]);
107	}
108	}
109	}
110
111	// Normalize density and distribution functions
112	d.p /= d.c[n - 1];
113	d.c /= d.c[n - 1];
114	}
115
116	distribution_.push_back(std::move(d));	1,159,064✔
117	} // incoming energies	1,159,064✔
118	}	59,993✔
119
120	void KalbachMann::sample(	3,201,388✔
121	double E_in, double& E_out, double& mu, uint64_t* seed) const
122	{
123	// Find energy bin and calculate interpolation factor -- if the energy is
124	// outside the range of the tabulated energies, choose the first or last bins
125	auto n_energy_in = energy_.size();	3,201,388✔
126	int i;
127	double r;
128	if (E_in < energy_[0]) {	3,201,388✔
129	i = 0;	×
130	r = 0.0;	×
131	} else if (E_in > energy_[n_energy_in - 1]) {	3,201,388✔
132	i = n_energy_in - 2;	×
133	r = 1.0;	×
134	} else {
135	i = lower_bound_index(energy_.begin(), energy_.end(), E_in);	3,201,388✔
136	r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);	3,201,388✔
137	}
138
139	// Sample between the ith and [i+1]th bin
140	int l = r > prn(seed) ? i + 1 : i;	3,201,388✔
141
142	// Interpolation for energy E1 and EK
143	int n_energy_out = distribution_[i].e_out.size();	3,201,388✔
144	int n_discrete = distribution_[i].n_discrete;	3,201,388✔
145	double E_i_1 = distribution_[i].e_out[n_discrete];	3,201,388✔
146	double E_i_K = distribution_[i].e_out[n_energy_out - 1];	3,201,388✔
147
148	n_energy_out = distribution_[i + 1].e_out.size();	3,201,388✔
149	n_discrete = distribution_[i + 1].n_discrete;	3,201,388✔
150	double E_i1_1 = distribution_[i + 1].e_out[n_discrete];	3,201,388✔
151	double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];	3,201,388✔
152
153	double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);	3,201,388✔
154	double E_K = E_i_K + r * (E_i1_K - E_i_K);	3,201,388✔
155
156	// Determine outgoing energy bin
157	n_energy_out = distribution_[l].e_out.size();	3,201,388✔
158	n_discrete = distribution_[l].n_discrete;	3,201,388✔
159	double r1 = prn(seed);	3,201,388✔
160	double c_k = distribution_[l].c[0];	3,201,388✔
161	int k = 0;	3,201,388✔
162	int end = n_energy_out - 2;	3,201,388✔
163
164	// Discrete portion
165	for (int j = 0; j < n_discrete; ++j) {	3,201,388✔
166	k = j;	×
167	c_k = distribution_[l].c[k];	×
168	if (r1 < c_k) {	×
169	end = j;	×
170	break;	×
171	}
172	}
173
174	// Continuous portion
175	double c_k1;
176	for (int j = n_discrete; j < end; ++j) {	87,777,381✔
177	k = j;	87,669,599✔
178	c_k1 = distribution_[l].c[k + 1];	87,669,599✔
179	if (r1 < c_k1)	87,669,599✔
180	break;	3,093,606✔
181	k = j + 1;	84,575,993✔
182	c_k = c_k1;	84,575,993✔
183	}
184
185	double E_l_k = distribution_[l].e_out[k];	3,201,388✔
186	double p_l_k = distribution_[l].p[k];	3,201,388✔
187	double km_r, km_a;
188	if (distribution_[l].interpolation == Interpolation::histogram) {	3,201,388✔
189	// Histogram interpolation
190	if (p_l_k > 0.0 && k >= n_discrete) {	3,181,060✔
191	E_out = E_l_k + (r1 - c_k) / p_l_k;	3,181,060✔
192	} else {
193	E_out = E_l_k;	×
194	}
195
196	// Determine Kalbach-Mann parameters
197	km_r = distribution_[l].r[k];	3,181,060✔
198	km_a = distribution_[l].a[k];	3,181,060✔
199
200	} else {
201	// Linear-linear interpolation
202	double E_l_k1 = distribution_[l].e_out[k + 1];	20,328✔
203	double p_l_k1 = distribution_[l].p[k + 1];	20,328✔
204
205	double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);	20,328✔
206	if (frac == 0.0) {	20,328✔
207	E_out = E_l_k + (r1 - c_k) / p_l_k;	×
208	} else {
209	E_out =	20,328✔
210	E_l_k +	20,328✔
211	(std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -	20,328✔
212	p_l_k) /	20,328✔
213	frac;
214	}
215
216	// Determine Kalbach-Mann parameters
217	km_r = distribution_[l].r[k] +	20,328✔
218	(E_out - E_l_k) / (E_l_k1 - E_l_k) *	40,656✔
219	(distribution_[l].r[k + 1] - distribution_[l].r[k]);	20,328✔
220	km_a = distribution_[l].a[k] +	20,328✔
221	(E_out - E_l_k) / (E_l_k1 - E_l_k) *	40,656✔
222	(distribution_[l].a[k + 1] - distribution_[l].a[k]);	20,328✔
223	}
224
225	// Now interpolate between incident energy bins i and i + 1
226	if (k >= n_discrete) {	3,201,388✔
227	if (l == i) {	3,201,388✔
228	E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);	1,599,190✔
229	} else {
230	E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);	1,602,198✔
231	}
232	}
233
234	// Sampled correlated angle from Kalbach-Mann parameters
235	if (prn(seed) > km_r) {	3,201,388✔
236	double T = uniform_distribution(-1., 1., seed) * std::sinh(km_a);	3,132,862✔
237	mu = std::log(T + std::sqrt(T * T + 1.0)) / km_a;	3,132,862✔
238	} else {
239	double r1 = prn(seed);	68,526✔
240	mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;	68,526✔
241	}
242	}	3,201,388✔
NEW 243	void KalbachMann::get_pdf(double det_pos[4], double E_in, double& E_out,	×
244	uint64_t* seed, Particle& p, std::vector<double>& mu_cm,
245	std::vector<double>& Js, std::vector<Particle>& ghost_particles,
246	std::vector<double>& pdfs_lab) const
247	{
248	// Find energy bin and calculate interpolation factor -- if the energy is
249	// outside the range of the tabulated energies, choose the first or last bins
NEW 250	const auto& nuc {data::nuclides[p.event_nuclide()]};	×
251
252	// double E_out;
NEW 253	auto n_energy_in = energy_.size();	×
254	int i;
255	double r;
NEW 256	if (E_in < energy_[0]) {	×
NEW 257	i = 0;	×
NEW 258	r = 0.0;	×
NEW 259	} else if (E_in > energy_[n_energy_in - 1]) {	×
NEW 260	i = n_energy_in - 2;	×
NEW 261	r = 1.0;	×
262	} else {
NEW 263	i = lower_bound_index(energy_.begin(), energy_.end(), E_in);	×
NEW 264	r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]);	×
265	}
266
267	// Sample between the ith and [i+1]th bin
NEW 268	int l = r > prn(seed) ? i + 1 : i;	×
269
270	// Interpolation for energy E1 and EK
NEW 271	int n_energy_out = distribution_[i].e_out.size();	×
NEW 272	int n_discrete = distribution_[i].n_discrete;	×
NEW 273	double E_i_1 = distribution_[i].e_out[n_discrete];	×
NEW 274	double E_i_K = distribution_[i].e_out[n_energy_out - 1];	×
275
NEW 276	n_energy_out = distribution_[i + 1].e_out.size();	×
NEW 277	n_discrete = distribution_[i + 1].n_discrete;	×
NEW 278	double E_i1_1 = distribution_[i + 1].e_out[n_discrete];	×
NEW 279	double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];	×
280
NEW 281	double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);	×
NEW 282	double E_K = E_i_K + r * (E_i1_K - E_i_K);	×
283
284	// Determine outgoing energy bin
NEW 285	n_energy_out = distribution_[l].e_out.size();	×
NEW 286	n_discrete = distribution_[l].n_discrete;	×
NEW 287	double r1 = prn(seed);	×
NEW 288	double c_k = distribution_[l].c[0];	×
NEW 289	int k = 0;	×
NEW 290	int end = n_energy_out - 2;	×
291
292	// Discrete portion
NEW 293	for (int j = 0; j < n_discrete; ++j) {	×
NEW 294	k = j;	×
NEW 295	c_k = distribution_[l].c[k];	×
NEW 296	if (r1 < c_k) {	×
NEW 297	end = j;	×
NEW 298	break;	×
299	}
300	}
301
302	// Continuous portion
303	double c_k1;
NEW 304	for (int j = n_discrete; j < end; ++j) {	×
NEW 305	k = j;	×
NEW 306	c_k1 = distribution_[l].c[k + 1];	×
NEW 307	if (r1 < c_k1)	×
NEW 308	break;	×
NEW 309	k = j + 1;	×
NEW 310	c_k = c_k1;	×
311	}
312
NEW 313	double E_l_k = distribution_[l].e_out[k];	×
NEW 314	double p_l_k = distribution_[l].p[k];	×
315	double km_r, km_a;
NEW 316	if (distribution_[l].interpolation == Interpolation::histogram) {	×
317	// Histogram interpolation
NEW 318	if (p_l_k > 0.0 && k >= n_discrete) {	×
NEW 319	E_out = E_l_k + (r1 - c_k) / p_l_k;	×
320	} else {
NEW 321	E_out = E_l_k;	×
322	}
323	// Determine Kalbach-Mann parameters
NEW 324	km_r = distribution_[l].r[k];	×
NEW 325	km_a = distribution_[l].a[k];	×
326
327	} else {
328	// Linear-linear interpolation
NEW 329	double E_l_k1 = distribution_[l].e_out[k + 1];	×
NEW 330	double p_l_k1 = distribution_[l].p[k + 1];	×
331
NEW 332	double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);	×
NEW 333	if (frac == 0.0) {	×
NEW 334	E_out = E_l_k + (r1 - c_k) / p_l_k;	×
335	} else {
NEW 336	E_out =	×
NEW 337	E_l_k +	×
NEW 338	(std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -	×
NEW 339	p_l_k) /	×
340	frac;
341	}
342	// Linear-linear interpolation
343	// Determine Kalbach-Mann parameters
NEW 344	km_r = distribution_[l].r[k] +	×
NEW 345	(E_out - E_l_k) / (E_l_k1 - E_l_k) *	×
NEW 346	(distribution_[l].r[k + 1] - distribution_[l].r[k]);	×
NEW 347	km_a = distribution_[l].a[k] +	×
NEW 348	(E_out - E_l_k) / (E_l_k1 - E_l_k) *	×
NEW 349	(distribution_[l].a[k + 1] - distribution_[l].a[k]);	×
350	}
351
352	// Now interpolate between incident energy bins i and i + 1
NEW 353	if (k >= n_discrete) {	×
NEW 354	if (l == i) {	×
NEW 355	E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);	×
356	} else {
NEW 357	E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);	×
358	}
359	}
360	/* function to transform pdf CM to lab frame - Future implementation
361	get_pdf_to_point_elastic(det_pos, p, mu_cm, Js, ghost_particles, E_out /
362	1e6); for (std::size_t i = 0; i < mu_cm.size(); ++i) {
363	// Assuming Js.size() is the same as mu_cm.size()
364	double mu_c = mu_cm[i];
365	double derivative = Js[i];
366	double pdf_cm = km_a / (2 * std::sinh(km_a)) *
367	(std::cosh(km_a * mu_c) +
368	km_r * std::sinh(km_a * mu_c)); // center of mass
369	pdfs_lab.push_back(pdf_cm / std::abs(derivative));
370	}
371	*/
372
373	// https://docs.openmc.org/en/v0.8.0/methods/physics.html#equation-KM-pdf-angle
374	// double pdf_mu = km_a / (2 * std::sinh(km_a)) * (std::cosh(km_a * mymu) +
375	// km_r * std::sinh(km_a * mymu)); // center of mass return pdf_mu;
376
NEW 377	const auto& rx {nuc->reactions_[p.event_index_mt()]};	×
NEW 378	if (!rx->scatter_in_cm_) {	×
NEW 379	fatal_error("didn't implement lab");	×
380	}
381	}
382
383	} // namespace openmc

openmc-dev / openmc / 17271733860

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