Energy in Simulation Exponentially Increasing

[electroflow]

Hi all,
I am new in the field of PIC simulation and have run into my first big roadblock and was wondering whether someone knows what is going on. In my simulations I always run into a point where the thermal energy of the electrons suddenly grows exponentially and drowns out everything else. Do you have a recommendation which parameter to adjust to control for this? I have attached a simple example of ionized hydrogen at solid densities and the plot of the mean energy per particle for the proton and electron. Here the energy increase isn't quite exponential but nevertheless it is constantly increasing, while the simulation should conserve energy(?)

Inputfile

my_constants.rho = 6e28
my_constants.eV = 1.6e-19
my_constants.c = 299792458

my_constants.energy = 10 * eV
my_constants.th_e = (2 * energy / m_e) ** 0.5 / c
my_constants.th_p = (2 * energy / m_p) ** 0.5 / c

amr.n_cell =  16 16

amr.blocking_factor = 16 
geometry.dims = 2
geometry.prob_lo = 0.0 0.0    
geometry.prob_hi = 0.00125e-6 0.00125e-6 
amr.max_level = 0

stop_time = 2e-10

boundary.field_lo = pec periodic 
boundary.field_hi = pec periodic
boundary.particle_lo = reflecting periodic
boundary.particle_hi = reflecting periodic

particles.species_names = e p

algo.particle_shape = 1

e.species_type = electron
e.injection_style = "NRandomPerCell"
e.num_particles_per_cell = 1
e.momentum_distribution_type = "gaussian"

e.uz_th = th_e
e.uy_th = th_e
e.ux_th = th_e

e.profile = constant
e.density = rho

p.species_type = proton
p.injection_style = "NRandomPerCell"
p.num_particles_per_cell = 1
p.momentum_distribution_type = "gaussian"

p.uz_th = th_p
p.uy_th = th_p
p.ux_th = th_p

p.profile = constant
p.density = rho

diagnostics.diags_names = diag particles
diag.intervals = 100000
diag.diag_type = Full
diag.fields_to_plot = rho

particles.intervals = 100000
particles.diag_type = Full
particles.format = openpmd
particles.write_species = 1
particles.species = e p

evaluation

import openpmd_api as io
from matplotlib import pyplot as plt
import numpy

toev = lambda x : numpy.array(x) / 1.6e-19

ser = io.Series("diags/particles/openpmd_%T.h5", io.Access.read_only)

fullen = {"e":[], "p":[]}

for j in range(len(ser.iterations)):
    i = ser.iterations[j * 100000]
    for part in fullen:
        frm = i.particles[part].to_df()    
        totw = frm["weighting"].sum()
        fullen[part].append((frm["weighting"] * 0.5 * (frm["momentum_x"] ** 2 + frm["momentum_y"] ** 2 + (frm["momentum_z"]) ** 2) / frm["mass"]).sum() / totw)

plt.plot(toev(fullen["e"]), label = "e")
plt.plot(toev(fullen["p"]), label = "p")
plt.ylabel("mean energy in eV")
plt.xlabel("timestep in $10^5$")
plt.legend()

[electroflow]

I have found the different algorithm options and this problem seems to be easily mended by setting algo.field_gathering = energy-conserving. For reference, I have run an example similar to the one above with the different possible settings for current_deposition, field_gathering, and particle_pusher.

[ax3l]

That's right - the PIC algorithm is generally not strictly energy conserving and convergence tests of the resolution and solvers you pick need to be done for the quantities of interest you are looking for.

Great that the energy-conserving gather helps so much in your case!

Things that I generally check to limit/avoid numerical heating on physical time scales:

• make sure the plasma frequency (Langmuir waves) are well resolved
• make sure that the Debye length is resolved
• make sure the shortest wavelength is resolved (e.g., if a laser pulse is in the sim)
• ensure energy-converving gather is used (it is only truly energy-conserving for the limit of infinitesimal small time steps)