From 4c83f43512dd33eaf5bb6fab4829d772e8c4e254 Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Sat, 5 Sep 2020 12:40:00 +0200
Subject: [PATCH] Compare force term against conditions from Schiller
---
lbmpy/forcemodels.py | 2 +-
lbmpy_tests/test_force.py | 18 +++++++++++++++++-
2 files changed, 18 insertions(+), 2 deletions(-)
diff --git a/lbmpy/forcemodels.py b/lbmpy/forcemodels.py
index 06060422..87924699 100644
--- a/lbmpy/forcemodels.py
+++ b/lbmpy/forcemodels.py
@@ -175,7 +175,7 @@ class Schiller:
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, 3) * (2 + omega_bulk)
+ (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):
diff --git a/lbmpy_tests/test_force.py b/lbmpy_tests/test_force.py
index 4757064c..e7021aa4 100644
--- a/lbmpy_tests/test_force.py
+++ b/lbmpy_tests/test_force.py
@@ -120,7 +120,7 @@ def test_modes(stencil, force_model):
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)
+ sp.Rational(1,method.dim) * (2 + omega_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)}
@@ -133,6 +133,22 @@ def test_modes(stencil, force_model):
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 + omega_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+omega_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":
--
GitLab