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

openmc-dev / openmc / 19830489469

01 Dec 2025 04:50PM UTC coverage: 82.12% (+0.05%) from 82.071%
19830489469

Pull #3645

github

web-flow
Merge 2329adead into ef22558f4
Pull Request #3645: Ray tracing volume calculation and geometry testing

17085 of 23668 branches covered (72.19%)

Branch coverage included in aggregate %.

418 of 429 new or added lines in 4 files covered. (97.44%)

120 existing lines in 8 files now uncovered.

55142 of 64285 relevant lines covered (85.78%)

42020645.37 hits per line

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

93.2
/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/geometry.h"
9
#include "openmc/hdf5_interface.h"
10
#include "openmc/material.h"
11
#include "openmc/message_passing.h"
12
#include "openmc/mgxs_interface.h"
13
#include "openmc/nuclide.h"
14
#include "openmc/openmp_interface.h"
15
#include "openmc/output.h"
16
#include "openmc/plot.h"
17
#include "openmc/random_dist.h"
18
#include "openmc/random_lcg.h"
19
#include "openmc/settings.h"
20
#include "openmc/timer.h"
21
#include "openmc/xml_interface.h"
22

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

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

31
namespace openmc {
32

33
//==============================================================================
34
// Global variables
35
//==============================================================================
36

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

41
#ifdef OPENMC_MPI
42
static MPI_Datatype mpi_vol_results;  //!< MPI struct for CalcResults
43
static MPI_Datatype mpi_volume_tally; //!< MPI struct for VolTally
44
#endif
45

46
//==============================================================================
47
// VolumeCalculation implementation
48
//==============================================================================
49

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

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

72
  // Determine volume of bounding box
73
  const Position d {upper_right_ - lower_left_};
358✔
74
  volume_sample_ = d.x * d.y * d.z;
358✔
75

76
  if (check_for_node(node, "threshold")) {
358✔
77
    pugi::xml_node threshold_node = node.child("threshold");
128✔
78

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

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

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

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

131
  // Type of volume estimator
132
  std::string tmp = get_node_value(node, "estimator_type");
358✔
133
  if (tmp == "hit") {
358✔
134
    mode_ = EstMode::REJECTION;
262✔
135
  } else if (tmp == "ray") {
96!
136
    mode_ = EstMode::RAYTRACE;
96✔
137
  } else {
NEW
138
    fatal_error(
×
NEW
139
      fmt::format("Invalid volume calculation mode '{}' provided.", tmp));
×
140
  }
141
}
358✔
142

143
void VolumeCalculation::execute(CalcResults& master_results) const
344✔
144
{
145
  // Check to make sure domain IDs are valid
146
  for (auto uid : domain_ids_) {
982✔
147
    switch (domain_type_) {
665!
148
    case TallyDomain::CELL:
365✔
149
      if (model::cell_map.find(uid) == model::cell_map.end()) {
365✔
150
        throw std::runtime_error {fmt::format(
18!
151
          "Cell {} in volume calculation does not exist in geometry.", uid)};
18✔
152
      }
153
      break;
356✔
154
    case TallyDomain::MATERIAL:
248✔
155
      if (model::material_map.find(uid) == model::material_map.end()) {
248✔
156
        throw std::runtime_error {fmt::format(
18!
157
          "Material {} in volume calculation does not exist in geometry.",
158
          uid)};
18✔
159
      }
160
      break;
239✔
161
    case TallyDomain::UNIVERSE:
52✔
162
      if (model::universe_map.find(uid) == model::universe_map.end()) {
52✔
163
        throw std::runtime_error {fmt::format(
18!
164
          "Universe {} in volume calculation does not exist in geometry.",
165
          uid)};
18✔
166
      }
167
    }
168
  }
169

170
  // Shared data that is collected from all threads
171
  const int n_domains = domain_ids_.size();
317✔
172

173
  // Divide work over MPI processes
174
  uint64_t min_samples = n_samples_ / mpi::n_procs;
317✔
175
  uint64_t remainder = n_samples_ % mpi::n_procs;
317✔
176
  uint64_t i_start, i_end;
177
  if (mpi::rank < remainder) {
317!
178
    i_start = (min_samples + 1) * mpi::rank;
×
179
    i_end = i_start + min_samples + 1;
×
180
  } else {
181
    i_start =
317✔
182
      (min_samples + 1) * remainder + (mpi::rank - remainder) * min_samples;
317✔
183
    i_end = i_start + min_samples;
317✔
184
  }
185

186
  while (true) {
187

188
#pragma omp parallel
5,867✔
189
    {
190
      // Temporary variables that are private to each thread
191
      CalcResults results(*this);
4,562✔
192
      results.sampling_time.start();
4,562✔
193

194
      // Sample locations and count scores
195
#pragma omp for
196
      for (size_t i = i_start; i < i_end; i++) {
3,359,612✔
197
        uint64_t id = master_results.iterations * n_samples_ + i;
3,355,050✔
198
        uint64_t seed = init_seed(id, STREAM_VOLUME);
3,355,050✔
199

200
        Position r {uniform_distribution(lower_left_.x, upper_right_.x, &seed),
3,355,050✔
201
          uniform_distribution(lower_left_.y, upper_right_.y, &seed),
3,355,050✔
202
          uniform_distribution(lower_left_.z, upper_right_.z, &seed)};
3,355,050✔
203

204
        results.n_samples++;
3,355,050✔
205

206
        switch (mode_) {
3,355,050!
207
        case EstMode::REJECTION: {
3,044,050✔
208
          constexpr double flt3 = 1. / std::sqrt(3.);
3,044,050✔
209
          Direction u {flt3, flt3, flt3};
3,044,050✔
210

211
          // Create zero-length ray, it is a bit excessive due to undemanded
212
          // internal ray variables initialization
213
          VolEstRay ray_hit(r, u, 0., 1., *this, results);
3,044,050✔
214
          ray_hit.score_hit();
3,044,050✔
215
        } break;
3,044,050✔
216
        case EstMode::RAYTRACE: {
311,000✔
217
          Direction u = isotropic_direction(&seed);
311,000✔
218

219
          // Compute lengths of the bounding box's chord segments
220
          const auto chord_len = get_box_chord(r, u);
311,000✔
221
          const double ch_len_tot = -chord_len.first + chord_len.second;
311,000✔
222

223
          r += chord_len.first * u;
311,000✔
224
          VolEstRay ray(r, u, ch_len_tot, 1. / ch_len_tot, *this, results);
311,000✔
225
          ray.trace(); // Trace from a boundary to another
311,000✔
226
        }
311,000✔
227
        }
228

229
        // This passing across all tallies after each sample can be
230
        // inefficient for the case of large number of domains and small
231
        // size of batch, but the batch size is currently assigned to 1
232
        // for keep the input format being unchanged
233
        const double batch_size_1 = 1. / static_cast<double>(1);
3,355,050✔
234
        for (auto& vt : results.vol_tallies) {
10,335,200✔
235
          for (auto& vol_tally : vt) {
21,880,615✔
236
            vol_tally.finalize_batch(batch_size_1);
14,900,465✔
237
          }
238
        }
239
      } // sample/batch loop
240

241
      results.sampling_time.stop();
4,562✔
242
      results.cost = results.sampling_time.elapsed();
4,562✔
243

244
      // At this point, each thread has its own volume tallies lists and we
245
      // now need to reduce them. OpenMP is not nearly smart enough to do this
246
      // on its own, so we have to manually reduce them
247
      reduce_results(results, master_results);
4,562✔
248
    } // omp parallel
4,562✔
249

250
    // bump iteration counter
251
    master_results.iterations++;
10,429✔
252

253
#ifdef OPENMC_MPI
254
    master_results.collect_MPI(); // collect results to master process
6,505✔
255
#endif
256
    // Process volume estimation results in master process for the trigger state
257
    // determination
258
    bool stop_calc =
10,429✔
259
      mpi::master && (trigger_type_ == TriggerMetric::not_active ||
17,469!
260
                       master_results.iterations == max_iterations_);
7,040✔
261

262
    if (!stop_calc) {
10,429✔
263
      // Compute current trigger state across totals (0th elements) only
264
      for (const auto& vt : master_results.vol_tallies) {
10,484✔
265
        stop_calc = vt[0].trigger_state(
10,388✔
266
          trigger_type_, threshold_cnd_, master_results.n_samples);
10,388✔
267
        if (!stop_calc)
10,388✔
268
          break;
10,152✔
269
      }
270
    }
271

272
#ifdef OPENMC_MPI
273
    // Send the state of calculation continuation just obtained in master
274
    // process to all processes
275
    MPI_Bcast(&stop_calc, 1, MPI_CXX_BOOL, 0, mpi::intracomm);
6,505✔
276
#endif
277

278
    if (!stop_calc)
10,429✔
279
      continue; // while loop
10,112✔
280
    // No trigger is applied or the trigger condition is satisfied, we're
281
    // done
282

283
    if (mpi::master) {
317✔
284
      // Normalize all results on the bounding primitive volume and compute
285
      // stddev
286
      for (auto& vt : master_results.vol_tallies) {
725✔
287
        for (auto& vol_tally : vt) {
1,488✔
288
          vol_tally.score_acc = get_tally_results(
1,010✔
289
            master_results.n_samples, volume_sample_, vol_tally);
1,010✔
290
        }
291
      }
292

293
      // Compute nuclides
294
      for (int i_domain = 0; i_domain < n_domains; ++i_domain) {
725✔
295
        // Create 2D array to store atoms/uncertainty for each nuclide. Later,
296
        // this is compressed into vectors storing only those nuclides that are
297
        // non-zero
298
        auto n_nuc =
299
          settings::run_CE ? data::nuclides.size() : data::mg.nuclides_.size();
478✔
300
        xt::xtensor<double, 2> atoms({n_nuc, 2}, 0.0);
478✔
301

302
        for (int j = 0; j < master_results.vol_tallies[i_domain].size(); ++j) {
1,488✔
303
          const int i_material = master_results.vol_tallies[i_domain][j].index;
1,010✔
304
          if (i_material == MATERIAL_VOID || i_material == _INDEX_TOTAL)
1,010✔
305
            continue;
489✔
306

307
          const auto& mat = model::materials[i_material];
521✔
308
          for (int k = 0; k < mat->nuclide_.size(); ++k) {
2,040✔
309
            auto& volume = master_results.vol_tallies[i_domain][j].score_acc;
1,519✔
310
            // Collect calculated nuclide amounts N [atoms] and stddev as
311
            // N = V [cm^3] * \ro [atoms/b-cm] * 1.e24 [b-cm/cm^3]
312
            atoms(mat->nuclide_[k], 0) +=
1,519✔
313
              volume[0] * mat->atom_density(k) * 1.0e24;
1,519✔
314
            atoms(mat->nuclide_[k], 1) +=
1,519✔
315
              volume[1] * mat->atom_density(k) * 1.0e24;
1,519✔
316
          }
317
        }
318

319
        // Get reference to result for this domain
320
        auto& result {master_results.nuc_results[i_domain]};
478✔
321
        // Convert full arrays to vectors
322
        for (int j = 0; j < n_nuc; ++j) {
3,440✔
323
          if (atoms(j, 0) > 0.0) {
2,962✔
324
            result.nuclides.push_back(j);
1,399✔
325
            result.atoms.push_back(atoms(j, 0));
1,399✔
326
            result.uncertainty.push_back(atoms(j, 1));
1,399✔
327
          }
328
        }
329
      } // end domains loop
478✔
330
    }
331

332
    return;
317✔
333

334
  } // end while
10,112✔
335
}
336

337
void VolumeCalculation::show_vol_stat(
873✔
338
  const std::string label, const std::string units, const double value) const
339
{
340
  fmt::print("{0:<20} = {2:10.4e} {1:<}\n", label, units, value);
873✔
341
}
873✔
342

343
void VolumeCalculation::show_volume(const std::string domain_type,
478✔
344
  const int domain_id, const std::string region_name, const double mean,
345
  const double stddev) const
346
{
347
  fmt::print(" {0:>9}{1:>6}: {2:10.4e} +/- {3:10.4e} cm^3", domain_type,
394✔
348
    domain_id, mean, stddev);
349
  if (!region_name.empty()) {
478✔
350
    fmt::print("   //{:<}", region_name);
80✔
351
  }
352
  fmt::print("\n");
478✔
353
}
478✔
354

355
void VolumeCalculation::show_results(const CalcResults& results) const
247✔
356
{
357
  // Show tracing statistics
358
  write_message(5, " ");
247✔
359
  show_vol_stat(
247✔
360
    "Total sample size", "", static_cast<double>(results.n_samples));
247✔
361
  show_vol_stat("Running cost", "thread-sec", results.cost);
247✔
362

363
  switch (mode_) {
247!
364
  case EstMode::REJECTION:
181✔
365
    show_vol_stat("Cost of hitting", "thread-sec/hit",
181✔
366
      static_cast<double>(results.cost) /
181✔
367
        static_cast<double>(results.n_samples));
181✔
368
    break;
181✔
369
  case EstMode::RAYTRACE:
66✔
370
    show_vol_stat("Rays traced", "", results.n_rays);
66✔
371
    show_vol_stat("Average ray length", "segments/ray",
66✔
372
      static_cast<double>(results.n_segs) /
66✔
373
        static_cast<double>(results.n_rays));
66✔
374
    show_vol_stat("Cost of tracing", "thread-sec/segment",
66✔
375
      static_cast<double>(results.cost) / static_cast<double>(results.n_segs));
66✔
376
    if (results.n_errors != 0)
66!
NEW
377
      show_vol_stat("Error rate", "errors/ray",
×
NEW
378
        static_cast<double>(results.n_errors) /
×
NEW
379
          static_cast<double>(results.n_rays));
×
380
  }
381
  write_message(5, " ");
247✔
382

383
  std::string domain_type;
247✔
384
  if (domain_type_ == TallyDomain::CELL) {
247✔
385
    domain_type = "Cell";
112✔
386
  } else if (domain_type_ == TallyDomain::MATERIAL) {
135✔
387
    domain_type = "Material";
102✔
388
  } else {
389
    domain_type = "Universe";
33✔
390
  }
391

392
  // Display domain volumes
393
  for (int j = 0; j < domain_ids_.size(); j++) {
725✔
394
    std::string region_name {""};
478✔
395
    if (domain_type_ == TallyDomain::CELL) {
478✔
396
      int cell_idx = model::cell_map[domain_ids_[j]];
266✔
397
      region_name = model::cells[cell_idx]->name();
266✔
398
    } else if (domain_type_ == TallyDomain::MATERIAL) {
212✔
399
      int mat_idx = model::material_map[domain_ids_[j]];
179✔
400
      region_name = model::materials[mat_idx]->name();
179✔
401
    }
402
    if (region_name.size())
478✔
403
      region_name.insert(0, " "); // prepend space for formatting
80✔
404

405
    show_volume(domain_type, domain_ids_[j], region_name,
478✔
406
      results.vol_tallies[j][0].score_acc[0],
478✔
407
      results.vol_tallies[j][0].score_acc[1]);
478✔
408
  }
478✔
409
  write_message(4, " "); // Blank line afer results printed
247✔
410
}
247✔
411

412
void VolumeCalculation::to_hdf5(
247✔
413
  const std::string& filename, const CalcResults& results) const
414
{
415
  // Create HDF5 file
416
  hid_t file_id = file_open(filename, 'w');
247✔
417

418
  // Write header info
419
  write_attribute(file_id, "filetype", "volume");
247✔
420
  write_attribute(file_id, "version", VERSION_VOLUME);
247✔
421
  write_attribute(file_id, "openmc_version", VERSION);
247✔
422
#ifdef GIT_SHA1
423
  write_attribute(file_id, "git_sha1", GIT_SHA1);
424
#endif
425

426
  // Write current date and time
427
  write_attribute(file_id, "date_and_time", time_stamp());
247✔
428

429
  // Write basic metadata
430
  write_attribute(file_id, "samples", n_samples_);
247✔
431
  write_attribute(file_id, "lower_left", lower_left_);
247✔
432
  write_attribute(file_id, "upper_right", upper_right_);
247✔
433
  // Write estimator info
434
  std::string estimator_str;
247✔
435
  switch (mode_) {
247!
436
  case EstMode::REJECTION:
181✔
437
    estimator_str = "hit";
181✔
438
    break;
181✔
439
  case EstMode::RAYTRACE:
66✔
440
    estimator_str = "ray";
66✔
441
    break;
66✔
442
  }
443
  write_attribute(file_id, "estimator_type", estimator_str);
247✔
444
  // Write trigger info
445
  if (trigger_type_ != TriggerMetric::not_active) {
247✔
446
    write_attribute(file_id, "iterations", results.iterations);
88✔
447
    write_attribute(file_id, "threshold", threshold_);
88✔
448
    std::string trigger_str;
88✔
449
    switch (trigger_type_) {
88!
450
    case TriggerMetric::variance:
22✔
451
      trigger_str = "variance";
22✔
452
      break;
22✔
453
    case TriggerMetric::standard_deviation:
22✔
454
      trigger_str = "std_dev";
22✔
455
      break;
22✔
456
    case TriggerMetric::relative_error:
44✔
457
      trigger_str = "rel_err";
44✔
458
      break;
44✔
459
    default:
×
460
      break;
×
461
    }
462
    write_attribute(file_id, "trigger_type", trigger_str);
88✔
463
    // check max_iterations on default value
464
    if (max_iterations_ !=
88✔
465
        std::numeric_limits<decltype(max_iterations_)>::max())
88✔
466
      write_attribute(file_id, "max_iterations", max_iterations_);
22✔
467
  } else {
88✔
468
    write_attribute(file_id, "iterations", 1);
159✔
469
  }
470

471
  if (domain_type_ == TallyDomain::CELL) {
247✔
472
    write_attribute(file_id, "domain_type", "cell");
112✔
473
  } else if (domain_type_ == TallyDomain::MATERIAL) {
135✔
474
    write_attribute(file_id, "domain_type", "material");
102✔
475
  } else if (domain_type_ == TallyDomain::UNIVERSE) {
33!
476
    write_attribute(file_id, "domain_type", "universe");
33✔
477
  }
478

479
  for (int i = 0; i < domain_ids_.size(); ++i) {
725✔
480
    hid_t group_id =
481
      create_group(file_id, fmt::format("domain_{}", domain_ids_[i]));
956✔
482

483
    // Write volume for domain
484
    const auto& result {results.nuc_results[i]};
478✔
485
    write_dataset(group_id, "volume", results.vol_tallies[i][0].score_acc);
478✔
486

487
    // Create array of nuclide names from the vector
488
    auto n_nuc = result.nuclides.size();
478✔
489

490
    vector<std::string> nucnames;
478✔
491
    for (int i_nuc : result.nuclides) {
1,877✔
492
      nucnames.push_back(settings::run_CE ? data::nuclides[i_nuc]->name_
1,883✔
493
                                          : data::mg.nuclides_[i_nuc].name);
484✔
494
    }
495

496
    // Create array of total # of atoms with uncertainty for each nuclide
497
    xt::xtensor<double, 2> atom_data({n_nuc, 2});
478✔
498
    xt::view(atom_data, xt::all(), 0) = xt::adapt(result.atoms);
478✔
499
    xt::view(atom_data, xt::all(), 1) = xt::adapt(result.uncertainty);
478✔
500

501
    // Write results
502
    write_dataset(group_id, "nuclides", nucnames);
478✔
503
    write_dataset(group_id, "atoms", atom_data);
478✔
504

505
    close_group(group_id);
478✔
506
  }
478✔
507

508
  file_close(file_id);
247✔
509
}
247✔
510

511
void VolumeCalculation::check_hit(const int32_t i_material,
10,687,121✔
512
  const double contrib, vector<VolTally>& vol_tallies) const
513
{
514
  // Contribute to entire domain result tally
515
  vol_tallies[0].score += contrib;
10,687,121✔
516

517
  // Check if this material was previously hit and if so, contribute score
518
  for (int j = 1; j < vol_tallies.size(); j++) {
11,335,821✔
519
    if (vol_tallies[j].index == i_material) {
11,288,587✔
520
      vol_tallies[j].score += contrib;
10,639,887✔
521
      return;
10,639,887✔
522
    }
523
  }
524

525
  // The material was not previously hit, append an entry to the material
526
  // indices and hits lists
527
  vol_tallies.push_back(VolTally(i_material, contrib));
47,234✔
528
}
529

530
void VolumeCalculation::reduce_results(
16,296✔
531
  const CalcResults& local_results, CalcResults& results) const
532
{
533
  auto n_threads = num_threads();
16,296✔
534

535
  // Collect scalar variables
536
#pragma omp for ordered schedule(static, 1)
11,734!
537
  for (int i = 0; i < n_threads; ++i) {
9,124✔
538
#pragma omp ordered
539
    {
540
      results.n_samples += local_results.n_samples;
16,296✔
541
      results.n_rays += local_results.n_rays;
16,296✔
542
      results.n_segs += local_results.n_segs;
16,296✔
543
      results.n_errors += local_results.n_errors;
16,296✔
544
      results.cost += local_results.cost;
16,296!
545
    }
546
  }
547

548
  // Collect vectored domain-wise results
549
  for (int i_domain = 0; i_domain < domain_ids_.size(); ++i_domain) {
64,048✔
550
    const vector<VolTally>& local_vol_tall =
551
      local_results.vol_tallies[i_domain];
47,752✔
552
    vector<VolTally>& vol_tall = results.vol_tallies[i_domain];
47,752✔
553

554
#pragma omp for ordered schedule(static, 1)
34,388!
555
    for (int i = 0; i < n_threads; ++i) {
26,728✔
556
#pragma omp ordered
557
      for (int j = 0; j < local_vol_tall.size(); ++j) {
142,738✔
558
        // Check if this material has been added to the master list and if
559
        // so, accumulate scores
560
        const auto ind {local_vol_tall[j].index};
94,986✔
561
        const auto it = std::find_if(vol_tall.begin(), vol_tall.end(),
94,986✔
562
          [ind](const VolTally& vt) { return vt.index == ind; });
3,340,387✔
563
        if (it == vol_tall.end()) {
94,986✔
564
          vol_tall.push_back(local_vol_tall[j]);
13,172✔
565
        } else {
566
          vol_tall[it - vol_tall.begin()].append_tally(local_vol_tall[j]);
81,814✔
567
        }
568
      }
569
    }
570
  }
571
}
16,296✔
572

573
#ifdef OPENMC_MPI
574
void VolumeCalculation::initialize_MPI_struct() const
197✔
575
{
576
  // This code is a slightly modified replica of initialize_mpi() from
577
  // initialize.cpp. It works under GCC in the Release configuration, but not
578
  // sure that the adress offsets of structure's memebrs should be necessary the
579
  // same everywhere as using an optimizing compiler.
580
  CalcResults cr(*this);
197✔
581
  MPI_Aint cr_disp[5], cr_d;
582
  MPI_Get_address(&cr, &cr_d);
197✔
583
  MPI_Get_address(&cr.n_errors, &cr_disp[0]);
197✔
584
  MPI_Get_address(&cr.n_rays, &cr_disp[1]);
197✔
585
  MPI_Get_address(&cr.n_segs, &cr_disp[2]);
197✔
586
  MPI_Get_address(&cr.n_samples, &cr_disp[3]);
197✔
587
  MPI_Get_address(&cr.cost, &cr_disp[4]);
197✔
588
  for (int i = 0; i < 5; i++) {
1,182✔
589
    cr_disp[i] -= cr_d;
985✔
590
  }
591

592
  int cr_blocks[] {1, 1, 1, 1, 1};
197✔
593
  MPI_Datatype cr_types[] {
197✔
594
    MPI_UINT64_T, MPI_UINT64_T, MPI_UINT64_T, MPI_UINT64_T, MPI_DOUBLE};
595
  MPI_Type_create_struct(5, cr_blocks, cr_disp, cr_types, &mpi_vol_results);
197✔
596
  MPI_Type_commit(&mpi_vol_results);
197✔
597

598
  VolTally vt;
197✔
599
  MPI_Aint vt_disp[3], vt_d;
600
  MPI_Get_address(&vt, &vt_d);
197✔
601
  MPI_Get_address(&vt.score, &vt_disp[0]);
197✔
602
  MPI_Get_address(&vt.score_acc, &vt_disp[1]);
197✔
603
  MPI_Get_address(&vt.index, &vt_disp[2]);
197✔
604
  for (int i = 0; i < 3; i++) {
788✔
605
    vt_disp[i] -= vt_d;
591✔
606
  }
607

608
  int vt_blocks[] {1, 2, 1};
197✔
609
  MPI_Datatype vt_types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT};
197✔
610
  MPI_Type_create_struct(3, vt_blocks, vt_disp, vt_types, &mpi_volume_tally);
197✔
611
  MPI_Type_commit(&mpi_volume_tally);
197✔
612
}
197✔
613

