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

openmc-dev / openmc / 25588358527

09 May 2026 01:53AM UTC coverage: 81.388% (+0.1%) from 81.266%
25588358527

push

github

web-flow
Implement DecaySpectrum distribution type and utilize in R2S (#3930)

Co-authored-by: Copilot <copilot@github.com>

17751 of 25650 branches covered (69.2%)

Branch coverage included in aggregate %.

235 of 276 new or added lines in 7 files covered. (85.14%)

1 existing line in 1 file now uncovered.

58747 of 68342 relevant lines covered (85.96%)

47204330.4 hits per line

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

81.15
/src/source.cpp
1
#include "openmc/source.h"
2

3
#if defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
4
#define HAS_DYNAMIC_LINKING
5
#endif
6

7
#include <utility> // for move
8

9
#ifdef HAS_DYNAMIC_LINKING
10
#include <dlfcn.h> // for dlopen, dlsym, dlclose, dlerror
11
#endif
12

13
#include "openmc/tensor.h"
14
#include <fmt/core.h>
15

16
#include "openmc/bank.h"
17
#include "openmc/capi.h"
18
#include "openmc/cell.h"
19
#include "openmc/container_util.h"
20
#include "openmc/error.h"
21
#include "openmc/file_utils.h"
22
#include "openmc/geometry.h"
23
#include "openmc/hdf5_interface.h"
24
#include "openmc/material.h"
25
#include "openmc/mcpl_interface.h"
26
#include "openmc/memory.h"
27
#include "openmc/message_passing.h"
28
#include "openmc/mgxs_interface.h"
29
#include "openmc/nuclide.h"
30
#include "openmc/random_lcg.h"
31
#include "openmc/search.h"
32
#include "openmc/settings.h"
33
#include "openmc/simulation.h"
34
#include "openmc/state_point.h"
35
#include "openmc/string_utils.h"
36
#include "openmc/xml_interface.h"
37

38
namespace openmc {
39

40
std::atomic<int64_t> source_n_accept {0};
41
std::atomic<int64_t> source_n_reject {0};
42

43
namespace {
44

45
void validate_particle_type(ParticleType type, const std::string& context)
78,885✔
46
{
47
  if (type.is_transportable())
78,885!
48
    return;
78,885✔
49

50
  fatal_error(
×
51
    fmt::format("Unsupported source particle type '{}' (PDG {}) in {}.",
×
52
      type.str(), type.pdg_number(), context));
×
53
}
54

55
} // namespace
56

57
//==============================================================================
58
// Global variables
59
//==============================================================================
60

61
namespace model {
62

63
vector<unique_ptr<Source>> external_sources;
64

65
vector<unique_ptr<Source>> adjoint_sources;
66

67
DiscreteIndex external_sources_probability;
68

69
} // namespace model
70

71
//==============================================================================
72
// Source implementation
73
//==============================================================================
74

75
Source::Source(pugi::xml_node node)
4,349✔
76
{
77
  // Check for source strength
78
  if (check_for_node(node, "strength")) {
4,349✔
79
    strength_ = std::stod(get_node_value(node, "strength"));
8,182✔
80
    if (strength_ < 0.0) {
4,091!
81
      fatal_error("Source strength is negative.");
×
82
    }
83
  }
84

85
  // Check for additional defined constraints
86
  read_constraints(node);
4,349✔
87
}
4,349✔
88

89
unique_ptr<Source> Source::create(pugi::xml_node node)
4,349✔
90
{
91
  // if the source type is present, use it to determine the type
92
  // of object to create
93
  if (check_for_node(node, "type")) {
4,349✔
94
    std::string source_type = get_node_value(node, "type");
3,998✔
95
    if (source_type == "independent") {
3,998✔
96
      return make_unique<IndependentSource>(node);
3,881✔
97
    } else if (source_type == "file") {
117✔
98
      return make_unique<FileSource>(node);
15✔
99
    } else if (source_type == "compiled") {
102✔
100
      return make_unique<CompiledSourceWrapper>(node);
12✔
101
    } else if (source_type == "mesh") {
90!
102
      return make_unique<MeshSource>(node);
90✔
103
    } else {
104
      fatal_error(fmt::format("Invalid source type '{}' found.", source_type));
×
105
    }
106
  } else {
3,993✔
107
    // support legacy source format
108
    if (check_for_node(node, "file")) {
351✔
109
      return make_unique<FileSource>(node);
12✔
110
    } else if (check_for_node(node, "library")) {
339!
111
      return make_unique<CompiledSourceWrapper>(node);
×
112
    } else {
113
      return make_unique<IndependentSource>(node);
339✔
114
    }
115
  }
116
}
117

118
void Source::read_constraints(pugi::xml_node node)
4,349✔
119
{
120
  // Check for constraints node. For backwards compatibility, if no constraints
121
  // node is given, still try searching for domain constraints from top-level
122
  // node.
123
  pugi::xml_node constraints_node = node.child("constraints");
4,349✔
124
  if (constraints_node) {
4,349✔
125
    node = constraints_node;
844✔
126
  }
127

128
  // Check for domains to reject from
129
  if (check_for_node(node, "domain_type")) {
4,349✔
130
    std::string domain_type = get_node_value(node, "domain_type");
222✔
131
    if (domain_type == "cell") {
222✔
132
      domain_type_ = DomainType::CELL;
45✔
133
    } else if (domain_type == "material") {
177✔
134
      domain_type_ = DomainType::MATERIAL;
27✔
135
    } else if (domain_type == "universe") {
150!
136
      domain_type_ = DomainType::UNIVERSE;
150✔
137
    } else {
138
      fatal_error(
×
139
        std::string("Unrecognized domain type for constraint: " + domain_type));
×
140
    }
141

142
    auto ids = get_node_array<int>(node, "domain_ids");
222✔
143
    domain_ids_.insert(ids.begin(), ids.end());
222✔
144
  }
222✔
145

146
  if (check_for_node(node, "time_bounds")) {
4,349✔
147
    auto ids = get_node_array<double>(node, "time_bounds");
5✔
148
    if (ids.size() != 2) {
5!
149
      fatal_error("Time bounds must be represented by two numbers.");
×
150
    }
151
    time_bounds_ = std::make_pair(ids[0], ids[1]);
5✔
152
  }
5✔
153
  if (check_for_node(node, "energy_bounds")) {
4,349✔
154
    auto ids = get_node_array<double>(node, "energy_bounds");
5✔
155
    if (ids.size() != 2) {
5!
156
      fatal_error("Energy bounds must be represented by two numbers.");
×
157
    }
158
    energy_bounds_ = std::make_pair(ids[0], ids[1]);
5✔
159
  }
5✔
160

161
  if (check_for_node(node, "fissionable")) {
4,349✔
162
    only_fissionable_ = get_node_value_bool(node, "fissionable");
617✔
163
  }
164

165
  // Check for how to handle rejected particles
166
  if (check_for_node(node, "rejection_strategy")) {
4,349!
167
    std::string rejection_strategy = get_node_value(node, "rejection_strategy");
×
168
    if (rejection_strategy == "kill") {
×
169
      rejection_strategy_ = RejectionStrategy::KILL;
×
170
    } else if (rejection_strategy == "resample") {
×
171
      rejection_strategy_ = RejectionStrategy::RESAMPLE;
×
172
    } else {
173
      fatal_error(std::string(
×
174
        "Unrecognized strategy source rejection: " + rejection_strategy));
175
    }
176
  }
×
177
}
4,349✔
178

179
void check_rejection_fraction(int64_t n_reject, int64_t n_accept)
15,567,505✔
180
{
181
  // Don't check unless we've hit a minimum number of total sites rejected
182
  if (n_reject < EXTSRC_REJECT_THRESHOLD)
15,567,505✔
183
    return;
184

185
  // Compute fraction of accepted sites and compare against minimum
186
  double fraction = static_cast<double>(n_accept) / n_reject;
607,094✔
187
  if (fraction <= settings::source_rejection_fraction) {
607,094✔
188
    fatal_error(fmt::format(
4✔
189
      "Too few source sites satisfied the constraints (minimum source "
190
      "rejection fraction = {}). Please check your source definition or "
191
      "set a lower value of Settings.source_rejection_fraction.",
192
      settings::source_rejection_fraction));
193
  }
194
}
195

196
SourceSite Source::sample_with_constraints(uint64_t* seed) const
15,567,505✔
197
{
198
  bool accepted = false;
15,567,505✔
199
  int64_t n_local_reject = 0;
15,567,505✔
200
  SourceSite site {};
15,567,505✔
201

202
  while (!accepted) {
47,319,313✔
203
    // Sample a source site without considering constraints yet
204
    site = this->sample(seed);
16,184,303✔
205

206
    if (constraints_applied()) {
16,184,303✔
207
      accepted = true;
208
    } else {
209
      // Check whether sampled site satisfies constraints
210
      accepted = satisfies_spatial_constraints(site.r) &&
17,731,563✔
211
                 satisfies_energy_constraints(site.E) &&
1,232,131✔
212
                 satisfies_time_constraints(site.time);
310,183✔
213
      if (!accepted) {
616,798✔
214
        ++n_local_reject;
616,798✔
215

216
        // Check per-particle rejection limit
217
        if (n_local_reject >= MAX_SOURCE_REJECTIONS_PER_SAMPLE) {
616,798!
218
          fatal_error("Exceeded maximum number of source rejections per "
×
219
                      "sample. Please check your source definition.");
220
        }
221

222
        // For the "kill" strategy, accept particle but set weight to 0 so that
223
        // it is terminated immediately
224
        if (rejection_strategy_ == RejectionStrategy::KILL) {
616,798!
225
          accepted = true;
×
226
          site.wgt = 0.0;
×
227
        }
228
      }
229
    }
230
  }
231

232
  // Flush local rejection count, update accept counter, and check overall
233
  // rejection fraction
234
  if (n_local_reject > 0) {
15,567,505✔
235
    source_n_reject += n_local_reject;
8,691✔
236
  }
237
  ++source_n_accept;
15,567,505✔
238
  check_rejection_fraction(source_n_reject, source_n_accept);
15,567,505✔
239

240
  return site;
15,567,501✔
241
}
242

243
bool Source::satisfies_energy_constraints(double E) const
15,582,517✔
244
{
245
  return E > energy_bounds_.first && E < energy_bounds_.second;
15,582,517!
246
}
247

248
bool Source::satisfies_time_constraints(double time) const
310,183✔
249
{
250
  return time > time_bounds_.first && time < time_bounds_.second;
310,183✔
251
}
252

253
bool Source::satisfies_spatial_constraints(Position r) const
17,899,014✔
254
{
255
  GeometryState geom_state;
17,899,014✔
256
  geom_state.r() = r;
17,899,014✔
257
  geom_state.u() = {0.0, 0.0, 1.0};
17,899,014✔
258

259
  // Reject particle if it's not in the geometry at all
260
  bool found = exhaustive_find_cell(geom_state);
17,899,014✔
261
  if (!found)
17,899,014✔
262
    return false;
263

264
  // Check the geometry state against specified domains
265
  bool accepted = true;
17,676,279✔
266
  if (!domain_ids_.empty()) {
17,676,279✔
267
    if (domain_type_ == DomainType::MATERIAL) {
1,091,051✔
268
      auto mat_index = geom_state.material();
100,000✔
269
      if (mat_index == MATERIAL_VOID) {
100,000!
270
        accepted = false;
271
      } else {
272
        accepted = contains(domain_ids_, model::materials[mat_index]->id());
200,000✔
273
      }
274
    } else {
275
      for (int i = 0; i < geom_state.n_coord(); i++) {
1,912,538✔
276
        auto id =
991,051✔
277
          (domain_type_ == DomainType::CELL)
278
            ? model::cells[geom_state.coord(i).cell()].get()->id_
991,051!
279
            : model::universes[geom_state.coord(i).universe()].get()->id_;
×
280
        if ((accepted = contains(domain_ids_, id)))
1,982,102✔
281
          break;
282
      }
283
    }
284
  }
285

286
  // Check if spatial site is in fissionable material
287
  if (accepted && only_fissionable_) {
17,676,279✔
288
    // Determine material
289
    auto mat_index = geom_state.material();
489,990✔
290
    if (mat_index == MATERIAL_VOID) {
489,990!
291
      accepted = false;
292
    } else {
293
      accepted = model::materials[mat_index]->fissionable();
489,990✔
294
    }
295
  }
296

297
  return accepted;
298
}
17,899,014✔
299

300
//==============================================================================
301
// IndependentSource implementation
302
//==============================================================================
303

304
IndependentSource::IndependentSource(
939✔
305
  UPtrSpace space, UPtrAngle angle, UPtrDist energy, UPtrDist time)
939✔
306
  : space_ {std::move(space)}, angle_ {std::move(angle)},
939✔
307
    energy_ {std::move(energy)}, time_ {std::move(time)}
939✔
308
{}
939✔
309

310
IndependentSource::IndependentSource(pugi::xml_node node) : Source(node)
4,220✔
311
{
312
  // Check for particle type
313
  if (check_for_node(node, "particle")) {
4,220✔
314
    auto temp_str = get_node_value(node, "particle", false, true);
3,881✔
315
    particle_ = ParticleType(temp_str);
3,881✔
316
    if (particle_ == ParticleType::photon() ||
3,881✔
317
        particle_ == ParticleType::electron() ||
3,881✔
318
        particle_ == ParticleType::positron()) {
3,791!
319
      settings::photon_transport = true;
90✔
320
    }
321
  }
3,881✔
322
  validate_particle_type(particle_, "IndependentSource");
4,220✔
323

324
  // Check for external source file
325
  if (check_for_node(node, "file")) {
4,220!
326

327
  } else {
328

329
    // Spatial distribution for external source
330
    if (check_for_node(node, "space")) {
4,220✔
331
      space_ = SpatialDistribution::create(node.child("space"));
3,248✔
332
    } else {
333
      // If no spatial distribution specified, make it a point source
334
      space_ = UPtrSpace {new SpatialPoint()};
972✔
335
    }
336

337
    // For backwards compatibility, check for only fissionable setting on box
338
    // source
339
    auto space_box = dynamic_cast<SpatialBox*>(space_.get());
4,220!
340
    if (space_box) {
4,220✔
341
      if (!only_fissionable_) {
1,709✔
342
        only_fissionable_ = space_box->only_fissionable();
1,092✔
343
      }
344
    }
345

346
    // Determine external source angular distribution
347
    if (check_for_node(node, "angle")) {
4,220✔
348
      angle_ = UnitSphereDistribution::create(node.child("angle"));
1,534✔
349
    } else {
350
      angle_ = UPtrAngle {new Isotropic()};
2,686✔
351
    }
352

353
    // Determine external source energy distribution
354
    if (check_for_node(node, "energy")) {
4,220✔
355
      pugi::xml_node node_dist = node.child("energy");
2,147✔
356
      energy_ = distribution_from_xml(node_dist);
2,147✔
357

358
      // For decay photon sources, use the absolute photon emission rate in
359
      // [photons/s] as the source strength
360
      if (dynamic_cast<DecaySpectrum*>(energy_.get())) {
2,147!
361
        if (strength_ != 1.0) {
25!
NEW
362
          warning(fmt::format(
×
363
            "Source strength of {} is ignored because the source uses a "
364
            "DecaySpectrum energy distribution. The source strength will be "
365
            "set from the DecaySpectrum emission rate.",
NEW
366
            strength_));
×
367
        }
368
        strength_ = energy_->integral();
25✔
369
      }
370
    } else {
371
      // Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
372
      energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)};
2,073✔
373
    }
374

375
    // Determine external source time distribution
376
    if (check_for_node(node, "time")) {
4,220✔
377
      pugi::xml_node node_dist = node.child("time");
17✔
378
      time_ = distribution_from_xml(node_dist);
17✔
379
    } else {
380
      // Default to a Constant time T=0
381
      double T[] {0.0};
4,203✔
382
      double p[] {1.0};
4,203✔
383
      time_ = UPtrDist {new Discrete {T, p, 1}};
4,203✔
384
    }
385
  }
386
}
4,220✔
387

