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

openmc-dev / openmc / 30830380864

03 Aug 2026 04:03PM UTC coverage: 81.48% (+0.05%) from 81.429%
30830380864

Pull #3971

github

web-flow
Merge 13d0297e2 into 5982acdf8
Pull Request #3971: Delta tracking

18827 of 27241 branches covered (69.11%)

Branch coverage included in aggregate %.

597 of 645 new or added lines in 20 files covered. (92.56%)

60852 of 70549 relevant lines covered (86.25%)

50385682.17 hits per line

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

90.93
/src/eigenvalue.cpp
1
#include "openmc/eigenvalue.h"
2

3
#include "openmc/tensor.h"
4

5
#include "openmc/array.h"
6
#include "openmc/bank.h"
7
#include "openmc/capi.h"
8
#include "openmc/constants.h"
9
#include "openmc/error.h"
10
#include "openmc/hdf5_interface.h"
11
#include "openmc/ifp.h"
12
#include "openmc/math_functions.h"
13
#include "openmc/mesh.h"
14
#include "openmc/message_passing.h"
15
#include "openmc/random_lcg.h"
16
#include "openmc/search.h"
17
#include "openmc/settings.h"
18
#include "openmc/simulation.h"
19
#include "openmc/tallies/tally.h"
20
#include "openmc/timer.h"
21

22
#include <algorithm> // for min
23
#include <cmath>     // for sqrt, abs, pow
24
#include <iterator>  // for back_inserter
25
#include <limits>    //for infinity
26
#include <string>
27

28
namespace openmc {
29

30
//==============================================================================
31
// Global variables
32
//==============================================================================
33

34
namespace simulation {
35

36
double keff_generation;
37
array<double, 2> k_sum;
38
vector<double> entropy;
39
tensor::Tensor<double> source_frac;
40

41
} // namespace simulation
42

43
//==============================================================================
44
// Non-member functions
45
//==============================================================================
46

47
void calculate_generation_keff()
102,873✔
48
{
49
  const auto& gt = simulation::global_tallies;
102,873✔
50

51
  // Get keff for this generation by subtracting off the starting value
52
  if (settings::delta_tracking) {
102,873✔
53
    simulation::keff_generation =
1,500✔
54
      gt(GlobalTally::K_COLLISION, TallyResult::VALUE) -
1,500✔
55
      simulation::keff_generation;
56
  } else {
57
    simulation::keff_generation =
101,373✔
58
      gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) -
101,373✔
59
      simulation::keff_generation;
60
  }
61

62
  double keff_reduced;
102,873✔
63
#ifdef OPENMC_MPI
64
  if (settings::solver_type != SolverType::RANDOM_RAY) {
46,632✔
65
    // Combine values across all processors
66
    MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE,
43,352✔
67
      MPI_SUM, mpi::intracomm);
68
  } else {
69
    // If using random ray, MPI parallelism is provided by domain replication.
70
    // As such, all fluxes will be reduced at the end of each transport sweep,
71
    // such that all ranks have identical scalar flux vectors, and will all
72
    // independently compute the same value of k. Thus, there is no need to
73
    // perform any additional MPI reduction here.
74
    keff_reduced = simulation::keff_generation;
3,280✔
75
  }
76
#else
77
  keff_reduced = simulation::keff_generation;
56,241✔
78
#endif
79

80
  // Normalize single batch estimate of k
81
  // TODO: This should be normalized by total_weight, not by n_particles
82
  if (settings::solver_type != SolverType::RANDOM_RAY) {
102,873✔
83
    keff_reduced /= settings::n_particles;
96,023✔
84
  }
85

86
  simulation::k_generation.push_back(keff_reduced);
102,873✔
87
}
102,873✔
88

