23097215950

Committed 14 Mar 2026 09:58PM UTC coverage: 81.435% (-0.1%) from 81.58%

Build # 23097215950

Build Type

push

github

Committed by

web-flow

Commit Message

Implement angular PDF evaluation for angle-energy distributions (#3550)

Co-authored-by: Your Name <you@example.com>
Co-authored-by: GuySten <62616591+GuySten@users.noreply.github.com>
Co-authored-by: GuySten <guyste@post.bgu.ac.il>
Co-authored-by: Eliezer214 <110336440+Eliezer214@users.noreply.github.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Run Details

17596 of 25365 branches covered (69.37%)

Branch coverage included in aggregate %.

67 of 181 new or added lines in 12 files covered. (37.02%)

3 existing lines in 1 file now uncovered.

58033 of 67505 relevant lines covered (85.97%)

44708275.34 hits per line

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

59.55

/src/distribution_multi.cpp

#include "openmc/distribution_multi.h"

#include <algorithm> // for move, clamp
#include <cmath>     // for sqrt, sin, cos, max

#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_dist.h"
#include "openmc/random_lcg.h"
#include "openmc/xml_interface.h"

namespace openmc {

unique_ptr<UnitSphereDistribution> UnitSphereDistribution::create(
  pugi::xml_node node)
{
  // Check for type of angular distribution
  std::string type;
  if (check_for_node(node, "type"))
    type = get_node_value(node, "type", true, true);
  if (type == "isotropic") {
    return UPtrAngle {new Isotropic(node)};
  } else if (type == "monodirectional") {
    return UPtrAngle {new Monodirectional(node)};
  } else if (type == "mu-phi") {
    return UPtrAngle {new PolarAzimuthal(node)};
  } else {
    fatal_error(fmt::format(
      "Invalid angular distribution for external source: {}", type));
  }
}

//==============================================================================
// UnitSphereDistribution implementation
//==============================================================================

UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node)
{
  // Read reference directional unit vector
  if (check_for_node(node, "reference_uvw")) {
    auto u_ref = get_node_array<double>(node, "reference_uvw");
    if (u_ref.size() != 3)
      fatal_error("Angular distribution reference direction must have "
                  "three parameters specified.");
    u_ref_ = Direction(u_ref.data());
    u_ref_ /= u_ref_.norm();
  }
}

//==============================================================================
// PolarAzimuthal implementation
//==============================================================================

PolarAzimuthal::PolarAzimuthal(Direction u, UPtrDist mu, UPtrDist phi)
  : UnitSphereDistribution {u}, mu_ {std::move(mu)}, phi_ {std::move(phi)}
{}

PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
  : UnitSphereDistribution {node}
{
  // Read reference directional unit vector
  if (check_for_node(node, "reference_vwu")) {
    auto v_ref = get_node_array<double>(node, "reference_vwu");
    if (v_ref.size() != 3)
      fatal_error("Angular distribution reference v direction must have "
                  "three parameters specified.");
    v_ref_ = Direction(v_ref.data());
    v_ref_ /= v_ref_.norm();
  }
  w_ref_ = u_ref_.cross(v_ref_);
  if (check_for_node(node, "mu")) {
    pugi::xml_node node_dist = node.child("mu");
    mu_ = distribution_from_xml(node_dist);
  } else {
    mu_ = UPtrDist {new Uniform(-1., 1.)};
  }

  if (check_for_node(node, "phi")) {
    pugi::xml_node node_dist = node.child("phi");
    phi_ = distribution_from_xml(node_dist);
  } else {
    phi_ = UPtrDist {new Uniform(0.0, 2.0 * PI)};
  }
}

std::pair<Direction, double> PolarAzimuthal::sample(uint64_t* seed) const
{
  return sample_impl(seed, false);
}

std::pair<Direction, double> PolarAzimuthal::sample_as_bias(
  uint64_t* seed) const
{
  return sample_impl(seed, true);
}

std::pair<Direction, double> PolarAzimuthal::sample_impl(
  uint64_t* seed, bool return_pdf) const
{
  // Sample cosine of polar angle
  auto [mu, mu_wgt] = mu_->sample(seed);

  // Sample azimuthal angle
  auto [phi, phi_wgt] = phi_->sample(seed);

  // Compute either the PDF value or the importance weight
  double weight =
    return_pdf ? (mu_->evaluate(mu) * phi_->evaluate(phi)) : (mu_wgt * phi_wgt);

  if (mu == 1.0)
    return {u_ref_, weight};
  if (mu == -1.0)
    return {-u_ref_, weight};

  double f = std::sqrt(1 - mu * mu);
  return {mu * u_ref_ + f * std::cos(phi) * v_ref_ + f * std::sin(phi) * w_ref_,
    weight};
}

double PolarAzimuthal::evaluate(Direction u) const
{
  double mu = std::clamp(u.dot(u_ref_), -1.0, 1.0);
  double phi = 0.0;
  double sin_theta_sq = std::max(0.0, 1.0 - mu * mu);
  if (sin_theta_sq > 0.0) {
    double sin_theta = std::sqrt(sin_theta_sq);
    double cos_phi = u.dot(v_ref_) / sin_theta;
    double sin_phi = u.dot(w_ref_) / sin_theta;
    phi = std::atan2(sin_phi, cos_phi);
    if (phi < 0.0)
      phi += 2.0 * PI;
  }
  return mu_->evaluate(mu) * phi_->evaluate(phi);
}

//==============================================================================
// Isotropic implementation
//==============================================================================

Isotropic::Isotropic(pugi::xml_node node) : UnitSphereDistribution {node}
{
  if (check_for_node(node, "bias")) {
    pugi::xml_node bias_node = node.child("bias");
    std::string bias_type = get_node_value(bias_node, "type", true, true);
    if (bias_type != "mu-phi") {
      openmc::fatal_error(
        "Isotropic distributions may only be biased by a PolarAzimuthal.");
    }
    auto bias = std::make_unique<PolarAzimuthal>(bias_node);
    if (bias->mu()->bias() || bias->phi()->bias()) {
      openmc::fatal_error(
        "Attempted to bias Isotropic distribution with a biased PolarAzimuthal "
        "distribution. Please ensure bias distributions are unbiased.");
    }
    this->set_bias(std::move(bias));
  }
}

Direction isotropic_direction(uint64_t* seed)
{
  double phi = uniform_distribution(0., 2.0 * PI, seed);
  double mu = uniform_distribution(-1., 1., seed);
  return {mu, std::sqrt(1.0 - mu * mu) * std::cos(phi),
    std::sqrt(1.0 - mu * mu) * std::sin(phi)};
}

std::pair<Direction, double> Isotropic::sample(uint64_t* seed) const
{
  if (bias()) {
    auto [val, eval] = bias()->sample_as_bias(seed);
    return {val, 1.0 / (4.0 * PI * eval)};
  } else {
    return {isotropic_direction(seed), 1.0};
  }
}

double Isotropic::evaluate(Direction u) const
{
  return 1.0 / (4.0 * PI);
}

//==============================================================================
// Monodirectional implementation
//==============================================================================

std::pair<Direction, double> Monodirectional::sample(uint64_t* seed) const
{
  return {u_ref_, 1.0};
}

} // namespace openmc

