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

openmc-dev / openmc / 20158005568

12 Dec 2025 06:07AM UTC coverage: 82.149% (-0.02%) from 82.166%
20158005568

Pull #3680

github

web-flow
Merge c09cb8305 into 5c4121efd
Pull Request #3680: Add a command-line argument for output verbosity

17011 of 23573 branches covered (72.16%)

Branch coverage included in aggregate %.

8 of 15 new or added lines in 3 files covered. (53.33%)

11 existing lines in 1 file now uncovered.

55098 of 64205 relevant lines covered (85.82%)

43482942.69 hits per line

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

86.53
/src/state_point.cpp
1
#include "openmc/state_point.h"
2

3
#include <algorithm>
4
#include <cstdint> // for int64_t
5
#include <string>
6

7
#include "xtensor/xbuilder.hpp" // for empty_like
8
#include "xtensor/xview.hpp"
9
#include <fmt/core.h>
10

11
#include "openmc/bank.h"
12
#include "openmc/bank_io.h"
13
#include "openmc/capi.h"
14
#include "openmc/constants.h"
15
#include "openmc/eigenvalue.h"
16
#include "openmc/error.h"
17
#include "openmc/file_utils.h"
18
#include "openmc/hdf5_interface.h"
19
#include "openmc/mcpl_interface.h"
20
#include "openmc/mesh.h"
21
#include "openmc/message_passing.h"
22
#include "openmc/mgxs_interface.h"
23
#include "openmc/nuclide.h"
24
#include "openmc/output.h"
25
#include "openmc/settings.h"
26
#include "openmc/simulation.h"
27
#include "openmc/tallies/derivative.h"
28
#include "openmc/tallies/filter.h"
29
#include "openmc/tallies/filter_mesh.h"
30
#include "openmc/tallies/tally.h"
31
#include "openmc/timer.h"
32
#include "openmc/vector.h"
33

