Schiller force model to be used instead of Guo on non-SRT

This is supposed to be equivalent to Guo, but confirmed to be working for non-SRT. Also includes two test cases:

• one that confirms momentum conservation for all available force models
• one that checks that the moments of the force look right, but only for Schiller and Simple force models

Fixes #11 (closed)

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 force_moments[dim+1+num_stresses:] == [0] * (len(stencil)-(dim+1+num_stresses))
+    elif force_model == "simple":
+        # All other moments should be zero
+        assert force_moments[dim+1:] == [0] * (len(stencil)-(dim+1))