EK density continuity broken for for boundaries in the ghostlayer

[itischler]

I found a bug, where the global ek density increases over time. This happens when the sample script is run with multiple mpi processes, put it works fine when only using 1 process.

import numpy as np
import espressomd
import espressomd.electrokinetics
import espressomd.shapes

box_l = 20
time_step = 0.2

# Parameters
kT = 1.0
density = 0.2
valency = 1.0
diffusion = 0.2
ext_efield = [0.01, 0., 0.]
sphere_radius = 4.0


# Set-up the system
system = espressomd.System(box_l=[box_l] * 3)
system.time_step = time_step
system.cell_system.skin = 1.0


# Set-up the EK
lattice = espressomd.electrokinetics.LatticeWalberla(agrid=1.0, n_ghost_layers=1)
ek_solver = espressomd.electrokinetics.EKNone(lattice=lattice)
system.ekcontainer = espressomd.electrokinetics.EKContainer(
    solver=ek_solver, tau=system.time_step)


# Set-up the EK species
ek_species_pos = espressomd.electrokinetics.EKSpecies(
    lattice=lattice,
    single_precision=False,
    kT=kT,
    tau=system.time_step,
    density=density,
    valency=valency,
    diffusion=diffusion,
    advection=False,
    friction_coupling=False,
    ext_efield=ext_efield,
    thermalized=False
)
system.ekcontainer.add(ek_species_pos)

# Set-up the boundary
sphere = espressomd.shapes.Sphere(center=[box_l/2]*3, radius=sphere_radius)

ek_species_pos.add_boundary_from_shape(
    shape=sphere, value=[0.0,0.0,0.0], boundary_type=espressomd.electrokinetics.FluxBoundary)
    

for i in range(20):
    system.integrator.run(100)
    print(np.mean(ek_species_pos[:,:,:].density))