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

openmc-dev / openmc / 21116084326

18 Jan 2026 05:53PM UTC coverage: 82.04% (-0.004%) from 82.044%
21116084326

Pull #3732

github

web-flow
Merge 8a0f48648 into 5847b0de2
Pull Request #3732: Volume Calculation enhancement including refactoring and real-valued scoring implementation

17295 of 24049 branches covered (71.92%)

Branch coverage included in aggregate %.

302 of 306 new or added lines in 5 files covered. (98.69%)

96 existing lines in 3 files now uncovered.

55771 of 65012 relevant lines covered (85.79%)

43326460.86 hits per line

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

94.12
/src/volume_calc.cpp
1
#include "openmc/volume_calc.h"
2

3
#include "openmc/capi.h"
4
#include "openmc/cell.h"
5
#include "openmc/constants.h"
6
#include "openmc/distribution_multi.h"
7
#include "openmc/error.h"
8
#include "openmc/finalize.h"
9
#include "openmc/geometry.h"
10
#include "openmc/hdf5_interface.h"
11
#include "openmc/initialize.h"
12
#include "openmc/material.h"
13
#include "openmc/message_passing.h"
14
#include "openmc/mgxs_interface.h"
15
#include "openmc/nuclide.h"
16
#include "openmc/openmp_interface.h"
17
#include "openmc/output.h"
18
#include "openmc/random_dist.h"
19
#include "openmc/random_lcg.h"
20
#include "openmc/settings.h"
21
#include "openmc/timer.h"
22
#include "openmc/xml_interface.h"
23

24
#include "xtensor/xadapt.hpp"
25
#include "xtensor/xview.hpp"
26
#include <fmt/core.h>
27

28
#include <algorithm> // for copy
29
#include <cmath>     // for pow, sqrt
30
#include <unordered_set>
31

