Picongpu: Current deposition accuracy study: Statistical analysis of simple test case

[steindev]

In reaction to issue #3377 I set up a simple test case to infer random and (possible) systematic inaccuracies in the current deposition in dependence of different methods used in a PIC simulation.

A parameter scan surveyed the dependence of the inaccuracies on the following configurations:

• startPosition: OnePosition, Random
• Particles per cell: 10, 1, 100
• Particle shape: PCS, TSC, CIC
• Current deposition method: EmZ, Esirkepov

TL;DR

• Inaccuracies are independent of the current deposition scheme: they are identical for the tested schemes (EmZ and Esirkepov)
• Inaccuracies are random: they are normal distributed with an expectation value of zero
• Inaccuracies are unexpectedly large: the distribution of deviations for a (supposedly) well-sampled setup shows that differences up to 12.1% appear at high frequency (2σ)
• It would be very informative if this test would be repeated with a different code @HighIander, @wmles, @n01r

Test setup

Current deposition is tested by initializing a homogeneous particle density resembling a square wire, where all particles move with the same velocity in the same direction such that they form a constant and homogeneous current density.

Particle (density) params

• square wire edge length = 1m in (x,y) direction, extends infinitely along z
• center of the wire's cross section is in the center of the simulation domain's (x,y) plane
• particle number density = 1.1960e+11/m^3
• initial particle velocity = 5.2187e+07m/s = 1.7408e-01c along z (chosen such that a particle crosses a cell in 10 time steps)
• particles start position within a cell: Random and OnePosition (all in lower left corner)
• Particles forming the current are positrons
• there are no other particles in the simulation

Grid params

• Δx = Δy = Δz = 7.8125e-03m
• Δt = 0.995 * Δx / c / sqrt(3.) = 1.4970e-11s

Particle & Field solver params

• Particle shape: CIC, TSC, PCS
• Particle pusher: Free
• homogenous macro-particle density of (10, 1, 100) macro-particles/cell
• minimum macro-particle weighting: 1
• particles are initially randomly distributed within a cell or all together at the lower left corner
• Current solver: EmZ, Esirkepov
• no current interpolation
• Field solver: None
• absorbing boundaries at simulation domain walls normal x-y plane
• periodic boundaries at remaining simulation domain walls parallel to x-y plane

Simulation params

• 3D simulation
• simulation domain length along x and y = 2m = 256Δx
• simulation domain length along z = 0.5m = 64Δz
• simulation time = 10Δt (but 1Δt would suffice)

Analysis and results overview

I study the particle number density and current density, provided by PIConGPU ADIOS1 sim_Data output, in a 20Δz x 100Δy x 100Δx sub volume which is cut out of the center of the total simulation volume.
Using a smaller sub-volume for the analysis ensures independence of density fluctuations at the wire surface due to the macro particle shape.
Within this sub volume, I regard the value of the number and current density at every cell as a measurement of the respective density and calculate their deviation from the expectation values.
Expectation values are the above given particle density and j = 1 A/m^2, respectively.
Thus there are 200000 measurements of deviation of the number and current density from their expectation value for each simulation / parameter set.

These measurements of deviation for the particle and current densities are statistically analyzed.
Before explaining the analysis, an example of such a measurement is shown.

Analysis method

Statistical analysis includes calculating several descriptive parameters of the distribution of deviations.
First, an estimator for the mean µ of the deviation of the measurements from the expectation is calculated.
If this value is significantly different from zero it will reveal a systematic error in our current deposition.

Second, the variance σ^2 of the distribution of deviations is calculated. This value characterizes the scattering of the deviations.
Its square root σ (i.e. the standard deviation of the distribution of deviations) is an estimator for the bandwidth in which we have to expect random inaccuracies within a simulation.
This value should be as close as possible to zero and values of σ on the order of 5% and larger are a sign for bad sampling, in my opinion.

