Segfault when using uniform distributions with space-charge effects

[DrakeWhu]

Hi WarpX team.

I am having problems with a simulation trying to model plasmonic oscillations of a graphene bilayer. I am trying to add space-charge effects but I keep getting segfaults. This is similar to an already closed issue #5678 but in my case it seems to run deeper. I tried changing every analytic distribution to a uniform one but I am still having issues. The script is this one:

from pywarpx import picmi
import numpy as np

# Constantes físicas
c = picmi.constants.c
q_e = picmi.constants.q_e

# Número de pasos de la simulación
max_steps = 5000

# =======================================================
# DOMINIO Y MALLADO (ESCALA NANOMÉTRICA)
# =======================================================
nx = 256
ny = 256
nz = 256

# Dominio físico reducido: 100 nm en x, y y z.
xmin = -50e-9    # -50 nm
xmax =  50e-9    # 50 nm
ymin = -50e-9    # -50 nm
ymax =  50e-9    # 50 nm
zmin = 0         # 0 nm
zmax = 200e-9    # 100 nm

beta = 0.25

grid = picmi.Cartesian3DGrid(
    number_of_cells=[nx, ny, nz],
    lower_bound=[xmin, ymin, zmin],
    upper_bound=[xmax, ymax, zmax],
    lower_boundary_conditions=["neumann", "neumann", "neumann"],
    upper_boundary_conditions=["neumann", "neumann", "neumann"],
    lower_boundary_conditions_particles=["open", "open", "open"],
    upper_boundary_conditions_particles=["open", "open", "open"],
    moving_window_velocity=[0.0, 0.0, beta * c],
    warpx_start_moving_window_step=0,
    warpx_max_grid_size=64,
    warpx_blocking_factor=32,
    pml_cells=[10, 10, 10]
)

# =======================================================
# PARÁMETROS DEL PLASMA Y CAPAS (ESCALA NANOMÉTRICA)
# =======================================================
plasma_surface_density = 1.53e20  # [m^-2] valor apropiado para el plasma

# Ahora definimos el grosor de cada capa de 3 nm y la separación entre capas de 50 nm.
layer_thickness = 3e-09
layer_separation = 5e-08

plasma_density = plasma_surface_density / (layer_thickness)  # [m^-3] valor apropiado para el plasma

n_layers = 2
# Las 2 capas se colocan centradas en y = ± (layer_separation/2)
y_centers = [-layer_separation/2, layer_separation/2]  # [-25e-9, 25e-9]

uniform_distribution_top = picmi.UniformDistribution(
    density=plasma_density / 4,
    lower_bound=[xmin, y_centers[0] - layer_thickness/2, zmin],
    upper_bound=[xmax, y_centers[0] + layer_thickness/2, zmax],
    fill_in=True
)

uniform_distribution_bottom = picmi.UniformDistribution(
    density=plasma_density / 4,
    lower_bound=[xmin, y_centers[1] - layer_thickness/2, zmin],
    upper_bound=[xmax, y_centers[1] + layer_thickness/2, zmax],
    fill_in=True
)

uniform_distribution_electrons = picmi.UniformDistribution(
    density=plasma_density,
    lower_bound=[xmin, y_centers[0] - layer_thickness/2 - 2e-9, zmin],
    upper_bound=[xmax, y_centers[1] + layer_thickness/2 + 2e-9, zmax],
    rms_velocity=[95361] * 3,
    fill_in=True
)

carbon_ions_top = picmi.Species(
    particle_type="C",
    name="carbon_ions_top",
    charge_state=4,
    initial_distribution=uniform_distribution_top
)

carbon_ions_bottom = picmi.Species(
    particle_type="C",
    name="carbon_ions_bottom",
    charge_state=4,
    initial_distribution=uniform_distribution_bottom
)

electrons = picmi.Species(
    particle_type="electron",
    name="electrons",
    initial_distribution=uniform_distribution_electrons
)
# =======================================================
# PARÁMETROS DEL BEAM (CON DIMENSIONES RESOLUCIONABLES)
# =======================================================

# Ahora definimos un beam con dimensiones de 3 nm en x, y y z
beam_size = 3e-9       # 3 nm en x e y
beam_thickness = 3e-9    # 3 nm en z

# Carga total deseada: la de un protón (1.6e-19 C).
q_tot = 1.6e-19  # [C]

# Calcular el volumen del beam.
beam_volume = beam_size * beam_size * beam_thickness  # [m^3]

# Definir la densidad tal que, al integrar sobre el volumen del beam, se obtenga 1 partícula física.
beam_density = 1.0 / beam_volume   # [partículas/m^3]

beam_distribution = picmi.GaussianBunchDistribution(
    n_physical_particles = q_tot / q_e,
    rms_bunch_size=[beam_size / 2, beam_size / 2, beam_thickness / 2],
    centroid_position=[0.0, 0.0, zmax - beam_thickness*2],
    centroid_velocity=[0.0, 0.0, beta * c],
    rms_velocity=[0, 0, 0]
)

