24471518122

Committed 15 Apr 2026 06:33PM UTC coverage: 81.192% (-0.1%) from 81.34%

Build # 24471518122

Build Type

Pull #3923

github

Committed by

web-flow

Commit Message

Merge 9fc69e91b into fd1bc26af

Pull Request Pull Request #3923: Increase numerical robustness of cylindrical and spherical distance to boundary checks

Coverage Stats

17131 of 24547 branches covered (69.79%)

Branch coverage included in aggregate %.

12 of 12 new or added lines in 1 file covered. (100.0%)

301 existing lines in 23 files now uncovered.

57289 of 67112 relevant lines covered (85.36%)

41016206.68 hits per line

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

97.96

/src/random_lcg.cpp

#include "openmc/random_lcg.h"

#include <cmath>

namespace openmc {

// Starting seed
int64_t master_seed {1};

// LCG parameters
constexpr uint64_t prn_mult {6364136223846793005ULL}; // multiplication
constexpr uint64_t prn_add {1442695040888963407ULL};  // additive factor, c
uint64_t prn_stride {DEFAULT_STRIDE}; // stride between particles

//==============================================================================
// PRN
//==============================================================================

// 64 bit implementation of the PCG-RXS-M-XS 64-bit state / 64-bit output
// geneator Adapted from: https://github.com/imneme/pcg-c, in particular
// https://github.com/imneme/pcg-c/blob/83252d9c23df9c82ecb42210afed61a7b42402d7/include/pcg_variants.h#L188-L192
// @techreport{oneill:pcg2014,
//    title = "PCG: A Family of Simple Fast Space-Efficient Statistically Good
//    Algorithms for Random Number Generation", author = "Melissa E. O'Neill",
//    institution = "Harvey Mudd College",
//    address = "Claremont, CA",
//    number = "HMC-CS-2014-0905",
//    year = "2014",
//    month = Sep,
//    xurl = "https://www.cs.hmc.edu/tr/hmc-cs-2014-0905.pdf",
//}
double prn(uint64_t* seed)
{
  // Advance the LCG
  *seed = (prn_mult * (*seed) + prn_add);

  // Permute the output
  uint64_t word =
    ((*seed >> ((*seed >> 59u) + 5u)) ^ *seed) * 12605985483714917081ull;
  uint64_t result = (word >> 43u) ^ word;

  // Convert output from unsigned integer to double
  return ldexp(result, -64);
}

//==============================================================================
// FUTURE_PRN
//==============================================================================

double future_prn(int64_t n, uint64_t seed)
{
  uint64_t fseed = future_seed(static_cast<uint64_t>(n), seed);
  return prn(&fseed);
}

//==============================================================================
// INIT_SEED
//==============================================================================

uint64_t init_seed(int64_t id, int offset)
{
  return future_seed(
    static_cast<uint64_t>(id) * prn_stride, master_seed + offset);
}

//==============================================================================
// INIT_PARTICLE_SEEDS
//==============================================================================

void init_particle_seeds(int64_t id, uint64_t* seeds)
{
  for (int i = 0; i < N_STREAMS; i++) {
    seeds[i] =
      future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + i);
  }
}

//==============================================================================
// ADVANCE_PRN_SEED
//==============================================================================

void advance_prn_seed(int64_t n, uint64_t* seed)
{
  *seed = future_seed(static_cast<uint64_t>(n), *seed);
}

//==============================================================================
// FUTURE_SEED
//==============================================================================

uint64_t future_seed(uint64_t n, uint64_t seed)
{
  // The algorithm here to determine the parameters used to skip ahead is
  // described in F. Brown, "Random Number Generation with Arbitrary Stride,"
  // Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in
  // O(log2(N)) operations instead of O(N). Basically, it computes parameters G
  // and C which can then be used to find x_N = G*x_0 + C mod 2^M.

  // Initialize constants
  uint64_t g {prn_mult};
  uint64_t c {prn_add};
  uint64_t g_new {1};
  uint64_t c_new {0};

  while (n > 0) {
    // Check if the least significant bit is 1.
    if (n & 1) {
      g_new *= g;
      c_new = c_new * g + c;
    }
    c *= (g + 1);
    g *= g;

    // Move bits right, dropping least significant bit.
    n >>= 1;
  }

  // With G and C, we can now find the new seed.
  return g_new * seed + c_new;
}

//==============================================================================
//                               API FUNCTIONS
//==============================================================================

extern "C" int64_t openmc_get_seed()
{
  return master_seed;
}

extern "C" void openmc_set_seed(int64_t new_seed)
{
  master_seed = new_seed;
}

extern "C" uint64_t openmc_get_stride()
{
  return prn_stride;
}

extern "C" void openmc_set_stride(uint64_t new_stride)
{
  prn_stride = new_stride;
}

} // namespace openmc

