Bug w/ Vay Deposition

[EZoni]

If we run WarpX (development branch) with this input file and look at the raw Jz data with guard cells (raw field jz_fp), we get

which shows some spurious signal at the low and high end of the domain along Z (note that the plot includes guard cells, for Z < 0 and Z > 128).

However, if we comment out the call to WarpXSumGuardCells inside WarpX::ApplyFilterandSumBoundaryJ (note that filtering is OFF in this test),

https://github.com/ECP-WarpX/WarpX/blob/5959723aa414d812b37151deffcba9cd05002b07/Source/Parallelization/WarpXComm.cpp#L992

and rerun the exact same test, we get

where the spurious signal is not there anymore.

@jlvay and I were thinking that @WeiqunZhang and/or @atmyers could be able to help us debug this, since it might be related to something happening in the guard cells. Will probably follow up on Slack and/or offline.

Note that the test uses algo.current_deposition = vay, while the issue does not occur with other current deposition schemes, such as algo.current_deposition = direct.

[NeilZaim]

Hi @EZoni. I'm not sure if this is related to this issue but it looks similar so I'm putting it here:

I tried to run the CI tests again in PR #2336, which reduces the number of guard cells, and the Vay deposition test benchmarks changed. I investigated the issue and the difference comes from the fact that with my PR ng_depos_J is decreased from 4 to 3 in the WarpXSumGuardCells call that you mentioned above.

So I wondered what is the correct value that should be put for ng_depos_J and I'm not sure what is the answer. Looking at the code, it looks like the J fields that we are summing are coming from an inverse Fourier transform, so in general they should be nonzero everywhere in the guard cells. So I thought that it might make sense to use j[idim]->nGrowVect() (i.e. the total number of allocated guard cells for the J field) instead of ng_depos_J in the WarpXSumGuardCells call in order to reduce the truncation errors. I tried to do that and indeed it changed the results. However, when I did that the charge conservation assert of the 2D Vay deposition CI test (see below) failed, although it was succeeding when using ng_depos_J = 3 or 4, so maybe using the total number of allocated guard cells in the current sum does not produce the best results physically.

https://github.com/ECP-WarpX/WarpX/blob/2191ca66c11ff71c62015234b7024fb85d36d2df/Examples/Tests/VayDeposition/analysis.py#L32-L37

Do you know what value of ng_depos_J is correct from the algorithm point of view here?

[EZoni]

@NeilZaim Thank you for reporting that. I tried running the CI tests with Vay deposition on the branch of #2336, so on commit 2903322 but with the following code changes,

--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -986,10 +986,10 @@ void WarpX::ApplyFilterandSumBoundaryJ (
             ng_depos_J.min(ng);
             MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), j[idim]->nComp(), ng);
             bilinear_filter.ApplyStencil(jf, *j[idim], lev);
-            WarpXSumGuardCells(*(j[idim]), jf, period, ng_depos_J, 0, (j[idim])->nComp());
+            WarpXSumGuardCells(*(j[idim]), jf, period, ng, 0, (j[idim])->nComp());
         } else {
             ng_depos_J.min(ng);
-            WarpXSumGuardCells(*(j[idim]), period, ng_depos_J, 0, (j[idim])->nComp());
+            WarpXSumGuardCells(*(j[idim]), period, ng, 0, (j[idim])->nComp());
         }
     }
 }

namely with ng instead of ng_depos_J (I think ng could be better than j[idim]->nGrowVect(), because it is initially assigned as j[idim]->nGrowVect() but then also modified to take into account the stencil of the bilinear filter and things like that).

I do see that the checksum values change, which I think is expected, but I do not see that the checks on charge conservation fail. I ran the tests Langmuir_multi_2d_psatd_Vay_deposition, Langmuir_multi_2d_psatd_Vay_deposition_nodal, Langmuir_multi_psatd_Vay_deposition and Langmuir_multi_psatd_Vay_deposition_nodal, and the checks on charge conservation seemed to pass in all cases. Could you double check that you still see a failure in the checks on charge conservation instead?

Also, note that the change from ng to ng_depos_J.min(ng) had been introduced by myself in #2323, so I guess we thought at the time that summing only a few of the guard cells would be enough from an algorithmic point of view. But I agree with you that those fields come from an inverse FFT and do supposedly have non-zero values in all the guard cells.

