Test continuity equation

[MaxThevenet]

Let us add a test checking divE-rho = 0 to machine precision for Esirkepov deposition. This could potentially be the source of discrepancy between single and double precision.

[EZoni]

I would be willing to help with this, if that is OK with everyone.

[MaxThevenet]

Sounds great, thanks! Note that this test could probably re-use existing example input files. Also, @mrowan137 has started looking into the discrepancy between single and double precision, for which this test should be relevant. Could you two coordinate to avoid duplicating the efforts?

[EZoni]

I think the effort of including a check on charge conservation is quite minor compared to the whole bug fixing, hence I think it probably makes more sense that both tasks are completed by @mrowan137, who is also investigating the bug. I expect that to be easier and faster than synchronizing separate efforts.

@mrowan137 @MaxThevenet Please, feel free to comment on this, if you do not agree.

On my side, I can provide some input on how I tested this before. If the goal is to assess how well charge is conserved, namely how accurately Gauss law is preserved, I would suggest the following:

1. add rho and divE to the output fields for the test case under study;
2. plot the absolute and relative errors of rho/epsilon_0 vs divE;
3. measure both errors with respect to some norm, for example the infinity norm (maximum in space).

I think that this can be done simply in the Python analysis script that is used to analyze the output data of the test case under study. It is done, for example, in the Python scripts Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py and Examples/Tests/Langmuir/analysis_langmuir_multi.py, where you find code blocks like:

# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E) when
# current correction (psatd.do_current_correction=1) is applied
if current_correction:
    rho  = data['rho' ].to_ndarray()
    divE = data['divE'].to_ndarray()
    Linf_norm = np.amax( np.abs( rho/epsilon_0 - divE ) ) / np.amax( np.abs( rho/epsilon_0 ) )
    print("error: " + str(Linf_norm))
    print("tolerance: 1.e-9")
    assert( Linf_norm < 1.e-9 )

I would suggest to add such checks for every existing test using Esirkepov deposition, rather than adding a new test whose only purpose is to check charge conservation. This could provide more robust safety checks for future development as well, as it would test one property (charge conservation) on a variety of test cases, rather than one property on one test case only. What do you think?

[MaxThevenet]

Yeah sorry I don't mean you should work together on the test, just decide together on who does it

[ax3l]

Not sure if helpful, but we usually plotted this in the past:

• https://github.com/ComputationalRadiationPhysics/picongpu/pull/1102#discussion_r39504877
• https://github.com/ComputationalRadiationPhysics/picongpu/pull/1582#issuecomment-251942671
• https://github.com/ComputationalRadiationPhysics/picongpu/blob/master/src/tools/bin/plot_chargeConservation.py
• https://github.com/ComputationalRadiationPhysics/picongpu/blob/master/src/tools/bin/plot_chargeConservation_overTime.py
• https://picongpu.readthedocs.io/en/0.4.3/usage/plugins/chargeConservation.html#usage-plugins-chargeconservation
• https://github.com/ComputationalRadiationPhysics/picongpu/pull/2275#issuecomment-331383352

For the first time steps one could see low-level spatial reduction groups at the machine precision level (in our case: supercell current deposition borders). We then tried to put those visual tests into some more quantitative curves.

[EZoni]

@MaxThevenet @mrowan137

Charge conservation with Esirkepov deposition seems well preserved for the following three test cases (input files followed by the infinity norm of the relative error of div(E)-rho/eps0, normalized with respect to max|rho/eps0|, at the last iteration):

input file	norm of error
inputs_2d_multi_rt.txt (similar to Examples/Tests/Langmuir/inputs_2d_multi_rt)	5.4604389e-10
inputs_3d_multi_rt.txt (similar to Examples/Tests/Langmuir/inputs_3d_multi_rt)	5.4511505e-10
inputs_2d.txt (similar to Examples/Modules/nci_corrector/inputs_2d)	2.2423738e-08

I'm currently looking at other test cases where results do not seem to be as good (e.g. tests in RZ geometry). Discrepancies are most likely due to diagnostics issues rather than to the deposition itself.

Note With the current version of the code div(E) at iteration i has to be compared with rho at iteration i+1, as the E and rho diagnostics are not synchronized in time when the FDTD solver is used. This is described in issue #1037. Automated unit tests to check charge conservation can be added only after this issue is solved.