# Usamos la especie "proton" para que la carga elemental sea 1.6e-19 C (positiva)
beam = picmi.Species(
    particle_type="proton",
    name="beam",
    mass = 1,
    initial_distribution=beam_distribution
)

# =======================================================
# SOLVER, DIAGNÓSTICOS Y CONFIGURACIÓN DE LA SIMULACIÓN
# =======================================================
solver = picmi.ElectromagneticSolver(
    grid=grid, method="Yee", cfl=1.0, divE_cleaning=0
)

diag_field_list = ["B", "E", "J", "rho"]

particle_diag = picmi.ParticleDiagnostic(
    name="diag1",
    period=5000,
    write_dir=".",
    warpx_file_prefix="3D",
    warpx_format="openpmd",
    warpx_openpmd_backend="h5"
)

field_diag = picmi.FieldDiagnostic(
    name="diag1",
    grid=grid,
    period=5000,
    data_list=diag_field_list,
    write_dir=".",
    warpx_file_prefix="3D",
    warpx_format="openpmd",
    warpx_openpmd_backend="h5"
)

sim = picmi.Simulation(
    solver=solver,
    max_steps=max_steps,
    verbose=1,
    particle_shape="cubic",
    warpx_use_filter=0,
    warpx_serialize_initial_conditions=1,
    warpx_do_dynamic_scheduling=0
)

# Se añade la especie del plasma (Carbono) con su layout
sim.add_species(
    carbon_ions_top,
    initialize_self_field=True,
    layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[1, 1, 1])
)

sim.add_species(
    carbon_ions_bottom,
    initialize_self_field=True,
    layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[1, 1, 1])
)

sim.add_species(
    electrons,
    layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[1, 1, 1])
)

# Se añade el beam con el layout definido
sim.add_species(beam, layout=picmi.PseudoRandomLayout(grid=grid, n_macroparticles=1000))

# Se añaden los diagnósticos
sim.add_diagnostic(particle_diag)
sim.add_diagnostic(field_diag)

# Se escribe el input file, se inicializan los inputs y se corre la simulación
sim.write_input_file(file_name="inputs_3d_picmi")
sim.initialize_inputs()
sim.initialize_warpx()
sim.step(max_steps)

The simulation log shows that it starts writing into the first dump, but the moment it deposits the currents, it has the segmentation fault. The backtrace shows that the problem is on the current deposition indeed:

===== TinyProfilers ======
REG::WarpX::Evolve()
WarpX::Evolve()
WarpX::Evolve::step
WarpX::OneStep_nosub()
PhysicalParticleContainer::Evolve()
WarpXParticleContainer::DepositCurrent::CurrentDeposition

Turning the initialize_self_field option off makes the segfault go away.

Any idea what may be happening? I thought that not using analytic expressions for the plasmas would solve the issue but it is not working and I am a bit puzzled right now.

On the last issue, I said it worked for me because the condition didn't produce a plasma. I changed the 'and' to an 'or' and the ions where properly set but the segfault reapeared.

If you need any extra info I can share it too.

[RemiLehe]

Thanks for reporting this.

Could you answer a few quick questions:

• Is the initial plasma expected to be charge neutral?
• What is the typical amplitude of the E field that you get after the initial Poisson solve (when using initialize_self_field)
• Would you be able to replace this line:
  https://github.com/ECP-WarpX/WarpX/blob/development/Source/Particles/WarpXParticleContainer.cpp#L407
  with AMREX_ALWAYS_ASSERT_WITH_MESSAGE(, and then recompile the code, and see if you get a clearer error message?

Thanks!

[DrakeWhu]

Hi @RemiLehe

The plasma is indeed expected to be charge neutral, as it is supposed to be two graphene layers without any doping. I ran the script again without recompiling in order to check the physical amplitude of the E field. Not sure where I can see it but I'm gonna share the simulation.log, warpx_used_inputs and inputs_3d_picmi, as well as one Backtrace file (I usually simulate using 32 parallel CPU cores so I always get 32 Backtraces). Here are those files:

simulation.log

inputs_3d_picmi.txt

warpx_used_inputs.txt

Backtrace.0.0.txt

I am going to attempt a new compìlation with the change you suggested. It may take a while to compile but the moment I can simulate I'll share the log too

[DrakeWhu]

Hi again,

I recompiled WarpX with the change @RemiLehe suggested. I ran the same script and segfault happened again as expected. The error message is way clearer now. It can be seen on the simulation.log:

simulation_debug_run.log

It looks like the particle shape does not fit within the tile (I am using CPU) used for current deposition. Not sure about what that means but you guys will probably be able to help me.

Thanks!

EDIT 1: I tried running the simulation with quadratic shape just in case it worked, but it still crashed. I also tried linear and it still had the segfault.

[RemiLehe]

Thanks for doing these tests and confirming that the plasma should be neutral.

Looking at your input script, it seems that the plasma is not initially neutral:

• positive ions are initialized between -40 nm and -10 nm, and then between +10 nm and +40 nm
• electrons are initialized between -40 nm and +40 nm, meaning there is a negatively charged volume between -10 nm and 10 nm.

Could you try to correct this and re-run the simulation?