Direction Splitting: Self-Consistent J & EMF

[ax3l]

This pull request introduces the correct handling for the binding of electro-magnetic fields to currents in the directional splitting (DS) scheme.

Rationale

In the directional splitting cell only the electro-magnetic fields are centered¹. Both B and E are furthermore defined at the same time t.²

The particles are still moved between the updates of B and E making their current density J be defined at t-1/2.

Let's have a look at the cell for our DS implementation³, E^t and B^t:

Compared to the Yee cell³, E^t and B^t-1/2:

Nevertheless the source terms³ (only J printed) are still at:

the last image can be derived from either analytically checking the DS scheme for a plane wave equation⁴ or by performing the substitution in the DS operators F^+-, G^+- in 2 or 3D and solving for J self-consistently, example for 1D in⁵ (too long for Markdown, I need LaTeX for that).

Fun fact: we can now still use our current current solvers since the position of J stays the same¹.

Nevertheless, the result is that J can not simply be subtracted from E. But simply put, we can just interpolate the 4 closed J values along the plane that is perpendicular to the J component (e.g.: x-y plane for E_z: J_z^ijk, J_z^i+1jk, J_z^ij+1k, J_z^i+1j+1k) and use the averaged results of each component for E (long reason: see derivation; E_x(J_x), E_y(J_y), E_z(J_z)).

Since the E and B field are calculated at the same point in time (and space) (and are not staggered in those as in Yee's scheme), the influence of J should propagate to both immediately. We now take Faraday's Law and split E^t due to superposition to E^t,free(E^t-1,B^t-1) and the delta E^t that is introduced from J^t-1/2. Now place the correct spatial averaged J's in again (I need LaTeX or will scan my notebook an other day for that) and one comes to the result that the change delta B is basically a specifically spatially averaged curl of J (dependencies: B_x(J_y,J_z), B_y(J_x,J_z), B_z(J_x,J_y) ).

Implementation

Based on the (for that reason) introduced current interpolation/filters from #888 (merged).

See src/picongpu/include/fields/currentInterpolation/NoneDS/NoneDS.hpp for the main implementation.

I renamed AllCenteredCell to EMFCenteredCell to reflect the actual distribution of quantities. The inherent assumption that the position of J depends on the position of E of a NumericalCellType was removed.

The missing event dependency in DS's after_current implementation was added.

Current Status

I just published this rather untested branch (1,5 days from scratch to result) since other developers were interested in the result. The KHI test is not yet correct.

Needs double checking: is the trilinear interpolation correct for the EMFCenteredCell? (needs -1 and +1 now for NGP!)⁶

Compile check: just a CFL assert due to my [WIP] commit that changes a file temporarily.

I did not cross-check my derivations (yet) completely and also need some more minutes to cross-check my constants. Furthermore, there are still several issues in the propagation of the fields #889 #890 #891 #893 ...

Thx a lot to @psychocoderHPC for event system and implementation support and for debugging the implementation of the vacuum field solver with @Heikman . Thx to @stetie for great general discussions and literature - why are all derivations 1D only? ;)

Footnotes

¹ One can place them also at some other point in the cell, e.g., the corners but we will see that this definition helps us in our implementation.
² Which btw. reduces some pusher artifacts that can lead to spurious emittance growth.
³ Examples taken and modified under CC-BY-SA 4.0 from my diploma thesis
⁴ Talk by Yasuhiko Sentoku, University of Nevada - Reno: Numerical Dispersion Free Maxwell Solver for multi-dimensional PIC
⁵ R. Lichters, R.E.W. Pfund, J. Meyer-ter-Vehn, MPI für Quantenoptik - Garching, LPIC++ A Parallel One-dimensional Relativistic Electromagnetic Particle-In-Cell Code for Simulating
Laser-Plasma-Interaction, chapter 5, https://books.google.de/books?id=FvkMHAAACAAJ (update: there is a download link including doc/book; code under GPLv2+; last update 07/2014 with v1.3.2)
⁶ General but independent side note: we could implement two additional interpolation methods... Galerkin/"energy conserving" and "momentum conserving" via interpolation to grid nodes first; the latter might be the same as trilinear for DS (probably not, since rho is not at E and B positions).

CC'ing @psychocoderHPC @Heikman @DrThomasKluge @bussmann @PrometheusPi @stetie

[ax3l]

@psychocoderHPC are you sure that the FieldToParticleInterpolation (currently AssignedTrilienarInterpolation) takes into account that it needs +1 cells in the lower direction as soon as the fields are not defined on the lower side of the cell?
Update: looks ok since the shift to the coordinate system of the fields should account for that.

[ax3l]

according to Sentoku's derivation,⁴ this should be 0.25.
changing this reduces the heating but is still too high (from the 10% to the few-% level in 2k steps) with TSC & Esirkepov and the KHI still shows artifacts.

[ax3l]

omg I am stupid, my Friday derivation also contains 0.25 xD

[ax3l]

micro optimization: actually a scalar such as GetComponentsType<Memory::ValueType>::type::create(0.5) would be enough (currently, it's a vector with all elements 0.5).

[bussmann]

I'd like to see the Galerkin Methods in PIConGPU. But first let's finish DS.

[psychocoderHPC]

here is a small update with results (debug pictures) of this pull request.

KHI el/pos with DS

const double DELTA_T_SI = 1.79e-16;
const double CELL_WIDTH_SI = DELTA_T_SI*SPEED_OF_LIGHT_SI;

LWFA with DS

const double CELL_WIDTH_SI = 0.4430e-7;
const double DELTA_T_SI = CELL_WIDTH_SI/ SPEED_OF_LIGHT_SI;

[psychocoderHPC]

@PrometheusPi and me checked with a plane wave in vacuum why the laser in DS is not as strong as with Yee (pictures in my previous post).
The problem in DS is that the laser amplitude is smaller by a factor of two. This originates from the missing B-field initialization in the laser initialization. The result in DS is that the wave is spitted in a forward and backward propagated wave. ( @PrometheusPi can explain this better)

possible solutions:

• add B-field while the laser injection (hard for laser with focus)
• add a factor of two to the laser amplitude @KlausSteiniger^^
• for DS only add the laser to the positive moving wave

Update: This problem is moved to a separate issue #916

[ax3l]

thx! Pls post this in the according pull for the solver. This pull only fixes the current interpolation (and the wrong laser init for DS does not explain a wrong khi).

Furthermore, using an antenna is probably easier than e+b init and should be implemented. Moving antennas also open the way to boosted frame sims.

[ax3l]

@psychocoderHPC this implementation is actually ready for a merge but does not implement all the missing points for a full DS.

someone or I will have to take a look if the field to particle interpolation is done with the correct fields in space and time, but this is totally independent of this problem.

same for the cross-check of the field solver (not only the #916 issue...)

[psychocoderHPC]

well done and big thanks for the nice description!