diff --git a/doc/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb b/doc/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb
index 20e1653b9692f9488ad55a36113cf52d0e70be9c..6c2ff255e9826e28b1b5baebca6190c480a38996 100644
--- a/doc/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb
+++ b/doc/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb
@@ -419,7 +419,7 @@
],
"source": [
"method = create_lb_method(method='mrt', weighted=True, stencil='D2Q9', force=(1e-6, 0),\n",
- " relaxation_rates=[0, 0, ω, 1.9, 1.9])\n",
+ " force_model='luo', relaxation_rates=[0, 0, ω, 1.9, 1.9])\n",
"method"
]
},
diff --git a/doc/sphinx/lbmpy.bib b/doc/sphinx/lbmpy.bib
index 9620dc345ce71d66b4fca4719f64ea83f654732b..bceb010e1ee47c31a8109031f90491bcd5b818d0 100644
--- a/doc/sphinx/lbmpy.bib
+++ b/doc/sphinx/lbmpy.bib
@@ -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},
diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 4eca94b03d4db460eda7dfecd5d5e29555c909b8..4b3cc76760f9474031016d5feee9f172d585c415 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -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,))
diff --git a/lbmpy/forcemodels.py b/lbmpy/forcemodels.py
index 017927c37c8cb27f962497ab24770f35a11329de..ad57d27c9b95011386b11528cdf44d426985e155 100644
--- a/lbmpy/forcemodels.py
+++ b/lbmpy/forcemodels.py
@@ -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)]
diff --git a/lbmpy/relaxationrates.py b/lbmpy/relaxationrates.py
index 70c9bdd5c99149546d55c2e01cc28321183c520b..adcdfcb65bdf6f1aac690fdc610935abafb0c8f1 100644
--- a/lbmpy/relaxationrates.py
+++ b/lbmpy/relaxationrates.py
@@ -1,6 +1,6 @@
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.
diff --git a/lbmpy_tests/test_entropic.py b/lbmpy_tests/test_entropic.py
index 284570c1d4c874f8d11737a70346ade24ce94eba..47f6e841288e93f26037f5a38a29d2895f609de8 100644
--- a/lbmpy_tests/test_entropic.py
+++ b/lbmpy_tests/test_entropic.py
@@ -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)
diff --git a/lbmpy_tests/test_fluctuation.py b/lbmpy_tests/test_fluctuation.py
index 63078aa12ef5373ae5d350f17c6143894c9a5f33..ded0db9149789bc95af08829a7e68b171752252e 100644
--- a/lbmpy_tests/test_fluctuation.py
+++ b/lbmpy_tests/test_fluctuation.py
@@ -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)
diff --git a/lbmpy_tests/test_force.py b/lbmpy_tests/test_force.py
new file mode 100644
index 0000000000000000000000000000000000000000..79872b7b17d7927af3b3bfdbcd2724ebdcd6d2b7
--- /dev/null
+++ b/lbmpy_tests/test_force.py
@@ -0,0 +1,157 @@
+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))
diff --git a/lbmpy_tests/test_split_optimization.py b/lbmpy_tests/test_split_optimization.py
index 89e73133c487b37af3716a956911d7f46af2ae87..dd1093806086c5adc7875a72d69cc515bca4b282 100644
--- a/lbmpy_tests/test_split_optimization.py
+++ b/lbmpy_tests/test_split_optimization.py
@@ -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)
diff --git a/lbmpy_tests/test_srt_trt_simplifications.py b/lbmpy_tests/test_srt_trt_simplifications.py
index ad0a96247593198ceb2e0db1dd536b5d9f8c2b83..ca830c4d3ce5f36a1fdc459be6bd598194c3bdb8 100644
--- a/lbmpy_tests/test_srt_trt_simplifications.py
+++ b/lbmpy_tests/test_srt_trt_simplifications.py
@@ -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])