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

openmc-dev / openmc / 20278909785

16 Dec 2025 06:41PM UTC coverage: 81.834% (-0.09%) from 81.92%
20278909785

Pull #3493

github

web-flow
Merge 9dc7c7a58 into bbfa18d72
Pull Request #3493: Implement vector fitting to replace external `vectfit` package

17020 of 23572 branches covered (72.2%)

Branch coverage included in aggregate %.

188 of 207 new or added lines in 3 files covered. (90.82%)

3101 existing lines in 56 files now uncovered.

54979 of 64410 relevant lines covered (85.36%)

41388074.54 hits per line

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

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

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

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

62
  // Parse command-line arguments
63
  int err = parse_command_line(argc, argv);
6,888✔
64
  if (err)
6,858✔
65
    return err;
10✔
66

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

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

88
#endif
89

90
  // Start total and initialization timer
91
  simulation::time_total.start();
6,848✔
92
  simulation::time_initialize.start();
6,848✔
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);
6,848✔
105
  openmc::openmc_set_stride(DEFAULT_STRIDE);
6,848✔
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);
6,848✔
113
  if (std::setlocale(LC_ALL, "C") == NULL) {
6,848!
114
    fatal_error("Cannot set locale to C.");
×
115
  }
116

117
  // Read XML input files
118
  if (!read_model_xml())
6,848✔
119
    read_separate_xml_files();
1,234✔
120

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

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

129
  // Check for particle restart run
130
  if (settings::particle_restart_run)
6,765✔
131
    settings::run_mode = RunMode::PARTICLE;
38✔
132

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

137
  return 0;
6,765✔
138
}
6,765✔
139

140
namespace openmc {
141

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

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

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

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

176
  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
3,310✔
177
  MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
3,310✔
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);
3,310✔
180
  MPI_Type_commit(&mpi::source_site);
3,310✔
181

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

204
  int blocksc[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
3,310✔
205
  MPI_Datatype typesc[] = {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
3,310✔
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(
3,310✔
210
    16, blocksc, dispc, typesc, &mpi::collision_track_site);
211
  MPI_Type_commit(&mpi::collision_track_site);
3,310✔
212
}
3,310✔
213
#endif // OPENMC_MPI
214

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

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

229
      } else if (arg == "-q" || arg == "--verbosity") {
651!
230
        i += 1;
×
231
        settings::verbosity = std::stoi(argv[i]);
×
232
        if (settings::verbosity > 10 || settings::verbosity < 1) {
×
UNCOV
233
          auto msg = fmt::format("Invalid verbosity: {}.", settings::verbosity);
×
234
          strcpy(openmc_err_msg, msg.c_str());
×
235
          return OPENMC_E_INVALID_ARGUMENT;
×
236
        }
×
237

238
      } else if (arg == "-e" || arg == "--event") {
651!
239
        settings::event_based = true;
187✔
240
      } else if (arg == "-r" || arg == "--restart") {
464!
241
        i += 1;
96✔
242
        // Check what type of file this is
243
        hid_t file_id = file_open(argv[i], 'r', true);
96✔
244
        std::string filetype;
96✔
245
        read_attribute(file_id, "filetype", filetype);
96✔
246
        file_close(file_id);
96✔
247

248
        // Set path and flag for type of run
249
        if (filetype == "statepoint") {
96✔
250
          settings::path_statepoint = argv[i];
58✔
251
          settings::path_statepoint_c = settings::path_statepoint.c_str();
58✔
252
          settings::restart_run = true;
58✔
253
        } else if (filetype == "particle restart") {
38!
254
          settings::path_particle_restart = argv[i];
38✔
255
          settings::particle_restart_run = true;
38✔
256
        } else {
257
          auto msg =
UNCOV
258
            fmt::format("Unrecognized file after restart flag: {}.", filetype);
×
UNCOV
259
          strcpy(openmc_err_msg, msg.c_str());
×
UNCOV
260
          return OPENMC_E_INVALID_ARGUMENT;
×
UNCOV
261
        }
×
262

263
        // If its a restart run check for additional source file
264
        if (settings::restart_run && i + 1 < argc) {
96!
265
          // Check if it has extension we can read
UNCOV
266
          if (ends_with(argv[i + 1], ".h5")) {
×
267

268
            // Check file type is a source file
UNCOV
269
            file_id = file_open(argv[i + 1], 'r', true);
×
UNCOV
270
            read_attribute(file_id, "filetype", filetype);
×
UNCOV
271
            file_close(file_id);
×
UNCOV
272
            if (filetype != "source") {
×
273
              std::string msg {
UNCOV
274
                "Second file after restart flag must be a source file"};
×
UNCOV
275
              strcpy(openmc_err_msg, msg.c_str());
×
UNCOV
276
              return OPENMC_E_INVALID_ARGUMENT;
×
277
            }
×
278

279
            // It is a source file
UNCOV
280
            settings::path_sourcepoint = argv[i + 1];
×
UNCOV
281
            i += 1;
×
282

283
          } else {
284
            // Source is in statepoint file
UNCOV
285
            settings::path_sourcepoint = settings::path_statepoint;
×
286
          }
287

UNCOV
288
        } else {
×
289
          // Source is assumed to be in statepoint file
290
          settings::path_sourcepoint = settings::path_statepoint;
96✔
291
        }
292

293
      } else if (arg == "-g" || arg == "--geometry-debug") {
464!
UNCOV
294
        settings::check_overlaps = true;
×
295
      } else if (arg == "-c" || arg == "--volume") {
368✔
296
        settings::run_mode = RunMode::VOLUME;
284✔
297
      } else if (arg == "-s" || arg == "--threads") {
84!
298
        // Read number of threads
299
        i += 1;
22✔
300

301
#ifdef _OPENMP
302
        // Read and set number of OpenMP threads
303
        int n_threads = std::stoi(argv[i]);
13✔
304
        if (n_threads < 1) {
13!
305
          std::string msg {"Number of threads must be positive."};
×
306
          strcpy(openmc_err_msg, msg.c_str());
307
          return OPENMC_E_INVALID_ARGUMENT;
308
        }
309
        omp_set_num_threads(n_threads);
13✔
310
#else
311
        if (mpi::master) {
9✔
312
          warning("Ignoring number of threads specified on command line.");
8✔
313
        }
314
#endif
315

316
      } else if (arg == "-?" || arg == "-h" || arg == "--help") {
62!
UNCOV
317
        print_usage();
×
UNCOV
318
        return OPENMC_E_UNASSIGNED;
×
319

320
      } else if (arg == "-v" || arg == "--version") {
62!
321
        print_version();
10✔
322
        print_build_info();
10✔
323
        return OPENMC_E_UNASSIGNED;
10✔
324

325
      } else if (arg == "-t" || arg == "--track") {
52!
326
        settings::write_all_tracks = true;
52✔
327

328
      } else {
UNCOV
329
        fmt::print(stderr, "Unknown option: {}\n", argv[i]);
×
UNCOV
330
        print_usage();
×
UNCOV
331
        return OPENMC_E_UNASSIGNED;
×
332
      }
