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

openmc-dev / openmc / 12776996362

14 Jan 2025 09:49PM UTC coverage: 84.938% (+0.2%) from 84.729%
12776996362

Pull #3133

github

web-flow
Merge 0495246d9 into 549cc0973
Pull Request #3133: Kinetics parameters using Iterated Fission Probability

318 of 330 new or added lines in 10 files covered. (96.36%)

1658 existing lines in 66 files now uncovered.

50402 of 59340 relevant lines covered (84.94%)

33987813.96 hits per line

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

76.56
/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,191✔
16
{
17
  // Check for type of spatial distribution and read
18
  std::string type;
6,191✔
19
  if (check_for_node(node, "type"))
6,191✔
20
    type = get_node_value(node, "type", true, true);
6,191✔
21
  if (type == "cartesian") {
6,191✔
22
    return UPtrSpace {new CartesianIndependent(node)};
17✔
23
  } else if (type == "cylindrical") {
6,174✔
24
    return UPtrSpace {new CylindricalIndependent(node)};
51✔
25
  } else if (type == "spherical") {
6,123✔
26
    return UPtrSpace {new SphericalIndependent(node)};
65✔
27
  } else if (type == "mesh") {
6,058✔
28
    return UPtrSpace {new MeshSpatial(node)};
7✔
29
  } else if (type == "cloud") {
6,051✔
30
    return UPtrSpace {new PointCloud(node)};
12✔
31
  } else if (type == "box") {
6,039✔
32
    return UPtrSpace {new SpatialBox(node)};
3,360✔
33
  } else if (type == "fission") {
2,679✔
34
    return UPtrSpace {new SpatialBox(node, true)};
67✔
35
  } else if (type == "point") {
2,612✔
36
    return UPtrSpace {new SpatialPoint(node)};
2,612✔
37
  } else {
UNCOV
38
    fatal_error(fmt::format(
×
39
      "Invalid spatial distribution for external source: {}", type));
40
  }
41
}
6,190✔
42

43
//==============================================================================
44
// CartesianIndependent implementation
45
//==============================================================================
46

47
CartesianIndependent::CartesianIndependent(pugi::xml_node node)
17✔
48
{
49
  // Read distribution for x coordinate
50
  if (check_for_node(node, "x")) {
17✔
51
    pugi::xml_node node_dist = node.child("x");
17✔
52
    x_ = distribution_from_xml(node_dist);
17✔
53
  } else {
54
    // If no distribution was specified, default to a single point at x=0
55
    double x[] {0.0};
×
UNCOV
56
    double p[] {1.0};
×
UNCOV
57
    x_ = UPtrDist {new Discrete {x, p, 1}};
×
58
  }
59

60
  // Read distribution for y coordinate
61
  if (check_for_node(node, "y")) {
17✔
62
    pugi::xml_node node_dist = node.child("y");
17✔
63
    y_ = distribution_from_xml(node_dist);
17✔
64
  } else {
65
    // If no distribution was specified, default to a single point at y=0
66
    double x[] {0.0};
×
UNCOV
67
    double p[] {1.0};
×
UNCOV
68
    y_ = UPtrDist {new Discrete {x, p, 1}};
×
69
  }
70

71
  // Read distribution for z coordinate
72
  if (check_for_node(node, "z")) {
17✔
73
    pugi::xml_node node_dist = node.child("z");
17✔
74
    z_ = distribution_from_xml(node_dist);
17✔
75
  } else {
76
    // If no distribution was specified, default to a single point at z=0
77
    double x[] {0.0};
×
UNCOV
78
    double p[] {1.0};
×
UNCOV
79
    z_ = UPtrDist {new Discrete {x, p, 1}};
×
80
  }
81
}
17✔
82

83
Position CartesianIndependent::sample(uint64_t* seed) const
3,480✔
84
{
85
  return {x_->sample(seed), y_->sample(seed), z_->sample(seed)};
3,480✔
86
}
87

88
//==============================================================================
89
// CylindricalIndependent implementation
90
//==============================================================================
91