89
void synchronize_bank()
96,023✔
90
{
91
  simulation::time_bank.start();
96,023✔
92

93
  // In order to properly understand the fission bank algorithm, you need to
94
  // think of the fission and source bank as being one global array divided
95
  // over multiple processors. At the start, each processor has a random amount
96
  // of fission bank sites -- each processor needs to know the total number of
97
  // sites in order to figure out the probability for selecting
98
  // sites. Furthermore, each proc also needs to know where in the 'global'
99
  // fission bank its own sites starts in order to ensure reproducibility by
100
  // skipping ahead to the proper seed.
101

102
#ifdef OPENMC_MPI
103
  int64_t start = 0;
43,352✔
104
  int64_t n_bank = simulation::fission_bank.size();
43,352✔
105
  MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
43,352✔
106

107
  // While we would expect the value of start on rank 0 to be 0, the MPI
108
  // standard says that the receive buffer on rank 0 is undefined and not
109
  // significant
110
  if (mpi::rank == 0)
43,352✔
111
    start = 0;
30,077✔
112

113
  int64_t finish = start + simulation::fission_bank.size();
43,352✔
114
  int64_t total = finish;
43,352✔
115
  MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm);
43,352✔
116

117
#else
118
  int64_t start = 0;
52,671✔
119
  int64_t finish = simulation::fission_bank.size();
52,671!
120
  int64_t total = finish;
52,671✔
121
#endif
122

123
  // If there are not that many particles per generation, it's possible that no
124
  // fission sites were created at all on a single processor. Rather than add
125
  // extra logic to treat this circumstance, we really want to ensure the user
126
  // runs enough particles to avoid this in the first place.
127

128
  if (simulation::fission_bank.size() == 0) {
96,023!
129
    fatal_error(
×
130
      "No fission sites banked on MPI rank " + std::to_string(mpi::rank));
×
131
  }
132

133
  simulation::time_bank_sample.start();
96,023✔
134

135
  // Allocate temporary source bank -- we don't really know how many fission
136
  // sites were created, so overallocate by a factor of 3
137
  int64_t index_temp = 0;
96,023✔
138

139
  vector<SourceSite> temp_sites(3 * simulation::work_per_rank);
96,023✔
140

141
  // Temporary banks for IFP
142
  vector<vector<int>> temp_delayed_groups;
96,023✔
143
  vector<vector<double>> temp_lifetimes;
96,023✔
144
  if (settings::ifp_on) {
96,023✔
145
    resize_ifp_data(
1,480✔
146
      temp_delayed_groups, temp_lifetimes, 3 * simulation::work_per_rank);
147
  }
148

149
  // ==========================================================================
150
  // SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
151

152
  // We use Uniform Combing method to exactly get the targeted particle size
153
  // [https://doi.org/10.1080/00295639.2022.2091906]
154

155
  // Make sure all processors use the same random number seed.
156
  int64_t id = simulation::total_gen + overall_generation();
96,023✔
157
  uint64_t seed = init_seed(id, STREAM_TRACKING);
96,023✔
158

159
  // Comb specification
160
  double teeth_distance = static_cast<double>(total) / settings::n_particles;
96,023✔
161
  double teeth_offset = prn(&seed) * teeth_distance;
96,023✔
162

163
  // First and last hitting tooth
164
  int64_t end = start + simulation::fission_bank.size();
96,023✔
165
  int64_t tooth_start = std::ceil((start - teeth_offset) / teeth_distance);
96,023✔
166
  int64_t tooth_end = std::floor((end - teeth_offset) / teeth_distance) + 1;
96,023✔
167

168
  // Locally comb particles in fission_bank
169
  double tooth = tooth_start * teeth_distance + teeth_offset;
96,023✔
170
  for (int64_t i = tooth_start; i < tooth_end; i++) {
151,133,023✔
171
    int64_t idx = std::floor(tooth) - start;
151,037,000✔
172
    temp_sites[index_temp] = simulation::fission_bank[idx];
151,037,000✔
173
    if (settings::ifp_on) {
151,037,000✔
174
      copy_ifp_data_from_fission_banks(
1,320,000✔
175
        idx, temp_delayed_groups[index_temp], temp_lifetimes[index_temp]);
1,320,000✔
176
    }
177
    ++index_temp;
151,037,000✔
178

179
    // Next tooth
180
    tooth += teeth_distance;
151,037,000✔
181
  }
