Warwick-Plasma/epoch

Possible deviation when interpolating the fields felt by the particles

XiangyanAn opened this issue · 4 comments

I just found that the EPOCH codes are one-cell shifted when calculating the nearest boundary cell of the particles and interpolating the fields felt by the particles. Such a deviation would result in a numerical increase of a conserved quantity when the electrons interact with a circular-polarized laser. This is confirmed by reading the source code and checking the fields of the particles after setting the fields as E_{x,y,z}=B_{x,y,z}=x. The numerical increase is fixed by adding one to the variables cell_x2, cell_y2, and cell_z2.

I wrote a detailed description in the attached document.

epoch_interpolation.pdf

Hey @XiangyanAn,

Thanks for the detailed write-up, this is definitely something worth checking up.

Your pdf file mentions a few different input decks - one with a circularly polarised beam, and another with the linear fields. Could you also send these if you still have them? I can have a go at recreating them if not.

Cheers,
Stuart

I'm still going through your report, but from a partial read I can see an EPOCH bug in the fields block. I followed your method of setting $E$ and $B$ to $x$, but in a 1D simulation. Input deck attached at the end.

In this simulation, we have 10 cells ranging from 0m to 10m, such that we have:

  • Cell edges (m) = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] (11 edges)
  • Cell centres (m) = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5] (10 cells)

Using the developer documentation stagger guide, or Fig 3b in your report, we expect the 10 cells of Ey, Ez and Bx to match the cell centre numbers, and this is what we observe. You acknowledge this agreement in figures 2b, 2c, 2e.

For staggered field-points, the Yee grid suggests a stagger in the positive direction. If cell with index ix has edges at 3m and 4m, then Ex(ix) should be evaluated on the high $x$ cell-edge, or at 4m in this case. Hence, I expect the $x$-stagger fields to be at positions:

  • Expected $x$ stagger fields = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

We observe:

  • Ex = [0, 1, 2, 3, 4, 5, 6, 7, 8, 0]
  • By = [0, 1, 2, 3, 4, 5, 6, 7, 8, 8.6780]
  • Bz = [0, 1, 2, 3, 4, 5, 6, 7, 8, 8.6780]

Even ignoring the strange value on the high $x$ simulation edge, we can see that these fields are being treated as if they evaluated on the low $x$ cell-edge, which is inconsistent with how the rest of the code treats grid stagger. This is why you were able to "correct" the problem by adjusting the interpolation. I would instead argue that it is not the interpolation which is broken, but the fields block.

With my understanding of the problem, this issue should only arise if you're using the fields block. When setting up your circular polarisation in equations 1 and 2, did you use the fields block, or the laser block?

Cheers,
Stuart

begin:control
    nx = 10
    t_end = 1.0e-15
    x_min = 0
    x_max = 10
    stdout_frequency = 100
end:control

begin:boundaries
    bc_x_min = simple_laser
    bc_x_max = open
end:boundaries

begin:fields
    ex = x
    ey = x
    ez = x
    bx = x
    by = x
    bz = x
end:fields

begin:output
    dt_snapshot = t_end
    ex = always
    ey = always
    ez = always
    bx = always
    by = always
    bz = always
end:output

I'll also take a look at this tomorrow.

I'm still going through your report, but from a partial read I can see an EPOCH bug in the fields block. I followed your method of setting E and B to x, but in a 1D simulation. Input deck attached at the end.

In this simulation, we have 10 cells ranging from 0m to 10m, such that we have:

  • Cell edges (m) = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] (11 edges)
  • Cell centres (m) = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5] (10 cells)

Using the developer documentation stagger guide, or Fig 3b in your report, we expect the 10 cells of Ey, Ez and Bx to match the cell centre numbers, and this is what we observe. You acknowledge this agreement in figures 2b, 2c, 2e.

For staggered field-points, the Yee grid suggests a stagger in the positive direction. If cell with index ix has edges at 3m and 4m, then Ex(ix) should be evaluated on the high x cell-edge, or at 4m in this case. Hence, I expect the x-stagger fields to be at positions:

  • Expected x stagger fields = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

We observe:

  • Ex = [0, 1, 2, 3, 4, 5, 6, 7, 8, 0]
  • By = [0, 1, 2, 3, 4, 5, 6, 7, 8, 8.6780]
  • Bz = [0, 1, 2, 3, 4, 5, 6, 7, 8, 8.6780]

Even ignoring the strange value on the high x simulation edge, we can see that these fields are being treated as if they evaluated on the low x cell-edge, which is inconsistent with how the rest of the code treats grid stagger. This is why you were able to "correct" the problem by adjusting the interpolation. I would instead argue that it is not the interpolation which is broken, but the fields block.

With my understanding of the problem, this issue should only arise if you're using the fields block. When setting up your circular polarisation in equations 1 and 2, did you use the fields block, or the laser block?

Cheers, Stuart

begin:control
    nx = 10
    t_end = 1.0e-15
    x_min = 0
    x_max = 10
    stdout_frequency = 100
end:control

begin:boundaries
    bc_x_min = simple_laser
    bc_x_max = open
end:boundaries

begin:fields
    ex = x
    ey = x
    ez = x
    bx = x
    by = x
    bz = x
end:fields

begin:output
    dt_snapshot = t_end
    ex = always
    ey = always
    ez = always
    bx = always
    by = always
    bz = always
end:output

Thanks for your reminder. Since I wanted a simplest problem in physics to make a test, I did use the fields block rather than the laser block to make the plane wave.

I have read the codes of 'fields.F90' and the developer documentation again and made some more tests. I now understand that the Yee grid suggests a stagger in the positive direction for staggered field-points. This means that the index of Ex, By, and Bz are 1 less than that of the corresponding xb grid, while I took them as the same in the previously written document.

I think you are right. It is not the interpolation which is broken, but the fields block.