388
SourceSite IndependentSource::sample(uint64_t* seed) const
16,295,355✔
389
{
390
  SourceSite site {};
16,295,355✔
391
  site.particle = particle_;
16,295,355✔
392
  double r_wgt = 1.0;
16,295,355✔
393
  double E_wgt = 1.0;
16,295,355✔
394

395
  // Repeat sampling source location until a good site has been accepted
396
  bool accepted = false;
16,295,355✔
397
  int64_t n_local_reject = 0;
16,295,355✔
398

399
  while (!accepted) {
33,272,421✔
400

401
    // Sample spatial distribution
402
    auto [r, r_wgt_temp] = space_->sample(seed);
16,977,066✔
403
    site.r = r;
16,977,066✔
404
    r_wgt = r_wgt_temp;
16,977,066✔
405

406
    // Check if sampled position satisfies spatial constraints
407
    accepted = satisfies_spatial_constraints(site.r);
16,977,066✔
408

409
    // Check for rejection
410
    if (!accepted) {
16,977,066✔
411
      ++n_local_reject;
681,711✔
412
      if (n_local_reject >= MAX_SOURCE_REJECTIONS_PER_SAMPLE) {
681,711!
413
        fatal_error("Exceeded maximum number of source rejections per "
×
414
                    "sample. Please check your source definition.");
415
      }
416
    }
417
  }
418

419
  // Sample angle
420
  auto [u, u_wgt] = angle_->sample(seed);
16,295,355✔
421
  site.u = u;
16,295,355✔
422

423
  site.wgt = r_wgt * u_wgt;
16,295,355✔
424

425
  // Sample energy and time for neutron and photon sources
426
  if (settings::solver_type != SolverType::RANDOM_RAY) {
16,295,355✔
427
    // Check for monoenergetic source above maximum particle energy
428
    auto p = particle_.transport_index();
15,262,355✔
429
    auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());
15,262,355!
430
    auto decay_spectrum = dynamic_cast<DecaySpectrum*>(energy_.get());
15,262,355!
431
    if (energy_ptr) {
15,262,355✔
432
      auto energies =
7,935,455✔
433
        tensor::Tensor<double>(energy_ptr->x().data(), energy_ptr->x().size());
7,935,455✔
434
      if ((energies > data::energy_max[p]).any()) {
23,806,365!
435
        fatal_error("Source energy above range of energies of at least "
×
436
                    "one cross section table");
437
      }
438
    }
7,935,455✔
439

440
    while (true) {
15,262,355✔
441
      // Sample energy spectrum. For decay photon sources, also get the parent
442
      // nuclide index to store in the source site for tallying purposes.
443
      if (decay_spectrum) {
15,262,355✔
444
        auto sample = decay_spectrum->sample_with_parent(seed);
175,000✔
445
        site.E = sample.energy;
175,000✔
446
        E_wgt = sample.weight;
175,000✔
447
        site.parent_nuclide = sample.parent_nuclide;
175,000✔
448
      } else {
449
        auto [E, E_wgt_temp] = energy_->sample(seed);
15,087,355✔
450
        site.E = E;
15,087,355✔
451
        E_wgt = E_wgt_temp;
15,087,355✔
452
      }
453

454
      // Resample if energy falls above maximum particle energy
455
      if (site.E < data::energy_max[p] &&
30,524,710!
456
          (satisfies_energy_constraints(site.E)))
15,262,355✔
457
        break;
458

459
      ++n_local_reject;
×
460
      if (n_local_reject >= MAX_SOURCE_REJECTIONS_PER_SAMPLE) {
×
461
        fatal_error("Exceeded maximum number of source rejections per "
×
462
                    "sample. Please check your source definition.");
463
      }
464
    }
465

466
    // Sample particle creation time
467
    auto [time, time_wgt] = time_->sample(seed);
15,262,355✔
468
    site.time = time;
15,262,355✔
469

470
    site.wgt *= (E_wgt * time_wgt);
15,262,355✔
471
  }
