From efb97424f9bca9f1713efa87657589df4efa9fc5 Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Fri, 4 Sep 2020 23:22:48 +0200
Subject: [PATCH] =?UTF-8?q?Test=20the=20simple=20force=20model=E2=80=99s?=
=?UTF-8?q?=20moments?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
lbmpy_tests/test_force.py | 58 +++++++++++++++++++++++----------------
1 file changed, 35 insertions(+), 23 deletions(-)
diff --git a/lbmpy_tests/test_force.py b/lbmpy_tests/test_force.py
index 6e013400..4757064c 100644
--- a/lbmpy_tests/test_force.py
+++ b/lbmpy_tests/test_force.py
@@ -82,7 +82,8 @@ def test_total_momentum(method, force_model, omega):
@pytest.mark.parametrize("stencil", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
-def test_stress(stencil):
+@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])
@@ -98,31 +99,42 @@ def test_stress(stencil):
stencil=stencil,
relaxation_rates=[omega_s, omega_b, omega_o, omega_e, omega_o, omega_e],
compressible=True,
- force_model="schiller",
+ force_model=force_model,
force=F)
force_moments = sp.simplify(method.moment_matrix * sp.Matrix(method.force_model(method)))
- # The momentum modes should contain the force
- assert list(force_moments[1:dim+1]) == F
+ # The mass mode should be zero
+ assert force_moments[0] == 0
- # 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)
+ # The momentum moments should contain the force
+ assert list(force_moments[1:dim+1]) == F
- 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)
- 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
+ if force_model == "schiller":
+ num_stresses = (dim*dim-dim)//2+dim
+
+ # 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 + 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)
+
+ 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
+
+ # 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))
--
GitLab