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

openmc-dev / openmc / 21139931957

19 Jan 2026 01:50PM UTC coverage: 81.835% (-0.2%) from 82.058%
21139931957

Pull #2693

github

web-flow
Merge 1258e263b into 5847b0de2
Pull Request #2693: Add reactivity control to coupled transport-depletion analyses

17180 of 23968 branches covered (71.68%)

Branch coverage included in aggregate %.

78 of 84 new or added lines in 4 files covered. (92.86%)

345 existing lines in 27 files now uncovered.

55616 of 64987 relevant lines covered (85.58%)

50201535.92 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 "xtensor/xbuilder.hpp"
4
#include "xtensor/xmath.hpp"
5
#include "xtensor/xtensor.hpp"
6
#include "xtensor/xview.hpp"
7

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

25
#include <algorithm> // for min
26
#include <cmath>     // for sqrt, abs, pow
27
#include <iterator>  // for back_inserter
28
#include <limits>    //for infinity
29
#include <string>
30

31
namespace openmc {
32

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

37
namespace simulation {
38

39
double keff_generation;
40
array<double, 2> k_sum;
41
vector<double> entropy;
42
xt::xtensor<double, 1> source_frac;
43

44
} // namespace simulation
45

46
//==============================================================================
47
// Non-member functions
48
//==============================================================================
49

50
void calculate_generation_keff()
164,557✔
51
{
52
  const auto& gt = simulation::global_tallies;
164,557✔
53

54
  // Get keff for this generation by subtracting off the starting value
55
  simulation::keff_generation =
164,557✔
56
    gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) -
164,557✔
57
    simulation::keff_generation;
58

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

77
  // Normalize single batch estimate of k
78
  // TODO: This should be normalized by total_weight, not by n_particles
79
  if (settings::solver_type != SolverType::RANDOM_RAY) {
164,557✔
80
    keff_reduced /= settings::n_particles;
160,957✔
81
  }
82

83
  simulation::k_generation.push_back(keff_reduced);
164,557✔
84
}
164,557✔
85

86
void synchronize_bank()
160,957✔
87
{
88
  simulation::time_bank.start();
160,957✔
89

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

99
#ifdef OPENMC_MPI
100
  int64_t start = 0;
89,015✔
101
  int64_t n_bank = simulation::fission_bank.size();
89,015✔
102
  MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
89,015✔
103

104
  // While we would expect the value of start on rank 0 to be 0, the MPI
105
  // standard says that the receive buffer on rank 0 is undefined and not
106
  // significant
107
  if (mpi::rank == 0)
89,015✔
108
    start = 0;
73,055✔
109

110
  int64_t finish = start + simulation::fission_bank.size();
89,015✔
111
  int64_t total = finish;
89,015✔
112
  MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm);
89,015✔
113

114
#else
115
  int64_t start = 0;
71,942✔
116
  int64_t finish = simulation::fission_bank.size();
71,942✔
117
  int64_t total = finish;
71,942✔
118
#endif
119

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

125
  if (simulation::fission_bank.size() == 0) {
160,957!
126
    fatal_error(
×
127
      "No fission sites banked on MPI rank " + std::to_string(mpi::rank));
×
128
  }
129

130
  simulation::time_bank_sample.start();
160,957✔
131

132
  // Allocate temporary source bank -- we don't really know how many fission
133
  // sites were created, so overallocate by a factor of 3
134
  int64_t index_temp = 0;
160,957✔
135

136
  vector<SourceSite> temp_sites(3 * simulation::work_per_rank);
160,957✔
137

138
  // Temporary banks for IFP
139
  vector<vector<int>> temp_delayed_groups;
160,957✔
140
  vector<vector<double>> temp_lifetimes;
160,957✔
141
  if (settings::ifp_on) {
160,957✔
142
    resize_ifp_data(
1,400✔
143
      temp_delayed_groups, temp_lifetimes, 3 * simulation::work_per_rank);
144
  }
145

146
  // ==========================================================================
147
  // SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
148