333

334
      last_flag = i;
892✔
335
    }
336
  }
995✔
337

338
  // Determine directory where XML input files are
339
  if (argc > 1 && last_flag < argc - 1) {
6,878✔
340
    settings::path_input = std::string(argv[last_flag + 1]);
93✔
341

342
    // check that the path is either a valid directory or file
343
    if (!dir_exists(settings::path_input) &&
175✔
344
        !file_exists(settings::path_input)) {
82✔
345
      fatal_error(fmt::format(
30✔
346
        "The path specified to the OpenMC executable '{}' does not exist.",
347
        settings::path_input));
348
    }
349

350
    // Add slash at end of directory if it isn't there
351
    if (!ends_with(settings::path_input, "/") &&
126!
352
        dir_exists(settings::path_input)) {
63✔
353
      settings::path_input += "/";
11✔
354
    }
355
  }
356

357
  return 0;
6,848✔
358
}
359

360
bool read_model_xml()
6,848✔
361
{
362
  std::string model_filename = settings::path_input;
6,848✔
363

364
  // if the current filename is a directory, append the default model filename
365
  if (model_filename.empty() || dir_exists(model_filename))
6,848✔
366
    model_filename += "model.xml";
6,796✔
367

368
  // if this file doesn't exist, stop here
369
  if (!file_exists(model_filename))
6,848✔
370
    return false;
1,234✔
371

372
  // try to process the path input as an XML file
373
  pugi::xml_document doc;
5,614✔
374
  if (!doc.load_file(model_filename.c_str())) {
5,614!
UNCOV
375
    fatal_error(fmt::format(
×
376
      "Error reading from single XML input file '{}'", model_filename));
377
  }
378

379
  pugi::xml_node root = doc.document_element();
5,614✔
380

381
  // Read settings
382
  if (!check_for_node(root, "settings")) {
5,614!
UNCOV
383
    fatal_error("No <settings> node present in the model.xml file.");
×
384
  }
385
  auto settings_root = root.child("settings");
5,614✔
386

387
  // Verbosity
388
  if (check_for_node(settings_root, "verbosity") && settings::verbosity == -1) {
5,614!
389
    settings::verbosity = std::stoi(get_node_value(settings_root, "verbosity"));
28✔
390
  } else if (settings::verbosity == -1) {
5,586!
391
    settings::verbosity = 7;
5,586✔
392
  }
393

394
  // To this point, we haven't displayed any output since we didn't know what
395
  // the verbosity is. Now that we checked for it, show the title if necessary
396
  if (mpi::master) {
5,614✔
397
    if (settings::verbosity >= 2)
4,926✔
398
      title();
4,906✔
399
  }
400

401
  write_message(
5,614✔
402
    fmt::format("Reading model XML file '{}' ...", model_filename), 5);
10,085✔
403

404
  read_settings_xml(settings_root);
5,614✔
405

406
  // If other XML files are present, display warning
407
  // that they will be ignored
408
  auto other_inputs = {"materials.xml", "geometry.xml", "settings.xml",
5,566✔
409
    "tallies.xml", "plots.xml"};
5,566✔
410
  for (const auto& input : other_inputs) {
32,875✔
411
    if (file_exists(settings::path_input + input)) {
27,432✔
412
      warning((fmt::format("Other XML file input(s) are present. These files "
123✔
413
                           "may be ignored in favor of the {} file.",
414
        model_filename)));
415
      break;
123✔
416
    }
417
  }
418

419
  // Read data from chain file
420
  read_chain_file_xml();
5,566✔
421

422
  // Read materials and cross sections
423
  if (!check_for_node(root, "materials")) {
5,566!
UNCOV
424
    fatal_error(fmt::format(
×
425
      "No <materials> node present in the {} file.", model_filename));
426
  }
427

428
  if (settings::run_mode != RunMode::PLOTTING) {
5,566✔
429
    read_cross_sections_xml(root.child("materials"));
5,345✔
430
  }
431
  read_materials_xml(root.child("materials"));
5,566✔
432

433
  // Read geometry
434
  if (!check_for_node(root, "geometry")) {
5,566!
UNCOV
435
    fatal_error(fmt::format(
×
436
      "No <geometry> node present in the {} file.", model_filename));
437
  }
438
  read_geometry_xml(root.child("geometry"));
5,566✔
439

440
  // Final geometry setup and assign temperatures
441
  finalize_geometry();
5,564✔
442

443
  // Finalize cross sections having assigned temperatures
444
  finalize_cross_sections();
5,564✔
445

446
  // Compute cell density multipliers now that material densities
447
  // have been finalized (from geometry_aux.h)
448
  finalize_cell_densities();
5,564✔
449

450
  if (check_for_node(root, "tallies"))
5,564✔
451
    read_tallies_xml(root.child("tallies"));
3,271✔
452

453
  // Initialize distribcell_filters
454
  prepare_distribcell();
5,548✔
455

456
  if (check_for_node(root, "plots")) {
5,548✔
457
    read_plots_xml(root.child("plots"));
425✔
458
  } else {
459
    // When no <plots> element is present in the model.xml file, check for a
460
    // regular plots.xml file
461
    std::string filename = settings::path_input + "plots.xml";
5,123✔
462
    if (file_exists(filename)) {
5,123!
UNCOV
463
      read_plots_xml();
×
464
    }
465
  }
5,123✔
466

467
  finalize_variance_reduction();
5,540✔
468

469
  return true;
5,540✔
470
}
6,774✔
471