182

183
  // At this point, the sampling of source sites is done and now we need to
184
  // figure out where to send source sites. Since it is possible that one
185
  // processor's share of the source bank spans more than just the immediate
186
  // neighboring processors, we have to perform an ALLGATHER to determine the
187
  // indices for all processors
188

189
#ifdef OPENMC_MPI
190
  // First do an exclusive scan to get the starting indices for
191
  start = 0;
43,352✔
192
  MPI_Exscan(&index_temp, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
43,352✔
193
  finish = start + index_temp;
43,352✔
194

195
  // TODO: protect for MPI_Exscan at rank 0
196

197
  // Allocate space for bank_position if this hasn't been done yet
198
  std::vector<int64_t> bank_position(mpi::n_procs);
43,352✔
199
  MPI_Allgather(&start, 1, MPI_INT64_T, bank_position.data(), 1, MPI_INT64_T,
43,352✔
200
    mpi::intracomm);
201
#else
202
  start = 0;
52,671✔
203
  finish = index_temp;
52,671✔
204
#endif
205

206
  simulation::time_bank_sample.stop();
96,023✔
207
  simulation::time_bank_sendrecv.start();
96,023✔
208

209
#ifdef OPENMC_MPI
210
  // ==========================================================================
211
  // SEND BANK SITES TO NEIGHBORS
212

213
  // IFP number of generation
214
  int ifp_n_generation;
43,352✔
215
  if (settings::ifp_on) {
43,352✔
216
    broadcast_ifp_n_generation(
640✔
217
      ifp_n_generation, temp_delayed_groups, temp_lifetimes);
218
  }
219

220
  int64_t index_local = 0;
43,352✔
221
  vector<MPI_Request> requests;
43,352!
222

223
  // IFP send buffers
224
  vector<int> send_delayed_groups;
43,352✔
225
  vector<double> send_lifetimes;
43,352✔
226

227
  if (start < settings::n_particles) {
43,352!
228
    // Determine the index of the processor which has the first part of the
229
    // source_bank for the local processor
230
    int neighbor = upper_bound_index(
43,352✔
231
      simulation::work_index.begin(), simulation::work_index.end(), start);
43,352✔
232

233
    // Resize IFP send buffers
234
    if (settings::ifp_on && mpi::n_procs > 1) {
43,352✔
235
      resize_ifp_data(send_delayed_groups, send_lifetimes,
320✔
236
        ifp_n_generation * 3 * simulation::work_per_rank);
320✔
237
    }
238

239
    while (start < finish) {
63,142✔
240
      // Determine the number of sites to send
241
      int64_t n =
56,409✔
242
        std::min(simulation::work_index[neighbor + 1], finish) - start;
56,409✔
243

244
      // Initiate an asynchronous send of source sites to the neighboring
245
      // process
246
      if (neighbor != mpi::rank) {
56,409✔
247
        requests.emplace_back();
13,057✔
248
        MPI_Isend(&temp_sites[index_local], static_cast<int>(n),
13,057✔
249
          mpi::source_site, neighbor, mpi::rank, mpi::intracomm,
250
          &requests.back());
13,057✔
251

252
        if (settings::ifp_on) {
13,057✔
253
          // Send IFP data
254
          if (is_beta_effective_or_both())
160!
255
            send_ifp_info(index_local, n, ifp_n_generation, neighbor, requests,
160✔
256
              temp_delayed_groups, send_delayed_groups);
257
          if (is_generation_time_or_both())
160!
258
            send_ifp_info(index_local, n, ifp_n_generation, neighbor, requests,
160✔
259
              temp_lifetimes, send_lifetimes);
260
        }
261
      }
262

263
      // Increment all indices
264
      start += n;
56,409✔
265
      index_local += n;
56,409✔
266
      ++neighbor;
56,409✔
267

268
      // Check for sites out of bounds -- this only happens in the rare
269
      // circumstance that a processor close to the end has so many sites that
270
      // it would exceed the bank on the last processor
271
      if (neighbor > mpi::n_procs - 1)
56,409✔
272
        break;
273
    }
274
  }
275

276
  // ==========================================================================
277
  // RECEIVE BANK SITES FROM NEIGHBORS OR TEMPORARY BANK
278

279
  start = simulation::work_index[mpi::rank];
43,352✔
280
  index_local = 0;
43,352✔
281

282
  // IFP receive buffers
283
  vector<int> recv_delayed_groups;
43,352✔
284
  vector<double> recv_lifetimes;
43,352✔
285
  vector<DeserializationInfo> deserialization_info;
43,352✔
286

287
  // Determine what process has the source sites that will need to be stored at
288
  // the beginning of this processor's source bank.
289

290
  int neighbor;
43,352✔
291
  if (start >= bank_position[mpi::n_procs - 1]) {
43,352✔
292
    neighbor = mpi::n_procs - 1;
293
  } else {
294
    neighbor =
19,817✔
295
      upper_bound_index(bank_position.begin(), bank_position.end(), start);
19,817✔
296
  }
297

298
  // Resize IFP receive buffers
299
  if (settings::ifp_on && mpi::n_procs > 1) {
43,352✔
300
    resize_ifp_data(recv_delayed_groups, recv_lifetimes,
320✔
301
      ifp_n_generation * simulation::work_per_rank);
320✔
302
  }
303

304
  while (start < simulation::work_index[mpi::rank + 1]) {
99,761✔
305
    // Determine how many sites need to be received
306
    int64_t n;
56,409✔
307
    if (neighbor == mpi::n_procs - 1) {
56,409✔
308
      n = simulation::work_index[mpi::rank + 1] - start;
36,592✔
309
    } else {
310
      n = std::min(bank_position[neighbor + 1],
19,817✔
311
            simulation::work_index[mpi::rank + 1]) -
19,817✔
312
          start;
313
    }
314

315
    if (neighbor != mpi::rank) {
56,409✔
316
      // If the source sites are not on this processor, initiate an
317
      // asynchronous receive for the source sites
318

319
      requests.emplace_back();
13,057✔
320
      MPI_Irecv(&simulation::source_bank[index_local], static_cast<int>(n),
13,057✔
321
        mpi::source_site, neighbor, neighbor, mpi::intracomm, &requests.back());
13,057✔
322

323
      if (settings::ifp_on) {
13,057✔
324
        // Receive IFP data
325
        if (is_beta_effective_or_both())
160!
326
          receive_ifp_data(index_local, n, ifp_n_generation, neighbor, requests,
160✔
327
            recv_delayed_groups, deserialization_info);
328
        if (is_generation_time_or_both())
160!
329
          receive_ifp_data(index_local, n, ifp_n_generation, neighbor, requests,
160✔
330
            recv_lifetimes, deserialization_info);
331
      }
332

333
    } else {
334
      // If the source sites are on this processor, we can simply copy them
335
      // from the temp_sites bank
336

337
      index_temp = start - bank_position[mpi::rank];
43,352✔
338
      std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n],
43,352✔
339
        &simulation::source_bank[index_local]);
43,352✔
340

341
      if (settings::ifp_on) {
43,352✔
342
        copy_partial_ifp_data_to_source_banks(
640✔
343
          index_temp, n, index_local, temp_delayed_groups, temp_lifetimes);
344
      }
345
    }