Third, an "error Δµ of the mean value µ" is calculated from the variance σ^2.
It is an estimator for the width of an interval around the mean value µ in which we have to expect the mean if we repeat the simulation.
That is, Δµ is an estimator for the variance of µ.

Results

The parameters µ, Δµ, and σ for every parameter set in the parameter scan are given in the following table

| START POSITION | PPC | PARTICLE SHAPE | CURRENT SOLVER | mean of particle density deviation µ±Δµ [%] | standard deviation of particle density deviation σ [%] | mean of current density deviation µ±Δµ [%] | standard deviation of current density deviation σ [%] |
| -- | -- | -- | -- | -- | -- | -- | -- |
| OnePosition | 10u | PCS | EmZ | 9.9813e-05 ± 1.0530e-07| 0.00002 | 8.0344e-05 ± 1.0424e-07 | 0.00002 |
| OnePosition | 10u | PCS | Esirkepov | 9.9817e-05 ± 1.0500e-07| 0.00002 | 1.0494e-04 ± 9.1846e-08 | 0.00002 |
| OnePosition | 10u | TSC | EmZ | 6.7043e-05 ± 5.9915e-08| 0.00001 | 2.0729e-05 ± 1.0273e-07 | 0.00002 |
| OnePosition | 10u | TSC | Esirkepov | 6.7045e-05 ± 6.0002e-08| 0.00001 | 2.1023e-05 ± 1.0195e-07 | 0.00002 |
| OnePosition | 10u | CIC | EmZ | 4.7684e-05 ± 0.0000e+00| 0.00000 | 0.0000e+00 ± 0.0000e+00 | 0.00000 |
| OnePosition | 10u | CIC | Esirkepov | 4.7684e-05 ± 0.0000e+00| 0.00000 | 1.1921e-06 ± 1.5994e-08 | 0.00000 |
| OnePosition | 1u | PCS | EmZ | 4.7624e-05 ± 4.2410e-08| 0.00001 | 7.9638e-05 ± 3.4484e-08 | 0.00001 |
| OnePosition | 1u | PCS | Esirkepov | 4.7609e-05 ± 4.2403e-08| 0.00001 | 1.1694e-04 ± 3.3851e-08 | 0.00001 |
| OnePosition | 1u | TSC | EmZ | 3.1704e-05 ± 3.3207e-08| 0.00001 | 1.0430e-05 ± 2.9997e-08 | 0.00001 |
| OnePosition | 1u | TSC | Esirkepov | 3.1711e-05 ± 3.3212e-08| 0.00001 | 1.0298e-05 ± 2.7510e-08 | 0.00001 |
| OnePosition | 1u | CIC | EmZ | 3.5763e-05 ± 0.0000e+00| 0.00000 | 1.1921e-05 ± 0.0000e+00 | 0.00000 |
| OnePosition | 1u | CIC | Esirkepov | 3.5763e-05 ± 0.0000e+00| 0.00000 | 1.1921e-05 ± 0.0000e+00 | 0.00000 |
| OnePosition | 100u | PCS | EmZ | -6.3475e-04 ± 1.2782e-06| 0.00029 | 3.3827e-05 ± 6.4411e-07 | 0.00014 |
| OnePosition | 100u | PCS | Esirkepov | -6.3454e-04 ± 1.2780e-06| 0.00029 | 4.8274e-05 ± 6.4448e-07 | 0.00014 |
| OnePosition | 100u | TSC | EmZ | 6.3971e-04 ± 1.8364e-06| 0.00041 | -2.4019e-04 ± 1.5583e-06 | 0.00035 |
| OnePosition | 100u | TSC | Esirkepov | 6.3966e-04 ± 1.8364e-06| 0.00041 | -2.3996e-04 ± 1.5402e-06 | 0.00034 |
| OnePosition | 100u | CIC | EmZ | 3.5763e-05 ± 0.0000e+00| 0.00000 | 1.3113e-04 ± 0.0000e+00 | 0.00000 |
| OnePosition | 100u | CIC | Esirkepov | 3.5763e-05 ± 0.0000e+00| 0.00000 | 1.3113e-04 ± 0.0000e+00 | 0.00000 |
| Random | 10u | PCS | EmZ | 6.0102e-03 ± 2.6092e-02| 5.83427 | 1.3480e-03 ± 2.7120e-02 | 6.06419 |
| Random | 10u | PCS | Esirkepov | 6.0101e-03 ± 2.6092e-02| 5.83427 | 1.3476e-03 ± 2.7120e-02 | 6.06419 |
| Random | 10u | TSC | EmZ | 9.1019e-03 ± 3.7086e-02| 8.29270 | 1.2750e-03 ± 3.8993e-02 | 8.71920 |
| Random | 10u | TSC | Esirkepov | 9.1019e-03 ± 3.7086e-02| 8.29270 | 1.2724e-03 ± 3.8993e-02 | 8.71920 |
| Random | 10u | CIC | EmZ | 1.2382e-02 ± 5.8439e-02| 13.06741 | 6.0885e-03 ± 6.3792e-02 | 14.26432 |
| Random | 10u | CIC | Esirkepov | 1.2382e-02 ± 5.8439e-02| 13.06741 | 6.0885e-03 ± 6.3792e-02 | 14.26432 |
| Random | 1u | PCS | EmZ | -2.3932e-03 ± 8.2873e-02| 18.53093 | -7.8465e-04 ± 8.6052e-02 | 19.24189 |
| Random | 1u | PCS | Esirkepov | -2.3933e-03 ± 8.2873e-02| 18.53093 | -7.8549e-04 ± 8.6052e-02 | 19.24189 |
| Random | 1u | TSC | EmZ | 2.7939e-03 ± 1.1798e-01| 26.38209 | 2.1950e-03 ± 1.2403e-01 | 27.73409 |
| Random | 1u | TSC | Esirkepov | 2.7939e-03 ± 1.1798e-01| 26.38209 | 2.1923e-03 ± 1.2403e-01 | 27.73409 |
| Random | 1u | CIC | EmZ | 8.5195e-03 ± 1.8599e-01| 41.58897 | -5.6566e-04 ± 2.0287e-01 | 45.36219 |
| Random | 1u | CIC | Esirkepov | 8.5195e-03 ± 1.8599e-01| 41.58897 | -5.6566e-04 ± 2.0287e-01 | 45.36219 |
| Random | 100u | PCS | EmZ | 1.9071e-04 ± 8.2538e-03| 1.84560 | -7.2185e-04 ± 8.5799e-03 | 1.91853 |
| Random | 100u | PCS | Esirkepov | 1.9069e-04 ± 8.2538e-03| 1.84560 | -7.1633e-04 ± 8.5799e-03 | 1.91853 |
| Random | 100u | TSC | EmZ | 1.0060e-03 ± 1.1739e-02| 2.62483 | -5.9165e-04 ± 1.2343e-02 | 2.75995 |
| Random | 100u | TSC | Esirkepov | 1.0060e-03 ± 1.1739e-02| 2.62483 | -5.9375e-04 ± 1.2343e-02 | 2.75995 |
| Random | 100u | CIC | EmZ | 2.2539e-03 ± 1.8482e-02| 4.13281 | 5.3286e-04 ± 2.0175e-02 | 4.51124 |
| Random | 100u | CIC | Esirkepov | 2.2539e-03 ± 1.8482e-02| 4.13281 | 5.3286e-04 ± 2.0175e-02 | 4.51124 |

