Commit 235979c0 authored by Markus Holzer's avatar Markus Holzer
Browse files

Merge branch 'schiller' into 'master'

Schiller force model to be used instead of Guo on non-SRT

Closes #11

See merge request pycodegen/lbmpy!34
parents 0153b244 47e9a1a9
......@@ -242,11 +242,11 @@
%% Cell type:code id: tags:
``` python
method = create_lb_method(method='mrt', weighted=True, stencil='D2Q9', force=(1e-6, 0),
relaxation_rates=[0, 0, ω, 1.9, 1.9])
force_model='luo', relaxation_rates=[0, 0, ω, 1.9, 1.9])
method
```
%%%% Output: execute_result
......
......@@ -41,6 +41,13 @@
publisher={APS}
}
@phdthesis{schiller2008thermal,
title={Thermal fluctuations and boundary conditions in the lattice Boltzmann method},
author={Schiller, Ulf Daniel},
year={2008},
school={Johannes Gutenberg Universit{\"a}t Mainz}
}
@article{Wohrwag2018,
archivePrefix = {arXiv},
......
......@@ -41,8 +41,8 @@ General:
compressible.
- ``equilibrium_order=2``: order in velocity, at which the equilibrium moment/cumulant approximation is
truncated. Order 2 is sufficient to approximate Navier-Stokes
- ``force_model=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'`` ``'buick'``, or an instance of a
class implementing the same methods as the classes in :mod:`lbmpy.forcemodels`. For details, see
- ``force_model=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'``, ``'buick'``, ``'schiller'``, or
an instance of a class implementing the same methods as the classes in :mod:`lbmpy.forcemodels`. For details, see
:mod:`lbmpy.forcemodels`
- ``force=(0,0,0)``: either constant force or a symbolic expression depending on field value
- ``maxwellian_moments=True``: way to compute equilibrium moments/cumulants, if False the standard
......@@ -387,7 +387,7 @@ def create_lb_method(**params):
no_force_model = 'force_model' not in params or params['force_model'] == 'none' or params['force_model'] is None
if not force_is_zero and no_force_model:
params['force_model'] = 'guo'
params['force_model'] = 'schiller'
if 'force_model' in params:
force_model = force_model_from_string(params['force_model'], params['force'][:dim])
......@@ -479,6 +479,7 @@ def force_model_from_string(force_model_name, force_values):
'buick': forcemodels.Buick,
'silva': forcemodels.Buick,
'edm': forcemodels.EDM,
'schiller': forcemodels.Schiller,
}
if force_model_name.lower() not in force_model_dict:
raise ValueError("Unknown force model %s" % (force_model_name,))
......
......@@ -7,7 +7,7 @@ Get started:
------------
This module offers different models to introduce a body force in the lattice Boltzmann scheme.
If you don't know which model to choose, use :class:`lbmpy.forcemodels.Guo`.
If you don't know which model to choose, use :class:`lbmpy.forcemodels.Schiller`.
For incompressible collision models the :class:`lbmpy.forcemodels.Buick` model can be better.
......@@ -88,7 +88,7 @@ second moment nonzero :class:`Luo` :class:`Guo`
import sympy as sp
from lbmpy.relaxationrates import get_shear_relaxation_rate
from lbmpy.relaxationrates import get_bulk_relaxation_rate, get_shear_relaxation_rate
class Simple:
......@@ -148,6 +148,7 @@ class Guo:
luo = Luo(self._force)
shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
assert len(set(lb_method.relaxation_rates)) == 1, "Guo only works for SRT, use Schiller instead"
correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
return [correction_factor * t for t in luo(lb_method)]
......@@ -158,6 +159,35 @@ class Guo:
return default_velocity_shift(density, self._force)
class Schiller:
r"""
Force model by Schiller :cite:`schiller2008thermal`, equation 4.67
Equivalent to Guo but not restricted to SRT.
"""
def __init__(self, force):
self._force = force
def __call__(self, lb_method):
u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)
force = sp.Matrix(self._force)
uf = u.dot(force) * sp.eye(len(force))
omega = get_shear_relaxation_rate(lb_method)
omega_bulk = get_bulk_relaxation_rate(lb_method)
G = (u * force.transpose() + force * u.transpose() - uf * sp.Rational(2, lb_method.dim)) * sp.Rational(1, 2) * \
(2 - omega) + uf * sp.Rational(1, lb_method.dim) * (2 - omega_bulk)
result = []
for direction, w_i in zip(lb_method.stencil, lb_method.weights):
direction = sp.Matrix(direction)
tr = sp.trace(G * (direction * direction.transpose() - sp.Rational(1, 3) * sp.eye(len(force))))
result.append(3 * w_i * (force.dot(direction) + sp.Rational(3, 2) * tr))
return result
def macroscopic_velocity_shift(self, density):
return default_velocity_shift(density, self._force)
class Buick:
r"""
This force model :cite:`buick2000gravity` has a force term with zero second moment.
......@@ -171,6 +201,7 @@ class Buick:
simple = Simple(self._force)
shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
assert len(set(lb_method.relaxation_rates)) == 1, "Buick only works for SRT"
correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
return [correction_factor * t for t in simple(lb_method)]
......
import sympy as sp
from lbmpy.moments import is_shear_moment
from lbmpy.moments import is_bulk_moment, is_shear_moment
def relaxation_rate_from_lattice_viscosity(nu):
......@@ -46,6 +46,18 @@ def get_shear_relaxation_rate(method):
"Can not determine their relaxation rate automatically")
def get_bulk_relaxation_rate(method):
"""
The bulk moment is x^2 + y^2 + z^2, plus a constant for orthogonalization.
If this moment does not exist, the bulk relaxation is part of the shear relaxation.
The bulk relaxation rate determines the bulk viscosity in hydrodynamic LBM schemes.
"""
for moment, relax_info in method.relaxation_info_dict.items():
if is_bulk_moment(moment, method.dim):
return relax_info.relaxation_rate
return get_shear_relaxation_rate(method)
def relaxation_rate_scaling(omega, level_scale_factor):
"""Computes adapted omega for refinement.
......
......@@ -10,13 +10,13 @@ from lbmpy.stencils import get_stencil
def test_entropic_methods():
sc_kbc = create_lid_driven_cavity((20, 20), method='trt-kbc-n4', relaxation_rate=1.9999,
entropic_newton_iterations=3, entropic=True, compressible=True,
force=(-1e-10, 0))
force=(-1e-10, 0), force_model="luo")
sc_srt = create_lid_driven_cavity((40, 40), relaxation_rate=1.9999, lid_velocity=0.05, compressible=True,
force=(-1e-10, 0))
force=(-1e-10, 0), force_model="luo")
sc_entropic = create_lid_driven_cavity((40, 40), method='entropic-srt', relaxation_rate=1.9999,
lid_velocity=0.05, compressible=True, force=(-1e-10, 0))
lid_velocity=0.05, compressible=True, force=(-1e-10, 0), force_model="luo")
sc_srt.run(1000)
sc_kbc.run(1000)
......
......@@ -5,7 +5,7 @@ from lbmpy.scenarios import create_channel
def test_fluctuating_generation_pipeline():
ch = create_channel((10, 10), stencil='D2Q9', method='mrt', weighted=True, relaxation_rates=[1.5] * 5, force=1e-5,
fluctuating={'temperature': 1e-9}, kernel_params={'time_step': 1, 'seed': 312},
force_model='luo', fluctuating={'temperature': 1e-9}, kernel_params={'time_step': 1, 'seed': 312},
optimization={'cse_global': True})
ch.run(10)
......
from pystencils.session import *
from lbmpy.session import *
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
import lbmpy.forcemodels
from lbmpy.moments import is_bulk_moment
import pytest
from contextlib import ExitStack as does_not_raise
force_models = [fm.lower() for fm in dir(lbmpy.forcemodels) if fm[0].isupper()]
def test_force_models_available():
assert 'guo' in force_models
assert 'luo' in force_models
@pytest.mark.parametrize("method", ["srt", "trt"])
@pytest.mark.parametrize("force_model", force_models)
@pytest.mark.parametrize("omega", [0.5, 1.5])
def test_total_momentum(method, force_model, omega):
L = (16, 16)
stencil = get_stencil("D2Q9")
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)
expectation = does_not_raise()
skip = False
if force_model in ['guo', 'buick'] and method != 'srt':
expectation = pytest.raises(AssertionError)
skip = True
with expectation:
collision = create_lb_update_rule(method=method,
stencil=stencil,
relaxation_rate=omega,
compressible=True,
force_model=force_model,
force=F,
kernel_type='collide_only',
optimization={'symbolic_field': src})
if skip:
return
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])
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 = 20
init()
time_loop(t)
total = np.sum(dh.gather_array(u.name), axis=(0,1))
assert np.allclose(total/np.prod(L)/F/t, 1)
@pytest.mark.parametrize("stencil", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
@pytest.mark.parametrize("force_model", ["simple", "schiller"])
def test_modes(stencil, force_model):
"""check Schiller's force term in mode space"""
stencil = get_stencil(stencil)
dim = len(stencil[0])
omega_s = sp.Symbol("omega_s")
omega_b = sp.Symbol("omega_b")
omega_o = sp.Symbol("omega_o")
omega_e = sp.Symbol("omega_e")
F = [sp.Symbol(f"F_{i}") for i in range(dim)]
method = create_lb_method(method="mrt", weighted=True,
stencil=stencil,
relaxation_rates=[omega_s, omega_b, omega_o, omega_e, omega_o, omega_e],
compressible=True,
force_model=force_model,
force=F)
force_moments = sp.simplify(method.moment_matrix * sp.Matrix(method.force_model(method)))
# The mass mode should be zero
assert force_moments[0] == 0
# The momentum moments should contain the force
assert list(force_moments[1:dim+1]) == F
if force_model == "schiller":
num_stresses = (dim*dim-dim)//2+dim
lambda_s, lambda_b = -omega_s, -omega_b
# The stress moments should match eq. 47 from https://doi.org/10.1023/A:1010414013942
u = method.first_order_equilibrium_moment_symbols
def traceless(m):
tr = sp.simplify(sp.Trace(m))
return m - tr/m.shape[0]*sp.eye(m.shape[0])
C = sp.Rational(1,2) * (2 + lambda_s) * (traceless(sp.Matrix(u) * sp.Matrix(F).transpose()) + \
traceless(sp.Matrix(F) * sp.Matrix(u).transpose())) + \
sp.Rational(1,method.dim) * (2 + lambda_b) * sp.Matrix(u).dot(F) * sp.eye(method.dim)
subs = {sp.Symbol(chr(ord("x")+i)) * sp.Symbol(chr(ord("x")+j)) : C[i,j]
for i in range(dim) for j in range(dim)}
for force_moment, moment in zip(force_moments[dim+1:dim+1+num_stresses],
method.moments[dim+1:dim+1+num_stresses]):
ref = moment.subs(subs)
diff = sp.simplify(ref - force_moment)
if is_bulk_moment(moment, dim):
assert diff == 0 or isinstance(diff, sp.Rational) # difference should be zero or a constant
else:
assert diff == 0 # difference should be zero
ff = sp.Matrix(method.force_model(method))
# Check eq. 4.53a from schiller2008thermal
assert sp.simplify(sum(ff)) == 0
# Check eq. 4.53b from schiller2008thermal
assert [sp.simplify(sum(ff[i] * stencil[i][j] for i in range(len(stencil)))) for j in range(dim)] == F
# Check eq. 4.61a from schiller2008thermal
ref = (2 + lambda_s)/2 * (traceless(sp.Matrix(u) * sp.Matrix(F).transpose()) + \
traceless(sp.Matrix(F) * sp.Matrix(u).transpose()))
s = sp.zeros(dim)
for i in range(0, len(stencil)):
s += ff[i] * traceless(sp.Matrix(stencil[i]) * sp.Matrix(stencil[i]).transpose())
assert sp.simplify(s-ref) == sp.zeros(dim)
# Check eq. 4.61b from schiller2008thermal
assert sp.simplify(sum(ff[i] * stencil[i][a]**2 for i in range(len(stencil)) for a in range(dim))
- (2+lambda_b)*sp.Matrix(u).dot(F)) == 0
# All other moments should be zero
assert list(force_moments[dim+1+num_stresses:]) == [0] * (len(stencil)-(dim+1+num_stresses))
elif force_model == "simple":
# All other moments should be zero
assert list(force_moments[dim+1:]) == [0] * (len(stencil)-(dim+1))
......@@ -17,6 +17,7 @@ def test_split_number_of_operations():
common_params = {'stencil': stencil,
'method': method,
'compressible': compressible,
'force_model': 'luo',
'force': (1e-6, 1e-5, 1e-7)
}
ast_with_splitting = create_lb_ast(optimization={'split': True}, **common_params)
......
......@@ -4,7 +4,7 @@ known acceptable values.
"""
import sympy as sp
from lbmpy.forcemodels import Guo
from lbmpy.forcemodels import Luo
from lbmpy.methods import create_srt, create_trt, create_trt_with_magic_number
from lbmpy.methods.momentbasedsimplifications import cse_in_opposing_directions
from lbmpy.simplificationfactory import create_simplification_strategy
......@@ -55,13 +55,13 @@ def test_simplifications_trt_d2q9_compressible():
def test_simplifications_trt_d3q19_force_incompressible():
o1, o2 = sp.symbols("omega_1 omega_2")
force_model = Guo([sp.Rational(1, 3), sp.Rational(1, 2), sp.Rational(1, 5)])
force_model = Luo([sp.Rational(1, 3), sp.Rational(1, 2), sp.Rational(1, 5)])
method = create_trt(get_stencil("D3Q19"), o1, o2, compressible=False, force_model=force_model)
check_method(method, [268, 281, 0], [241, 175, 1])
def test_simplifications_trt_d3q19_force_compressible():
o1, o2 = sp.symbols("omega_1 omega_2")
force_model = Guo([sp.Rational(1, 3), sp.Rational(1, 2), sp.Rational(1, 5)])
force_model = Luo([sp.Rational(1, 3), sp.Rational(1, 2), sp.Rational(1, 5)])
method = create_trt_with_magic_number(get_stencil("D3Q19"), o1, compressible=False, force_model=force_model)
check_method(method, [270, 284, 1], [243, 178, 1])
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment