diff --git a/lbmpy/forcemodels.py b/lbmpy/forcemodels.py index e8895124ca7c6a6fa191973f6dcaa22f5d797ab6..060604225ff784900938f71259ea9b48efee6724 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 028dda313fab5aa7a5176186606f814106393c03..6e0134001fada55418b60a2024520000fd072b19 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