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

openmc-dev / openmc / 25218672962

01 May 2026 02:43PM UTC coverage: 81.452% (+0.1%) from 81.326%
25218672962

Pull #3863

github

web-flow
Merge bf4da9a8a into 368ea069c
Pull Request #3863: Shared Secondary Particle Bank

17907 of 25811 branches covered (69.38%)

Branch coverage included in aggregate %.

415 of 427 new or added lines in 17 files covered. (97.19%)

1149 existing lines in 27 files now uncovered.

58842 of 68415 relevant lines covered (86.01%)

48738312.94 hits per line

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

81.73
/src/output.cpp
1
#include "openmc/output.h"
2

3
#include <algorithm> // for transform, max
4
#include <cstdio>    // for stdout
5
#include <cstring>   // for strlen
6
#include <ctime>     // for time, localtime
7
#include <fstream>
8
#include <iomanip> // for setw, setprecision, put_time
9
#include <ios>     // for fixed, scientific, left
10
#include <iostream>
11
#include <sstream>
12
#include <unordered_map>
13
#include <utility> // for pair
14

15
#include <fmt/core.h>
16
#include <fmt/ostream.h>
17
#ifdef _OPENMP
18
#include <omp.h>
19
#endif
20
#include "openmc/tensor.h"
21

22
#include "openmc/capi.h"
23
#include "openmc/cell.h"
24
#include "openmc/constants.h"
25
#include "openmc/eigenvalue.h"
26
#include "openmc/error.h"
27
#include "openmc/geometry.h"
28
#include "openmc/lattice.h"
29
#include "openmc/math_functions.h"
30
#include "openmc/message_passing.h"
31
#include "openmc/mgxs_interface.h"
32
#include "openmc/nuclide.h"
33
#include "openmc/plot.h"
34
#include "openmc/random_ray/flat_source_domain.h"
35
#include "openmc/reaction.h"
36
#include "openmc/settings.h"
37
#include "openmc/simulation.h"
38
#include "openmc/surface.h"
39
#include "openmc/tallies/derivative.h"
40
#include "openmc/tallies/filter.h"
41
#include "openmc/tallies/tally.h"
42
#include "openmc/tallies/tally_scoring.h"
43
#include "openmc/timer.h"
44