34
namespace openmc {
35

36
extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
7,679✔
37
{
38
  simulation::time_statepoint.start();
7,679✔
39

40
  // If a nullptr is passed in, we assume that the user
41
  // wants a default name for this, of the form like output/statepoint.20.h5
42
  std::string filename_;
7,679✔
43
  if (filename) {
7,679✔
44
    filename_ = filename;
797✔
45
  } else {
46
    // Determine width for zero padding
47
    int w = std::to_string(settings::n_max_batches).size();
6,882✔
48

49
    // Set filename for state point
50
    filename_ = fmt::format("{0}statepoint.{1:0{2}}.h5", settings::path_output,
12,481✔
51
      simulation::current_batch, w);
6,882✔
52
  }
53

54
  // If a file name was specified, ensure it has .h5 file extension
55
  const auto extension = get_file_extension(filename_);
7,679✔
56
  if (extension != "h5") {
7,679!
57
    warning("openmc_statepoint_write was passed a file extension differing "
×
58
            "from .h5, but an hdf5 file will be written.");
59
  }
60

61
  // Determine whether or not to write the source bank
62
  bool write_source_ = write_source ? *write_source : true;
7,679!
63

64
  // Write message
65
  write_message("Creating state point " + filename_ + "...", 5);
7,679✔
66

67
  hid_t file_id;
68
  if (mpi::master) {
7,679✔
69
    // Create statepoint file
70
    file_id = file_open(filename_, 'w');
6,544✔
71

72
    // Write file type
73
    write_attribute(file_id, "filetype", "statepoint");
6,544✔
74

75
    // Write revision number for state point file
76
    write_attribute(file_id, "version", VERSION_STATEPOINT);
6,544✔
77

78
    // Write OpenMC version
79
    write_attribute(file_id, "openmc_version", VERSION);
6,544✔
80
#ifdef GIT_SHA1
81
    write_attribute(file_id, "git_sha1", GIT_SHA1);
82
#endif
83

84
    // Write current date and time
85
    write_attribute(file_id, "date_and_time", time_stamp());
6,544✔
86

87
    // Write path to input
88
    write_attribute(file_id, "path", settings::path_input);
6,544✔
89

90
    // Write out random number seed
91
    write_dataset(file_id, "seed", openmc_get_seed());
6,544✔
92

93
    // Write out random number stride
94
    write_dataset(file_id, "stride", openmc_get_stride());
6,544✔
95

96
    // Write run information
97
    write_dataset(file_id, "energy_mode",
6,544✔
98
      settings::run_CE ? "continuous-energy" : "multi-group");
99
    switch (settings::run_mode) {
6,544!
100
    case RunMode::FIXED_SOURCE:
2,518✔
101
      write_dataset(file_id, "run_mode", "fixed source");
2,518✔
102
      break;
2,518✔
103
    case RunMode::EIGENVALUE:
4,026✔
104
      write_dataset(file_id, "run_mode", "eigenvalue");
4,026✔
105
      break;
4,026✔
106
    default:
×
107
      break;
×
108
    }
109
    write_attribute(file_id, "photon_transport", settings::photon_transport);
6,544✔
110
    write_dataset(file_id, "n_particles", settings::n_particles);
6,544✔
111
    write_dataset(file_id, "n_batches", settings::n_batches);
6,544✔
112

113
    // Write out current batch number
114
    write_dataset(file_id, "current_batch", simulation::current_batch);
6,544✔
115

116
    // Indicate whether source bank is stored in statepoint
117
    write_attribute(file_id, "source_present", write_source_);
6,544✔
118

119
    // Write out information for eigenvalue run
120
    if (settings::run_mode == RunMode::EIGENVALUE)
6,544✔
121
      write_eigenvalue_hdf5(file_id);
4,026✔
122

123
    hid_t tallies_group = create_group(file_id, "tallies");
6,544✔
124

125
    // Write meshes
126
    meshes_to_hdf5(tallies_group);
6,544✔
127

128
    // Write information for derivatives
129
    if (!model::tally_derivs.empty()) {
6,544✔
130
      hid_t derivs_group = create_group(tallies_group, "derivatives");
11✔
131
      for (const auto& deriv : model::tally_derivs) {
66✔
132
        hid_t deriv_group =
133
          create_group(derivs_group, "derivative " + std::to_string(deriv.id));
55✔
134
        write_dataset(deriv_group, "material", deriv.diff_material);
55✔
135
        if (deriv.variable == DerivativeVariable::DENSITY) {
55✔
136
          write_dataset(deriv_group, "independent variable", "density");
22✔
137
        } else if (deriv.variable == DerivativeVariable::NUCLIDE_DENSITY) {
33✔
138
          write_dataset(deriv_group, "independent variable", "nuclide_density");
22✔
139
          write_dataset(
22✔
140
            deriv_group, "nuclide", data::nuclides[deriv.diff_nuclide]->name_);
22✔
141
        } else if (deriv.variable == DerivativeVariable::TEMPERATURE) {
11!
142
          write_dataset(deriv_group, "independent variable", "temperature");
11✔
143
        } else {
144
          fatal_error("Independent variable for derivative " +
×
145
                      std::to_string(deriv.id) +
×
146
                      " not defined in state_point.cpp");
147
        }
148
        close_group(deriv_group);
55✔
149
      }
150
      close_group(derivs_group);
11✔
151
    }
152

153
    // Write information for filters
154
    hid_t filters_group = create_group(tallies_group, "filters");
6,544✔
155
    write_attribute(filters_group, "n_filters", model::tally_filters.size());
6,544✔
156
    if (!model::tally_filters.empty()) {
6,544✔
157
      // Write filter IDs
158
      vector<int32_t> filter_ids;
4,059✔
159
      filter_ids.reserve(model::tally_filters.size());
4,059✔
160
      for (const auto& filt : model::tally_filters)
14,282✔
161
        filter_ids.push_back(filt->id());
10,223✔
162
      write_attribute(filters_group, "ids", filter_ids);
4,059✔
163

164
      // Write info for each filter
165
      for (const auto& filt : model::tally_filters) {
14,282✔
166
        hid_t filter_group =
167
          create_group(filters_group, "filter " + std::to_string(filt->id()));
10,223✔
168
        filt->to_statepoint(filter_group);
10,223✔
169
        close_group(filter_group);
10,223✔
170
      }
171
    }
4,059✔
172
    close_group(filters_group);
6,544✔
173

174
    // Write information for tallies
175
    write_attribute(tallies_group, "n_tallies", model::tallies.size());
6,544✔
176
    if (!model::tallies.empty()) {
6,544✔
177
      // Write tally IDs
178
      vector<int32_t> tally_ids;
4,587✔
179
      tally_ids.reserve(model::tallies.size());
4,587✔
180
      for (const auto& tally : model::tallies)
25,641✔
181
        tally_ids.push_back(tally->id_);
21,054✔
182
      write_attribute(tallies_group, "ids", tally_ids);
4,587✔
183

184
      // Write all tally information except results
185
      for (const auto& tally : model::tallies) {
25,641✔
186
        hid_t tally_group =
187
          create_group(tallies_group, "tally " + std::to_string(tally->id_));
21,054✔
188

189
        write_dataset(tally_group, "name", tally->name_);
21,054✔
190

191
        if (tally->writable_) {
21,054✔
192
          write_attribute(tally_group, "internal", 0);
20,016✔
193
        } else {
194
          write_attribute(tally_group, "internal", 1);
1,038✔
195
          close_group(tally_group);
1,038✔
196
          continue;
1,038✔
197
        }
198

199
        if (tally->multiply_density()) {
20,016✔
200
          write_attribute(tally_group, "multiply_density", 1);
19,961✔
201
        } else {
202
          write_attribute(tally_group, "multiply_density", 0);
55✔
203
        }
204

205
        if (tally->higher_moments()) {
20,016✔
206
          write_attribute(tally_group, "higher_moments", 1);
11✔
207
        } else {
208
          write_attribute(tally_group, "higher_moments", 0);
20,005✔
209
        }
210

211
        if (tally->estimator_ == TallyEstimator::ANALOG) {
20,016✔
212
          write_dataset(tally_group, "estimator", "analog");
7,198✔
213
        } else if (tally->estimator_ == TallyEstimator::TRACKLENGTH) {
12,818✔
214
          write_dataset(tally_group, "estimator", "tracklength");
11,954✔
215
        } else if (tally->estimator_ == TallyEstimator::COLLISION) {
864!
216
          write_dataset(tally_group, "estimator", "collision");
864✔
217
        }
218

219
        write_dataset(tally_group, "n_realizations", tally->n_realizations_);
20,016✔
220

221
        // Write the ID of each filter attached to this tally
222
        write_dataset(tally_group, "n_filters", tally->filters().size());
20,016✔
223
        if (!tally->filters().empty()) {
20,016✔
224
          vector<int32_t> filter_ids;
18,839✔
225
          filter_ids.reserve(tally->filters().size());
18,839✔
226
          for (auto i_filt : tally->filters())
56,885✔
227
            filter_ids.push_back(model::tally_filters[i_filt]->id());
38,046✔
228
          write_dataset(tally_group, "filters", filter_ids);
18,839✔
229
        }
18,839✔
230

231
        // Write the nuclides this tally scores
232
        vector<std::string> nuclides;
20,016✔
233
        for (auto i_nuclide : tally->nuclides_) {
47,798✔
234
          if (i_nuclide == -1) {
27,782✔
235
            nuclides.push_back("total");
17,387✔
236
          } else {
237
            if (settings::run_CE) {
10,395✔
238
              nuclides.push_back(data::nuclides[i_nuclide]->name_);
10,285✔
239
            } else {
240
              nuclides.push_back(data::mg.nuclides_[i_nuclide].name);
110✔
241
            }
242
          }
243
        }
244
        write_dataset(tally_group, "nuclides", nuclides);
20,016✔
245

246
        if (tally->deriv_ != C_NONE)
20,016✔
247
          write_dataset(
220✔
248
            tally_group, "derivative", model::tally_derivs[tally->deriv_].id);
220✔
249

250
        // Write the tally score bins
251
        vector<std::string> scores;
20,016✔
252
        for (auto sc : tally->scores_)
48,475✔
253
          scores.push_back(reaction_name(sc));
28,459✔
254
        write_dataset(tally_group, "n_score_bins", scores.size());
20,016✔
255
        write_dataset(tally_group, "score_bins", scores);
20,016✔
256

257
        close_group(tally_group);
20,016✔
258
      }
20,016✔
259
    }
4,587✔
260

261
    if (settings::reduce_tallies) {
6,544✔
262
      // Write global tallies
263
      write_dataset(file_id, "global_tallies", simulation::global_tallies);
6,533✔
264

265
      // Write tallies
266
      if (model::active_tallies.size() > 0) {
6,533✔
267
        // Indicate that tallies are on
268
        write_attribute(file_id, "tallies_present", 1);
4,356✔
269

270
        // Write all tally results
271
        for (const auto& tally : model::tallies) {
25,179✔
272
          if (!tally->writable_)
20,823✔
273
            continue;
818✔
274

275
          // Write results for each bin
276
          std::string name = "tally " + std::to_string(tally->id_);
20,005✔
277
          hid_t tally_group = open_group(tallies_group, name.c_str());
20,005✔
278
          auto& results = tally->results_;
20,005✔
279
          write_tally_results(tally_group, results.shape()[0],
20,005✔
280
            results.shape()[1], results.shape()[2], results.data());
20,005✔
281
          close_group(tally_group);
20,005✔
282
        }
20,005✔
283
      } else {
284
        // Indicate tallies are off
285
        write_attribute(file_id, "tallies_present", 0);
2,177✔
286
      }
287
    }
288

289
    close_group(tallies_group);
6,544✔
290
  }
291

292
  // Check for the no-tally-reduction method
293
  if (!settings::reduce_tallies) {
7,679✔
294
    // If using the no-tally-reduction method, we need to collect tally
295
    // results before writing them to the state point file.
296
    write_tally_results_nr(file_id);
16✔
297

298
  } else if (mpi::master) {
7,663✔
299
    // Write number of global realizations
300
    write_dataset(file_id, "n_realizations", simulation::n_realizations);
6,533✔
301
  }
302

303
  if (mpi::master) {
7,679✔
304
    // Write out the runtime metrics.
305
    using namespace simulation;
306
    hid_t runtime_group = create_group(file_id, "runtime");
6,544✔
307
    write_dataset(
6,544✔
308
      runtime_group, "total initialization", time_initialize.elapsed());
309
    write_dataset(
6,544✔
310
      runtime_group, "reading cross sections", time_read_xs.elapsed());
311
    write_dataset(runtime_group, "simulation",
6,544✔
312
      time_inactive.elapsed() + time_active.elapsed());
6,544✔
313
    write_dataset(runtime_group, "transport", time_transport.elapsed());
6,544✔
314
    if (settings::run_mode == RunMode::EIGENVALUE) {
6,544✔
315
      write_dataset(runtime_group, "inactive batches", time_inactive.elapsed());
4,026✔
316
    }
317
    write_dataset(runtime_group, "active batches", time_active.elapsed());
6,544✔
318
    if (settings::run_mode == RunMode::EIGENVALUE) {
6,544✔
319
      write_dataset(
4,026✔
320
        runtime_group, "synchronizing fission bank", time_bank.elapsed());
321
      write_dataset(
4,026✔
322
        runtime_group, "sampling source sites", time_bank_sample.elapsed());
323
      write_dataset(
4,026✔
324
        runtime_group, "SEND-RECV source sites", time_bank_sendrecv.elapsed());
325
    }
326
    write_dataset(
6,544✔
327
      runtime_group, "accumulating tallies", time_tallies.elapsed());
328
    write_dataset(runtime_group, "total", time_total.elapsed());
6,544✔
329
    write_dataset(
6,544✔
330
      runtime_group, "writing statepoints", time_statepoint.elapsed());
331
    close_group(runtime_group);
6,544✔
332

333
    file_close(file_id);
6,544✔
334
  }
335

336
#ifdef PHDF5
337
  bool parallel = true;
4,114✔
338
#else
339
  bool parallel = false;
3,565✔
340
#endif
341

342
  // Write the source bank if desired
343
  if (write_source_) {
7,679✔
344
    if (mpi::master || parallel)
3,876!
345
      file_id = file_open(filename_, 'a', true);
3,876✔
346
    write_source_bank(file_id, simulation::source_bank, simulation::work_index);
3,876✔
347
    if (mpi::master || parallel)
3,876!
348
      file_close(file_id);
3,876✔
349
  }
350

351
#if defined(OPENMC_LIBMESH_ENABLED) || defined(OPENMC_DAGMC_ENABLED)
352
  // write unstructured mesh tally files
353
  write_unstructured_mesh_results();
2,289✔
354
#endif
355

356
  simulation::time_statepoint.stop();
7,679✔
357

358
  return 0;
7,679✔
359
}
7,679✔
360

361
void restart_set_keff()
65✔
362
{
363
  if (simulation::restart_batch > settings::n_inactive) {
65!
364
    for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) {
309✔
365
      simulation::k_sum[0] += simulation::k_generation[i];
244✔
366
      simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
244✔
367
    }
368
    int n = settings::gen_per_batch * simulation::n_realizations;
65✔
369
    simulation::keff = simulation::k_sum[0] / n;
65✔
370
  } else {
371
    simulation::keff = simulation::k_generation.back();
×
372
  }
373
}
65✔
374

