Guo and Buick force is wrong when applied to non-SRT LBs

As a simple test, we apply a constant force to a fluid and measure the resulting total momentum. It should be force * number of cells * number of time steps. We do this for different relaxation and force models.

from pystencils.session import *
from lbmpy.session import *
import lbmpy
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter

L = (40, 40)
stencil = get_stencil("D2Q9")
eta = 0.15
omega = lbmpy.relaxationrates.relaxation_rate_from_lattice_viscosity(eta)
F = [2e-4, -3e-4]

dh = ps.create_data_handling(L, periodicity=True, default_target='cpu')
src = dh.add_array('src', values_per_cell=len(stencil))
dst = dh.add_array_like('dst', 'src')
ρ = dh.add_array('rho')
u = dh.add_array('u', values_per_cell=dh.dim)

collision = create_lb_update_rule(method="trt",
                                  stencil=stencil,
                                  relaxation_rate=omega, 
                                  compressible=True,
                                  force_model='guo', 
                                  force=F,
                                  kernel_type='collide_only',
                                  optimization={'symbolic_field': src})

stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
                                               {'density': ρ, 'velocity': u})

opts = {'cpu_openmp': True, 
        'cpu_vectorize_info': None,
        'target': dh.default_target}

stream_kernel = ps.create_kernel(stream, **opts).compile()
collision_kernel = ps.create_kernel(collision, **opts).compile()

def init():
    dh.fill(ρ.name, 1)
    dh.fill(u.name, 0)
    
    setter = macroscopic_values_setter(collision.method, velocity=(0,)*dh.dim, 
                                   pdfs=src.center_vector, density=ρ.center)
    kernel = ps.create_kernel(setter, ghost_layers=0).compile()
    dh.run_kernel(kernel)

sync_pdfs = dh.synchronization_function([src.name]) # needed before stream, but after collision

def time_loop(steps):
    dh.all_to_gpu()
    for i in range(steps):
        dh.run_kernel(collision_kernel)
        sync_pdfs()
        dh.run_kernel(stream_kernel)
        dh.swap(src.name, dst.name)
    dh.all_to_cpu()

t = 100
init()
time_loop(t)
total = np.sum(dh.gather_array(u.name), axis=(0,1))
print(total/np.prod(L)/F/t)

We see that for SRT or Luo, the script always prints 1.0, so all force applied has been converted into momentum. For TRT with Guo, it produces an omega-dependent number, which means that force has not been applied correctly.

Method	Momentum per Force
SRT Guo	1.0
SRT Luo	1.0
SRT Buick	1.0
TRT Guo	*
TRT Luo	1.0
TRT Buick	*

omega	*
2	0
1.8	0.22903226
1.5	0.55769231
1.2	0.87058824
1	1.07142857
0.8	1.26666667
0.5	1.55
0.25	1.77822581
0.1	1.92
0	2