Simulation results not corresponding with deck inputs (injector block related)

[DrakeWhu]

Hi EPOCH team,

I am running a simulation where I am injecting a bunch of charged particles inside a hollow plasma tube. The simulations runs succesfully. Here's is the input.deck:

begin:constant
dif_x = 0.01 * micro
dif_y = 0.02 * micro
dt = 0.8 * dif_x/c
n_steps_max = 4000
end:constant

begin:control
nx = 1500
ny = 400

t_end = dt * n_steps_max
x_min = 0
x_max = nx * dif_x
y_min = -ny*dif_y/2
y_max = ny * dif_y/2
stdout_frequency = 100
end:control

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

begin:constant # Tube parameters
plasma_density = 1e28
omega_p = sqrt(plasma_density*qe^2/(me*epsilon0))
lambda_p = 2*pi/omega_p
r_e = 4e-6
r_i = 1e-6
tube_structure = if(((abs(y) gt r_i) or (abs(y) eq r_i)) and ((abs(y) lt r_e) or (abs(y) eq r_e)), plasma_density, 0)
end:constant

begin:constant # Beam params
beta = 0.99
db_mass = 2.0e8
gamma = 1/sqrt(1-beta^2)
beam_density = 2.0e8

drift_p = beta * c * db_mass * gamma # Relativistic momentum
db_charge = -2.0e8
y_fwhm = 3.0e-6
t_fwhm = y_fwhm / c
w_r = y_fwhm / (2.0*sqrt(loge(2)))
w_t = w_r / c
t_hw01m = 0.5 * t_fwhm * sqrt(loge(10)/loge(2))

sprofile_beam = gauss(y, 0, 200 * dif_y)
tprofile_beam = gauss(time, 0, 200 * dif_y / c)
end:constant

begin:species
name = plasma_electrons
charge = -1.0
mass = 1.0
number_density = tube_structure
temp = 0
npart_per_cell = 4
end:species

begin:species
name = driver_beam
charge = db_charge
mass = db_mass
end:species

begin:injector
boundary = x_min
species = driver_beam
t_start = 0
t_end =  40 * dt
number_density =  beam_density * sprofile_beam * tprofile_beam
drift_x = drift_p
npart_per_cell = 10
use_flux_maxwellian = T
end:injector

begin:window
move_window = T
window_v_x = c
window_start_time = (x_max-x_min)/c
bc_x_min_after_move = open
bc_x_max_after_move = open
end:window

begin:output
dt_snapshot = dt * 100

particles = always
px = always
vx = always

ex = always
ey = always
bz = always
number_density = always + species
end:output

The problem I am facing is that changing beta should change the speed of the bunch of particles but it does not, the bunch still is traveling with the moving window, which moves at c. This is even more drastic when setting not relativistic speeds like 0.1c. I'd like to know why this is not working.

I also have questions about the injector block. I am trying to give the injected particle bunch a gaussian shape but apparently it is not working either. The bunch is detectable when you ask for the 1D plot of the middle y direction along the x axis. There we see the bunch is not gaussian in time. Haven't checked the spatial one. One other thing I'd like is for this bunch to have constant charge, so the density/spatial/temporal profile should be adjusted so that charge is constant. I intend to do this because charge is one of the most important parameters of this system and I'd like to change directly the total charge of the bunch if possible.

Lastly, I'd also like to plot a phase space of the particles that cross a certain plane after a given time, that way I can infer if I have acceleration or not.

I asked a bunch of unrelated questions, help with any of them will be greatly appreciated, but I believe the most important part is that I can ensure my colleagues that the results of the simulation reflect the physics we are studying, and the beta and charge of the bunch is something really important to get right for this scan.

For some context, this EPOCH installation was compiled using the HC solver to yield realistic relativistic results, and I believe I am using also intel prefetch.

Thanks, Juan

[Status-Mirror]

Hey Juan,

Sorry for the delay in getting back to you. When you calculate drift_p, you seem to be using $\gamma m c \beta$, but your mass $m$ is in units of electron masses. What you're actually calculating is the momentum for a particle with mass $2\times 10^8$ kg, then assigning this momentum to a particle of mass $2\times 10^8 m_e$. I'm not surprised that your particle is relativistic for various different $\beta$ values.

It's not obvious to me why your profiles wouldn't be working. Could the background plasma be influencing this? It would be better to test your injector by trying to run it in a simulation with a vacuum background.