614
void VolumeCalculation::delete_MPI_struct() const
185✔
615
{
616
  MPI_Type_free(&mpi_vol_results);
185✔
617
  MPI_Type_free(&mpi_volume_tally);
185✔
618
}
185✔
619

620
void VolumeCalculation::CalcResults::collect_MPI()
6,505✔
621
{
622
  vector<int> domain_sizes(vol_tallies.size());
6,505✔
623

624
  if (!mpi::master) {
6,505✔
625
    // n_domain + 2 MPI messages will be send in total to node mpi::master
626

627
    // To determination of an unique tag for each MPI message (as below)
628
    int mpi_offset = mpi::rank * (vol_tallies.size() + 2) + 2;
3,230✔
629

630
    // Pass root data of the struct
631
    MPI_Send(
3,230✔
632
      (void*)this, 1, mpi_vol_results, 0, mpi_offset - 2, mpi::intracomm);
633

634
    // Pass sizes of domain-wise data
635
    for (int i_domain = 0; i_domain < vol_tallies.size(); i_domain++) {
12,740✔
636
      domain_sizes[i_domain] = vol_tallies[i_domain].size();
9,510✔
637
    }
638

639
    MPI_Send(domain_sizes.data(), domain_sizes.size(), MPI_INT, 0,
3,230✔
640
      mpi_offset - 1, mpi::intracomm);
641

642
    // Pass domain-wise data of struct
643
    for (int i_domain = 0; i_domain < vol_tallies.size(); i_domain++) {
12,740✔
644
      MPI_Send(vol_tallies[i_domain].data(), domain_sizes[i_domain],
9,510✔
645
        mpi_volume_tally, 0, mpi_offset + i_domain, mpi::intracomm);
646
    }
647

648
    this->reset(); // Delete passed to main process data
3,230✔
649

650
  } else {
651

652
    // n_domain + 2 MPI messages will be recieved in total on node mpi::master
653

654
    for (int i_proc = 1; i_proc < mpi::n_procs; i_proc++) {
6,505✔
655

656
      CalcResults res_buff(*this); // temporary storage for recived data
3,230✔
657

658
      // To determination of an unique tag for each MPI message (as above)
659
      int mpi_offset = i_proc * (vol_tallies.size() + 2) + 2;
3,230✔
660

661
      // Pass root data of struct
662
      MPI_Recv(&res_buff, 1, mpi_vol_results, i_proc, mpi_offset - 2,
3,230✔
663
        mpi::intracomm, MPI_STATUS_IGNORE);
664

665
      // Pass sizes of domain-wise data
666
      MPI_Recv(domain_sizes.data(), domain_sizes.size(), MPI_INT, i_proc,
3,230✔
667
        mpi_offset - 1, mpi::intracomm, MPI_STATUS_IGNORE);
668

669
      // Pass domain-wise data of struct
670
      for (int i_domain = 0; i_domain < vol_tallies.size(); i_domain++) {
12,740✔
671
        res_buff.vol_tallies[i_domain].resize(domain_sizes[i_domain]);
9,510✔
672
        MPI_Recv(res_buff.vol_tallies[i_domain].data(), domain_sizes[i_domain],
9,510✔
673
          mpi_volume_tally, i_proc, mpi_offset + i_domain, mpi::intracomm,
674
          MPI_STATUS_IGNORE);
675
      }
676

677
      this->append(res_buff);
3,230✔
678
    }
3,230✔
679
  }
680
}
6,505✔
681
#endif
682