[EZoni]

Regarding why WarpXSumGuardCells is called using ng_depos_J.min(ng) rather than ng, the reason is that SyncCurrent is usually called right after the current deposition (where only at most ng_depos_J guard cells are written), before any other operation in Fourier space.

However, after #3012 this is not the case anymore for Vay deposition: now SyncCurrent is called after the FFT of D and the inverse FFT of J (computed from D in Fourier space), as you pointed out. Hence, I think I agree that now, for the case of Vay deposition, it would make more sense to call WarpXSumGuardCells with ng instead of ng_depos_J.min(ng).

However, this does not fix the issue reported here, per se, it just moves a bit the longitudinal position of the spurious signal mentioned in the description of the issue.

[EZoni]

Suggestion by @jlvay: the bilinear filter needs to be applied to D, not J, while the sum over guard cells needs to be applied to J (after inversion in Fourier space). This can be tried locally, by adding some flag that triggers either filtering or sum over guard cells in SyncCurrent.

[NeilZaim]

Thanks for your replies @EZoni!

I have checked again with your modifications and I still see the charge conservation error with the VayDeposition2D test. Running ./run_test.sh VayDeposition2D locally I get:

working on test: VayDeposition2D
   re-making clean...
   building...

configuring VayDeposition2D build...
   mkdir /tmp/ci-OIrlI1zRvD/warpx/builddir
   cmake  -DAMReX_ASSERTIONS=ON -DAMReX_TESTING=ON -DWarpX_LIB=ON -DWarpX_DIMS=2 -DWarpX_PSATD=ON -S /tmp/ci-OIrlI1zRvD/warpx/ -B /tmp/ci-OIrlI1zRvD/warpx/builddir 

building VayDeposition2D...
   cmake --build /tmp/ci-OIrlI1zRvD/warpx/builddir -j 8 --  
   Compilation time: 240.433 s
   run & test directory: /tmp/ci-OIrlI1zRvD/rt-WarpX/WarpX-tests/2022-06-23/VayDeposition2D/
   copying files to run directory...
   path to input file: Examples/Tests/VayDeposition/inputs_2d
   running the test...
   mpiexec -n 2 ./VayDeposition2D.ex inputs_2d  diag1.file_prefix=VayDeposition2D_plt   amrex.abort_on_unused_inputs=1 amrex.fpe_trap_invalid=1 amrex.fpe_trap_zero=1 amrex.fpe_trap_overflow=1 warpx.always_warn_immediately=1 warpx.abort_on_warning_threshold=low
   Execution time: 1.216 s
   WARNING: unable to open the job_info file
   doing the analysis...
   WARNING: analysis failed...
./analysis.py  VayDeposition2D_plt000050
Error on charge conservation:
error_rel = 0.010440988968585644
tolerance = 0.001
/tmp/ci-OIrlI1zRvD/py-venv/lib/python3.8/site-packages/yt/utilities/logger.py:4: VisibleDeprecationWarning: The configuration file /home/zaimn/.config/yt/ytrc is deprecated in favor of /home/zaimn/.config/yt/yt.toml. Currently, both are present. Please manually remove the deprecated one to silence this warning.
Deprecated since v4.0.0. This feature will be removed in v4.1.0
  from yt.config import ytcfg
Traceback (most recent call last):
  File "./analysis.py", line 37, in <module>
    assert( error_rel < tolerance )
AssertionError

   Analysis time: 4.227 s
   WARNING: unable to copy analysis image
   archiving the output...
   creating problem test report ...
   VayDeposition2D FAILED

You don't obtain the same thing with this specific test?

Other than that, I agree that generally it is a good idea to use j[idim]->nGrowVect(), except for the specific case of the Vay deposition for which the current density has been obtained with an inverse fft and is nonzero everywhere in the guard cells.

[EZoni]

The suggestion mentioned above in https://github.com/ECP-WarpX/WarpX/issues/3189#issuecomment-1163776204 is now implemented in #3208. Note that the plots shown in the description of this issue are not really of interest, because having the bilinear filter turned off is problematic in itself. I will close the issue. We will see if #3208 makes any difference for the changes reported in #2336 (not sure at this point).