149
  // We use Uniform Combing method to exactly get the targeted particle size
150
  // [https://doi.org/10.1080/00295639.2022.2091906]
151

152
  // Make sure all processors use the same random number seed.
153
  int64_t id = simulation::total_gen + overall_generation();
160,957✔
154
  uint64_t seed = init_seed(id, STREAM_TRACKING);
160,957✔
155

156
  // Comb specification
157
  double teeth_distance = static_cast<double>(total) / settings::n_particles;
160,957✔
158
  double teeth_offset = prn(&seed) * teeth_distance;
160,957✔
159

160
  // First and last hitting tooth
161
  int64_t end = start + simulation::fission_bank.size();
160,957✔
162
  int64_t tooth_start = std::ceil((start - teeth_offset) / teeth_distance);
160,957✔
163
  int64_t tooth_end = std::floor((end - teeth_offset) / teeth_distance) + 1;
160,957✔
164

165
  // Locally comb particles in fission_bank
166
  double tooth = tooth_start * teeth_distance + teeth_offset;
160,957✔
167
  for (int64_t i = tooth_start; i < tooth_end; i++) {
212,915,457✔
168
    int64_t idx = std::floor(tooth) - start;
212,754,500✔
169
    temp_sites[index_temp] = simulation::fission_bank[idx];
212,754,500✔
170
    if (settings::ifp_on) {
212,754,500✔
171
      copy_ifp_data_from_fission_banks(
1,200,000✔
172
        idx, temp_delayed_groups[index_temp], temp_lifetimes[index_temp]);
1,200,000✔
173
    }
174
    ++index_temp;
212,754,500✔
175

176
    // Next tooth
177
    tooth += teeth_distance;
212,754,500✔
178
  }
179

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

186
#ifdef OPENMC_MPI
187
  // First do an exclusive scan to get the starting indices for
188
  start = 0;