346

347
    // Increment all indices
348
    start += n;
56,409✔
349
    index_local += n;
56,409✔
350
    ++neighbor;
56,409✔
351
  }
352

353
  // Since we initiated a series of asynchronous ISENDs and IRECVs, now we have
354
  // to ensure that the data has actually been communicated before moving on to
355
  // the next generation
356

357
  int n_request = requests.size();
43,352✔
358
  MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
43,352✔
359

360
  if (settings::ifp_on) {
43,352✔
361
    if (is_beta_effective_or_both())
640!
362
      deserialize_ifp_info(ifp_n_generation, recv_delayed_groups,
640✔
363
        simulation::ifp_source_delayed_group_bank, deserialization_info);
364
    if (is_generation_time_or_both())
640!
365
      deserialize_ifp_info(ifp_n_generation, recv_lifetimes,
640✔
366
        simulation::ifp_source_lifetime_bank, deserialization_info);
367
  }
368

369
#else
370
  std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles,
52,671✔
371
    simulation::source_bank.begin());
372
  if (settings::ifp_on) {
52,671✔
373
    copy_complete_ifp_data_to_source_banks(temp_delayed_groups, temp_lifetimes);
840✔
374
  }
375
#endif
376

377
  simulation::time_bank_sendrecv.stop();