If you keep your spatial and temporal profiles fixed, then the total injected charge can be modified by just changing the injector number density. If you want to change your profiles, I don't think EPOCH can modify your number density automatically to keep a fixed total charge - the code doesn't have a capability to integrate profiles over space and time.

EPOCH has particle probes which can output the properties of particles as they pass a simulation point/line/plane in 1D/2D/3D. If you're interested in what these particles do later, run the code with particle ID switched on (see compilation options), output the ID of probe-particles, and search for these particles in later outputs.

Hope this helps, Stuart

[DrakeWhu]

Hi Stuart,

I did what you suggested and set the mass taking into account that it is in units of electron masses. I am still debugging this part, gonna launch some simulations to keep testing.

About the spatial and temporal profile of the beam, I managed to change the spatial profile succesfully, but for some reason the time profile is still messed up. I ran a simulation with vacuum background:

input.deck:

begin:constant
dif_x = 0.01 * micro
dif_y = 0.02 * micro
dt = 0.8 * dif_x/c
n_steps_max = 4000
end:constant

begin:control
nx = 1500
ny = 400

t_end = dt * n_steps_max
x_min = 0
x_max = nx * dif_x
y_min = -ny*dif_y/2
y_max = ny * dif_y/2
stdout_frequency = 100
end:control

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

begin:constant # Beam params
beta = 0.99
db_mass = 2.0e8
gamma = 1/sqrt(1-beta^2)
beam_density = 2.0e8

drift_p = beta * c * db_mass * gamma # Relativistic momentum
db_charge = -2.0e8
y_fwhm = 3.0e-6
t_fwhm = y_fwhm / c
w_r = y_fwhm / (2.0*sqrt(loge(2)))
w_t = w_r / c
t_hw01m = 0.5 * t_fwhm * sqrt(loge(10)/loge(2))

sprofile_beam = gauss(y, 0, 10 * dif_y)
tprofile_beam = gauss(time, 0, -t_fwhm)
end:constant

begin:species
name = driver_beam
charge = db_charge
mass = db_mass
end:species

begin:injector
boundary = x_min
species = driver_beam
t_start = 0
t_end =  40 * dt
number_density =  beam_density * sprofile_beam * tprofile_beam
drift_x = drift_p
npart_per_cell = 10
use_flux_maxwellian = T
end:injector

begin:window
move_window = T
window_v_x = c
window_start_time = (x_max-x_min)/c
bc_x_min_after_move = open
bc_x_max_after_move = open
end:window

begin:output
dt_snapshot = dt * 100

particles = always
px = always
vx = always

ex = always
ey = always
bz = always
number_density = always + species
end:output

I got these animations out of the simulation:

This is the derived number density in 2D, where the spatial profile is appreciable:

When we check the time profile along the middle row, we get this shape:

Where we can see the shape of the beam does not match a gaussian, or at least it is a really thight gaussian. Changing the third input in the gaussian of the tprofile_beam doesn't seem to be doing anything.

What I expect to see is typical gaussian bunch with an oval/circular shape in the 2D derived number density, and a gaussian pulse in the 1D derived number density along the middle axis.

Thanks for the help and don't worry about the late response, Juan

[DrakeWhu]

Hi again,

I solved the problem with the bunch. I only had to adjust the parameters of the gaussian. The sprofile and tprofile are now properly set. Here's the updated deck:

begin:constant
dif_x = 0.01 * micro
dif_y = 0.02 * micro
dt = 0.8 * dif_x/c
n_steps_max = 4000
end:constant

begin:control
nx = 1500
ny = 400

t_end = dt * n_steps_max
x_min = 0
x_max = nx * dif_x
y_min = -ny*dif_y/2
y_max = ny * dif_y/2
stdout_frequency = 100
end:control

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

begin:constant # Beam params
beta = 0.99
db_mass = 2.0e8
gamma = 1/sqrt(1-beta^2)
beam_density = 2.0e16

drift_p = beta * c * db_mass * gamma # Relativistic momentum
db_charge = -2.0e8
y_fwhm = 3.0e-6
t_fwhm = y_fwhm / c
w_r = y_fwhm / (2.0*sqrt(loge(2)))
w_t = w_r / c
t_hw01m = 0.5 * t_fwhm * sqrt(loge(10)/loge(2))

sprofile_beam = gauss(y, 0, 10 * dif_y)
tprofile_beam = gauss(time, t_fwhm, 0.1*t_fwhm)
end:constant