89,015✔
189
  MPI_Exscan(&index_temp, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
89,015✔
190
  finish = start + index_temp;
89,015✔
191

192
  // TODO: protect for MPI_Exscan at rank 0
193

194
  // Allocate space for bank_position if this hasn't been done yet
195
  int64_t bank_position[mpi::n_procs];
89,015✔
196
  MPI_Allgather(
89,015✔
197
    &start, 1, MPI_INT64_T, bank_position, 1, MPI_INT64_T, mpi::intracomm);
198
#else
199
  start = 0;
71,942✔
200
  finish = index_temp;
71,942✔
201
#endif
202

203
  simulation::time_bank_sample.stop();
160,957✔
204
  simulation::time_bank_sendrecv.start();
160,957✔
205

206
#ifdef OPENMC_MPI
207
  // ==========================================================================
208
  // SEND BANK SITES TO NEIGHBORS
209

210
  // IFP number of generation
211
  int ifp_n_generation;
212
  if (settings::ifp_on) {
89,015✔
213
    broadcast_ifp_n_generation(
800✔
214
      ifp_n_generation, temp_delayed_groups, temp_lifetimes);
215
  }
216

217
  int64_t index_local = 0;
89,015✔
218
  vector<MPI_Request> requests;
89,015✔
219

220
  // IFP send buffers
221
  vector<int> send_delayed_groups;
89,015✔
222
  vector<double> send_lifetimes;
89,015✔
223

224
  if (start < settings::n_particles) {
89,015!
225
    // Determine the index of the processor which has the first part of the
226
    // source_bank for the local processor
227
    int neighbor = upper_bound_index(
89,015✔
228
      simulation::work_index.begin(), simulation::work_index.end(), start);
89,015✔
229

230
    // Resize IFP send buffers
231
    if (settings::ifp_on && mpi::n_procs > 1) {
89,015✔
232
      resize_ifp_data(send_delayed_groups, send_lifetimes,
400✔
233
        ifp_n_generation * 3 * simulation::work_per_rank);
400✔
234
    }
235

236
    while (start < finish) {
112,840✔
237
      // Determine the number of sites to send
238
      int64_t n =
239
        std::min(simulation::work_index[neighbor + 1], finish) - start;
104,724✔
240

241
      // Initiate an asynchronous send of source sites to the neighboring
242
      // process
243
      if (neighbor != mpi::rank) {
104,724✔
244
        requests.emplace_back();
15,709✔
245
        MPI_Isend(&temp_sites[index_local], static_cast<int>(n),
15,709✔
246
          mpi::source_site, neighbor, mpi::rank, mpi::intracomm,
247
          &requests.back());
15,709✔
248

249
        if (settings::ifp_on) {
15,709✔
250
          // Send IFP data
251
          if (is_beta_effective_or_both())
200!
252
            send_ifp_info(index_local, n, ifp_n_generation, neighbor, requests,
200✔
253
              temp_delayed_groups, send_delayed_groups);
254
          if (is_generation_time_or_both())
200!
255
            send_ifp_info(index_local, n, ifp_n_generation, neighbor, requests,
200✔
256
              temp_lifetimes, send_lifetimes);
257
        }
258
      }
259

260
      // Increment all indices
261
      start += n;
104,724✔
262
      index_local += n;
104,724✔
263
      ++neighbor;
104,724✔
264

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

273
  // ==========================================================================
274
  // RECEIVE BANK SITES FROM NEIGHBORS OR TEMPORARY BANK
275

276
  start = simulation::work_index[mpi::rank];
89,015✔
277
  index_local = 0;
89,015✔
278

279
  // IFP receive buffers
280
  vector<int> recv_delayed_groups;
89,015✔
281
  vector<double> recv_lifetimes;
89,015✔
282
  vector<DeserializationInfo> deserialization_info;
89,015✔
283

284
  // Determine what process has the source sites that will need to be stored at
285
  // the beginning of this processor's source bank.
286

287
  int neighbor;
288
  if (start >= bank_position[mpi::n_procs - 1]) {
89,015✔
289
    neighbor = mpi::n_procs - 1;
65,211✔
290
  } else {
291
    neighbor =
23,804✔
292
      upper_bound_index(bank_position, bank_position + mpi::n_procs, start);
23,804✔
293
  }
294

295
  // Resize IFP receive buffers
296
  if (settings::ifp_on && mpi::n_procs > 1) {
89,015✔
297
    resize_ifp_data(recv_delayed_groups, recv_lifetimes,
400✔
298
      ifp_n_generation * simulation::work_per_rank);
400✔
299
  }
300

301
  while (start < simulation::work_index[mpi::rank + 1]) {
193,739✔
302
    // Determine how many sites need to be received
303
    int64_t n;
304
    if (neighbor == mpi::n_procs - 1) {
104,724✔
305
      n = simulation::work_index[mpi::rank + 1] - start;
80,920✔
306
    } else {
307
      n = std::min(bank_position[neighbor + 1],
23,804✔
308
            simulation::work_index[mpi::rank + 1]) -
47,608✔
309
          start;
310
    }
311

312
    if (neighbor != mpi::rank) {
104,724✔
313
      // If the source sites are not on this processor, initiate an
314
      // asynchronous receive for the source sites
315

316
      requests.emplace_back();
15,709✔
317
      MPI_Irecv(&simulation::source_bank[index_local], static_cast<int>(n),
15,709✔
318
        mpi::source_site, neighbor, neighbor, mpi::intracomm, &requests.back());
15,709✔
319

320
      if (settings::ifp_on) {
15,709✔
321
        // Receive IFP data
322
        if (is_beta_effective_or_both())
200!
323
          receive_ifp_data(index_local, n, ifp_n_generation, neighbor, requests,
200✔
324
            recv_delayed_groups, deserialization_info);
325
        if (is_generation_time_or_both())
200!
326
          receive_ifp_data(index_local, n, ifp_n_generation, neighbor, requests,
200✔
327
            recv_lifetimes, deserialization_info);
328
      }
329

330
    } else {
331
      // If the source sites are on this processor, we can simply copy them
332
      // from the temp_sites bank
333

334
      index_temp = start - bank_position[mpi::rank];
89,015✔
335
      std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n],
89,015✔
336
        &simulation::source_bank[index_local]);
89,015✔
337

338
      if (settings::ifp_on) {
89,015✔
339
        copy_partial_ifp_data_to_source_banks(
800✔
340
          index_temp, n, index_local, temp_delayed_groups, temp_lifetimes);
341
      }
342
    }
343

344
    // Increment all indices
345
    start += n;
104,724✔
346
    index_local += n;
104,724✔
347
    ++neighbor;
104,724✔
348
  }
349

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

354
  int n_request = requests.size();
89,015✔
355
  MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
89,015✔
356

357
  if (settings::ifp_on) {
89,015✔
358
    if (is_beta_effective_or_both())
800!
359
      deserialize_ifp_info(ifp_n_generation, recv_delayed_groups,
800✔
360
        simulation::ifp_source_delayed_group_bank, deserialization_info);
361
    if (is_generation_time_or_both())
800!
362
      deserialize_ifp_info(ifp_n_generation, recv_lifetimes,
800✔
363
        simulation::ifp_source_lifetime_bank, deserialization_info);
364
  }
365

366
#else
367
  std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles,
71,942✔
368
    simulation::source_bank.begin());
