• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 18299973882

07 Oct 2025 02:14AM UTC coverage: 81.92% (-3.3%) from 85.194%
18299973882

push

github

web-flow
Switch to using coveralls github action for reporting (#3594)

16586 of 23090 branches covered (71.83%)

Branch coverage included in aggregate %.

53703 of 62712 relevant lines covered (85.63%)

43428488.21 hits per line

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

72.4
/src/distribution_spatial.cpp
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,663✔
16
{
17
  // Check for type of spatial distribution and read
18
  std::string type;
6,663✔
19
  if (check_for_node(node, "type"))
6,663!
20
    type = get_node_value(node, "type", true, true);
6,663✔
21
  if (type == "cartesian") {
6,663✔
22
    return UPtrSpace {new CartesianIndependent(node)};
16✔
23
  } else if (type == "cylindrical") {
6,647✔
24
    return UPtrSpace {new CylindricalIndependent(node)};
48✔
25
  } else if (type == "spherical") {
6,599✔
26
    return UPtrSpace {new SphericalIndependent(node)};
63✔
27
  } else if (type == "mesh") {
6,536✔
28
    return UPtrSpace {new MeshSpatial(node)};
46✔
29
  } else if (type == "cloud") {
6,490✔
30
    return UPtrSpace {new PointCloud(node)};
13✔
31
  } else if (type == "box") {
6,477✔
32
    return UPtrSpace {new SpatialBox(node)};
3,663✔
33
  } else if (type == "fission") {
2,814✔
34
    return UPtrSpace {new SpatialBox(node, true)};
26✔
35
  } else if (type == "point") {
2,788!
36
    return UPtrSpace {new SpatialPoint(node)};
2,788✔
37
  } else {
38
    fatal_error(fmt::format(
×
39
      "Invalid spatial distribution for external source: {}", type));
40
  }
41
}
6,662✔
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
}
48✔
141

142
Position CylindricalIndependent::sample(uint64_t* seed) const
3,542✔
143
{
144
  double r = r_->sample(seed);
3,542✔
145
  double phi = phi_->sample(seed);
3,542✔
146
  double x = r * cos(phi) + origin_.x;
3,542✔
147
  double y = r * sin(phi) + origin_.y;
3,542✔
148
  double z = z_->sample(seed) + origin_.z;
3,542✔
149
  return {x, y, z};
3,542✔
150
}
151

152
//==============================================================================
153
// SphericalIndependent implementation
154
//==============================================================================
155

156
SphericalIndependent::SphericalIndependent(pugi::xml_node node)
63✔
157
{
158
  // Read distribution for r-coordinate
159
  if (check_for_node(node, "r")) {
63!
160
    pugi::xml_node node_dist = node.child("r");
63✔
161
    r_ = distribution_from_xml(node_dist);
63✔
162
  } else {
163
    // If no distribution was specified, default to a single point at r=0
164
    double x[] {0.0};
×
165
    double p[] {1.0};
×
166
    r_ = make_unique<Discrete>(x, p, 1);
×
167
  }
168

169
  // Read distribution for cos_theta-coordinate
170
  if (check_for_node(node, "cos_theta")) {
63!
171
    pugi::xml_node node_dist = node.child("cos_theta");
63✔
172
    cos_theta_ = distribution_from_xml(node_dist);
63✔
173
  } else {
174
    // If no distribution was specified, default to a single point at
175
    // cos_theta=0
176
    double x[] {0.0};
×
177
    double p[] {1.0};
×
178
    cos_theta_ = make_unique<Discrete>(x, p, 1);
×
179
  }
180

181
  // Read distribution for phi-coordinate
182
  if (check_for_node(node, "phi")) {
63!
183
    pugi::xml_node node_dist = node.child("phi");
63✔
184
    phi_ = distribution_from_xml(node_dist);
63✔
185
  } else {
186
    // If no distribution was specified, default to a single point at phi=0
187
    double x[] {0.0};
×
188
    double p[] {1.0};
×
189
    phi_ = make_unique<Discrete>(x, p, 1);
×
190
  }
191

192
  // Read sphere center coordinates
193
  if (check_for_node(node, "origin")) {
63!
194
    auto origin = get_node_array<double>(node, "origin");
63✔
195
    if (origin.size() == 3) {
63!
196
      origin_ = origin;
63✔
197
    } else {
198
      fatal_error("Origin for spherical source distribution must be length 3");
×
199
    }
200
  } else {
63✔
201
    // If no coordinates were specified, default to (0, 0, 0)
202
    origin_ = {0.0, 0.0, 0.0};
×
203
  }
204
}
63✔
205

