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

openmc-dev / openmc / 13688481331

06 Mar 2025 12:35AM UTC coverage: 85.002% (-0.04%) from 85.042%
13688481331

Pull #3333

github

web-flow
Merge 5654c92ef into ced892912
Pull Request #3333: Random Ray Source Region Mesh Subdivision (Cell-Under-Voxel Geometry)

560 of 646 new or added lines in 11 files covered. (86.69%)

361 existing lines in 16 files now uncovered.

51473 of 60555 relevant lines covered (85.0%)

36077311.57 hits per line

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

86.12
/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 "xtensor/xview.hpp"
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
//==============================================================================
48

49
void title()
5,227✔
50
{
51
  fmt::print("                                %%%%%%%%%%%%%%%\n"
4,285✔
52
             "                           %%%%%%%%%%%%%%%%%%%%%%%%\n"
53
             "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
54
             "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
55
             "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
56
             "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
57
             "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n"
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\n");
74

75
  // Write version information
76
  fmt::print(
4,285✔
77
    "                 | The OpenMC Monte Carlo Code\n"
78
    "       Copyright | 2011-2025 MIT, UChicago Argonne LLC, and contributors\n"
79
    "         License | https://docs.openmc.org/en/latest/license.html\n"
80
    "         Version | {}.{}.{}{}{}\n",
81
    VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE, VERSION_DEV ? "-dev" : "",
10,454✔
82
    VERSION_COMMIT_COUNT);
83
  fmt::print("     Commit Hash | {}\n", VERSION_COMMIT_HASH);
4,285✔
84

85
  // Write the date and time
86
  fmt::print("       Date/Time | {}\n", time_stamp());
10,454✔
87

88
#ifdef OPENMC_MPI
89
  // Write number of processors
90
  fmt::print("   MPI Processes | {}\n", mpi::n_procs);
1,936✔
91
#endif
92

93
#ifdef _OPENMP
94
  // Write number of OpenMP threads
95
  fmt::print("  OpenMP Threads | {}\n", omp_get_max_threads());
5,718✔
96
#endif
97
  fmt::print("\n");
5,227✔
98
  std::fflush(stdout);
5,227✔
99
}
5,227✔
100

101
//==============================================================================
102

103
std::string header(const char* msg)
33,224✔
104
{
105
  // Determine how many times to repeat the '=' character.
106
  int n_prefix = (63 - strlen(msg)) / 2;
33,224✔
107
  int n_suffix = n_prefix;
33,224✔
108
  if ((strlen(msg) % 2) == 0)
33,224✔
109
    ++n_suffix;
7,554✔
110

111
  // Convert to uppercase.
112
  std::string upper(msg);
33,224✔
113
  std::transform(upper.begin(), upper.end(), upper.begin(), ::toupper);
33,224✔
114

115
  // Add ===>  <=== markers.
116
  std::stringstream out;
33,224✔
117
  out << ' ';
33,224✔
118
  for (int i = 0; i < n_prefix; i++)
855,448✔
119
    out << '=';
822,224✔
120
  out << ">     " << upper << "     <";
33,224✔
121
  for (int i = 0; i < n_suffix; i++)
863,002✔
122
    out << '=';
829,778✔
123

124
  return out.str();
66,448✔
125
}
33,224✔
126

127
std::string header(const std::string& msg)
18,358✔
128
{
129
  return header(msg.c_str());
18,358✔
130
}
131

132
void header(const char* msg, int level)
14,866✔
133
{
134
  auto out = header(msg);
14,866✔
135

136
  // Print header based on verbosity level.
137
  if (settings::verbosity >= level) {
14,866✔
138
    fmt::print("\n{}\n\n", out);
11,871✔
139
    std::fflush(stdout);
14,547✔
140
  }
141
}
14,866✔
142

143
//==============================================================================
144

145
std::string time_stamp()
15,932✔
146
{
147
  std::stringstream ts;
15,932✔
148
  std::time_t t = std::time(nullptr); // get time now
15,932✔
149
  ts << std::put_time(std::localtime(&t), "%Y-%m-%d %H:%M:%S");
15,932✔
150
  return ts.str();
31,864✔
151
}
15,932✔
152

153
//==============================================================================
154

155
void print_particle(Particle& p)
33✔
156
{
157
  // Display particle type and ID.
158
  switch (p.type()) {
33✔
159
  case ParticleType::neutron:
33✔
160
    fmt::print("Neutron ");
27✔
161
    break;
33✔
162
  case ParticleType::photon:
×
163
    fmt::print("Photon ");
×
164
    break;
×
165
  case ParticleType::electron:
×
166
    fmt::print("Electron ");
×
167
    break;
×
168
  case ParticleType::positron:
×
169
    fmt::print("Positron ");
×
170
    break;
×
171
  default:
×
172
    fmt::print("Unknown Particle ");
×
173
  }
174
  fmt::print("{}\n", p.id());
60✔
175

176
  // Display particle geometry hierarchy.
177
  for (auto i = 0; i < p.n_coord(); i++) {
66✔
178
    fmt::print("  Level {}\n", i);
27✔
179

180
    if (p.coord(i).cell != C_NONE) {
33✔
181
      const Cell& c {*model::cells[p.coord(i).cell]};
33✔
182
      fmt::print("    Cell             = {}\n", c.id_);
66✔
183
    }
184

185
    if (p.coord(i).universe != C_NONE) {
33✔
186
      const Universe& u {*model::universes[p.coord(i).universe]};
33✔
187
      fmt::print("    Universe         = {}\n", u.id_);
66✔
188
    }
189

190
    if (p.coord(i).lattice != C_NONE) {
33✔
191
      const Lattice& lat {*model::lattices[p.coord(i).lattice]};
×
192
      fmt::print("    Lattice          = {}\n", lat.id_);
×
193
      fmt::print("    Lattice position = ({},{},{})\n", p.coord(i).lattice_i[0],
×
194
        p.coord(i).lattice_i[1], p.coord(i).lattice_i[2]);
×
195
    }
196

197
    fmt::print("    r = {}\n", p.coord(i).r);
60✔
198
    fmt::print("    u = {}\n", p.coord(i).u);
66✔
199
  }
200

201
  // Display miscellaneous info.
202
  if (p.surface() != SURFACE_NONE) {
33✔
203
    // Surfaces identifiers are >= 1, but indices are >= 0 so we need -1
204
    const Surface& surf {*model::surfaces[p.surface_index()]};
×
205
    fmt::print("  Surface = {}\n", (p.surface() > 0) ? surf.id_ : -surf.id_);
×
206
  }
207
  fmt::print("  Weight = {}\n", p.wgt());
60✔
208
  if (settings::run_CE) {
33✔
209
    fmt::print("  Energy = {}\n", p.E());
66✔
210
  } else {
211
    fmt::print("  Energy Group = {}\n", p.g());
×
212
  }
213
  fmt::print("  Delayed Group = {}\n\n", p.delayed_group());
60✔
214
}
33✔
215

216
//==============================================================================
217

218
void print_plot()
264✔
219
{
220
  header("PLOTTING SUMMARY", 5);
264✔
221
  if (settings::verbosity < 5)
264✔
222
    return;
×
223

224
  for (const auto& pl : model::plots) {
792✔
225
    fmt::print("Plot ID: {}\n", pl->id());
960✔
226
    fmt::print("Plot file: {}\n", pl->path_plot());
960✔
227
    fmt::print("Universe depth: {}\n", pl->level());
960✔
228
    pl->print_info(); // prints type-specific plot info
528✔
229
    fmt::print("\n");
528✔
230
  }
231
}
232

233
//==============================================================================
234

235
void print_overlap_check()
×
236
{
237
#ifdef OPENMC_MPI
238
  vector<int64_t> temp(model::overlap_check_count);
239
  MPI_Reduce(temp.data(), model::overlap_check_count.data(),
240
    model::overlap_check_count.size(), MPI_INT64_T, MPI_SUM, 0, mpi::intracomm);
241
#endif
242

243
  if (mpi::master) {
×
244
    header("cell overlap check summary", 1);
×
245
    fmt::print(" Cell ID      No. Overlap Checks\n");
×
246

247
    vector<int32_t> sparse_cell_ids;
×
248
    for (int i = 0; i < model::cells.size(); i++) {
×
249
      fmt::print(
×
250
        " {:8} {:17}\n", model::cells[i]->id_, model::overlap_check_count[i]);
×
251
      if (model::overlap_check_count[i] < 10) {
×
252
        sparse_cell_ids.push_back(model::cells[i]->id_);
×
253
      }
254
    }
255

256
    fmt::print("\n There were {} cells with less than 10 overlap checks\n",
×
257
      sparse_cell_ids.size());
×
258
    for (auto id : sparse_cell_ids) {
×
259
      fmt::print(" {}", id);
×
260
    }
261
    fmt::print("\n");
×
262
  }
263
}
264

265
//==============================================================================
266

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

288
//==============================================================================
289

290
void print_version()
11✔
291
{
292
  if (mpi::master) {
11✔
293
    fmt::print("OpenMC version {}.{}.{}{}{}\n", VERSION_MAJOR, VERSION_MINOR,
9✔
294
      VERSION_RELEASE, VERSION_DEV ? "-dev" : "", VERSION_COMMIT_COUNT);
22✔
295
    fmt::print("Commit hash: {}\n", VERSION_COMMIT_HASH);
11✔
296
    fmt::print("Copyright (c) 2011-2025 MIT, UChicago Argonne LLC, and "
11✔
297
               "contributors\nMIT/X license at "
298
               "<https://docs.openmc.org/en/latest/license.html>\n");
299
  }
300
}
11✔
301

302
//==============================================================================
303

304
void print_build_info()
11✔
305
{
306
  const std::string n("no");
11✔
307
  const std::string y("yes");
11✔
308

309
  std::string mpi(n);
11✔
310
  std::string phdf5(n);
11✔
311
  std::string dagmc(n);
11✔
312
  std::string libmesh(n);
11✔
313
  std::string png(n);
11✔
314
  std::string profiling(n);
11✔
315
  std::string coverage(n);
11✔
316
  std::string mcpl(n);
11✔
317
  std::string uwuw(n);
11✔
318

319
#ifdef PHDF5
320
  phdf5 = y;
5✔
321
#endif
322
#ifdef OPENMC_MPI
323
  mpi = y;
5✔
324
#endif
325
#ifdef DAGMC
326
  dagmc = y;
1✔
327
#endif
328
#ifdef LIBMESH
329
  libmesh = y;
2✔
330
#endif
331
#ifdef OPENMC_MCPL
332
  mcpl = y;
11✔
333
#endif
334
#ifdef USE_LIBPNG
335
  png = y;
11✔
336
#endif
337
#ifdef PROFILINGBUILD
338
  profiling = y;
339
#endif
340
#ifdef COVERAGEBUILD
341
  coverage = y;
11✔
342
#endif
343
#ifdef OPENMC_UWUW
344
  uwuw = y;
1✔
345
#endif
346

347
  // Wraps macro variables in quotes
348
#define STRINGIFY(x) STRINGIFY2(x)
349
#define STRINGIFY2(x) #x
350

351
  if (mpi::master) {
11✔
352
    fmt::print("Build type:            {}\n", STRINGIFY(BUILD_TYPE));
11✔
353
    fmt::print("Compiler ID:           {} {}\n", STRINGIFY(COMPILER_ID),
11✔
354
      STRINGIFY(COMPILER_VERSION));
355
    fmt::print("MPI enabled:           {}\n", mpi);
11✔
356
    fmt::print("Parallel HDF5 enabled: {}\n", phdf5);
11✔
357
    fmt::print("PNG support:           {}\n", png);
11✔
358
    fmt::print("DAGMC support:         {}\n", dagmc);
11✔
359
    fmt::print("libMesh support:       {}\n", libmesh);
11✔
360
    fmt::print("MCPL support:          {}\n", mcpl);
11✔
361
    fmt::print("Coverage testing:      {}\n", coverage);
11✔
362
    fmt::print("Profiling flags:       {}\n", profiling);
11✔
363
    fmt::print("UWUW support:          {}\n", uwuw);
11✔
364
  }
365
}
11✔
366

367
//==============================================================================
368

369
void print_columns()
2,626✔
370
{
371
  if (settings::entropy_on) {
2,626✔
372
    fmt::print("  Bat./Gen.      k       Entropy         Average k \n"
308✔
373
               "  =========   ========   ========   ====================\n");
374
  } else {
375
    fmt::print("  Bat./Gen.      k            Average k\n"
2,318✔
376
               "  =========   ========   ====================\n");
377
  }
378
}
2,626✔
379

380
//==============================================================================
381

382
void print_generation()
51,646✔
383
{
384
  // Determine overall generation index and number of active generations
385
  int idx = overall_generation() - 1;
51,646✔
386
  int n = simulation::current_batch > settings::n_inactive
103,292✔
387
            ? settings::gen_per_batch * simulation::n_realizations +
51,646✔
388
                simulation::current_gen
389
            : 0;
390

391
  // write out batch/generation and generation k-effective
392
  auto batch_and_gen = std::to_string(simulation::current_batch) + "/" +
103,292✔
393
                       std::to_string(simulation::current_gen);
103,292✔
394
  fmt::print("  {:>9}   {:8.5f}", batch_and_gen, simulation::k_generation[idx]);
93,900✔
395

396
  // write out entropy info
397
  if (settings::entropy_on) {
51,646✔
398
    fmt::print("   {:8.5f}", simulation::entropy[idx]);
18,370✔
399
  }
400

401
  if (n > 1) {
51,646✔
402
    fmt::print("   {:8.5f} +/-{:8.5f}", simulation::keff, simulation::keff_std);
39,047✔
403
  }
404
  fmt::print("\n");
42,254✔
405
  std::fflush(stdout);
51,646✔
406
}
51,646✔
407

408
//==============================================================================
409

410
void show_time(const char* label, double secs, int indent_level)
54,154✔
411
{
412
  int width = 33 - indent_level * 2;
54,154✔
413
  fmt::print("{0:{1}} {2:<{3}} = {4:>10.4e} seconds\n", "", 2 * indent_level,
98,482✔
414
    label, width, secs);
415
}
54,154✔
416

417
void show_rate(const char* label, double particles_per_sec)
6,037✔
418
{
419
  fmt::print(" {:<33} = {:.6} particles/second\n", label, particles_per_sec);
6,037✔
420
}
6,037✔
421

422
void print_runtime()
4,208✔
423
{
424
  using namespace simulation;
425

426
  // display header block
427
  header("Timing Statistics", 6);
4,208✔
428
  if (settings::verbosity < 6)
4,208✔
UNCOV
429
    return;
×
430

431
  // display time elapsed for various sections
432
  show_time("Total time for initialization", time_initialize.elapsed());
4,208✔
433
  show_time("Reading cross sections", time_read_xs.elapsed(), 1);
4,208✔
434
  show_time("Total time in simulation",
4,208✔
435
    time_inactive.elapsed() + time_active.elapsed());
4,208✔
436
  show_time("Time in transport only", time_transport.elapsed(), 1);
4,208✔
437
  if (settings::event_based) {
4,208✔
438
    show_time("Particle initialization", time_event_init.elapsed(), 2);
127✔
439
    show_time("XS lookups", time_event_calculate_xs.elapsed(), 2);
127✔
440
    show_time("Advancing", time_event_advance_particle.elapsed(), 2);
127✔
441
    show_time("Surface crossings", time_event_surface_crossing.elapsed(), 2);
127✔
442
    show_time("Collisions", time_event_collision.elapsed(), 2);
127✔
443
    show_time("Particle death", time_event_death.elapsed(), 2);
127✔
444
  }
445
  if (settings::run_mode == RunMode::EIGENVALUE) {
4,208✔
446
    show_time("Time in inactive batches", time_inactive.elapsed(), 1);
2,538✔
447
  }
448
  show_time("Time in active batches", time_active.elapsed(), 1);
4,208✔
449
  if (settings::run_mode == RunMode::EIGENVALUE) {
4,208✔
450
    show_time("Time synchronizing fission bank", time_bank.elapsed(), 1);
2,538✔
451
    show_time("Sampling source sites", time_bank_sample.elapsed(), 2);
2,538✔
452
    show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2);
2,538✔
453
  }
454
  show_time("Time accumulating tallies", time_tallies.elapsed(), 1);
4,208✔
455
  show_time("Time writing statepoints", time_statepoint.elapsed(), 1);
4,208✔
456
  show_time("Total time for finalization", time_finalize.elapsed());
4,208✔
457
  show_time("Total time elapsed", time_total.elapsed());
4,208✔
458

459
  // Calculate particle rate in active/inactive batches
460
  int n_active = simulation::current_batch - settings::n_inactive;
4,208✔
461
  double speed_inactive = 0.0;
4,208✔
462
  double speed_active;
463
  if (settings::restart_run) {
4,208✔
464
    if (simulation::restart_batch < settings::n_inactive) {
33✔
UNCOV
465
      speed_inactive = (settings::n_particles *
×
UNCOV
466
                         (settings::n_inactive - simulation::restart_batch) *
×
UNCOV
467
                         settings::gen_per_batch) /
×
UNCOV
468
                       time_inactive.elapsed();
×
UNCOV
469
      speed_active =
×
470
        (settings::n_particles * n_active * settings::gen_per_batch) /
×
471
        time_active.elapsed();
×
472
    } else {
473
      speed_active = (settings::n_particles *
66✔
474
                       (settings::n_batches - simulation::restart_batch) *
33✔
475
                       settings::gen_per_batch) /
33✔
476
                     time_active.elapsed();
33✔
477
    }
478
  } else {
479
    if (settings::n_inactive > 0) {
4,175✔
480
      speed_inactive = (settings::n_particles * settings::n_inactive *
3,658✔
481
                         settings::gen_per_batch) /
1,829✔
482
                       time_inactive.elapsed();
1,829✔
483
    }
484
    speed_active =
4,175✔
485
      (settings::n_particles * n_active * settings::gen_per_batch) /
4,175✔
486
      time_active.elapsed();
4,175✔
487
  }
488

489
  // display calculation rate
490
  if (!(settings::restart_run &&
4,208✔
491
        (simulation::restart_batch >= settings::n_inactive)) &&
33✔
492
      settings::n_inactive > 0) {
4,175✔
493
    show_rate("Calculation Rate (inactive)", speed_inactive);
1,829✔
494
  }
495
  show_rate("Calculation Rate (active)", speed_active);
4,208✔
496
}
497

498
//==============================================================================
499

500
std::pair<double, double> mean_stdev(const double* x, int n)
13,744,652✔
501
{
502
  double mean = x[static_cast<int>(TallyResult::SUM)] / n;
13,744,652✔
503
  double stdev =
504
    n > 1 ? std::sqrt(std::max(0.0,
27,481,626✔
505
              (x[static_cast<int>(TallyResult::SUM_SQ)] / n - mean * mean) /
27,481,626✔
506
                (n - 1)))
13,736,974✔
507
          : 0.0;
13,744,652✔
508
  return {mean, stdev};
13,744,652✔
509
}
510

511
//==============================================================================
512

513
void print_results()
4,208✔
514
{
515
  // display header block for results
516
  header("Results", 4);
4,208✔
517
  if (settings::verbosity < 4)
4,208✔
UNCOV
518
    return;
×
519

520
  // Calculate t-value for confidence intervals
521
  int n = simulation::n_realizations;
4,208✔
522
  double alpha, t_n1, t_n3;
523
  if (settings::confidence_intervals) {
4,208✔
524
    alpha = 1.0 - CONFIDENCE_LEVEL;
11✔
525
    t_n1 = t_percentile(1.0 - alpha / 2.0, n - 1);
11✔
526
    t_n3 = t_percentile(1.0 - alpha / 2.0, n - 3);
11✔
527
  } else {
528
    t_n1 = 1.0;
4,197✔
529
    t_n3 = 1.0;
4,197✔
530
  }
531

532
  // write global tallies
533
  const auto& gt = simulation::global_tallies;
4,208✔
534
  double mean, stdev;
535
  if (n > 1) {
4,208✔
536
    if (settings::run_mode == RunMode::EIGENVALUE) {
4,084✔
537
      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_COLLISION, 0), n);
2,505✔
538
      fmt::print(" k-effective (Collision)     = {:.5f} +/- {:.5f}\n", mean,
2,049✔
539
        t_n1 * stdev);
2,505✔
540
      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_TRACKLENGTH, 0), n);
2,505✔
541
      fmt::print(" k-effective (Track-length)  = {:.5f} +/- {:.5f}\n", mean,
2,049✔
542
        t_n1 * stdev);
2,505✔
543
      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_ABSORPTION, 0), n);
2,505✔
544
      fmt::print(" k-effective (Absorption)    = {:.5f} +/- {:.5f}\n", mean,
2,049✔
545
        t_n1 * stdev);
2,505✔
546
      if (n > 3) {
2,505✔
547
        double k_combined[2];
548
        openmc_get_keff(k_combined);
2,450✔
549
        fmt::print(" Combined k-effective        = {:.5f} +/- {:.5f}\n",
2,450✔
550
          k_combined[0], k_combined[1]);
551
      }
552
    }
553
    std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::LEAKAGE, 0), n);
4,084✔
554
    fmt::print(
3,330✔
555
      " Leakage Fraction            = {:.5f} +/- {:.5f}\n", mean, t_n1 * stdev);
8,168✔
556
  } else {
557
    if (mpi::master)
124✔
558
      warning("Could not compute uncertainties -- only one "
124✔
559
              "active batch simulated!");
560

561
    if (settings::run_mode == RunMode::EIGENVALUE) {
124✔
562
      fmt::print(" k-effective (Collision)    = {:.5f}\n",
27✔
563
        gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n);
33✔
564
      fmt::print(" k-effective (Track-length) = {:.5f}\n",
27✔
565
        gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n);
33✔
566
      fmt::print(" k-effective (Absorption)   = {:.5f}\n",
27✔
567
        gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n);
66✔
568
    }
569
    fmt::print(" Leakage Fraction           = {:.5f}\n",
100✔
570
      gt(GlobalTally::LEAKAGE, TallyResult::SUM) / n);
248✔
571
  }