369
  if (settings::ifp_on) {
71,942✔
370
    copy_complete_ifp_data_to_source_banks(temp_delayed_groups, temp_lifetimes);
600✔
371
  }
372
#endif
373

374
  simulation::time_bank_sendrecv.stop();
160,957✔
375
  simulation::time_bank.stop();
160,957✔
376
}
249,972✔
377

378
void calculate_average_keff()
164,557✔
379
{
380
  // Determine overall generation and number of active generations
381
  int i = overall_generation() - 1;
164,557✔
382
  int n;
383
  if (simulation::current_batch > settings::n_inactive) {
164,557✔
384
    n = settings::gen_per_batch * simulation::n_realizations +
145,222✔
385
        simulation::current_gen;
386
  } else {
387
    n = 0;
19,335✔
388
  }
389

390
  if (n <= 0) {
164,557✔
391
    // For inactive generations, use current generation k as estimate for next
392
    // generation
393
    simulation::keff = simulation::k_generation[i];
19,335✔
394
  } else {
395
    // Sample mean of keff
396
    simulation::k_sum[0] += simulation::k_generation[i];
145,222✔
397
    simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
145,222✔
398

399
    // Determine mean
400
    simulation::keff = simulation::k_sum[0] / n;
145,222✔
401

402
    if (n > 1) {
145,222✔
403
      double t_value;
404
      if (settings::confidence_intervals) {
141,060✔
405
        // Calculate t-value for confidence intervals
406
        double alpha = 1.0 - CONFIDENCE_LEVEL;
105✔
407
        t_value = t_percentile(1.0 - alpha / 2.0, n - 1);
105✔
408
      } else {
409
        t_value = 1.0;
140,955✔
410
      }
411

412
      // Standard deviation of the sample mean of k
413
      simulation::keff_std =
141,060✔
414
        t_value *
141,060✔
415
        std::sqrt(
141,060✔
416
          (simulation::k_sum[1] / n - std::pow(simulation::keff, 2)) / (n - 1));
141,060✔
417

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

431
int openmc_get_keff(double* k_combined)
12,370✔
432
{
433
  k_combined[0] = 0.0;
12,370✔
434
  k_combined[1] = 0.0;
12,370✔
435

436
  // Special case for n <=3. Notice that at the end,
437
  // there is a N-3 term in a denominator.
438
  if (simulation::n_realizations <= 3 ||
12,370✔
439
      settings::solver_type == SolverType::RANDOM_RAY) {
9,770✔
440
    k_combined[0] = simulation::keff;
2,770✔
441
    k_combined[1] = simulation::keff_std;
2,770✔
442
    if (simulation::n_realizations <= 1) {
2,770✔
443
      k_combined[1] = std::numeric_limits<double>::infinity();
1,900✔
444
    }
445
    return 0;
2,770✔
446
  }
447

448
  // Initialize variables
449
  int64_t n = simulation::n_realizations;
9,600✔
450

451
  // Copy estimates of k-effective and its variance (not variance of the mean)
452
  const auto& gt = simulation::global_tallies;
9,600✔
453

454
  array<double, 3> kv {};
9,600✔
455
  xt::xtensor<double, 2> cov = xt::zeros<double>({3, 3});
9,600✔
456
  kv[0] = gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n;
9,600✔
457
  kv[1] = gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n;
9,600✔
458
  kv[2] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n;
9,600✔
459
  cov(0, 0) =
19,200✔
460
    (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) /
9,600✔
461
    (n - 1);
9,600✔
462
  cov(1, 1) =
19,200✔
463
    (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n * kv[1] * kv[1]) /
9,600✔
464
    (n - 1);
9,600✔
465
  cov(2, 2) =
19,200✔
466
    (gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM_SQ) - n * kv[2] * kv[2]) /
9,600✔
467
    (n - 1);
9,600✔
468

469
  // Calculate covariances based on sums with Bessel's correction
470
  cov(0, 1) = (simulation::k_col_abs - n * kv[0] * kv[1]) / (n - 1);
9,600✔
471
  cov(0, 2) = (simulation::k_col_tra - n * kv[0] * kv[2]) / (n - 1);
9,600✔
472
  cov(1, 2) = (simulation::k_abs_tra - n * kv[1] * kv[2]) / (n - 1);
9,600✔
473
  cov(1, 0) = cov(0, 1);
9,600✔
474
  cov(2, 0) = cov(0, 2);
9,600✔
475
  cov(2, 1) = cov(1, 2);
9,600✔
476

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

484
  // First we will identify if there are any matching estimators
485
  int i, j;
486
  bool use_three = false;
9,600✔
487
  if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) &&
9,623✔
488
      (std::abs(cov(0, 0) - cov(1, 1)) / cov(0, 0) < FP_REL_PRECISION)) {
23✔
489
    // 0 and 1 match, so only use 0 and 2 in our comparisons
490
    i = 0;
20✔
491
    j = 2;
20✔
492

493
  } else if ((std::abs(kv[0] - kv[2]) / kv[0] < FP_REL_PRECISION) &&
9,580!
UNCOV
494
             (std::abs(cov(0, 0) - cov(2, 2)) / cov(0, 0) < FP_REL_PRECISION)) {
×
495
    // 0 and 2 match, so only use 0 and 1 in our comparisons
UNCOV
496
    i = 0;
×
UNCOV
497
    j = 1;
×
498

499
  } else if ((std::abs(kv[1] - kv[2]) / kv[1] < FP_REL_PRECISION) &&
9,582!
500
             (std::abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) {
2!
501
    // 1 and 2 match, so only use 0 and 1 in our comparisons
UNCOV
502
    i = 0;
×
UNCOV
503
    j = 1;
×
504

505
  } else {
506
    // No two estimators match, so set boolean to use all three estimators.
507
    use_three = true;
9,580✔
508
  }
509

510
  if (use_three) {
9,600✔
511
    // Use three estimators as derived in the paper by Urbatsch
512

513
    // Initialize variables
514
    double g = 0.0;
9,580✔
515
    array<double, 3> S {};
9,580✔
516

517
    for (int l = 0; l < 3; ++l) {
38,320✔
518
      // Permutations of estimates
519
      int k;
520
      switch (l) {
28,740!
521
      case 0:
9,580✔
522
        // i = collision, j = absorption, k = tracklength
523
        i = 0;
9,580✔
524
        j = 1;
9,580✔
525
        k = 2;
9,580✔
526
        break;
9,580✔
527
      case 1:
9,580✔
528
        // i = absortion, j = tracklength, k = collision
529
        i = 1;
9,580✔
530
        j = 2;
9,580✔
531
        k = 0;
9,580✔
532
        break;
9,580✔
533
      case 2:
9,580✔
534
        // i = tracklength, j = collision, k = absorption
535
        i = 2;
9,580✔
536
        j = 0;
9,580✔
537
        k = 1;
9,580✔
538
        break;
9,580✔
539
      }
540

541
      // Calculate weighting
542
      double f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) +
28,740✔
543
                 cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k));
28,740✔
544

545
      // Add to S sums for variance of combined estimate
546
      S[0] += f * cov(0, l);
28,740✔
547
      S[1] += (cov(j, j) + cov(k, k) - 2.0 * cov(j, k)) * kv[l] * kv[l];
28,740✔
548
      S[2] += (cov(k, k) + cov(i, j) - cov(j, k) - cov(i, k)) * kv[l] * kv[j];
28,740✔
549

550
      // Add to sum for combined k-effective
551
      k_combined[0] += f * kv[l];
28,740✔
552
      g += f;
28,740✔
553
    }
