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

openmc-dev / openmc / 18538152141

15 Oct 2025 06:02PM UTC coverage: 81.97% (-3.2%) from 85.194%
18538152141

Pull #3417

github

web-flow
Merge 4604e1321 into e9077b137
Pull Request #3417: Addition of a collision tracking feature

16794 of 23357 branches covered (71.9%)

Branch coverage included in aggregate %.

480 of 522 new or added lines in 13 files covered. (91.95%)

457 existing lines in 53 files now uncovered.

54128 of 63165 relevant lines covered (85.69%)

42776927.64 hits per line

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

78.45
/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)
7,728✔
46
{
47
  using namespace openmc;
48

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

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

62
  // Parse command-line arguments
63
  int err = parse_command_line(argc, argv);
7,728✔
64
  if (err)
7,695✔
65
    return err;
11✔
66

67
#ifdef OPENMC_LIBMESH_ENABLED
68
  const int n_threads = num_threads();
1,417✔
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()) {
1,417!
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);
863✔
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);
554✔
83
#endif
84

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

88
#endif
89

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

94
#ifdef _OPENMP
95
  // If OMP_SCHEDULE is not set, default to a static schedule
96
  char* envvar = std::getenv("OMP_SCHEDULE");
4,251✔
97
  if (!envvar) {
4,251!
98
    omp_set_schedule(omp_sched_static, 0);
4,251✔
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);
7,684✔
105
  openmc::openmc_set_stride(DEFAULT_STRIDE);
7,684✔
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);
7,684✔
113
  if (std::setlocale(LC_ALL, "C") == NULL) {
7,684!
114
    fatal_error("Cannot set locale to C.");
×
115
  }
116

117
  // Read XML input files
118
  if (!read_model_xml())
7,684✔
119
    read_separate_xml_files();
1,637✔
120

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

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

129
  // Check for particle restart run
130
  if (settings::particle_restart_run)
7,591✔
131
    settings::run_mode = RunMode::PARTICLE;
43✔
132

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

137
  return 0;
7,591✔
138
}
7,591✔
139

140
namespace openmc {
141

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

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

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

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

176
  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
4,372✔
177
  MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
4,372✔
178
    MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
179
  MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site);
4,372✔
180
  MPI_Type_commit(&mpi::source_site);
4,372✔
181

182
  CollisionTrackSite bc;
4,372✔
183
  MPI_Aint dispc[16];
184
  MPI_Get_address(&bc.r, &dispc[0]);             // double
4,372✔
185
  MPI_Get_address(&bc.u, &dispc[1]);             // double
4,372✔
186
  MPI_Get_address(&bc.E, &dispc[2]);             // double
4,372✔
187
  MPI_Get_address(&bc.dE, &dispc[3]);            // double
4,372✔
188
  MPI_Get_address(&bc.time, &dispc[4]);          // double
4,372✔
189
  MPI_Get_address(&bc.wgt, &dispc[5]);           // double
4,372✔
190
  MPI_Get_address(&bc.event_mt, &dispc[6]);      // int
4,372✔
191
  MPI_Get_address(&bc.delayed_group, &dispc[7]); // int
4,372✔
192
  MPI_Get_address(&bc.cell_id, &dispc[8]);       // int
4,372✔
193
  MPI_Get_address(&bc.nuclide_id, &dispc[9]);    // int
4,372✔
194
  MPI_Get_address(&bc.material_id, &dispc[10]);  // int
4,372✔
195
  MPI_Get_address(&bc.universe_id, &dispc[11]);  // int
4,372✔
196
  MPI_Get_address(&bc.n_collision, &dispc[12]);  // int
4,372✔
197
  MPI_Get_address(&bc.particle, &dispc[13]);     // int
4,372✔
198
  MPI_Get_address(&bc.parent_id, &dispc[14]);    // int64_t
4,372✔
199
  MPI_Get_address(&bc.progeny_id, &dispc[15]);   // int64_t
4,372✔
200
  for (int i = 15; i >= 0; --i) {
74,324✔
201
    dispc[i] -= dispc[0];
69,952✔
202
  }
203