375
void load_state_point()
65✔
376
{
377
  write_message(
65✔
378
    fmt::format("Loading state point {}...", settings::path_statepoint_c), 5);
118✔
379
  openmc_statepoint_load(settings::path_statepoint.c_str());
65✔
380
}
65✔
381

382
void statepoint_version_check(hid_t file_id)
65✔
383
{
384
  // Read revision number for state point file and make sure it matches with
385
  // current version
386
  array<int, 2> version_array;
387
  read_attribute(file_id, "version", version_array);
65✔
388
  if (version_array != VERSION_STATEPOINT) {
65!
389
    fatal_error(
×
390
      "State point version does not match current version in OpenMC.");
391
  }
392
}
65✔
393

394
extern "C" int openmc_statepoint_load(const char* filename)
65✔
395
{
396
  // Open file for reading
397
  hid_t file_id = file_open(filename, 'r', true);
65✔
398

399
  // Read filetype
400
  std::string word;
65✔
401
  read_attribute(file_id, "filetype", word);
65✔
402
  if (word != "statepoint") {
65!
403
    fatal_error("OpenMC tried to restart from a non-statepoint file.");
×
404
  }
405

406
  statepoint_version_check(file_id);
65✔
407

408
  // Read and overwrite random number seed
409
  int64_t seed;
410
  read_dataset(file_id, "seed", seed);
65✔
411
  openmc_set_seed(seed);
65✔
412

413
  // Read and overwrite random number stride
414
  uint64_t stride;
415
  read_dataset(file_id, "stride", stride);
65✔
416
  openmc_set_stride(stride);
65✔
417

418
  // It is not impossible for a state point to be generated from a CE run but
419
  // to be loaded in to an MG run (or vice versa), check to prevent that.
420
  read_dataset(file_id, "energy_mode", word);
65✔
421
  if (word == "multi-group" && settings::run_CE) {
65!
422
    fatal_error("State point file is from multigroup run but current run is "
×
423
                "continous energy.");
424
  } else if (word == "continuous-energy" && !settings::run_CE) {
65!
425
    fatal_error("State point file is from continuous-energy run but current "
×
426
                "run is multigroup!");
427
  }
428

429
  // Read and overwrite run information except number of batches
430
  read_dataset(file_id, "run_mode", word);
65✔
431
  if (word == "fixed source") {
65!
432
    settings::run_mode = RunMode::FIXED_SOURCE;
×
433
  } else if (word == "eigenvalue") {
65!
434
    settings::run_mode = RunMode::EIGENVALUE;
65✔
435
  }
436
  read_attribute(file_id, "photon_transport", settings::photon_transport);
65✔
437
  read_dataset(file_id, "n_particles", settings::n_particles);
65✔
438
  int temp;
439
  read_dataset(file_id, "n_batches", temp);
65✔
440

441
  // Take maximum of statepoint n_batches and input n_batches
442
  settings::n_batches = std::max(settings::n_batches, temp);
65✔
443

444
  // Read batch number to restart at
445
  read_dataset(file_id, "current_batch", simulation::restart_batch);
65✔
446

447
  if (settings::restart_run &&
65!
448
      simulation::restart_batch >= settings::n_max_batches) {
65✔
449
    warning(fmt::format(
11✔
450
      "The number of batches specified for simulation ({}) is smaller "
451
      "than or equal to the number of batches in the restart statepoint file "
452
      "({})",
453
      settings::n_max_batches, simulation::restart_batch));
454
  }
455

456
  // Logical flag for source present in statepoint file
457
  bool source_present;
458
  read_attribute(file_id, "source_present", source_present);
65✔
459

460
  // Read information specific to eigenvalue run
461
  if (settings::run_mode == RunMode::EIGENVALUE) {
65!
462
    read_dataset(file_id, "n_inactive", temp);
65✔
463
    read_eigenvalue_hdf5(file_id);
65✔
464

465
    // Take maximum of statepoint n_inactive and input n_inactive
466
    settings::n_inactive = std::max(settings::n_inactive, temp);
65✔
467

468
    // Check to make sure source bank is present
469
    if (settings::path_sourcepoint == settings::path_statepoint &&
130!
470
        !source_present) {
65!
471
      fatal_error("Source bank must be contained in statepoint restart file");
×
472
    }
473
  }
474

475
  // Read number of realizations for global tallies
476
  read_dataset(file_id, "n_realizations", simulation::n_realizations);
65✔
477

478
  // Set k_sum, keff, and current_batch based on whether restart file is part
479
  // of active cycle or inactive cycle
480
  if (settings::run_mode == RunMode::EIGENVALUE) {
65!
481
    restart_set_keff();
65✔
482
  }
483

484
  // Set current batch number
485
  simulation::current_batch = simulation::restart_batch;
65✔
486

487
  // Read tallies to master. If we are using Parallel HDF5, all processes
488
  // need to be included in the HDF5 calls.
489
#ifdef PHDF5
490
  if (true) {
491
#else
492
  if (mpi::master) {
30!
493
#endif
494
    // Read global tally data
495
    read_dataset_lowlevel(file_id, "global_tallies", H5T_NATIVE_DOUBLE, H5S_ALL,
65✔
496
      false, simulation::global_tallies.data());
65✔
497

498
    // Check if tally results are present
499
    bool present;
500
    read_attribute(file_id, "tallies_present", present);
65✔
501

502
    // Read in sum and sum squared
503
    if (present) {
65!
504
      hid_t tallies_group = open_group(file_id, "tallies");
65✔
505

506
      for (auto& tally : model::tallies) {
223✔
507
        // Read sum, sum_sq, and N for each bin
508
        std::string name = "tally " + std::to_string(tally->id_);
158✔
509
        hid_t tally_group = open_group(tallies_group, name.c_str());
158✔
510

511
        int internal = 0;
158✔
512
        if (attribute_exists(tally_group, "internal")) {
158!
513
          read_attribute(tally_group, "internal", internal);
158✔
514
        }
515
        if (internal) {
158!
516
          tally->writable_ = false;
×
517
        } else {
518
          auto& results = tally->results_;
158✔
519
          read_tally_results(tally_group, results.shape()[0],
316✔
520
            results.shape()[1], results.shape()[2], results.data());
158✔
521

522
          read_dataset(tally_group, "n_realizations", tally->n_realizations_);
158✔
523
          close_group(tally_group);
158✔
524
        }
525
      }
158✔
526
      close_group(tallies_group);
65✔
527
    }
528
  }
529

530
  // Read source if in eigenvalue mode
531
  if (settings::run_mode == RunMode::EIGENVALUE) {
65!
532

533
    // Check if source was written out separately
534
    if (!source_present) {
65!
535

536
      // Close statepoint file
537
      file_close(file_id);
×
538

539
      // Write message
540
      write_message(
×
541
        "Loading source file " + settings::path_sourcepoint + "...", 5);
×
542

543
      // Open source file
544
      file_id = file_open(settings::path_sourcepoint.c_str(), 'r', true);
×
545
    }
546

547
    // Read source
548
    read_source_bank(file_id, simulation::source_bank, true);
65✔
549
  }
550

551
  // Close file
552
  file_close(file_id);
65✔
553

554
  return 0;
65✔
555
}
65✔
556

557
hid_t h5banktype(bool memory)
10,288✔
558
{
559
  // Create compound type for position
560
  hid_t postype = H5Tcreate(H5T_COMPOUND, sizeof(struct Position));
10,288✔
561
  H5Tinsert(postype, "x", HOFFSET(Position, x), H5T_NATIVE_DOUBLE);
10,288✔
562
  H5Tinsert(postype, "y", HOFFSET(Position, y), H5T_NATIVE_DOUBLE);
10,288✔
563
  H5Tinsert(postype, "z", HOFFSET(Position, z), H5T_NATIVE_DOUBLE);
10,288✔
564

565
  // Create bank datatype
566
  //
567
  // If you make changes to the compound datatype here, make sure you update:
568
  // - openmc/source.py
569
  // - openmc/statepoint.py
570
  // - docs/source/io_formats/statepoint.rst
571
  // - docs/source/io_formats/source.rst
572
  auto n = sizeof(SourceSite);
10,288✔
573
  if (!memory)
10,288✔
574
    n = 2 * sizeof(struct Position) + 3 * sizeof(double) + 3 * sizeof(int);
5,080✔
575
  hid_t banktype = H5Tcreate(H5T_COMPOUND, n);
10,288✔
576
  H5Tinsert(banktype, "r", HOFFSET(SourceSite, r), postype);
10,288✔
577
  H5Tinsert(banktype, "u", HOFFSET(SourceSite, u), postype);
10,288✔
578
  H5Tinsert(banktype, "E", HOFFSET(SourceSite, E), H5T_NATIVE_DOUBLE);
10,288✔
579
  H5Tinsert(banktype, "time", HOFFSET(SourceSite, time), H5T_NATIVE_DOUBLE);
10,288✔
580
  H5Tinsert(banktype, "wgt", HOFFSET(SourceSite, wgt), H5T_NATIVE_DOUBLE);
10,288✔
581
  H5Tinsert(banktype, "delayed_group", HOFFSET(SourceSite, delayed_group),
10,288✔
582
    H5T_NATIVE_INT);
10,288✔
583
  H5Tinsert(banktype, "surf_id", HOFFSET(SourceSite, surf_id), H5T_NATIVE_INT);
10,288✔
584
  H5Tinsert(
10,288✔
585
    banktype, "particle", HOFFSET(SourceSite, particle), H5T_NATIVE_INT);
10,288✔
586

587
  H5Tclose(postype);
10,288✔
588
  return banktype;
10,288✔
589
}
590

591
void write_source_point(std::string filename, span<SourceSite> source_bank,
1,242✔
592
  const vector<int64_t>& bank_index, bool use_mcpl)
593
{
594
  std::string ext = use_mcpl ? "mcpl" : "h5";
1,242✔
595
  write_message("Creating source file {}.{} with {} particles ...", filename,
1,242✔
596
    ext, source_bank.size(), 5);
1,242✔
597

598
  // Dispatch to appropriate function based on file type
599
  if (use_mcpl) {
1,242✔
600
    filename.append(".mcpl");
38✔
601
    write_mcpl_source_point(filename.c_str(), source_bank, bank_index);
38✔
602
  } else {
603
    filename.append(".h5");
1,204✔
604
    write_h5_source_point(filename.c_str(), source_bank, bank_index);
1,204✔
605
  }
606
}
1,242✔
607

608
void write_h5_source_point(const char* filename, span<SourceSite> source_bank,
1,204✔
609
  const vector<int64_t>& bank_index)
610
{
611
  // When using parallel HDF5, the file is written to collectively by all
612
  // processes. With MPI-only, the file is opened and written by the master
613
  // (note that the call to write_source_bank is by all processes since slave
614
  // processes need to send source bank data to the master.
615
#ifdef PHDF5
616
  bool parallel = true;
587✔
617
#else
618
  bool parallel = false;
617✔
619
#endif
620

621
  if (!filename)
1,204!
UNCOV
622
    fatal_error("write_source_point filename needs a nonempty name.");
×
623

624
  std::string filename_(filename);
1,204✔
625
  const auto extension = get_file_extension(filename_);
1,204✔
626
  if (extension != "h5") {
1,204!
UNCOV
627
    warning("write_source_point was passed a file extension differing "
×
628
            "from .h5, but an hdf5 file will be written.");
629
  }
630

631
  hid_t file_id;
632
  if (mpi::master || parallel) {
1,204!
633
    file_id = file_open(filename_.c_str(), 'w', true);
1,204✔
634
    write_attribute(file_id, "filetype", "source");
1,204✔
635
  }
636

637
  // Get pointer to source bank and write to file
638
  write_source_bank(file_id, source_bank, bank_index);
1,204✔
639

640
  if (mpi::master || parallel)
1,204!
641
    file_close(file_id);
1,204✔
642
}
1,204✔
643

644
void write_source_bank(hid_t group_id, span<SourceSite> source_bank,
5,080✔
645
  const vector<int64_t>& bank_index)
646
{
647
  hid_t membanktype = h5banktype(true);
5,080✔
648
  hid_t filebanktype = h5banktype(false);
5,080✔
649

650
#ifdef OPENMC_MPI
651
  write_bank_dataset("source_bank", group_id, source_bank, bank_index,
2,736✔
652
    membanktype, filebanktype, mpi::source_site);
653
#else
654
  write_bank_dataset("source_bank", group_id, source_bank, bank_index,
2,344✔
655
    membanktype, filebanktype);
656
#endif
657

658
  H5Tclose(membanktype);
5,080✔
659
  H5Tclose(filebanktype);
5,080✔
660
}
5,080✔
661

662
// Determine member names of a compound HDF5 datatype
663
std::string dtype_member_names(hid_t dtype_id)
256✔
664
{
665
  int nmembers = H5Tget_nmembers(dtype_id);
256✔
666
  std::string names;
256✔
667
  for (int i = 0; i < nmembers; i++) {
2,259✔
668
    char* name = H5Tget_member_name(dtype_id, i);
2,003✔
669
    names = names.append(name);
2,003✔
670
    H5free_memory(name);
2,003✔
671
    if (i < nmembers - 1)
2,003✔
672
      names += ", ";
1,747✔
673
  }
674
  return names;
256✔
UNCOV
675
}
×
676

677
void read_source_bank(
128✔
678
  hid_t group_id, vector<SourceSite>& sites, bool distribute)
679
{
680
  hid_t banktype = h5banktype(true);
128✔
681

682
  // Open the dataset
683
  hid_t dset = H5Dopen(group_id, "source_bank", H5P_DEFAULT);
128✔
684

685
  // Make sure number of members matches
686
  hid_t dtype = H5Dget_type(dset);
128✔
687
  auto file_member_names = dtype_member_names(dtype);
128✔
688
  auto bank_member_names = dtype_member_names(banktype);
128✔
689
  if (file_member_names != bank_member_names) {
128✔
690
    fatal_error(fmt::format(
9✔
691
      "Source site attributes in file do not match what is "
692
      "expected for this version of OpenMC. File attributes = ({}). Expected "
693
      "attributes = ({})",
694
      file_member_names, bank_member_names));
695
  }
696

697
  hid_t dspace = H5Dget_space(dset);
119✔
698
  hsize_t n_sites;
699
  H5Sget_simple_extent_dims(dspace, &n_sites, nullptr);
119✔
700

701
  // Make sure vector is big enough in case where we're reading entire source on
702
  // each process
703
  if (!distribute)
119✔
704
    sites.resize(n_sites);
54✔
705

706
  hid_t memspace;
707
  if (distribute) {
119✔
708
    if (simulation::work_index[mpi::n_procs] > n_sites) {
65!
UNCOV
709
      fatal_error("Number of source sites in source file is less "
×
710
                  "than number of source particles per generation.");
711
    }
712

713
    // Create another data space but for each proc individually
714
    hsize_t n_sites_local = simulation::work_per_rank;
65✔
715
    memspace = H5Screate_simple(1, &n_sites_local, nullptr);
65✔
716

717
    // Select hyperslab for each process
718
    hsize_t offset = simulation::work_index[mpi::rank];
65✔
719
    H5Sselect_hyperslab(
65✔
720
      dspace, H5S_SELECT_SET, &offset, nullptr, &n_sites_local, nullptr);
721
  } else {
722
    memspace = H5S_ALL;
54✔
723
  }
724

725
#ifdef PHDF5
726
  // Read data in parallel
727
  hid_t plist = H5Pcreate(H5P_DATASET_XFER);
65✔
728
  H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
65✔
729
  H5Dread(dset, banktype, memspace, dspace, plist, sites.data());
65✔
730
  H5Pclose(plist);
65✔
731
#else
732
  H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, sites.data());
54✔
733
#endif
734

735
  // Close all ids
736
  H5Sclose(dspace);
119✔
737
  if (distribute)
119✔
738
    H5Sclose(memspace);
65✔
739
  H5Dclose(dset);
119✔
740
  H5Tclose(banktype);
119✔
741
}
119✔
742

