Matching Staggered Field with Charge on the Grid

[rmcg15]

Hello,

I've been dealing with the problem of initialising the fields of a charge/current on the grid and attempting to resolve artefacts arising from the initialisation. My method is to calculate the initial fields of a distribution externally and load the data files for the fields such that deposited charge begins the simulation with the correct static fields.

At present it appears that I have a mismatch between the charge and fields owing to their relative positioning. See in the two plots below. A small, relativistic bunch of charge is deposited in the middle of the grid
`
begin:species
name = distribution
charge = -1.0
mass = 1.0
fraction = 1.0
drift_px = drift_p
number_density = dens * gauss(x,0,sqrt(2)*wdthx) * gauss(y,0,sqrt(2)*wdthy)

bc_x_min = open
bc_x_max = open
bc_y_min = open
bc_y_max = open

end:species
`
The fields are loaded externally, and when propagating I see the following artefact.

If I somewhat arbitrarily shift the grid that my field is evaluated on by a value of the order of my grid dimensions then the artefact is reduced (but in the other hand the transverse field looks cruddier).

Clearly this indicates that the fields and charge are mismatched wrt to eachother and this is the cause of the artefacts. However I find that the documentation is not very clear on how the charge and field are staggered in space and time. Among other things there's also a mention of ghost cells, do I need to account for this in my data file?

Thanks

[TomGoffrey]

I'm not familiar with the code that reads fields from file (so another developer might add to this), but I would expect that you do need to take into account of the fact that fields are staggered in space using the standard Yee grid. Your observations about the changes you see from shifting the field seems to back that up. Ghostcell values will be set by EPOCH, so there's no need for you to do anything about those.

[rmcg15]

Hi Tom,

Thanks for your response. Currently when solving for fields I offset Ex by dx/2 (where dx is the cell dimension in x), Ey by dy/2 and Bz by dx/2 and dy/2. It is my understanding that EPOCH uses a leapfrog algorithm such that the fields and charge are evaluated such that there's no off-set in time (so it calculates values such as Ex[t+dt], Bz[d+dt] as opposed to Ex[t+dt], Bz[d+dt/2]).