204
  int blocksc[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
4,372✔
205
  MPI_Datatype typesc[] = {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
4,372✔
206
    MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT,
207
    MPI_INT, MPI_INT, MPI_INT, MPI_INT64_T, MPI_INT64_T};
208

209
  MPI_Type_create_struct(
4,372✔
210
    16, blocksc, dispc, typesc, &mpi::collision_track_site);
211
  MPI_Type_commit(&mpi::collision_track_site);
4,372✔
212
}
4,372✔
213
#endif // OPENMC_MPI
214

215
int parse_command_line(int argc, char* argv[])
7,728✔
216
{
217
  int last_flag = 0;
7,728✔
218
  for (int i = 1; i < argc; ++i) {
8,787✔
219
    std::string arg {argv[i]};
1,070✔
220
    if (arg[0] == '-') {
1,070✔
221
      if (arg == "-p" || arg == "--plot") {
962!
222
        settings::run_mode = RunMode::PLOTTING;
276✔
223
        settings::check_overlaps = true;
276✔
224

225
      } else if (arg == "-n" || arg == "--particles") {
686!
226
        i += 1;
×
227
        settings::n_particles = std::stoll(argv[i]);
×
228

229
      } else if (arg == "-e" || arg == "--event") {
686!
230
        settings::event_based = true;
180✔
231
      } else if (arg == "-r" || arg == "--restart") {
506!
232
        i += 1;
111✔
233
        // Check what type of file this is
234
        hid_t file_id = file_open(argv[i], 'r', true);
111✔
235
        std::string filetype;
111✔
236
        read_attribute(file_id, "filetype", filetype);
111✔
237
        file_close(file_id);
111✔
238

239
        // Set path and flag for type of run
240
        if (filetype == "statepoint") {
111✔
241
          settings::path_statepoint = argv[i];
68✔
242
          settings::path_statepoint_c = settings::path_statepoint.c_str();
68✔
243
          settings::restart_run = true;
68✔
244
        } else if (filetype == "particle restart") {
43!
245
          settings::path_particle_restart = argv[i];
43✔
246
          settings::particle_restart_run = true;
43✔
247
        } else {
248
          auto msg =
249
            fmt::format("Unrecognized file after restart flag: {}.", filetype);
×
250
          strcpy(openmc_err_msg, msg.c_str());
×
251
          return OPENMC_E_INVALID_ARGUMENT;
×
UNCOV
252
        }
×
253

254
        // If its a restart run check for additional source file
255
        if (settings::restart_run && i + 1 < argc) {
111!
256
          // Check if it has extension we can read
257
          if (ends_with(argv[i + 1], ".h5")) {
×
258

259
            // Check file type is a source file
260
            file_id = file_open(argv[i + 1], 'r', true);
×
261
            read_attribute(file_id, "filetype", filetype);
×
262
            file_close(file_id);
×
263
            if (filetype != "source") {
×
264
              std::string msg {
265
                "Second file after restart flag must be a source file"};
×
266
              strcpy(openmc_err_msg, msg.c_str());
×
267
              return OPENMC_E_INVALID_ARGUMENT;
×
UNCOV
268
            }
×
269

270
            // It is a source file
271
            settings::path_sourcepoint = argv[i + 1];
×
272
            i += 1;
×
273

274
          } else {
275
            // Source is in statepoint file
276
            settings::path_sourcepoint = settings::path_statepoint;
×
277
          }
278

279
        } else {
×
280
          // Source is assumed to be in statepoint file
281
          settings::path_sourcepoint = settings::path_statepoint;
111✔
282
        }
283

284
      } else if (arg == "-g" || arg == "--geometry-debug") {
506!
285
        settings::check_overlaps = true;
×
286
      } else if (arg == "-c" || arg == "--volume") {
395✔
287
        settings::run_mode = RunMode::VOLUME;
300✔
288
      } else if (arg == "-s" || arg == "--threads") {
95!
289
        // Read number of threads
290
        i += 1;
25✔
291

292
#ifdef _OPENMP
293
        // Read and set number of OpenMP threads
294
        int n_threads = std::stoi(argv[i]);
13✔
295
        if (n_threads < 1) {
13!
296
          std::string msg {"Number of threads must be positive."};
×
297
          strcpy(openmc_err_msg, msg.c_str());
298
          return OPENMC_E_INVALID_ARGUMENT;
299
        }
300
        omp_set_num_threads(n_threads);
13✔
301
#else
302
        if (mpi::master) {
12✔
303
          warning("Ignoring number of threads specified on command line.");
10✔
304
        }
305
#endif
306

307
      } else if (arg == "-?" || arg == "-h" || arg == "--help") {
70!
308
        print_usage();
×
309
        return OPENMC_E_UNASSIGNED;
×
310

311
      } else if (arg == "-v" || arg == "--version") {
70!
312
        print_version();
11✔
313
        print_build_info();
11✔
314
        return OPENMC_E_UNASSIGNED;
11✔
315

316
      } else if (arg == "-t" || arg == "--track") {
59!
317
        settings::write_all_tracks = true;
59✔
318

319
      } else {
320
        fmt::print(stderr, "Unknown option: {}\n", argv[i]);
×
321
        print_usage();
×
322
        return OPENMC_E_UNASSIGNED;
×
323
      }
324

325
      last_flag = i;
951✔
326
    }
327
  }
1,070✔
328

329
  // Determine directory where XML input files are
330
  if (argc > 1 && last_flag < argc - 1) {
7,717✔
331
    settings::path_input = std::string(argv[last_flag + 1]);
108✔
332

333
    // check that the path is either a valid directory or file
334
    if (!dir_exists(settings::path_input) &&
204✔
335
        !file_exists(settings::path_input)) {
96✔
336
      fatal_error(fmt::format(
33✔
337
        "The path specified to the OpenMC executable '{}' does not exist.",
338
        settings::path_input));
339
    }
340

341
    // Add slash at end of directory if it isn't there
342
    if (!ends_with(settings::path_input, "/") &&
150!
343
        dir_exists(settings::path_input)) {
75✔
344
      settings::path_input += "/";
12✔
345
    }
346
  }
347

348
  return 0;
7,684✔
349
}
350