96,023✔
378
  simulation::time_bank.stop();
96,023✔
379
}
96,023✔
380

381
void calculate_average_keff()
102,873✔
382
{
383
  // Determine overall generation and number of active generations
384
  int i = overall_generation() - 1;
102,873✔
385
  int n;
102,873✔
386
  if (simulation::current_batch > settings::n_inactive) {
102,873✔
387
    n = settings::gen_per_batch * simulation::n_realizations +
77,836✔
388
        simulation::current_gen;
389
  } else {
390
    n = 0;
391
  }
392

393
  if (n <= 0) {
77,836!
394
    // For inactive generations, use current generation k as estimate for next
395
    // generation
396
    simulation::keff = simulation::k_generation[i];
25,037✔
397
  } else {
398
    // Sample mean of keff
399
    simulation::k_sum[0] += simulation::k_generation[i];
77,836✔
400
    simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
77,836✔
401

402
    // Determine mean
403
    simulation::keff = simulation::k_sum[0] / n;
77,836✔
404

405
    if (n > 1) {
77,836✔
406
      double t_value;
73,235✔
407
      if (settings::confidence_intervals) {
73,235✔
408
        // Calculate t-value for confidence intervals
409
        double alpha = 1.0 - CONFIDENCE_LEVEL;
105✔
410
        t_value = t_percentile(1.0 - alpha / 2.0, n - 1);
105✔
411
      } else {
412
        t_value = 1.0;
413
      }
414

415
      // Standard deviation of the sample mean of k
416
      simulation::keff_std =
146,470✔
417
        t_value *
73,235✔
418
        std::sqrt(
73,235✔
419
          (simulation::k_sum[1] / n - std::pow(simulation::keff, 2)) / (n - 1));
73,235✔
420

421
      // In some cases (such as an infinite medium problem), random ray
422
      // may estimate k exactly and in an unvarying manner between iterations.
423
      // In this case, the floating point roundoff between the division and the
424
      // power operations may cause an extremely small negative value to occur
425
      // inside the sqrt operation, leading to NaN. If this occurs, we check for
426
      // it and set the std dev to zero.
427
      if (!std::isfinite(simulation::keff_std)) {
73,235!
428
        simulation::keff_std = 0.0;
×
429
      }
430
    }
431
  }
432
}
102,873✔
433

434
int openmc_get_keff(double* k_combined)
13,652✔
435
{
436
  k_combined[0] = 0.0;
13,652✔
437
  k_combined[1] = 0.0;
13,652✔
438

439
  // Special case for n <=3. Notice that at the end,
440
  // there is a N-3 term in a denominator.
441
  if (simulation::n_realizations <= 3 ||
13,652✔
442
      settings::solver_type == SolverType::RANDOM_RAY) {
10,715✔
443
    k_combined[0] = simulation::keff;
3,223✔
444
    k_combined[1] = simulation::keff_std;
3,223✔
445
    if (simulation::n_realizations <= 1) {
3,223✔
446
      k_combined[1] = std::numeric_limits<double>::infinity();
2,167✔
447
    }
448
    return 0;
3,223✔
449
  }
450

451
  // Initialize variables
452
  int64_t n = simulation::n_realizations;
10,429✔
453

454
  // Copy estimates of k-effective and its variance (not variance of the mean)
455
  const auto& gt = simulation::global_tallies;
10,429✔
456

457
  array<double, 3> kv {};
10,429✔
458
  tensor::Tensor<double> cov = tensor::zeros<double>({3, 3});
10,429✔
459
  kv[0] = gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n;
10,429✔
460
  kv[1] = gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n;
10,429✔
461
  kv[2] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n;
10,429✔
462
  cov(0, 0) =
10,429✔
463
    (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) /
10,429✔
464
    (n - 1);
10,429✔
465
  cov(1, 1) =
10,429✔
466
    (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n * kv[1] * kv[1]) /
10,429✔
467
    (n - 1);
468
  cov(2, 2) =
10,429✔
469
    (gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM_SQ) - n * kv[2] * kv[2]) /
