Picongpu: When generating particles, evaluate density at centers of cells ?

[sbastrakov]

This is a side note found while looking at #3377. It is definitely not that reason of it, but can also cause some issues.

Currently inside CreateDensity we eventually go over each cell and estimate density at cell corners to compute the number of real particles in the corresponding cell (and then apply our logic to convert to number of macroparticles, but this is irrelevant for this topic). The fact that we estimate density at cell corners to compute the number of real particles seems to have the following two downsides:

• Implementors of density functors with cell-aligned sharp edges have to take care of rounding off (without necessarily knowing how to, unless they know exactly our scheme!). Which is not trivial, and if we estimated at centers more naive implementations would have worked as their authors expected. Of course, for smooth densities like Gaussian this is not a real issue.
• When computing number of real particles in a cell as a product of cell volume and density at a point. The exact number is integral of density over the cell. With the current way we approximate it by a side-rectangle rule, since the density is estimated at a corner. Estimating at center would also be a rectangle rule and same computational complexity, but would give a better central-rectangle rule.

[sbastrakov]

IF we decide to change, however we would need to manually go over all density profiles and see if they implicitly assumed the currently used corner evaluation scheme.

[steindev]

I agree with sbastrakov in both points and support changing the current method to determining the density for particle generation at the centers of cells.

[bussmann]

IF we decide to change, however we would need to manually go over all density profiles and see if they implicitly assumed the currently used corner evaluation scheme.

These assumptions must be made explicit!

[pordyna]

Does CreateDensity use picongpu::traits::FieldPosition? I think we should avoid the situation when the density in the output is different form the one generated by CreateDensity. So we may need to consider changing FieldPosition to the center of the cell, for FieldTmp, as well. Though this would affect all scalar fields.

[steindev]

There is a subtle difference in meaning here. During a simulation and when writing output, the position were density will be computed from the distribution of the macro particles will still be the cell corner. And this is the place were it needs to be in order to be self-consistent with our definition of the Yee cell, i.e.with the positions of the E- and B-field components.

But the question here just regards macro particle generation at the beginning of a simulation (or during sliding into new volumes during a moving window simulation). The better approximation for the number of macro particles to be generated in a cell from a given density distribution is obtained by evaluating the given density at the cell center, not at the cell corner, as explained by @sbastrakov in point two of its initial post before.

However, I agreee with @pordyna that this may lead to confusion, as the density evaluated by PIConGPU for macro particle initialization may be different from the density field output. More precisely, this difference is to be expected when initializing with startPosition::OnePosition (and startPosition::QuietStart?) if the starting position is the cell corner.
But I think we can easily remedy this issue by setting the default start position of a macro particle to the cell center instead of the cell corner.

[sbastrakov]

Sure, regardless of the scheme being used the functors that depend on it have to document it. My point was that they seem to be really old code and do not do it currently. So that needs to be fixed regardless of where do we evaluate density.

[sbastrakov]

I think it's even more complicated than what @steindev explained. As CreateDensity is eventually used (together with other things like min weighting and desired number of macroparticles per cell) to determine the number of macroparticles per cell. However, it currently does not influence the distribution of macroparticle centers inside the cell. While for the output we have different options: simply count to the cell where the center is, or make weighted contributions according to particles' assignment function values. Also, of course the output is openPMD annotated so at least it self-describes the spatial positions of data points.

[sbastrakov]

So upon a further look the changes required are actually rather small. #3414 extends comments of what's the requirements for density functors and what's the meaning of their return value. Namely, it is a single value to approximate density in the given cell. Then naturally the best constant-value-approximation is the density value at cell center. However, implementors are free to choose what they want as long as it corresponds to the cell. Then #3415 only changes the evaluation point for user-provided functors with FreeFormulaImpl, as these functors are evaluatable at any point already. Then we can go through existing profiles and see where a half-cell shift is appropriate. But at worst they are just slightly less accurate then possible, since density should not have a subcell level fluctuations anyways.

[sbastrakov]

The change was done for the free functor density. Other density profiles may be updated on a case by case basis if necessary.