From 464159115e8a39c15d1c971e5ee37e479ca75754 Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Fri, 4 Sep 2020 19:31:29 +0200
Subject: [PATCH] fix force test

---
 lbmpy/forcemodels.py      | 2 +-
 lbmpy_tests/test_force.py | 8 ++++++--
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/lbmpy/forcemodels.py b/lbmpy/forcemodels.py
index e8895124..06060422 100644
--- a/lbmpy/forcemodels.py
+++ b/lbmpy/forcemodels.py
@@ -174,7 +174,7 @@ class Schiller:
         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, 3)) * sp.Rational(1, 2) * \
+        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)
 
         result = []
diff --git a/lbmpy_tests/test_force.py b/lbmpy_tests/test_force.py
index 028dda31..6e013400 100644
--- a/lbmpy_tests/test_force.py
+++ b/lbmpy_tests/test_force.py
@@ -2,6 +2,7 @@ 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
@@ -103,7 +104,7 @@ def test_stress(stencil):
     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
+    assert list(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
@@ -121,4 +122,7 @@ def test_stress(stencil):
                                     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
+        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
-- 
GitLab