1D diode at space charge limit

[AlexGlock]

Hello everybody,
I recently started with WarpX and I am now working on a very simple 1D model of a diode. I really like the picmi python interface and would prefer to work with it, so I started with the capacitive discharge python example and came up with an picmi_input file for the diode (as attached). It does run and gives me the desired diagnostic outputs - however, there seem to be no particles in the results. I strongly suspect the UniformFluxDistribution to be faulty but so far I didnt manage to spot the error. Maybe someone else has done something like this and knows the common pitfalls.

from pywarpx import picmi
import numpy as np
import scipy as sp
import math

# --------------------------------------------------
# Input   0 o----------------o L 
#       emitter (Phi_0)   collector (Phi_L)
# --------------------------------------------------
gap = 5e-3      # gap distance [m]
Phi_0 = 0.0     # emitter potential [V]
Phi_L = 10.0    # collector potential [V]

t_end = 1e-7    # end time  [s]
nz = 128        # number of cells
energy = 1      # initial electron energy [eV]

# --------------------------------------------------
# Constants
# --------------------------------------------------
gap_voltage = Phi_L - Phi_0
m_e = sp.constants.m_e 
q_e = sp.constants.e 
eps_0 = sp.constants.epsilon_0
rms_velocity = np.sqrt(2 * energy * q_e / m_e)

dz = gap / nz
vzfinal = np.sqrt(2.0 * abs(gap_voltage) * q_e / m_e)
dt = 0.4 * (dz / vzfinal)
timesteps = math.trunc(t_end/dt)

# Child–Langmuir current density
J_CL = ((4.0 / 9.0) * eps_0 *
        np.sqrt(2.0 * abs(q_e) / m_e) *
        abs(gap_voltage)**(3.0 / 2.0) /
        gap**2)

# --------------------------------------------------
# Grid (1D)
# --------------------------------------------------
grid = picmi.Cartesian1DGrid(
    number_of_cells=[nz],
    lower_bound=[0.0],
    upper_bound=[gap],
    lower_boundary_conditions=['dirichlet'],
    upper_boundary_conditions=['dirichlet'],
    lower_boundary_conditions_particles=['absorbing'],
    upper_boundary_conditions_particles=['absorbing'],
    warpx_potential_lo_z=Phi_0,
    warpx_potential_hi_z=Phi_L
)

# --------------------------------------------------
# Electrostatic solver
# --------------------------------------------------
solver = picmi.ElectrostaticSolver(
    grid=grid
    #method='Multigrid'
)

# --------------------------------------------------
# Electron species
# --------------------------------------------------
electrons = picmi.Species(
    particle_type='electron',
    name='electrons',
    initial_distribution=picmi.UniformDistribution(
        density=0.0,

    )
)

# --------------------------------------------------
# Continuous injection (flux at cathode)
# --------------------------------------------------
electron_flux = picmi.UniformFluxDistribution(
    flux_normal_axis='z',
    surface_flux_position=0.0,
    flux_direction=1,
    flux=str(J_CL / q_e),
    rms_velocity=[0., 0., rms_velocity]
)

# --------------------------------------------------
# Simulation
# --------------------------------------------------
sim = picmi.Simulation(
    solver=solver,
    time_step_size=dt,
    max_time=t_end,
    particle_shape='linear'
)

sim.add_species(
    species=electrons,
    layout=picmi.GriddedLayout(
        n_macroparticle_per_cell=[15]
    )
)

# --------------------------------------------------
# Diagnostics
# --------------------------------------------------
field_diag = picmi.FieldDiagnostic(
    name='diag1',
    grid=grid,
    period=timesteps,
    data_list=['Ez', 'rho', 'Jz', 'phi'],
    write_dir='diags',
    warpx_format='openpmd'
)

particle_diag = picmi.ParticleDiagnostic(
    name='diag1',
    species=electrons,
    period=timesteps,
    data_list=['position', 'momentum'],
    write_dir='diags',
    warpx_format='openpmd'
)

sim.add_diagnostic(field_diag)
sim.add_diagnostic(particle_diag)

# --------------------------------------------------
# Run
# --------------------------------------------------
sim.step()

[AlexGlock]

I solved the issue here is the running example in case someone is interested

from pywarpx import picmi
import numpy as np
import scipy as sp
import math

