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

openmc-dev / openmc / 27444237628

12 Jun 2026 09:31PM UTC coverage: 81.304% (+0.02%) from 81.281%
27444237628

Pull #3971

github

web-flow
Merge 9487507b3 into 02eb999af
Pull Request #3971: Delta tracking

18480 of 26797 branches covered (68.96%)

Branch coverage included in aggregate %.

592 of 644 new or added lines in 20 files covered. (91.93%)

31 existing lines in 3 files now uncovered.

59819 of 69507 relevant lines covered (86.06%)

49514105.74 hits per line

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

90.91
/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,653✔
48
{
49
  const auto& gt = simulation::global_tallies;
102,653✔
50

51
  // Get keff for this generation by subtracting off the starting value
52
  if (settings::delta_tracking) {
102,653✔
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,153✔
58
      gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) -
101,153✔
59
      simulation::keff_generation;
60
  }
61

62
  double keff_reduced;
102,653✔
63
#ifdef OPENMC_MPI
64
  if (settings::solver_type != SolverType::RANDOM_RAY) {
46,552✔
65
    // Combine values across all processors
66
    MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE,
43,272✔
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,101✔
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,653✔
83
    keff_reduced /= settings::n_particles;
95,803✔
84
  }
85

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

89
void synchronize_bank()
95,803✔
90
{
91
  simulation::time_bank.start();
95,803✔
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,272✔
104
  int64_t n_bank = simulation::fission_bank.size();
43,272✔
105
  MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
43,272✔
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,272✔
111
    start = 0;
29,997✔
112

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

117
#else
118
  int64_t start = 0;
52,531✔
119
  int64_t finish = simulation::fission_bank.size();
52,531!
120
  int64_t total = finish;
52,531✔
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) {
95,803!
129
    fatal_error(
×
130
      "No fission sites banked on MPI rank " + std::to_string(mpi::rank));
×
131
  }
132

133
  simulation::time_bank_sample.start();
95,803✔
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;
95,803✔
138

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

141
  // Temporary banks for IFP
142
  vector<vector<int>> temp_delayed_groups;
95,803✔
143
  vector<vector<double>> temp_lifetimes;
