17445237309

Committed 03 Sep 2025 08:25PM UTC coverage: 84.864% (-0.3%) from 85.209%

Build # 17445237309

Build Type

Pull #3460

github

Committed by

web-flow

Commit Message

Merge 5e7a2eae9 into 591856472

Pull Request Pull Request #3460: Source biasing capabilities

Run Details

343 of 638 new or added lines in 11 files covered. (53.76%)

4 existing lines in 3 files now uncovered.

53134 of 62611 relevant lines covered (84.86%)

38459716.64 hits per line

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

65.71

/src/distribution_multi.cpp

#include "openmc/distribution_multi.h"

#include <algorithm> // for move
#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());
  }
}

//==============================================================================
// 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}
{
  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
{
  // Sample cosine of polar angle
  auto [mu, mu_wgt] = mu_->sample(seed);

  // Sample azimuthal angle
  auto [phi, phi_wgt] = phi_->sample(seed);
  if (mu == 1.0)
    return {u_ref_, mu_wgt * phi_wgt};

  return {rotate_angle(u_ref_, mu, &phi, seed), mu_wgt * phi_wgt};
}

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

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

  double pdf_evaluation = mu_->evaluate(mu) * phi_->evaluate(phi);

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

  return {rotate_angle(u_ref_, mu, &phi, seed), pdf_evaluation};
}

//==============================================================================
// 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};
  }
}

//==============================================================================
// 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
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,169✔
16	pugi::xml_node node)
17	{
18	// Check for type of angular distribution
19	std::string type;	3,169✔
20	if (check_for_node(node, "type"))	3,169✔
21	type = get_node_value(node, "type", true, true);	3,169✔
22	if (type == "isotropic") {	3,169✔
23	return UPtrAngle {new Isotropic(node)};	255✔
24	} else if (type == "monodirectional") {	2,914✔
25	return UPtrAngle {new Monodirectional(node)};	2,876✔
26	} else if (type == "mu-phi") {	38✔
27	return UPtrAngle {new PolarAzimuthal(node)};	38✔
28	} else {
29	fatal_error(fmt::format(	×
30	"Invalid angular distribution for external source: {}", type));
31	}
32	}	3,169✔
33
34	//==============================================================================
35	// UnitSphereDistribution implementation
36	//==============================================================================
37
38	UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node)	3,169✔
39	{
40	// Read reference directional unit vector
41	if (check_for_node(node, "reference_uvw")) {	3,169✔
42	auto u_ref = get_node_array<double>(node, "reference_uvw");	2,914✔
43	if (u_ref.size() != 3)	2,914✔
44	fatal_error("Angular distribution reference direction must have "	×
45	"three parameters specified.");
46	u_ref_ = Direction(u_ref.data());	2,914✔
47	}	2,914✔
48	}	3,169✔
49
50	//==============================================================================
51	// PolarAzimuthal implementation
52	//==============================================================================
53
54	PolarAzimuthal::PolarAzimuthal(Direction u, UPtrDist mu, UPtrDist phi)	×
55	: UnitSphereDistribution {u}, mu_ {std::move(mu)}, phi_ {std::move(phi)}	×
56	{}	×
57
58	PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)	38✔
59	: UnitSphereDistribution {node}	38✔
60	{
61	if (check_for_node(node, "mu")) {	38✔
62	pugi::xml_node node_dist = node.child("mu");	38✔
63	mu_ = distribution_from_xml(node_dist);	38✔
64	} else {
65	mu_ = UPtrDist {new Uniform(-1., 1.)};	×
66	}
67
68	if (check_for_node(node, "phi")) {	38✔
69	pugi::xml_node node_dist = node.child("phi");	38✔
70	phi_ = distribution_from_xml(node_dist);	38✔
71	} else {
72	phi_ = UPtrDist {new Uniform(0.0, 2.0 * PI)};	×
73	}
74	}	38✔
75
76	std::pair<Direction, double> PolarAzimuthal::sample(uint64_t* seed) const	25,300✔
77	{
78	// Sample cosine of polar angle
79	auto [mu, mu_wgt] = mu_->sample(seed);	25,300✔
80
81	// Sample azimuthal angle
82	auto [phi, phi_wgt] = phi_->sample(seed);	25,300✔
83	if (mu == 1.0)	25,300✔
84	return {u_ref_, mu_wgt * phi_wgt};	759✔
85
86	return {rotate_angle(u_ref_, mu, &phi, seed), mu_wgt * phi_wgt};	24,541✔
87	}
88
NEW 89	std::pair<Direction, double> PolarAzimuthal::sample_as_bias(	×
90	uint64_t* seed) const
91	{
92	// Sample cosine of polar angle
NEW 93	auto [mu, mu_wgt] = mu_->sample(seed);	×
94
95	// Sample azimuthal angle
NEW 96	auto [phi, phi_wgt] = phi_->sample(seed);	×
97
NEW 98	double pdf_evaluation = mu_->evaluate(mu) * phi_->evaluate(phi);	×
99
NEW 100	if (mu == 1.0)	×
NEW 101	return {u_ref_, pdf_evaluation};	×
102
NEW 103	return {rotate_angle(u_ref_, mu, &phi, seed), pdf_evaluation};	×
104	}
105
106	//==============================================================================
107	// Isotropic implementation
108	//==============================================================================
109
110	Isotropic::Isotropic(pugi::xml_node node) : UnitSphereDistribution {node}	255✔
111	{
112	if (check_for_node(node, "bias")) {	255✔
NEW 113	pugi::xml_node bias_node = node.child("bias");	×
NEW 114	std::string bias_type = get_node_value(bias_node, "type", true, true);	×
NEW 115	if (bias_type != "mu-phi") {	×
NEW 116	openmc::fatal_error(	×
117	"Isotropic distributions may only be biased by a PolarAzimuthal.");
118	}
NEW 119	auto bias = std::make_unique<PolarAzimuthal>(bias_node);	×
NEW 120	if (bias->mu()->bias() \|\| bias->phi()->bias()) {	×
NEW 121	openmc::fatal_error(	×
122	"Attempted to bias Isotropic distribution with a biased PolarAzimuthal "
123	"distribution. Please ensure bias distributions are unbiased.");
124	}
NEW 125	this->set_bias(std::move(bias));	×
126	}
127	}	255✔
128
129	Direction isotropic_direction(uint64_t* seed)	111,314,611✔
130	{
131	double phi = uniform_distribution(0., 2.0 * PI, seed);	111,314,611✔
132	double mu = uniform_distribution(-1., 1., seed);	111,314,611✔
133	return {mu, std::sqrt(1.0 - mu * mu) * std::cos(phi),	111,314,611✔
134	std::sqrt(1.0 - mu * mu) * std::sin(phi)};	111,314,611✔
135	}
136
137	std::pair<Direction, double> Isotropic::sample(uint64_t* seed) const	21,614,969✔
138	{
139	if (bias()) {	21,614,969✔
NEW 140	auto [val, eval] = bias()->sample_as_bias(seed);	×
NEW 141	return {val, 1.0 / (4.0 * PI * eval)};	×
142	} else {
143	return {isotropic_direction(seed), 1.0};	21,614,969✔
144	}
145	}
146
147	//==============================================================================
148	// Monodirectional implementation
149	//==============================================================================
150
151	std::pair<Direction, double> Monodirectional::sample(uint64_t* seed) const	12,466,214✔
152	{
153	return {u_ref_, 1.0};	12,466,214✔
154	}
155
156	} // namespace openmc

openmc-dev / openmc / 17445237309

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