LB GPU does not apply velocity bc under unclear circumstances

[RudolfWeeber]

Below script simulates a slab geometry with v=0 boundary conditions on top and bottom.
An external force density is applied to the fluid/
Also, a particle is advected.

The resulting flow profile is Poisseullew'ish for LB CPU.
For LB GPU, it is completeley flat unless the WCA interaction is
disabled. Then, it is poisseulle like as well.

import espressomd
import espressomd.lb
from espressomd.interactions import HarmonicBond
import espressomd.lbboundaries
import espressomd.shapes as shapes
from espressomd.observables import *
from espressomd.virtual_sites import VirtualSitesRelative

import sys
from time import time
import argparse
import cPickle
import gzip

import numpy as np

kT=0.
gamma=1.


# box
agrid=1.
l_xy=64
l_z=25

rasp_center_type=1
rasp_beed_type=2
wall_type=3


s=espressomd.System(box_l=(l_xy,l_xy,l_z))
s.periodicity=1,1,1

s.time_step=0.01






center =s.part.add(pos=s.box_l/2,type=rasp_center_type)

# Boundaries
wall_shape1 = espressomd.shapes.Wall(normal=[0, 0, 1], dist=agrid)
wall_shape2 = espressomd.shapes.Wall(
        normal=[0, 0, -1], dist=-(s.box_l[2] - agrid))

s.constraints.add(shape=wall_shape1,particle_type=wall_type)
s.constraints.add(shape=wall_shape2,particle_type=wall_type)

wall1 = espressomd.lbboundaries.LBBoundary(shape=wall_shape1,velocity=[0,0,0])
wall2 = espressomd.lbboundaries.LBBoundary(shape=wall_shape2,velocity=[0,0,0])
s.lbboundaries.add(wall1)
s.lbboundaries.add(wall2)


# Interactions
#rasp wall
s.non_bonded_inter[rasp_center_type,rasp_center_type].wca.set_params(sigma=4,epsilon=1)

# lb
lbf=espressomd.lb.LBFluidGPU(kT=0,seed=0,agrid=agrid,dens=5,visc=3,tau=s.time_step,ext_force_density=(1,0,0))
s.actors.add(lbf)
s.thermostat.set_lb(LB_fluid=lbf,gamma=20, seed=1,act_on_virtual=True)

s.cell_system.skin=0.3


lb_vel_profile=espressomd.observables.LBVelocityProfile(
        n_x_bins=1,
        n_y_bins=1,
        n_z_bins=int(s.box_l[2]),
        min_x=0,
        min_y=0,
        min_z=0,
        max_x=s.box_l[0],
        max_y=s.box_l[1],
        max_z=s.box_l[2],
        sampling_offset_x=0.*agrid,
        sampling_offset_y=0.*agrid,
        sampling_offset_z=0.*agrid,
        allow_empty_bins=True)

measurement_rounds=10000
measurement_interval=10

for i in range(measurement_rounds):
    tick=time()
    s.integrator.run(measurement_interval)
    tock=time()
    if i%100==0:
      print np.array(lb_vel_profile.calculate()).reshape((int(s.box_l[2]),3))
      print center.pos,center.v
      print

[KaiSzuttor]

your gamma is quite large... are you sure that the LB is stable? Also, we have a python test that checks a very similar system...

[RudolfWeeber]

setting a short range ia changes the cell geometry.
cells_re_init() calls lb_init_gpu() which calls lb_init_gpu() which deletetes the boundary flasg on the lbgpu nodes via the
reset_boundaries() kernel.

[KaiSzuttor]

remove the kernel call?

[RudolfWeeber]

This is a regression introduced after 4.0.
In 4.0 the code in on_cell_structure_change() was
#ifdef LB
if (lattice_switch & LATTICE_LB) {
lb_init();
}
#endif

now it reads:
#ifdef LB
lb_lbfluid_init();
#endif
So, 4.0.2 does not have to wait for this.
I'll open a pr restoring the old behavior.