LB: Particle not coupled due to rounding error

[RudolfWeeber]

Under veyr specific circumstances, particles are not fully coupled to LB:

Trigger: a particle coordinate is -epsilon but epsilon + box_l is = box_l due to a rounding error. In this case, neither the original position nor the folded positoin (-eps+box_l) is in the local domain and no force is added to the particle.

Probability of this occurring: probably very low in a system with randomly placed particles, but likely, if particles are initialized with a zero in one of the positoin coordinates. this can easily propagate to -epsilon due small numerical errors in the LB forces

MNWE (but this might be compiler dependent):

import espressomd

import espressomd.interactions
import espressomd.lb

import numpy as np

system = espressomd.System(box_l=[50,50,40])

system.time_step =0.01
system.cell_system.skin = 0.4

sigma = 1





if hasattr(espressomd.lb, "LBFluidWalberla"):
    lbf = espressomd.lb.LBFluidWalberla(
      agrid=1,
      density=1,
      tau=system.time_step,
      kinematic_viscosity=1)
else:
    lbf = espressomd.lb.LBFluid(
      agrid=1,
      dens=1,
      tau=system.time_step,
      visc=1)
if hasattr(system, "actors"): 
    system.actors.add(lbf)
else:
    system.lb = lbf
system.thermostat.set_lb(LB_fluid=lbf, gamma=1,seed=1, act_on_virtual=False)

p1 = system.part.add(pos=(0,-1E-25,0),v=(-1,0,0))
p2 = system.part.add(pos=(30,0,0),v=(1,0,0))

#system.thermostat.set_langevin(kT=0,gamma=1,seed=2)
for i in range(500):
  system.integrator.run(1)
  print(system.distance(p1,p2), p1.v[0],p2.v[0],p1.f[0],p2.f[0])
  np.testing.assert_allclose(np.copy(p1.f), -np.copy(p2.f), atol=1E-5)

Likely, both, the dev version and 4.2 are affected.

[RudolfWeeber]

4.2.1 is also affected

[RudolfWeeber]

As an aside: the correct term for the error here is cancellation, not rounding, I think.