351
bool read_model_xml()
7,684✔
352
{
353
  std::string model_filename = settings::path_input;
7,684✔
354

355
  // if the current filename is a directory, append the default model filename
356
  if (model_filename.empty() || dir_exists(model_filename))
7,684✔
357
    model_filename += "model.xml";
7,621✔
358

359
  // if this file doesn't exist, stop here
360
  if (!file_exists(model_filename))
7,684✔
361
    return false;
1,637✔
362

363
  // try to process the path input as an XML file
364
  pugi::xml_document doc;
6,047✔
365
  if (!doc.load_file(model_filename.c_str())) {
6,047!
366
    fatal_error(fmt::format(
×
367
      "Error reading from single XML input file '{}'", model_filename));
368
  }
369

370
  pugi::xml_node root = doc.document_element();
6,047✔
371

372
  // Read settings
373
  if (!check_for_node(root, "settings")) {
6,047!
374
    fatal_error("No <settings> node present in the model.xml file.");
×
375
  }
376
  auto settings_root = root.child("settings");
6,047✔
377

378
  // Verbosity
379
  if (check_for_node(settings_root, "verbosity")) {
6,047✔
380
    settings::verbosity = std::stoi(get_node_value(settings_root, "verbosity"));
32✔
381
  }
382

383
  // To this point, we haven't displayed any output since we didn't know what
384
  // the verbosity is. Now that we checked for it, show the title if necessary
385
  if (mpi::master) {
6,047✔
386
    if (settings::verbosity >= 2)
5,237✔
387
      title();
5,215✔
388
  }
389

390
  write_message(
6,047✔
391
    fmt::format("Reading model XML file '{}' ...", model_filename), 5);
10,992✔
392

393
  read_settings_xml(settings_root);
6,047✔
394

395
  // If other XML files are present, display warning
396
  // that they will be ignored
397
  auto other_inputs = {"materials.xml", "geometry.xml", "settings.xml",
5,993✔
398
    "tallies.xml", "plots.xml"};
5,993✔
399
  for (const auto& input : other_inputs) {
34,965✔
400
    if (file_exists(settings::path_input + input)) {
29,192✔
401
      warning((fmt::format("Other XML file input(s) are present. These files "
220✔
402
                           "may be ignored in favor of the {} file.",
403
        model_filename)));
404
      break;
220✔
405
    }
406
  }
407

408
  // Read data from chain file
409
  read_chain_file_xml();
5,993✔
410

411
  // Read materials and cross sections
412
  if (!check_for_node(root, "materials")) {
5,993!
413
    fatal_error(fmt::format(
×
414
      "No <materials> node present in the {} file.", model_filename));
415
  }
416

417
  if (settings::run_mode != RunMode::PLOTTING) {
5,993✔
418
    read_cross_sections_xml(root.child("materials"));
5,750✔
419
  }
420
  read_materials_xml(root.child("materials"));
5,993✔
421

422
  // Read geometry
423
  if (!check_for_node(root, "geometry")) {
5,993!
424
    fatal_error(fmt::format(
×
425
      "No <geometry> node present in the {} file.", model_filename));
426
  }
427
  read_geometry_xml(root.child("geometry"));
5,993✔
428

429
  // Final geometry setup and assign temperatures
430
  finalize_geometry();
5,991✔
431

432
  // Finalize cross sections having assigned temperatures
433
  finalize_cross_sections();
5,991✔
434

435
  // Compute cell density multipliers now that material densities
436
  // have been finalized (from geometry_aux.h)
437
  finalize_cell_densities();
5,991✔
438

439
  if (check_for_node(root, "tallies"))
5,991✔
440
    read_tallies_xml(root.child("tallies"));
3,573✔
441

442
  // Initialize distribcell_filters
443
  prepare_distribcell();
5,973✔
444

445
  if (check_for_node(root, "plots")) {
5,973✔
446
    read_plots_xml(root.child("plots"));
472✔
447
  } else {
448
    // When no <plots> element is present in the model.xml file, check for a
449
    // regular plots.xml file
450
    std::string filename = settings::path_input + "plots.xml";
5,501✔
451
    if (file_exists(filename)) {
5,501!
452
      read_plots_xml();
×
453
    }
454
  }
5,501✔
455

456
  finalize_variance_reduction();
5,964✔
457

458
  return true;
5,964✔
459
}
7,601✔
460

