15666011990

Committed 15 Jun 2025 06:12PM UTC coverage: 85.117% (-0.05%) from 85.162%

Build # 15666011990

Build Type

Pull #3444

github

Committed by

web-flow

Commit Message

Merge d5e5b6009 into 23eab2c89

Pull Request Pull Request #3444: Add Ball spatial distribution

Coverage Stats

13 of 48 new or added lines in 2 files covered. (27.08%)

52386 of 61546 relevant lines covered (85.12%)

36738320.84 hits per line

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

71.01

/src/distribution_spatial.cpp

#include "openmc/distribution_spatial.h"

#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/mesh.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/xml_interface.h"

namespace openmc {

//==============================================================================
// SpatialDistribution implementation
//==============================================================================

unique_ptr<SpatialDistribution> SpatialDistribution::create(pugi::xml_node node)
{
  // Check for type of spatial distribution and read
  std::string type;
  if (check_for_node(node, "type"))
    type = get_node_value(node, "type", true, true);
  if (type == "cartesian") {
    return UPtrSpace {new CartesianIndependent(node)};
  } else if (type == "cylindrical") {
    return UPtrSpace {new CylindricalIndependent(node)};
  } else if (type == "spherical") {
    return UPtrSpace {new SphericalIndependent(node)};
  } else if (type == "mesh") {
    return UPtrSpace {new MeshSpatial(node)};
  } else if (type == "cloud") {
    return UPtrSpace {new PointCloud(node)};
  } else if (type == "box") {
    return UPtrSpace {new SpatialBox(node)};
  } else if (type == "fission") {
    return UPtrSpace {new SpatialBox(node, true)};
  } else if (type == "point") {
    return UPtrSpace {new SpatialPoint(node)};
  } else if (type == "ball") {
    return UPtrSpace {new SpatialBall(node)};
  } else {
    fatal_error(fmt::format(
      "Invalid spatial distribution for external source: {}", type));
  }
}

//==============================================================================
// CartesianIndependent implementation
//==============================================================================

CartesianIndependent::CartesianIndependent(pugi::xml_node node)
{
  // Read distribution for x coordinate
  if (check_for_node(node, "x")) {
    pugi::xml_node node_dist = node.child("x");
    x_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at x=0
    double x[] {0.0};
    double p[] {1.0};
    x_ = UPtrDist {new Discrete {x, p, 1}};
  }

  // Read distribution for y coordinate
  if (check_for_node(node, "y")) {
    pugi::xml_node node_dist = node.child("y");
    y_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at y=0
    double x[] {0.0};
    double p[] {1.0};
    y_ = UPtrDist {new Discrete {x, p, 1}};
  }

  // Read distribution for z coordinate
  if (check_for_node(node, "z")) {
    pugi::xml_node node_dist = node.child("z");
    z_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at z=0
    double x[] {0.0};
    double p[] {1.0};
    z_ = UPtrDist {new Discrete {x, p, 1}};
  }
}

Position CartesianIndependent::sample(uint64_t* seed) const
{
  return {x_->sample(seed), y_->sample(seed), z_->sample(seed)};
}

//==============================================================================
// CylindricalIndependent implementation
//==============================================================================

CylindricalIndependent::CylindricalIndependent(pugi::xml_node node)
{
  // Read distribution for r-coordinate
  if (check_for_node(node, "r")) {
    pugi::xml_node node_dist = node.child("r");
    r_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at r=0
    double x[] {0.0};
    double p[] {1.0};
    r_ = make_unique<Discrete>(x, p, 1);
  }

  // Read distribution for phi-coordinate
  if (check_for_node(node, "phi")) {
    pugi::xml_node node_dist = node.child("phi");
    phi_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at phi=0
    double x[] {0.0};
    double p[] {1.0};
    phi_ = make_unique<Discrete>(x, p, 1);
  }

  // Read distribution for z-coordinate
  if (check_for_node(node, "z")) {
    pugi::xml_node node_dist = node.child("z");
    z_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at z=0
    double x[] {0.0};
    double p[] {1.0};
    z_ = make_unique<Discrete>(x, p, 1);
  }

  // Read cylinder center coordinates
  if (check_for_node(node, "origin")) {
    auto origin = get_node_array<double>(node, "origin");
    if (origin.size() == 3) {
      origin_ = origin;
    } else {
      fatal_error(
        "Origin for cylindrical source distribution must be length 3");
    }
  } else {
    // If no coordinates were specified, default to (0, 0, 0)
    origin_ = {0.0, 0.0, 0.0};
  }
}

Position CylindricalIndependent::sample(uint64_t* seed) const
{
  double r = r_->sample(seed);
  double phi = phi_->sample(seed);
  double x = r * cos(phi) + origin_.x;
  double y = r * sin(phi) + origin_.y;
  double z = z_->sample(seed) + origin_.z;
  return {x, y, z};
}

//==============================================================================
// SphericalIndependent implementation
//==============================================================================

SphericalIndependent::SphericalIndependent(pugi::xml_node node)
{
  // Read distribution for r-coordinate
  if (check_for_node(node, "r")) {
    pugi::xml_node node_dist = node.child("r");
    r_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at r=0
    double x[] {0.0};
    double p[] {1.0};
    r_ = make_unique<Discrete>(x, p, 1);
  }

  // Read distribution for cos_theta-coordinate
  if (check_for_node(node, "cos_theta")) {
    pugi::xml_node node_dist = node.child("cos_theta");
    cos_theta_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at
    // cos_theta=0
    double x[] {0.0};
    double p[] {1.0};
    cos_theta_ = make_unique<Discrete>(x, p, 1);
  }

  // Read distribution for phi-coordinate
  if (check_for_node(node, "phi")) {
    pugi::xml_node node_dist = node.child("phi");
    phi_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at phi=0
    double x[] {0.0};
    double p[] {1.0};
    phi_ = make_unique<Discrete>(x, p, 1);
  }

  // Read sphere center coordinates
  if (check_for_node(node, "origin")) {
    auto origin = get_node_array<double>(node, "origin");
    if (origin.size() == 3) {
      origin_ = origin;
    } else {
      fatal_error("Origin for spherical source distribution must be length 3");
    }
  } else {
    // If no coordinates were specified, default to (0, 0, 0)
    origin_ = {0.0, 0.0, 0.0};
  }
}

Position SphericalIndependent::sample(uint64_t* seed) const
{
  double r = r_->sample(seed);
  double cos_theta = cos_theta_->sample(seed);
  double phi = phi_->sample(seed);
  // sin(theta) by sin**2 + cos**2 = 1
  double x = r * std::sqrt(1 - cos_theta * cos_theta) * cos(phi) + origin_.x;
  double y = r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y;
  double z = r * cos_theta + origin_.z;
  return {x, y, z};
}

//==============================================================================
// MeshSpatial implementation
//==============================================================================

MeshSpatial::MeshSpatial(pugi::xml_node node)
{

  if (get_node_value(node, "type", true, true) != "mesh") {
    fatal_error(fmt::format(
      "Incorrect spatial type '{}' for a MeshSpatial distribution"));
  }

  // No in-tet distributions implemented, could include distributions for the
  // barycentric coords Read in unstructured mesh from mesh_id value
  int32_t mesh_id = std::stoi(get_node_value(node, "mesh_id"));
  // Get pointer to spatial distribution
  mesh_idx_ = model::mesh_map.at(mesh_id);

  const auto mesh_ptr = model::meshes.at(mesh_idx_).get();

  check_element_types();

  size_t n_bins = this->n_sources();
  std::vector<double> strengths(n_bins, 1.0);

  // Create cdfs for sampling for an element over a mesh
  // Volume scheme is weighted by the volume of each tet
  // File scheme is weighted by an array given in the xml file
  if (check_for_node(node, "strengths")) {
    strengths = get_node_array<double>(node, "strengths");
    if (strengths.size() != n_bins) {
      fatal_error(
        fmt::format("Number of entries in the source strengths array {} does "
                    "not match the number of entities in mesh {} ({}).",
          strengths.size(), mesh_id, n_bins));
    }
  }

  if (get_node_value_bool(node, "volume_normalized")) {
    for (int i = 0; i < n_bins; i++) {
      strengths[i] *= this->mesh()->volume(i);
    }
  }

  elem_idx_dist_.assign(strengths);
}

MeshSpatial::MeshSpatial(int32_t mesh_idx, span<const double> strengths)
  : mesh_idx_(mesh_idx)
{
  check_element_types();
  elem_idx_dist_.assign(strengths);
}

void MeshSpatial::check_element_types() const
{
  const auto umesh_ptr = dynamic_cast<const UnstructuredMesh*>(this->mesh());
  if (umesh_ptr) {
    // ensure that the unstructured mesh contains only linear tets
    for (int bin = 0; bin < umesh_ptr->n_bins(); bin++) {
      if (umesh_ptr->element_type(bin) != ElementType::LINEAR_TET) {
        fatal_error(
          "Mesh specified for source must contain only linear tetrahedra.");
      }
    }
  }
}

int32_t MeshSpatial::sample_element_index(uint64_t* seed) const
{
  return elem_idx_dist_.sample(seed);
}

std::pair<int32_t, Position> MeshSpatial::sample_mesh(uint64_t* seed) const
{
  // Sample the CDF defined in initialization above
  int32_t elem_idx = this->sample_element_index(seed);
  return {elem_idx, mesh()->sample_element(elem_idx, seed)};
}

Position MeshSpatial::sample(uint64_t* seed) const
{
  return this->sample_mesh(seed).second;
}

//==============================================================================
// PointCloud implementation
//==============================================================================

PointCloud::PointCloud(pugi::xml_node node)
{
  if (check_for_node(node, "coords")) {
    point_cloud_ = get_node_position_array(node, "coords");
  } else {
    fatal_error("No coordinates were provided for the PointCloud "
                "spatial distribution");
  }

  std::vector<double> strengths;

  if (check_for_node(node, "strengths"))
    strengths = get_node_array<double>(node, "strengths");
  else
    strengths.resize(point_cloud_.size(), 1.0);

  if (strengths.size() != point_cloud_.size()) {
    fatal_error(
      fmt::format("Number of entries for the strengths array {} does "
                  "not match the number of spatial points provided {}.",
        strengths.size(), point_cloud_.size()));
  }

  point_idx_dist_.assign(strengths);
}

PointCloud::PointCloud(
  std::vector<Position> point_cloud, span<const double> strengths)
{
  point_cloud_.assign(point_cloud.begin(), point_cloud.end());
  point_idx_dist_.assign(strengths);
}

Position PointCloud::sample(uint64_t* seed) const
{
  int32_t index = point_idx_dist_.sample(seed);
  return point_cloud_[index];
}

//==============================================================================
// SpatialBox implementation
//==============================================================================

SpatialBox::SpatialBox(pugi::xml_node node, bool fission)
  : only_fissionable_ {fission}
{
  // Read lower-right/upper-left coordinates
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() != 6)
    openmc::fatal_error("Box/fission spatial source must have six "
                        "parameters specified.");

