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

openmc-dev / openmc / 23011177584

12 Mar 2026 03:57PM UTC coverage: 80.895% (-0.7%) from 81.566%
23011177584

Pull #3863

github

web-flow
Merge 00a6e11ae into ba94c5823
Pull Request #3863: Shared Secondary Particle Bank

16504 of 23690 branches covered (69.67%)

Branch coverage included in aggregate %.

261 of 367 new or added lines in 17 files covered. (71.12%)

609 existing lines in 46 files now uncovered.

56448 of 66491 relevant lines covered (84.9%)

25290933.62 hits per line

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

71.24
/src/initialize.cpp
1
#include "openmc/initialize.h"
2

3
#include <clocale>
4
#include <cstddef>
5
#include <cstdlib> // for getenv
6
#include <cstring>
7
#include <string>
8

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

14
#include "openmc/capi.h"
15
#include "openmc/chain.h"
16
#include "openmc/constants.h"
17
#include "openmc/cross_sections.h"
18
#include "openmc/error.h"
19
#include "openmc/file_utils.h"
20
#include "openmc/geometry_aux.h"
21
#include "openmc/hdf5_interface.h"
22
#include "openmc/material.h"
23
#include "openmc/memory.h"
24
#include "openmc/message_passing.h"
25
#include "openmc/mgxs_interface.h"
26
#include "openmc/nuclide.h"
27
#include "openmc/openmp_interface.h"
28
#include "openmc/output.h"
29
#include "openmc/plot.h"
30
#include "openmc/random_lcg.h"
31
#include "openmc/settings.h"
32
#include "openmc/simulation.h"
33
#include "openmc/string_utils.h"
34
#include "openmc/summary.h"
35
#include "openmc/tallies/tally.h"
36
#include "openmc/thermal.h"
37
#include "openmc/timer.h"
38
#include "openmc/vector.h"
39
#include "openmc/weight_windows.h"
40

41
#ifdef OPENMC_LIBMESH_ENABLED
42
#include "libmesh/libmesh.h"
43
#endif
44

45
int openmc_init(int argc, char* argv[], const void* intracomm)
2,711✔
46
{
47
  using namespace openmc;
2,711✔
48

49
#ifdef OPENMC_MPI
50
  // Check if intracomm was passed
51
  MPI_Comm comm;
52
  if (intracomm) {
53
    comm = *static_cast<const MPI_Comm*>(intracomm);
54
  } else {
55
    comm = MPI_COMM_WORLD;
56
  }
57

58
  // Initialize MPI for C++
59
  initialize_mpi(comm);
60
#endif
61

62
  // Parse command-line arguments
63
  int err = parse_command_line(argc, argv);
2,711✔
64
  if (err)
2,699✔
65
    return err;
66

67
#ifdef OPENMC_LIBMESH_ENABLED
68
  const int n_threads = num_threads();
69
  // initialize libMesh if it hasn't been initialized already
70
  // (if initialized externally, the libmesh_init object needs to be provided
71
  // also)
72
  if (!settings::libmesh_init && !libMesh::initialized()) {
73
#ifdef OPENMC_MPI
74
    // pass command line args, empty MPI communicator, and number of threads.
75
    // Because libMesh was not initialized, we assume that OpenMC is the primary
76
    // application and that its main MPI comm should be used.
77
    settings::libmesh_init =
78
      make_unique<libMesh::LibMeshInit>(argc, argv, comm, n_threads);
79
#else
80
    // pass command line args, empty MPI communicator, and number of threads
81
    settings::libmesh_init =
82
      make_unique<libMesh::LibMeshInit>(argc, argv, 0, n_threads);
83
#endif
84

85
    settings::libmesh_comm = &(settings::libmesh_init->comm());
86
  }
87

88
#endif
89

90
  // Start total and initialization timer
91
  simulation::time_total.start();
2,695✔
92
  simulation::time_initialize.start();
2,695✔
93

94
#ifdef _OPENMP
95
  // If OMP_SCHEDULE is not set, default to a static schedule
96
  char* envvar = std::getenv("OMP_SCHEDULE");
673✔
97
  if (!envvar) {
673!
98
    omp_set_schedule(omp_sched_static, 0);
673✔
99
  }
100
#endif
101

102
  // Initialize random number generator -- if the user specifies a seed and/or
103
  // stride, it will be re-initialized later
104
  openmc::openmc_set_seed(DEFAULT_SEED);
2,695✔
105
  openmc::openmc_set_stride(DEFAULT_STRIDE);
2,695✔
106

107
  // Copy previous locale and set locale to C. This is a workaround for an issue
108
  // whereby when openmc_init is called from the plotter, the Qt application
109
  // framework first calls std::setlocale, which affects how pugixml reads
110
  // floating point numbers due to a bug:
111
  // https://github.com/zeux/pugixml/issues/469
112
  std::string prev_locale = std::setlocale(LC_ALL, nullptr);
2,695✔
113
  if (std::setlocale(LC_ALL, "C") == NULL) {
2,695!
114
    fatal_error("Cannot set locale to C.");
×
115
  }
116

117
  // Read XML input files
118
  if (!read_model_xml())
2,695✔
119
    read_separate_xml_files();
428✔
120

121
  // Reset locale to previous state
122
  if (std::setlocale(LC_ALL, prev_locale.c_str()) == NULL) {
2,655!
123
    fatal_error("Cannot reset locale.");
×
124
  }
125

126
  // Write some initial output under the header if needed
127
  initial_output();
2,655✔
128

129
  // Check for particle restart run
130
  if (settings::particle_restart_run)
2,655✔
131
    settings::run_mode = RunMode::PARTICLE;
16✔
132

133
  // Stop initialization timer
134
  simulation::time_initialize.stop();
2,655✔
135
  simulation::time_total.stop();
2,655✔
136

137
  return 0;
2,655✔
138
}
2,659✔
139