683
VolumeCalculation::CalcResults::CalcResults(const VolumeCalculation& vol_calc)
16,837✔
684
{
685
  n_samples = 0;
16,837✔
686
  n_rays = 0;
16,837✔
687
  n_segs = 0;
16,837✔
688
  n_errors = 0;
16,837✔
689
  iterations = 0;
16,837✔
690
  cost = 0.;
16,837✔
691
  for (int i = 0; i < vol_calc.domain_ids_.size(); i++) {
65,652✔
692
    vector<VolTally> vt_vect; // Tally group for a domain
48,815✔
693
    vt_vect.push_back(
48,815✔
694
      VolTally(_INDEX_TOTAL)); // Zero-element tally for entire domain totals
48,815✔
695
    vol_tallies.push_back(vt_vect);
48,815✔
696
    nuc_results.push_back(NuclResult());
48,815✔
697
  }
48,815✔
698
}
16,837✔
699

700
void VolumeCalculation::CalcResults::reset()
3,230✔
701
{
702
  n_samples = 0;
3,230✔
703
  n_rays = 0;
3,230✔
704
  n_segs = 0;
3,230✔
705
  n_errors = 0;
3,230✔
706
  cost = 0.;
3,230✔
707
  for (auto& vt : vol_tallies) {
12,740✔
708
    std::fill(vt.begin(), vt.end(), VolTally());
9,510!
709
  }
710
}
3,230✔
711

