Picongpu: dummy unbounding of the initially neutral ions removes the species

[hightower8083]

Faced a kinda wobbly behaviour in a case, if one starts with the initially neutral ions, but leaves in speciesInitialization.param the unbounding manupulator

ManipulateDerive<
            manipulators::binary::UnboundElectronsTimesWeighting

Since ions are set to neutral (particle[ boundElectrons_ ] = protonNumber), this manipulator should do nothing, but actually it somewhat removes all species. To make it more puzzling, this happens after 0-th raw output -- so one sees all the densities correctly, but next moment all particles just disappear.

This of course does not affect any simulation designed from the start, but in case if someone looks at the difference between pre-ionized and neutral cases (and is not careful enough to remove this manipulator) user can get really confused... took me some time to identify this issue, so I could fix it in PoGit.

Therefore, I don't suggest any actions on code (error message?), but it can be useful to keep it mentioned here.

[psychocoderHPC]

Could you please post the full Initialization pipeline using InitPipeline = ... from the file speciesInitialization.param

[hightower8083]

Hey Rene, sorry I should have been more specific.

The sortest way to reproduce the issue is to take a FoilLCT example, go to particle.param and in TwiceIonizedImpl replace

particle[ boundElectrons_ ] = protonNumber - 2._X;

by

particle[ boundElectrons_ ] = protonNumber - 0._X;

I look at species density fields in hdf5 plugin, and I see that in the initial record step=0, all the species are present as defined, but in the next output all species fields become 0 (laser is still here), and particle containers become empty.
Also the per-step computational time remains constant and laser goes as in vacuum, so I assume that all particles are really gone from the simulation (and not from plugin's range)

[sbastrakov]

After a discussion with @psychocoderHPC , we think the issue is as follows. Due to your change the UnboundElectronsTimesWeighting called here will create electrons with zero weights. Note that these are not all electrons: the previously created electrons here should be okay. Zero weights are currently not supported and they generate NaNs first in momentum due to how pushers are currently done, and then propagate towards currents and the EM field. You can see this in hdf5 output and particle output.

We are not yet fully sure if the particles are removed or "just" have NaN data.

This is of course an unfortunate and unexpected behaviour. Currently not sure how we can fix it or warn about it.

[hightower8083]

Thank you @sbastrakov for investigation. I had an impression that all species (not only derived ones) get corrupted and this happens after the initialisation -- at 0th iteration output looks fine, and next moment all's gone. But I don't have a data at hand for a quick check.

In PoGit I haven't yet implemented much of manipulators (e.g. Derive), but instead initialise species independently -- each has its own _initialPosition_, _densityProfile_, pusher, interpolator etc (hence its probably sub-optimal). I've patched this issue on my side, by suppressing UnboundElectronsTimesWeighting manipulator for initial_charge=0 species. Therefore, this behaviour does not really bother me, but it would be good to unravel, what is actually happening, since its kinda alarming.

BTW, I've noticed, that recent issues (#2943, #2944) affect functionality covered by PoGit. I'm keeping an eye on the component: user input issues, but it would be better to consider possible CI integration (mentioned by @ax3l in #2918)

[sbastrakov]

Yes, eventually all species get corrupted, as field values become NaN, which then sets momentum and position to NaN as well. It's the derived species that cause this behaviour in the first place.

[ax3l]

@hightower8083 the question is what you would like the manipulator UnboundElectronsTimesWeighting to do if operated on a neutral ion.

Here is its implementation: https://github.com/ComputationalRadiationPhysics/picongpu/blob/0.4.3/include/picongpu/particles/manipulators/binary/UnboundElectronsTimesWeighting.def#L72-L75

If it's called as linked above, falling silently back to "multiply by one" (aka no manipulation) might be unexpected, too... no?
We could set the multiMask to invalid and erase it xD

[hightower8083]

@ax3l I'm actually unsure of what should be a proper solution. In this particular case, ManipulateDerive needs to exit, rather than derive-without-manipulation. Probably a general flow should be producing the particles with zero weights and consequently removing them via MIN_WEIGHTING.
A more temporal patch could be to raise an error, in this particular case.

[ax3l]

Hm, so you as well suggest the two approaches (exit == "multiply-by-one") or erase ;-)
Multiply-by-one/do-nothing/exit is easiest to implement but removing might be more obvious.

[hightower8083]

@ax3l U've a bit lost me with how's "multiply-by-one"==exit ...
Since the whole unbounding operation for initially neutral ions is illegal, the proper solution seems to just abort the compilation and suggest user to remove the illegal manipulator. Does that make sense?

[sbastrakov]

I think the idea was that currently we first set the weighting of a derived macroparticle to be the same as for the original one, and then here multiply it by freeElectrons. So if before this line we add a conditional exit like if( freeElectrons == 0 ) return; it would effectively mean multiply-by-one.

About compile-time checks, not sure if we have data at that moment: while protonNumber is available at compile time, but I think boundElectrons is only accessible at run time.

[hightower8083]

OK, I think that my confusion comes from the fact that this operation does two things -- derives (creates) new particles and set proper weights. Since electrons should not be created in the first place, this cannot be fixed by tweaking the manipulator, except if there is a way to create and immediately remove these particles (e.g. by treating zero-weights) -- though its not a smartest thing to do anyway =)

So on my side, in pogit, I've made neutral ions to be a particular case with a truncated operator, which is not beautiful, but makes sense since its indeed a particular and frequent case.

On PIConGPU side, I think its important that users don't get mislead -- if we make it just skip multiplication, and someone will do this, he'll start with extra electrons in the box and will get false results. IMO, crushing the simulation with some message is a more healthy way.

[ax3l]

I see, you want to be able to specify a filter during ManipulateDerive, which tells if an individual particle shall be derived at all.

You can do this today via ManipulateDerive -> Manipulate (set multiMask of selected particles to invalid) -> FillAllGaps docs. The first two steps can potentially be one as well.

[hightower8083]

Hey thanx for an update! Indeed i see there are new options in particles class, and once I have time I could try improve my solution in pogit and enable more features.

However, for a moment I'm doing less picongpu cases and trying to get a hold of our cpu part of local cluster. Also looking at the pace at which code evolves, I've kinda put on hold the pogit support at least till next release (and if anyone will actually start using it).

[sbastrakov]

@hightower8083 thanks for keeping us updated on the pogit plans. I was actually a little reluctant to make any changes to the contents of .param files, so in #3022 went for a soft change.

[hightower8083]

@sbastrakov , yes I'm more or less keeping an eye on the evolution and i support a more clear moving window definition -- ideally a support for multiple entries with start, stop, speed arguments.

On the other hand backward compatibility is an inevitable burden to be carried, especially since your inputs are significantly entangled with the code structure. For this I'm still convinced that adding a more simple, user-oriented input processor (sorta PICMI) would relief this task letting devs to modify the code freely with only necessary adjustments to the input processor. But that's gonna take a compatibility-breaking release anyway, so guess its not for the near future..

[ax3l]

Don't worry about compatibility breaking, this is fine and can be documented in changelogs. Move fast and break things (but document it).