Possible wrong charge (and current) deposition at the edges in 2D-RZ

[Yin-YinjianZhao]

@dpgrote @RemiLehe @ax3l

Recently I'm using 2D-RZ to carry out some simulations, and noticed an error occurred at the edges (r=0 and r=r_max) for the rho computation. Although I'm actually using EM solver but only outputting rho, which means the "wrong" rho is not actually used in this simulation since only the current density J is used to solve the Maxwell equations here, but it is worth checking if rho and J are computed correctly in the code.

The wrong result or the big error occurs at the edges as shown in this animation: The plotted rho_e values are averaged over z and plotted over r. The picmi.py script is here: picmi.txt And the plotting script plot_for+gif.ipynb is here: plot_for_gif.txt BTW, I'm using Dave's new branch in this PR #3837 which just fixed the reflection BC.

And usually a correction volume term needs to be used as described in J.P. Verboncoeur's paper: https://doi.org/10.1006/jcph.2001.6923 Or described in a paper of mine, in which I generalized the method to 3D cylindrical charge and current deposition: https://doi.org/10.1007/s40571-022-00513-6

I also noticed the volume correction is related to the PR #3820

Although this error is usually small and may not affects the simulated physics a lot, but it should be fixed if what I said was right.

[RemiLehe]

Thanks for raising this issue.

The above is a bit surprising, since, by default, WarpX does already apply the Verboncoeur correction: https://github.com/ECP-WarpX/WarpX/blob/development/Source/FieldSolver/WarpXPushFieldsEM.cpp#L1395 https://github.com/ECP-WarpX/WarpX/blob/development/Source/WarpX.cpp#L163

Are you suggesting that we should apply a different correction, or do you think that, with your particular input script, this correction is not applied?

[Yin-YinjianZhao]

@RemiLehe The Verboncoeur correction is perfect for this. So I believe in the simulation the correct values are computed. Since my input script gives no more than a random but uniform distribution, and it shows the errors at the edges, I wonder if the Verboncoeur correction is applied in the rho diagnostic, do you know?

[Yin-YinjianZhao]

In other words, does the current deposition and the charge deposition in WarpX's PIC loop both apply Verboncoeur's correction? Or only the current deposition applies? And does the diagnostic outputs use the same functions as the ones in the PIC loop?

[Yin-YinjianZhao]

And recently I'm giving a distribution such that n(r)=3n0 - 3n0*r/R, but again, there is an error on the axis.

picmi.txt

[RemiLehe]

Thanks for doing the above test. Could you also do a test with a constant density (e.g. n(r)=3n0) and compare the results for different number of particles per cell? In my understanding the Verboncoeur correction is calculated for an infinite number of macroparticles per cell, so I would expect the on-axis spike to decrease as you increase the number of particles per cell.

[RemiLehe]

Also: I don't quite remember: isn't the Verboncoeur specifically computed for linear shape factors (i.e. order 1 shape factors)? (See equation 3, 4, 5 in the Verboncoeur paper.)

However, if I understand correctly, the above tests are for order 3 shape factors. Is that correct?

If the above is correct, I think that means that we should apply a correction that depends on the shape factor.

[dpgrote]

The correction that is applied in WarpX is only good for order 1 shape factor and also assumes a uniform distribution. One can use the same procedure and calculate the correction for higher order shape factors (which modifies rho at grid cells further away from the axis with higher order), but that is not implemented in WarpX.

Can you try turning off the correction, setting boundary.verboncoeur_axis_correction = 0 and see what happens?

[Yin-YinjianZhao]

Thank you for your help! @RemiLehe @dpgrote I guess Verboncoeur's correction is only good for 1 shape factor, and I forgot about that. After using 1 shape factor, and increasing the number of particles per cell to 64, I got good result.

Feel free to close this issue, I will use shape factor=1 for the moment.

[RemiLehe]

@dpgrote @Yin-YinjianZhao Do you know of a paper that gives the Verboncoeur correction for shape 2 and 3? Also: @dpgrote Do you remember whether the correction for shape 2 and 3 where implemented in Warp?

[dpgrote]

I have not seen it published anywhere for higher order shape factors. I remember that it is relatively easy to derive, something like the integral of the shape factor between the grid cells of a particle on the axis. I had written a script to calculate the correction many years ago but can't find it now. It was implemented for 2nd order shape in Warp in on place.

[Yin-YinjianZhao]

However, when I try to use a constant uniform density profile along r, I found the error at r=0 can be fixed by shape factor=1, but there is still an error at the bigger edge. This may needs to be double checked in the code if the bigger edge is handled correctly.

picmi.txt

[dpgrote]

Currently, there is no correction applied to the boundary at the upper radius. A Verboncoeur type correction can be derived at the upper boundary, but it is not applied in WarpX. This should be relatively straight forward to implement, with the same assumptions of uniform density and space factor of 1.

[Yin-YinjianZhao]

@dpgrote @RemiLehe

It could be something like Eq.(12) and Eq.(13) in 3D cylindrical.

Or like the following code in 2D RZ, where i=0 is the axis, i1 and i2 are two intermediate boundaries we don't need to care about, and i3 is the bigger r boundary.

subroutine sub_scatter_2d_rz_n(i1,i2,i3,k1,k2,k3,npml,den,par,w,dr,dz)

    implicit none
    integer :: i1,i2,i3,k1,k2,k3,npml
    real,dimension(0:i3+1,0-1-npml:k3+1+npml) :: den
    real,dimension(:,:) :: par
    real :: w,dr,dz

    integer :: ip,i,k
    real :: r1,r2,z1,z2,rp,zp
    real :: vi0,vi1,vk0,vk1
    real :: dr2
    real,parameter :: pi2 = 6.283185307179586d0

    dr2 = dr*dr

    do ip = 1, size(par,2)

        rp = sqrt(par(1,ip)**2+par(2,ip)**2)
        zp = par(3,ip)

        i = floor(rp/dr)
        k = floor(zp/dz)

        r1 = (  rp -  i   *dr ) / dr
        r2 = ( -rp + (i+1)*dr ) / dr
        z1 = (  zp -  k   *dz ) / dz
        z2 = ( -zp + (k+1)*dz ) / dz

        if ( i==0 ) then
            vi0 = dr2/6
        else if ( i==i2 .and. k<k2 ) then
            vi0 = (3*i2+1)*dr2/6
        else
            vi0 = i*dr2
        end if

        if ( i+1==i3 ) then
            vi1 = (3*i3-1)*dr2/6
        else if ( i+1==i1 .and. k<k2 ) then
            vi1 = (3*i1-1)*dr2/6
        else
            vi1 = (i+1)*dr2
        end if

        if ( k==0 ) then
            vk0 = 0.5*dz
        else if ( k==k1 .and. i<=i1 ) then
            vk0 = 0.5*dz
        else if ( k==k2 .and. i>i1 .and. i<i2-1 ) then
            vk0 = 0.5*dz
        else
            vk0 = dz
        end if

        if ( k+1==k3 ) then
            vk1 = 0.5*dz
        else
            vk1 = dz
        end if

        den(i  ,k  ) = den(i  ,k  ) + w*r2*z2/(vi0*vk0*pi2)
        den(i+1,k  ) = den(i+1,k  ) + w*r1*z2/(vi1*vk0*pi2)
        den(i  ,k+1) = den(i  ,k+1) + w*r2*z1/(vi0*vk1*pi2)
        den(i+1,k+1) = den(i+1,k+1) + w*r1*z1/(vi1*vk1*pi2)

    end do

end

The above is for the (charge) density, for the current density, this may be helpful, which I just wrote recently for my own code.