472

473
  // Flush local rejection count into global counter
474
  if (n_local_reject > 0) {
16,295,355✔
475
    source_n_reject += n_local_reject;
160,195✔
476
  }
477

478
  return site;
16,295,355✔
479
}
480

481
//==============================================================================
482
// FileSource implementation
483
//==============================================================================
484

485
FileSource::FileSource(pugi::xml_node node) : Source(node)
27✔
486
{
487
  auto path = get_node_value(node, "file", false, true);
27✔
488
  load_sites_from_file(path);
27✔
489
}
22✔
490

491
FileSource::FileSource(const std::string& path)
12✔
492
{
493
  load_sites_from_file(path);
12✔
494
}
12✔
495

496
void FileSource::load_sites_from_file(const std::string& path)
39✔
497
{
498
  // If MCPL file, use the dedicated file reader
499
  if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) {
66!
500
    sites_ = mcpl_source_sites(path);
12✔
501
  } else {
502
    // Check if source file exists
503
    if (!file_exists(path)) {
27!
504
      fatal_error(fmt::format("Source file '{}' does not exist.", path));
×
505
    }
506

507
    write_message(6, "Reading source file from {}...", path);
27✔
508

509
    // Open the binary file
510
    hid_t file_id = file_open(path, 'r', true);
27✔
511

512
    // Check to make sure this is a source file
513
    std::string filetype;
27✔
514
    read_attribute(file_id, "filetype", filetype);
27✔
515
    if (filetype != "source" && filetype != "statepoint") {
27!
516
      fatal_error("Specified starting source file not a source file type.");
×
517
    }
518

519
    // Read in the source particles
520
    read_source_bank(file_id, sites_, false);
27✔
521

522
    // Close file
523
    file_close(file_id);
22✔
524
  }
22✔
525

526
  // Make sure particles in source file have valid types
527
  for (const auto& site : this->sites_) {
74,039✔
528
    validate_particle_type(site.particle, "FileSource");
148,010✔
529
  }
530
}
34✔
531

532
SourceSite FileSource::sample(uint64_t* seed) const
130,162✔
533
{
534
  // Sample a particle randomly from list
535
  size_t i_site = sites_.size() * prn(seed);
130,162✔
536
  return sites_[i_site];
130,162✔
537
}
538

539
//==============================================================================
540
// CompiledSourceWrapper implementation
541
//==============================================================================
542

543
CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) : Source(node)
12✔
544
{
545
  // Get shared library path and parameters
546
  auto path = get_node_value(node, "library", false, true);
12✔
547
  std::string parameters;
12✔
548
  if (check_for_node(node, "parameters")) {
12✔
549
    parameters = get_node_value(node, "parameters", false, true);
6✔
550
  }
551
  setup(path, parameters);
12✔
552
}
12✔
553

554
void CompiledSourceWrapper::setup(
12✔
555
  const std::string& path, const std::string& parameters)
556
{
557
#ifdef HAS_DYNAMIC_LINKING
558
  // Open the library
559
  shared_library_ = dlopen(path.c_str(), RTLD_LAZY);
12✔
560
  if (!shared_library_) {
12!
561
    fatal_error("Couldn't open source library " + path);
×
562
  }
563

564
  // reset errors
565
  dlerror();
12✔
566

567
  // get the function to create the custom source from the library
568
  auto create_compiled_source = reinterpret_cast<create_compiled_source_t*>(
12✔
569
    dlsym(shared_library_, "openmc_create_source"));
12✔
570

571
  // check for any dlsym errors
572
  auto dlsym_error = dlerror();
12✔
573
  if (dlsym_error) {
12!
574
    std::string error_msg = fmt::format(
×
575
      "Couldn't open the openmc_create_source symbol: {}", dlsym_error);
×
576
    dlclose(shared_library_);
×
577
    fatal_error(error_msg);
×
578
  }
×
579

580
  // create a pointer to an instance of the custom source
581
  compiled_source_ = create_compiled_source(parameters);
12✔
582

583
#else
584
  fatal_error("Custom source libraries have not yet been implemented for "
585
              "non-POSIX systems");
586
#endif
587
}
12✔
588

589
CompiledSourceWrapper::~CompiledSourceWrapper()
24✔
590
{
591
  // Make sure custom source is cleared before closing shared library
592
  if (compiled_source_.get())
12!
593
    compiled_source_.reset();
12✔
594

595
#ifdef HAS_DYNAMIC_LINKING
596
  dlclose(shared_library_);
12✔
597
#else
598
  fatal_error("Custom source libraries have not yet been implemented for "
599
              "non-POSIX systems");
600
#endif
601
}
24✔
602

603
//==============================================================================
604
// MeshElementSpatial implementation
605
//==============================================================================
606

607
std::pair<Position, double> MeshElementSpatial::sample(uint64_t* seed) const
697,023✔
608
{
609
  return {model::meshes[mesh_index_]->sample_element(elem_index_, seed), 1.0};
697,023✔
610
}
611

612
//==============================================================================
613
// MeshSource implementation
614
//==============================================================================
615

616
MeshSource::MeshSource(pugi::xml_node node) : Source(node)
90✔
617
{
618
  int32_t mesh_id = stoi(get_node_value(node, "mesh"));
180✔
619
  int32_t mesh_idx = model::mesh_map.at(mesh_id);
90✔
620
  const auto& mesh = model::meshes[mesh_idx];
90✔
621

622
  std::vector<double> strengths;
90✔
623
  // read all source distributions and populate strengths vector for MeshSpatial
624
  // object
625
  for (auto source_node : node.children("source")) {
750✔
626
    auto src = Source::create(source_node);
660✔
627
    if (auto ptr = dynamic_cast<IndependentSource*>(src.get())) {
660!
628
      src.release();
660✔
629
      sources_.emplace_back(ptr);
660✔
630
    } else {
631
      fatal_error(
×
632
        "The source assigned to each element must be an IndependentSource.");
633
    }
634
    strengths.push_back(sources_.back()->strength());
660✔
635
  }
660✔
636

637
  // Set spatial distributions for each mesh element
638
  for (int elem_index = 0; elem_index < sources_.size(); ++elem_index) {
750✔
639
    sources_[elem_index]->set_space(
660✔
640
      std::make_unique<MeshElementSpatial>(mesh_idx, elem_index));
1,320✔
641
  }
642

643
  // Make sure sources use valid particle types
644
  for (const auto& src : sources_) {
750✔
645
    validate_particle_type(src->particle_type(), "MeshSource");
1,320✔
646
  }
647

648
  // the number of source distributions should either be one or equal to the
649
  // number of mesh elements
650
  if (sources_.size() > 1 && sources_.size() != mesh->n_bins()) {
90!
651
    fatal_error(fmt::format("Incorrect number of source distributions ({}) for "
×
652
                            "mesh source with {} elements.",
653
      sources_.size(), mesh->n_bins()));
×
654
  }
655

656
  space_ = std::make_unique<MeshSpatial>(mesh_idx, strengths);
90✔
657
}
90✔
658

659
SourceSite MeshSource::sample(uint64_t* seed) const
691,786✔
660
{
661
  // Sample a mesh element based on the relative strengths
662
  int32_t element = space_->sample_element_index(seed);
691,786✔
663

664
  // Sample the distribution for the specific mesh element; note that the
665
  // spatial distribution has been set for each element using MeshElementSpatial
666
  return source(element)->sample_with_constraints(seed);
1,383,572!
667
}
668

669
//==============================================================================
670
// Non-member functions
671
//==============================================================================
672

673
void initialize_source()
1,715✔
674
{
675
  write_message("Initializing source particles...", 5);
1,715✔
676

677
// Generation source sites from specified distribution in user input
678
#pragma omp parallel for
679
  for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
1,290,380✔
680
    // initialize random number seed
681
    int64_t id = simulation::total_gen * settings::n_particles +
1,288,665✔
682
                 simulation::work_index[mpi::rank] + i + 1;
1,288,665✔
683
    uint64_t seed = init_seed(id, STREAM_SOURCE);
1,288,665✔
684

685
    // sample external source distribution
686
    simulation::source_bank[i] = sample_external_source(&seed);
1,288,665✔
687
  }
688

689
  // Write out initial source
690
  if (settings::write_initial_source) {
1,715!
691
    write_message("Writing out initial source...", 5);
×
692
    std::string filename = settings::path_output + "initial_source.h5";
×
693
    hid_t file_id = file_open(filename, 'w', true);
×
694
    write_source_bank(file_id, simulation::source_bank, simulation::work_index);
×
695
    file_close(file_id);
×
696
  }
