openmc/tests/cpp_unit_tests/test_distribution.cpp at develop · openmc-dev/openmc · GitHub

History

219 lines (182 loc) · 6.72 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

#include "openmc/distribution.h"

#include "openmc/distribution_spatial.h"

#include "openmc/position.h"

#include "openmc/random_lcg.h"

#include <catch2/catch_test_macros.hpp>

#include <catch2/matchers/catch_matchers_floating_point.hpp>

#include <cmath>

#include <pugixml.hpp>

TEST_CASE("Test alias method sampling of a discrete distribution")

{

constexpr int n_samples = 1000000;

double x[5] = {-1.6, 1.1, 20.3, 4.7, 0.9};

double p[5] = {0.2, 0.1, 0.65, 0.02, 0.03};

// Initialize distribution

openmc::Discrete dist(x, p, 5);

uint64_t seed = openmc::init_seed(0, 0);

// Calculate expected distribution mean

double mean = 0.0;

for (size_t i = 0; i < 5; i++) {

mean += x[i] * p[i];

}

// Sample distribution and calculate mean, standard deviation, and number of

// x[0] sampled

double dist_mean = 0.0;

double std = 0.0;

int counter = 0;

for (size_t i = 0; i < n_samples; i++) {

auto sample = dist.sample(&seed).first;

std += sample * sample / n_samples;

dist_mean += sample;

if (sample == x[0])

counter++;

}

dist_mean /= n_samples;

std -= dist_mean * dist_mean;

// Require sampled distribution mean is within 4 standard deviations of the

// expected mean

REQUIRE(std::abs(dist_mean - mean) < 4 * std);

// Require counter of number of x[0] is within the 95% confidence interval

// assuming a Poisson distribution of 200,000

REQUIRE(std::abs((double)counter / n_samples - p[0]) <

1.96 * std::sqrt(p[0] / n_samples));

}

TEST_CASE("Test alias sampling method for pugixml constructor")