554

555
    // Complete calculations of S sums
556
    for (auto& S_i : S) {
38,320✔
557
      S_i *= (n - 1);
28,740✔
558
    }
559
    S[0] *= (n - 1) * (n - 1);
9,580✔
560

561
    // Calculate combined estimate of k-effective
562
    k_combined[0] /= g;
9,580✔
563

564
    // Calculate standard deviation of combined estimate
565
    g *= (n - 1) * (n - 1);
9,580✔
566
    k_combined[1] =
9,580✔
567
      std::sqrt(S[0] / (g * n * (n - 3)) * (1 + n * ((S[1] - 2 * S[2]) / g)));
9,580✔
568

569
  } else {
570
    // Use only two estimators
571
    // These equations are derived analogously to that done in the paper by
572
    // Urbatsch, but are simpler than for the three estimators case since the
573
    // block matrices of the three estimator equations reduces to scalars here
574

575
    // Store the commonly used term
576
    double f = kv[i] - kv[j];
20✔
577
    double g = cov(i, i) + cov(j, j) - 2.0 * cov(i, j);
20✔
578

579
    // Calculate combined estimate of k-effective
580
    k_combined[0] = kv[i] - (cov(i, i) - cov(i, j)) / g * f;
20✔
581

582
    // Calculate standard deviation of combined estimate
583
    k_combined[1] = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) *
20✔
584
                    (g + n * f * f) / (n * (n - 2) * g * g);
20✔
585
    k_combined[1] = std::sqrt(k_combined[1]);
20✔
586
  }
587
  return 0;
9,600✔
588
}
9,600✔
589

590
void shannon_entropy()
7,000✔
591
{
592
  // Get source weight in each mesh bin
593
  bool sites_outside;
594
  xt::xtensor<double, 1> p =
595
    simulation::entropy_mesh->count_sites(simulation::fission_bank.data(),
7,000✔
596
      simulation::fission_bank.size(), &sites_outside);
14,000✔
597

598
  // display warning message if there were sites outside entropy box
599
  if (sites_outside) {
7,000!
UNCOV
600
    if (mpi::master)
×
UNCOV
601
      warning("Fission source site(s) outside of entropy box.");
×
602
  }
603

604
  if (mpi::master) {
7,000✔
605
    // Normalize to total weight of bank sites
606
    p /= xt::sum(p);
6,950✔
607

608
    // Sum values to obtain Shannon entropy
609
    double H = 0.0;
6,950✔
610
    for (auto p_i : p) {
533,050✔
611
      if (p_i > 0.0) {
526,100✔
612
        H -= p_i * std::log2(p_i);
412,450✔
613
      }
614
    }
615

616
    // Add value to vector
617
    simulation::entropy.push_back(H);
6,950✔
618
  }
619
}
7,000✔
620