472
void read_separate_xml_files()
1,234✔
473
{
474
  read_settings_xml();
1,234✔
475
  if (settings::run_mode != RunMode::PLOTTING) {
1,225✔
476
    read_cross_sections_xml();
1,195✔
477
  }
478

479
  // Read data from chain file
480
  read_chain_file_xml();
1,225✔
481

482
  read_materials_xml();
1,225✔
483
  read_geometry_xml();
1,225✔
484

485
  // Final geometry setup and assign temperatures
486
  finalize_geometry();
1,225✔
487

488
  // Finalize cross sections having assigned temperatures
489
  finalize_cross_sections();
1,225✔
490

491
  // Compute cell density multipliers now that material densities
492
  // have been finalized (from geometry_aux.h)
493
  finalize_cell_densities();
1,225✔
494

495
  read_tallies_xml();
1,225✔
496

497
  // Initialize distribcell_filters
498
  prepare_distribcell();
1,225✔
499

500
  // Read the plots.xml regardless of plot mode in case plots are requested
501
  // via the API
502
  read_plots_xml();
1,225✔
503

504
  finalize_variance_reduction();
1,225✔
505
}
1,225✔
506

507
void initial_output()
6,765✔
508
{
509
  // write initial output
510
  if (settings::run_mode == RunMode::PLOTTING) {
6,765✔
511
    // Read plots.xml if it exists
512
    if (mpi::master && settings::verbosity >= 5)
243!
513
      print_plot();
223✔
514

515
  } else {
516
    // Write summary information
517
    if (mpi::master && settings::output_summary)
6,522✔
518
      write_summary();
5,360✔
519

520
    // Warn if overlap checking is on
521
    if (mpi::master && settings::check_overlaps) {
6,522!
UNCOV
522
      warning("Cell overlap checking is ON.");
×
523
    }
524
  }
525
}
6,765✔
526

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