10,429✔
470
    (n - 1);
471

472
  // Calculate covariances based on sums with Bessel's correction
473
  cov(0, 1) = (simulation::k_col_abs - n * kv[0] * kv[1]) / (n - 1);
10,429✔
474
  cov(0, 2) = (simulation::k_col_tra - n * kv[0] * kv[2]) / (n - 1);
10,429✔
475
  cov(1, 2) = (simulation::k_abs_tra - n * kv[1] * kv[2]) / (n - 1);
10,429✔
476
  cov(1, 0) = cov(0, 1);
10,429✔
477
  cov(2, 0) = cov(0, 2);
10,429✔
478
  cov(2, 1) = cov(1, 2);
10,429✔
479

480
  // Check to see if two estimators are the same; this is guaranteed to happen
481
  // in MG-mode with survival biasing when the collision and absorption
482
  // estimators are the same, but can theoretically happen at anytime.
483
  // If it does, the standard estimators will produce floating-point
484
  // exceptions and an expression specifically derived for the combination of
485
  // two estimators (vice three) should be used instead.
486

487
  // If delta tracking is enabled, use the collision and absorption estimators
488
  // only. Otherwise, we will identify if there are any matching estimators. If
489
  // none match, all three estimates are used.
490
  int i, j;
10,429✔
491
  bool use_three = false;
10,429✔
492
  if (settings::delta_tracking) {
10,429✔
493
    i = 0;
494
    j = 1;
495
  } else {
496
    if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) &&
10,209!
497
        (std::abs(cov(0, 0) - cov(1, 1)) / cov(0, 0) < FP_REL_PRECISION)) {
22!
498
      // 0 and 1 match, so only use 0 and 2 in our comparisons
499
      i = 0;
500
      j = 2;
501

502
    } else if ((std::abs(kv[0] - kv[2]) / kv[0] < FP_REL_PRECISION) &&
10,187!
NEW
503
               (std::abs(cov(0, 0) - cov(2, 2)) / cov(0, 0) <
×
504
                 FP_REL_PRECISION)) {
505
      // 0 and 2 match, so only use 0 and 1 in our comparisons
506
      i = 0;
507
      j = 1;
508

509
    } else if ((std::abs(kv[1] - kv[2]) / kv[1] < FP_REL_PRECISION) &&
10,187!
NEW
510
               (std::abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) <
×
511
                 FP_REL_PRECISION)) {
512
      // 1 and 2 match, so only use 0 and 1 in our comparisons
513
      i = 0;
514
      j = 1;
515

516
    } else {
517
      // No two estimators match, so set boolean to use all three estimators.
518
      use_three = true;
10,187✔
519
    }
520
  }
521

522
  if (use_three) {
10,187✔
523
    // Use three estimators as derived in the paper by Urbatsch
524

525
    // Initialize variables
526
    double g = 0.0;
10,187✔
527
    array<double, 3> S {};
10,187✔
528

529
    for (int l = 0; l < 3; ++l) {
40,748✔
530
      // Permutations of estimates
531
      int k;
30,561✔
532
      switch (l) {
30,561✔
533
      case 0:
534
        // i = collision, j = absorption, k = tracklength
535
        i = 0;
536
        j = 1;
537
        k = 2;
538
        break;
539
      case 1:
10,187✔
540
        // i = absortion, j = tracklength, k = collision
541
        i = 1;
10,187✔
542
        j = 2;
10,187✔
543
        k = 0;
10,187✔
544
        break;
10,187✔
545
      case 2:
10,187✔
546
        // i = tracklength, j = collision, k = absorption
547
        i = 2;
10,187✔
548
        j = 0;
10,187✔
549
        k = 1;
10,187✔
550
        break;
10,187✔
551
      }
552

553
      // Calculate weighting
554
      double f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) +
30,561✔
555
                 cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k));