begin:species
name = driver_beam
charge = db_charge
mass = db_mass
end:species

begin:injector
boundary = x_min
species = driver_beam
t_start = 0
t_end =  4000 * dt
number_density =  beam_density * sprofile_beam * tprofile_beam
drift_x = drift_p
npart_per_cell = 10
use_flux_maxwellian = T
end:injector

begin:window
move_window = T
window_v_x = c
window_start_time = (x_max-x_min)/c
bc_x_min_after_move = open
bc_x_max_after_move = open
end:window

begin:output
dt_snapshot = dt * 100

particles = always
px = always
vx = always

ex = always
ey = always
bz = always
number_density = always + species
end:output

This yields the expected result. I'm showing here the last frame of the simulation:

Thanks for your help with this issue. I also tried the same with the plasma tube and, as @Status-Mirror suggested, it seems to mess up the bunch. I ran a couple of simulations and all show like if it was a point charge traveling at the speed of light. For example, here's one of the decks where this happens:

begin:constant
dif_x = 0.01 * micro
dif_y = 0.02 * micro
dt = 0.8 * dif_x/c
n_steps_max = 4000
end:constant

begin:control
nx = 1500
ny = 400

t_end = dt * n_steps_max
x_min = 0
x_max = nx * dif_x
y_min = -ny*dif_y/2
y_max = ny * dif_y/2
stdout_frequency = 100
end:control

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

begin:constant # Tube parameters
plasma_density = 1e28
omega_p = sqrt(plasma_density*qe^2/(me*epsilon0))
lambda_p = 2*pi/omega_p
r_e = 3e-6
r_i = 1e-6
tube_structure = if(((abs(y) gt r_i) or (abs(y) eq r_i)) and ((abs(y) lt r_e) or (abs(y) eq r_e)), plasma_density, 0)
end:constant

begin:constant # Beam params
beta  = 0.1
db_mass = 1.0e8
gamma = sqrt(1-beta^2)
db_dens  = 100000000.0

drift_p = beta * c * db_mass * gamma # Relativistic momentum
db_charge = 1.0e8
y_fwhm = 1.0e-6
t_fwhm = y_fwhm / c
w_r = y_fwhm / (2.0*sqrt(loge(2)))
w_t = w_r / c
t_hw01m = 0.5 * t_fwhm * sqrt(loge(10)/loge(2))

sprofile_beam = gauss(y, 0, 10* dif_y)
tprofile_beam = gauss(time, t_fwhm, 0.1 * t_fwhm)
end:constant

begin:species
name = plasma_electrons
charge = -1.0
mass = 1.0
number_density = tube_structure
temp = 0
npart_per_cell = 4
end:species

begin:species
name = driver_beam
charge = db_charge
mass = db_mass
end:species

begin:injector
boundary = x_min
species = driver_beam
t_start = 0
t_end =  40 * dt
number_density =  db_dens * sprofile_beam * tprofile_beam
drift_x = drift_p
npart_per_cell = 10
end:injector

begin:window
move_window = T
window_v_x = c
window_start_time = (x_max-x_min)/(c)
bc_x_min_after_move = open
bc_x_max_after_move = open
end:window

begin:output
dt_snapshot = dt * 100

particles = always
px = always
vx = always

ex = always
ey = always
bz = always
number_density = always + species
end:output

I'm not sure what's happening, could it be related to the issue regarding decks with both a moving window and an injector block?

Thanks for the help, Juan

[Status-Mirror]

Hey Juan,

It's not clear to me that you're seeing something unphysical. Electron bunches are known to have modified profiles when travelling through plasma - this is covered in one of our papers.

If you think the moving window is affecting things, you could try switching this off to see if the pulse changes in an equivalent case, or look at the moving window for an electron bunch in a vacuum. Alternatively, you could run without the injector to see if the background plasma is loaded correctly with the moving window.

Hope this helps, Stuart

[DrakeWhu]

Hi Stuart,

Thanks for the clarification. Accelerator physics is not my forte so some things I see that seem unphysical for me may not be in reality.

Right now I am running a parameter scan to gather more data. Maybe if there is a problem it will show up. I am running 900 different combinations changing plasma density, beta, driver beam mass, charge, density and y_fwhm. I'll share what I find when simulations are finished. Should take a couple hours so for now I just have to wait.

If this scan yields sensible results, I'll close the issue.

Thanks for your kind help, Juan

[Status-Mirror]

Closing the issue for now, but feel free to re-open later.

Cheers, Stuart