  lower_left_ = Position {params[0], params[1], params[2]};
  upper_right_ = Position {params[3], params[4], params[5]};
}

Position SpatialBox::sample(uint64_t* seed) const
{
  Position xi {prn(seed), prn(seed), prn(seed)};
  return lower_left_ + xi * (upper_right_ - lower_left_);
}

//==============================================================================
// SpatialBall implementation
//==============================================================================

SpatialBall::SpatialBall(pugi::xml_node node)
{
  // Read origin coordinates and radius
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() != 4)
    openmc::fatal_error("Ball spatial source must have four "
                        "parameters specified.");

  origin_ = Position {params[0], params[1], params[2]};
  radius_ = params[3];
}

Position SpatialBall::sample(uint64_t* seed) const
{
  double u = 2.0 * prn(seed) - 1.0;
  double phi = 2 * PI * prn(seed);
  double r = radius_ * std::cbrt(prn(seed));
  Position xi {std::cos(phi) * std::sqrt(1 - u * u),
    std::sin(phi) * std::sqrt(1 - u * u), u};
  return origin_ + xi * r;
}

//==============================================================================
// SpatialPoint implementation
//==============================================================================

SpatialPoint::SpatialPoint(pugi::xml_node node)
{
  // Read location of point source
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() != 3)
    openmc::fatal_error("Point spatial source must have three "
                        "parameters specified.");

  // Set position
  r_ = Position {params.data()};
}