572
  fmt::print("\n");
3,430✔
573
  std::fflush(stdout);
4,208✔
574
}
575

576
//==============================================================================
577

578
const std::unordered_map<int, const char*> score_names = {
579
  {SCORE_FLUX, "Flux"},
580
  {SCORE_TOTAL, "Total Reaction Rate"},
581
  {SCORE_SCATTER, "Scattering Rate"},
582
  {SCORE_NU_SCATTER, "Scattering Production Rate"},
583
  {SCORE_ABSORPTION, "Absorption Rate"},
584
  {SCORE_FISSION, "Fission Rate"},
585
  {SCORE_NU_FISSION, "Nu-Fission Rate"},
586
  {SCORE_KAPPA_FISSION, "Kappa-Fission Rate"},
587
  {SCORE_EVENTS, "Events"},
588
  {SCORE_DECAY_RATE, "Decay Rate"},
589
  {SCORE_DELAYED_NU_FISSION, "Delayed-Nu-Fission Rate"},
590
  {SCORE_PROMPT_NU_FISSION, "Prompt-Nu-Fission Rate"},
591
  {SCORE_INVERSE_VELOCITY, "Flux-Weighted Inverse Velocity"},
592
  {SCORE_FISS_Q_PROMPT, "Prompt fission power"},
593
  {SCORE_FISS_Q_RECOV, "Recoverable fission power"},
594
  {SCORE_CURRENT, "Current"},
595
  {SCORE_PULSE_HEIGHT, "pulse-height"},
596
};
597