712
void VolumeCalculation::CalcResults::append(const CalcResults& other)
3,230✔
713
{
714
  n_samples += other.n_samples;
3,230✔
715
  n_rays += other.n_rays;
3,230✔
716
  n_segs += other.n_segs;
3,230✔
717
  n_errors += other.n_errors;
3,230✔
718
  cost += other.cost;
3,230✔
719

720
  // The domain-wise vectors this.vol_tallies and other.vol_tallies are
721
  // conformed each to other by definition
722
  for (auto id = 0; id < vol_tallies.size(); id++) {
12,740✔
723
    // Merging current domain vector from other.vol_tallies into this via
724
    // pair-wise comparisons
725
    for (const auto& vt_other : other.vol_tallies[id]) {
1,366,310✔
726
      bool already_appended = false;
1,356,800✔
727
      for (auto& vt : vol_tallies[id]) {
3,373,240✔
728
        if (vt.index == vt_other.index) {
3,373,200✔
729
          vt.append_tally(vt_other);
1,356,760✔
730
          already_appended = true;
1,356,760✔
731
          break;
1,356,760✔
732
        }
733
      }
734
      if (!already_appended)
1,356,800✔
735
        vol_tallies[id].push_back(vt_other);
40!
736
    }
737
  }
738
}
3,230✔
739