621
void ufs_count_sites()
150✔
622
{
623
  if (simulation::current_batch == 1 && simulation::current_gen == 1) {
150!
624
    // On the first generation, just assume that the source is already evenly
625
    // distributed so that effectively the production of fission sites is not
626
    // biased
627

628
    std::size_t n = simulation::ufs_mesh->n_bins();
15✔
629
    double vol_frac = simulation::ufs_mesh->volume_frac_;
15✔
630
    simulation::source_frac = xt::xtensor<double, 1>({n}, vol_frac);
15✔
631

632
  } else {
15✔
633
    // count number of source sites in each ufs mesh cell
634
    bool sites_outside;
635
    simulation::source_frac =
636
      simulation::ufs_mesh->count_sites(simulation::source_bank.data(),
270✔
637
        simulation::source_bank.size(), &sites_outside);
270✔
638

639
    // Check for sites outside of the mesh
640
    if (mpi::master && sites_outside) {
135!
UNCOV
641
      fatal_error("Source sites outside of the UFS mesh!");
×
642
    }
643

644
#ifdef OPENMC_MPI
645
    // Send source fraction to all processors
646
    int n_bins = simulation::ufs_mesh->n_bins();
90✔
647
    MPI_Bcast(
90✔
648
      simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
90✔
649
#endif
650

651
    // Normalize to total weight to get fraction of source in each cell
652
    double total = xt::sum(simulation::source_frac)();
135✔
653
    simulation::source_frac /= total;
135✔
654

655
    // Since the total starting weight is not equal to n_particles, we need to
656
    // renormalize the weight of the source sites
657
    for (int i = 0; i < simulation::work_per_rank; ++i) {
90,135✔
658
      simulation::source_bank[i].wgt *= settings::n_particles / total;
90,000✔
659
    }
660
  }
661
}
150✔
662

663
double ufs_get_weight(const Particle& p)
85,620✔
664
{
665
  // Determine indices on ufs mesh for current location
666
  int mesh_bin = simulation::ufs_mesh->get_bin(p.r());
85,620✔
667
  if (mesh_bin < 0) {
85,620!
UNCOV
668
    p.write_restart();
×
UNCOV
669
    fatal_error("Source site outside UFS mesh!");
×
670
  }
671

672
  if (simulation::source_frac(mesh_bin) != 0.0) {
85,620✔
673
    return simulation::ufs_mesh->volume_frac_ /
46,450✔
674
           simulation::source_frac(mesh_bin);
46,450✔
675
  } else {
676
    return 1.0;
39,170✔
677
  }
678
}
679

680
void write_eigenvalue_hdf5(hid_t group)
4,154✔
681
{
682
  write_dataset(group, "n_inactive", settings::n_inactive);
4,154✔
683
  write_dataset(group, "generations_per_batch", settings::gen_per_batch);
4,154✔
684
  write_dataset(group, "k_generation", simulation::k_generation);
4,154✔
685
  if (settings::entropy_on) {
4,154✔
686
    write_dataset(group, "entropy", simulation::entropy);
420✔
687
  }
688
  write_dataset(group, "k_col_abs", simulation::k_col_abs);
4,154✔
689
  write_dataset(group, "k_col_tra", simulation::k_col_tra);
4,154✔
690
  write_dataset(group, "k_abs_tra", simulation::k_abs_tra);
4,154✔
691
  array<double, 2> k_combined;
692
  openmc_get_keff(k_combined.data());
4,154✔
693
  write_dataset(group, "k_combined", k_combined);
4,154✔
694
}
4,154✔
695

696
void read_eigenvalue_hdf5(hid_t group)
60✔
697
{
698
  read_dataset(group, "generations_per_batch", settings::gen_per_batch);
60✔
699
  int n = simulation::restart_batch * settings::gen_per_batch;
60✔
700
  simulation::k_generation.resize(n);
60✔
701
  read_dataset(group, "k_generation", simulation::k_generation);
60✔
702
  if (settings::entropy_on) {
60✔
703
    read_dataset(group, "entropy", simulation::entropy);
10✔
704
  }
705
  read_dataset(group, "k_col_abs", simulation::k_col_abs);
60✔
706
  read_dataset(group, "k_col_tra", simulation::k_col_tra);
60✔
707
  read_dataset(group, "k_abs_tra", simulation::k_abs_tra);
60✔
708
}
60✔
709

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