45
namespace openmc {
46

47
#ifdef OPENMC_ENABLE_STRICT_FP
48
const bool STRICT_FP_ENABLED = true;
49
#else
50
const bool STRICT_FP_ENABLED = false;
51
#endif
52

53
//==============================================================================
54

55
void title()
7,340✔
56
{
57
  fmt::print("                                %%%%%%%%%%%%%%%\n"
7,340✔
58
             "                           %%%%%%%%%%%%%%%%%%%%%%%%\n"
59
             "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
60
             "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
61
             "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
62
             "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
63
             "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n"
64
             "                                     %%%%%%%%%%%%%%%%%%%%%%%%\n"
65
             "                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%\n"
66
             "                ##################     %%%%%%%%%%%%%%%%%%%%%%%\n"
67
             "                ###################     %%%%%%%%%%%%%%%%%%%%%%%\n"
68
             "                ####################     %%%%%%%%%%%%%%%%%%%%%%\n"
69
             "                #####################     %%%%%%%%%%%%%%%%%%%%%\n"
70
             "                ######################     %%%%%%%%%%%%%%%%%%%%\n"
71
             "                #######################     %%%%%%%%%%%%%%%%%%\n"
72
             "                 #######################     %%%%%%%%%%%%%%%%%\n"
73
             "                 ######################     %%%%%%%%%%%%%%%%%\n"
74
             "                  ####################     %%%%%%%%%%%%%%%%%\n"
75
             "                    #################     %%%%%%%%%%%%%%%%%\n"
76
             "                     ###############     %%%%%%%%%%%%%%%%\n"
77
             "                       ############     %%%%%%%%%%%%%%%\n"
78
             "                          ########     %%%%%%%%%%%%%%\n"
79
             "                                      %%%%%%%%%%%\n\n");
80

81
  // Write version information
82
  fmt::print(
14,680✔
83
    "                 | The OpenMC Monte Carlo Code\n"
84
    "       Copyright | 2011-2026 MIT, UChicago Argonne LLC, and contributors\n"
85
    "         License | https://docs.openmc.org/en/latest/license.html\n"
86
    "         Version | {}.{}.{}{}{}\n",
87
    VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE, VERSION_DEV ? "-dev" : "",
7,340✔
88
    VERSION_COMMIT_COUNT);
89
  fmt::print("     Commit Hash | {}\n", VERSION_COMMIT_HASH);
7,340✔
90

91
  // Write the date and time
92
  fmt::print("       Date/Time | {}\n", time_stamp());
8,659✔
93

94
#ifdef OPENMC_MPI
95
  // Write number of processors
96
  fmt::print("   MPI Processes | {}\n", mpi::n_procs);
2,706✔
97
#endif
98

99
#ifdef _OPENMP
100
  // Write number of OpenMP threads
101
  fmt::print("  OpenMP Threads | {}\n", omp_get_max_threads());
4,012✔
102
#endif
103
  fmt::print("\n");
7,340✔
104
  std::fflush(stdout);
7,340✔
105
}
7,340✔
106

107
//==============================================================================
108

109
std::string header(const char* msg)
39,815✔
110
{
111
  // Determine how many times to repeat the '=' character.
112
  int n_prefix = (63 - strlen(msg)) / 2;
39,815✔
113
  int n_suffix = n_prefix;
39,815✔
114
  if ((strlen(msg) % 2) == 0)
39,815✔
115
    ++n_suffix;
7,877✔
116

117
  // Convert to uppercase.
118
  std::string upper(msg);
39,815✔
119
  std::transform(upper.begin(), upper.end(), upper.begin(), ::toupper);
39,815✔
120

121
  // Add ===>  <=== markers.
122
  std::stringstream out;
39,815✔
123
  out << ' ';
39,815✔
124
  for (int i = 0; i < n_prefix; i++)
1,009,749✔
125
    out << '=';
969,934✔
126
  out << ">     " << upper << "     <";
39,815✔
127
  for (int i = 0; i < n_suffix; i++)
1,017,626✔
128
    out << '=';
977,811✔
129

130
  return out.str();
39,815✔
131
}
39,815✔
132

133
std::string header(const std::string& msg)
19,697✔
134
{
135
  return header(msg.c_str());
19,697✔
136
}
137

138
void header(const char* msg, int level)
20,118✔
139
{
140
  auto out = header(msg);
20,118✔
141

142
  // Print header based on verbosity level.
143
  if (settings::verbosity >= level) {
20,118✔
144
    fmt::print("\n{}\n\n", out);
19,682✔
145
    std::fflush(stdout);
19,682✔
146
  }
147
}
20,118✔
148

149
//==============================================================================
150

151
std::string time_stamp()
22,263✔
152
{
153
  std::stringstream ts;
22,263✔
154
  std::time_t t = std::time(nullptr); // get time now
22,263✔
155
  ts << std::put_time(std::localtime(&t), "%Y-%m-%d %H:%M:%S");
22,263✔
156
  return ts.str();
44,526✔
157
}
22,263✔
158

159
//==============================================================================
160

161
void print_particle(Particle& p)
44✔
162
{
163
  // Display particle type and ID.
164
  switch (p.type().pdg_number()) {
44!
165
  case PDG_NEUTRON:
44✔
166
    fmt::print("Neutron ");
44✔
167
    break;
44✔
168
  case PDG_PHOTON:
×
169
    fmt::print("Photon ");
×
170
    break;
×
171
  case PDG_ELECTRON:
×
172
    fmt::print("Electron ");
×
173
    break;
×
174
  case PDG_POSITRON:
×
175
    fmt::print("Positron ");
×
176
    break;
×
177
  default:
×
178
    fmt::print("Particle {} ", p.type().str());
×
179
  }
180
  fmt::print("{}\n", p.id());
44✔
181

182
  // Display particle geometry hierarchy.
183
  for (auto i = 0; i < p.n_coord(); i++) {
88✔
184
    fmt::print("  Level {}\n", i);
44✔
185

186
    if (p.coord(i).cell() != C_NONE) {
44!
187
      const Cell& c {*model::cells[p.coord(i).cell()]};
44✔
188
      fmt::print("    Cell             = {}\n", c.id_);
44✔
189
    }
190

191
    if (p.coord(i).universe() != C_NONE) {
44!
192
      const Universe& u {*model::universes[p.coord(i).universe()]};
44✔
193
      fmt::print("    Universe         = {}\n", u.id_);
44✔
194
    }
195

196
    if (p.coord(i).lattice() != C_NONE) {
44!
197
      const Lattice& lat {*model::lattices[p.coord(i).lattice()]};
×
198
      fmt::print("    Lattice          = {}\n", lat.id_);
×
199
      fmt::print("    Lattice position = ({},{},{})\n",
×
200
        p.coord(i).lattice_index()[0], p.coord(i).lattice_index()[1],
×
201
        p.coord(i).lattice_index()[2]);
×
202
    }
203

204
    fmt::print("    r = {}\n", p.coord(i).r());
44✔
205
    fmt::print("    u = {}\n", p.coord(i).u());
44✔
206
  }
207

208
  // Display miscellaneous info.
209
  if (p.surface() != SURFACE_NONE) {
44!
210
    // Surfaces identifiers are >= 1, but indices are >= 0 so we need -1
211
    const Surface& surf {*model::surfaces[p.surface_index()]};
×
212
    fmt::print("  Surface = {}\n", (p.surface() > 0) ? surf.id_ : -surf.id_);
×
213
  }
214
  fmt::print("  Weight = {}\n", p.wgt());
44✔
215
  if (settings::run_CE) {
44!
216
    fmt::print("  Energy = {}\n", p.E());
44✔
217
  } else {
218
    fmt::print("  Energy Group = {}\n", p.g());
×
219
  }
220
  fmt::print("  Delayed Group = {}\n\n", p.delayed_group());
44✔
221
}
44✔
222

223
//==============================================================================
224

225
void print_plot()
88✔
226
{
227
  header("PLOTTING SUMMARY", 5);
88✔
228
  if (settings::verbosity < 5)
88!
229
    return;
230

231
  for (const auto& pl : model::plots) {
330✔
232
    fmt::print("Plot ID: {}\n", pl->id());
242✔
233
    fmt::print("Plot file: {}\n", pl->path_plot());
242✔
234
    fmt::print("Universe depth: {}\n", pl->level());
242✔
235
    pl->print_info(); // prints type-specific plot info
242✔
236
    fmt::print("\n");
242✔
237
  }
238
}
239

240
//==============================================================================
241

242
void print_overlap_check()
×
243
{
244
#ifdef OPENMC_MPI
245
  vector<int64_t> temp(model::overlap_check_count);
246
  MPI_Reduce(temp.data(), model::overlap_check_count.data(),
×
247
    model::overlap_check_count.size(), MPI_INT64_T, MPI_SUM, 0, mpi::intracomm);
×
248
#endif
249

250
  if (mpi::master) {
×
251
    header("cell overlap check summary", 1);
×
252
    fmt::print(" Cell ID      No. Overlap Checks\n");
×
253

254
    vector<int32_t> sparse_cell_ids;
×
255
    for (int i = 0; i < model::cells.size(); i++) {
×
256
      fmt::print(
×
257
        " {:8} {:17}\n", model::cells[i]->id_, model::overlap_check_count[i]);
×
258
      if (model::overlap_check_count[i] < 10) {
×
259
        sparse_cell_ids.push_back(model::cells[i]->id_);
×
260
      }
261
    }
262

263
    fmt::print("\n There were {} cells with less than 10 overlap checks\n",
×
264
      sparse_cell_ids.size());
×
265
    for (auto id : sparse_cell_ids) {
×
266
      fmt::print(" {}", id);
×
267
    }
268
    fmt::print("\n");
×
269
  }
×
270
}
×
271

272
//==============================================================================
273

274
void print_usage()
×
275
{
276
  if (mpi::master) {
×
277
    fmt::print(
×
278
      "Usage: openmc [options] [path]\n\n"
279
      "Options:\n"
280
      "  -c, --volume           Run in stochastic volume calculation mode\n"
281
      "  -g, --geometry-debug   Run with geometry debugging on\n"
282
      "  -n, --particles        Number of particles per generation\n"
283
      "  -p, --plot             Run in plotting mode\n"
284
      "  -r, --restart          Restart a previous run from a state point\n"
285
      "                         or a particle restart file\n"
286
      "  -s, --threads          Number of OpenMP threads\n"
287
      "  -t, --track            Write tracks for all particles (up to "
288
      "max_tracks)\n"
289
      "  -e, --event            Run using event-based parallelism\n"
290
      "  -q, --verbosity        Output verbosity\n"
291
      "  -v, --version          Show version information\n"
292
      "  -h, --help             Show this message\n");
293
  }
294
}
×
295

296
//==============================================================================
297

298
void print_version()
11✔
299
{
300
  if (mpi::master) {
11!
301
    fmt::print("OpenMC version {}.{}.{}{}{}\n", VERSION_MAJOR, VERSION_MINOR,
22✔
302
      VERSION_RELEASE, VERSION_DEV ? "-dev" : "", VERSION_COMMIT_COUNT);
11✔
303
    fmt::print("Commit hash: {}\n", VERSION_COMMIT_HASH);
11✔
304
    fmt::print("Copyright (c) 2011-2026 MIT, UChicago Argonne LLC, and "
11✔
305
               "contributors\nMIT/X license at "
306
               "<https://docs.openmc.org/en/latest/license.html>\n");
307
  }
308
}
11✔
309

310
//==============================================================================
311

312
void print_build_info()
11✔
313
{
314
  const std::string n("no");
11✔
315
  const std::string y("yes");
11✔
316

317
  std::string mpi(n);
11✔
318
  std::string phdf5(n);
11✔
319
  std::string dagmc(n);
11✔
320
  std::string libmesh(n);
11✔
321
  std::string png(n);
11✔
322
  std::string profiling(n);
11✔
323
  std::string coverage(n);
11✔
324
  std::string uwuw(n);
11✔
325
  std::string strict_fp(n);
11✔
326

327
#ifdef PHDF5
328
  phdf5 = y;
4✔
329
#endif
330
#ifdef OPENMC_MPI
331
  mpi = y;
4✔
332
#endif
333
#ifdef OPENMC_DAGMC_ENABLED
334
  dagmc = y;
1✔
335
#endif
336
#ifdef OPENMC_LIBMESH_ENABLED
337
  libmesh = y;
2✔
338
#endif
339
#ifdef USE_LIBPNG
340
  png = y;
11✔
341
#endif
342
#ifdef PROFILINGBUILD
343
  profiling = y;
344
#endif
345
#ifdef COVERAGEBUILD
346
  coverage = y;
11✔
347
#endif
348
#ifdef OPENMC_UWUW_ENABLED
349
  uwuw = y;
1✔
350
#endif
351
#ifdef OPENMC_ENABLE_STRICT_FP
352
  strict_fp = y;
11✔
353
#endif
354

355
  // Wraps macro variables in quotes
356
#define STRINGIFY(x) STRINGIFY2(x)
357
#define STRINGIFY2(x) #x
358

359
  if (mpi::master) {
11!
360
    fmt::print("Build type:            {}\n", STRINGIFY(BUILD_TYPE));
11✔
361
    fmt::print("Compiler ID:           {} {}\n", STRINGIFY(COMPILER_ID),
11✔
362
      STRINGIFY(COMPILER_VERSION));
363
    fmt::print("MPI enabled:           {}\n", mpi);
11✔
364
    fmt::print("Parallel HDF5 enabled: {}\n", phdf5);
11✔
365
    fmt::print("PNG support:           {}\n", png);
11✔
366
    fmt::print("DAGMC support:         {}\n", dagmc);
11✔
367
    fmt::print("libMesh support:       {}\n", libmesh);
11✔
368
    fmt::print("Coverage testing:      {}\n", coverage);
11✔
369
    fmt::print("Profiling flags:       {}\n", profiling);
11✔
370
    fmt::print("UWUW support:          {}\n", uwuw);
11✔
371
    fmt::print("Strict FP:             {}\n", strict_fp);
13✔
372
  }
373
}
11✔
374

375
//==============================================================================
376

377
void print_columns()
3,395✔
378
{
379
  if (settings::entropy_on) {
3,395✔
380
    fmt::print("  Bat./Gen.      k       Entropy         Average k \n"
495✔
381
               "  =========   ========   ========   ====================\n");
382
  } else {
383
    fmt::print("  Bat./Gen.      k            Average k\n"
2,900✔
384
               "  =========   ========   ====================\n");
385
  }
386
}
3,395✔
387

388
//==============================================================================
389

390
void print_generation()
75,918✔
391
{
392
  // Determine overall generation index and number of active generations
393
  int idx = overall_generation() - 1;
75,918✔
394
  int n = simulation::current_batch > settings::n_inactive
151,836✔
395
            ? settings::gen_per_batch * simulation::n_realizations +
75,918✔
396
                simulation::current_gen
397
            : 0;
398

399
  // write out batch/generation and generation k-effective
400
  auto batch_and_gen = std::to_string(simulation::current_batch) + "/" +
151,836✔
401
                       std::to_string(simulation::current_gen);
151,836✔
402
  fmt::print("  {:>9}   {:8.5f}", batch_and_gen, simulation::k_generation[idx]);
75,918✔
403

404
  // write out entropy info
405
  if (settings::entropy_on) {
75,918✔
406
    fmt::print("   {:8.5f}", simulation::entropy[idx]);
13,255✔
407
  }
408

409
  if (n > 1) {
75,918✔
410
    fmt::print("   {:8.5f} +/-{:8.5f}", simulation::keff, simulation::keff_std);
53,806✔
411
  }
412
  fmt::print("\n");
75,918✔
413
  std::fflush(stdout);
75,918✔
414
}
75,918✔
415

416
//==============================================================================
417

418
void show_time(const char* label, double secs, int indent_level)
71,867✔
419
{
420
  int width = 33 - indent_level * 2;
71,867✔
421
  fmt::print("{0:{1}} {2:<{3}} = {4:>10.4e} seconds\n", "", 2 * indent_level,
71,867✔
422
    label, width, secs);
423
}
71,867✔
424

425
void show_rate(const char* label, double particles_per_sec)
7,829✔
426
{
427
  fmt::print(" {:<33} = {:.6} particles/second\n", label, particles_per_sec);
7,829✔
428
}
7,829✔
429

430
void print_runtime()
5,594✔
431
{
432
  using namespace simulation;
5,594✔
433

434
  // display header block
435
  header("Timing Statistics", 6);
5,594✔
436
  if (settings::verbosity < 6)
5,594!
437
    return;
438

439
  // display time elapsed for various sections
440
  show_time("Total time for initialization", time_initialize.elapsed());
5,594✔
441
  show_time("Reading cross sections", time_read_xs.elapsed(), 1);
5,594✔
442
  show_time("Total time in simulation",
5,594✔
443
    time_inactive.elapsed() + time_active.elapsed());
5,594✔
444
  show_time("Time in transport only", time_transport.elapsed(), 1);
5,594✔
445
  if (settings::event_based) {
5,594✔
446
    show_time("Particle initialization", time_event_init.elapsed(), 2);
159✔
447
    show_time("XS lookups", time_event_calculate_xs.elapsed(), 2);
159✔
448
    show_time("Advancing", time_event_advance_particle.elapsed(), 2);
159✔
449
    show_time("Surface crossings", time_event_surface_crossing.elapsed(), 2);
159✔
450
    show_time("Collisions", time_event_collision.elapsed(), 2);
159✔
451
    show_time("Particle death", time_event_death.elapsed(), 2);
159✔
452
  }
453
  if (settings::run_mode == RunMode::EIGENVALUE) {
5,594✔
454
    show_time("Time in inactive batches", time_inactive.elapsed(), 1);
3,120✔
455
  }
456
  show_time("Time in active batches", time_active.elapsed(), 1);
5,594✔
457
  if (settings::run_mode == RunMode::EIGENVALUE) {
5,594✔
458
    show_time("Time synchronizing fission bank", time_bank.elapsed(), 1);
3,120✔
459
    show_time("Sampling source sites", time_bank_sample.elapsed(), 2);
3,120✔
460
    show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2);
3,120✔
461
  }
462
  show_time("Time accumulating tallies", time_tallies.elapsed(), 1);
5,594✔
463
  show_time("Time writing statepoints", time_statepoint.elapsed(), 1);
5,594✔
464
  show_time("Total time for finalization", time_finalize.elapsed());
5,594✔
465
  show_time("Total time elapsed", time_total.elapsed());
5,594✔
466

467
  // Calculate particle rate in active/inactive batches
468
  int n_active = simulation::current_batch - settings::n_inactive;
5,594✔
469
  double speed_inactive = 0.0;
5,594✔
470
  double speed_active;
5,594✔
471
  if (settings::restart_run) {
5,594✔
472
    if (simulation::restart_batch < settings::n_inactive) {
33!
UNCOV
473
      speed_inactive = (settings::n_particles *
×
UNCOV
474
                         (settings::n_inactive - simulation::restart_batch) *
×
UNCOV
475
                         settings::gen_per_batch) /
×
UNCOV
476
                       time_inactive.elapsed();
×
UNCOV
477
      speed_active =
×
478
        (settings::n_particles * n_active * settings::gen_per_batch) /
×
479
        time_active.elapsed();
×
480
    } else {
481
      speed_active = (settings::n_particles *
33✔
482
                       (settings::n_batches - simulation::restart_batch) *
33✔
483
                       settings::gen_per_batch) /
33✔
484
                     time_active.elapsed();
33✔
485
    }
486
  } else {
487
    if (settings::n_inactive > 0) {
5,561✔
488
      speed_inactive = (settings::n_particles * settings::n_inactive *
2,235✔
489
                         settings::gen_per_batch) /
2,235✔
490
                       time_inactive.elapsed();
2,235✔
491
    }
492
    speed_active =
5,561✔
493
      (settings::n_particles * n_active * settings::gen_per_batch) /
5,561✔
494
      time_active.elapsed();
5,561✔
495
  }
496

497
  // display calculation rate
498
  if (!(settings::restart_run &&
5,594✔
499
        (simulation::restart_batch >= settings::n_inactive)) &&
33!
500
      settings::n_inactive > 0) {
5,561✔
501
    show_rate("Calculation Rate (inactive)", speed_inactive);
2,235✔
502
  }
503
  show_rate("Calculation Rate (active)", speed_active);
5,594✔
504

505
  // Display track rate when weight windows are enabled
506
  if (settings::weight_windows_on) {
5,594✔
507
    double speed_tracks =
269✔
508
      simulation::simulation_tracks_completed / time_active.elapsed();
269✔
509
    fmt::print(
269✔
510
      " {:<33} = {:.6} tracks/second\n", "Track Rate (active)", speed_tracks);
511
  }
512
}
513

