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

openmc-dev / openmc / 19929505944

04 Dec 2025 12:47PM UTC coverage: 82.063% (-0.005%) from 82.068%
19929505944

Pull #3659

github

web-flow
Merge 5d6b8bf87 into ad5a876be
Pull Request #3659: Delayed neutron fraction by isotope with IFP method

16971 of 23544 branches covered (72.08%)

Branch coverage included in aggregate %.

72 of 72 new or added lines in 11 files covered. (100.0%)

162 existing lines in 7 files now uncovered.

54950 of 64097 relevant lines covered (85.73%)

43444833.37 hits per line

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

90.28
/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()
85,873✔
51
{
52
  const auto& gt = simulation::global_tallies;
85,873✔
53

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

59
  double keff_reduced;
60
#ifdef OPENMC_MPI
61
  if (settings::solver_type != SolverType::RANDOM_RAY) {
48,335✔
62
    // Combine values across all processors
63
    MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE,
46,035✔
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,300✔
72
  }
73
#else
74
  keff_reduced = simulation::keff_generation;
37,538✔
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) {
85,873✔
80
    keff_reduced /= settings::n_particles;
82,193✔
81
  }
82

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

86
void synchronize_bank()
82,193✔
87
{
88
  simulation::time_bank.start();
82,193✔
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;
46,035✔
101
  int64_t n_bank = simulation::fission_bank.size();
46,035✔
102
  MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
46,035✔
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)
46,035✔
108
    start = 0;
30,100✔
109

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

114
#else
115
  int64_t start = 0;
36,158✔
116
  int64_t finish = simulation::fission_bank.size();
36,158✔
117
  int64_t total = finish;
36,158✔
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) {
82,193!
126
    fatal_error(
×
127
      "No fission sites banked on MPI rank " + std::to_string(mpi::rank));
×
128
  }
129

130
  simulation::time_bank_sample.start();
82,193✔
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;
82,193✔
135

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

138
  // Temporary banks for IFP
139
  vector<vector<int>> temp_delayed_groups;
82,193✔
140
  vector<vector<int>> temp_ancestors;
82,193✔
141
  vector<vector<double>> temp_lifetimes;
82,193✔
142
  if (settings::ifp_on) {
82,193✔
143
    resize_ifp_data(temp_delayed_groups, temp_ancestors, temp_lifetimes,
2,720✔
144
      3 * simulation::work_per_rank);
145
  }
146

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

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

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

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

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

166
  // Locally comb particles in fission_bank
167
  double tooth = tooth_start * teeth_distance + teeth_offset;
82,193✔
168
  for (int64_t i = tooth_start; i < tooth_end; i++) {
140,952,893✔
169
    int64_t idx = std::floor(tooth) - start;
140,870,700✔
170
    temp_sites[index_temp] = simulation::fission_bank[idx];
140,870,700✔
171
    if (settings::ifp_on) {
140,870,700✔
172
      copy_ifp_data_from_fission_banks(idx, temp_delayed_groups[index_temp],
2,420,000✔
173
        temp_ancestors[index_temp], temp_lifetimes[index_temp]);
2,420,000✔
174
    }
175
    ++index_temp;
140,870,700✔
176

177
    // Next tooth
178
    tooth += teeth_distance;
140,870,700✔
179
  }
180

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

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

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

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

204
  simulation::time_bank_sample.stop();
82,193✔
205
  simulation::time_bank_sendrecv.start();
82,193✔
206

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

211
  // IFP number of generation
212
  int ifp_n_generation;
213
  if (settings::ifp_on) {
46,035✔
214
    broadcast_ifp_n_generation(
1,400✔
215
      ifp_n_generation, temp_delayed_groups, temp_ancestors, temp_lifetimes);
216
  }
217

218
  int64_t index_local = 0;
46,035✔
219
  vector<MPI_Request> requests;
46,035✔
220

221
  // IFP send buffers
222
  vector<int> send_delayed_groups;
46,035✔
223
  vector<int> send_ancestors;
46,035✔
224
  vector<double> send_lifetimes;
46,035✔
225

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

232
    // Resize IFP send buffers
233
    if (settings::ifp_on && mpi::n_procs > 1) {
46,035✔
234
      resize_ifp_data(send_delayed_groups, send_ancestors, send_lifetimes,
600✔
235
        ifp_n_generation * 3 * simulation::work_per_rank);
600✔
236
    }
237

238
    while (start < finish) {
69,825✔
239
      // Determine the number of sites to send
240
      int64_t n =
241
        std::min(simulation::work_index[neighbor + 1], finish) - start;
61,729✔
242

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

251
        if (settings::ifp_on) {
15,694✔
252
          // Send IFP data
253
          send_ifp_info(index_local, n, ifp_n_generation, neighbor, requests,
295✔
254
            temp_delayed_groups, send_delayed_groups, temp_ancestors,
255
            send_ancestors, temp_lifetimes, send_lifetimes);
256
        }
257
      }
258

259
      // Increment all indices
260
      start += n;
61,729✔
261
      index_local += n;
61,729✔
262
      ++neighbor;
61,729✔
263

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

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

275
  start = simulation::work_index[mpi::rank];
46,035✔
276
  index_local = 0;
46,035✔
277

278
  // IFP receive buffers
279
  vector<int> recv_delayed_groups;
46,035✔
280
  vector<int> recv_ancestors;
46,035✔
281
  vector<double> recv_lifetimes;
46,035✔
282
  vector<DeserializationInfo> deserialization_info;
46,035✔
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]) {
46,035✔
289
    neighbor = mpi::n_procs - 1;
22,261✔
290
  } else {
291
    neighbor =
23,774✔
292
      upper_bound_index(bank_position, bank_position + mpi::n_procs, start);
23,774✔
293
  }
294

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

301
  while (start < simulation::work_index[mpi::rank + 1]) {
107,764✔
302
    // Determine how many sites need to be received
303
    int64_t n;
304
    if (neighbor == mpi::n_procs - 1) {
61,729✔
305
      n = simulation::work_index[mpi::rank + 1] - start;
37,955✔
306
    } else {
307
      n = std::min(bank_position[neighbor + 1],
23,774✔
308
            simulation::work_index[mpi::rank + 1]) -
47,548✔
309
          start;
310
    }
311

312
    if (neighbor != mpi::rank) {
61,729✔
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,694✔
317
      MPI_Irecv(&simulation::source_bank[index_local], static_cast<int>(n),
15,694✔
318
        mpi::source_site, neighbor, neighbor, mpi::intracomm, &requests.back());
15,694✔
319

320
      if (settings::ifp_on) {
15,694✔
321
        // Receive IFP data
322
        receive_ifp_data(index_local, n, ifp_n_generation, neighbor, requests,
295✔
323
          recv_delayed_groups, recv_ancestors, recv_lifetimes,
324
          deserialization_info);
325
      }
326

327
    } else {
328
      // If the source sites are on this processor, we can simply copy them
329
      // from the temp_sites bank
330

331
      index_temp = start - bank_position[mpi::rank];
46,035✔
332
      std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n],
46,035✔
333
        &simulation::source_bank[index_local]);
46,035✔
334

335
      if (settings::ifp_on) {
46,035✔
336
        copy_partial_ifp_data_to_source_banks(index_temp, n, index_local,
1,400✔
337
          temp_delayed_groups, temp_ancestors, temp_lifetimes);
338
      }
339
    }
340

341
    // Increment all indices
342
    start += n;
61,729✔
343
    index_local += n;
61,729✔
344
    ++neighbor;
61,729✔
345
  }
346

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

351
  int n_request = requests.size();
46,035✔
352
  MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
46,035✔
353

354
  if (settings::ifp_on) {
46,035✔
355
    deserialize_ifp_info(ifp_n_generation, deserialization_info,
1,400✔
356
      recv_delayed_groups, recv_ancestors, recv_lifetimes);
357
  }
358

359
#else
360
  std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles,
36,158✔
361
    simulation::source_bank.begin());
362
  if (settings::ifp_on) {
36,158✔
363
    copy_complete_ifp_data_to_source_banks(
1,320✔
364
      temp_delayed_groups, temp_ancestors, temp_lifetimes);
365
  }
366
#endif
367

368
  simulation::time_bank_sendrecv.stop();
82,193✔
369
  simulation::time_bank.stop();
82,193✔
370
}
128,228✔
371

372
void calculate_average_keff()
85,873✔
373
{
374
  // Determine overall generation and number of active generations
375
  int i = overall_generation() - 1;
85,873✔
376
  int n;
377
  if (simulation::current_batch > settings::n_inactive) {
85,873✔
378
    n = settings::gen_per_batch * simulation::n_realizations +
67,643✔
379
        simulation::current_gen;
380
  } else {
381
    n = 0;
18,230✔
382
  }
383

384
  if (n <= 0) {
85,873✔
385
    // For inactive generations, use current generation k as estimate for next
386
    // generation
387
    simulation::keff = simulation::k_generation[i];
18,230✔
388
  } else {
389
    // Sample mean of keff
390
    simulation::k_sum[0] += simulation::k_generation[i];
67,643✔
391
    simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
67,643✔
392

393
    // Determine mean
394
    simulation::keff = simulation::k_sum[0] / n;
67,643✔
395

396
    if (n > 1) {
67,643✔
397
      double t_value;
398
      if (settings::confidence_intervals) {
63,547✔
399
        // Calculate t-value for confidence intervals
400
        double alpha = 1.0 - CONFIDENCE_LEVEL;
112✔
401
        t_value = t_percentile(1.0 - alpha / 2.0, n - 1);
112✔
402
      } else {
403
        t_value = 1.0;
63,435✔
404
      }
405

406
      // Standard deviation of the sample mean of k
407
      simulation::keff_std =
63,547✔
408
        t_value *
63,547✔
409
        std::sqrt(
63,547✔
410
          (simulation::k_sum[1] / n - std::pow(simulation::keff, 2)) / (n - 1));
63,547✔
411

412
      // In some cases (such as an infinite medium problem), random ray
413
      // may estimate k exactly and in an unvarying manner between iterations.
414
      // In this case, the floating point roundoff between the division and the
415
      // power operations may cause an extremely small negative value to occur
416
      // inside the sqrt operation, leading to NaN. If this occurs, we check for
417
      // it and set the std dev to zero.
418
      if (!std::isfinite(simulation::keff_std)) {
63,547!
UNCOV
419
        simulation::keff_std = 0.0;
×
420
      }
421
    }
422
  }
423
}
85,873✔
424

425
int openmc_get_keff(double* k_combined)
12,659✔
426
{
427
  k_combined[0] = 0.0;
12,659✔
428
  k_combined[1] = 0.0;
12,659✔
429

430
  // Special case for n <=3. Notice that at the end,
431
  // there is a N-3 term in a denominator.
432
  if (simulation::n_realizations <= 3 ||
12,659✔
433
      settings::solver_type == SolverType::RANDOM_RAY) {
9,799✔
434
    k_combined[0] = simulation::keff;
3,036✔
435
    k_combined[1] = simulation::keff_std;
3,036✔
436
    if (simulation::n_realizations <= 1) {
3,036✔
437
      k_combined[1] = std::numeric_limits<double>::infinity();
2,090✔
438
    }
439
    return 0;
3,036✔
440
  }
441

442
  // Initialize variables
443
  int64_t n = simulation::n_realizations;
9,623✔
444

445
  // Copy estimates of k-effective and its variance (not variance of the mean)
446
  const auto& gt = simulation::global_tallies;
9,623✔
447

448
  array<double, 3> kv {};
9,623✔
449
  xt::xtensor<double, 2> cov = xt::zeros<double>({3, 3});
9,623✔
450
  kv[0] = gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n;
9,623✔
451
  kv[1] = gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n;
9,623✔
452
  kv[2] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n;
9,623✔
453
  cov(0, 0) =
19,246✔
454
    (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) /
9,623✔
455
    (n - 1);
9,623✔
456
  cov(1, 1) =
19,246✔
457
    (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n * kv[1] * kv[1]) /
9,623✔
458
    (n - 1);
9,623✔
459
  cov(2, 2) =
19,246✔
460
    (gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM_SQ) - n * kv[2] * kv[2]) /
9,623✔
461
    (n - 1);
9,623✔
462

463
  // Calculate covariances based on sums with Bessel's correction
464
  cov(0, 1) = (simulation::k_col_abs - n * kv[0] * kv[1]) / (n - 1);
9,623✔
465
  cov(0, 2) = (simulation::k_col_tra - n * kv[0] * kv[2]) / (n - 1);
9,623✔
466
  cov(1, 2) = (simulation::k_abs_tra - n * kv[1] * kv[2]) / (n - 1);
9,623✔
467
  cov(1, 0) = cov(0, 1);
9,623✔
468
  cov(2, 0) = cov(0, 2);
9,623✔
469
  cov(2, 1) = cov(1, 2);
9,623✔
470

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

478
  // First we will identify if there are any matching estimators
479
  int i, j;
480
  bool use_three = false;
9,623✔
481
  if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) &&
9,645✔
482
      (std::abs(cov(0, 0) - cov(1, 1)) / cov(0, 0) < FP_REL_PRECISION)) {
22!
483
    // 0 and 1 match, so only use 0 and 2 in our comparisons
484
    i = 0;
22✔
485
    j = 2;
22✔
486

487
  } else if ((std::abs(kv[0] - kv[2]) / kv[0] < FP_REL_PRECISION) &&
9,601!
UNCOV
488
             (std::abs(cov(0, 0) - cov(2, 2)) / cov(0, 0) < FP_REL_PRECISION)) {
×
489
    // 0 and 2 match, so only use 0 and 1 in our comparisons
UNCOV
490
    i = 0;
×
UNCOV
491
    j = 1;
×
492

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

499
  } else {
500
    // No two estimators match, so set boolean to use all three estimators.
501
    use_three = true;
9,601✔
502
  }
503

504
  if (use_three) {
9,623✔
505
    // Use three estimators as derived in the paper by Urbatsch
506

507
    // Initialize variables
508
    double g = 0.0;
9,601✔
509
    array<double, 3> S {};
9,601✔
510

511
    for (int l = 0; l < 3; ++l) {
38,404✔
512
      // Permutations of estimates
513
      int k;
514
      switch (l) {
28,803!
515
      case 0:
9,601✔
516
        // i = collision, j = absorption, k = tracklength
517
        i = 0;
9,601✔
518
        j = 1;
9,601✔
519
        k = 2;
9,601✔
520
        break;
9,601✔
521
      case 1:
9,601✔
522
        // i = absortion, j = tracklength, k = collision
523
        i = 1;
9,601✔
524
        j = 2;
9,601✔
525
        k = 0;
9,601✔
526
        break;
9,601✔
527
      case 2:
9,601✔
528
        // i = tracklength, j = collision, k = absorption
529
        i = 2;
9,601✔
530
        j = 0;
9,601✔
531
        k = 1;
9,601✔
532
        break;
9,601✔
533
      }
534

535
      // Calculate weighting
536
      double f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) +
28,803✔
537
                 cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k));
28,803✔
538

539
      // Add to S sums for variance of combined estimate
540
      S[0] += f * cov(0, l);
28,803✔
541
      S[1] += (cov(j, j) + cov(k, k) - 2.0 * cov(j, k)) * kv[l] * kv[l];
28,803✔
542
      S[2] += (cov(k, k) + cov(i, j) - cov(j, k) - cov(i, k)) * kv[l] * kv[j];
28,803✔
543

544
      // Add to sum for combined k-effective
545
      k_combined[0] += f * kv[l];
28,803✔
546
      g += f;
28,803✔
547
    }
548

549
    // Complete calculations of S sums
550
    for (auto& S_i : S) {
38,404✔
551
      S_i *= (n - 1);
28,803✔
552
    }
553
    S[0] *= (n - 1) * (n - 1);
9,601✔
554

555
    // Calculate combined estimate of k-effective
556
    k_combined[0] /= g;
9,601✔
557

558
    // Calculate standard deviation of combined estimate
559
    g *= (n - 1) * (n - 1);
9,601✔
560
    k_combined[1] =
9,601✔
561
      std::sqrt(S[0] / (g * n * (n - 3)) * (1 + n * ((S[1] - 2 * S[2]) / g)));
9,601✔
562

563
  } else {
564
    // Use only two estimators
565
    // These equations are derived analogously to that done in the paper by
566
    // Urbatsch, but are simpler than for the three estimators case since the
567
    // block matrices of the three estimator equations reduces to scalars here
568

569
    // Store the commonly used term
570
    double f = kv[i] - kv[j];
22✔
571
    double g = cov(i, i) + cov(j, j) - 2.0 * cov(i, j);
22✔
572

573
    // Calculate combined estimate of k-effective
574
    k_combined[0] = kv[i] - (cov(i, i) - cov(i, j)) / g * f;
22✔
575

576
    // Calculate standard deviation of combined estimate
577
    k_combined[1] = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) *
22✔
578
                    (g + n * f * f) / (n * (n - 2) * g * g);
22✔
579
    k_combined[1] = std::sqrt(k_combined[1]);
22✔
580
  }
581
  return 0;
9,623✔
582
}
9,623✔
583

584
void shannon_entropy()
7,695✔
585
{
586
  // Get source weight in each mesh bin
587
  bool sites_outside;
588
  xt::xtensor<double, 1> p =
589
    simulation::entropy_mesh->count_sites(simulation::fission_bank.data(),
7,695✔
590
      simulation::fission_bank.size(), &sites_outside);
15,390✔
591

592
  // display warning message if there were sites outside entropy box
593
  if (sites_outside) {
7,695!
UNCOV
594
    if (mpi::master)
×
UNCOV
595
      warning("Fission source site(s) outside of entropy box.");
×
596
  }
597

598
  if (mpi::master) {
7,695✔
599
    // Normalize to total weight of bank sites
600
    p /= xt::sum(p);
7,645✔
601

602
    // Sum values to obtain Shannon entropy
603
    double H = 0.0;
7,645✔
604
    for (auto p_i : p) {
586,355✔
605
      if (p_i > 0.0) {
578,710✔
606
        H -= p_i * std::log2(p_i);
453,695✔
607
      }
608
    }
609

610
    // Add value to vector
611
    simulation::entropy.push_back(H);
7,645✔
612
  }
613
}
7,695✔
614

615
void ufs_count_sites()
160✔
616
{
617
  if (simulation::current_batch == 1 && simulation::current_gen == 1) {
160!
618
    // On the first generation, just assume that the source is already evenly
619
    // distributed so that effectively the production of fission sites is not
620
    // biased
621

622
    std::size_t n = simulation::ufs_mesh->n_bins();
16✔
623
    double vol_frac = simulation::ufs_mesh->volume_frac_;
16✔
624
    simulation::source_frac = xt::xtensor<double, 1>({n}, vol_frac);
16✔
625

626
  } else {
16✔
627
    // count number of source sites in each ufs mesh cell
628
    bool sites_outside;
629
    simulation::source_frac =
630
      simulation::ufs_mesh->count_sites(simulation::source_bank.data(),
288✔
631
        simulation::source_bank.size(), &sites_outside);
288✔
632

633
    // Check for sites outside of the mesh
634
    if (mpi::master && sites_outside) {
144!
UNCOV
635
      fatal_error("Source sites outside of the UFS mesh!");
×
636
    }
637

638
#ifdef OPENMC_MPI
639
    // Send source fraction to all processors
640
    int n_bins = simulation::ufs_mesh->n_bins();
90✔
641
    MPI_Bcast(
90✔
642
      simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
90✔
643
#endif
644

645
    // Normalize to total weight to get fraction of source in each cell
646
    double total = xt::sum(simulation::source_frac)();
144✔
647
    simulation::source_frac /= total;
144✔
648

649
    // Since the total starting weight is not equal to n_particles, we need to
650
    // renormalize the weight of the source sites
651
    for (int i = 0; i < simulation::work_per_rank; ++i) {
99,144✔
652
      simulation::source_bank[i].wgt *= settings::n_particles / total;
99,000✔
653
    }
654
  }
655
}
160✔
656

657
double ufs_get_weight(const Particle& p)
94,182✔
658
{
659
  // Determine indices on ufs mesh for current location
660
  int mesh_bin = simulation::ufs_mesh->get_bin(p.r());
94,182✔
661
  if (mesh_bin < 0) {
94,182!
UNCOV
662
    p.write_restart();
×
UNCOV
663
    fatal_error("Source site outside UFS mesh!");
×
664
  }
665

666
  if (simulation::source_frac(mesh_bin) != 0.0) {
94,182✔
667
    return simulation::ufs_mesh->volume_frac_ /
51,095✔
668
           simulation::source_frac(mesh_bin);
51,095✔
669
  } else {
670
    return 1.0;
43,087✔
671
  }
672
}
673

674
void write_eigenvalue_hdf5(hid_t group)
4,048✔
675
{
676
  write_dataset(group, "n_inactive", settings::n_inactive);
4,048✔
677
  write_dataset(group, "generations_per_batch", settings::gen_per_batch);
4,048✔
678
  write_dataset(group, "k_generation", simulation::k_generation);
4,048✔
679
  if (settings::entropy_on) {
4,048✔
680
    write_dataset(group, "entropy", simulation::entropy);
451✔
681
  }
682
  write_dataset(group, "k_col_abs", simulation::k_col_abs);
4,048✔
683
  write_dataset(group, "k_col_tra", simulation::k_col_tra);
4,048✔
684
  write_dataset(group, "k_abs_tra", simulation::k_abs_tra);
4,048✔
685
  array<double, 2> k_combined;
686
  openmc_get_keff(k_combined.data());
4,048✔
687
  write_dataset(group, "k_combined", k_combined);
4,048✔
688
}
4,048✔
689

690
void read_eigenvalue_hdf5(hid_t group)
65✔
691
{
692
  read_dataset(group, "generations_per_batch", settings::gen_per_batch);
65✔
693
  int n = simulation::restart_batch * settings::gen_per_batch;
65✔
694
  simulation::k_generation.resize(n);
65✔
695
  read_dataset(group, "k_generation", simulation::k_generation);
65✔
696
  if (settings::entropy_on) {
65✔
697
    read_dataset(group, "entropy", simulation::entropy);
11✔
698
  }
699
  read_dataset(group, "k_col_abs", simulation::k_col_abs);
65✔
700
  read_dataset(group, "k_col_tra", simulation::k_col_tra);
65✔
701
  read_dataset(group, "k_abs_tra", simulation::k_abs_tra);
65✔
702
}
65✔
703

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

© 2025 Coveralls, Inc