×
697
}
1,715✔
698

699
SourceSite sample_external_source(uint64_t* seed)
14,875,719✔
700
{
701
  // Sample from among multiple source distributions
702
  int i = 0;
14,875,719✔
703
  int n_sources = model::external_sources.size();
14,875,719✔
704
  if (n_sources > 1) {
14,875,719✔
705
    if (settings::uniform_source_sampling) {
1,616,500✔
706
      i = prn(seed) * n_sources;
1,000✔
707
    } else {
708
      i = model::external_sources_probability.sample(seed);
1,615,500✔
709
    }
710
  }
711

712
  // Sample source site from i-th source distribution
713
  SourceSite site {model::external_sources[i]->sample_with_constraints(seed)};
14,875,719✔
714

715
  // For uniform source sampling, multiply the weight by the ratio of the actual
716
  // probability of sampling source i to the biased probability of sampling
717
  // source i, which is (strength_i / total_strength) / (1 / n)
718
  if (n_sources > 1 && settings::uniform_source_sampling) {
14,875,715✔
719
    double total_strength = model::external_sources_probability.integral();
1,000✔
720
    site.wgt *=
2,000✔
721
      model::external_sources[i]->strength() * n_sources / total_strength;
1,000✔
722
  }
723

724
  // If running in MG, convert site.E to group
725
  if (!settings::run_CE) {
14,875,715✔
726
    site.E = lower_bound_index(data::mg.rev_energy_bins_.begin(),
792,150✔
727
      data::mg.rev_energy_bins_.end(), site.E);
728
    site.E = data::mg.num_energy_groups_ - site.E - 1.;
792,150✔
729
  }
730

731
  return site;
14,875,715✔
732
}
733