743
void write_unstructured_mesh_results()
2,289✔
744
{
745

746
  for (auto& tally : model::tallies) {
10,839✔
747

748
    vector<std::string> tally_scores;
8,550✔
749
    for (auto filter_idx : tally->filters()) {
25,034✔
750
      auto& filter = model::tally_filters[filter_idx];
16,484✔
751
      if (filter->type() != FilterType::MESH)
16,484!
752
        continue;
16,469✔
753

754
      // check if the filter uses an unstructured mesh
755
      auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
1,922!
756
      auto mesh_idx = mesh_filter->mesh();
1,922!
757
      auto umesh =
758
        dynamic_cast<UnstructuredMesh*>(model::meshes[mesh_idx].get());
1,922!
759

760
      if (!umesh)
1,922✔
761
        continue;
1,887✔
762

763
      if (!umesh->output_)
35!
UNCOV
764
        continue;
×
765

766
      if (umesh->library() == "moab") {
35!
767
        if (mpi::master)
20✔
768
          warning(fmt::format(
10!
769
            "Output for a MOAB mesh (mesh {}) was "
770
            "requested but will not be written. Please use the Python "
771
            "API to generated the desired VTK tetrahedral mesh.",
772
            umesh->id_));
10✔
773
        continue;
20✔
774
      }
775

776
      // if this tally has more than one filter, print
777
      // warning and skip writing the mesh
778
      if (tally->filters().size() > 1) {
15!
UNCOV
779
        warning(fmt::format("Skipping unstructured mesh writing for tally "
×
780
                            "{}. More than one filter is present on the tally.",
UNCOV
781
          tally->id_));
×
UNCOV
782
        break;
×
783
      }
784

785
      int n_realizations = tally->n_realizations_;
15✔
786

787
      for (int score_idx = 0; score_idx < tally->scores_.size(); score_idx++) {
30✔
788
        for (int nuc_idx = 0; nuc_idx < tally->nuclides_.size(); nuc_idx++) {
30✔
789
          // combine the score and nuclide into a name for the value
790
          auto score_str = fmt::format("{}_{}", tally->score_name(score_idx),
30!
791
            tally->nuclide_name(nuc_idx));
30!
792
          // add this score to the mesh
793
          // (this is in a separate loop because all variables need to be added
794
          //  to libMesh's equation system before any are initialized, which
795
          //  happens in set_score_data)
796
          umesh->add_score(score_str);
15!
797
        }
15✔
798
      }
799

800
      for (int score_idx = 0; score_idx < tally->scores_.size(); score_idx++) {
30✔
801
        for (int nuc_idx = 0; nuc_idx < tally->nuclides_.size(); nuc_idx++) {
30✔
802
          // combine the score and nuclide into a name for the value
803
          auto score_str = fmt::format("{}_{}", tally->score_name(score_idx),
30!
804
            tally->nuclide_name(nuc_idx));
30!
805

806
          // index for this nuclide and score
807
          int nuc_score_idx = score_idx + nuc_idx * tally->scores_.size();
15✔
808

809
          // construct result vectors
810
          vector<double> mean_vec(umesh->n_bins()),
15!
811
            std_dev_vec(umesh->n_bins());
15!
812
          for (int j = 0; j < tally->results_.shape()[0]; j++) {
146,799✔
813
            // get the volume for this bin
814
            double volume = umesh->volume(j);
146,784!
815
            // compute the mean
816
            double mean = tally->results_(j, nuc_score_idx, TallyResult::SUM) /
146,784!
817
                          n_realizations;
146,784✔
818
            mean_vec.at(j) = mean / volume;
146,784!
819

820
            // compute the standard deviation
821
            double sum_sq =
822
              tally->results_(j, nuc_score_idx, TallyResult::SUM_SQ);
146,784!
823
            double std_dev {0.0};
146,784✔
824
            if (n_realizations > 1) {
146,784!
825
              std_dev = sum_sq / n_realizations - mean * mean;
146,784✔
826
              std_dev = std::sqrt(std_dev / (n_realizations - 1));
146,784✔
827
            }
828
            std_dev_vec[j] = std_dev / volume;
146,784✔
829
          }
830
#ifdef OPENMC_MPI
831
          MPI_Bcast(
10!
832
            mean_vec.data(), mean_vec.size(), MPI_DOUBLE, 0, mpi::intracomm);
10✔
833
          MPI_Bcast(std_dev_vec.data(), std_dev_vec.size(), MPI_DOUBLE, 0,
10!
834
            mpi::intracomm);
835
#endif
836
          // set the data for this score
837
          umesh->set_score_data(score_str, mean_vec, std_dev_vec);
15!
838
        }
15✔
839
      }
840

841
      // Generate a file name based on the tally id
842
      // and the current batch number
843
      size_t batch_width {std::to_string(settings::n_max_batches).size()};
15!
844
      std::string filename = fmt::format("tally_{0}.{1:0{2}}", tally->id_,
15✔
UNCOV
845
        simulation::current_batch, batch_width);
×
846

847
      // Write the unstructured mesh and data to file
848
      umesh->write(filename);
15!
849

850
      // remove score data added for this mesh write
851
      umesh->remove_scores();
15!
852
    }
15✔
853
  }
8,550✔
854
}
2,289✔
855