# --------------------------------------------------
# Input   0 o----------------o L 
#       emitter (Phi_0)   collector (Phi_L)
# --------------------------------------------------
average_time = True
diag_name = 'diag_uniform_0.2eV'
gap = 5e-3             # electrode gap [m]
phi_0 = 0.0            # cathode potential [V]
phi_L = 10.0           # anode potential [V]

nz = 256               # number of cells
t_end = 1.5e-7         # final simulation time [s]
cfl_multipl = 0.95     # run with dt = <X> * dt_CFL
max_steps = 25000      # total number of steps
intervals = 250        # export diagnostics every <X> steps

initial_energy = 0.2    # initial electron energy
# 0.2eV -> 0.1 beta | 5ev -> beta: 0.5 | 20eV -> beta 1.0

# --- Physical constants
q_e = picmi.constants.q_e               # electron charge [C]
m_e = picmi.constants.m_e               # electron mass [kg]
epsilon0 = picmi.constants.ep0          # vacuum permittivity [F/m]
clight = picmi.constants.c              # speed of light [m/s]

# --- Derived quantities
dz = gap / nz
gap_voltage = phi_L - phi_0
vz_0 = math.sqrt(initial_energy * q_e / m_e)  # thermal velocity
vzfinal = math.sqrt(2.0 * abs(gap_voltage) * q_e / m_e)
dt = cfl_multipl * dz / (2 * math.sqrt(vz_0**2 + 2 * q_e * phi_L / m_e))
beta = math.sqrt(m_e * vz_0**2 / (2 * q_e * phi_L))

# Child-Langmuir current density
J_CL = (4.0/9.0) * epsilon0 * math.sqrt(2 * abs(q_e)/m_e) * abs(gap_voltage)**(1.5) / gap**2
J_SCL = (beta + math.sqrt(beta**2 + 1))**3 * J_CL

# --------------------------------------------------
# Grid (1D)
# --------------------------------------------------
grid = picmi.Cartesian1DGrid(
    number_of_cells=[nz],
    lower_bound=[0.0],
    upper_bound=[gap],
    lower_boundary_conditions=['dirichlet'],
    upper_boundary_conditions=['dirichlet'],
    lower_boundary_conditions_particles=['absorbing'],
    upper_boundary_conditions_particles=['absorbing'],
    warpx_potential_lo_z=phi_0,
    warpx_potential_hi_z=phi_L
)


# --------------------------------------------------
# Electrostatic solver
# --------------------------------------------------
solver = picmi.ElectrostaticSolver(
    grid=grid,
    warpx_dt_update_interval=100,
    warpx_cfl=cfl_multipl
)


# --------------------------------------------------
# Particles
# --------------------------------------------------
electron_injector = picmi.AnalyticFluxDistribution(
    # inject at cathode surface z=0
    flux=0.995*(J_SCL/q_e),
    flux_normal_axis='z',
    surface_flux_position=0.0,
    flux_direction=1,
    directed_velocity= [0.,0.,vz_0]
)

electrons = picmi.Species(
    #particle_type='electron',
    name='electrons',
    charge=-q_e,
    mass=m_e,
    warpx_save_particles_at_zhi = True,
    warpx_save_particles_at_zlo = True,
    initial_distribution=electron_injector
)


# --------------------------------------------------
# Simulation
# --------------------------------------------------
sim = picmi.Simulation(
    solver=solver,
    max_steps=max_steps,
    particle_shape='linear',
    max_time=t_end,
    #warpx_used_inputs_file=diag_name+"/warpx_used_inputs"
)

sim.add_species(
    species=electrons,
    initialize_self_field=True, 
    layout=picmi.PseudoRandomLayout(
            n_macroparticles_per_cell=100,
            grid=grid
        )
    )

#sim.write_input_file(diag_name+"/input")

# --------------------------------------------------
# Diagnostics
# --------------------------------------------------
field_diag = picmi.FieldDiagnostic(
    name='fields',
    grid=grid,
    period=intervals,
    data_list=['Ez', 'rho', 'Jz', 'phi'],
    write_dir=diag_name,
    warpx_format='openpmd',
)

particle_diag = picmi.ParticleDiagnostic(
    name='fluxes',
    period=intervals,
    species=[electrons],
    data_list=['position','momentum'],
    write_dir=diag_name,
    warpx_format='openpmd',
)

sim.add_diagnostic(field_diag)
sim.add_diagnostic(particle_diag)

# --- Run simulation
sim.step()