140
namespace openmc {
141

142
#ifdef OPENMC_MPI
143
void initialize_mpi(MPI_Comm intracomm)
144
{
145
  mpi::intracomm = intracomm;
146

147
  // Initialize MPI
148
  int flag;
149
  MPI_Initialized(&flag);
150
  if (!flag)
151
    MPI_Init(nullptr, nullptr);
152

153
  // Determine number of processes and rank for each
154
  MPI_Comm_size(intracomm, &mpi::n_procs);
155
  MPI_Comm_rank(intracomm, &mpi::rank);
156
  mpi::master = (mpi::rank == 0);
157

158
  // Create bank datatype
159
  SourceSite b;
160
  MPI_Aint disp[14];
161
  MPI_Get_address(&b.r, &disp[0]);
162
  MPI_Get_address(&b.u, &disp[1]);
163
  MPI_Get_address(&b.E, &disp[2]);
164
  MPI_Get_address(&b.time, &disp[3]);
165
  MPI_Get_address(&b.wgt, &disp[4]);
166
  MPI_Get_address(&b.delayed_group, &disp[5]);
167
  MPI_Get_address(&b.surf_id, &disp[6]);
168
  MPI_Get_address(&b.particle, &disp[7]);
169
  MPI_Get_address(&b.parent_nuclide, &disp[8]);
170
  MPI_Get_address(&b.parent_id, &disp[9]);
171
  MPI_Get_address(&b.progeny_id, &disp[10]);
172
  MPI_Get_address(&b.wgt_born, &disp[11]);
173
  MPI_Get_address(&b.wgt_ww_born, &disp[12]);
174
  MPI_Get_address(&b.n_split, &disp[13]);
175
  for (int i = 13; i >= 0; --i) {
176
    disp[i] -= disp[0];
177
  }
178

179
  // Block counts for each field
180
  int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
181

182
  // Types for each field
183
  MPI_Datatype types[] = {
184
    MPI_DOUBLE,  // r (3 doubles)
185
    MPI_DOUBLE,  // u (3 doubles)
186
    MPI_DOUBLE,  // E
187
    MPI_DOUBLE,  // time
188
    MPI_DOUBLE,  // wgt
189
    MPI_INT,     // delayed_group
190
    MPI_INT,     // surf_id
191
    MPI_INT,     // particle (enum)
192
    MPI_INT,     // parent_nuclide
193
    MPI_INT64_T, // parent_id
194
    MPI_INT64_T, // progeny_id
195
    MPI_DOUBLE,  // wgt_born
196
    MPI_DOUBLE,  // wgt_ww_born
197
    MPI_INT64_T  // n_split
198
  };
199

200
  MPI_Type_create_struct(14, blocks, disp, types, &mpi::source_site);
201
  MPI_Type_commit(&mpi::source_site);
202

203
  CollisionTrackSite bc;
204
  MPI_Aint dispc[16];
205
  MPI_Get_address(&bc.r, &dispc[0]);             // double
206
  MPI_Get_address(&bc.u, &dispc[1]);             // double
207
  MPI_Get_address(&bc.E, &dispc[2]);             // double
208
  MPI_Get_address(&bc.dE, &dispc[3]);            // double
209
  MPI_Get_address(&bc.time, &dispc[4]);          // double
210
  MPI_Get_address(&bc.wgt, &dispc[5]);           // double
211
  MPI_Get_address(&bc.event_mt, &dispc[6]);      // int
212
  MPI_Get_address(&bc.delayed_group, &dispc[7]); // int
213
  MPI_Get_address(&bc.cell_id, &dispc[8]);       // int
214
  MPI_Get_address(&bc.nuclide_id, &dispc[9]);    // int
215
  MPI_Get_address(&bc.material_id, &dispc[10]);  // int
216
  MPI_Get_address(&bc.universe_id, &dispc[11]);  // int
217
  MPI_Get_address(&bc.n_collision, &dispc[12]);  // int
218
  MPI_Get_address(&bc.particle, &dispc[13]);     // int
219
  MPI_Get_address(&bc.parent_id, &dispc[14]);    // int64_t
220
  MPI_Get_address(&bc.progeny_id, &dispc[15]);   // int64_t
221
  for (int i = 15; i >= 0; --i) {
222
    dispc[i] -= dispc[0];
223
  }
224

225
  int blocksc[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
226
  MPI_Datatype typesc[] = {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
227
    MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT,
228
    MPI_INT, MPI_INT, MPI_INT, MPI_INT64_T, MPI_INT64_T};
229

230
  MPI_Type_create_struct(
231
    16, blocksc, dispc, typesc, &mpi::collision_track_site);
232
  MPI_Type_commit(&mpi::collision_track_site);
233
}
234
#endif // OPENMC_MPI
235

236
int parse_command_line(int argc, char* argv[])
2,711✔
237
{
238
  int last_flag = 0;
2,711✔
239
  for (int i = 1; i < argc; ++i) {
3,099✔
240
    std::string arg {argv[i]};
392✔
241
    if (arg[0] == '-') {
392✔
242
      if (arg == "-p" || arg == "--plot") {
668!
243
        settings::run_mode = RunMode::PLOTTING;
44✔
244
        settings::check_overlaps = true;
44✔
245

246
      } else if (arg == "-n" || arg == "--particles") {
624!
247
        i += 1;
×
248
        settings::n_particles = std::stoll(argv[i]);
×
249

250
      } else if (arg == "-q" || arg == "--verbosity") {
624!
251
        i += 1;
×
252
        settings::verbosity = std::stoi(argv[i]);
×
253
        if (settings::verbosity > 10 || settings::verbosity < 1) {
×
254
          auto msg = fmt::format("Invalid verbosity: {}.", settings::verbosity);
×
255
          strcpy(openmc_err_msg, msg.c_str());
×
256
          return OPENMC_E_INVALID_ARGUMENT;
×
257
        }
×
258

259
      } else if (arg == "-e" || arg == "--event") {
624!
UNCOV
260
        settings::event_based = true;
×
261
      } else if (arg == "-r" || arg == "--restart") {
588!
262
        i += 1;
36✔
263
        // Check what type of file this is
264
        hid_t file_id = file_open(argv[i], 'r', true);
36✔
265
        std::string filetype;
36✔
266
        read_attribute(file_id, "filetype", filetype);
36✔
267
        file_close(file_id);
36✔
268

269
        // Set path and flag for type of run
270
        if (filetype == "statepoint") {
36✔
271
          settings::path_statepoint = argv[i];
20✔
272
          settings::path_statepoint_c = settings::path_statepoint.c_str();
20✔
273
          settings::restart_run = true;
20✔
274
        } else if (filetype == "particle restart") {
16!
275
          settings::path_particle_restart = argv[i];
16✔
276
          settings::particle_restart_run = true;
16✔
277
        } else {
278
          auto msg =
×
279
            fmt::format("Unrecognized file after restart flag: {}.", filetype);
×
280
          strcpy(openmc_err_msg, msg.c_str());
×
281
          return OPENMC_E_INVALID_ARGUMENT;
×
282
        }
×
283

284
        // If its a restart run check for additional source file
285
        if (settings::restart_run && i + 1 < argc) {
36!
286
          // Check if it has extension we can read
287
          if (ends_with(argv[i + 1], ".h5")) {
×
288

289
            // Check file type is a source file
290
            file_id = file_open(argv[i + 1], 'r', true);
×
291
            read_attribute(file_id, "filetype", filetype);
×
292
            file_close(file_id);
×
293
            if (filetype != "source") {
×
294
              std::string msg {
×
295
                "Second file after restart flag must be a source file"};
×
296
              strcpy(openmc_err_msg, msg.c_str());
×
297
              return OPENMC_E_INVALID_ARGUMENT;
×
298
            }
×
299

300
            // It is a source file
301
            settings::path_sourcepoint = argv[i + 1];
36!
302
            i += 1;
303

304
          } else {
305
            // Source is in statepoint file
306
            settings::path_sourcepoint = settings::path_statepoint;
×
307
          }
308

309
        } else {
310
          // Source is assumed to be in statepoint file
311
          settings::path_sourcepoint = settings::path_statepoint;
72✔
312
        }
313

314
      } else if (arg == "-g" || arg == "--geometry-debug") {
588!
315
        settings::check_overlaps = true;
×
316
      } else if (arg == "-c" || arg == "--volume") {
336✔
317
        settings::run_mode = RunMode::VOLUME;
248✔
318
      } else if (arg == "-s" || arg == "--threads") {
48!
319
        // Read number of threads
320
        i += 1;
8✔
321

322
#ifdef _OPENMP
323
        // Read and set number of OpenMP threads
324
        int n_threads = std::stoi(argv[i]);
4✔
325
        if (n_threads < 1) {
2!
326
          std::string msg {"Number of threads must be positive."};
×
327
          strcpy(openmc_err_msg, msg.c_str());
328
          return OPENMC_E_INVALID_ARGUMENT;
329
        }
330
        omp_set_num_threads(n_threads);
2✔
331
#else
332
        if (mpi::master) {
6!
333
          warning("Ignoring number of threads specified on command line.");
12✔
334
        }
335
#endif
336

337
      } else if (arg == "-?" || arg == "-h" || arg == "--help") {
60!
338
        print_usage();
×
339
        return OPENMC_E_UNASSIGNED;
×
340

341
      } else if (arg == "-v" || arg == "--version") {
36!
342
        print_version();
4✔
343
        print_build_info();
4✔
344
        return OPENMC_E_UNASSIGNED;
4✔
345

346
      } else if (arg == "-t" || arg == "--track") {
16!
347
        settings::write_all_tracks = true;
16✔
348

349
      } else {
350
        fmt::print(stderr, "Unknown option: {}\n", argv[i]);
×
351
        print_usage();
×
352
        return OPENMC_E_UNASSIGNED;
×
353
      }
354

355
      last_flag = i;
356
    }
357
  }
392✔
358

359
  // Determine directory where XML input files are
360
  if (argc > 1 && last_flag < argc - 1) {
2,707✔
361
    settings::path_input = std::string(argv[last_flag + 1]);
36✔
362

363
    // check that the path is either a valid directory or file
364
    if (!dir_exists(settings::path_input) &&
68✔
365
        !file_exists(settings::path_input)) {
32✔
366
      fatal_error(fmt::format(
12✔
367
        "The path specified to the OpenMC executable '{}' does not exist.",
368
        settings::path_input));
369
    }
370

371
    // Add slash at end of directory if it isn't there
372
    if (!ends_with(settings::path_input, "/") &&
72!
373
        dir_exists(settings::path_input)) {
24✔
374
      settings::path_input += "/";
2,699✔
375
    }
376
  }
377

378
  return 0;
379
}
380

381
// TODO: Pulse-height tallies require per-history scoring across the full
382
// particle tree (parent + all descendants). The shared secondary bank
383
// transports each secondary as an independent Particle, breaking this
384
// assumption. A proper fix would defer pulse-height scoring: save
385
// (root_source_id, cell, pht_storage) per particle, then aggregate by
386
// root_source_id after all secondary generations complete before scoring
387
// into the histogram. For now, disable shared secondary when pulse-height
388
// tallies are present.
389
static void check_pulse_height_compatibility()
2,659✔
390
{
391
  if (settings::use_shared_secondary_bank) {
2,659✔
392
    for (const auto& t : model::tallies) {
148✔
393
      if (t->type_ == TallyType::PULSE_HEIGHT) {
84✔
394
        settings::use_shared_secondary_bank = false;
8✔
395
        warning("Pulse-height tallies are not yet compatible with the shared "
8✔
396
                "secondary bank. Disabling shared secondary bank.");
397
        break;
8✔
398
      }
399
    }
400
  }
401
}
2,659✔
402

403
bool read_model_xml()
2,695✔
404
{
405
  std::string model_filename = settings::path_input;
2,695✔
406

407
  // if the current filename is a directory, append the default model filename
408
  if (model_filename.empty() || dir_exists(model_filename))
2,695✔
409
    model_filename += "model.xml";
2,675✔
410

411
  // if this file doesn't exist, stop here
412
  if (!file_exists(model_filename))
2,695✔
413
    return false;
414

415
  // try to process the path input as an XML file
416
  pugi::xml_document doc;
2,267✔
417
  if (!doc.load_file(model_filename.c_str())) {
2,267!
418
    fatal_error(fmt::format(
×
419
      "Error reading from single XML input file '{}'", model_filename));
420
  }
421

422
  pugi::xml_node root = doc.document_element();
2,267✔
423

424
  // Read settings
425
  if (!check_for_node(root, "settings")) {
2,267!
426
    fatal_error("No <settings> node present in the model.xml file.");
×
427
  }
428
  auto settings_root = root.child("settings");
2,267✔
429

430
  // Verbosity
431
  if (check_for_node(settings_root, "verbosity") && settings::verbosity == -1) {
2,267!
432
    settings::verbosity = std::stoi(get_node_value(settings_root, "verbosity"));
16✔
433
  } else if (settings::verbosity == -1) {
2,259!
434
    settings::verbosity = 7;
2,259✔
435
  }
436

437
  // To this point, we haven't displayed any output since we didn't know what
438
  // the verbosity is. Now that we checked for it, show the title if necessary
439
  if (mpi::master) {
2,267!
440
    if (settings::verbosity >= 2)
2,267✔
441
      title();
2,259✔
442
  }
443

444
  write_message(
2,267✔
445
    fmt::format("Reading model XML file '{}' ...", model_filename), 5);
2,267✔
446

447
  read_settings_xml(settings_root);
2,267✔
448

449
  // If other XML files are present, display warning
450
  // that they will be ignored
451
  auto other_inputs = {"materials.xml", "geometry.xml", "settings.xml",
2,243✔
452
    "tallies.xml", "plots.xml"};
2,243✔
453
  for (const auto& input : other_inputs) {
13,186✔
454
    if (file_exists(settings::path_input + input)) {
11,007✔
455
      warning((fmt::format("Other XML file input(s) are present. These files "
64✔
456
                           "may be ignored in favor of the {} file.",
457
        model_filename)));
458
      break;
64✔
459
    }
460
  }
461

462
  // Read data from chain file
463
  read_chain_file_xml();
2,243✔
464

465
  // Read materials and cross sections
466
  if (!check_for_node(root, "materials")) {
2,243!
467
    fatal_error(fmt::format(
×
468
      "No <materials> node present in the {} file.", model_filename));
469
  }
470

471
  if (settings::run_mode != RunMode::PLOTTING) {
2,243✔
472
    read_cross_sections_xml(root.child("materials"));
2,211✔
473
  }
474
  read_materials_xml(root.child("materials"));
2,243✔
475

476
  // Read geometry
477
  if (!check_for_node(root, "geometry")) {
2,243!
478
    fatal_error(fmt::format(
×
479
      "No <geometry> node present in the {} file.", model_filename));
480
  }
481
  read_geometry_xml(root.child("geometry"));
2,243✔
482

483
  // Final geometry setup and assign temperatures
484
  finalize_geometry();
2,243✔
485

486
  // Finalize cross sections having assigned temperatures
487
  finalize_cross_sections();
2,243✔
488

489
  // Compute cell density multipliers now that material densities
490
  // have been finalized (from geometry_aux.h)
491
  finalize_cell_densities();
2,243✔
492

493
  if (check_for_node(root, "tallies"))
2,243✔
494
    read_tallies_xml(root.child("tallies"));
1,328✔
495

496
  check_pulse_height_compatibility();
2,235✔
497

498
  // Initialize distribcell_filters
499
  prepare_distribcell();
2,235✔
500

501
  if (check_for_node(root, "plots")) {
2,235✔
502
    read_plots_xml(root.child("plots"));
124✔
503
  } else {
504
    // When no <plots> element is present in the model.xml file, check for a
505
    // regular plots.xml file
506
    std::string filename = settings::path_input + "plots.xml";
2,111✔
507
    if (file_exists(filename)) {
2,111!
508
      read_plots_xml();
×
509
    }
510
  }
2,111✔
511

512
  finalize_variance_reduction();
2,231✔
513

514
  return true;
2,231✔
515
}
4,890✔
516

517
void read_separate_xml_files()
428✔
518
{
519
  read_settings_xml();
428✔
520
  if (settings::run_mode != RunMode::PLOTTING) {
424✔
521
    read_cross_sections_xml();
412✔
522
  }
523

524
  // Read data from chain file
525
  read_chain_file_xml();
424✔
526

527
  read_materials_xml();
424✔
528
  read_geometry_xml();
424✔
529

530
  // Final geometry setup and assign temperatures
531
  finalize_geometry();
424✔
532

533
  // Finalize cross sections having assigned temperatures
534
  finalize_cross_sections();
424✔
535

536
  // Compute cell density multipliers now that material densities
537
  // have been finalized (from geometry_aux.h)
538
  finalize_cell_densities();
424✔
539

540
  read_tallies_xml();
424✔
541

542
  check_pulse_height_compatibility();
424✔
543

544
  // Initialize distribcell_filters
545
  prepare_distribcell();
424✔
546

547
  // Read the plots.xml regardless of plot mode in case plots are requested
548
  // via the API
549
  read_plots_xml();
424✔
550

551
  finalize_variance_reduction();
424✔
552
}
424✔
553

554
void initial_output()
2,655✔
555
{
556
  // write initial output
557
  if (settings::run_mode == RunMode::PLOTTING) {
2,655✔
558
    // Read plots.xml if it exists
559
    if (mpi::master && settings::verbosity >= 5)
40!
560
      print_plot();
32✔
561

562
  } else {
563
    // Write summary information
564
    if (mpi::master && settings::output_summary)
2,615!
565
      write_summary();
2,491✔
566

567
    // Warn if overlap checking is on
568
    if (mpi::master && settings::check_overlaps) {
2,615!
569
      warning("Cell overlap checking is ON.");
×
570
    }
571
  }
572
}
2,655✔
573

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