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!