92
CylindricalIndependent::CylindricalIndependent(pugi::xml_node node)
51✔
93
{
94
  // Read distribution for r-coordinate
95
  if (check_for_node(node, "r")) {
51✔
96
    pugi::xml_node node_dist = node.child("r");
51✔
97
    r_ = distribution_from_xml(node_dist);
51✔
98
  } else {
99
    // If no distribution was specified, default to a single point at r=0
100
    double x[] {0.0};
×
UNCOV
101
    double p[] {1.0};
×
UNCOV
102
    r_ = make_unique<Discrete>(x, p, 1);
×
103
  }
104

105
  // Read distribution for phi-coordinate
106
  if (check_for_node(node, "phi")) {
51✔
107
    pugi::xml_node node_dist = node.child("phi");
51✔
108
    phi_ = distribution_from_xml(node_dist);
51✔
109
  } else {
110
    // If no distribution was specified, default to a single point at phi=0
111
    double x[] {0.0};
×
UNCOV
112
    double p[] {1.0};
×
UNCOV
113
    phi_ = make_unique<Discrete>(x, p, 1);
×
114
  }
115

116
  // Read distribution for z-coordinate
117
  if (check_for_node(node, "z")) {
51✔
118
    pugi::xml_node node_dist = node.child("z");
51✔
119
    z_ = distribution_from_xml(node_dist);
51✔
120
  } else {
121
    // If no distribution was specified, default to a single point at z=0
122
    double x[] {0.0};
×
UNCOV
123
    double p[] {1.0};
×
UNCOV
124
    z_ = make_unique<Discrete>(x, p, 1);
×
125
  }
126

127
  // Read cylinder center coordinates
128
  if (check_for_node(node, "origin")) {
51✔
129
    auto origin = get_node_array<double>(node, "origin");
51✔
130
    if (origin.size() == 3) {
51✔
131
      origin_ = origin;
51✔
132
    } else {
UNCOV
133
      fatal_error(
×
134
        "Origin for cylindrical source distribution must be length 3");
135
    }
136
  } else {
51✔
137
    // If no coordinates were specified, default to (0, 0, 0)
UNCOV
138
    origin_ = {0.0, 0.0, 0.0};
×
139
  }
140
}
51✔
141

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

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

156
SphericalIndependent::SphericalIndependent(pugi::xml_node node)
65✔
157
{
158
  // Read distribution for r-coordinate
159
  if (check_for_node(node, "r")) {
65✔
160
    pugi::xml_node node_dist = node.child("r");
65✔
161
    r_ = distribution_from_xml(node_dist);
65✔
162
  } else {
163
    // If no distribution was specified, default to a single point at r=0
164
    double x[] {0.0};
×
UNCOV
165
    double p[] {1.0};
×
UNCOV
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")) {
65✔
171
    pugi::xml_node node_dist = node.child("cos_theta");
65✔
172
    cos_theta_ = distribution_from_xml(node_dist);
65✔
173
  } else {
174
    // If no distribution was specified, default to a single point at
175
    // cos_theta=0
176
    double x[] {0.0};
×
UNCOV
177
    double p[] {1.0};
×
UNCOV
178
    cos_theta_ = make_unique<Discrete>(x, p, 1);
×
179
  }
180

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

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

206
Position SphericalIndependent::sample(uint64_t* seed) const
182,496✔
207
{
208
  double r = r_->sample(seed);
182,496✔
209
  double cos_theta = cos_theta_->sample(seed);
182,496✔
210
  double phi = phi_->sample(seed);
182,496✔
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,496✔
213
  double y = r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y;
182,496✔
214
  double z = r * cos_theta + origin_.z;
182,496✔
215
  return {x, y, z};
182,496✔
216
}
217

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

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

225
  if (get_node_value(node, "type", true, true) != "mesh") {
7✔
UNCOV
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"));
7✔
233
  // Get pointer to spatial distribution
234
  mesh_idx_ = model::mesh_map.at(mesh_id);
7✔
235

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

238
  check_element_types();
7✔
239

240
  size_t n_bins = this->n_sources();
7✔
241
  std::vector<double> strengths(n_bins, 1.0);
7✔
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")) {
7✔
247
    strengths = get_node_array<double>(node, "strengths");
4✔
248
    if (strengths.size() != n_bins) {
4✔
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")) {
6✔
257
    for (int i = 0; i < n_bins; i++) {
36,003✔
258
      strengths[i] *= this->mesh()->volume(i);
36,000✔
259
    }
260
  }
261

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

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

272
void MeshSpatial::check_element_types() const
226✔
273
{
274
  const auto umesh_ptr = dynamic_cast<const UnstructuredMesh*>(this->mesh());
226✔
275
  if (umesh_ptr) {
226✔
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✔
UNCOV
279
        fatal_error(
×
280
          "Mesh specified for source must contain only linear tetrahedra.");
281
      }
282
    }
283
  }
284
}
226✔
285

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

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

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

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

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

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

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

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

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

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

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

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

350
SpatialBox::SpatialBox(pugi::xml_node node, bool fission)
3,427✔
351
  : only_fissionable_ {fission}
3,427✔
352
{
353
  // Read lower-right/upper-left coordinates
354
  auto params = get_node_array<double>(node, "parameters");
3,427✔
355
  if (params.size() != 6)
3,427✔
UNCOV
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,427✔
360
  upper_right_ = Position {params[3], params[4], params[5]};
3,427✔
361
}
3,427✔
362

363
Position SpatialBox::sample(uint64_t* seed) const
3,857,664✔
364
{
365
  Position xi {prn(seed), prn(seed), prn(seed)};
3,857,664✔
366
  return lower_left_ + xi * (upper_right_ - lower_left_);
7,715,328✔
367
}
368

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

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

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

385
Position SpatialPoint::sample(uint64_t* seed) const
12,190,025✔
386
{
387
  return r_;
12,190,025✔
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