Discussion

• Current deposition by EmZ and Esirkepov produce identical results
• The mean of the density deviations is always ~0, which tells us that there are no systematic errors in the current deposition
• As expected and hoped for, the OnePosition start position initialization yields exact results as the particle sits directly in the lower left corner of the cell
• The width of the deviations distribution becomes smaller with higher order particle shapes and more particles
• For supposedly good sampling, i.e. PCS shape and 10 particles per cell, the full width at half maximum of the deviation distribution is (in my opinion) very large with a value of 14% for random initialization.
• I am especially surprised, that at better sampling, i.e. PCS shape and 100 particles per cell, the full width at half maximum still approaches 4.5% for random initialization. My (obviously ill-founded) feeling is, that it should become close to or below 1%.

I would love the see results of such a simulation from a different code in order to rate these results obtained with PIConGPU.
(And rebase my feelings :smile:)

cc'ing: @sbastrakov, @psychocoderHPC

[steindev]

@psychocoderHPC Since you were asking for it: The results do not change if FieldToParticleInterpolationNative and EsirkepovNative is used. See table below

| START POSITION | PPC | PARTICLE SHAPE| INTERPOLATION SCHEME | CURRENT SOLVER | Number dens µ±Δµ [%] | σ [%] | Current dens µ±Δµ [%] | σ [%] |
| -- | -- | -- | -- | -- | -- | -- | -- | -- |
| Random | 10u | PCS | FieldToParticleInterpolation | EsirkepovNative | -1.6126e-03 ± 2.6023e-02| 5.81897 | -1.1687e-03 ± 2.7130e-02 | 6.06647 |
| Random | 10u | PCS | FieldToParticleInterpolation | Esirkepov | -1.6126e-03 ± 2.6023e-02| 5.81897 | -1.1728e-03 ± 2.7130e-02 | 6.06647 |
| Random | 10u | PCS | FieldToParticleInterpolationNative | EsirkepovNative | -1.6126e-03 ± 2.6023e-02| 5.81897 | -1.1688e-03 ± 2.7130e-02 | 6.06647 |
| Random | 10u | PCS | FieldToParticleInterpolationNative | Esirkepov | -1.6126e-03 ± 2.6023e-02| 5.81897 | -1.1728e-03 ± 2.7130e-02 | 6.06647 |
| Random | 10u | TSC | FieldToParticleInterpolation | EsirkepovNative | -2.7108e-03 ± 3.6916e-02| 8.25476 | 9.3237e-04 ± 3.9008e-02 | 8.72245 |
| Random | 10u | TSC | FieldToParticleInterpolation | Esirkepov | -2.7108e-03 ± 3.6916e-02| 8.25476 | 9.2838e-04 ± 3.9008e-02 | 8.72245 |
| Random | 10u | TSC | FieldToParticleInterpolationNative | EsirkepovNative | -2.7108e-03 ± 3.6916e-02| 8.25476 | 9.3238e-04 ± 3.9008e-02 | 8.72245 |
| Random | 10u | TSC | FieldToParticleInterpolationNative | Esirkepov | -2.7108e-03 ± 3.6916e-02| 8.25476 | 9.2839e-04 ± 3.9008e-02 | 8.72245 |
| Random | 10u | CIC | FieldToParticleInterpolation | EsirkepovNative | -5.1717e-03 ± 5.7932e-02| 12.95401 | 9.7750e-03 ± 6.3843e-02 | 14.27572 |
| Random | 10u | CIC | FieldToParticleInterpolation | Esirkepov | -5.1717e-03 ± 5.7932e-02| 12.95401 | 9.7724e-03 ± 6.3843e-02 | 14.27572 |
| Random | 10u | CIC | FieldToParticleInterpolationNative | EsirkepovNative | -5.1717e-03 ± 5.7932e-02| 12.95401 | 9.7750e-03 ± 6.3843e-02 | 14.27572 |
| Random | 10u | CIC | FieldToParticleInterpolationNative | Esirkepov | -5.1717e-03 ± 5.7932e-02| 12.95401 | 9.7723e-03 ± 6.3843e-02 | 14.27572 |
| Random | 100u | PCS | FieldToParticleInterpolation | EsirkepovNative | -1.5216e-03 ± 8.2285e-03| 1.83994 | -1.0551e-03 ± 8.5819e-03 | 1.91896 |
| Random | 100u | PCS | FieldToParticleInterpolation | Esirkepov | -1.5216e-03 ± 8.2285e-03| 1.83994 | -1.0589e-03 ± 8.5819e-03 | 1.91896 |
| Random | 100u | PCS | FieldToParticleInterpolationNative | EsirkepovNative | -1.5216e-03 ± 8.2285e-03| 1.83994 | -1.0551e-03 ± 8.5819e-03 | 1.91896 |
| Random | 100u | PCS | FieldToParticleInterpolationNative | Esirkepov | -1.5216e-03 ± 8.2285e-03| 1.83994 | -1.0589e-03 ± 8.5819e-03 | 1.91896 |
| Random | 100u | TSC | FieldToParticleInterpolation | EsirkepovNative | -1.7148e-03 ± 1.1680e-02| 2.61175 | -4.8502e-04 ± 1.2345e-02 | 2.76043 |
| Random | 100u | TSC | FieldToParticleInterpolation | Esirkepov | -1.7148e-03 ± 1.1680e-02| 2.61175 | -4.8839e-04 ± 1.2345e-02 | 2.76043 |
| Random | 100u | TSC | FieldToParticleInterpolationNative | EsirkepovNative | -1.7148e-03 ± 1.1680e-02| 2.61175 | -4.8501e-04 ± 1.2345e-02 | 2.76043 |
| Random | 100u | TSC | FieldToParticleInterpolationNative | Esirkepov | -1.7148e-03 ± 1.1680e-02| 2.61175 | -4.8843e-04 ± 1.2345e-02 | 2.76043 |
| Random | 100u | CIC | FieldToParticleInterpolation | EsirkepovNative | -2.0976e-03 ± 1.8325e-02| 4.09755 | -2.9008e-04 ± 2.0214e-02 | 4.52006 |
| Random | 100u | CIC | FieldToParticleInterpolation | Esirkepov | -2.0976e-03 ± 1.8325e-02| 4.09755 | -2.9277e-04 ± 2.0214e-02 | 4.52006 |
| Random | 100u | CIC | FieldToParticleInterpolationNative | EsirkepovNative | -2.0976e-03 ± 1.8325e-02| 4.09755 | -2.9009e-04 ± 2.0214e-02 | 4.52006 |
| Random | 100u | CIC | FieldToParticleInterpolationNative | Esirkepov | -2.0977e-03 ± 1.8325e-02| 4.09755 | -2.9280e-04 ± 2.0214e-02 | 4.52006 |

