Weight windows while using survival biasing

[teade]

Bug Description

Turning on survival biasing significantly reduces the effectiveness of weight windows. First image is the neutron flux across a block with a weight window applied and survival biasing turned off. The second image is the neutron flux across the same block with a weight window applied and survival biasing turned on.

Steps to Reproduce

Having looked at the source code this appears to be caused by this if statement in the physics.cpp file:

  // Play russian roulette if survival biasing is turned on
  if (settings::survival_biasing) {
    if (p.wgt() < settings::weight_cutoff) {
      russian_roulette(p, settings::weight_survive);
    }
  }

By default the settings::weight_cutoff is set to 0.25, so any particle which is split by the weight window to be below this value is rouletted. A work around a user can employ is to reduce the settings::weight_cutoff value. However, there is also a weight_window.cutoff value which limits the particles weight in the weight window routines and it would be better from a users point of view not to have two weight cutoff values which need to be set for different variance reduction games (and make weight windows less effective if not set properly).

Environment

Linux, OpenMC Version 0.14.0-dev

[yrrepy]

It is standard to disable survival biasing when using weight windows.
The two VR schemes clash.

https://link.springer.com/chapter/10.1007/978-3-031-04129-7_4
Section 4.1.2.3 Iteration 2
end of paragraph 1

[eepeterson]

It should definitely be feasible to apply the weight cutoff relative to the local weight window lower bound to get the benefits of both, no?

[teade]

It is standard to disable survival biasing when using weight windows. The two VR schemes clash.

https://link.springer.com/chapter/10.1007/978-3-031-04129-7_4 Section 4.1.2.3 Iteration 2 end of paragraph 1

I regularly run with both survival biasing/implicit capture and weight windows in MCNP. It seems to work OK. Does anyone else have any experience of this?

Trying to understand what is written in the book you link to; it appears the author is saying that particles which reach a part of the mesh without a weight window value set are rouletted quickly as they will be below the weight cutoff on the cut:n card. However, in the MCNP source code there appears to be an if statement in the colidn routine which means that if a WWP card is specified it skips the survival bias rouletting. I haven't played around with this too much, so there might be further conditions which mean the survival bias rouletting is activated if the particle leaves the weight window mesh.

However, in general, if your weight window goes across the whole area you are interested in, I don't see why you can't have survival biasing and ww on at the same time and just leave the rouletting to the weight window routine (and skip the roulette function in physics.cpp).

[yrrepy]

I've used weight windows extensively in MCNP, always turning survival biasing off. Via the cut or phys card.
Not sure about the source itself having further disabling mechanisms (wrt WWP)

For your problem you may be seeing a difference between the two and greater degradation due to different survival biasing parameters.
MCNP defaults to 0.25, 0.5.
OpenMC defaults to 0.25, 1.0

[pshriwise]

Thanks for noting this @teade. @eepeterson noticed some odd interactivity between these two forms of variance reduction.

I'll break the survival biasing method into two processes:

• Reduced particle weight by $1 - \frac{\sigma_a}{\sigma_t}$ when absorption events occur
• Rouletting of particles below the survival biasing weight cutoff, $w_c$

As you noted above, particle weights are likely to be much lower than $w_c$ and as a result of the survival biasing roulette occuring for every particle reaction and terminating particles with probability $1 - \frac{w_c}{w_s}$, it's highly likely that particles will be terminated earlier than is desirable when combining these two variance reduction techniques. Keeping the absorption reactions from terminating particles is generally desirable in problems that require variance reduction, however.

Ideally, we would be able to keep the survival biasing without the collision event rouletting that comes with it. This can be done by enabling survival biasing but setting the weight cutoff for rouletting to a very low value or zero. I'm inclined to say that OpenMC should do this automatically when weight windows are present and leave the rouletting to the weight windows, but there are particles that move through a weight window mesh region and exit to consider. In this case, they would likely continue on without being terminated until reaching a vacuum boundary, re-entering the weight window region, or reaching the minimum energy threshold.

In summary, when weight window mesh(es) are present and survival biasing is enabled my proposed options are:

1. Display a warning about both being enabled and the weight cutoff parameter
2. Do 1 and also set the weight cutoff parameter to zero (this change would be included in the warning)

I'm leaning toward 1. Option 2 is cleaner in some ways, but less flexible. What might be nice is to provide some guidance on how to set the $w_c$ parameter to avoid interference with the weight windows variance reduction based on weight window lower bounds perhaps.

Thoughts?

[teade]

Thanks @pshriwise; that is a good summary. However, I am not sure just displaying a warning will be sufficient; users don't always pay much attention to warnings. I would propose a slight alteration to 2 might be better.

As you point out the tricky bit is if the particle leaves the weight window without a weight cutoff it may lead to long run times. Which is not always ideal. I wonder if, while the particle is in the weight window mesh, $w_c$ could be set to zero, then a suitable $w_c$ value could be set on the fly if/when a particle leaves the weight window. Using @eepeterson suggestion, maybe this could be made relative to the lower weight window bound ($w_l$) of the last voxel the particle was in i.e. $w_c \times w_l$ ?

[pshriwise]

I think as a short-term measure a warning about this would be useful, no? I'm primarily for that b/c survival biasing isn't enabled by default.

Hrmmm I like it @teade, though this would require we change the $w_c$ parameter from a global value to a member of the Particle class as the value would be associated each particle as it leaves the mesh. Not impossible but more complicated than our current setup for survival biasing and I think that data attached to the particle would go unused more often than not while increasing the size of the object (something we try to avoid when we can). Another idea I'll pitch is setting $w_c$ to the lowest value of the weight window lower bounds for the problem would be effective "enough"?

Either way, I'd like to setup a contrived problem to understand how long these histories can be when they leave the mesh w/ survivial biasing on and $w_c = 0$. I don't have a good understanding of that right now.

[pshriwise]

@paulromano and I ended up returning to this issue today.

After a lot of discussion, the proposition is that that problems using both survival biasing and weight window variance reduction at the same time will rely on weight windows to roulette particles (no global rouletting) using the weight window lower bounds only while in the mesh.

If a particle exits the weight window mesh, the global roulette based on the survival biasing parameters will be applied. The assumption here is that if a particle exits the weight window mesh, it isn't going to contribute regions of interest within the mesh. I'll note that if the survival_normalization option is turned on it will still apply.

[pshriwise]

Thankfully the places where weight adjustments are made to particles are already quite close together in the code, making the introduction of a more generalized function for weight adjustment relatively easy to implement

Global rouletting occurs at the end of sample_neutron_reaction after which we immediately call apply_weight_windows

openmc/src/physics.cpp

Lines 49 to 66 in 7382b5d

	switch (p.type()) {
	case ParticleType::neutron:
	sample_neutron_reaction(p);
	break;
	case ParticleType::photon:
	sample_photon_reaction(p);
	break;
	case ParticleType::electron:
	sample_electron_reaction(p);
	break;
	case ParticleType::positron:
	sample_positron_reaction(p);
	break;
	}
	if (settings::weight_window_checkpoint_collision)
	apply_weight_windows(p);

[GuySten]

I've opened PR #3751 to try to fix this bug.