206
Position SphericalIndependent::sample(uint64_t* seed) const
182,288✔
207
{
208
  double r = r_->sample(seed);
182,288✔
209
  double cos_theta = cos_theta_->sample(seed);
182,288✔
210
  double phi = phi_->sample(seed);
182,288✔
211
  // sin(theta) by sin**2 + cos**2 = 1
212
  double x = r * std::sqrt(1 - cos_theta * cos_theta) * cos(phi) + origin_.x;
182,288✔
213
  double y = r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y;
182,288✔
214
  double z = r * cos_theta + origin_.z;
182,288✔
215
  return {x, y, z};
182,288✔
216
}
217

218
//==============================================================================
219
// MeshSpatial implementation
220
//==============================================================================
221

222
MeshSpatial::MeshSpatial(pugi::xml_node node)
46✔
223
{
224

225
  if (get_node_value(node, "type", true, true) != "mesh") {
46!
226
    fatal_error(fmt::format(
×
227
      "Incorrect spatial type '{}' for a MeshSpatial distribution"));
228
  }
229

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

236
  const auto mesh_ptr = model::meshes.at(mesh_idx_).get();
46✔
237

238
  check_element_types();
46✔
239

240
  size_t n_bins = this->n_sources();
46✔
241
  std::vector<double> strengths(n_bins, 1.0);
46✔
242

243
  // Create cdfs for sampling for an element over a mesh
244
  // Volume scheme is weighted by the volume of each tet
245
  // File scheme is weighted by an array given in the xml file
246
  if (check_for_node(node, "strengths")) {
46✔
247
    strengths = get_node_array<double>(node, "strengths");
43✔
248
    if (strengths.size() != n_bins) {
43✔
249
      fatal_error(
1✔
250
        fmt::format("Number of entries in the source strengths array {} does "
1!
251
                    "not match the number of entities in mesh {} ({}).",
252
          strengths.size(), mesh_id, n_bins));
1✔
253
    }
254
  }
255

256
  if (get_node_value_bool(node, "volume_normalized")) {
45✔
257
    for (int i = 0; i < n_bins; i++) {
37,602✔
258
      strengths[i] *= this->mesh()->volume(i);
37,560✔
259
    }
260
  }
261

262
  elem_idx_dist_.assign(strengths);
45✔
263
}
45✔
264

265
MeshSpatial::MeshSpatial(int32_t mesh_idx, span<const double> strengths)
205✔
266
  : mesh_idx_(mesh_idx)
205✔
267
{
268
  check_element_types();
205✔
269
  elem_idx_dist_.assign(strengths);
205✔
270
}
205✔
271

272
void MeshSpatial::check_element_types() const
251✔
273
{
274
  const auto umesh_ptr = dynamic_cast<const UnstructuredMesh*>(this->mesh());
251!
275
  if (umesh_ptr) {
251✔
276
    // ensure that the unstructured mesh contains only linear tets
277
    for (int bin = 0; bin < umesh_ptr->n_bins(); bin++) {
120,010✔
278
      if (umesh_ptr->element_type(bin) != ElementType::LINEAR_TET) {
120,000!
279
        fatal_error(
×
280
          "Mesh specified for source must contain only linear tetrahedra.");
281
      }
282
    }
283
  }
284
}
251✔
285