514
//==============================================================================
515

516
std::pair<double, double> mean_stdev(const double* x, int n)
13,837,297✔
517
{
518
  double mean = x[static_cast<int>(TallyResult::SUM)] / n;
13,837,297✔
519
  double stdev =
13,837,297✔
520
    n > 1 ? std::sqrt(std::max(0.0,
27,666,388✔
521
              (x[static_cast<int>(TallyResult::SUM_SQ)] / n - mean * mean) /
27,666,388✔
522
                (n - 1)))
13,829,091✔
523
          : 0.0;
13,837,297✔
524
  return {mean, stdev};
13,837,297✔
525
}
526

527
//==============================================================================
528

529
void print_results()
5,594✔
530
{
531
  // display header block for results
532
  header("Results", 4);
5,594✔
533
  if (settings::verbosity < 4)
5,594!
UNCOV
534
    return;
×
535

536
  // Calculate t-value for confidence intervals
537
  int n = simulation::n_realizations;
5,594✔
538
  double alpha, t_n1, t_n3;
5,594✔
539
  if (settings::confidence_intervals) {
5,594✔
540
    alpha = 1.0 - CONFIDENCE_LEVEL;
11✔
541
    t_n1 = t_percentile(1.0 - alpha / 2.0, n - 1);
11✔
542
    t_n3 = t_percentile(1.0 - alpha / 2.0, n - 3);
11✔
543
  } else {
544
    t_n1 = 1.0;
545
    t_n3 = 1.0;
546
  }
547

548
  // write global tallies
549
  const auto& gt = simulation::global_tallies;
5,594✔
550
  double mean, stdev;
5,594✔
551
  if (n > 1) {
5,594✔
552
    if (settings::run_mode == RunMode::EIGENVALUE) {
5,298✔
553
      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_COLLISION, 0), n);
3,087✔
554
      fmt::print(" k-effective (Collision)     = {:.5f} +/- {:.5f}\n", mean,
6,174✔
555
        t_n1 * stdev);
3,087✔
556
      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_TRACKLENGTH, 0), n);