32
namespace openmc {
33

34
//==============================================================================
35
// Global variables
36
//==============================================================================
37

38
namespace model {
39
vector<VolumeCalculation> volume_calcs;
40
} // namespace model
41

42
//==============================================================================
43
// VolumeCalculation implementation
44
//==============================================================================
45

46
VolumeCalculation::VolumeCalculation(pugi::xml_node node)
358✔
47
{
48
  // Read domain type (cell, material or universe)
49
  std::string domain_type = get_node_value(node, "domain_type");
358✔
50
  if (domain_type == "cell") {
358✔
51
    domain_type_ = TallyDomain::CELL;
158✔
52
  } else if (domain_type == "material") {
200✔
53
    domain_type_ = TallyDomain::MATERIAL;
148✔
54
  } else if (domain_type == "universe") {
52!
55
    domain_type_ = TallyDomain::UNIVERSE;
52✔
56
  } else {
57
    fatal_error(std::string("Unrecognized domain type for stochastic "
×
58
                            "volume calculation: " +
59
                            domain_type));
60
  }
61

62
  // Read domain IDs, bounding corodinates and number of samples
63
  domain_ids_ = get_node_array<int>(node, "domain_ids");
358✔
64
  lower_left_ = get_node_array<double>(node, "lower_left");
358✔
65
  upper_right_ = get_node_array<double>(node, "upper_right");
358✔
66
  n_samples_ = std::stoull(get_node_value(node, "samples"));
358✔
67

68
  // Determine volume of bounding box
69
  const Position d {upper_right_ - lower_left_};
358✔
70
  volume_sample_ = d.x * d.y * d.z;
358✔
71

72
  if (check_for_node(node, "threshold")) {
358✔
73
    pugi::xml_node threshold_node = node.child("threshold");
128✔
74

75
    threshold_ = std::stod(get_node_value(threshold_node, "threshold"));
128✔
76
    if (threshold_ <= 0.0) {
128!
77
      fatal_error(fmt::format("Invalid error threshold {} provided for a "
×
78
                              "volume calculation.",
79
        threshold_));
×
80
    }
81

82
    // Prior calculation of the trigger volume fractions-related values t'
83
    // deriviated from the given volumes-related values t to prevent excessive
84
    // repeated computations during Monte Carlo execution. The values of t' are
85
    // computed via the sample mean \bar{x} and the adjusted sample variance
86
    // s^2.
87
    std::string tmp = get_node_value(threshold_node, "type");
128✔
88
    if (tmp == "variance") {
128✔
89
      trigger_type_ = TriggerMetric::variance;
32✔
90
      // Condition: s^2 < t' / N
91
      // t' = s^2 = t / (V_b)^2, in sq. volume fraction units
92
      threshold_cnd_ = threshold_ / std::pow(volume_sample_, 2);
32✔
93
      ;
94
    } else if (tmp == "std_dev") {
96✔
95
      trigger_type_ = TriggerMetric::standard_deviation;
32✔
96
      // Condition: s^2 < t' / N
97
      // t' = s^2 = (t / V_b)^2, in sq. volume fraction units
98
      threshold_cnd_ = std::pow(threshold_ / volume_sample_, 2);
32✔
99
    } else if (tmp == "rel_err") {
64!
100
      trigger_type_ = TriggerMetric::relative_error;
64✔
101
      // Condition: s^2 / \bar{x}^2 < t' / N
102
      // t' = s^2 / \bar{x}^2 = t^2, in relative units
103
      threshold_cnd_ = threshold_ * threshold_;
64✔
104
    } else {
105
      fatal_error(fmt::format(
×
106
        "Invalid volume calculation trigger type '{}' provided.", tmp));
107
    }
108

109
    if (check_for_node(threshold_node, "max_iterations")) {
128✔
110
      max_iterations_ =
32✔
111
        std::stoi(get_node_value(threshold_node, "max_iterations"));
32✔
112
      if (max_iterations_ <= 0) {
32!
NEW
113
        fatal_error(fmt::format(
×
NEW
114
          "Invalid error max_iterations {} provided.", max_iterations_));
×
115
      }
116
    }
117
  }
128✔
118

119
  // Ensure there are no duplicates by copying elements to a set and then
120
  // comparing the length with the original vector
121
  std::unordered_set<int> unique_ids(domain_ids_.cbegin(), domain_ids_.cend());
358✔
122
  if (unique_ids.size() != domain_ids_.size()) {
358!
123
    throw std::runtime_error {"Domain IDs for a volume calculation "
×
124
                              "must be unique."};
×
125
  }
126
}
358✔
127

128
void VolumeCalculation::execute(CalcResults& master_results) const
344✔
129
{
130
  // Check to make sure domain IDs are valid
131
  for (auto uid : domain_ids_) {
982✔
132
    switch (domain_type_) {
665!
133
    case TallyDomain::CELL:
365✔
134
      if (model::cell_map.find(uid) == model::cell_map.end()) {
365✔
135
        throw std::runtime_error {fmt::format(
18!
136
          "Cell {} in volume calculation does not exist in geometry.", uid)};
18✔
137
      }
138
      break;
356✔
139
    case TallyDomain::MATERIAL:
248✔
140
      if (model::material_map.find(uid) == model::material_map.end()) {
248✔
141
        throw std::runtime_error {fmt::format(
18!
142
          "Material {} in volume calculation does not exist in geometry.",
143
          uid)};
18✔
144
      }
145
      break;
239✔
146
    case TallyDomain::UNIVERSE:
52✔
147
      if (model::universe_map.find(uid) == model::universe_map.end()) {
52✔
148
        throw std::runtime_error {fmt::format(
18!
149
          "Universe {} in volume calculation does not exist in geometry.",
150
          uid)};
18✔
151
      }
152
    }
153
  }
154

155
  // Shared data that is collected from all threads
156
  const int n_domains = domain_ids_.size();
317✔
157

158
  // Divide work over MPI processes
159
  uint64_t min_samples = n_samples_ / mpi::n_procs;
317✔
160
  uint64_t remainder = n_samples_ % mpi::n_procs;
317✔
161
  uint64_t i_start, i_end;
162
  if (mpi::rank < remainder) {
317!
163
    i_start = (min_samples + 1) * mpi::rank;
×
164
    i_end = i_start + min_samples + 1;
×
165
  } else {
166
    i_start =
317✔
167
      (min_samples + 1) * remainder + (mpi::rank - remainder) * min_samples;
317✔
168
    i_end = i_start + min_samples;
317✔
169
  }
170

171
  while (true) {
172

173
#pragma omp parallel
11,735✔
174
    {
175
      // Temporary variables that are private to each thread
176
      CalcResults results(*this);
9,126✔
177
      results.sampling_time.start();
9,126✔
178
      Particle p;
9,126✔
179

180
      // Sample locations and count scores
181
#pragma omp for
182
      for (size_t i = i_start; i < i_end; i++) {
3,690,176✔
183
        uint64_t id = master_results.iterations * n_samples_ + i;
3,681,050✔
184
        uint64_t seed = init_seed(id, STREAM_VOLUME);
3,681,050✔
185

186
        p.n_coord() = 1;
3,681,050✔
187
        p.r() = {uniform_distribution(lower_left_.x, upper_right_.x, &seed),
7,362,100✔
188
          uniform_distribution(lower_left_.y, upper_right_.y, &seed),
3,681,050✔
189
          uniform_distribution(lower_left_.z, upper_right_.z, &seed)};
3,681,050✔
190
        constexpr double sqrt3_1 = 1. / std::sqrt(3.);
3,681,050✔
191
        p.u() = {sqrt3_1, sqrt3_1, sqrt3_1};
3,681,050✔
192

193
        // TO REVIEWER: THE SWITCH IS TRANSFERED TO score_hit()
194
        if (exhaustive_find_cell(p))
3,681,050✔
195
          this->score_hit(p, results);
2,733,235✔
196

197
        results.n_samples++;
3,681,050✔
198

199
        // This passing across all tallies after each sample can be
200
        // inefficient for the case of large number of domains and small
201
        // size of batch, but the batch size is currently assigned to 1
202
        // for keep the input format being unchanged
203
        const double batch_size_1 = 1. / static_cast<double>(1);
3,681,050✔
204
        for (auto& vt : results.vol_tallies) {
11,636,200✔
205
          for (auto& vol_tally : vt) {
24,517,551✔
206
            vol_tally.finalize_batch(batch_size_1);
16,562,401✔
207
          }
208
        }
209
      } // sample/batch loop
210

211
      results.sampling_time.stop();
9,126✔
212
      results.cost = results.sampling_time.elapsed();
9,126✔
213

214
      // At this point, each thread has its own volume tallies lists and we
215
      // now need to reduce them. OpenMP is not nearly smart enough to do
216
      // this on its own, so we have to manually reduce them
217
      reduce_results(results, master_results);
9,126✔
218
    } // omp parallel
9,126✔
219

220
    // bump iteration counter
221
    master_results.iterations++;
20,861✔
222

223
    // TO REVIEWER: THE MPI-RELATED PART IS MOVED TO collect_MPI()
224
#ifdef OPENMC_MPI
225
    master_results.collect_MPI(); // collect results to master process
13,025✔
226
#endif
227
    // Process volume estimation results in master process for the trigger state
228
    // determination
229
    bool stop_calc =
20,861✔
230
      mpi::master && (trigger_type_ == TriggerMetric::not_active ||
35,073!
231
                       master_results.iterations == max_iterations_);
14,212✔
232

233
    if (!stop_calc) {
20,861✔
234
      // Compute current trigger state across totals (0th elements) only
235
      for (const auto& vt : master_results.vol_tallies) {
20,916✔
236
        stop_calc = vt[0].trigger_state(
20,820✔
237
          trigger_type_, threshold_cnd_, master_results.n_samples);
20,820✔
238
        if (!stop_calc)
20,820✔
239
          break;
20,584✔
240
      }
241
    }
242

243
#ifdef OPENMC_MPI
244
    // Send the state of calculation continuation just obtained in master
245
    // process to all processes
246
    MPI_Bcast(&stop_calc, 1, MPI_CXX_BOOL, 0, mpi::intracomm);
13,025✔
247
#endif
248

249
    if (!stop_calc)
20,861✔
250
      continue; // while loop
20,544✔
251
    // No trigger is applied or the trigger condition is satisfied, we're
252
    // done
253

254
    if (mpi::master) {
317✔
255
      // Normalize all results on the bounding primitive volume and compute
256
      // stddev
257
      for (auto& vt : master_results.vol_tallies) {
725✔
258
        for (auto& vol_tally : vt) {
1,488✔
259
          vol_tally.score_acc = get_tally_results(
1,010✔
260
            master_results.n_samples, volume_sample_, vol_tally);
1,010✔
261
        }
262
      }
263

264
      // Compute nuclides
265
      for (int i_domain = 0; i_domain < n_domains; ++i_domain) {
725✔
266
        // Create 2D array to store atoms/uncertainty for each nuclide. Later,
267
        // this is compressed into vectors storing only those nuclides that are
268
        // non-zero
269
        auto n_nuc =
270
          settings::run_CE ? data::nuclides.size() : data::mg.nuclides_.size();
478✔
271
        xt::xtensor<double, 2> atoms({n_nuc, 2}, 0.0);
478✔
272

273
        for (int j = 0; j < master_results.vol_tallies[i_domain].size(); ++j) {
1,488✔
274
          const int i_material = master_results.vol_tallies[i_domain][j].index;
1,010✔
275
          if (i_material == MATERIAL_VOID || i_material == _INDEX_TOTAL)
1,010✔
276
            continue;
489✔
277

278
          const auto& mat = model::materials[i_material];
521✔
279
          for (int k = 0; k < mat->nuclide_.size(); ++k) {
2,040✔
280
            auto& volume = master_results.vol_tallies[i_domain][j].score_acc;
1,519✔
281
            // Collect calculated nuclide amounts N [atoms] and stddev as
282
            // N = V [cm^3] * \rho [atoms/b-cm] * 1.e24 [b-cm/cm^3]
283
            atoms(mat->nuclide_[k], 0) +=
1,519✔
284
              volume[0] * mat->atom_density(k) * 1.0e24;
1,519✔
285
            atoms(mat->nuclide_[k], 1) +=
1,519✔
286
              volume[1] * mat->atom_density(k) * 1.0e24;
1,519✔
287
          }
288
        }
289

290
        // Get reference to result for this domain
291
        auto& result {master_results.nuc_results[i_domain]};
478✔
292
        // Convert full arrays to vectors
293
        for (int j = 0; j < n_nuc; ++j) {
3,440✔
294
          if (atoms(j, 0) > 0.0) {
2,962✔
295
            result.nuclides.push_back(j);
1,399✔
296
            result.atoms.push_back(atoms(j, 0));
1,399✔
297
            result.uncertainty.push_back(atoms(j, 1));
1,399✔
298
          }
299
        }
300
      } // end domains loop
478✔
301
    }
302

303
    return;
317✔
304

305
  } // end while
20,544✔
306
}
307

308
void VolumeCalculation::show_vol_stat(
741✔
309
  const std::string label, const std::string units, const double value) const
310
{
311
  fmt::print("{0:<20} = {2:10.4e} {1:<}\n", label, units, value);
741✔
312
}
741✔
313

314
void VolumeCalculation::show_volume(const std::string domain_type,
478✔
315
  const int domain_id, const std::string region_name, const double mean,
316
  const double stddev) const
317
{
318
  fmt::print(" {0:>9}{1:>6}: {2:10.4e} +/- {3:10.4e} cm^3", domain_type,
394✔
319
    domain_id, mean, stddev);
320
  if (!region_name.empty()) {
478✔
321
    fmt::print("   //{:<}", region_name);
80✔
322
  }
323
  fmt::print("\n");
478✔
324
}
478✔
325

326
// TO REVIEWER: I/O BLOCK BEGINS, I WOULD PREFFER TO PUT IT AT THE END OF FILE
327
// BUT KEEPT THE FILE STRUCTURE UNCHANGED
328
void VolumeCalculation::show_results(const CalcResults& results) const
247✔
329
{
330
  // Show tracing statistics
331
  write_message(5, " ");
247✔
332
  show_vol_stat(
247✔
333
    "Total sample size", "", static_cast<double>(results.n_samples));
247✔
334
  show_vol_stat("Running cost", "thread-sec", results.cost);
247✔
335
  show_vol_stat("Cost of hitting", "thread-sec/hit",
247✔
336
    static_cast<double>(results.cost) / static_cast<double>(results.n_samples));
247✔
337
  write_message(5, " ");
247✔
338

339
  std::string domain_type;
247✔
340
  if (domain_type_ == TallyDomain::CELL) {
247✔
341
    domain_type = "Cell";
112✔
342
  } else if (domain_type_ == TallyDomain::MATERIAL) {
135✔
343
    domain_type = "Material";
102✔
344
  } else {
345
    domain_type = "Universe";
33✔
346
  }
347

348
  // Display domain volumes
349
  for (int j = 0; j < domain_ids_.size(); j++) {
725✔
350
    std::string region_name {""};
478✔
351
    if (domain_type_ == TallyDomain::CELL) {
478✔
352
      int cell_idx = model::cell_map[domain_ids_[j]];
266✔
353
      region_name = model::cells[cell_idx]->name();
266✔
354
    } else if (domain_type_ == TallyDomain::MATERIAL) {
212✔
355
      int mat_idx = model::material_map[domain_ids_[j]];
179✔
356
      region_name = model::materials[mat_idx]->name();
179✔
357
    }
358
    if (region_name.size())
478✔
359
      region_name.insert(0, " "); // prepend space for formatting
80✔
360

361
    show_volume(domain_type, domain_ids_[j], region_name,
478✔
362
      results.vol_tallies[j][0].score_acc[0],
478✔
363
      results.vol_tallies[j][0].score_acc[1]);
478✔
364
  }
478✔
365
  write_message(4, " "); // Blank line afer results printed
247✔
366
}
247✔
367

368
void VolumeCalculation::to_hdf5(
247✔
369
  const std::string& filename, const CalcResults& results) const
370
{
371
  // Create HDF5 file
372
  hid_t file_id = file_open(filename, 'w');
247✔
373

374
  // Write header info
375
  write_attribute(file_id, "filetype", "volume");
247✔
376
  write_attribute(file_id, "version", VERSION_VOLUME);
247✔
377
  write_attribute(file_id, "openmc_version", VERSION);
247✔
378
#ifdef GIT_SHA1
379
  write_attribute(file_id, "git_sha1", GIT_SHA1);
380
#endif
381

382
  // Write current date and time
383
  write_attribute(file_id, "date_and_time", time_stamp());
247✔
384

385
  // Write basic metadata
386
  write_attribute(file_id, "samples", n_samples_);
247✔
387
  write_attribute(file_id, "lower_left", lower_left_);
247✔
388
  write_attribute(file_id, "upper_right", upper_right_);
247✔
389
  // Write trigger info
390
  if (trigger_type_ != TriggerMetric::not_active) {
247✔
391
    write_attribute(file_id, "iterations", results.iterations);
88✔
392
    write_attribute(file_id, "threshold", threshold_);
88✔
393
    std::string trigger_str;
88✔
394
    switch (trigger_type_) {
88!
395
    case TriggerMetric::variance:
22✔
396
      trigger_str = "variance";
22✔
397
      break;
22✔
398
    case TriggerMetric::standard_deviation:
22✔
399
      trigger_str = "std_dev";
22✔
400
      break;
22✔
401
    case TriggerMetric::relative_error:
44✔
402
      trigger_str = "rel_err";
44✔
403
      break;
44✔
404
    default:
×
405
      break;
×
406
    }
407
    write_attribute(file_id, "trigger_type", trigger_str);
88✔
408
    // check max_iterations on default value
409
    if (max_iterations_ !=
88✔
410
        std::numeric_limits<decltype(max_iterations_)>::max())
88✔
411
      write_attribute(file_id, "max_iterations", max_iterations_);
22✔
412
  } else {
88✔
413
    write_attribute(file_id, "iterations", 1);
159✔
414
  }
415

416
  if (domain_type_ == TallyDomain::CELL) {
247✔
417
    write_attribute(file_id, "domain_type", "cell");
112✔
418
  } else if (domain_type_ == TallyDomain::MATERIAL) {
135✔
419
    write_attribute(file_id, "domain_type", "material");
102✔
420
  } else if (domain_type_ == TallyDomain::UNIVERSE) {
33!
421
    write_attribute(file_id, "domain_type", "universe");
33✔
422
  }
423

424
  for (int i = 0; i < domain_ids_.size(); ++i) {
725✔
425
    hid_t group_id =
426
      create_group(file_id, fmt::format("domain_{}", domain_ids_[i]));
956✔
427

428
    // Write volume for domain
429
    const auto& result {results.nuc_results[i]};
478✔
430
    write_dataset(group_id, "volume", results.vol_tallies[i][0].score_acc);
478✔
431

432
    // Create array of nuclide names from the vector
433
    auto n_nuc = result.nuclides.size();
478✔
434

435
    vector<std::string> nucnames;
478✔
436
    for (int i_nuc : result.nuclides) {
1,877✔
437
      nucnames.push_back(settings::run_CE ? data::nuclides[i_nuc]->name_
1,883✔
438
                                          : data::mg.nuclides_[i_nuc].name);
484✔
439
    }
440

441
    // Create array of total # of atoms with uncertainty for each nuclide
442
    xt::xtensor<double, 2> atom_data({n_nuc, 2});
478✔
443
    xt::view(atom_data, xt::all(), 0) = xt::adapt(result.atoms);
478✔
444
    xt::view(atom_data, xt::all(), 1) = xt::adapt(result.uncertainty);
478✔
445

446
    // Write results
447
    write_dataset(group_id, "nuclides", nucnames);
478✔
448
    write_dataset(group_id, "atoms", atom_data);
478✔
449

450
    close_group(group_id);
478✔
451
  }
478✔
452

453
  file_close(file_id);
247✔
454
}
247✔
455

456
// TO REVIEWER: THIS IS THE SAME check_hit()
457
void VolumeCalculation::check_hit(const int32_t i_material,
10,974,551✔
458
  const double contrib, vector<VolTally>& vol_tallies) const
459
{
460
  // Contribute to entire domain result tally
461
  vol_tallies[0].score += contrib;
10,974,551✔
462

463
  // Check if this material was previously hit and if so, contribute score
464
  for (int j = 1; j < vol_tallies.size(); j++) {
11,623,251✔
465
    if (vol_tallies[j].index == i_material) {
11,538,587✔
466
      vol_tallies[j].score += contrib;
10,889,887✔
467
      return;
10,889,887✔
468
    }
469
  }
470

471
  // The material was not previously hit, append an entry to the material
472
  // indices and hits lists
473
  vol_tallies.push_back(VolTally(i_material, contrib));
84,664✔
474
}
475

476
// TO REVIEWER: THIS IS AN EQUIVALENT OF THE reduce_indicies_hits() TEMPLATE
477
// FROM volume_calc.h
478
void VolumeCalculation::reduce_results(
32,596✔
479
  const CalcResults& local_results, CalcResults& results) const
480
{
481
  auto n_threads = num_threads();
32,596✔
482

483
  // Collect scalar variables
484
#pragma omp for ordered schedule(static, 1)
23,470!
485
  for (int i = 0; i < n_threads; ++i) {
18,252✔
486
#pragma omp ordered
487
    {
488
      results.n_samples += local_results.n_samples;
32,596✔
489
      results.cost += local_results.cost;
32,596!
490
    }
491
  }
492

493
  // Collect vectored domain-wise results
494
  for (int i_domain = 0; i_domain < domain_ids_.size(); ++i_domain) {
129,098✔
495
    const vector<VolTally>& local_vol_tall =
496
      local_results.vol_tallies[i_domain];
96,502✔
497
    vector<VolTally>& vol_tall = results.vol_tallies[i_domain];
96,502✔
498

499
#pragma omp for ordered schedule(static, 1)
69,488!
500
    for (int i = 0; i < n_threads; ++i) {
54,028✔
501
#pragma omp ordered
502
      for (int j = 0; j < local_vol_tall.size(); ++j) {
277,668✔
503
        // Check if this material has been added to the master list and if
504
        // so, accumulate scores
505
        const auto ind {local_vol_tall[j].index};
181,166✔
506
        const auto it = std::find_if(vol_tall.begin(), vol_tall.end(),
181,166✔
507
          [ind](const VolTally& vt) { return vt.index == ind; });
13,879,145✔
508
        if (it == vol_tall.end()) {
181,166✔
509
          vol_tall.push_back(local_vol_tall[j]);
26,162✔
510
        } else {
511
          vol_tall[it - vol_tall.begin()].append_tally(local_vol_tall[j]);
155,004✔
512
        }
513
      }
514
    }
515
  }
516
}
32,596✔
517

518
#ifdef OPENMC_MPI
519
// TO REVIEWER: THIS IS AN EQUIVALENT OF THE RELATED MPI-INTERCHANGE BLOCK OF
520
// execute()
521
void VolumeCalculation::CalcResults::collect_MPI()
13,025✔
522
{
523
  vector<int> domain_sizes(vol_tallies.size());
13,025✔
524

525
  if (!mpi::master) {
13,025✔
526
    // n_domain + 2 MPI messages will be send in total to node mpi::master
527

528
    // To determination of an unique tag for each MPI message (as below)
529
    int mpi_offset = mpi::rank * (vol_tallies.size() + 2) + 2;
6,490✔
530

531
    // Pass root data of the struct
532
    MPI_Send(
6,490✔
533
      (void*)this, 1, mpi::volume_results, 0, mpi_offset - 2, mpi::intracomm);
534

535
    // Pass sizes of domain-wise data
536
    for (int i_domain = 0; i_domain < vol_tallies.size(); i_domain++) {
25,750✔
537
      domain_sizes[i_domain] = vol_tallies[i_domain].size();
19,260✔
538
    }
539

540
    MPI_Send(domain_sizes.data(), domain_sizes.size(), MPI_INT, 0,
6,490✔
541
      mpi_offset - 1, mpi::intracomm);
542

543
    // Pass domain-wise data of struct
544
    for (int i_domain = 0; i_domain < vol_tallies.size(); i_domain++) {
25,750✔
545
      MPI_Send(vol_tallies[i_domain].data(), domain_sizes[i_domain],
19,260✔
546
        mpi::volume_tally, 0, mpi_offset + i_domain, mpi::intracomm);
547
    }
548

549
    this->reset(); // Delete passed to main process data
6,490✔
550

551
  } else {
552

553
    // n_domain + 2 MPI messages will be recieved in total from node mpi::master
554

555
    for (int i_proc = 1; i_proc < mpi::n_procs; i_proc++) {
13,025✔
556

557
      CalcResults res_buff(*this); // temporary storage for recived data
6,490✔
558

559
      // To determination of an unique tag for each MPI message (as above)
560
      int mpi_offset = i_proc * (vol_tallies.size() + 2) + 2;
6,490✔
561

562
      // Pass root data of struct
563
      MPI_Recv(&res_buff, 1, mpi::volume_results, i_proc, mpi_offset - 2,
6,490✔
564
        mpi::intracomm, MPI_STATUS_IGNORE);
565

566
      // Pass sizes of domain-wise data
567
      MPI_Recv(domain_sizes.data(), domain_sizes.size(), MPI_INT, i_proc,
6,490✔
568
        mpi_offset - 1, mpi::intracomm, MPI_STATUS_IGNORE);
569

570
      // Pass domain-wise data of struct
571
      for (int i_domain = 0; i_domain < vol_tallies.size(); i_domain++) {
25,750✔
572
        res_buff.vol_tallies[i_domain].resize(domain_sizes[i_domain]);
19,260✔
573
        MPI_Recv(res_buff.vol_tallies[i_domain].data(), domain_sizes[i_domain],
19,260✔
574
          mpi::volume_tally, i_proc, mpi_offset + i_domain, mpi::intracomm,
575
          MPI_STATUS_IGNORE);
576
      }
577

578
      this->append(res_buff);
6,490✔
579
    }
6,490✔
580
  }
581
}
13,025✔
582
#endif
583

584
// TO REVIEWER: STATISTICS AND TALLY MANIPULATION BLOCK BEGIN
585
VolumeCalculation::CalcResults::CalcResults(const VolumeCalculation& vol_calc)
32,940✔
586
{
587
  n_samples = 0;
32,940✔
588
  iterations = 0;
32,940✔
589
  cost = 0.;
32,940✔
590
  for (int i = 0; i < vol_calc.domain_ids_.size(); i++) {
130,107✔
591
    vector<VolTally> vt_vect; // Tally group for a domain
97,167✔
592
    vt_vect.push_back(
97,167✔
593
      VolTally(_INDEX_TOTAL)); // Zero-element tally for entire domain totals
97,167✔
594
    vol_tallies.push_back(vt_vect);
97,167✔
595
    nuc_results.push_back(NuclResult());
97,167✔
596
  }
97,167✔
597
}
32,940✔
598

599
void VolumeCalculation::CalcResults::reset()
6,490✔
600
{
601
  n_samples = 0;
6,490✔
602
  cost = 0.;
6,490✔
603
  for (auto& vt : vol_tallies) {
25,750✔
604
    std::fill(vt.begin(), vt.end(), VolTally());
19,260!
605
  }
606
}
6,490✔
607

608
// TO REVIEWER: THIS IS A PART OF THE RELATED TO MPI-INTERCHANGE
609
// BLOCK OF execute()
610
void VolumeCalculation::CalcResults::append(const CalcResults& other)
6,490✔
611
{
612
  n_samples += other.n_samples;
6,490✔
613
  cost += other.cost;
6,490✔
614

615
  // The domain-wise vectors this.vol_tallies and other.vol_tallies are
616
  // conformed each to other by definition
617
  for (auto id = 0; id < vol_tallies.size(); id++) {
25,750✔
618
    // Merging current domain vector from other.vol_tallies into this via
619
    // pair-wise comparisons
620
    for (const auto& vt_other : other.vol_tallies[id]) {
5,738,470✔
621
      bool already_appended = false;
5,719,210✔
622
      for (auto& vt : vol_tallies[id]) {
14,259,870✔
623
        if (vt.index == vt_other.index) {
14,259,810✔
624
          vt.append_tally(vt_other);
5,719,150✔
625
          already_appended = true;
5,719,150✔
626
          break;
5,719,150✔
627
        }
628
      }
629
      if (!already_appended)
5,719,210✔
630
        vol_tallies[id].push_back(vt_other);
60!
631
    }
632
  }
633
}
6,490✔
634

635
inline void VolumeCalculation::VolTally::finalize_batch(
60,122,709✔
636
  const double batch_size_1)
637
{
638
  if (score != 0.) {
60,122,709✔
639
    score *= batch_size_1;
21,949,102✔
640
    score_acc[0] += score;
21,949,102✔
641
    score_acc[1] += score * score;
21,949,102✔
642
    score = 0.;
21,949,102✔
643
  }
644
}
60,122,709✔
645

646
inline void VolumeCalculation::VolTally::assign_tally(const VolTally& vol_tally)
647
{
648
  score = vol_tally.score;
649
  score_acc = vol_tally.score_acc;
650
  index = vol_tally.index;
651
}
652

653
inline void VolumeCalculation::VolTally::append_tally(const VolTally& vol_tally)
5,874,154✔
654
{
655
  score += vol_tally.score;
5,874,154✔
656
  score_acc[0] += vol_tally.score_acc[0];
5,874,154✔
657
  score_acc[1] += vol_tally.score_acc[1];
5,874,154✔
658
}
5,874,154✔
659

660
// TO REVIEWER: IF THIS APPROACH WILL BE FOUND OBFUSCATED, IT CAN BE ADJUSTED TO
661
// LITERAL FORMULAE TRANSLATION WITH REPEATING MULTIPLY CALCULATIONS OF
662
// THRESHOLD SQUARE ROOTS AND DIVISIONS FOR EACH VOLUME DOMAIN IN THE END OF
663
// EACH BATCH
664
inline bool VolumeCalculation::VolTally::trigger_state(
20,820✔
665
  const TriggerMetric trigger_type, const double threshold,
666
  const size_t& n_samples) const
667
{
668
  // For sample contribution to volume fraction limited by 1, the maximal
669
  // allowed n_samples value is around 1.e102, but this is still much larger
670
  // than the size_t limit equal to ~1.8e19
671
  const double ns1 = static_cast<double>(n_samples - 1);
20,820✔
672
  const double ns = static_cast<double>(n_samples);
20,820✔
673
  const double mean_xi_sq = score_acc[0] * score_acc[0];
20,820✔
674

675
  // Adjusted sample variance: s^2 = (\bar{x^2} - \bar{x}^2) / (N-1)
676
  // \bar{x}^2 = mean_xi_sq / N^2, \bar{\x^2} = score_acc[1] / N, N = ns
677
  switch (trigger_type) {
20,820✔
678
  case TriggerMetric::variance:
20,184✔
679
  case TriggerMetric::standard_deviation:
680
    // Condition: s^2 / N < t'
681
    // Equivalent implementation:
682
    // N^2 * (\bar{x^2} - \bar{x}^2) < t' * (N-1) * N^2
683
    return score_acc[1] * ns - mean_xi_sq < threshold * ns1 * ns * ns;
20,184✔
684
  case TriggerMetric::relative_error:
576✔
685
    // Condition: (s^2 / \mu^2) / N < t'
686
    // Equivalent implementation:
687
    // N^2 * (\bar{x^2} - \bar{x}^2) < t' * (N-1) * (N * \bar{x})^2
688
    return score_acc[1] * ns - mean_xi_sq < threshold * ns1 * mean_xi_sq;
576✔
689
  default:
60✔
690
    return true;
60✔
691
  }
692
}
693

694
array<double, 2> VolumeCalculation::get_tally_results(const size_t& n_samples,
1,010✔
695
  const double coeff_norm, const VolTally& vol_tally) const
696
{
697
  array<double, 2> volume;
698
  const double ns_1 = 1. / static_cast<double>(n_samples);
1,010✔
699
  volume[0] = vol_tally.score_acc[0] * ns_1;
1,010✔
700
  volume[1] = vol_tally.score_acc[1] * ns_1;
1,010✔
701
  volume[1] = std::sqrt(
2,020✔
702
    (volume[1] - volume[0] * volume[0]) / static_cast<double>(n_samples - 1));
1,010✔
703
  volume[0] *= coeff_norm;
1,010✔
704
  volume[1] *= coeff_norm;
1,010✔
705
  return volume;
1,010✔
706
}
707

708
// TO REVIEWER: THIS IS A FULL EQUIVALENT OF THE SWITCH FROM THE HISTORY LOOP
709
void VolumeCalculation::score_hit(const Particle& p, CalcResults& results) const
11,013,117✔
710
{
711
  const auto n_domains = domain_ids_.size();
11,013,117✔
712
  const auto id_mat = p.material();
11,013,117✔
713
  const auto score = 1.; // Floating-point score value
11,013,117✔
714

715
  switch (domain_type_) {
11,013,117!
716
  case VolumeCalculation::TallyDomain::MATERIAL:
4,689,853✔
717
    if (id_mat != MATERIAL_VOID) {
4,689,853✔
718
      for (auto i_domain = 0; i_domain < n_domains; i_domain++) {
7,653,891✔
719
        if (model::materials[id_mat]->id_ == domain_ids_[i_domain]) {
7,644,739✔
720
          this->check_hit(id_mat, score, results.vol_tallies[i_domain]);
4,675,718✔
721
          break;
4,675,718✔
722
        }
723
      }
724
    }
725
    break;
4,689,853✔
726
  case VolumeCalculation::TallyDomain::CELL:
4,684,429✔
727
    for (auto level = 0; level < p.n_coord_last(); ++level) {
9,368,858✔
728
      for (auto i_domain = 0; i_domain < n_domains; i_domain++) {
5,190,583✔
729
        if (model::cells[p.coord(level).cell()]->id_ == domain_ids_[i_domain]) {
5,166,152✔
730
          this->check_hit(id_mat, score, results.vol_tallies[i_domain]);
4,659,998✔
731
          break;
4,659,998✔
732
        }
733
      }
734
    }
735
    break;
4,684,429✔
736
  case VolumeCalculation::TallyDomain::UNIVERSE:
1,638,835✔
737
    for (auto level = 0; level < p.n_coord_last(); ++level) {
3,277,670✔
738
      for (auto i_domain = 0; i_domain < n_domains; ++i_domain) {
1,638,835!
739
        if (model::universes[model::cells[p.coord(level).cell()]->universe_]
1,638,835✔
740
              ->id_ == domain_ids_[i_domain]) {
1,638,835!
741
          this->check_hit(id_mat, score, results.vol_tallies[i_domain]);
1,638,835✔
742
          break;
1,638,835✔
743
        }
744
      }
745
    }
746
  }
747
}
11,013,117✔
748

749
void free_memory_volume()
8,050✔
750
{
751
  openmc::model::volume_calcs.clear();
8,050✔
752
}
8,050✔
753

754
} // namespace openmc
755