740
inline VolumeCalculation::VolTally::VolTally(const int i_material,
1,440,406✔
741
  const double contrib, const double score_acc_, const double score2_acc_)
1,440,406✔
742
{
743
  score = contrib;
1,440,406✔
744
  score_acc[0] = score_acc_;
1,440,406✔
745
  score_acc[1] = score2_acc_;
1,440,406✔
746
  index = i_material;
1,440,406✔
747
}
1,440,406✔
748

749
inline void VolumeCalculation::VolTally::finalize_batch(
56,700,445✔
750
  const double batch_size_1)
751
{
752
  if (score != 0.) {
56,700,445✔
753
    score *= batch_size_1;
21,374,198✔
754
    score_acc[0] += score;
21,374,198✔
755
    score_acc[1] += score * score;
21,374,198✔
756
    score = 0.;
21,374,198✔
757
  }
758
}
56,700,445✔
759

760
inline void VolumeCalculation::VolTally::assign_tally(const VolTally& vol_tally)
761
{
762
  score = vol_tally.score;
763
  score_acc = vol_tally.score_acc;
764
  index = vol_tally.index;
765
}
766

767
inline void VolumeCalculation::VolTally::append_tally(const VolTally& vol_tally)
1,438,574✔
768
{
769
  score += vol_tally.score;
1,438,574✔
770
  score_acc[0] += vol_tally.score_acc[0];
1,438,574✔
771
  score_acc[1] += vol_tally.score_acc[1];
1,438,574✔
772
}
1,438,574✔
773