286
int32_t MeshSpatial::sample_element_index(uint64_t* seed) const
2,344,199✔
287
{
288
  return elem_idx_dist_.sample(seed);
2,344,199✔
289
}
290

291
std::pair<int32_t, Position> MeshSpatial::sample_mesh(uint64_t* seed) const
601,590✔
292
{
293
  // Sample the CDF defined in initialization above
294
  int32_t elem_idx = this->sample_element_index(seed);
601,590✔
295
  return {elem_idx, mesh()->sample_element(elem_idx, seed)};
601,590✔
296
}
297

298
Position MeshSpatial::sample(uint64_t* seed) const
601,590✔
299
{
300
  return this->sample_mesh(seed).second;
601,590✔
301
}
302

303
//==============================================================================
304
// PointCloud implementation
305
//==============================================================================
306

307
PointCloud::PointCloud(pugi::xml_node node)
13✔
308
{
309
  if (check_for_node(node, "coords")) {
13!
310
    point_cloud_ = get_node_position_array(node, "coords");
13✔
311
  } else {
312
    fatal_error("No coordinates were provided for the PointCloud "
×
313
                "spatial distribution");
314
  }
315

316
  std::vector<double> strengths;
13✔
317

318
  if (check_for_node(node, "strengths"))
13!
319
    strengths = get_node_array<double>(node, "strengths");
13✔
320
  else
321
    strengths.resize(point_cloud_.size(), 1.0);
×
322

323
  if (strengths.size() != point_cloud_.size()) {
13!
324
    fatal_error(
×
325
      fmt::format("Number of entries for the strengths array {} does "
×
326
                  "not match the number of spatial points provided {}.",
327
        strengths.size(), point_cloud_.size()));
×
328
  }
329

330
  point_idx_dist_.assign(strengths);
13✔
331
}
13✔
332

333
PointCloud::PointCloud(
×
334
  std::vector<Position> point_cloud, span<const double> strengths)
×
335
{
336
  point_cloud_.assign(point_cloud.begin(), point_cloud.end());
×
337
  point_idx_dist_.assign(strengths);
×
338
}
×
339

340
Position PointCloud::sample(uint64_t* seed) const
650,000✔
341
{
342
  int32_t index = point_idx_dist_.sample(seed);
650,000✔
343
  return point_cloud_[index];
650,000✔
344
}
345

346
//==============================================================================
347
// SpatialBox implementation
348
//==============================================================================
349

350
SpatialBox::SpatialBox(pugi::xml_node node, bool fission)
3,689✔
351
  : only_fissionable_ {fission}
3,689✔
352
{
353
  // Read lower-right/upper-left coordinates
354
  auto params = get_node_array<double>(node, "parameters");
3,689✔
355
  if (params.size() != 6)
3,689!
356
    openmc::fatal_error("Box/fission spatial source must have six "
×
357
                        "parameters specified.");
358

359
  lower_left_ = Position {params[0], params[1], params[2]};
3,689✔
360
  upper_right_ = Position {params[3], params[4], params[5]};
3,689✔
361
}
3,689✔
362

363
Position SpatialBox::sample(uint64_t* seed) const
5,243,300✔
364
{
365
  Position xi {prn(seed), prn(seed), prn(seed)};
5,243,300✔
366
  return lower_left_ + xi * (upper_right_ - lower_left_);
10,486,600✔
367
}
368

369
//==============================================================================
370
// SpatialPoint implementation
371
//==============================================================================
372

373
SpatialPoint::SpatialPoint(pugi::xml_node node)
2,788✔
374
{
375
  // Read location of point source
376
  auto params = get_node_array<double>(node, "parameters");
2,788✔
377
  if (params.size() != 3)
2,788!
378
    openmc::fatal_error("Point spatial source must have three "
×
379
                        "parameters specified.");
380

381
  // Set position
382
  r_ = Position {params.data()};
2,788✔
383
}
2,788✔
384

385
Position SpatialPoint::sample(uint64_t* seed) const
25,623,612✔
386
{
387
  return r_;
25,623,612✔
388
}
389

390
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc