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

openmc-dev / openmc / 20733102133

06 Jan 2026 12:03AM UTC coverage: 81.932% (-0.2%) from 82.174%
20733102133

Pull #3702

github

web-flow
Merge b36a045c7 into 60ddafa9b
Pull Request #3702: Random Ray Kinetic Simulation Mode

17551 of 24380 branches covered (71.99%)

Branch coverage included in aggregate %.

869 of 1117 new or added lines in 21 files covered. (77.8%)

483 existing lines in 19 files now uncovered.

55915 of 65287 relevant lines covered (85.64%)

50695599.69 hits per line

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

75.06
/src/settings.cpp
1
#include "openmc/settings.h"
2
#include "openmc/random_ray/flat_source_domain.h"
3

4
#include <cmath>  // for ceil, pow
5
#include <limits> // for numeric_limits
6
#include <string>
7

8
#include <fmt/core.h>
9
#ifdef _OPENMP
10
#include <omp.h>
11
#endif
12

13
#include "openmc/capi.h"
14
#include "openmc/collision_track.h"
15
#include "openmc/constants.h"
16
#include "openmc/container_util.h"
17
#include "openmc/distribution.h"
18
#include "openmc/distribution_multi.h"
19
#include "openmc/distribution_spatial.h"
20
#include "openmc/eigenvalue.h"
21
#include "openmc/error.h"
22
#include "openmc/file_utils.h"
23
#include "openmc/mcpl_interface.h"
24
#include "openmc/mesh.h"
25
#include "openmc/message_passing.h"
26
#include "openmc/output.h"
27
#include "openmc/plot.h"
28
#include "openmc/random_lcg.h"
29
#include "openmc/random_ray/random_ray.h"
30
#include "openmc/reaction.h"
31
#include "openmc/simulation.h"
32
#include "openmc/source.h"
33
#include "openmc/string_utils.h"
34
#include "openmc/tallies/trigger.h"
35
#include "openmc/volume_calc.h"
36
#include "openmc/weight_windows.h"
37
#include "openmc/xml_interface.h"
38