95,803✔
144
  if (settings::ifp_on) {
95,803✔
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();
95,803✔
157
  uint64_t seed = init_seed(id, STREAM_TRACKING);
95,803✔
158

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

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

168
  // Locally comb particles in fission_bank
169
  double tooth = tooth_start * teeth_distance + teeth_offset;
95,803✔
170
  for (int64_t i = tooth_start; i < tooth_end; i++) {
151,099,803✔
171
    int64_t idx = std::floor(tooth) - start;
151,004,000✔
172
    temp_sites[index_temp] = simulation::fission_bank[idx];
151,004,000✔
173
    if (settings::ifp_on) {
151,004,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,004,000✔
178

179
    // Next tooth
180
    tooth += teeth_distance;
151,004,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,272✔
192
  MPI_Exscan(&index_temp, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
43,272✔
193
  finish = start + index_temp;
43,272✔
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,272✔
199
  MPI_Allgather(&start, 1, MPI_INT64_T, bank_position.data(), 1, MPI_INT64_T,
43,272✔
200
    mpi::intracomm);
201
#else
202
  start = 0;
52,531✔
203
  finish = index_temp;
52,531✔
204
#endif
205

206
  simulation::time_bank_sample.stop();
95,803✔
207
  simulation::time_bank_sendrecv.start();
95,803✔
208

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

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

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

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

227
  if (start < settings::n_particles) {
43,272!
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,272✔
231
      simulation::work_index.begin(), simulation::work_index.end(), start);
43,272✔
232

233
    // Resize IFP send buffers
234
    if (settings::ifp_on && mpi::n_procs > 1) {
43,272✔
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,054✔
240
      // Determine the number of sites to send
241
      int64_t n =
56,333✔
242
        std::min(simulation::work_index[neighbor + 1], finish) - start;
56,333✔
243

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

252
        if (settings::ifp_on) {
13,061✔
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,333✔
265
      index_local += n;
56,333✔
266
      ++neighbor;
56,333✔
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,333✔
272
        break;
273
    }
274
  }
275

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

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

282
  // IFP receive buffers
283
  vector<int> recv_delayed_groups;
43,272✔
284
  vector<double> recv_lifetimes;
43,272✔
285
  vector<DeserializationInfo> deserialization_info;
43,272✔
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,272✔
291
  if (start >= bank_position[mpi::n_procs - 1]) {
43,272✔
292
    neighbor = mpi::n_procs - 1;
293
  } else {
294
    neighbor =
19,829✔
295
      upper_bound_index(bank_position.begin(), bank_position.end(), start);
19,829✔
296
  }
297

298
  // Resize IFP receive buffers
299
  if (settings::ifp_on && mpi::n_procs > 1) {
43,272✔
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,605✔
305
    // Determine how many sites need to be received
306
    int64_t n;
56,333✔
307
    if (neighbor == mpi::n_procs - 1) {
56,333✔
308
      n = simulation::work_index[mpi::rank + 1] - start;
36,504✔
309
    } else {
310
      n = std::min(bank_position[neighbor + 1],
19,829✔
311
            simulation::work_index[mpi::rank + 1]) -
19,829✔
312
          start;
313
    }
314

315
    if (neighbor != mpi::rank) {
56,333✔
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,061✔
320
      MPI_Irecv(&simulation::source_bank[index_local], static_cast<int>(n),
13,061✔
321
        mpi::source_site, neighbor, neighbor, mpi::intracomm, &requests.back());
13,061✔
322

323
      if (settings::ifp_on) {
13,061✔
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,272✔
338
      std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n],
43,272✔
339
        &simulation::source_bank[index_local]);
43,272✔
340

341
      if (settings::ifp_on) {
43,272✔
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,333✔
349
    index_local += n;
56,333✔
350
    ++neighbor;
56,333✔
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,272✔
358
  MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
43,272✔
359

360
  if (settings::ifp_on) {
43,272✔
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,531✔
371
    simulation::source_bank.begin());
372
  if (settings::ifp_on) {
52,531✔
373
    copy_complete_ifp_data_to_source_banks(temp_delayed_groups, temp_lifetimes);
840✔
374
  }
375
#endif
376

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

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

393
  if (n <= 0) {
77,616!
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,616✔
400
    simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
77,616✔
401

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

405
    if (n > 1) {
77,616✔
406
      double t_value;
73,037✔
407
      if (settings::confidence_intervals) {
73,037✔
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,074✔
417
        t_value *
73,037✔
418
        std::sqrt(
73,037✔
419
          (simulation::k_sum[1] / n - std::pow(simulation::keff, 2)) / (n - 1));
73,037✔
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,037!
428
        simulation::keff_std = 0.0;
×
429
      }
430
    }
431
  }
432
}
102,653✔
433

434
int openmc_get_keff(double* k_combined)
13,564✔
435
{
436
  k_combined[0] = 0.0;
13,564✔
437
  k_combined[1] = 0.0;
13,564✔
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,564✔
442
      settings::solver_type == SolverType::RANDOM_RAY) {
10,627✔
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,341✔
453

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

457
  array<double, 3> kv {};
10,341✔
458
  tensor::Tensor<double> cov = tensor::zeros<double>({3, 3});
10,341✔
459
  kv[0] = gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n;
10,341✔
460
  kv[1] = gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n;
10,341✔
461
  kv[2] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n;
10,341✔
462
  cov(0, 0) =
10,341✔
463
    (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) /
10,341✔
464
    (n - 1);
10,341✔
465
  cov(1, 1) =
10,341✔
466
    (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n * kv[1] * kv[1]) /
10,341✔
467
    (n - 1);
468
  cov(2, 2) =
10,341✔
469
    (gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM_SQ) - n * kv[2] * kv[2]) /
10,341✔
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,341✔
474
  cov(0, 2) = (simulation::k_col_tra - n * kv[0] * kv[2]) / (n - 1);
10,341✔
475
  cov(1, 2) = (simulation::k_abs_tra - n * kv[1] * kv[2]) / (n - 1);
10,341✔
476
  cov(1, 0) = cov(0, 1);
10,341✔
477
  cov(2, 0) = cov(0, 2);
10,341✔
478
  cov(2, 1) = cov(1, 2);
10,341✔
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,341✔
491
  bool use_three = false;
10,341✔
492
  if (settings::delta_tracking) {
10,341✔
493
    i = 0;
494
    j = 1;
495
  } else {
496
    if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) &&
10,121!
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,099!
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,099!
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,099✔
519
    }
520
  }
521

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

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

529
    for (int l = 0; l < 3; ++l) {
40,396✔
530
      // Permutations of estimates
531
      int k;
30,297✔
532
      switch (l) {
30,297✔
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,099✔
540
        // i = absortion, j = tracklength, k = collision
541
        i = 1;
10,099✔
542
        j = 2;
10,099✔
543
        k = 0;
10,099✔
544
        break;
10,099✔
545
      case 2:
10,099✔
546
        // i = tracklength, j = collision, k = absorption
547
        i = 2;
10,099✔
548
        j = 0;
10,099✔
549
        k = 1;
10,099✔
550
        break;
10,099✔
551
      }
552

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

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

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

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

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

576
    // Calculate standard deviation of combined estimate
577
    g *= (n - 1) * (n - 1);
10,099✔
578
    k_combined[1] =
10,099✔
579
      std::sqrt(S[0] / (g * n * (n - 3)) * (1 + n * ((S[1] - 2 * S[2]) / g)));
10,099✔
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 the commonly used term
588
    double f = kv[i] - kv[j];
242✔
589
    double g = cov(i, i) + cov(j, j) - 2.0 * cov(i, j);
242✔
590

591
    // Calculate combined estimate of k-effective
592
    k_combined[0] = kv[i] - (cov(i, i) - cov(i, j)) / g * f;
242✔
593

594
    // Calculate standard deviation of combined estimate
595
    k_combined[1] = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) *
242✔
596
                    (g + n * f * f) / (n * (n - 2) * g * g);
242✔
597
    k_combined[1] = std::sqrt(k_combined[1]);
242✔
598
  }
599
  return 0;
10,341✔
600
}
13,564✔
601

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

© 2026 Coveralls, Inc