598
//! Create an ASCII output file showing all tally results.
599

600
void write_tallies()
4,802✔
601
{
602
  if (model::tallies.empty())
4,802✔
603
    return;
1,510✔
604

605
  // Set filename for tallies_out
606
  std::string filename = fmt::format("{}tallies.out", settings::path_output);
2,688✔
607

608
  // Open the tallies.out file.
609
  std::ofstream tallies_out;
3,292✔
610
  tallies_out.open(filename, std::ios::out | std::ios::trunc);
3,292✔
611

612
  // Loop over each tally.
613
  for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
21,650✔
614
    const auto& tally {*model::tallies[i_tally]};
18,358✔
615

616
    // Write header block.
617
    std::string tally_header("TALLY " + std::to_string(tally.id_));
18,358✔
618
    if (!tally.name_.empty())
18,358✔
619
      tally_header += ": " + tally.name_;
1,721✔
620
    fmt::print(tallies_out, "{}\n\n", header(tally_header));
36,716✔
621

622
    if (!tally.writable_) {
18,358✔
623
      fmt::print(tallies_out, " Internal\n\n");
330✔
624
      continue;
330✔
625
    }
626

627
    // Calculate t-value for confidence intervals
628
    double t_value = 1;
18,028✔
629
    if (settings::confidence_intervals) {
18,028✔
630
      auto alpha = 1 - CONFIDENCE_LEVEL;
11✔
631
      t_value = t_percentile(1 - alpha * 0.5, tally.n_realizations_ - 1);
11✔
632
    }
633

634
    // Write derivative information.
635
    if (tally.deriv_ != C_NONE) {
18,028✔
636
      const auto& deriv {model::tally_derivs[tally.deriv_]};
220✔
637
      switch (deriv.variable) {
220✔
638
      case DerivativeVariable::DENSITY:
88✔
639
        fmt::print(tallies_out, " Density derivative Material {}\n",
88✔
640
          deriv.diff_material);
88✔
641
        break;
88✔
642
      case DerivativeVariable::NUCLIDE_DENSITY:
88✔
643
        fmt::print(tallies_out,
88✔
644
          " Nuclide density derivative Material {} Nuclide {}\n",
645
          deriv.diff_material, data::nuclides[deriv.diff_nuclide]->name_);
88✔
646
        break;
88✔
647
      case DerivativeVariable::TEMPERATURE:
44✔
648
        fmt::print(tallies_out, " Temperature derivative Material {}\n",
44✔
649
          deriv.diff_material);
44✔
650
        break;
44✔
UNCOV
651
      default:
×
UNCOV
652
        fatal_error(fmt::format("Differential tally dependent variable for "
×
653
                                "tally {} not defined in output.cpp",
UNCOV
654
          tally.id_));
×
655
      }
656
    }
657

658
    // Initialize Filter Matches Object
659
    vector<FilterMatch> filter_matches;
36,056✔
660
    // Allocate space for tally filter matches
661
    filter_matches.resize(model::tally_filters.size());
18,028✔
662

663
    // Loop over all filter bin combinations.
664
    auto filter_iter = FilterBinIter(tally, false, &filter_matches);
18,028✔
665
    auto end = FilterBinIter(tally, true, &filter_matches);
18,028✔
666
    for (; filter_iter != end; ++filter_iter) {
9,450,719✔
667
      auto filter_index = filter_iter.index_;
9,432,691✔
668

669
      // Print info about this combination of filter bins.  The stride check
670
      // prevents redundant output.
671
      int indent = 0;
9,432,691✔
672
      for (auto i = 0; i < tally.filters().size(); ++i) {
29,776,007✔
673
        if (filter_index % tally.strides(i) == 0) {
20,343,316✔
674
          auto i_filt = tally.filters(i);
14,209,430✔
675
          const auto& filt {*model::tally_filters[i_filt]};
14,209,430✔
676
          auto& match {filter_matches[i_filt]};
14,209,430✔
677
          fmt::print(tallies_out, "{0:{1}}{2}\n", "", indent + 1,
28,418,860✔
678
            filt.text_label(match.i_bin_));
28,418,860✔
679
        }
680
        indent += 2;
20,343,316✔
681
      }
682

683
      // Loop over all nuclide and score combinations.
684
      int score_index = 0;
9,432,691✔
685
      for (auto i_nuclide : tally.nuclides_) {
20,293,622✔
686
        // Write label for this nuclide bin.
687
        if (i_nuclide == -1) {
10,860,931✔
688
          fmt::print(tallies_out, "{0:{1}}Total Material\n", "", indent + 1);
16,048,194✔
689
        } else {
690
          if (settings::run_CE) {
2,836,834✔
691
            fmt::print(tallies_out, "{0:{1}}{2}\n", "", indent + 1,
2,836,460✔
692
              data::nuclides[i_nuclide]->name_);
2,836,460✔
693
          } else {
694
            fmt::print(tallies_out, "{0:{1}}{2}\n", "", indent + 1,
374✔
695
              data::mg.nuclides_[i_nuclide].name);
374✔
696
          }
697
        }
698

699
        // Write the score, mean, and uncertainty.
700
        indent += 2;
10,860,931✔
701
        for (auto score : tally.scores_) {
24,593,984✔
702
          std::string score_name =
703
            score > 0 ? reaction_name(score) : score_names.at(score);
13,733,053✔
704
          double mean, stdev;
705
          std::tie(mean, stdev) =
13,733,053✔
706
            mean_stdev(&tally.results_(filter_index, score_index, 0),
13,733,053✔
707
              tally.n_realizations_);
27,466,106✔
708
          fmt::print(tallies_out, "{0:{1}}{2:<36} {3:.6} +/- {4:.6}\n", "",
13,733,053✔
709
            indent + 1, score_name, mean, t_value * stdev);
13,733,053✔
710
          score_index += 1;
13,733,053✔
711
        }
13,733,053✔
712
        indent -= 2;
10,860,931✔
713
      }
714
    }
715
  }
18,358✔
716
}
3,292✔
717

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