30,561✔
556

557
      // Add to S sums for variance of combined estimate
558
      S[0] += f * cov(0, l);
30,561✔
559
      S[1] += (cov(j, j) + cov(k, k) - 2.0 * cov(j, k)) * kv[l] * kv[l];
30,561✔
560
      S[2] += (cov(k, k) + cov(i, j) - cov(j, k) - cov(i, k)) * kv[l] * kv[j];
30,561✔
561

562
      // Add to sum for combined k-effective
563
      k_combined[0] += f * kv[l];
30,561✔
564
      g += f;
30,561✔
565
    }
566

567
    // Complete calculations of S sums
568
    for (auto& S_i : S) {
40,748✔
569
      S_i *= (n - 1);
30,561✔
570
    }
571
    S[0] *= (n - 1) * (n - 1);
10,187✔
572

573
    // Calculate combined estimate of k-effective
574
    k_combined[0] /= g;
10,187✔
575

576
    // Calculate standard deviation of combined estimate
577
    g *= (n - 1) * (n - 1);
10,187✔
578
    k_combined[1] =
10,187✔
579
      std::sqrt(S[0] / (g * n * (n - 3)) * (1 + n * ((S[1] - 2 * S[2]) / g)));
10,187✔
580

581
  } else {
582
    // Use only two estimators
583
    // These equations are derived analogously to that done in the paper by
584
    // Urbatsch, but are simpler than for the three estimators case since the
585
    // block matrices of the three estimator equations reduces to scalars here
586

587
    // Store common terms.
588
    const double f = cov(i, i) + cov(j, j) - 2.0 * cov(i, j);
242✔
589
    const double g = cov(i, i) - cov(i, j);
242✔
590
    const double h = kv[i] - kv[j];
242✔
591

592
    // Eq. 36 in LA-12658-MS
593
    k_combined[0] = kv[i] - g * h / f;
242✔
594

595
    // This is re-derived from the \Sigma matrix based on Eq. 40 in LA-12658-MS
596
    k_combined[1] =
242✔
597
      (cov(i, i) - g * g / f) * ((n - 1.0) / n + h * h / f) / (n - 2.0);
242✔
598
    k_combined[1] = std::sqrt(k_combined[1]);
242✔
599
  }
600
  return 0;
10,429✔
601
}
13,652✔
602

603
void shannon_entropy()
7,685✔
604
{
605
  // Get source weight in each mesh bin
606
  bool sites_outside;
7,685✔
607
  tensor::Tensor<double> p =
7,685✔
608
    simulation::entropy_mesh->count_sites(simulation::fission_bank.data(),
7,685✔
609
      simulation::fission_bank.size(), &sites_outside);
7,685✔
610

611
  // display warning message if there were sites outside entropy box
612
  if (sites_outside) {
7,685!
613
    if (mpi::master)
×
614
      warning("Fission source site(s) outside of entropy box.");
×
615
  }
616

617
  if (mpi::master) {
7,685✔
618
    // Normalize to total weight of bank sites
619
    p /= p.sum();
7,645✔
620

621
    // Sum values to obtain Shannon entropy
622
    double H = 0.0;
7,645✔
623
    for (auto p_i : p) {
586,355✔
624
      if (p_i > 0.0) {
578,710✔
625
        H -= p_i * std::log2(p_i);
453,695✔
626
      }
627
    }
628

629
    // Add value to vector
630
    simulation::entropy.push_back(H);
7,645✔
631
  }
632
}
7,685✔
633

