Time dependent MC

[ChasingNeutrons]

This PR adds kinetic transport capabilities. This is heavily based on work done by @martin220199

I have implemented a kinetic physics package which tracks neutrons and precursors/delayed neutrons through time. This required adding a new particle type (precursors), tweaking transport operators to track time, modifying collision operators to produce delayed neutrons, implementing combing in the dungeon, and a map for time. I also fixed one or two silly things, like getting our shakes -> seconds conversion wrong! I have included a prompt only benchmark problem, Stacy 30, which agrees very well with Serpent results. I am a bit less certain on problems with delayed neutrons, though some very very preliminary runs suggest they agree.

The algorithm is the standard one for time-dependent without feedback: the outer loop is batches, the inner loop is time, i.e., we follow a batch of neutrons throughout all time and then repeat from the start.

It is still a bit primitive in several senses. The source of neutrons is arbitrary, but can only be introduced at time t = 0. There is no critical source sampling implemented yet. There is no grouped precursor sampling: this may affect variance significantly, but is somewhat complicated to implement due to the different spectra of different precursor groups (we would need to communicate the spectrum from the collision operator, where precursors are produced, to wherever particles are sampled subsequently). It is also not yet reproducible! Any obvious tips on that would be appreciated.

I would also appreciate any comments from @martin220199 who knows this stuff better than I do.

I would also appreciate any thoughts on future work, e.g., incorporating multiphysics. The Serpent approach seems reasonable: do one step at a time (with multiple batches), before communicating results, wiping tallies, and starting from the next step. This may require a different physics package.

[martin220199]

LGTM.

1. I agree the "sequential" Serpent (and TRIPOLI) approach seems sensible for incorporating statistically robust feedback, and this is what I had in mind when implementing kinetic to set up for dynamic. Do you have any thoughts on how to properly handle different simulation (for population control) and tally grid structures when communicating results and wiping tallies? This discussion is perhaps easier to have in person, but say we have simulation time grid t_s, and tally grid t_T. I currently only support t_s = n*t_T, for n >= 1, i.e., each simulation time step covers a multiple of tally time steps, and they are aligned to ensure statistics are computed properly batch-wise. Otherwise, if you have incomplete tally intervals for which you compute statistics, you would have to keep the accumulations of the first and second moments for those tallies, to then infer the total accumulations of the first and second order moments for the corresponding complete tally steps, to then compute mean and STD. So wiping becomes less trivial if we want to be more flexible, I think.

2. Is there a reason why sources cannot simply be introduced at any given input time t_i? That is, the user inputs t_i -> map to corresponding timestep -> call fixedSource % generate and particles are initialised with p % t = t_i? Moreover, for more complicated sources (distributed in time), I think that in the bufferloop, before calling the transport operator, you can check if the particle has fate == aged_fate, which shouldn't have been set unless it crossed a time boundary, so particles that are not born yet will not be simulated.

[ChasingNeutrons]

Thanks @martin220199 !

Honestly, I don't have any good ideas which aren't very complicated, constraining, or which don't radically change tallies. My view (again following Serpent) is that in the dynamicPhysicsPackage (or whatever it is called), it simulates one tally time step at a time. This may have several population control points, but essentially it won't really be tallying over time - it will produce a tally output for that step, which another solver will receive to update the system's physics. Then all tallies will be wiped before commencing the next step. I am not sure how valid higher statistical moments in such simulations are due to correlations from physics feedback anyway. I am hoping to at least have some means of addressing multiphysics in a few upcoming PRs, so maybe this will look a bit clearer then. Does that sort of answer things?

That is a good point on time-dependent sources - I hadn't thought of that approach! Very elegant. Perhaps the neatest way is to make a follow-up PR which add's time dependence in the source and adds these tweaks to the physicsPackage.

[martin220199]

Ok, that sounds good!

[valeriaRaffuzzi]

I am confused about some points, clarification would be good... and I might have missed some because it's a big one!

This will be a great feature to have :)

[valeriaRaffuzzi]

It might be more readable if you did the switch after the transport loop: during transport, store all the precursors being generated into nextPrecursors, and swap this and next at the end.

[ChasingNeutrons]

I think I need to talk through this as I'm not sure if the result is more readable. That means, at the start of the cycle, checking whether nextPrecursors is empty, for example, which is a bit confusing. I am probably not thinking it through correctly though.

[valeriaRaffuzzi]

I think that flipping the precursors should be done at the end of the timeStep, rather than half way through: during a step, you are saving precursors meant for the next step inside 'thisPrecursors', due to flipping before the step is over. This way, the notions of 'this' and 'next' don't mean anything anymore since you're using them arbitrarily. This is a matter of aesthetics and consistency.

[valeriaRaffuzzi]

I don't understand this: above (lines 209-212) you are forcing decay and adding pNew to the queue to simulate it. Here, you are forcing decay and simulating pNew, while keeping the precursor too for the next cycle?
pNew has an adjusted weight, but p itself is unchanged: it seems to me like this is duplicated.

[ChasingNeutrons]

This is the standard for forced precursor decay. Is the complaint that this is unclear? The operation is literally finding the particle weight at some time within the remaining interval corresponding to the precursor decaying then. The precursor will also decay in weight as time advances. Can you suggest any comments to make it clearer? I think most of the relevant description is maybe inside the forcedPrecursorDecay subroutine.

[valeriaRaffuzzi]

Great work, feel free to go ahead!