19041149218

Committed 03 Nov 2025 04:12PM UTC coverage: 81.294% (-0.6%) from 81.872%

Build # 19041149218

Build Type

Pull #3550

github

Committed by

web-flow

Commit Message

Merge cc7c3a5ec into fd964bc9b

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

Run Details

16664 of 23490 branches covered (70.94%)

Branch coverage included in aggregate %.

2 of 387 new or added lines in 13 files covered. (0.52%)

3 existing lines in 1 file now uncovered.

54148 of 63616 relevant lines covered (85.12%)

43073383.89 hits per line

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

29.87

/src/secondary_nbody.cpp

#include "openmc/secondary_nbody.h"

#include <cmath> // for log

#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/random_dist.h"
#include "openmc/random_lcg.h"

namespace openmc {

//==============================================================================
// NBodyPhaseSpace implementation
//==============================================================================

NBodyPhaseSpace::NBodyPhaseSpace(hid_t group)
{
  read_attribute(group, "n_particles", n_bodies_);
  read_attribute(group, "total_mass", mass_ratio_);
  read_attribute(group, "atomic_weight_ratio", A_);
  read_attribute(group, "q_value", Q_);
}

void NBodyPhaseSpace::sample(
  double E_in, double& E_out, double& mu, uint64_t* seed) const
{
  // By definition, the distribution of the angle is isotropic for an N-body
  // phase space distribution
  mu = uniform_distribution(-1., 1., seed);

  // Determine E_max parameter
  double Ap = mass_ratio_;
  double E_max = (Ap - 1.0) / Ap * (A_ / (A_ + 1.0) * E_in + Q_);

  // x is essentially a Maxwellian distribution
  double x = maxwell_spectrum(1.0, seed);

  double y;
  double r1, r2, r3, r4, r5, r6;
  switch (n_bodies_) {
  case 3:
    y = maxwell_spectrum(1.0, seed);
    break;
  case 4:
    r1 = prn(seed);
    r2 = prn(seed);
    r3 = prn(seed);
    y = -std::log(r1 * r2 * r3);
    break;
  case 5:
    r1 = prn(seed);
    r2 = prn(seed);
    r3 = prn(seed);
    r4 = prn(seed);
    r5 = prn(seed);
    r6 = prn(seed);
    y = -std::log(r1 * r2 * r3 * r4) -
        std::log(r5) * std::pow(std::cos(PI / 2.0 * r6), 2);
    break;
  default:
    fatal_error("N-body phase space with >5 bodies.");
  }

  // Now determine v and E_out
  double v = x / (x + y);
  E_out = E_max * v;
}

double NBodyPhaseSpace::sample_energy_and_pdf(
  double E_in, double mu, double& E_out, uint64_t* seed) const
{
  // Determine E_max parameter
  double Ap = mass_ratio_;
  double E_max = (Ap - 1.0) / Ap * (A_ / (A_ + 1.0) * E_in + Q_);

  // x is essentially a Maxwellian distribution
  double x = maxwell_spectrum(1.0, seed);

  double y;
  double r1, r2, r3, r4, r5, r6;
  switch (n_bodies_) {
  case 3:
    y = maxwell_spectrum(1.0, seed);
    break;
  case 4:
    r1 = prn(seed);
    r2 = prn(seed);
    r3 = prn(seed);
    y = -std::log(r1 * r2 * r3);
    break;
  case 5:
    r1 = prn(seed);
    r2 = prn(seed);
    r3 = prn(seed);
    r4 = prn(seed);
    r5 = prn(seed);
    r6 = prn(seed);
    y = -std::log(r1 * r2 * r3 * r4) -
        std::log(r5) * std::pow(std::cos(PI / 2.0 * r6), 2);
    break;
  default:
    fatal_error("N-body phase space with >5 bodies.");
  }

  // Now determine v and E_out
  double v = x / (x + y);
  E_out = E_max * v;
  return 0.5;
}

} // namespace openmc

1	#include "openmc/secondary_nbody.h"
2
3	#include <cmath> // for log
4
5	#include "openmc/constants.h"
6	#include "openmc/hdf5_interface.h"
7	#include "openmc/math_functions.h"
8	#include "openmc/random_dist.h"
9	#include "openmc/random_lcg.h"
10
11	namespace openmc {
12
13	//==============================================================================
14	// NBodyPhaseSpace implementation
15	//==============================================================================
16
17	NBodyPhaseSpace::NBodyPhaseSpace(hid_t group)	735✔
18	{
19	read_attribute(group, "n_particles", n_bodies_);	735✔
20	read_attribute(group, "total_mass", mass_ratio_);	735✔
21	read_attribute(group, "atomic_weight_ratio", A_);	735✔
22	read_attribute(group, "q_value", Q_);	735✔
23	}	735✔
24
25	void NBodyPhaseSpace::sample(	440✔
26	double E_in, double& E_out, double& mu, uint64_t* seed) const
27	{
28	// By definition, the distribution of the angle is isotropic for an N-body
29	// phase space distribution
30	mu = uniform_distribution(-1., 1., seed);	440✔
31
32	// Determine E_max parameter
33	double Ap = mass_ratio_;	440✔
34	double E_max = (Ap - 1.0) / Ap * (A_ / (A_ + 1.0) * E_in + Q_);	440✔
35
36	// x is essentially a Maxwellian distribution
37	double x = maxwell_spectrum(1.0, seed);	440✔
38
39	double y;
40	double r1, r2, r3, r4, r5, r6;
41	switch (n_bodies_) {	440!
42	case 3:	440✔
43	y = maxwell_spectrum(1.0, seed);	440✔
44	break;	440✔
45	case 4:	×
46	r1 = prn(seed);	×
47	r2 = prn(seed);	×
48	r3 = prn(seed);	×
49	y = -std::log(r1 * r2 * r3);	×
50	break;	×
51	case 5:	×
52	r1 = prn(seed);	×
53	r2 = prn(seed);	×
54	r3 = prn(seed);	×
55	r4 = prn(seed);	×
56	r5 = prn(seed);	×
57	r6 = prn(seed);	×
58	y = -std::log(r1 * r2 * r3 * r4) -	×
59	std::log(r5) * std::pow(std::cos(PI / 2.0 * r6), 2);	×
60	break;	×
61	default:	×
NEW 62	fatal_error("N-body phase space with >5 bodies.");	×
63	}
64
65	// Now determine v and E_out
66	double v = x / (x + y);	440✔
67	E_out = E_max * v;	440✔
68	}	440✔
69
NEW 70	double NBodyPhaseSpace::sample_energy_and_pdf(	×
71	double E_in, double mu, double& E_out, uint64_t* seed) const
72	{
73	// Determine E_max parameter
NEW 74	double Ap = mass_ratio_;	×
NEW 75	double E_max = (Ap - 1.0) / Ap * (A_ / (A_ + 1.0) * E_in + Q_);	×
76
77	// x is essentially a Maxwellian distribution
NEW 78	double x = maxwell_spectrum(1.0, seed);	×
79
80	double y;
81	double r1, r2, r3, r4, r5, r6;
NEW 82	switch (n_bodies_) {	×
NEW 83	case 3:	×
NEW 84	y = maxwell_spectrum(1.0, seed);	×
NEW 85	break;	×
NEW 86	case 4:	×
NEW 87	r1 = prn(seed);	×
NEW 88	r2 = prn(seed);	×
NEW 89	r3 = prn(seed);	×
NEW 90	y = -std::log(r1 * r2 * r3);	×
NEW 91	break;	×
NEW 92	case 5:	×
NEW 93	r1 = prn(seed);	×
NEW 94	r2 = prn(seed);	×
NEW 95	r3 = prn(seed);	×
NEW 96	r4 = prn(seed);	×
NEW 97	r5 = prn(seed);	×
NEW 98	r6 = prn(seed);	×
NEW 99	y = -std::log(r1 * r2 * r3 * r4) -	×
NEW 100	std::log(r5) * std::pow(std::cos(PI / 2.0 * r6), 2);	×
NEW 101	break;	×
NEW 102	default:	×
NEW 103	fatal_error("N-body phase space with >5 bodies.");	×
104	}
105
106	// Now determine v and E_out
NEW 107	double v = x / (x + y);	×
NEW 108	E_out = E_max * v;	×
NEW 109	return 0.5;	×
110	}
111
112	} // namespace openmc

openmc-dev / openmc / 19041149218

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