3,087✔
557
      fmt::print(" k-effective (Track-length)  = {:.5f} +/- {:.5f}\n", mean,
6,174✔
558
        t_n1 * stdev);
3,087✔
559
      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_ABSORPTION, 0), n);
3,087✔
560
      fmt::print(" k-effective (Absorption)    = {:.5f} +/- {:.5f}\n", mean,
6,174✔
561
        t_n1 * stdev);
3,087✔
562
      if (n > 3) {
3,087✔
563
        double k_combined[2];
3,021✔
564
        openmc_get_keff(k_combined);
3,021✔
565
        fmt::print(" Combined k-effective        = {:.5f} +/- {:.5f}\n",
3,021✔
566
          k_combined[0], k_combined[1]);
567
      }
568
    }
569
    std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::LEAKAGE, 0), n);
5,298✔
570
    fmt::print(
9,622✔
571
      " Leakage Fraction            = {:.5f} +/- {:.5f}\n", mean, t_n1 * stdev);
5,298✔
572
  } else {
573
    if (mpi::master)
296!
574
      warning("Could not compute uncertainties -- only one "
592✔
575
              "active batch simulated!");
576

577
    if (settings::run_mode == RunMode::EIGENVALUE) {
296✔
578
      fmt::print(" k-effective (Collision)    = {:.5f}\n",
66✔
579
        gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n);
33✔
580
      fmt::print(" k-effective (Track-length) = {:.5f}\n",
66✔
581
        gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n);
33✔
582
      fmt::print(" k-effective (Absorption)   = {:.5f}\n",
60✔
583
        gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n);
33✔
584
    }
585
    fmt::print(" Leakage Fraction           = {:.5f}\n",
537✔
586
      gt(GlobalTally::LEAKAGE, TallyResult::SUM) / n);
296✔
587
  }
588
  fmt::print("\n");
5,594✔
589
  std::fflush(stdout);
5,594✔
590
}
591

592
//==============================================================================
593

594
const std::unordered_map<int, const char*> score_names = {
595
  {SCORE_FLUX, "Flux"},
596
  {SCORE_TOTAL, "Total Reaction Rate"},
597
  {SCORE_SCATTER, "Scattering Rate"},
598
  {SCORE_NU_SCATTER, "Scattering Production Rate"},
599
  {SCORE_ABSORPTION, "Absorption Rate"},
600
  {SCORE_FISSION, "Fission Rate"},
601
  {SCORE_NU_FISSION, "Nu-Fission Rate"},
602
  {SCORE_KAPPA_FISSION, "Kappa-Fission Rate"},
603
  {SCORE_EVENTS, "Events"},
604
  {SCORE_DECAY_RATE, "Decay Rate"},
605
  {SCORE_DELAYED_NU_FISSION, "Delayed-Nu-Fission Rate"},
606
  {SCORE_PROMPT_NU_FISSION, "Prompt-Nu-Fission Rate"},
607
  {SCORE_INVERSE_VELOCITY, "Flux-Weighted Inverse Velocity"},
608
  {SCORE_FISS_Q_PROMPT, "Prompt fission power"},
609
  {SCORE_FISS_Q_RECOV, "Recoverable fission power"},
610
  {SCORE_CURRENT, "Current"},
611
  {SCORE_PULSE_HEIGHT, "pulse-height"},
612
  {SCORE_IFP_TIME_NUM, "IFP lifetime numerator"},
613
  {SCORE_IFP_BETA_NUM, "IFP delayed fraction numerator"},
614
  {SCORE_IFP_DENOM, "IFP common denominator"},
615
};
616