756
//==============================================================================
757
// OPENMC_CALCULATE_VOLUMES runs each of the stochastic volume calculations
758
// that the user has specified and writes results to HDF5 files
759
//==============================================================================
760

761
int openmc_calculate_volumes()
128✔
762
{
763
  using namespace openmc;
764

765
  if (mpi::master) {
128✔
766
    header("STOCHASTIC VOLUME CALCULATION", 3);
118✔
767
  }
768
  Timer time_volume;
128✔
769
  time_volume.start();
128✔
770

771
  for (int i = 0; i < model::volume_calcs.size(); ++i) {
445✔
772
    write_message(4, "Running volume calculation {}", i + 1);
344✔
773

774
    // Run volume calculation
775
    const auto& vol_calc {model::volume_calcs[i]};
344✔
776
    VolumeCalculation::CalcResults results(vol_calc);
344✔
777

778
    try {
779
      vol_calc.execute(results);
344✔
780
    } catch (const std::exception& e) {
27!
781
      set_errmsg(e.what());
27✔
782
      return OPENMC_E_UNASSIGNED;
27✔
783
    }
27✔
784

785
    if (mpi::master) {
317✔
786

787
      // Output volume calculation results and statistics
788
      vol_calc.show_results(results);
247✔
789

790
      // Write volumes to HDF5 file
791
      std::string filename =
792
        fmt::format("{}volume_{}.h5", settings::path_output, i + 1);
450✔
793
      vol_calc.to_hdf5(filename, results);
247✔
794
    }
247✔
795
  }
344✔
796

797
  // Show elapsed time
798
  time_volume.stop();
101✔
799
  write_message(6, "Elapsed time: {} sec", time_volume.elapsed());
101✔
800

801
  return 0;
101✔
802
}
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