1	#include "openmc/random_lcg.h"
2
3	#include <cmath>
4
5	namespace openmc {
6
7	// Starting seed
8	int64_t master_seed {1};
9
10	// LCG parameters
11	constexpr uint64_t prn_mult {6364136223846793005ULL}; // multiplication
12	constexpr uint64_t prn_add {1442695040888963407ULL}; // additive factor, c
13	uint64_t prn_stride {DEFAULT_STRIDE}; // stride between particles
14
15	//==============================================================================
16	// PRN
17	//==============================================================================
18
19	// 64 bit implementation of the PCG-RXS-M-XS 64-bit state / 64-bit output
20	// geneator Adapted from: https://github.com/imneme/pcg-c, in particular
21	// https://github.com/imneme/pcg-c/blob/83252d9c23df9c82ecb42210afed61a7b42402d7/include/pcg_variants.h#L188-L192
22	// @techreport{oneill:pcg2014,
23	// title = "PCG: A Family of Simple Fast Space-Efficient Statistically Good
24	// Algorithms for Random Number Generation", author = "Melissa E. O'Neill",
25	// institution = "Harvey Mudd College",
26	// address = "Claremont, CA",
27	// number = "HMC-CS-2014-0905",
28	// year = "2014",
29	// month = Sep,
30	// xurl = "https://www.cs.hmc.edu/tr/hmc-cs-2014-0905.pdf",
31	//}
32	double prn(uint64_t* seed)	2,147,483,647✔
33	{
34	// Advance the LCG
35	seed = (prn_mult (*seed) + prn_add);	2,147,483,647✔
36
37	// Permute the output
38	uint64_t word =	2,147,483,647✔
39	((seed >> ((seed >> 59u) + 5u)) ^ seed) 12605985483714917081ull;	2,147,483,647✔
40	uint64_t result = (word >> 43u) ^ word;	2,147,483,647✔
41
42	// Convert output from unsigned integer to double
43	return ldexp(result, -64);	2,147,483,647✔
44	}
45
46	//==============================================================================
47	// FUTURE_PRN
48	//==============================================================================
49
50	double future_prn(int64_t n, uint64_t seed)	213,437,329✔
51	{
52	uint64_t fseed = future_seed(static_cast<uint64_t>(n), seed);	213,437,329✔
53	return prn(&fseed);	213,437,329✔
54	}
55
56	//==============================================================================
57	// INIT_SEED
58	//==============================================================================
59
60	uint64_t init_seed(int64_t id, int offset)	105,583,056✔
61	{
62	return future_seed(	105,583,056✔
63	static_cast<uint64_t>(id) * prn_stride, master_seed + offset);	105,583,056✔
64	}
65
66	//==============================================================================
67	// INIT_PARTICLE_SEEDS
68	//==============================================================================
69
70	void init_particle_seeds(int64_t id, uint64_t* seeds)	211,954,738✔
71	{
72	for (int i = 0; i < N_STREAMS; i++) {	1,059,773,690✔
73	seeds[i] =	847,818,952✔
74	future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + i);	847,818,952✔
75	}
76	}	211,954,738✔
77
78	//==============================================================================
79	// ADVANCE_PRN_SEED
80	//==============================================================================
81
82	void advance_prn_seed(int64_t n, uint64_t* seed)	816,419,997✔
83	{
84	seed = future_seed(static_cast<uint64_t>(n), seed);	816,419,997✔
85	}	816,419,997✔
86
87	//==============================================================================
88	// FUTURE_SEED
89	//==============================================================================
90
91	uint64_t future_seed(uint64_t n, uint64_t seed)	1,983,259,334✔
92	{
93	// The algorithm here to determine the parameters used to skip ahead is
94	// described in F. Brown, "Random Number Generation with Arbitrary Stride,"
95	// Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in
96	// O(log2(N)) operations instead of O(N). Basically, it computes parameters G
97	// and C which can then be used to find x_N = G*x_0 + C mod 2^M.
98
99	// Initialize constants
100	uint64_t g {prn_mult};	1,983,259,334✔
101	uint64_t c {prn_add};	1,983,259,334✔
102	uint64_t g_new {1};	1,983,259,334✔
103	uint64_t c_new {0};	1,983,259,334✔
104
UNCOV 105	while (n > 0) {	✔
106	// Check if the least significant bit is 1.
107	if (n & 1) {	2,147,483,647✔
108	g_new *= g;	2,147,483,647✔
109	c_new = c_new * g + c;	2,147,483,647✔
110	}
111	c *= (g + 1);	2,147,483,647✔
112	g *= g;	2,147,483,647✔
113
114	// Move bits right, dropping least significant bit.
115	n >>= 1;	2,147,483,647✔
116	}
117
118	// With G and C, we can now find the new seed.
119	return g_new * seed + c_new;	1,983,259,334✔
120	}
121
122	//==============================================================================
123	// API FUNCTIONS
124	//==============================================================================
125
126	extern "C" int64_t openmc_get_seed()	5,776✔
127	{
128	return master_seed;	5,776✔
129	}
130
131	extern "C" void openmc_set_seed(int64_t new_seed)	14,034✔
132	{
133	master_seed = new_seed;	14,034✔
134	}	14,034✔
135
136	extern "C" uint64_t openmc_get_stride()	5,767✔
137	{
138	return prn_stride;	5,767✔
139	}
140
141	extern "C" void openmc_set_stride(uint64_t new_stride)	13,580✔
142	{
143	prn_stride = new_stride;	13,580✔
144	}	13,580✔
145
146	} // namespace openmc

openmc-dev / openmc / 24471518122

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