617
//! Create an ASCII output file showing all tally results.
618

619
void write_tallies()
6,280✔
620
{
621
  if (model::tallies.empty())
6,280✔
622
    return;
2,056✔
623

624
  // Set filename for tallies_out
625
  std::string filename = fmt::format("{}tallies.out", settings::path_output);
4,224✔
626

627
  // Open the tallies.out file.
628
  std::ofstream tallies_out;
4,224✔
629
  tallies_out.open(filename, std::ios::out | std::ios::trunc);
4,224✔
630

631
  // Loop over each tally.
632
  for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
23,921✔
633
    const auto& tally {*model::tallies[i_tally]};
19,697✔
634

635
    // Write header block.
636
    std::string tally_header("TALLY " + std::to_string(tally.id_));
19,697✔
637
    if (!tally.name_.empty())
19,697✔
638
      tally_header += ": " + tally.name_;
4,632✔
639
    fmt::print(tallies_out, "{}\n\n", header(tally_header));
19,697✔
640

641
    if (!tally.writable_) {
19,697✔
642
      fmt::print(tallies_out, " Internal\n\n");
579✔
643
      continue;
579✔
644
    }
645

646
    // Calculate t-value for confidence intervals
647
    double t_value = 1;
19,118✔
648
    if (settings::confidence_intervals) {
19,118✔
649
      auto alpha = 1 - CONFIDENCE_LEVEL;
11✔
650
      t_value = t_percentile(1 - alpha * 0.5, tally.n_realizations_ - 1);
11✔
651
    }
652

653
    // Write derivative information.
654
    if (tally.deriv_ != C_NONE) {
19,118✔
655
      const auto& deriv {model::tally_derivs[tally.deriv_]};
220!
656
      switch (deriv.variable) {
220!
657
      case DerivativeVariable::DENSITY:
88✔
658
        fmt::print(tallies_out, " Density derivative Material {}\n",
176✔
659
          deriv.diff_material);
88✔
660
        break;
88✔
661
      case DerivativeVariable::NUCLIDE_DENSITY:
88✔
662
        fmt::print(tallies_out,
176✔
663
          " Nuclide density derivative Material {} Nuclide {}\n",
664
          deriv.diff_material, data::nuclides[deriv.diff_nuclide]->name_);
88✔
665
        break;
88✔
666
      case DerivativeVariable::TEMPERATURE:
44✔
667
        fmt::print(tallies_out, " Temperature derivative Material {}\n",
88✔
668
          deriv.diff_material);
44✔
669
        break;
44✔
UNCOV
670
      default:
×
UNCOV
671
        fatal_error(fmt::format("Differential tally dependent variable for "
×
672
                                "tally {} not defined in output.cpp",
UNCOV
673
          tally.id_));
×
674
      }
675
    }
676

677
    // Initialize Filter Matches Object
678
    vector<FilterMatch> filter_matches;
38,236✔
679
    // Allocate space for tally filter matches
680
    filter_matches.resize(model::tally_filters.size());
19,118✔
681

682
    // Loop over all filter bin combinations.
683
    auto filter_iter = FilterBinIter(tally, false, &filter_matches);
19,118✔
684
    auto end = FilterBinIter(tally, true, &filter_matches);
19,118✔
685
    for (; filter_iter != end; ++filter_iter) {
9,537,930✔
686
      auto filter_index = filter_iter.index_;
687

688
      // Print info about this combination of filter bins.  The stride check
689
      // prevents redundant output.
690
      int indent = 0;
691
      for (auto i = 0; i < tally.filters().size(); ++i) {
30,061,934✔
692
        if (filter_index % tally.strides(i) == 0) {
20,543,122✔
693
          auto i_filt = tally.filters(i);
14,345,656✔
694
          const auto& filt {*model::tally_filters[i_filt]};
14,345,656✔
695
          auto& match {filter_matches[i_filt]};
14,345,656✔
696
          fmt::print(tallies_out, "{0:{1}}{2}\n", "", indent + 1,
14,345,656✔
697
            filt.text_label(match.i_bin_));
28,691,312✔
698
        }
699
        indent += 2;
20,543,122✔
700
      }
701

702
      // Loop over all nuclide and score combinations.
703
      int score_index = 0;
9,518,812✔
704
      for (auto i_nuclide : tally.nuclides_) {
20,465,589✔
705
        // Write label for this nuclide bin.
706
        if (i_nuclide == -1) {
10,946,777✔
707
          fmt::print(tallies_out, "{0:{1}}Total Material\n", "", indent + 1);
8,110,174✔
708
        } else {
709
          if (settings::run_CE) {
2,836,603✔
710
            fmt::print(tallies_out, "{0:{1}}{2}\n", "", indent + 1,
2,836,229✔
711
              data::nuclides[i_nuclide]->name_);
2,836,229✔
712
          } else {
713
            fmt::print(tallies_out, "{0:{1}}{2}\n", "", indent + 1,
374✔
714
              data::mg.nuclides_[i_nuclide].name);
374✔
715
          }
716
        }
717

718
        // Write the score, mean, and uncertainty.
719
        indent += 2;
10,946,777✔
720
        for (auto score : tally.scores_) {
24,769,515✔
721
          std::string score_name =
13,822,738✔
722
            score > 0 ? reaction_name(score) : score_names.at(score);
13,822,738✔
723
          double mean, stdev;
13,822,738✔
724
          std::tie(mean, stdev) =
13,822,738✔
725
            mean_stdev(&tally.results_(filter_index, score_index, 0),
27,645,476✔
726
              tally.n_realizations_);
13,822,738✔
727
          fmt::print(tallies_out, "{0:{1}}{2:<36} {3:.6} +/- {4:.6}\n", "",
27,645,476✔
728
            indent + 1, score_name, mean, t_value * stdev);
13,822,738✔
729
          score_index += 1;
13,822,738✔
730
        }
13,822,738✔
731
        indent -= 2;
10,946,777✔
732
      }
733
    }
734
  }
19,697✔
735
}
4,224✔
736

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