774
inline bool VolumeCalculation::VolTally::trigger_state(
10,388✔
775
  const TriggerMetric trigger_type, const double threshold,
776
  const size_t& n_samples) const
777
{
778
  // For sample contribution to volume fraction limited by 1, the maximal
779
  // allowed n_samples value is around 1.e102, but this is still much larger
780
  // than the size_t limit equal to ~1.8e19
781
  const double ns1 = static_cast<double>(n_samples - 1);
10,388✔
782
  const double ns = static_cast<double>(n_samples);
10,388✔
783
  const double mean_xi_sq = score_acc[0] * score_acc[0];
10,388✔
784

785
  // Adjusted sample variance: s^2 = (\bar{x^2} - \bar{x}^2) / (N-1)
786
  // \bar{x}^2 = mean_xi_sq / N^2, \bar{\x^2} = score_acc[1] / N, N = ns
787
  switch (trigger_type) {
10,388✔
788
  case TriggerMetric::variance:
9,848✔
789
  case TriggerMetric::standard_deviation:
790
    // Condition: s^2 / N < t'
791
    // Equivalent implementation:
792
    // N^2 * (\bar{x^2} - \bar{x}^2) < t' * (N-1) * N^2
793
    return score_acc[1] * ns - mean_xi_sq < threshold * ns1 * ns * ns;
9,848✔
794
  case TriggerMetric::relative_error:
480✔
795
    // Condition: (s^2 / \mu^2) / N < t'
796
    // Equivalent implementation:
797
    // N^2 * (\bar{x^2} - \bar{x}^2) < t' * (N-1) * (N * \bar{x})^2
798
    return score_acc[1] * ns - mean_xi_sq < threshold * ns1 * mean_xi_sq;
480✔
799
  default:
60✔
800
    return true;
60✔
801
  }
802
}
803

