16033463050

Committed 02 Jul 2025 06:57PM UTC coverage: 85.256% (+0.02%) from 85.241%

Build # 16033463050

Build Type

Pull #3474

github

Committed by

web-flow

Commit Message

Merge e9cd2e48b into eb74d497d

Pull Request Pull Request #3474: added option to have cylindrical sources with general axis

Run Details

57 of 65 new or added lines in 2 files covered. (87.69%)

151 existing lines in 7 files now uncovered.

52702 of 61816 relevant lines covered (85.26%)

36453600.67 hits per line

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

76.3

/src/distribution_spatial.cpp

#include "openmc/distribution_spatial.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 {
    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};
  }

  // Read cylinder z_dir
  if (check_for_node(node, "z_dir")) {
    auto z_dir = get_node_array<double>(node, "z_dir");
    if (z_dir.size() == 3) {
      z_dir_ = z_dir;
      z_dir_ /= z_dir_.norm();
    } else {
      fatal_error("z_dir for cylindrical source distribution must be length 3");
    }
  } else {
    // If no z_dir was specified, default to (0, 0, 1)
    z_dir_ = {0.0, 0.0, 1.0};
  }

  // Read cylinder r_dir
  if (check_for_node(node, "r_dir")) {
    auto r_dir = get_node_array<double>(node, "r_dir");
    if (r_dir.size() == 3) {
      r_dir_ = r_dir;
      r_dir_ /= r_dir_.norm();
    } else {
      fatal_error("r_dir for cylindrical source distribution must be length 3");
    }
  } else {
    // If no r_dir was specified, default to (1, 0, 0)
    r_dir_ = {1.0, 0.0, 0.0};
  }

  if (r_dir_.dot(z_dir_) > 1e-12)
    fatal_error("r_dir must be perpendicular to z_dir");

  auto phi_dir = z_dir_.cross(r_dir_);
  phi_dir /= phi_dir.norm();
  phi_dir_ = phi_dir;
}

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

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

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

openmc-dev / openmc / 16033463050

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