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

openmc-dev / openmc / 18656608910

20 Oct 2025 03:16PM UTC coverage: 81.844% (-3.4%) from 85.218%
18656608910

Pull #3454

github

web-flow
Merge 5eee478a5 into 055ea15a2
Pull Request #3454: Adding variance of variance and normality tests for tally statistics

16610 of 23118 branches covered (71.85%)

Branch coverage included in aggregate %.

188 of 312 new or added lines in 7 files covered. (60.26%)

2157 existing lines in 77 files now uncovered.

53896 of 63029 relevant lines covered (85.51%)

42554448.35 hits per line

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

76.72
/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,602✔
46
{
47
  using namespace openmc;
48

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

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

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

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

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

88
#endif
89

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

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

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

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

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

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

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

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

140
namespace openmc {
141

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

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

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

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

176
  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
4,299✔
177
  MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
4,299✔
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,299✔
180
  MPI_Type_commit(&mpi::source_site);
4,299✔
181
}
4,299✔
182
#endif // OPENMC_MPI
183

184
int parse_command_line(int argc, char* argv[])
7,602✔
185
{
186
  int last_flag = 0;
7,602✔
187
  for (int i = 1; i < argc; ++i) {
8,634✔
188
    std::string arg {argv[i]};
1,043✔
189
    if (arg[0] == '-') {
1,043✔
190
      if (arg == "-p" || arg == "--plot") {
933!
191
        settings::run_mode = RunMode::PLOTTING;
276✔
192
        settings::check_overlaps = true;
276✔
193

194
      } else if (arg == "-n" || arg == "--particles") {
657!
195
        i += 1;
×
196
        settings::n_particles = std::stoll(argv[i]);
×
197

198
      } else if (arg == "-e" || arg == "--event") {
657!
199
        settings::event_based = true;
174✔
200
      } else if (arg == "-r" || arg == "--restart") {
483!
201
        i += 1;
108✔
202
        // Check what type of file this is
203
        hid_t file_id = file_open(argv[i], 'r', true);
108✔
204
        std::string filetype;
108✔
205
        read_attribute(file_id, "filetype", filetype);
108✔
206
        file_close(file_id);
108✔
207

208
        // Set path and flag for type of run
209
        if (filetype == "statepoint") {
108✔
210
          settings::path_statepoint = argv[i];
65✔
211
          settings::path_statepoint_c = settings::path_statepoint.c_str();
65✔
212
          settings::restart_run = true;
65✔
213
        } else if (filetype == "particle restart") {
43!
214
          settings::path_particle_restart = argv[i];
43✔
215
          settings::particle_restart_run = true;
43✔
216
        } else {
217
          auto msg =
218
            fmt::format("Unrecognized file after restart flag: {}.", filetype);
×
219
          strcpy(openmc_err_msg, msg.c_str());
×
220
          return OPENMC_E_INVALID_ARGUMENT;
×
UNCOV
221
        }
×
222

223
        // If its a restart run check for additional source file
224
        if (settings::restart_run && i + 1 < argc) {
108!
225
          // Check if it has extension we can read
226
          if (ends_with(argv[i + 1], ".h5")) {
×
227

228
            // Check file type is a source file
229
            file_id = file_open(argv[i + 1], 'r', true);
×
230
            read_attribute(file_id, "filetype", filetype);
×
231
            file_close(file_id);
×
232
            if (filetype != "source") {
×
233
              std::string msg {
234
                "Second file after restart flag must be a source file"};
×
235
              strcpy(openmc_err_msg, msg.c_str());
×
236
              return OPENMC_E_INVALID_ARGUMENT;
×
UNCOV
237
            }
×
238

239
            // It is a source file
240
            settings::path_sourcepoint = argv[i + 1];
×
241
            i += 1;
×
242

243
          } else {
244
            // Source is in statepoint file
245
            settings::path_sourcepoint = settings::path_statepoint;
×
246
          }
247

248
        } else {
×
249
          // Source is assumed to be in statepoint file
250
          settings::path_sourcepoint = settings::path_statepoint;
108✔
251
        }
252

253
      } else if (arg == "-g" || arg == "--geometry-debug") {
483!
254
        settings::check_overlaps = true;
×
255
      } else if (arg == "-c" || arg == "--volume") {
375✔
256
        settings::run_mode = RunMode::VOLUME;
280✔
257
      } else if (arg == "-s" || arg == "--threads") {
95!
258
        // Read number of threads
259
        i += 1;
25✔
260

261
#ifdef _OPENMP
262
        // Read and set number of OpenMP threads
263
        int n_threads = std::stoi(argv[i]);
13✔
264
        if (n_threads < 1) {
13!
265
          std::string msg {"Number of threads must be positive."};
×
266
          strcpy(openmc_err_msg, msg.c_str());
267
          return OPENMC_E_INVALID_ARGUMENT;
268
        }
269
        omp_set_num_threads(n_threads);
13✔
270
#else
271
        if (mpi::master) {
12✔
272
          warning("Ignoring number of threads specified on command line.");
10✔
273
        }
274
#endif
275

276
      } else if (arg == "-?" || arg == "-h" || arg == "--help") {
70!
277
        print_usage();
×
278
        return OPENMC_E_UNASSIGNED;
×
279

280
      } else if (arg == "-v" || arg == "--version") {
70!
281
        print_version();
11✔
282
        print_build_info();
11✔
283
        return OPENMC_E_UNASSIGNED;
11✔
284

285
      } else if (arg == "-t" || arg == "--track") {
59!
286
        settings::write_all_tracks = true;
59✔
287

288
      } else {
289
        fmt::print(stderr, "Unknown option: {}\n", argv[i]);
×
290
        print_usage();
×
291
        return OPENMC_E_UNASSIGNED;
×
292
      }
293

294
      last_flag = i;
922✔
295
    }
296
  }
1,043✔
297

298
  // Determine directory where XML input files are
299
  if (argc > 1 && last_flag < argc - 1) {
7,591✔
300
    settings::path_input = std::string(argv[last_flag + 1]);
110✔
301

302
    // check that the path is either a valid directory or file
303
    if (!dir_exists(settings::path_input) &&
208✔
304
        !file_exists(settings::path_input)) {
98✔
305
      fatal_error(fmt::format(
33✔
306
        "The path specified to the OpenMC executable '{}' does not exist.",
307
        settings::path_input));
308
    }
309

310
    // Add slash at end of directory if it isn't there
311
    if (!ends_with(settings::path_input, "/") &&
154!
312
        dir_exists(settings::path_input)) {
77✔
313
      settings::path_input += "/";
12✔
314
    }
315
  }
316

317
  return 0;
7,558✔
318
}
319

320
bool read_model_xml()
7,558✔
321
{
322
  std::string model_filename = settings::path_input;
7,558✔
323

324
  // if the current filename is a directory, append the default model filename
325
  if (model_filename.empty() || dir_exists(model_filename))
7,558✔
326
    model_filename += "model.xml";
7,493✔
327

328
  // if this file doesn't exist, stop here
329
  if (!file_exists(model_filename))
7,558✔
330
    return false;
1,660✔
331

332
  // try to process the path input as an XML file
333
  pugi::xml_document doc;
5,898✔
334
  if (!doc.load_file(model_filename.c_str())) {
5,898!
335
    fatal_error(fmt::format(
×
336
      "Error reading from single XML input file '{}'", model_filename));
337
  }
338

339
  pugi::xml_node root = doc.document_element();
5,898✔
340

341
  // Read settings
342
  if (!check_for_node(root, "settings")) {
5,898!
343
    fatal_error("No <settings> node present in the model.xml file.");
×
344
  }
345
  auto settings_root = root.child("settings");
5,898✔
346

347
  // Verbosity
348
  if (check_for_node(settings_root, "verbosity")) {
5,898✔
349
    settings::verbosity = std::stoi(get_node_value(settings_root, "verbosity"));
32✔
350
  }
351

352
  // To this point, we haven't displayed any output since we didn't know what
353
  // the verbosity is. Now that we checked for it, show the title if necessary
354
  if (mpi::master) {
5,898✔
355
    if (settings::verbosity >= 2)
5,118✔
356
      title();
5,096✔
357
  }
358

359
  write_message(
5,898✔
360
    fmt::format("Reading model XML file '{}' ...", model_filename), 5);
10,696✔
361

362
  read_settings_xml(settings_root);
5,898✔
363

364
  // If other XML files are present, display warning
365
  // that they will be ignored
366
  auto other_inputs = {"materials.xml", "geometry.xml", "settings.xml",
5,844✔
367
    "tallies.xml", "plots.xml"};
5,844✔
368
  for (const auto& input : other_inputs) {
34,011✔
369
    if (file_exists(settings::path_input + input)) {
28,399✔
370
      warning((fmt::format("Other XML file input(s) are present. These files "
232✔
371
                           "may be ignored in favor of the {} file.",
372
        model_filename)));
373
      break;
232✔
374
    }
375
  }
376

377
  // Read data from chain file
378
  read_chain_file_xml();
5,844✔
379

380
  // Read materials and cross sections
381
  if (!check_for_node(root, "materials")) {
5,844!
382
    fatal_error(fmt::format(
×
383
      "No <materials> node present in the {} file.", model_filename));
384
  }
385

386
  if (settings::run_mode != RunMode::PLOTTING) {
5,844✔
387
    read_cross_sections_xml(root.child("materials"));
5,601✔
388
  }
389
  read_materials_xml(root.child("materials"));
5,844✔
390

391
  // Read geometry
392
  if (!check_for_node(root, "geometry")) {
5,844!
393
    fatal_error(fmt::format(
×
394
      "No <geometry> node present in the {} file.", model_filename));
395
  }
396
  read_geometry_xml(root.child("geometry"));
5,844✔
397

398
  // Final geometry setup and assign temperatures
399
  finalize_geometry();
5,842✔
400

401
  // Finalize cross sections having assigned temperatures
402
  finalize_cross_sections();
5,842✔
403

404
  // Compute cell density multipliers now that material densities
405
  // have been finalized (from geometry_aux.h)
406
  finalize_cell_densities();
5,842✔
407

408
  if (check_for_node(root, "tallies"))
5,842✔
409
    read_tallies_xml(root.child("tallies"));
3,576✔
410

411
  // Initialize distribcell_filters
412
  prepare_distribcell();
5,824✔
413

414
  if (check_for_node(root, "plots")) {
5,824✔
415
    read_plots_xml(root.child("plots"));
468✔
416
  } else {
417
    // When no <plots> element is present in the model.xml file, check for a
418
    // regular plots.xml file
419
    std::string filename = settings::path_input + "plots.xml";
5,356✔
420
    if (file_exists(filename)) {
5,356!
UNCOV
421
      read_plots_xml();
×
422
    }
423
  }
5,356✔
424

425
  finalize_variance_reduction();
5,815✔
426

427
  return true;
5,815✔
428
}
7,475✔
429

430
void read_separate_xml_files()
1,660✔
431
{
432
  read_settings_xml();
1,660✔
433
  if (settings::run_mode != RunMode::PLOTTING) {
1,650✔
434
    read_cross_sections_xml();
1,617✔
435
  }
436

437
  // Read data from chain file
438
  read_chain_file_xml();
1,650✔
439

440
  read_materials_xml();
1,650✔
441
  read_geometry_xml();
1,650✔
442

443
  // Final geometry setup and assign temperatures
444
  finalize_geometry();
1,650✔
445

446
  // Finalize cross sections having assigned temperatures
447
  finalize_cross_sections();
1,650✔
448

449
  // Compute cell density multipliers now that material densities
450
  // have been finalized (from geometry_aux.h)
451
  finalize_cell_densities();
1,650✔
452

453
  read_tallies_xml();
1,650✔
454

455
  // Initialize distribcell_filters
456
  prepare_distribcell();
1,650✔
457

458
  // Read the plots.xml regardless of plot mode in case plots are requested
459
  // via the API
460
  read_plots_xml();
1,650✔
461

462
  finalize_variance_reduction();
1,650✔
463
}
1,650✔
464

465
void initial_output()
7,465✔
466
{
467
  // write initial output
468
  if (settings::run_mode == RunMode::PLOTTING) {
7,465✔
469
    // Read plots.xml if it exists
470
    if (mpi::master && settings::verbosity >= 5)
267!
471
      print_plot();
267✔
472

473
  } else {
474
    // Write summary information
475
    if (mpi::master && settings::output_summary)
7,198✔
476
      write_summary();
5,822✔
477

478
    // Warn if overlap checking is on
479
    if (mpi::master && settings::check_overlaps) {
7,198!
UNCOV
480
      warning("Cell overlap checking is ON.");
×
481
    }
482
  }
483
}
7,465✔
484

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