1	#include "openmc/distribution_multi.h"
2
3	#include <algorithm> // for move, clamp
4	#include <cmath> // for sqrt, sin, cos, max
5
6	#include "openmc/constants.h"
7	#include "openmc/error.h"
8	#include "openmc/math_functions.h"
9	#include "openmc/random_dist.h"
10	#include "openmc/random_lcg.h"
11	#include "openmc/xml_interface.h"
12
13	namespace openmc {
14
15	unique_ptr<UnitSphereDistribution> UnitSphereDistribution::create(	3,402✔
16	pugi::xml_node node)
17	{
18	// Check for type of angular distribution
19	std::string type;	3,402✔
20	if (check_for_node(node, "type"))	3,402!
21	type = get_node_value(node, "type", true, true);	3,402✔
22	if (type == "isotropic") {	3,402✔
23	return UPtrAngle {new Isotropic(node)};	249✔
24	} else if (type == "monodirectional") {	3,153✔
25	return UPtrAngle {new Monodirectional(node)};	3,116✔
26	} else if (type == "mu-phi") {	37!
27	return UPtrAngle {new PolarAzimuthal(node)};	37✔
28	} else {
29	fatal_error(fmt::format(	×
30	"Invalid angular distribution for external source: {}", type));
31	}
32	}	3,402✔
33
34	//==============================================================================
35	// UnitSphereDistribution implementation
36	//==============================================================================
37
38	UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node)	3,402✔
39	{
40	// Read reference directional unit vector
41	if (check_for_node(node, "reference_uvw")) {	3,402✔
42	auto u_ref = get_node_array<double>(node, "reference_uvw");	3,153✔
43	if (u_ref.size() != 3)	3,153!
44	fatal_error("Angular distribution reference direction must have "	×
45	"three parameters specified.");
46	u_ref_ = Direction(u_ref.data());	3,153✔
47	u_ref_ /= u_ref_.norm();	3,153✔
48	}	3,153✔
49	}	3,402✔
50
51	//==============================================================================
52	// PolarAzimuthal implementation
53	//==============================================================================
54
55	PolarAzimuthal::PolarAzimuthal(Direction u, UPtrDist mu, UPtrDist phi)	×
56	: UnitSphereDistribution {u}, mu_ {std::move(mu)}, phi_ {std::move(phi)}	×
57	{}	×
58
59	PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)	37✔
60	: UnitSphereDistribution {node}	37✔
61	{
62	// Read reference directional unit vector
63	if (check_for_node(node, "reference_vwu")) {	37!
64	auto v_ref = get_node_array<double>(node, "reference_vwu");	37✔
65	if (v_ref.size() != 3)	37!
66	fatal_error("Angular distribution reference v direction must have "	×
67	"three parameters specified.");
68	v_ref_ = Direction(v_ref.data());	37✔
69	v_ref_ /= v_ref_.norm();	37✔
70	}	37✔
71	w_ref_ = u_ref_.cross(v_ref_);	37✔
72	if (check_for_node(node, "mu")) {	37!
73	pugi::xml_node node_dist = node.child("mu");	37✔
74	mu_ = distribution_from_xml(node_dist);	37✔
75	} else {
76	mu_ = UPtrDist {new Uniform(-1., 1.)};	×
77	}
78
79	if (check_for_node(node, "phi")) {	37!
80	pugi::xml_node node_dist = node.child("phi");	37✔
81	phi_ = distribution_from_xml(node_dist);	37✔
82	} else {
83	phi_ = UPtrDist {new Uniform(0.0, 2.0 * PI)};	×
84	}
85	}	37✔
86
87	std::pair<Direction, double> PolarAzimuthal::sample(uint64_t* seed) const	25,300✔
88	{
89	return sample_impl(seed, false);	25,300✔
90	}
91
92	std::pair<Direction, double> PolarAzimuthal::sample_as_bias(	×
93	uint64_t* seed) const
94	{
95	return sample_impl(seed, true);	×
96	}
97
98	std::pair<Direction, double> PolarAzimuthal::sample_impl(	25,300✔
99	uint64_t* seed, bool return_pdf) const
100	{
101	// Sample cosine of polar angle
102	auto [mu, mu_wgt] = mu_->sample(seed);	25,300✔
103
104	// Sample azimuthal angle
105	auto [phi, phi_wgt] = phi_->sample(seed);	25,300!
106
107	// Compute either the PDF value or the importance weight
108	double weight =	25,300✔
109	return_pdf ? (mu_->evaluate(mu) * phi_->evaluate(phi)) : (mu_wgt * phi_wgt);	25,300!
110
111	if (mu == 1.0)	25,300✔
112	return {u_ref_, weight};	759✔
113	if (mu == -1.0)	24,541✔
114	return {-u_ref_, weight};	1,628✔
115
116	double f = std::sqrt(1 - mu * mu);	22,913✔
117	return {mu * u_ref_ + f * std::cos(phi) * v_ref_ + f * std::sin(phi) * w_ref_,	22,913✔
118	weight};	22,913✔
119	}
120
NEW 121	double PolarAzimuthal::evaluate(Direction u) const	×
122	{
NEW 123	double mu = std::clamp(u.dot(u_ref_), -1.0, 1.0);	×
NEW 124	double phi = 0.0;	×
NEW 125	double sin_theta_sq = std::max(0.0, 1.0 - mu * mu);	×
NEW 126	if (sin_theta_sq > 0.0) {	×
NEW 127	double sin_theta = std::sqrt(sin_theta_sq);	×
NEW 128	double cos_phi = u.dot(v_ref_) / sin_theta;	×
NEW 129	double sin_phi = u.dot(w_ref_) / sin_theta;	×
NEW 130	phi = std::atan2(sin_phi, cos_phi);	×
NEW 131	if (phi < 0.0)	×
NEW 132	phi += 2.0 * PI;	×
133	}
NEW 134	return mu_->evaluate(mu) * phi_->evaluate(phi);	×
135	}
136
137	//==============================================================================
138	// Isotropic implementation
139	//==============================================================================
140
141	Isotropic::Isotropic(pugi::xml_node node) : UnitSphereDistribution {node}	249✔
142	{
143	if (check_for_node(node, "bias")) {	249!
144	pugi::xml_node bias_node = node.child("bias");	×
145	std::string bias_type = get_node_value(bias_node, "type", true, true);	×
146	if (bias_type != "mu-phi") {	×
147	openmc::fatal_error(	×
148	"Isotropic distributions may only be biased by a PolarAzimuthal.");
149	}
150	auto bias = std::make_unique<PolarAzimuthal>(bias_node);	×
151	if (bias->mu()->bias() \|\| bias->phi()->bias()) {	×
152	openmc::fatal_error(	×
153	"Attempted to bias Isotropic distribution with a biased PolarAzimuthal "
154	"distribution. Please ensure bias distributions are unbiased.");
155	}
156	this->set_bias(std::move(bias));	×
157	}	×
158	}	249✔
159
160	Direction isotropic_direction(uint64_t* seed)	114,849,917✔
161	{
162	double phi = uniform_distribution(0., 2.0 * PI, seed);	114,849,917✔
163	double mu = uniform_distribution(-1., 1., seed);	114,849,917✔
164	return {mu, std::sqrt(1.0 - mu * mu) * std::cos(phi),	114,849,917✔
165	std::sqrt(1.0 - mu * mu) * std::sin(phi)};	114,849,917✔
166	}
167
168	std::pair<Direction, double> Isotropic::sample(uint64_t* seed) const	21,574,774✔
169	{
170	if (bias()) {	21,574,774!
171	auto [val, eval] = bias()->sample_as_bias(seed);	×
172	return {val, 1.0 / (4.0 * PI * eval)};	×
173	} else {
174	return {isotropic_direction(seed), 1.0};	21,574,774✔
175	}
176	}
177
NEW 178	double Isotropic::evaluate(Direction u) const	×
179	{
NEW 180	return 1.0 / (4.0 * PI);	×
181	}
182
183	//==============================================================================
184	// Monodirectional implementation
185	//==============================================================================
186
187	std::pair<Direction, double> Monodirectional::sample(uint64_t* seed) const	14,343,114✔
188	{
189	return {u_ref_, 1.0};	14,343,114✔
190	}
191
192	} // namespace openmc

openmc-dev / openmc / 23097215950

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