Thermalized CPU LB produces unphysical currents

[christophlohrmann]

When using the thermalized CPU LB with boundaries that create corners, unphysical currents occur. See below an example. With the CPU implementation, there are currents in y and z direction without any external forces or velocity boundary conditions. For the GPU implementation, this does not happen.

from espressomd import System, lb, lbboundaries, shapes
import matplotlib.pyplot as plt
import numpy as np

impl = 'CPU' #'GPU'

box_l = [40,15,40]
LB_params = {'agrid':1,
             'tau':1,
             'dens':1,
             'visc':0.2,
             'kT':1e-8,
             'seed':41}
padding = 4
system = System(box_l = box_l)
system.cell_system.skin = 1
system.time_step = LB_params['tau']

if impl == 'GPU':
    lbf = lb.LBFluidGPU(**LB_params)
elif impl == 'CPU':
    lbf = lb.LBFluid(**LB_params)
    
system.actors.add(lbf)

top_wall = shapes.Wall(dist = -(box_l[2]-padding), normal = [0,0,-1])
bottom_wall = shapes.Wall(dist = padding, normal = [0,0,1])
left_wall = shapes.Wall(dist = padding, normal = [1,0,0])
right_wall = shapes.Wall(dist = -(box_l[0]-padding), normal = [-1,0,0])
for wall in [top_wall,bottom_wall, left_wall,right_wall]:
    noslip_wall = lbboundaries.LBBoundary(shape = wall)
    system.lbboundaries.add(noslip_wall)

system.integrator.run(100)

#lbf.print_vtk_velocity(f'./velocity_{impl}.vtk')  

probe_xs = np.linspace(0,box_l[0],num=40)
vels_over_time = []
n_samples = 100
for _ in range(n_samples):
    probe_vels = np.array([lbf.get_interpolated_velocity([x,0,box_l[2]/2.]) for x in probe_xs])
    vels_over_time.append(np.mean(probe_vels,axis = 0))
    system.integrator.run(10)
vels_over_time = np.array(vels_over_time)
mean_vel = np.mean(vels_over_time, axis = 0)
std_vel = np.std(vels_over_time, axis=0)
SEM_vel = std_vel/np.sqrt(n_samples)

print(f'velocity = {mean_vel} +- {SEM_vel}')
print(np.all(np.abs(mean_vel)<std_vel))
#I know we should compare to the standard error of the mean (divide by sqrt(n_samples)) but that takes forever to produce True on GPU. 
                          
#fig = plt.figure()
#for direction in [0,1,2]:
    #plt.plot(probe_xs, probe_vels[:,direction], label = f'v_{direction}')
#plt.xlabel('x')
#plt.legend()
#plt.show()