804
array<double, 2> VolumeCalculation::get_tally_results(const size_t& n_samples,
1,010✔
805
  const double coeff_norm, const VolTally& vol_tally) const
806
{
807
  array<double, 2> volume;
808
  const double ns_1 = 1. / static_cast<double>(n_samples);
1,010✔
809
  volume[0] = vol_tally.score_acc[0] * ns_1;
1,010✔
810
  volume[1] = vol_tally.score_acc[1] * ns_1;
1,010✔
811
  volume[1] = std::sqrt(
2,020✔
812
    (volume[1] - volume[0] * volume[0]) / static_cast<double>(n_samples - 1));
1,010✔
813
  volume[0] *= coeff_norm;
1,010✔
814
  volume[1] *= coeff_norm;
1,010✔
815
  return volume;
1,010✔
816
}
817

818
std::pair<double, double> VolumeCalculation::get_box_chord(
684,200✔
819
  const Position& r, const Direction& u) const
820
{
821
  // Compute distanses to each box plane orthogonal to an axis
822
  Direction u_1 = {1., 1., 1.};
684,200✔
823
  u_1 = u_1 / u;
684,200✔
824
  Position xi = (lower_left_ - r) * u_1;
684,200✔
825
  const array<double, 3> dist1 = {xi.x, xi.y, xi.z};
684,200✔
826
  xi = (upper_right_ - r) * u_1;
684,200✔
827
  const array<double, 3> dist2 = {xi.x, xi.y, xi.z};
684,200✔
828

829
  // Find the minimal forward (positive values) and backward (negative values)
830
  // distances across the computed ones (probably there is some STL
831
  // alternatives)
832
  std::pair<double, double> chord_lengths {std::minmax(dist1[0], dist2[0])};
684,200✔
833
  for (int i = 1; i < 3; i++) {
2,052,600✔
834
    const std::pair<double, double> dist_mm = std::minmax(dist1[i], dist2[i]);
1,368,400✔
835
    chord_lengths.first = std::max(chord_lengths.first, dist_mm.first);
1,368,400✔
836
    chord_lengths.second = std::min(chord_lengths.second, dist_mm.second);
1,368,400✔
837
  }
838
  return chord_lengths;
684,200✔
839
}
840

841
void VolEstRay::on_intersection()
1,371,656✔
842
{
843
  if (traversal_distance_ == 0.) {
1,371,656✔
844
    return; // No crossing model
620,972✔
845
  }
846

847
  results_.n_segs++;
750,684✔
848
  if (traversal_distance_last_ == 0.) {
750,684✔
849
    results_.n_rays++; // First segment of new ray
620,972✔
850
  }
851

852
  // At this point, current GeometryState parameters represent the cell behind
853
  // crossed surface, but the segment length is known for the previous
854
  // cell only, therefore we use below the last cell keept in GeometryState
855
  const auto score = (std::min(traversal_distance_, traversal_distance_max_) -
750,684✔
856
                       traversal_distance_last_) *
750,684✔
857
                     coeff_mult_;
750,684✔
858

859
  //----------------------------------------------------------------------------
860
  // Tracing error diagnostic
861
  // TODO: those can be implemented here clever diagnostic and guidance for
862
  // user to fix input mistakes
863
  //----------------------------------------------------------------------------
864

865
  if (traversal_distance_ >= traversal_distance_max_) {
750,684!
NEW
866
    stop();
×
867
  } else {
868
    traversal_distance_last_ = traversal_distance_;
750,684✔
869
  }
870

871
  // In a case of single-segment ray (leakage after 1st surface crossing),
872
  // current segment material ID seems to be contained in material(), but in
873
  // other cases it is in material_last(). Due to this unclear behavior, it is
874
  // used below a more fundamental way for material determination -- via the
875
  // stable last cell record.
876
  vol_scoring(VolumeCalculation::EstMode::RAYTRACE, score,
1,501,368✔
877
    (current_cell(VolumeCalculation::EstMode::RAYTRACE, n_coord_last() - 1)
750,684✔
878
        ->material(n_coord_last() - 1)));
750,684✔
879
}
880

