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 |
Edited by Michael Kuron