From c3dbce3ca9ecca9fc69c9d2db85b62795623cab8 Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Fri, 4 Sep 2020 16:55:36 +0200
Subject: [PATCH] Default to Schiller force model and test against Ladd&Verberg
---
lbmpy/creationfunctions.py | 6 +++---
lbmpy/forcemodels.py | 4 ++--
lbmpy_tests/test_force.py | 44 ++++++++++++++++++++++++++++++++++++++
3 files changed, 49 insertions(+), 5 deletions(-)
diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 52c9bbbb..4b3cc767 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])
diff --git a/lbmpy/forcemodels.py b/lbmpy/forcemodels.py
index 5b4a0b19..e8895124 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.
@@ -161,7 +161,7 @@ class Guo:
class Schiller:
r"""
- Force model by Schiller :cite:`schiller2008thermal`
+ Force model by Schiller :cite:`schiller2008thermal`, equation 4.67
Equivalent to Guo but not restricted to SRT.
"""
def __init__(self, force):
diff --git a/lbmpy_tests/test_force.py b/lbmpy_tests/test_force.py
index b42c03eb..028dda31 100644
--- a/lbmpy_tests/test_force.py
+++ b/lbmpy_tests/test_force.py
@@ -78,3 +78,47 @@ def test_total_momentum(method, force_model, omega):
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"])
+def test_stress(stencil):
+ """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="schiller",
+ force=F)
+
+ force_moments = sp.simplify(method.moment_matrix * sp.Matrix(method.force_model(method)))
+
+ # The momentum modes should contain the force
+ assert force_moments[1:dim+1] == F
+
+ # The stress modes 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 + omega_s) * (traceless(sp.Matrix(u) * sp.Matrix(F).transpose()) + \
+ traceless(sp.Matrix(F) * sp.Matrix(u).transpose())) + \
+ sp.Rational(1,3) * (2 + omega_b) * sp.Matrix(u).dot(F) * sp.eye(method.dim)
+
+ num_stresses = (dim*dim-dim)//2+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)
+ assert diff == 0 or isinstance(diff, sp.Rational) # difference should be zero or a constant
--
GitLab