If this correct then it indicates that perhaps something is going awry in my matlab script elsewhere (I've been adding some bells and whistles to it so perhaps something is going wrong there).

[Status-Mirror]

Hi @rmcg15,

There is a known bug with this file-reader, where I think the spatial stagger for loading fields is in the wrong direction. See issue 597 for a more complete discussion and workaround.

Unfortunately I haven't had the chance to fix this recently, as we're busy with updating the whole code to C++.

Hope this helps,
Stuart

[rmcg15]

Hi @rmcg15,

There is a known bug with this file-reader, where I think the spatial stagger for loading fields is in the wrong direction. See issue 597 for a more complete discussion and workaround.

Unfortunately I haven't had the chance to fix this recently, as we're busy with updating the whole code to C++.

Hope this helps, Stuart

Hi,

Firstly apologies for the late response, it took awhile for me to apply the suggested fix in the thread as I am not running it locally.

I tried applying the fix suggested by the author of the thread by increasing cell_x2, cell_y2 by one however from my testing there is apparently no change in performance, I obtain the exact same artefacts as prior. Is there something in the fields block I should be editing instead as implied in the following discussion in the thread?

Regards, Ryan

[rmcg15]

bump.

[Status-Mirror]

Good use of a bump! I've made some progress but I forgot to update this. I was going through the source-code to fix issue 597, where I learned this error was only related to using maths parser expressions in the fields block, and it shouldn't affect fields set from files.

You are to write your fields files assuming that everything is evaluated at the same point in time, $t=0$. This is also the assumed time when particles are loaded. EPOCH will perform a half-step in the magnetic field before the PIC-loop begins, to allow the standard leap-frog method.

If loading from a file, you are to give the code the values of the fields at their evaluation points (including grid-stagger). I'll include an example here.

Suppose we have 6 cells, from 0 to 6 microns, in a 1D simulation. The $E_y$ evaluation points are cell-centred, and the $E_x$ points lie on the high $x$ cell edge. Hence, the 6 evaluation points to load in your fields block correspond to:

• $E_x$ : (1, 2, 3, 4, 5, 6) microns
• $E_y$ : (0.5, 1.5, 2.5, 3.5, 4.5, 5.5) microns

With outflow-boundaries, any fields on the simulation edge will be assumed 0. This includes the $E_x$ point on x_min, and the value at x_max which you do give. This will be changed by the boundary routines before your simulation begins.

To test how the fields are loaded relative to particles, let us create a set of fields which are anti-symmetric about a point at 3 microns. This means that a particle positioned exactly at 3 microns will feel no acceleration from $E_x$ and $E_y$, as they will be balanced on both sides. To do this, we will use the fields block to setup the fields, and we'll use the particles_from_file block to position our test particle. If our test-particle has a net acceleration, we know we've got the stagger wrong.

To write the binary files, I used this MATLAB script:

% Simulation parameters
nx = 6;

% Setup anti-symmetric fields with stagger
ex_vals = 1.0e10 * [-1, -1, 0, 1, 1, 0];
ey_vals = 1.0e10 * [-1, -1, -1, 1, 1, 1];

% Write field files
ex_file = fopen('ex.bin','w');
ey_file = fopen('ey.bin','w');
for i = 1:nx
    fwrite(ex_file,ex_vals(i),'double');
    fwrite(ey_file,ey_vals(i),'double');
end
fclose(ex_file);
fclose(ey_file);

% Write files to load a particle at the zero-field point at 3um
x_file = fopen('x.bin', 'w');
w_file = fopen('w.bin', 'w');
fwrite(x_file,3.0e-6,'double');
fwrite(w_file,1,'double');
fclose(x_file);
fclose(w_file);

And I used this input deck, where I've given the absolute path to the data files on my system:

begin:control
    nx = 6
    x_min = 0
    x_max = 6e-6
    nsteps = 2
    stdout_frequency = 100
end:control

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

begin:fields
    ex = "/home/stu/Debugging/disc_685/ex.bin"
    ey = "/home/stu/Debugging/disc_685/ey.bin"
end:fields

begin:species
    name = Electron
    mass = 1.0
    charge = -1.0
end:species

begin:particles_from_file
    species = Electron
    x_data = "/home/stu/Debugging/disc_685/x.bin"
    w_data = "/home/stu/Debugging/disc_685/w.bin"
end:particles_from_file

begin:output
    dt_snapshot = t_end
    px = always
    ex = always
    ey = always
end:output

When running this, I find that in the 0001.sdf output file, the particle has 0 momentum. Hence, fields must be evaluated at their stagger co-ordinates to interact with particles correctly. If I ignore stagger and write ex_vals = 1.0e10 * [-1, -1, -1, 1, 1, 1];, then the particle gains momentum 5.9856e-24 kgm/s by the simulation end.

Hope this helps,
Stuart

[rmcg15]

Hi Stuart,

Thanks very much for your detailed response. I think your outlined procedure for obtaining the correct stagger should match up with what I've done. I did have a couple of subtle errors in my initialisation which has slightly reduced the severity of the artefacts I'm seeing but they are still present. I perhaps wonder if there's still a mismatch somewhere.

In an example, I have a 1000x1000 grid over a domain of 3000microns so dx = 3microns. Keeping to the x-direction, If I deposit charge according to the input deck I gave above then I assume my charge is deposited symmetrically on the grid according to distribution function f(x,y), we could say if the total width is 30microns, so it extends 15 microns either way along x, occupying 5 cells on one side and 5 cells on the other. If I understand the depositing correctly, the points at which the charge is deposited (i.e. points at which f(x,y) is sampled to give the initial number density) along the x axis would be [-13.5 -10.5 -7.5 -4.5 -1.5 1.5 4.5 7.5 10.5 13.5 ].

In my initialisation code in very basic terms, I recreate a grid with the spacing discussed above, and then for each component of the field I evaluate at the grid boundaries:
ExPoints = [-12, -9, -6,-3, 0, 3, 6, 9, 12, 15]
EyPoints = [-13.5 -10.5 -7.5 -4.5 -1.5 1.5 4.5 7.5 10.5 13.5 ]

• the rest of space etc. etc.

What I wonder is if where I expect the charge to be deposited on the grid and the actual position EPOCH deposits as a function f(x,y) are different. It does confuse me somewhat that in my data output I have a number-density grid that is 1000x1000 but the grid coordinates are 1001x1001.

It's possible the degree to which I can initialise this system correctly is limited however, when I play and shift the off-set of the fields and charge slightly, I tend to see a visible amount of anomalous radiation, so I might have gotten things as good as I can.

Thanks again

Ryan

[Status-Mirror]

Hi Ryan,

When loading macro-particles, the function given to your density key in the species block will evaluate that function at the $x$ and $y$ positions of the cell centre. This will set the number density for that cell.

In the default particle-loader, the code then randomly positions macro-particles at various points throughout that cell, and gives them a weight such that the sum of the weights is equal to the local number density multiplied by the cell volume. In a 2D code, the cell-size in $z$ is assumed to be 1m for the volume calculation - this actually doesn't matter, as the field-solver is only sensitive to current densities, not absolute current.

So your density won't give you the smooth analytic function, but will instead give you a blocky distribution, with structure on the order of the cell-size. But there are additional caveats, as the macro-particles also have shape functions which project themselves over neighbouring cells. A single macro-particle actually represents a cloud of real particles spread over multiple cells. This adds a second departure from a clean $f(x,y)$ function.

These issues are resolution based, and I would expect them to be reduced if you switch to a system with a smaller cell-size. You can try setting the particle positions and weights manually like in the particles from file demo, but I think it would be hard to get smooth functions due to the particle shapes.

Another issue could be due to the beam itself. If you have too few macro-particles representing this beam, then any currents would introduce noise, which would affect the field-solver. It's also been a good few years since I took an electrodynamics undergraduate module, so I can't comment on whether you have the right fields for your given electron bunch.

On your final question, your number density grid is 1000x1000 cells to match your 1000x1000 cell simulation, but EPOCH outputs the cell-edges for its grid, which is why there is an extra point to give 1001x1001.

Hope this helps,
Stuart

[rmcg15]

Hi Stuart,

Thanks again for your responsive, this is very informative. It sounds very much like my issues could as you say be resolution based, previously I'd had some issues with noise (especially with too few macroparticles) even when loading a "clean" distribution in which I had concluded meant there was some random sampling across the grid cells occurring which you describe.

I wonder then if I introduce some aspects of the weighting to my distribution then it might help smooth out these remaining issues . . .

Thanks! Ryan

[Status-Mirror]

Hi Ryan,

I'd advise caution in that approach. You may be able to set up your system so it has a perfect initial charge distribution, but the steps are still evaluated using cell-discretised EM fields and current densities. In my experience, noise in the current densities (and hence the fields) are unavoidable in PIC codes. I think your safest bet would be to reduce the cell-size in your current set-up.

If your bunch has width 30 microns in a simulation with a 3 micron cell size, that's only 10 cells across. The particle-shape overflow would be an extra cell on each side, so a 20% modification to your intended shape.

If your simulations are quite expensive in 2D, it might be worth testing the effect of cell-size on your anomalous fields in 1D simulations first.

Cheers,
Stuart

[rmcg15]

Hi Stuart,

Yes that's a good point. Ultimately in a low resolution system the difference between the distributions might be quite significant. It does appear that increasing the grid size helps reduce this error. What is ultimately frustrating is that I'd like to simulate quite relativistic electrons, I find that these errors tend to increase quickly as I increase the energy/velocity of the charges. So even if they appear to vanish for N=5000 for a lorentz factor of 5 they may appear again when I move to a lorentz factor of 50.

Without a Lorentz boost I may have to content myself with having the initialisation work as well as I can and studying lower energy cases.

Regards, Ryan