856
void write_tally_results_nr(hid_t file_id)
16✔
857
{
858
  // ==========================================================================
859
  // COLLECT AND WRITE GLOBAL TALLIES
860

861
  hid_t tallies_group;
862
  if (mpi::master) {
16✔
863
    // Write number of realizations
864
    write_dataset(file_id, "n_realizations", simulation::n_realizations);
11✔
865

866
    tallies_group = open_group(file_id, "tallies");
11✔
867
  }
868

869
  // Get global tallies
870
  auto& gt = simulation::global_tallies;
16✔
871

872
#ifdef OPENMC_MPI
873
  // Reduce global tallies
874
  xt::xtensor<double, 2> gt_reduced = xt::empty_like(gt);
10✔
875
  MPI_Reduce(gt.data(), gt_reduced.data(), gt.size(), MPI_DOUBLE, MPI_SUM, 0,
10✔
876
    mpi::intracomm);
877

878
  // Transfer values to value on master
879
  if (mpi::master) {
10✔
880
    if (simulation::current_batch == settings::n_max_batches ||
5!
881
        simulation::satisfy_triggers) {
882
      std::copy(gt_reduced.begin(), gt_reduced.end(), gt.begin());
5✔
883
    }
884
  }
885
#endif
886

887
  // Write out global tallies sum and sum_sq
888
  if (mpi::master) {
16✔
889
    write_dataset(file_id, "global_tallies", gt);
11✔
890
  }
891

892
  for (const auto& t : model::tallies) {
32✔
893
    // Skip any tallies that are not active
894
    if (!t->active_)
16!
UNCOV
895
      continue;
×
896
    if (!t->writable_)
16!
UNCOV
897
      continue;
×
898

899
    if (mpi::master && !attribute_exists(file_id, "tallies_present")) {
16!
900
      write_attribute(file_id, "tallies_present", 1);
11✔
901
    }
902

903
    // Get view of accumulated tally values
904
    auto values_view = xt::view(t->results_, xt::all(), xt::all(),
16✔
905
      xt::range(static_cast<int>(TallyResult::SUM),
16✔
906
        static_cast<int>(TallyResult::SUM_SQ) + 1));
16✔
907

908
    // Make copy of tally values in contiguous array
909
    xt::xtensor<double, 3> values = values_view;
16✔
910

911
    if (mpi::master) {
16✔
912
      // Open group for tally
913
      std::string groupname {"tally " + std::to_string(t->id_)};
11✔
914
      hid_t tally_group = open_group(tallies_group, groupname.c_str());
11✔
915

916
      // The MPI_IN_PLACE specifier allows the master to copy values into
917
      // a receive buffer without having a temporary variable
918
#ifdef OPENMC_MPI
919
      MPI_Reduce(MPI_IN_PLACE, values.data(), values.size(), MPI_DOUBLE,
5✔
920
        MPI_SUM, 0, mpi::intracomm);
921
#endif
922

923
      // At the end of the simulation, store the results back in the
924
      // regular TallyResults array
925
      if (simulation::current_batch == settings::n_max_batches ||
11!
926
          simulation::satisfy_triggers) {
927
        values_view = values;
11✔
928
      }
929

930
      // Put in temporary tally result
931
      xt::xtensor<double, 3> results_copy = xt::zeros_like(t->results_);
11✔
932
      auto copy_view = xt::view(results_copy, xt::all(), xt::all(),
11✔
933
        xt::range(static_cast<int>(TallyResult::SUM),
11✔
934
          static_cast<int>(TallyResult::SUM_SQ) + 1));
11✔
935
      copy_view = values;
11✔
936

937
      // Write reduced tally results to file
938
      auto shape = results_copy.shape();
11✔
939
      write_tally_results(
11✔
940
        tally_group, shape[0], shape[1], shape[2], results_copy.data());
11✔
941

942
      close_group(tally_group);
11✔
943
    } else {
11✔
944
      // Receive buffer not significant at other processors
945
#ifdef OPENMC_MPI
946
      MPI_Reduce(values.data(), nullptr, values.size(), MPI_DOUBLE, MPI_SUM, 0,
5✔
947
        mpi::intracomm);
948
#endif
949
    }
950
  }
16✔
951

952
  if (mpi::master) {
16✔
953
    if (!object_exists(file_id, "tallies_present")) {
11!
954
      // Indicate that tallies are off
955
      write_dataset(file_id, "tallies_present", 0);
11✔
956
    }
957

958
    close_group(tallies_group);
11✔
959
  }
960
}
16✔
961

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