subroutine sub_scatter_2d_rz_f_f(i1,i2,i3,k1,k2,k3,npml,Jr,Jz,q,w,dr,dz,dt,r1,r2,z1,z2,i,k)

    implicit none
    integer :: i1,i2,i3,k1,k2,k3,npml
    real,dimension(0:i3+1,0-1-npml:k3+1+npml) :: Jr,Jz
    real :: q,w,dr,dz,dt
    real :: r1,r2,z1,z2
    integer :: i,k

    real :: del_r,del_z,rmid,zmid,dr2
    real :: vi0,vi1,vk0,vk1
    real,parameter :: pi2 = 6.283185307179586d0
    real :: V1,V2,V,buf

    del_r = r2-r1
    del_z = z2-z1
    if ( abs(del_r)<tiny(1.0) .and. abs(del_z)<tiny(1.0) ) return

#if defined(DEBUG)
    if (i>=i3+1 .or. i<0-1) then
        write(*,*) "Yin: sub_scatter_2d_rz_f_f: i>=i3+1 or i<0-1."
        write(*,*) "i,i3,k,k3:",i,i3,k,k3
        write(*,*) "r1,r2,z1,z2",r1,r2,z1,z2
        stop
    end if
    if (k>=k3+1 .or. k<0-1) then
        write(*,*) "Yin: sub_scatter_2d_rz_f_f: k>=k3+1 or k<0-1."
        write(*,*) "i,i3,k,k3:",i,i3,k,k3
        write(*,*) "r1,r2,z1,z2",r1,r2,z1,z2
        stop
    end if
#endif

    rmid = 0.5*(r1+r2)
    zmid = 0.5*(z1+z2)

    dr2 = dr*dr

    if ( i==0 ) then
        vi0 = dr2/6
    else if ( i==i2 .and. k<k2 ) then
        vi0 = (3*i2+1)*dr2/6
    else
        vi0 = i*dr2
    end if

    if ( i+1==i3 ) then
        vi1 = (3*i3-1)*dr2/6
    else if ( i+1==i1 .and. k<k2 ) then
        vi1 = (3*i1-1)*dr2/6
    else
        vi1 = (i+1)*dr2
    end if

    if ( k==0 ) then
        vk0 = 0.5*dz
    else if ( k==k1 .and. i<=i1 ) then
        vk0 = 0.5*dz
    else if ( k==k2 .and. i>=i1 .and. i<i2 ) then
        vk0 = 0.5*dz
    else
        vk0 = dz
    end if

    if ( k+1==k3 ) then
        vk1 = 0.5*dz
    else
        vk1 = dz
    end if

    V1 =  ((k+1)*dz-zmid)*abs(del_r)
    V2 = -((k  )*dz-zmid)*abs(del_r)
    V  = V1+V2
    if ( V<tiny(1.0) ) V=1.0d0
    buf = q*w*del_r/dt/((dble(i)+0.5d0)*dr2)/pi2
    Jr(i,k)   = Jr(i,k)   + buf*V1/V/vk0
    Jr(i,k+1) = Jr(i,k+1) + buf*V2/V/vk1

    if ( i==0 ) then
        vi0 = dr2/6
    else if ( i==i2 .and. k<k2 ) then
        vi0 = (3*i2+1)*dr2/6
    else
        vi0 = i*dr2
    end if

    if ( i+1==i3 ) then
        vi1 = (3*i3-1)*dr2/6
    else if ( i+1==i1 .and. k<k2 ) then
        vi1 = (3*i1-1)*dr2/6
    else
        vi1 = (i+1)*dr2
    end if

    if ( k==0 ) then
        vk0 = 0.5*dz
    else if ( k==k1 .and. i<=i1 ) then
        vk0 = 0.5*dz
    else if ( k==k2 .and. i>i1 .and. i<i2-1 ) then
        vk0 = 0.5*dz
    else
        vk0 = dz
    end if

    if ( k+1==k3 ) then
        vk1 = 0.5*dz
    else
        vk1 = dz
    end if

    V1 =  ((i+1)*dr-rmid)*abs(del_z)
    V2 = -((i  )*dr-rmid)*abs(del_z)
    V  = V1+V2
    if ( V<tiny(1.0) ) V=1.0d0
    buf = q*w*del_z/dt/dz/pi2
    Jz(i,k)   = Jz(i,k)   + buf*V1/V/vi0
    Jz(i+1,k) = Jz(i+1,k) + buf*V2/V/vi1

end