[steindev]

To finish this study (and rebase my feelings) I wrote a little python script that evaluates particle density at grid nodes from a random initial particle distribution taking the particle shape into account.
Repeating particle distribution and density measurement 200k times resembles the above experiment in PIConGPU. I obtained the following values after statistical analysis

| PPC | PARTICLE SHAPE | Number dens µ±Δµ [%] | σ [%] |
| --- | ------ | -------------------- | ---------- |
| 1 | CIC | -0.0625 ± 1.8299e-01 | 40.91868 |
| 1 | TSC | -0.0210 ± 1.3432e-01 | 30.03420 |
| 1 | PQS | -0.1010 ± 1.0556e-01 | 23.60390 |
| 10 | CIC | -0.0115 ± 5.7723e-02 | 12.90736 |
| 10 | TSC | -0.0193 ± 4.2419e-02 | 9.48526 |
| 10 | PQS | 0.0185 ± 3.3433e-02 | 7.47594 |
| 100 | CIC | -0.0010 ± 1.8230e-02 | 4.07644 |
| 100 | TSC | -0.0018 ± 1.3425e-02 | 3.00196 |
| 100 | PQS | 0.0073 ± 1.0544e-02 | 2.35772 |

The error of the mean and the standard deviation are actually in quite good agreement with the results obtained from PIConGPU. :rocket:
As anticipated, only my initial expectations of what to obtain were wrong.
I will close the issue.

Big thanks to @sbastrakov for the idea!

Btw, I think the shapes of order m>=3 (PCS, P4S) are named incorrectly in PIConGPU.
See issue #3407.