461
void read_separate_xml_files()
1,637✔
462
{
463
  read_settings_xml();
1,637✔
464
  if (settings::run_mode != RunMode::PLOTTING) {
1,627✔
465
    read_cross_sections_xml();
1,594✔
466
  }
467

468
  // Read data from chain file
469
  read_chain_file_xml();
1,627✔
470

471
  read_materials_xml();
1,627✔
472
  read_geometry_xml();
1,627✔
473

474
  // Final geometry setup and assign temperatures
475
  finalize_geometry();
1,627✔
476

477
  // Finalize cross sections having assigned temperatures
478
  finalize_cross_sections();
1,627✔
479

480
  // Compute cell density multipliers now that material densities
481
  // have been finalized (from geometry_aux.h)
482
  finalize_cell_densities();
1,627✔
483

484
  read_tallies_xml();
1,627✔
485

486
  // Initialize distribcell_filters
487
  prepare_distribcell();
1,627✔
488

489
  // Read the plots.xml regardless of plot mode in case plots are requested
490
  // via the API
491
  read_plots_xml();
1,627✔
492

493
  finalize_variance_reduction();
1,627✔
494
}
1,627✔
495

496
void initial_output()
7,591✔
497
{
498
  // write initial output
499
  if (settings::run_mode == RunMode::PLOTTING) {
7,591✔
500
    // Read plots.xml if it exists
501
    if (mpi::master && settings::verbosity >= 5)
267!
502
      print_plot();
267✔
503

504
  } else {
505
    // Write summary information
506
    if (mpi::master && settings::output_summary)
7,324✔
507
      write_summary();
5,933✔
508

509
    // Warn if overlap checking is on
510
    if (mpi::master && settings::check_overlaps) {
7,324!
511
      warning("Cell overlap checking is ON.");
×
512
    }
513
  }
514
}
7,591✔
515

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

© 2025 Coveralls, Inc