734
void free_memory_source()
3,625✔
735
{
736
  model::external_sources.clear();
3,625✔
737
  model::adjoint_sources.clear();
3,625✔
738
  reset_source_rejection_counters();
3,625✔
739
}
3,625✔
740

741
void reset_source_rejection_counters()
6,813✔
742
{
743
  source_n_accept = 0;
6,813✔
744
  source_n_reject = 0;
6,813✔
745
}
6,813✔
746

747
//==============================================================================
748
// C API
749
//==============================================================================
750

751
extern "C" int openmc_sample_external_source(
165✔
752
  size_t n, uint64_t* seed, void* sites)
753
{
754
  if (!sites || !seed) {
165!
755
    set_errmsg("Received null pointer.");
×
756
    return OPENMC_E_INVALID_ARGUMENT;
×
757
  }
758

759
  if (model::external_sources.empty()) {
165!
760
    set_errmsg("No external sources have been defined.");
×
761
    return OPENMC_E_OUT_OF_BOUNDS;
×
762
  }
763

764
  auto sites_array = static_cast<SourceSite*>(sites);
165✔
765

766
  // Derive independent per-particle seeds from the base seed so that
767
  // each iteration has its own RNG state for thread-safe parallel sampling.
768
  uint64_t base_seed = *seed;
165✔
769

770
#pragma omp parallel for schedule(static)
771
  for (size_t i = 0; i < n; ++i) {
1,071,485✔
772
    uint64_t particle_seed = init_seed(base_seed + i, STREAM_SOURCE);
1,071,320✔
773
    sites_array[i] = sample_external_source(&particle_seed);
1,071,320✔
774
  }
775
  return 0;
776
}
777

778
} // 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

© 2026 Coveralls, Inc