634
void ufs_count_sites()
150✔
635
{
636
  if (simulation::current_batch == 1 && simulation::current_gen == 1) {
150!
637
    // On the first generation, just assume that the source is already evenly
638
    // distributed so that effectively the production of fission sites is not
639
    // biased
640

641
    std::size_t n = simulation::ufs_mesh->n_bins();
15✔
642
    double vol_frac = simulation::ufs_mesh->volume_frac_;
15✔
643
    simulation::source_frac = tensor::Tensor<double>({n}, vol_frac);
15✔
644

645
  } else {
15✔
646
    // count number of source sites in each ufs mesh cell
647
    bool sites_outside;
135✔
648
    simulation::source_frac =
135✔
649
      simulation::ufs_mesh->count_sites(simulation::source_bank.data(),
135✔
650
        simulation::source_bank.size(), &sites_outside);
135✔
651

652
    // Check for sites outside of the mesh
653
    if (mpi::master && sites_outside) {
135!
654
      fatal_error("Source sites outside of the UFS mesh!");
×
655
    }
656

657
#ifdef OPENMC_MPI
658
    // Send source fraction to all processors
659
    int n_bins = simulation::ufs_mesh->n_bins();
72✔
660
    MPI_Bcast(
72✔
661
      simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
72✔
662
#endif
663

664
    // Normalize to total weight to get fraction of source in each cell
665
    double total = simulation::source_frac.sum();
135✔
666
    simulation::source_frac /= total;
135✔
667

668
    // Since the total starting weight is not equal to n_particles, we need to
669
    // renormalize the weight of the source sites
670
    for (int i = 0; i < simulation::work_per_rank; ++i) {
99,135✔
671
      simulation::source_bank[i].wgt *= settings::n_particles / total;
99,000✔
672
    }
673
  }
674
}
150✔
675

676
double ufs_get_weight(const Particle& p)
94,182✔
677
{
678
  // Determine indices on ufs mesh for current location
679
  int mesh_bin = simulation::ufs_mesh->get_bin(p.r());
94,182✔
680
  if (mesh_bin < 0) {
94,182!
681
    p.write_restart();
×
682
    fatal_error("Source site outside UFS mesh!");
×
683
  }
684

685
  if (simulation::source_frac(mesh_bin) != 0.0) {
94,182✔
686
    return simulation::ufs_mesh->volume_frac_ /
51,095✔
687
           simulation::source_frac(mesh_bin);
51,095✔
688
  } else {
689
    return 1.0;
690
  }
691
}
692

693
void write_eigenvalue_hdf5(hid_t group)
4,739✔
694
{
695
  write_dataset(group, "n_inactive", settings::n_inactive);
4,739✔
696
  write_dataset(group, "generations_per_batch", settings::gen_per_batch);
4,739✔
697
  write_dataset(group, "k_generation", simulation::k_generation);
4,739✔
698
  if (settings::entropy_on) {
4,739✔
699
    write_dataset(group, "entropy", simulation::entropy);
550✔
700
  }
701
  write_dataset(group, "k_col_abs", simulation::k_col_abs);
4,739✔
702
  write_dataset(group, "k_col_tra", simulation::k_col_tra);
4,739✔
703
  write_dataset(group, "k_abs_tra", simulation::k_abs_tra);
4,739✔
704
  array<double, 2> k_combined;
4,739✔
705
  openmc_get_keff(k_combined.data());
4,739✔
706
  write_dataset(group, "k_combined", k_combined);
4,739✔
707
}
4,739✔
708

709
void read_eigenvalue_hdf5(hid_t group)
63✔
710
{
711
  read_dataset(group, "generations_per_batch", settings::gen_per_batch);
63✔
712
  int n = simulation::restart_batch * settings::gen_per_batch;
63✔
713
  simulation::k_generation.resize(n);
63✔
714
  read_dataset(group, "k_generation", simulation::k_generation);
63✔
715
  if (settings::entropy_on) {
63✔
716
    read_dataset(group, "entropy", simulation::entropy);
11✔
717
  }
718
  read_dataset(group, "k_col_abs", simulation::k_col_abs);
63✔
719
  read_dataset(group, "k_col_tra", simulation::k_col_tra);
63✔
720
  read_dataset(group, "k_abs_tra", simulation::k_abs_tra);
63✔
721
}
63✔
722

723
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE TRIAL · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc