Combing for fission site sampling, and delayed neutron emission time

[ilhamv]

Description

Currently, OpenMC uses uniform Splitting-Roulette with fix-up for fission site sampling in synchronize_bank(). The fix-up is needed to exactly produce the targeted sample size; however, the fix-up introduces bias (even though negligible if the particle size is sufficiently large). This PR replaces the current fission site sampling with uniform Combing, which is unbiased and exactly produces the targeted sample size. The implementation maintains the reproducibility and scalability of synchronize_bank().

This PR changes the simulation random number stream, so a revision of regression test answer keys would be needed. The plan is to group this PR along with other PRs that also change the simulation random number stream: #2898.

cc: @paulromano

Splitting-Roulette and the fix-up

For k-eigenvalue simulation, OpenMC currently uses the uniform Splitting-Roulette technique followed by a fix-up to ensure exactly settings::n_particles are sampled.

The uniform Splitting-Roulette:

openmc/src/eigenvalue.cpp

Lines 158 to 177 in cfe210d

	for (int64_t i = 0; i < simulation::fission_bank.size(); i++) {
	const auto& site = simulation::fission_bank[i];
	// If there are less than n_particles particles banked, automatically add
	// int(n_particles/total) sites to temp_sites. For example, if you need
	// 1000 and 300 were banked, this would add 3 source sites per banked site
	// and the remaining 100 would be randomly sampled.
	if (total < settings::n_particles) {
	for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
	temp_sites[index_temp] = site;
	++index_temp;
	}
	}
	// Randomly sample sites needed
	if (prn(&seed) < p_sample) {
	temp_sites[index_temp] = site;
	++index_temp;
	}
	}

The fix-up:

openmc/src/eigenvalue.cpp

Lines 200 to 223 in cfe210d

	// Now that the sampling is complete, we need to ensure that we have exactly
	// n_particles source sites. The way this is done in a reproducible manner is
	// to adjust only the source sites on the last processor.
	if (mpi::rank == mpi::n_procs - 1) {
	if (finish > settings::n_particles) {
	// If we have extra sites sampled, we will simply discard the extra
	// ones on the last processor
	index_temp = settings::n_particles - start;
	} else if (finish < settings::n_particles) {
	// If we have too few sites, repeat sites from the very end of the
	// fission bank
	sites_needed = settings::n_particles - finish;
	for (int i = 0; i < sites_needed; ++i) {
	int i_bank = simulation::fission_bank.size() - sites_needed + i;
	temp_sites[index_temp] = simulation::fission_bank[i_bank];
	++index_temp;
	}
	}
	// the last processor should not be sending sites to right
	finish = simulation::work_index[mpi::rank + 1];
	}

The uniform Splitting-Roulette is unbiased and introduces minimum variance to the sampling process [Ref. 1]. However, Splitting-Roulette does not exactly produce the targeted population size. Since OpenMC requires each generation to run exactly with settings::n_particles, the fix-up is performed. This fix-up, however, introduces bias, even though negligible if the particle size is sufficiently large.

(Note that the fix-up is unnecessary if OpenMC allows changing particle sizes in a simulation. Also, it is worth noting that the population control itself is not strictly necessary, but it helps ensure that the particle population does not wildly vary in the early generations, especially if we have a poor initial guess.)

Alternative techniques

Two other population control techniques, Combing and Duplicate-Discard [Ref 1], are also unbiased and introduce minimum variance in the sampling process. However, unlike Splitting-Roulette, they produce exactly the targeted size. Among the two, Combing is simpler to implement and faster to perform (see [Ref 2] for a study on their implementations in OpenMC).

Combing may suffer from unwanted behavior due to possible correlation in the particle order in the bank. However, other than in contrived demonstrations, such unwanted behavior has not been reported in practice in multiple implementations of Combing in time-dependent Monte Carlo applications. In static simulations like k-eigenvalue calculation, the unwanted behavior becomes less relevant (it may increase the tally variance, but such evidence is lacking, as shown in [Ref 2]).

Test

The OPR benchmark was run with 1M particles in 200/500 inactive/active batches on 8 MPI processes, each with 72 threads.

With the current Splitting-Roulette with fix-up:
k = 1.07473 +- 4.18987e-05
Total runtime = 277.428 s
Site sampling runtime = 16.5646 s

With the proposed Combing:
k = 1.07466 +- 4.1204e-05
Total runtime = 275.814 s
Site sampling runtime = 17.1973 s

For a more extensive test, see [Ref 2].

Future plan

This fission site sampling is essentially population control [Ref 1], which can also be used in time-censussed particles in time-dependent simulations. The plan is to make the procedure a function so that it can be used in both eigenvalue and time-dependent fixed source simulations.

If it is desirable to support different population control techniques that the user can choose from, we can create a population control class like the one in https://github.com/ilhamv/openmc/tree/population_control.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[paulromano]

@ilhamv This update looks great to me; thanks a lot for putting this together! I'm current in the process of merging #2898 into this one and then updating test results so that we only have to do that once. I'm also going to run a series of criticality benchmarks to make sure that we are getting results that are consistent with the current develop branch and will report the results here.

For anyone who is interested in this topic, I highly highly recommend reading [Ref. 1] that is linked above, which has got to be one of the most thorough and well-written articles I've seen in a while! Kudos @ilhamv 😄

[paulromano]

Alright, I've merged #2898 into this branch and updated all the test results. Hopefully they pass in CI but we'll see how it goes.

Any time we do a wholesale update of regression test results like this, I also run a series of ICSBEP benchmarks to make sure that we haven't fundamentally messed up the physics behavior. I've completed that and happy to share the below reuslts, showing the difference in k-eff between the current develop branch and this PR. As you can see, the differences are pretty much random, within uncertainty, and average out to zero—all of which is expected ✔️

[ilhamv]

Thanks, Paul! 😄 🙏🏽

It's great to see that the use of combing passes/is verified through the ICSBEP validation tests. This well supports our claim that the "fear of the correlation effect" in combing is not relevant in static, k-eigenvalue simulations.

===

I'm not familiar yet with the CIs. But it appears that many of them failed..?

[paulromano]

There's just one pesky test that is failing in CI (deplete_with_transfer_rates) that I'm trying to track down. Once that's resolved we're ready to merge!