Position SpatialPoint::sample(uint64_t* seed) const
{
  return r_;
}

} // namespace openmc

1	#include "openmc/distribution_spatial.h"
2
3	#include "openmc/constants.h"
4	#include "openmc/error.h"
5	#include "openmc/mesh.h"
6	#include "openmc/random_lcg.h"
7	#include "openmc/search.h"
8	#include "openmc/xml_interface.h"
9
10	namespace openmc {
11
12	//==============================================================================
13	// SpatialDistribution implementation
14	//==============================================================================
15
16	unique_ptr<SpatialDistribution> SpatialDistribution::create(pugi::xml_node node)	6,156✔
17	{
18	// Check for type of spatial distribution and read
19	std::string type;	6,156✔
20	if (check_for_node(node, "type"))	6,156✔
21	type = get_node_value(node, "type", true, true);	6,156✔
22	if (type == "cartesian") {	6,156✔
23	return UPtrSpace {new CartesianIndependent(node)};	16✔
24	} else if (type == "cylindrical") {	6,140✔
25	return UPtrSpace {new CylindricalIndependent(node)};	48✔
26	} else if (type == "spherical") {	6,092✔
27	return UPtrSpace {new SphericalIndependent(node)};	63✔
28	} else if (type == "mesh") {	6,029✔
29	return UPtrSpace {new MeshSpatial(node)};	40✔
30	} else if (type == "cloud") {	5,989✔
31	return UPtrSpace {new PointCloud(node)};	11✔
32	} else if (type == "box") {	5,978✔
33	return UPtrSpace {new SpatialBox(node)};	3,354✔
34	} else if (type == "fission") {	2,624✔
35	return UPtrSpace {new SpatialBox(node, true)};	26✔
36	} else if (type == "point") {	2,598✔
37	return UPtrSpace {new SpatialPoint(node)};	2,598✔
NEW 38	} else if (type == "ball") {	×
NEW 39	return UPtrSpace {new SpatialBall(node)};	×
40	} else {
41	fatal_error(fmt::format(	×
42	"Invalid spatial distribution for external source: {}", type));
43	}
44	}	6,155✔
45
46	//==============================================================================
47	// CartesianIndependent implementation
48	//==============================================================================
49
50	CartesianIndependent::CartesianIndependent(pugi::xml_node node)	16✔
51	{
52	// Read distribution for x coordinate
53	if (check_for_node(node, "x")) {	16✔
54	pugi::xml_node node_dist = node.child("x");	16✔
55	x_ = distribution_from_xml(node_dist);	16✔
56	} else {
57	// If no distribution was specified, default to a single point at x=0
58	double x[] {0.0};	×
59	double p[] {1.0};	×
60	x_ = UPtrDist {new Discrete {x, p, 1}};	×
61	}
62
63	// Read distribution for y coordinate
64	if (check_for_node(node, "y")) {	16✔
65	pugi::xml_node node_dist = node.child("y");	16✔
66	y_ = distribution_from_xml(node_dist);	16✔
67	} else {
68	// If no distribution was specified, default to a single point at y=0
69	double x[] {0.0};	×
70	double p[] {1.0};	×
71	y_ = UPtrDist {new Discrete {x, p, 1}};	×
72	}
73
74	// Read distribution for z coordinate
75	if (check_for_node(node, "z")) {	16✔
76	pugi::xml_node node_dist = node.child("z");	16✔
77	z_ = distribution_from_xml(node_dist);	16✔
78	} else {
79	// If no distribution was specified, default to a single point at z=0
80	double x[] {0.0};	×
81	double p[] {1.0};	×
82	z_ = UPtrDist {new Discrete {x, p, 1}};	×
83	}
84	}	16✔
85
86	Position CartesianIndependent::sample(uint64_t* seed) const	3,190✔
87	{
88	return {x_->sample(seed), y_->sample(seed), z_->sample(seed)};	3,190✔
89	}
90
91	//==============================================================================
92	// CylindricalIndependent implementation
93	//==============================================================================
94
95	CylindricalIndependent::CylindricalIndependent(pugi::xml_node node)	48✔
96	{
97	// Read distribution for r-coordinate
98	if (check_for_node(node, "r")) {	48✔
99	pugi::xml_node node_dist = node.child("r");	48✔
100	r_ = distribution_from_xml(node_dist);	48✔
101	} else {
102	// If no distribution was specified, default to a single point at r=0
103	double x[] {0.0};	×
104	double p[] {1.0};	×
105	r_ = make_unique<Discrete>(x, p, 1);	×
106	}
107
108	// Read distribution for phi-coordinate
109	if (check_for_node(node, "phi")) {	48✔
110	pugi::xml_node node_dist = node.child("phi");	48✔
111	phi_ = distribution_from_xml(node_dist);	48✔
112	} else {
113	// If no distribution was specified, default to a single point at phi=0
114	double x[] {0.0};	×
115	double p[] {1.0};	×
116	phi_ = make_unique<Discrete>(x, p, 1);	×
117	}
118
119	// Read distribution for z-coordinate
120	if (check_for_node(node, "z")) {	48✔
121	pugi::xml_node node_dist = node.child("z");	48✔
122	z_ = distribution_from_xml(node_dist);	48✔
123	} else {
124	// If no distribution was specified, default to a single point at z=0
125	double x[] {0.0};	×
126	double p[] {1.0};	×
127	z_ = make_unique<Discrete>(x, p, 1);	×
128	}
129
130	// Read cylinder center coordinates
131	if (check_for_node(node, "origin")) {	48✔
132	auto origin = get_node_array<double>(node, "origin");	48✔
133	if (origin.size() == 3) {	48✔
134	origin_ = origin;	48✔
135	} else {
136	fatal_error(	×
137	"Origin for cylindrical source distribution must be length 3");
138	}
139	} else {	48✔
140	// If no coordinates were specified, default to (0, 0, 0)
141	origin_ = {0.0, 0.0, 0.0};	×
142	}
143	}	48✔
144
145	Position CylindricalIndependent::sample(uint64_t* seed) const	3,542✔
146	{
147	double r = r_->sample(seed);	3,542✔
148	double phi = phi_->sample(seed);	3,542✔
149	double x = r * cos(phi) + origin_.x;	3,542✔
150	double y = r * sin(phi) + origin_.y;	3,542✔
151	double z = z_->sample(seed) + origin_.z;	3,542✔
152	return {x, y, z};	3,542✔
153	}
154
155	//==============================================================================
156	// SphericalIndependent implementation
157	//==============================================================================
158
159	SphericalIndependent::SphericalIndependent(pugi::xml_node node)	63✔
160	{
161	// Read distribution for r-coordinate
162	if (check_for_node(node, "r")) {	63✔
163	pugi::xml_node node_dist = node.child("r");	63✔
164	r_ = distribution_from_xml(node_dist);	63✔
165	} else {
166	// If no distribution was specified, default to a single point at r=0
167	double x[] {0.0};	×
168	double p[] {1.0};	×
169	r_ = make_unique<Discrete>(x, p, 1);	×
170	}
171
172	// Read distribution for cos_theta-coordinate
173	if (check_for_node(node, "cos_theta")) {	63✔
174	pugi::xml_node node_dist = node.child("cos_theta");	63✔
175	cos_theta_ = distribution_from_xml(node_dist);	63✔
176	} else {
177	// If no distribution was specified, default to a single point at
178	// cos_theta=0
179	double x[] {0.0};	×
180	double p[] {1.0};	×
181	cos_theta_ = make_unique<Discrete>(x, p, 1);	×
182	}
183
184	// Read distribution for phi-coordinate
185	if (check_for_node(node, "phi")) {	63✔
186	pugi::xml_node node_dist = node.child("phi");	63✔
187	phi_ = distribution_from_xml(node_dist);	63✔
188	} else {
189	// If no distribution was specified, default to a single point at phi=0
190	double x[] {0.0};	×
191	double p[] {1.0};	×
192	phi_ = make_unique<Discrete>(x, p, 1);	×
193	}
194
195	// Read sphere center coordinates
196	if (check_for_node(node, "origin")) {	63✔
197	auto origin = get_node_array<double>(node, "origin");	63✔
198	if (origin.size() == 3) {	63✔
199	origin_ = origin;	63✔
200	} else {
201	fatal_error("Origin for spherical source distribution must be length 3");	×
202	}
203	} else {	63✔
204	// If no coordinates were specified, default to (0, 0, 0)
205	origin_ = {0.0, 0.0, 0.0};	×
206	}
207	}	63✔
208
209	Position SphericalIndependent::sample(uint64_t* seed) const	182,288✔
210	{
211	double r = r_->sample(seed);	182,288✔
212	double cos_theta = cos_theta_->sample(seed);	182,288✔
213	double phi = phi_->sample(seed);	182,288✔
214	// sin(theta) by sin2 + cos2 = 1
215	double x = r * std::sqrt(1 - cos_theta * cos_theta) * cos(phi) + origin_.x;	182,288✔
216	double y = r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y;	182,288✔
217	double z = r * cos_theta + origin_.z;	182,288✔
218	return {x, y, z};	182,288✔
219	}
220
221	//==============================================================================
222	// MeshSpatial implementation
223	//==============================================================================
224
225	MeshSpatial::MeshSpatial(pugi::xml_node node)	40✔
226	{
227
228	if (get_node_value(node, "type", true, true) != "mesh") {	40✔
229	fatal_error(fmt::format(	×
230	"Incorrect spatial type '{}' for a MeshSpatial distribution"));
231	}
232
233	// No in-tet distributions implemented, could include distributions for the
234	// barycentric coords Read in unstructured mesh from mesh_id value
235	int32_t mesh_id = std::stoi(get_node_value(node, "mesh_id"));	40✔
236	// Get pointer to spatial distribution
237	mesh_idx_ = model::mesh_map.at(mesh_id);	40✔
238
239	const auto mesh_ptr = model::meshes.at(mesh_idx_).get();	40✔
240
241	check_element_types();	40✔
242
243	size_t n_bins = this->n_sources();	40✔
244	std::vector<double> strengths(n_bins, 1.0);	40✔
245
246	// Create cdfs for sampling for an element over a mesh
247	// Volume scheme is weighted by the volume of each tet
248	// File scheme is weighted by an array given in the xml file
249	if (check_for_node(node, "strengths")) {	40✔
250	strengths = get_node_array<double>(node, "strengths");	37✔
251	if (strengths.size() != n_bins) {	37✔
252	fatal_error(	1✔
253	fmt::format("Number of entries in the source strengths array {} does "	1✔
254	"not match the number of entities in mesh {} ({}).",
255	strengths.size(), mesh_id, n_bins));	1✔
256	}
257	}
258
259	if (get_node_value_bool(node, "volume_normalized")) {	39✔
260	for (int i = 0; i < n_bins; i++) {	37,356✔
261	strengths[i] *= this->mesh()->volume(i);	37,320✔
262	}
263	}
264
265	elem_idx_dist_.assign(strengths);	39✔
266	}	39✔
267
268	MeshSpatial::MeshSpatial(int32_t mesh_idx, span<const double> strengths)	201✔
269	: mesh_idx_(mesh_idx)	201✔
270	{
271	check_element_types();	201✔
272	elem_idx_dist_.assign(strengths);	201✔
273	}	201✔
274
275	void MeshSpatial::check_element_types() const	241✔
276	{
277	const auto umesh_ptr = dynamic_cast<const UnstructuredMesh*>(this->mesh());	241✔
278	if (umesh_ptr) {	241✔
279	// ensure that the unstructured mesh contains only linear tets
280	for (int bin = 0; bin < umesh_ptr->n_bins(); bin++) {	120,010✔
281	if (umesh_ptr->element_type(bin) != ElementType::LINEAR_TET) {	120,000✔
282	fatal_error(	×
283	"Mesh specified for source must contain only linear tetrahedra.");
284	}
285	}
286	}
287	}	241✔
288
289	int32_t MeshSpatial::sample_element_index(uint64_t* seed) const	788,670✔
290	{
291	return elem_idx_dist_.sample(seed);	788,670✔
292	}
293
294	std::pair<int32_t, Position> MeshSpatial::sample_mesh(uint64_t* seed) const	601,530✔
295	{
296	// Sample the CDF defined in initialization above
297	int32_t elem_idx = this->sample_element_index(seed);	601,530✔
298	return {elem_idx, mesh()->sample_element(elem_idx, seed)};	601,530✔
299	}
300
301	Position MeshSpatial::sample(uint64_t* seed) const	601,530✔
302	{
303	return this->sample_mesh(seed).second;	601,530✔
304	}
305
306	//==============================================================================
307	// PointCloud implementation
308	//==============================================================================
309
310	PointCloud::PointCloud(pugi::xml_node node)	11✔
311	{
312	if (check_for_node(node, "coords")) {	11✔
313	point_cloud_ = get_node_position_array(node, "coords");	11✔
314	} else {
315	fatal_error("No coordinates were provided for the PointCloud "	×
316	"spatial distribution");
317	}
318
319	std::vector<double> strengths;	11✔
320
321	if (check_for_node(node, "strengths"))	11✔
322	strengths = get_node_array<double>(node, "strengths");	11✔
323	else
324	strengths.resize(point_cloud_.size(), 1.0);	×
325
326	if (strengths.size() != point_cloud_.size()) {	11✔
327	fatal_error(	×
328	fmt::format("Number of entries for the strengths array {} does "	×
329	"not match the number of spatial points provided {}.",
330	strengths.size(), point_cloud_.size()));	×
331	}
332
333	point_idx_dist_.assign(strengths);	11✔
334	}	11✔
335
336	PointCloud::PointCloud(	×
337	std::vector<Position> point_cloud, span<const double> strengths)	×
338	{
339	point_cloud_.assign(point_cloud.begin(), point_cloud.end());	×
340	point_idx_dist_.assign(strengths);	×
341	}
342
343	Position PointCloud::sample(uint64_t* seed) const	550,000✔
344	{
345	int32_t index = point_idx_dist_.sample(seed);	550,000✔
346	return point_cloud_[index];	550,000✔
347	}
348
349	//==============================================================================
350	// SpatialBox implementation
351	//==============================================================================
352
353	SpatialBox::SpatialBox(pugi::xml_node node, bool fission)	3,380✔
354	: only_fissionable_ {fission}	3,380✔
355	{
356	// Read lower-right/upper-left coordinates
357	auto params = get_node_array<double>(node, "parameters");	3,380✔
358	if (params.size() != 6)	3,380✔
359	openmc::fatal_error("Box/fission spatial source must have six "	×
360	"parameters specified.");
361
362	lower_left_ = Position {params[0], params[1], params[2]};	3,380✔
363	upper_right_ = Position {params[3], params[4], params[5]};	3,380✔
364	}	3,380✔
365
366	Position SpatialBox::sample(uint64_t* seed) const	5,140,167✔
367	{
368	Position xi {prn(seed), prn(seed), prn(seed)};	5,140,167✔
369	return lower_left_ + xi * (upper_right_ - lower_left_);	10,280,334✔
370	}
371
372	//==============================================================================
373	// SpatialBall implementation
374	//==============================================================================
375
NEW 376	SpatialBall::SpatialBall(pugi::xml_node node)	×
377	{
378	// Read origin coordinates and radius
NEW 379	auto params = get_node_array<double>(node, "parameters");	×
NEW 380	if (params.size() != 4)	×
NEW 381	openmc::fatal_error("Ball spatial source must have four "	×
382	"parameters specified.");
383
NEW 384	origin_ = Position {params[0], params[1], params[2]};	×
NEW 385	radius_ = params[3];	×
386	}
387
NEW 388	Position SpatialBall::sample(uint64_t* seed) const	×
389	{
NEW 390	double u = 2.0 * prn(seed) - 1.0;	×
NEW 391	double phi = 2 * PI * prn(seed);	×
NEW 392	double r = radius_ * std::cbrt(prn(seed));	×
NEW 393	Position xi {std::cos(phi) * std::sqrt(1 - u * u),	×
NEW 394	std::sin(phi) * std::sqrt(1 - u * u), u};	×
NEW 395	return origin_ + xi * r;	×
396	}
397
398	//==============================================================================
399	// SpatialPoint implementation
400	//==============================================================================
401
402	SpatialPoint::SpatialPoint(pugi::xml_node node)	2,598✔
403	{
404	// Read location of point source
405	auto params = get_node_array<double>(node, "parameters");	2,598✔
406	if (params.size() != 3)	2,598✔
407	openmc::fatal_error("Point spatial source must have three "	×
408	"parameters specified.");
409
410	// Set position
411	r_ = Position {params.data()};	2,598✔
412	}	2,598✔
413
414	Position SpatialPoint::sample(uint64_t* seed) const	22,533,234✔
415	{
416	return r_;	22,533,234✔
417	}
418
419	} // namespace openmc

openmc-dev / openmc / 15666011990

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