39
namespace openmc {
40

41
//==============================================================================
42
// Global variables
43
//==============================================================================
44

45
namespace settings {
46

47
// Default values for boolean flags
48
bool assume_separate {false};
49
bool check_overlaps {false};
50
bool collision_track {false};
51
bool cmfd_run {false};
52
bool confidence_intervals {false};
53
bool create_delayed_neutrons {true};
54
bool create_fission_neutrons {true};
55
bool delayed_photon_scaling {true};
56
bool entropy_on {false};
57
bool event_based {false};
58
bool ifp_on {false};
59
bool kinetic_simulation {false};
60
bool legendre_to_tabular {true};
61
bool material_cell_offsets {true};
62
bool output_summary {true};
63
bool output_tallies {true};
64
bool particle_restart_run {false};
65
bool photon_transport {false};
66
bool reduce_tallies {true};
67
bool res_scat_on {false};
68
bool restart_run {false};
69
bool run_CE {true};
70
bool source_latest {false};
71
bool source_separate {false};
72
bool source_write {true};
73
bool source_mcpl_write {false};
74
bool surf_source_write {false};
75
bool surf_mcpl_write {false};
76
bool surf_source_read {false};
77
bool survival_biasing {false};
78
bool survival_normalization {false};
79
bool temperature_multipole {false};
80
bool trigger_on {false};
81
bool trigger_predict {false};
82
bool uniform_source_sampling {false};
83
bool ufs_on {false};
84
bool urr_ptables_on {true};
85
bool use_decay_photons {false};
86
bool weight_windows_on {false};
87
bool weight_window_checkpoint_surface {false};
88
bool weight_window_checkpoint_collision {true};
89
bool write_all_tracks {false};
90
bool write_initial_source {false};
91

92
std::string path_cross_sections;
93
std::string path_input;
94
std::string path_output;
95
std::string path_particle_restart;
96
std::string path_sourcepoint;
97
std::string path_statepoint;
98
const char* path_statepoint_c {path_statepoint.c_str()};
99
std::string weight_windows_file;
100

101
int32_t n_inactive {0};
102
int32_t max_lost_particles {10};
103
double rel_max_lost_particles {1.0e-6};
104
int32_t max_write_lost_particles {-1};
105
int32_t gen_per_batch {1};
106
int64_t n_particles {-1};
107

108
int64_t max_particles_in_flight {100000};
109
int max_particle_events {1000000};
110

111
ElectronTreatment electron_treatment {ElectronTreatment::TTB};
112
array<double, 4> energy_cutoff {0.0, 1000.0, 0.0, 0.0};
113
array<double, 4> time_cutoff {INFTY, INFTY, INFTY, INFTY};
114
int ifp_n_generation {-1};
115
IFPParameter ifp_parameter {IFPParameter::None};
116
int legendre_to_tabular_points {C_NONE};
117
int max_order {0};
118
int n_log_bins {8000};
119
int n_batches;
120
int n_max_batches;
121
int max_secondaries {10000};
122
int max_history_splits {10'000'000};
123
int max_tracks {1000};
124
ResScatMethod res_scat_method {ResScatMethod::rvs};
125
double res_scat_energy_min {0.01};
126
double res_scat_energy_max {1000.0};
127
vector<std::string> res_scat_nuclides;
128
RunMode run_mode {RunMode::UNSET};
129
SolverType solver_type {SolverType::MONTE_CARLO};
130
std::unordered_set<int> sourcepoint_batch;
131
std::unordered_set<int> statepoint_batch;
132
double source_rejection_fraction {0.05};
133
double free_gas_threshold {400.0};
134
std::unordered_set<int> source_write_surf_id;
135
CollisionTrackConfig collision_track_config {};
136
int64_t ssw_max_particles;
137
int64_t ssw_max_files;
138
int64_t ssw_cell_id {C_NONE};
139
SSWCellType ssw_cell_type {SSWCellType::None};
140
TemperatureMethod temperature_method {TemperatureMethod::NEAREST};
141
double temperature_tolerance {10.0};
142
double temperature_default {293.6};
143
array<double, 2> temperature_range {0.0, 0.0};
144
int trace_batch;
145
int trace_gen;
146
int64_t trace_particle;
147
vector<array<int, 3>> track_identifiers;
148
int trigger_batch_interval {1};
149
int verbosity {-1};
150
double weight_cutoff {0.25};
151
double weight_survive {1.0};
152

153
// Timestep variables for kinetic simulation
154
int n_timesteps;
155
double dt;
156

157
} // namespace settings
158

159
//==============================================================================
160
// Functions
161
//==============================================================================
162

163
void get_run_parameters(pugi::xml_node node_base)
7,226✔
164
{
165
  using namespace settings;
166
  using namespace pugi;
167

168
  // Check number of particles
169
  if (!check_for_node(node_base, "particles")) {
7,226!
170
    fatal_error("Need to specify number of particles.");
×
171
  }
172

173
  // Get number of particles if it wasn't specified as a command-line argument
174
  if (n_particles == -1) {
7,226!
175
    n_particles = std::stoll(get_node_value(node_base, "particles"));
7,226✔
176
  }
177

178
  // Get maximum number of in flight particles for event-based mode
179
  if (check_for_node(node_base, "max_particles_in_flight")) {
7,226!
180
    max_particles_in_flight =
×
181
      std::stoll(get_node_value(node_base, "max_particles_in_flight"));
×
182
  }
183

184
  // Get maximum number of events allowed per particle
185
  if (check_for_node(node_base, "max_particle_events")) {
7,226!
186
    max_particle_events =
×
187
      std::stoll(get_node_value(node_base, "max_particle_events"));
×
188
  }
189

190
  // Get number of basic batches
191
  if (check_for_node(node_base, "batches")) {
7,226!
192
    n_batches = std::stoi(get_node_value(node_base, "batches"));
7,226✔
193
  }
194
  if (!trigger_on)
7,226✔
195
    n_max_batches = n_batches;
7,080✔
196

197
  // Get max number of lost particles
198
  if (check_for_node(node_base, "max_lost_particles")) {
7,226✔
199
    max_lost_particles =
49✔
200
      std::stoi(get_node_value(node_base, "max_lost_particles"));
49✔
201
  }
202

203
  // Get relative number of lost particles
204
  if (check_for_node(node_base, "rel_max_lost_particles")) {
7,226!
205
    rel_max_lost_particles =
×
206
      std::stod(get_node_value(node_base, "rel_max_lost_particles"));
×
207
  }
208

209
  // Get relative number of lost particles
210
  if (check_for_node(node_base, "max_write_lost_particles")) {
7,226✔
211
    max_write_lost_particles =
16✔
212
      std::stoi(get_node_value(node_base, "max_write_lost_particles"));
16✔
213
  }
214

215
  // Get number of inactive batches
216
  if (run_mode == RunMode::EIGENVALUE ||
7,226✔
217
      solver_type == SolverType::RANDOM_RAY) {
2,688✔
218
    if (check_for_node(node_base, "inactive")) {
4,939✔
219
      n_inactive = std::stoi(get_node_value(node_base, "inactive"));
4,754✔
220
    }
221
    if (check_for_node(node_base, "generations_per_batch")) {
4,939✔
222
      gen_per_batch =
16✔
223
        std::stoi(get_node_value(node_base, "generations_per_batch"));
16✔
224
    }
225

226
    // Preallocate space for keff and entropy by generation
227
    int m = settings::n_max_batches * settings::gen_per_batch;
4,939✔
228
    simulation::k_generation.reserve(m);
4,939✔
229
    simulation::entropy.reserve(m);
4,939✔
230

231
    // Get the trigger information for keff
232
    if (check_for_node(node_base, "keff_trigger")) {
4,939✔
233
      xml_node node_keff_trigger = node_base.child("keff_trigger");
107✔
234

235
      if (check_for_node(node_keff_trigger, "type")) {
107!
236
        auto temp = get_node_value(node_keff_trigger, "type", true, true);
107✔
237
        if (temp == "std_dev") {
107!
238
          keff_trigger.metric = TriggerMetric::standard_deviation;
107✔
239
        } else if (temp == "variance") {
×
240
          keff_trigger.metric = TriggerMetric::variance;
×
241
        } else if (temp == "rel_err") {
×
242
          keff_trigger.metric = TriggerMetric::relative_error;
×
243
        } else {
244
          fatal_error("Unrecognized keff trigger type " + temp);
×
245
        }
246
      } else {
107✔
247
        fatal_error("Specify keff trigger type in settings XML");
×
248
      }
249

250
      if (check_for_node(node_keff_trigger, "threshold")) {
107!
251
        keff_trigger.threshold =
107✔
252
          std::stod(get_node_value(node_keff_trigger, "threshold"));
107✔
253
        if (keff_trigger.threshold <= 0) {
107!
254
          fatal_error("keff trigger threshold must be positive");
×
255
        }
256
      } else {
257
        fatal_error("Specify keff trigger threshold in settings XML");
×
258
      }
259
    }
260
  }
261

262
  // Kinetic variables
263
  if (check_for_node(node_base, "kinetic_simulation")) {
7,226✔
264
    kinetic_simulation = get_node_value_bool(node_base, "kinetic_simulation");
112✔
265
    if (solver_type != SolverType::RANDOM_RAY) {
112!
NEW
266
      fatal_error("Unsupported solver selected for kinetic simulation. Kinetic "
×
267
                  "simulations currently only support the random ray solver.");
268
    }
269
    if (run_mode != RunMode::EIGENVALUE) {
112!
NEW
270
      fatal_error(
×
271
        "Unsupported run mode selected for kinetic simulation. Kinetic "
272
        "simulations currently only support run mode based on an eigenvalue "
273
        "simulation establishing an initial condition.");
274
    }
275
  }
276

277
  // Get timestep parameters for kinetic simulations
278
  if (kinetic_simulation) {
7,226✔
279
    xml_node ts_node = node_base.child("timestep_parameters");
112✔
280
    if (check_for_node(ts_node, "n_timesteps")) {
112!
281
      n_timesteps = std::stoi(get_node_value(ts_node, "n_timesteps"));
112✔
282
    } else {
NEW
283
      fatal_error("Specify number of timesteps in settings XML");
×
284
    }
285
    if (check_for_node(ts_node, "timestep_units")) {
112!
286
      std::string units = get_node_value(ts_node, "timestep_units");
112✔
287
      if (check_for_node(ts_node, "dt")) {
112!
288
        dt = std::stod(get_node_value(ts_node, "dt"));
112✔
289
        double factor_to_seconds;
290
        if (units == "ms") {
112!
NEW
291
          factor_to_seconds = 1e-3;
×
292
        } else if (units == "s") {
112!
293
          factor_to_seconds = 1.0;
112✔
NEW
294
        } else if (units == "min") {
×
NEW
295
          factor_to_seconds = 1 / 60;
×
296
        } else {
NEW
297
          fatal_error("Invalid timestep unit, " + units);
×
298
        }
299
        dt *= factor_to_seconds;
112✔
300
      } else {
NEW
301
        fatal_error("Specify dt in settings XML");
×
302
      }
303
    } else {
112✔
NEW
304
      fatal_error("Specify timestep units in settings XML");
×
305
    }
306
  }
307

308
  // Random ray variables
309
  if (solver_type == SolverType::RANDOM_RAY) {
7,226✔
310
    xml_node random_ray_node = node_base.child("random_ray");
753✔
311
    if (check_for_node(random_ray_node, "distance_active")) {
753!
312
      RandomRay::distance_active_ =
753✔
313
        std::stod(get_node_value(random_ray_node, "distance_active"));
753✔
314
      if (RandomRay::distance_active_ <= 0.0) {
753!
315
        fatal_error("Random ray active distance must be greater than 0");
×
316
      }
317
    } else {
318
      fatal_error("Specify random ray active distance in settings XML");
×
319
    }
320
    if (check_for_node(random_ray_node, "distance_inactive")) {
753!
321
      RandomRay::distance_inactive_ =
753✔
322
        std::stod(get_node_value(random_ray_node, "distance_inactive"));
753✔
323
      if (RandomRay::distance_inactive_ < 0) {
753!
324
        fatal_error(
×
325
          "Random ray inactive distance must be greater than or equal to 0");
326
      }
327
    } else {
328
      fatal_error("Specify random ray inactive distance in settings XML");
×
329
    }
330
    if (check_for_node(random_ray_node, "source")) {
753!
331
      xml_node source_node = random_ray_node.child("source");
753✔
332
      // Get point to list of <source> elements and make sure there is at least
333
      // one
334
      RandomRay::ray_source_ = Source::create(source_node);
753✔
335
    } else {
336
      fatal_error("Specify random ray source in settings XML");
×
337
    }
338
    if (check_for_node(random_ray_node, "volume_estimator")) {
753✔
339
      std::string temp_str =
340
        get_node_value(random_ray_node, "volume_estimator", true, true);
129✔
341
      if (temp_str == "simulation_averaged") {
129✔
342
        FlatSourceDomain::volume_estimator_ =
32✔
343
          RandomRayVolumeEstimator::SIMULATION_AVERAGED;
344
      } else if (temp_str == "naive") {
97✔
345
        FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::NAIVE;
65✔
346
      } else if (temp_str == "hybrid") {
32!
347
        FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID;
32✔
348
      } else {
349
        fatal_error("Unrecognized volume estimator: " + temp_str);
×
350
      }
351
    }
129✔
352
    if (check_for_node(random_ray_node, "source_shape")) {
753✔
353
      std::string temp_str =
354
        get_node_value(random_ray_node, "source_shape", true, true);
368✔
355
      if (temp_str == "flat") {
368✔
356
        RandomRay::source_shape_ = RandomRaySourceShape::FLAT;
64✔
357
      } else if (temp_str == "linear") {
304✔
358
        RandomRay::source_shape_ = RandomRaySourceShape::LINEAR;
256✔
359
        if (settings::kinetic_simulation) {
256!
NEW
360
          fatal_error(
×
361
            "linear source shapes unimplemented for kinetic simulations.");
362
        }
363
      } else if (temp_str == "linear_xy") {
48!
364
        RandomRay::source_shape_ = RandomRaySourceShape::LINEAR_XY;
48✔
365
        if (settings::kinetic_simulation) {
48!
NEW
366
          fatal_error(
×
367
            "linear_xy source shapes unimplemented for kinetic simulations.");
368
        }
369
      } else {
370
        fatal_error("Unrecognized source shape: " + temp_str);
×
371
      }
372
    }
368✔
373
    if (check_for_node(random_ray_node, "volume_normalized_flux_tallies")) {
753✔
374
      FlatSourceDomain::volume_normalized_flux_tallies_ =
544✔
375
        get_node_value_bool(random_ray_node, "volume_normalized_flux_tallies");
544✔
376
    }
377
    if (check_for_node(random_ray_node, "adjoint")) {
753✔
378
      FlatSourceDomain::adjoint_ =
32✔
379
        get_node_value_bool(random_ray_node, "adjoint");
32✔
380
    }
381
    if (check_for_node(random_ray_node, "sample_method")) {
753✔
382
      std::string temp_str =
383
        get_node_value(random_ray_node, "sample_method", true, true);
16✔
384
      if (temp_str == "prng") {
16!
385
        RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
×
386
      } else if (temp_str == "halton") {
16!
387
        RandomRay::sample_method_ = RandomRaySampleMethod::HALTON;
16✔
388
      } else {
389
        fatal_error("Unrecognized sample method: " + temp_str);
×
390
      }
391
    }
16✔
392
    if (check_for_node(random_ray_node, "source_region_meshes")) {
753✔
393
      pugi::xml_node node_source_region_meshes =
394
        random_ray_node.child("source_region_meshes");
305✔
395
      for (pugi::xml_node node_mesh :
305✔
396
        node_source_region_meshes.children("mesh")) {
979✔
397
        int mesh_id = std::stoi(node_mesh.attribute("id").value());
369✔
398
        for (pugi::xml_node node_domain : node_mesh.children("domain")) {
738✔
399
          int domain_id = std::stoi(node_domain.attribute("id").value());
369✔
400
          std::string domain_type = node_domain.attribute("type").value();
369✔
401
          Source::DomainType type;
402
          if (domain_type == "material") {
369✔
403
            type = Source::DomainType::MATERIAL;
32✔
404
          } else if (domain_type == "cell") {
337✔
405
            type = Source::DomainType::CELL;
32✔
406
          } else if (domain_type == "universe") {
305!
407
            type = Source::DomainType::UNIVERSE;
305✔
408
          } else {
409
            throw std::runtime_error("Unknown domain type: " + domain_type);
×
410
          }
411
          FlatSourceDomain::mesh_domain_map_[mesh_id].emplace_back(
369✔
412
            type, domain_id);
413
        }
369✔
414
      }
415
    }
416
    if (check_for_node(random_ray_node, "diagonal_stabilization_rho")) {
753✔
417
      FlatSourceDomain::diagonal_stabilization_rho_ = std::stod(
32✔
418
        get_node_value(random_ray_node, "diagonal_stabilization_rho"));
64✔
419
      if (FlatSourceDomain::diagonal_stabilization_rho_ < 0.0 ||
32!
420
          FlatSourceDomain::diagonal_stabilization_rho_ > 1.0) {
32!
421
        fatal_error("Random ray diagonal stabilization rho factor must be "
×
422
                    "between 0 and 1");
423
      }
424
    }
425
    if (kinetic_simulation) {
753✔
426
      if (check_for_node(random_ray_node, "bd_order")) {
112✔
427
        static int n = std::stod(get_node_value(random_ray_node, "bd_order"));
48!
428
        if (n < 1 || n > 6) {
48!
NEW
429
          fatal_error("Specified BD order of " + std::to_string(n) +
×
430
                      ". BD order must be between 1 and 6");
431
        } else {
432
          RandomRay::bd_order_ = n;
48✔
433
        }
434
      }
435
      if (check_for_node(random_ray_node, "time_derivative_method")) {
112!
436
        std::string temp_str =
NEW
437
          get_node_value(random_ray_node, "time_derivative_method", true, true);
×
NEW
438
        if (temp_str == "isotropic") {
×
NEW
439
          RandomRay::time_method_ = RandomRayTimeMethod::ISOTROPIC;
×
NEW
440
        } else if (temp_str == "propagation") {
×
NEW
441
          RandomRay::time_method_ = RandomRayTimeMethod::PROPAGATION;
×
442
        } else {
NEW
443
          fatal_error("Unrecognized time derivative method: " + temp_str);
×
444
        }
NEW
445
      }
×
446
    }
447
  }
448
}
7,226✔
449

450
void read_settings_xml()
1,381✔
451
{
452
  using namespace settings;
453
  using namespace pugi;
454
  // Check if settings.xml exists
455
  std::string filename = settings::path_input + "settings.xml";
1,381✔
456
  if (!file_exists(filename)) {
1,381✔
457
    if (run_mode != RunMode::PLOTTING) {
22!
458
      fatal_error("Could not find any XML input files! In order to run OpenMC, "
×
459
                  "you first need a set of input files; at a minimum, this "
460
                  "includes settings.xml, geometry.xml, and materials.xml or a "
461
                  "single model XML file. Please consult the user's guide at "
462
                  "https://docs.openmc.org for further information.");
463
    } else {
464
      // The settings.xml file is optional if we just want to make a plot.
465
      return;
22✔
466
    }
467
  }
468

469
  // Parse settings.xml file
470
  xml_document doc;
1,359✔
471
  auto result = doc.load_file(filename.c_str());
1,359✔
472
  if (!result) {
1,359!
473
    fatal_error("Error processing settings.xml file.");
×
474
  }
475

476
  // Get root element
477
  xml_node root = doc.document_element();
1,359✔
478

479
  // Verbosity
480
  if (check_for_node(root, "verbosity") && verbosity == -1) {
1,359!
481
    verbosity = std::stoi(get_node_value(root, "verbosity"));
182✔
482
  } else if (verbosity == -1) {
1,177!
483
    verbosity = 7;
1,177✔
484
  }
485

486
  // To this point, we haven't displayed any output since we didn't know what
487
  // the verbosity is. Now that we checked for it, show the title if necessary
488
  if (mpi::master) {
1,359✔
489
    if (verbosity >= 2)
1,134✔
490
      title();
962✔
491
  }
492

493
  write_message("Reading settings XML file...", 5);
1,359✔
494

495
  read_settings_xml(root);
1,359✔
496
}
1,371✔
497

498
void read_settings_xml(pugi::xml_node root)
7,825✔
499
{
500
  using namespace settings;
501
  using namespace pugi;
502

503
  // Find if a multi-group or continuous-energy simulation is desired
504
  if (check_for_node(root, "energy_mode")) {
7,825✔
505
    std::string temp_str = get_node_value(root, "energy_mode", true, true);
1,287✔
506
    if (temp_str == "mg" || temp_str == "multi-group") {
1,287!
507
      run_CE = false;
1,287✔
508
    } else if (temp_str == "ce" || temp_str == "continuous-energy") {
×
509
      run_CE = true;
×
510
    }
511
  }
1,287✔
512

513
  // Check for user meshes and allocate
514
  read_meshes(root);
7,825✔
515

516
  // Look for deprecated cross_sections.xml file in settings.xml
517
  if (check_for_node(root, "cross_sections")) {
7,825!
518
    warning(
×
519
      "Setting cross_sections in settings.xml has been deprecated."
520
      " The cross_sections are now set in materials.xml and the "
521
      "cross_sections input to materials.xml and the OPENMC_CROSS_SECTIONS"
522
      " environment variable will take precendent over setting "
523
      "cross_sections in settings.xml.");
524
    path_cross_sections = get_node_value(root, "cross_sections");
×
525
  }
526

527
  if (!run_CE) {
7,825✔
528
    // Scattering Treatments
529
    if (check_for_node(root, "max_order")) {
1,287✔
530
      max_order = std::stoi(get_node_value(root, "max_order"));
16✔
531
    } else {
532
      // Set to default of largest int - 1, which means to use whatever is
533
      // contained in library. This is largest int - 1 because for legendre
534
      // scattering, a value of 1 is added to the order; adding 1 to the largest
535
      // int gets you the largest negative integer, which is not what we want.
536
      max_order = std::numeric_limits<int>::max() - 1;
1,271✔
537
    }
538
  }
539

540
  // Check for a trigger node and get trigger information
541
  if (check_for_node(root, "trigger")) {
7,825✔
542
    xml_node node_trigger = root.child("trigger");
162✔
543

544
    // Check if trigger(s) are to be turned on
545
    trigger_on = get_node_value_bool(node_trigger, "active");
162✔
546

547
    if (trigger_on) {
162✔
548
      if (check_for_node(node_trigger, "max_batches")) {
146!
549
        n_max_batches = std::stoi(get_node_value(node_trigger, "max_batches"));
146✔
550
      } else {
551
        fatal_error("<max_batches> must be specified with triggers");
×
552
      }
553

554
      // Get the batch interval to check triggers
555
      if (!check_for_node(node_trigger, "batch_interval")) {
146✔
556
        trigger_predict = true;
16✔
557
      } else {
558
        trigger_batch_interval =
130✔
559
          std::stoi(get_node_value(node_trigger, "batch_interval"));
130✔
560
        if (trigger_batch_interval <= 0) {
130!
561
          fatal_error("Trigger batch interval must be greater than zero");
×
562
        }
563
      }
564
    }
565
  }
566

567
  // Check run mode if it hasn't been set from the command line
568
  xml_node node_mode;
7,825✔
569
  if (run_mode == RunMode::UNSET) {
7,825✔
570
    if (check_for_node(root, "run_mode")) {
7,258✔
571
      std::string temp_str = get_node_value(root, "run_mode", true, true);
7,226✔
572
      if (temp_str == "eigenvalue") {
7,226✔
573
        run_mode = RunMode::EIGENVALUE;
4,506✔
574
      } else if (temp_str == "fixed source") {
2,720✔
575
        run_mode = RunMode::FIXED_SOURCE;
2,688✔
576
      } else if (temp_str == "plot") {
32!
577
        run_mode = RunMode::PLOTTING;
×
578
      } else if (temp_str == "particle restart") {
32!
579
        run_mode = RunMode::PARTICLE;
×
580
      } else if (temp_str == "volume") {
32!
581
        run_mode = RunMode::VOLUME;
32✔
582
      } else {
583
        fatal_error("Unrecognized run mode: " + temp_str);
×
584
      }
585

586
      // Assume XML specifies <particles>, <batches>, etc. directly
587
      node_mode = root;
7,226✔
588
    } else {
7,226✔
589
      warning("<run_mode> should be specified.");
32✔
590

591
      // Make sure that either eigenvalue or fixed source was specified
592
      node_mode = root.child("eigenvalue");
32✔
593
      if (node_mode) {
32!
594
        run_mode = RunMode::EIGENVALUE;
32✔
595
      } else {
596
        node_mode = root.child("fixed_source");
×
597
        if (node_mode) {
×
598
          run_mode = RunMode::FIXED_SOURCE;
×
599
        } else {
600
          fatal_error("<eigenvalue> or <fixed_source> not specified.");
×
601
        }
602
      }
603
    }
604
  }
605

606
  // Check solver type
607
  if (check_for_node(root, "random_ray")) {
7,825✔
608
    solver_type = SolverType::RANDOM_RAY;
753✔
609
    if (run_CE)
753!
610
      fatal_error("multi-group energy mode must be specified in settings XML "
×
611
                  "when using the random ray solver.");
612
  }
613

614
  if (run_mode == RunMode::EIGENVALUE || run_mode == RunMode::FIXED_SOURCE) {
7,825✔
615
    // Read run parameters
616
    get_run_parameters(node_mode);
7,226✔
617

618
    // Check number of active batches, inactive batches, max lost particles and
619
    // particles
620
    if (n_batches <= n_inactive) {
7,226!
621
      fatal_error("Number of active batches must be greater than zero.");
×
622
    } else if (n_inactive < 0) {
7,226!
623
      fatal_error("Number of inactive batches must be non-negative.");
×
624
    } else if (n_particles <= 0) {
7,226!
625
      fatal_error("Number of particles must be greater than zero.");
×
626
    } else if (max_lost_particles <= 0) {
7,226!
627
      fatal_error("Number of max lost particles must be greater than zero.");
×
628
    } else if (rel_max_lost_particles <= 0.0 || rel_max_lost_particles >= 1.0) {
7,226!
629
      fatal_error("Relative max lost particles must be between zero and one.");
×
630
    }
631

632
    // Check for user value for the number of generation of the Iterated Fission
633
    // Probability (IFP) method
634
    if (check_for_node(root, "ifp_n_generation")) {
7,226✔
635
      ifp_n_generation = std::stoi(get_node_value(root, "ifp_n_generation"));
85✔
636
      if (ifp_n_generation <= 0) {
85!
637
        fatal_error("'ifp_n_generation' must be greater than 0.");
×
638
      }
639
      // Avoid tallying 0 if IFP logs are not complete when active cycles start
640
      if (ifp_n_generation > n_inactive) {
85✔
641
        fatal_error("'ifp_n_generation' must be lower than or equal to the "
9✔
642
                    "number of inactive cycles.");
643
      }
644
    }
645
  }
646

647
  // Copy plotting random number seed if specified
648
  if (check_for_node(root, "plot_seed")) {
7,816!
649
    auto seed = std::stoll(get_node_value(root, "plot_seed"));
×
650
    model::plotter_seed = seed;
×
651
  }
652

653
  // Copy random number seed if specified
654
  if (check_for_node(root, "seed")) {
7,816✔
655
    auto seed = std::stoll(get_node_value(root, "seed"));
589✔
656
    openmc_set_seed(seed);
589✔
657
  }
658

659
  // Copy random number stride if specified
660
  if (check_for_node(root, "stride")) {
7,816✔
661
    auto stride = std::stoull(get_node_value(root, "stride"));
16✔
662
    openmc_set_stride(stride);
16✔
663
  }
664

665
  // Check for electron treatment
666
  if (check_for_node(root, "electron_treatment")) {
7,816✔
667
    auto temp_str = get_node_value(root, "electron_treatment", true, true);
70✔
668
    if (temp_str == "led") {
70✔
669
      electron_treatment = ElectronTreatment::LED;
11✔
670
    } else if (temp_str == "ttb") {
59!
671
      electron_treatment = ElectronTreatment::TTB;
59✔
672
    } else {
673
      fatal_error("Unrecognized electron treatment: " + temp_str + ".");
×
674
    }
675
  }
70✔
676

677
  // Check for photon transport
678
  if (check_for_node(root, "photon_transport")) {
7,816✔
679
    photon_transport = get_node_value_bool(root, "photon_transport");
199✔
680

681
    if (!run_CE && photon_transport) {
199!
682
      fatal_error("Photon transport is not currently supported in "
×
683
                  "multigroup mode");
684
    }
685
  }
686

687
  // Number of bins for logarithmic grid
688
  if (check_for_node(root, "log_grid_bins")) {
7,816✔
689
    n_log_bins = std::stoi(get_node_value(root, "log_grid_bins"));
16✔
690
    if (n_log_bins < 1) {
16!
691
      fatal_error("Number of bins for logarithmic grid must be greater "
×
692
                  "than zero.");
693
    }
694
  }
695

696
  // Number of OpenMP threads
697
  if (check_for_node(root, "threads")) {
7,816!
698
    if (mpi::master)
×
699
      warning("The <threads> element has been deprecated. Use "
×
700
              "the OMP_NUM_THREADS environment variable to set the number of "
701
              "threads.");
702
  }
703

704
  // ==========================================================================
705
  // EXTERNAL SOURCE
706

707
  // Get point to list of <source> elements and make sure there is at least one
708
  for (pugi::xml_node node : root.children("source")) {
15,077✔
709
    model::external_sources.push_back(Source::create(node));
7,271✔
710
  }
711

712
  // Check if the user has specified to read surface source
713
  if (check_for_node(root, "surf_source_read")) {
7,806✔
714
    surf_source_read = true;
32✔
715
    // Get surface source read node
716
    xml_node node_ssr = root.child("surf_source_read");
32✔
717

718
    std::string path = "surface_source.h5";
32✔
719
    // Check if the user has specified different file for surface source reading
720
    if (check_for_node(node_ssr, "path")) {
32!
721
      path = get_node_value(node_ssr, "path", false, true);
32✔
722
    }
723
    model::external_sources.push_back(make_unique<FileSource>(path));
32✔
724
  }
32✔
725

726
  // If no source specified, default to isotropic point source at origin with
727
  // Watt spectrum. No default source is needed in random ray mode.
728
  if (model::external_sources.empty() &&
9,972✔
729
      settings::solver_type != SolverType::RANDOM_RAY) {
2,166✔
730
    double T[] {0.0};
2,006✔
731
    double p[] {1.0};
2,006✔
732
    model::external_sources.push_back(make_unique<IndependentSource>(
2,006✔
733
      UPtrSpace {new SpatialPoint({0.0, 0.0, 0.0})},
4,012✔
734
      UPtrAngle {new Isotropic()}, UPtrDist {new Watt(0.988e6, 2.249e-6)},
4,012✔
735
      UPtrDist {new Discrete(T, p, 1)}));
4,012✔
736
  }
737

738
  // Build probability mass function for sampling external sources
739
  vector<double> source_strengths;
7,806✔
740
  for (auto& s : model::external_sources) {
17,105✔
741
    source_strengths.push_back(s->strength());
9,299✔
742
  }
743
  model::external_sources_probability.assign(source_strengths);
7,806✔
744

745
  // Check if we want to write out source
746
  if (check_for_node(root, "write_initial_source")) {
7,806!
747
    write_initial_source = get_node_value_bool(root, "write_initial_source");
×
748
  }
749

750
  // Get relative number of lost particles
751
  if (check_for_node(root, "source_rejection_fraction")) {
7,806✔
752
    source_rejection_fraction =
6✔
753
      std::stod(get_node_value(root, "source_rejection_fraction"));
6!
754
  }
755

756
  if (check_for_node(root, "free_gas_threshold")) {
7,806!
757
    free_gas_threshold = std::stod(get_node_value(root, "free_gas_threshold"));
×
758
  }
759

760
  // Survival biasing
761
  if (check_for_node(root, "survival_biasing")) {
7,806✔
762
    survival_biasing = get_node_value_bool(root, "survival_biasing");
177✔
763
  }
764

765
  // Probability tables
766
  if (check_for_node(root, "ptables")) {
7,806✔
767
    urr_ptables_on = get_node_value_bool(root, "ptables");
16✔
768
  }
769

770
  // Cutoffs
771
  if (check_for_node(root, "cutoff")) {
7,806✔
772
    xml_node node_cutoff = root.child("cutoff");
145✔
773
    if (check_for_node(node_cutoff, "weight")) {
145✔
774
      weight_cutoff = std::stod(get_node_value(node_cutoff, "weight"));
16✔
775
    }
776
    if (check_for_node(node_cutoff, "weight_avg")) {
145✔
777
      weight_survive = std::stod(get_node_value(node_cutoff, "weight_avg"));
16✔
778
    }
779
    if (check_for_node(node_cutoff, "survival_normalization")) {
145!
780
      survival_normalization =
×
781
        get_node_value_bool(node_cutoff, "survival_normalization");
×
782
    }
783
    if (check_for_node(node_cutoff, "energy_neutron")) {
145✔
784
      energy_cutoff[0] =
32✔
785
        std::stod(get_node_value(node_cutoff, "energy_neutron"));
16✔
786
    } else if (check_for_node(node_cutoff, "energy")) {
129!
787
      warning("The use of an <energy> cutoff is deprecated and should "
×
788
              "be replaced by <energy_neutron>.");
789
      energy_cutoff[0] = std::stod(get_node_value(node_cutoff, "energy"));
×
790
    }
791
    if (check_for_node(node_cutoff, "energy_photon")) {
145✔
792
      energy_cutoff[1] =
172✔
793
        std::stod(get_node_value(node_cutoff, "energy_photon"));
86✔
794
    }
795
    if (check_for_node(node_cutoff, "energy_electron")) {
145!
796
      energy_cutoff[2] =
×
797
        std::stof(get_node_value(node_cutoff, "energy_electron"));
×
798
    }
799
    if (check_for_node(node_cutoff, "energy_positron")) {
145!
800
      energy_cutoff[3] =
×
801
        std::stod(get_node_value(node_cutoff, "energy_positron"));
×
802
    }
803
    if (check_for_node(node_cutoff, "time_neutron")) {
145✔
804
      time_cutoff[0] = std::stod(get_node_value(node_cutoff, "time_neutron"));
27✔
805
    }
806
    if (check_for_node(node_cutoff, "time_photon")) {
145!
807
      time_cutoff[1] = std::stod(get_node_value(node_cutoff, "time_photon"));
×
808
    }
809
    if (check_for_node(node_cutoff, "time_electron")) {
145!
810
      time_cutoff[2] = std::stod(get_node_value(node_cutoff, "time_electron"));
×
811
    }
812
    if (check_for_node(node_cutoff, "time_positron")) {
145!
813
      time_cutoff[3] = std::stod(get_node_value(node_cutoff, "time_positron"));
×
814
    }
815
  }
816

817
  // Particle trace
818
  if (check_for_node(root, "trace")) {
7,806✔
819
    auto temp = get_node_array<int64_t>(root, "trace");
16✔
820
    if (temp.size() != 3) {
16!
821
      fatal_error("Must provide 3 integers for <trace> that specify the "
×
822
                  "batch, generation, and particle number.");
823
    }
824
    trace_batch = temp.at(0);
16✔
825
    trace_gen = temp.at(1);
16✔
826
    trace_particle = temp.at(2);
16✔
827
  }
16✔
828

829
  // Particle tracks
830
  if (check_for_node(root, "track")) {
7,806✔
831
    // Get values and make sure there are three per particle
832
    auto temp = get_node_array<int>(root, "track");
48✔
833
    if (temp.size() % 3 != 0) {
48!
834
      fatal_error(
×
835
        "Number of integers specified in 'track' is not "
836
        "divisible by 3.  Please provide 3 integers per particle to be "
837
        "tracked.");
838
    }
839

840
    // Reshape into track_identifiers
841
    int n_tracks = temp.size() / 3;
48✔
842
    for (int i = 0; i < n_tracks; ++i) {
192✔
843
      track_identifiers.push_back(
144✔
844
        {temp[3 * i], temp[3 * i + 1], temp[3 * i + 2]});
144✔
845
    }
846
  }
48✔
847

848
  // Shannon entropy
849
  if (solver_type == SolverType::RANDOM_RAY) {
7,806✔
850
    if (check_for_node(root, "entropy_mesh")) {
753!
851
      fatal_error("Random ray uses FSRs to compute the Shannon entropy. "
×
852
                  "No user-defined entropy mesh is supported.");
853
    }
854
    entropy_on = true;
753✔
855
  } else if (solver_type == SolverType::MONTE_CARLO) {
7,053!
856
    if (check_for_node(root, "entropy_mesh")) {
7,053✔
857
      int temp = std::stoi(get_node_value(root, "entropy_mesh"));
335✔
858
      if (model::mesh_map.find(temp) == model::mesh_map.end()) {
335!
859
        fatal_error(fmt::format(
×
860
          "Mesh {} specified for Shannon entropy does not exist.", temp));
861
      }
862

863
      auto* m = dynamic_cast<RegularMesh*>(
335!
864
        model::meshes[model::mesh_map.at(temp)].get());
335✔
865
      if (!m)
335!
866
        fatal_error("Only regular meshes can be used as an entropy mesh");
×
867
      simulation::entropy_mesh = m;
335✔
868

869
      // Turn on Shannon entropy calculation
870
      entropy_on = true;
335✔
871

872
    } else if (check_for_node(root, "entropy")) {
6,718!
873
      fatal_error(
×
874
        "Specifying a Shannon entropy mesh via the <entropy> element "
875
        "is deprecated. Please create a mesh using <mesh> and then reference "
876
        "it by specifying its ID in an <entropy_mesh> element.");
877
    }
878
  }
879
  // Uniform fission source weighting mesh
880
  if (check_for_node(root, "ufs_mesh")) {
7,806✔
881
    auto temp = std::stoi(get_node_value(root, "ufs_mesh"));
16✔
882
    if (model::mesh_map.find(temp) == model::mesh_map.end()) {
16!
883
      fatal_error(fmt::format("Mesh {} specified for uniform fission site "
×
884
                              "method does not exist.",
885
        temp));
886
    }
887

888
    auto* m =
889
      dynamic_cast<RegularMesh*>(model::meshes[model::mesh_map.at(temp)].get());
16!
890
    if (!m)
16!
891
      fatal_error("Only regular meshes can be used as a UFS mesh");
×
892
    simulation::ufs_mesh = m;
16✔
893

894
    // Turn on uniform fission source weighting
895
    ufs_on = true;
16✔
896

897
  } else if (check_for_node(root, "uniform_fs")) {
7,790!
898
    fatal_error(
×
899
      "Specifying a UFS mesh via the <uniform_fs> element "
900
      "is deprecated. Please create a mesh using <mesh> and then reference "
901
      "it by specifying its ID in a <ufs_mesh> element.");
902
  }
903

904
  // Check if the user has specified to write state points
905
  if (check_for_node(root, "state_point")) {
7,806✔
906

907
    // Get pointer to state_point node
908
    auto node_sp = root.child("state_point");
167✔
909

910
    // Determine number of batches at which to store state points
911
    if (check_for_node(node_sp, "batches")) {
167!
912
      // User gave specific batches to write state points
913
      auto temp = get_node_array<int>(node_sp, "batches");
167✔
914
      for (const auto& b : temp) {
512✔
915
        statepoint_batch.insert(b);
345✔
916
      }
917
    } else {
167✔
918
      // If neither were specified, write state point at last batch
919
      statepoint_batch.insert(n_batches);
×
920
    }
921
  } else {
922
    // If no <state_point> tag was present, by default write state point at
923
    // last batch only
924
    statepoint_batch.insert(n_batches);
7,639✔
925
  }
926

927
  // Check if the user has specified to write source points
928
  if (check_for_node(root, "source_point")) {
7,806✔
929
    // Get source_point node
930
    xml_node node_sp = root.child("source_point");
107✔
931

932
    // Determine batches at which to store source points
933
    if (check_for_node(node_sp, "batches")) {
107✔
934
      // User gave specific batches to write source points
935
      auto temp = get_node_array<int>(node_sp, "batches");
48✔
936
      for (const auto& b : temp) {
128✔
937
        sourcepoint_batch.insert(b);
80✔
938
      }
939
    } else {
48✔
940
      // If neither were specified, write source points with state points
941
      sourcepoint_batch = statepoint_batch;
59✔
942
    }
943

944
    // Check if the user has specified to write binary source file
945
    if (check_for_node(node_sp, "separate")) {
107✔
946
      source_separate = get_node_value_bool(node_sp, "separate");
75✔
947
    }
948
    if (check_for_node(node_sp, "write")) {
107!
949
      source_write = get_node_value_bool(node_sp, "write");
×
950
    }
951
    if (check_for_node(node_sp, "mcpl")) {
107✔
952
      source_mcpl_write = get_node_value_bool(node_sp, "mcpl");
27✔
953
    }
954
    if (check_for_node(node_sp, "overwrite_latest")) {
107✔
955
      source_latest = get_node_value_bool(node_sp, "overwrite_latest");
16✔
956
      source_separate = source_latest;
16✔
957
    }
958
  } else {
959
    // If no <source_point> tag was present, by default we keep source bank in
960
    // statepoint file and write it out at statepoints intervals
961
    source_separate = false;
7,699✔
962
    sourcepoint_batch = statepoint_batch;
7,699✔
963
  }
964

965
  // Check is the user specified to convert strength to statistical weight
966
  if (check_for_node(root, "uniform_source_sampling")) {
7,806✔
967
    uniform_source_sampling =
55✔
968
      get_node_value_bool(root, "uniform_source_sampling");
55✔
969
  }
970

971
  // Check if the user has specified to write surface source
972
  if (check_for_node(root, "surf_source_write")) {
7,806✔
973
    surf_source_write = true;
412✔
974
    // Get surface source write node
975
    xml_node node_ssw = root.child("surf_source_write");
412✔
976

977
    // Determine surface ids at which crossing particles are to be banked.
978
    // If no surfaces are specified, all surfaces in the model will be used
979
    // to bank source points.
980
    if (check_for_node(node_ssw, "surface_ids")) {
412✔
981
      auto temp = get_node_array<int>(node_ssw, "surface_ids");
202✔
982
      for (const auto& b : temp) {
994✔
983
        source_write_surf_id.insert(b);
792✔
984
      }
985
    }
202✔
986

987
    // Get maximum number of particles to be banked per surface
988
    if (check_for_node(node_ssw, "max_particles")) {
412✔
989
      ssw_max_particles = std::stoll(get_node_value(node_ssw, "max_particles"));
403✔
990
    } else {
991
      fatal_error("A maximum number of particles needs to be specified "
9✔
992
                  "using the 'max_particles' parameter to store surface "
993
                  "source points.");
994
    }
995

996
    // Get maximum number of surface source files to be created
997
    if (check_for_node(node_ssw, "max_source_files")) {
403✔
998
      ssw_max_files = std::stoll(get_node_value(node_ssw, "max_source_files"));
33✔
999
    } else {
1000
      ssw_max_files = 1;
370✔
1001
    }
1002

1003
    if (check_for_node(node_ssw, "mcpl")) {
403✔
1004
      surf_mcpl_write = get_node_value_bool(node_ssw, "mcpl");
11✔
1005
    }
1006
    // Get cell information
1007
    if (check_for_node(node_ssw, "cell")) {
403✔
1008
      ssw_cell_id = std::stoll(get_node_value(node_ssw, "cell"));
104✔
1009
      ssw_cell_type = SSWCellType::Both;
104✔
1010
    }
1011
    if (check_for_node(node_ssw, "cellfrom")) {
403✔
1012
      if (ssw_cell_id != C_NONE) {
90✔
1013
        fatal_error(
18✔
1014
          "'cell', 'cellfrom' and 'cellto' cannot be used at the same time.");
1015
      }
1016
      ssw_cell_id = std::stoll(get_node_value(node_ssw, "cellfrom"));
72✔
1017
      ssw_cell_type = SSWCellType::From;
72✔
1018
    }
1019
    if (check_for_node(node_ssw, "cellto")) {
385✔
1020
      if (ssw_cell_id != C_NONE) {
71✔
1021
        fatal_error(
18✔
1022
          "'cell', 'cellfrom' and 'cellto' cannot be used at the same time.");
1023
      }
1024
      ssw_cell_id = std::stoll(get_node_value(node_ssw, "cellto"));
53✔
1025
      ssw_cell_type = SSWCellType::To;
53✔
1026
    }
1027
  }
1028

1029
  // Check if the user has specified to write specific collisions
1030
  if (check_for_node(root, "collision_track")) {
7,761✔
1031
    settings::collision_track = true;
155✔
1032
    // Get collision track node
1033
    xml_node node_ct = root.child("collision_track");
155✔
1034
    collision_track_config = CollisionTrackConfig {};
155✔
1035

1036
    // Determine cell ids at which crossing particles are to be banked
1037
    if (check_for_node(node_ct, "cell_ids")) {
155✔
1038
      auto temp = get_node_array<int>(node_ct, "cell_ids");
81✔
1039
      for (const auto& b : temp) {
211✔
1040
        collision_track_config.cell_ids.insert(b);
130✔
1041
      }
1042
    }
81✔
1043
    if (check_for_node(node_ct, "reactions")) {
155✔
1044
      auto temp = get_node_array<std::string>(node_ct, "reactions");
65✔
1045
      for (const auto& b : temp) {
178✔
1046
        int reaction_int = reaction_type(b);
113✔
1047
        if (reaction_int > 0) {
113!
1048
          collision_track_config.mt_numbers.insert(reaction_int);
113✔
1049
        }
1050
      }
1051
    }
65✔
1052
    if (check_for_node(node_ct, "universe_ids")) {
155✔
1053
      auto temp = get_node_array<int>(node_ct, "universe_ids");
32✔
1054
      for (const auto& b : temp) {
64✔
1055
        collision_track_config.universe_ids.insert(b);
32✔
1056
      }
1057
    }
32✔
1058
    if (check_for_node(node_ct, "material_ids")) {
155✔
1059
      auto temp = get_node_array<int>(node_ct, "material_ids");
32✔
1060
      for (const auto& b : temp) {
80✔
1061
        collision_track_config.material_ids.insert(b);
48✔
1062
      }
1063
    }
32✔
1064
    if (check_for_node(node_ct, "nuclides")) {
155✔
1065
      auto temp = get_node_array<std::string>(node_ct, "nuclides");
32✔
1066
      for (const auto& b : temp) {
128✔
1067
        collision_track_config.nuclides.insert(b);
96✔
1068
      }
1069
    }
32✔
1070
    if (check_for_node(node_ct, "deposited_E_threshold")) {
155✔
1071
      collision_track_config.deposited_energy_threshold =
32✔
1072
        std::stod(get_node_value(node_ct, "deposited_E_threshold"));
32✔
1073
    }
1074
    // Get maximum number of particles to be banked per collision
1075
    if (check_for_node(node_ct, "max_collisions")) {
155!
1076
      collision_track_config.max_collisions =
155✔
1077
        std::stoll(get_node_value(node_ct, "max_collisions"));
155✔
1078
    } else {
1079
      warning("A maximum number of collisions needs to be specified. "
×
1080
              "By default the code sets 'max_collisions' parameter equals to "
1081
              "1000.");
1082
    }
1083
    // Get maximum number of collision_track files to be created
1084
    if (check_for_node(node_ct, "max_collision_track_files")) {
155!
1085
      collision_track_config.max_files =
×
1086
        std::stoll(get_node_value(node_ct, "max_collision_track_files"));
×
1087
    }
1088
    if (check_for_node(node_ct, "mcpl")) {
155✔
1089
      collision_track_config.mcpl_write = get_node_value_bool(node_ct, "mcpl");
22✔
1090
    }
1091
  }
1092

1093
  // If source is not separate and is to be written out in the statepoint
1094
  // file, make sure that the sourcepoint batch numbers are contained in the
1095
  // statepoint list
1096
  if (!source_separate) {
7,761✔
1097
    for (const auto& b : sourcepoint_batch) {
15,470✔
1098
      if (!contains(statepoint_batch, b)) {
7,800!
1099
        fatal_error(
×
1100
          "Sourcepoint batches are not a subset of statepoint batches.");
1101
      }
1102
    }
1103
  }
1104

1105
  // Check if the user has specified to not reduce tallies at the end of every
1106
  // batch
1107
  if (check_for_node(root, "no_reduce")) {
7,761✔
1108
    reduce_tallies = !get_node_value_bool(root, "no_reduce");
32✔
1109
  }
1110

1111
  // Check if the user has specified to use confidence intervals for
1112
  // uncertainties rather than standard deviations
1113
  if (check_for_node(root, "confidence_intervals")) {
7,761✔
1114
    confidence_intervals = get_node_value_bool(root, "confidence_intervals");
16✔
1115
  }
1116

1117
  // Check for output options
1118
  if (check_for_node(root, "output")) {
7,761✔
1119
    // Get pointer to output node
1120
    pugi::xml_node node_output = root.child("output");
672✔
1121

1122
    // Check for summary option
1123
    if (check_for_node(node_output, "summary")) {
672✔
1124
      output_summary = get_node_value_bool(node_output, "summary");
645✔
1125
    }
1126

1127
    // Check for ASCII tallies output option
1128
    if (check_for_node(node_output, "tallies")) {
672✔
1129
      output_tallies = get_node_value_bool(node_output, "tallies");
250✔
1130
    }
1131

1132
    // Set output directory if a path has been specified
1133
    if (check_for_node(node_output, "path")) {
672!
1134
      path_output = get_node_value(node_output, "path");
×
1135
      if (!ends_with(path_output, "/")) {
×
1136
        path_output += "/";
×
1137
      }
1138
    }
1139
  }
1140

1141
  // Resonance scattering parameters
1142
  if (check_for_node(root, "resonance_scattering")) {
7,761✔
1143
    xml_node node_res_scat = root.child("resonance_scattering");
16✔
1144

1145
    // See if resonance scattering is enabled
1146
    if (check_for_node(node_res_scat, "enable")) {
16!
1147
      res_scat_on = get_node_value_bool(node_res_scat, "enable");
16✔
1148
    } else {
1149
      res_scat_on = true;
×
1150
    }
1151

1152
    // Determine what method is used
1153
    if (check_for_node(node_res_scat, "method")) {
16!
1154
      auto temp = get_node_value(node_res_scat, "method", true, true);
16✔
1155
      if (temp == "rvs") {
16!
1156
        res_scat_method = ResScatMethod::rvs;
16✔
1157
      } else if (temp == "dbrc") {
×
1158
        res_scat_method = ResScatMethod::dbrc;
×
1159
      } else {
1160
        fatal_error(
×
1161
          "Unrecognized resonance elastic scattering method: " + temp + ".");
×
1162
      }
1163
    }
16✔
1164

1165
    // Minimum energy for resonance scattering
1166
    if (check_for_node(node_res_scat, "energy_min")) {
16!
1167
      res_scat_energy_min =
16✔
1168
        std::stod(get_node_value(node_res_scat, "energy_min"));
16✔
1169
    }
1170
    if (res_scat_energy_min < 0.0) {
16!
1171
      fatal_error("Lower resonance scattering energy bound is negative");
×
1172
    }
1173

1174
    // Maximum energy for resonance scattering
1175
    if (check_for_node(node_res_scat, "energy_max")) {
16!
1176
      res_scat_energy_max =
16✔
1177
        std::stod(get_node_value(node_res_scat, "energy_max"));
16✔
1178
    }
1179
    if (res_scat_energy_max < res_scat_energy_min) {
16!
1180
      fatal_error("Upper resonance scattering energy bound is below the "
×
1181
                  "lower resonance scattering energy bound.");
1182
    }
1183

1184
    // Get resonance scattering nuclides
1185
    if (check_for_node(node_res_scat, "nuclides")) {
16!
1186
      res_scat_nuclides =
1187
        get_node_array<std::string>(node_res_scat, "nuclides");
16✔
1188
    }
1189
  }
1190

1191
  // Get volume calculations
1192
  for (pugi::xml_node node_vol : root.children("volume_calc")) {
8,087✔
1193
    model::volume_calcs.emplace_back(node_vol);
326✔
1194
  }
1195

1196
  // Get temperature settings
1197
  if (check_for_node(root, "temperature_default")) {
7,761✔
1198
    temperature_default =
172✔
1199
      std::stod(get_node_value(root, "temperature_default"));
172✔
1200
  }
1201
  if (check_for_node(root, "temperature_method")) {
7,761✔
1202
    auto temp = get_node_value(root, "temperature_method", true, true);
325✔
1203
    if (temp == "nearest") {
325✔
1204
      temperature_method = TemperatureMethod::NEAREST;
140✔
1205
    } else if (temp == "interpolation") {
185!
1206
      temperature_method = TemperatureMethod::INTERPOLATION;
185✔
1207
    } else {
1208
      fatal_error("Unknown temperature method: " + temp);
×
1209
    }
1210
  }
325✔
1211
  if (check_for_node(root, "temperature_tolerance")) {
7,761✔
1212
    temperature_tolerance =
171✔
1213
      std::stod(get_node_value(root, "temperature_tolerance"));
171✔
1214
  }
1215
  if (check_for_node(root, "temperature_multipole")) {
7,761✔
1216
    temperature_multipole = get_node_value_bool(root, "temperature_multipole");
32✔
1217

1218
    // Multipole currently doesn't work with photon transport
1219
    if (temperature_multipole && photon_transport) {
32!
1220
      fatal_error("Multipole data cannot currently be used in conjunction with "
×
1221
                  "photon transport.");
1222
    }
1223
  }
1224
  if (check_for_node(root, "temperature_range")) {
7,761!
1225
    auto range = get_node_array<double>(root, "temperature_range");
×
1226
    temperature_range[0] = range.at(0);
×
1227
    temperature_range[1] = range.at(1);
×
1228
  }
×
1229

1230
  // Check for tabular_legendre options
1231
  if (check_for_node(root, "tabular_legendre")) {
7,761✔
1232
    // Get pointer to tabular_legendre node
1233
    xml_node node_tab_leg = root.child("tabular_legendre");
96✔
1234

1235
    // Check for enable option
1236
    if (check_for_node(node_tab_leg, "enable")) {
96!
1237
      legendre_to_tabular = get_node_value_bool(node_tab_leg, "enable");
96✔
1238
    }
1239

1240
    // Check for the number of points
1241
    if (check_for_node(node_tab_leg, "num_points")) {
96!
1242
      legendre_to_tabular_points =
×
1243
        std::stoi(get_node_value(node_tab_leg, "num_points"));
×
1244
      if (legendre_to_tabular_points <= 1 && !run_CE) {
×
1245
        fatal_error(
×
1246
          "The 'num_points' subelement/attribute of the "
1247
          "<tabular_legendre> element must contain a value greater than 1");
1248
      }
1249
    }
1250
  }
1251

1252
  // Check whether create delayed neutrons in fission
1253
  if (check_for_node(root, "create_delayed_neutrons")) {
7,761!
1254
    create_delayed_neutrons =
×
1255
      get_node_value_bool(root, "create_delayed_neutrons");
×
1256
  }
1257

1258
  // Check whether create fission sites
1259
  if (run_mode == RunMode::FIXED_SOURCE) {
7,761✔
1260
    if (check_for_node(root, "create_fission_neutrons")) {
2,642✔
1261
      create_fission_neutrons =
193✔
1262
        get_node_value_bool(root, "create_fission_neutrons");
193✔
1263
    }
1264
  }
1265

1266
  // Check whether to scale fission photon yields
1267
  if (check_for_node(root, "delayed_photon_scaling")) {
7,761!
1268
    delayed_photon_scaling =
×
1269
      get_node_value_bool(root, "delayed_photon_scaling");
×
1270
  }
1271

1272
  // Check whether to use event-based parallelism
1273
  if (check_for_node(root, "event_based")) {
7,761!
1274
    event_based = get_node_value_bool(root, "event_based");
×
1275
  }
1276

1277
  // Check whether material cell offsets should be generated
1278
  if (check_for_node(root, "material_cell_offsets")) {
7,761!
1279
    material_cell_offsets = get_node_value_bool(root, "material_cell_offsets");
×
1280
  }
1281

1282
  // Weight window information
1283
  for (pugi::xml_node node_ww : root.children("weight_windows")) {
7,853✔
1284
    variance_reduction::weight_windows.emplace_back(
92✔
1285
      std::make_unique<WeightWindows>(node_ww));
184✔
1286
  }
1287

1288
  // Enable weight windows by default if one or more are present
1289
  if (variance_reduction::weight_windows.size() > 0)
7,761✔
1290
    settings::weight_windows_on = true;
65✔
1291

1292
  // read weight windows from file
1293
  if (check_for_node(root, "weight_windows_file")) {
7,761!
1294
    weight_windows_file = get_node_value(root, "weight_windows_file");
×
1295
  }
1296

1297
  // read settings for weight windows value, this will override
1298
  // the automatic setting even if weight windows are present
1299
  if (check_for_node(root, "weight_windows_on")) {
7,761✔
1300
    weight_windows_on = get_node_value_bool(root, "weight_windows_on");
38✔
1301
  }
1302

1303
  if (check_for_node(root, "max_secondaries")) {
7,761!
1304
    settings::max_secondaries =
×
1305
      std::stoi(get_node_value(root, "max_secondaries"));
×
1306
  }
1307

1308
  if (check_for_node(root, "max_history_splits")) {
7,761✔
1309
    settings::max_history_splits =
215✔
1310
      std::stoi(get_node_value(root, "max_history_splits"));
215✔
1311
  }
1312

1313
  if (check_for_node(root, "max_tracks")) {
7,761✔
1314
    settings::max_tracks = std::stoi(get_node_value(root, "max_tracks"));
48✔
1315
  }
1316

1317
  // Create weight window generator objects
1318
  if (check_for_node(root, "weight_window_generators")) {
7,761✔
1319
    auto wwgs_node = root.child("weight_window_generators");
82✔
1320
    for (pugi::xml_node node_wwg :
82✔
1321
      wwgs_node.children("weight_windows_generator")) {
246✔
1322
      variance_reduction::weight_windows_generators.emplace_back(
82✔
1323
        std::make_unique<WeightWindowsGenerator>(node_wwg));
164✔
1324
    }
1325
    // if any of the weight windows are intended to be generated otf, make
1326
    // sure they're applied
1327
    for (const auto& wwg : variance_reduction::weight_windows_generators) {
82!
1328
      if (wwg->on_the_fly_) {
82!
1329
        settings::weight_windows_on = true;
82✔
1330
        break;
82✔
1331
      }
1332
    }
1333
  }
1334

1335
  // Set up weight window checkpoints
1336
  if (check_for_node(root, "weight_window_checkpoints")) {
7,761✔
1337
    xml_node ww_checkpoints = root.child("weight_window_checkpoints");
2!
1338
    if (check_for_node(ww_checkpoints, "collision")) {
2!
1339
      weight_window_checkpoint_collision =
2✔
1340
        get_node_value_bool(ww_checkpoints, "collision");
2!
1341
    }
1342
    if (check_for_node(ww_checkpoints, "surface")) {
2!
1343
      weight_window_checkpoint_surface =
2✔
1344
        get_node_value_bool(ww_checkpoints, "surface");
2!
1345
    }
1346
  }
1347

1348
  if (check_for_node(root, "use_decay_photons")) {
7,761✔
1349
    settings::use_decay_photons =
11✔
1350
      get_node_value_bool(root, "use_decay_photons");
11✔
1351
  }
1352
}
7,761✔
1353

1354
void free_memory_settings()
7,901✔
1355
{
1356
  settings::statepoint_batch.clear();
7,901✔
1357
  settings::sourcepoint_batch.clear();
7,901✔
1358
  settings::source_write_surf_id.clear();
7,901✔
1359
  settings::res_scat_nuclides.clear();
7,901✔
1360
}
7,901✔
1361

1362
//==============================================================================
1363
// C API functions
1364
//==============================================================================
1365

1366
extern "C" int openmc_set_n_batches(
44✔
1367
  int32_t n_batches, bool set_max_batches, bool add_statepoint_batch)
1368
{
1369
  if (settings::n_inactive >= n_batches) {
44✔
1370
    set_errmsg("Number of active batches must be greater than zero.");
11✔
1371
    return OPENMC_E_INVALID_ARGUMENT;
11✔
1372
  }
1373

1374
  if (!settings::trigger_on) {
33✔
1375
    // Set n_batches and n_max_batches to same value
1376
    settings::n_batches = n_batches;
11✔
1377
    settings::n_max_batches = n_batches;
11✔
1378
  } else {
1379
    // Set n_batches and n_max_batches based on value of set_max_batches
1380
    if (set_max_batches) {
22✔
1381
      settings::n_max_batches = n_batches;
11✔
1382
    } else {
1383
      settings::n_batches = n_batches;
11✔
1384
    }
1385
  }
1386

1387
  // Update size of k_generation and entropy
1388
  int m = settings::n_max_batches * settings::gen_per_batch;
33✔
1389
  simulation::k_generation.reserve(m);
33✔
1390
  simulation::entropy.reserve(m);
33✔
1391

1392
  // Add value of n_batches to statepoint_batch
1393
  if (add_statepoint_batch &&
55✔
1394
      !(contains(settings::statepoint_batch, n_batches)))
22!
1395
    settings::statepoint_batch.insert(n_batches);
22✔
1396

1397
  return 0;
33✔
1398
}
1399

1400
extern "C" int openmc_get_n_batches(int* n_batches, bool get_max_batches)
2,530✔
1401
{
1402
  *n_batches = get_max_batches ? settings::n_max_batches : settings::n_batches;
2,530✔
1403

1404
  return 0;
2,530✔
1405
}
1406

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