{

// XML doc node for Discrete contructor

pugi::xml_document doc;

pugi::xml_node energy = doc.append_child("energy");

pugi::xml_node parameters = energy.append_child("parameters");

parameters.append_child(pugi::node_pcdata)

.set_value("800 500000 30000 0.1 0.6 0.3");

// Initialize discrete distribution and seed

openmc::Discrete dist(energy);

uint64_t seed = openmc::init_seed(0, 0);

auto sample = dist.sample(&seed).first;

// Assertions

REQUIRE(dist.x().size() == 3);

REQUIRE(dist.prob().size() == 3);

REQUIRE(dist.alias().size() == 3);

openmc::vector<double> correct_x = {800, 500000, 30000};

openmc::vector<double> correct_prob = {0.3, 1.0, 0.9};

openmc::vector<size_t> correct_alias = {1, 0, 1};

for (size_t i = 0; i < 3; i++) {

REQUIRE(dist.x()[i] == correct_x[i]);

REQUIRE_THAT(

dist.prob()[i], Catch::Matchers::WithinAbs(correct_prob[i], 1e-12));

REQUIRE(dist.alias()[i] == correct_alias[i]);

}

TEST_CASE("Test sampling a large linear-linear tabular distribution")

{

constexpr int n_points = 10001;

constexpr int n_samples = 200000;

openmc::vector<double> x(n_points);

openmc::vector<double> p(n_points);

for (int i = 0; i < n_points; ++i) {

x[i] = static_cast<double>(i) / (n_points - 1);

p[i] = 2.0 * x[i];

}

openmc::Tabular dist(

x.data(), p.data(), n_points, openmc::Interpolation::lin_lin);

uint64_t seed = openmc::init_seed(0, 0);

double mean = 0.0;

for (int i = 0; i < n_samples; ++i) {

mean += dist.sample(&seed).first;

}

mean /= n_samples;

// The normalized PDF is 2x on [0, 1], which has a mean of 2/3.

REQUIRE_THAT(mean, Catch::Matchers::WithinAbs(2.0 / 3.0, 0.003));

}

TEST_CASE("Test construction of SpatialBox with parameters")

{

openmc::Position ll {-1, -2, -3};

openmc::Position ur {30, 15, 5};

openmc::SpatialBox box(ll, ur);

REQUIRE(box.lower_left() == openmc::Position {-1, -2, -3});

REQUIRE(box.upper_right() == openmc::Position {30, 15, 5});

REQUIRE_FALSE(box.only_fissionable());

}

TEST_CASE("Test Normal distribution")

{

// Test untruncated normal distribution

openmc::Normal normal_unbounded(0.0, 1.0);

// Check PDF at mean (should be 1/sqrt(2*pi) ≈ 0.3989)

REQUIRE_THAT(

normal_unbounded.evaluate(0.0), Catch::Matchers::WithinRel(0.3989, 0.001));

// Check that it's not truncated

REQUIRE_FALSE(normal_unbounded.is_truncated());

// Check accessors

REQUIRE(normal_unbounded.mean_value() == 0.0);

REQUIRE(normal_unbounded.std_dev() == 1.0);

REQUIRE(normal_unbounded.lower() == -openmc::INFTY);

REQUIRE(normal_unbounded.upper() == openmc::INFTY);

}

TEST_CASE("Test truncated Normal distribution")

{

// Create a truncated normal: mean=0, std=1, bounds=[-1, 1]

openmc::Normal normal_truncated(0.0, 1.0, -1.0, 1.0);

// Check that it's truncated

REQUIRE(normal_truncated.is_truncated());

// Check accessors

REQUIRE(normal_truncated.lower() == -1.0);

REQUIRE(normal_truncated.upper() == 1.0);

// PDF should be zero outside bounds

REQUIRE(normal_truncated.evaluate(-2.0) == 0.0);

REQUIRE(normal_truncated.evaluate(2.0) == 0.0);

// PDF inside bounds should be higher than untruncated (due to

// renormalization)

openmc::Normal normal_unbounded(0.0, 1.0);

REQUIRE(normal_truncated.evaluate(0.0) > normal_unbounded.evaluate(0.0));

// The truncated PDF at mean should be approximately 0.3989 / 0.6827 ≈ 0.584

// (0.6827 is the probability mass of N(0,1) in [-1,1])

REQUIRE_THAT(

normal_truncated.evaluate(0.0), Catch::Matchers::WithinRel(0.584, 0.01));

}

TEST_CASE("Test truncated Normal sampling")

{

constexpr int n_samples = 10000;

openmc::Normal normal_truncated(0.0, 1.0, -1.0, 1.0);

uint64_t seed = openmc::init_seed(0, 0);

// Sample and verify all samples are within bounds

for (int i = 0; i < n_samples; ++i) {

auto [x, w] = normal_truncated.sample(&seed);

REQUIRE(x >= -1.0);

REQUIRE(x <= 1.0);

REQUIRE(w == 1.0); // Unbiased sampling should have weight 1

}

TEST_CASE("Test one-sided truncated Normal")

{

// Test lower-bounded only (positive half-normal)

openmc::Normal lower_bounded(0.0, 1.0, 0.0, openmc::INFTY);

REQUIRE(lower_bounded.is_truncated());

REQUIRE(lower_bounded.evaluate(-1.0) == 0.0);

REQUIRE(lower_bounded.evaluate(1.0) > 0.0);

// PDF at 0 should be approximately 2 * 0.3989 ≈ 0.798 (half-normal)

REQUIRE_THAT(

lower_bounded.evaluate(0.0), Catch::Matchers::WithinRel(0.798, 0.01));

// Test upper-bounded only

openmc::Normal upper_bounded(0.0, 1.0, -openmc::INFTY, 0.0);

REQUIRE(upper_bounded.is_truncated());

REQUIRE(upper_bounded.evaluate(1.0) == 0.0);

REQUIRE(upper_bounded.evaluate(-1.0) > 0.0);

}

TEST_CASE("Test Normal XML constructor with truncation")

{

// XML doc node for truncated Normal

pugi::xml_document doc;

pugi::xml_node energy = doc.append_child("energy");

energy.append_child("type")

.append_child(pugi::node_pcdata)

.set_value("normal");

energy.append_child("parameters")

.append_child(pugi::node_pcdata)

.set_value("1.0e6 1.0e5 0.8e6 1.2e6");

openmc::Normal dist(energy);

REQUIRE(dist.mean_value() == 1.0e6);

REQUIRE(dist.std_dev() == 1.0e5);

REQUIRE(dist.lower() == 0.8e6);

REQUIRE(dist.upper() == 1.2e6);

REQUIRE(dist.is_truncated());

}

openmc/tests/cpp_unit_tests/test_distribution.cpp at develop · openmc-dev/openmc · GitHub

Uh oh!

Expand file tree

Breadcrumbs

test_distribution.cpp

Latest commit

History

Breadcrumbs

test_distribution.cpp

File metadata and controls