881
void VolEstRay::score_hit()
11,696,910✔
882
{
883
  if (exhaustive_find_cell(*this))
11,696,910✔
884
    vol_scoring(VolumeCalculation::EstMode::REJECTION, 1.,
9,975,003✔
885
      material()); // One hit score
9,975,003✔
886
}
11,696,910✔
887

888
void VolEstRay::vol_scoring(
10,725,687✔
889
  const VolumeCalculation::EstMode mode, const double score, const int id_mat)
890
{
891
  const auto n_domains = vol_calc_.domain_ids_.size();
10,725,687✔
892

893
  switch (vol_calc_.domain_type_) {
10,725,687!
894
  case VolumeCalculation::TallyDomain::MATERIAL:
4,689,875✔
895
    if (id_mat != MATERIAL_VOID) {
4,689,875✔
896
      for (auto i_domain = 0; i_domain < n_domains; i_domain++) {
7,652,901✔
897
        if (model::materials[id_mat]->id_ == vol_calc_.domain_ids_[i_domain]) {
7,643,749✔
898
          vol_calc_.check_hit(id_mat, score, results_.vol_tallies[i_domain]);
4,675,740✔
899
          break;
4,675,740✔
900
        }
901
      }
902
    }
903
    break;
4,689,875✔
904
  case VolumeCalculation::TallyDomain::CELL:
4,396,977✔
905
    for (auto level = 0; level < n_coord_last(); ++level) {
8,793,954✔
906
      for (auto i_domain = 0; i_domain < n_domains; i_domain++) {
4,964,291✔
907
        if (current_cell(mode, level)->id_ == vol_calc_.domain_ids_[i_domain]) {
4,939,860✔
908
          vol_calc_.check_hit(id_mat, score, results_.vol_tallies[i_domain]);
4,372,546✔
909
          break;
4,372,546✔
910
        }
911
      }
912
    }
913
    break;
4,396,977✔
914
  case VolumeCalculation::TallyDomain::UNIVERSE:
1,638,835✔
915
    for (auto level = 0; level < n_coord_last(); ++level) {
3,277,670✔
916
      for (auto i_domain = 0; i_domain < n_domains; ++i_domain) {
1,638,835!
917
        if (model::universes[current_cell(mode, level)->universe_]->id_ ==
3,277,670!
918
            vol_calc_.domain_ids_[i_domain]) {
1,638,835✔
919
          vol_calc_.check_hit(id_mat, score, results_.vol_tallies[i_domain]);
1,638,835✔
920
          break;
1,638,835✔
921
        }
922
      }
923
    }
924
  }
925
}
10,725,687✔
926

927
inline std::unique_ptr<Cell>& VolEstRay::current_cell(
7,329,379✔
928
  VolumeCalculation::EstMode mode, int level)
929
{
930
  switch (mode) {
7,329,379!
931
  case VolumeCalculation::EstMode::REJECTION:
5,594,635✔
932
    return model::cells[coord(level).cell()]; // Current position
5,594,635✔
933
  case VolumeCalculation::EstMode::RAYTRACE:
1,734,744✔
934
    return model::cells[cell_last(level)]; // Previous segment
1,734,744✔
935
  }
UNCOV
936
}
×
937

938
void free_memory_volume()
7,490✔
939
{
940
  openmc::model::volume_calcs.clear();
7,490✔
941
}
7,490✔
942

943
} // namespace openmc
944

945
//==============================================================================
946
// OPENMC_CALCULATE_VOLUMES runs each of the stochastic volume calculations
947
// that the user has specified and writes results to HDF5 files
948
//==============================================================================
949

950
int openmc_calculate_volumes()
128✔
951
{
952
  using namespace openmc;
953

954
  if (mpi::master) {
128✔
955
    header("STOCHASTIC VOLUME CALCULATION", 3);
118✔
956
  }
957
  Timer time_volume;
128✔
958
  time_volume.start();
128✔
959

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

963
    // Run volume calculation
964
    const auto& vol_calc {model::volume_calcs[i]};
344✔
965
    VolumeCalculation::CalcResults results(vol_calc);
344✔
966

967
#ifdef OPENMC_MPI
968
    vol_calc.initialize_MPI_struct();
197✔
969
#endif
970

971
    try {
972
      vol_calc.execute(results);
344✔
973
    } catch (const std::exception& e) {
27!
974
      set_errmsg(e.what());
27✔
975
      return OPENMC_E_UNASSIGNED;
27✔
976
    }
27✔
977

978
#ifdef OPENMC_MPI
979
    vol_calc.delete_MPI_struct();
185✔
980
#endif
981

982
    if (mpi::master) {
317✔
983

984
      // Output volume calculation results and statistics
985
      vol_calc.show_results(results);
247✔
986

987
      // Write volumes to HDF5 file
988
      std::string filename =
989
        fmt::format("{}volume_{}.h5", settings::path_output, i + 1);
450✔
990
      vol_calc.to_hdf5(filename, results);
247✔
991
    }
247✔
992
  }
344✔
993

994
  // Show elapsed time
995
  time_volume.stop();
101✔
996
  write_message(6, "Elapsed time: {} sec", time_volume.elapsed());
101✔
997

998
  return 0;
101✔
999
}
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