Schiller force model

lbmpy/creationfunctions.py

@@ -479,6 +479,7 @@ def force_model_from_string(force_model_name, force_values):
         'buick': forcemodels.Buick,
         'silva': forcemodels.Buick,
         'edm': forcemodels.EDM,
+        'schiller': forcemodels.Schiller,
     }
     if force_model_name.lower() not in force_model_dict:
         raise ValueError("Unknown force model %s" % (force_model_name,))

lbmpy/forcemodels.py

@@ -88,7 +88,7 @@ second moment nonzero   :class:`Luo`         :class:`Guo`
 
 import sympy as sp
 
-from lbmpy.relaxationrates import get_shear_relaxation_rate
+from lbmpy.relaxationrates import get_bulk_relaxation_rate, get_shear_relaxation_rate
 
 
 class Simple:
@@ -148,6 +148,7 @@ class Guo:
         luo = Luo(self._force)
 
         shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
+        assert len(set(lb_method.relaxation_rates)) == 1, "Guo only works for SRT, use Schiller instead"
         correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
         return [correction_factor * t for t in luo(lb_method)]
 
@@ -158,6 +159,35 @@ class Guo:
         return default_velocity_shift(density, self._force)
 
 
+class Schiller:
+    r"""
+    Force model by Schiller  :cite:`schiller2008thermal`
+    Equivalent to Guo but not restricted to SRT.
+    """
+    def __init__(self, force):
+        self._force = force
+
+    def __call__(self, lb_method):
+        u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)
+        force = sp.Matrix(self._force)
+        
+        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) * \
+            (2 + omega) + uf * sp.Rational(1, 3) * (2 + omega_bulk)
+
+        result = []
+        for direction, w_i in zip(lb_method.stencil, lb_method.weights):
+            direction = sp.Matrix(direction)
+            tr = sp.trace(G * (direction * direction.transpose() - sp.Rational(1, 3) * sp.eye(len(force))))
+            result.append(3 * w_i * (force.dot(direction) + sp.Rational(3, 2) * tr))
+        return result
+    
+    def macroscopic_velocity_shift(self, density):
+        return default_velocity_shift(density, self._force)
+
+
 class Buick:
     r"""
     This force model :cite:`buick2000gravity` has a force term with zero second moment. 

lbmpy/relaxationrates.py

 import sympy as sp
 
-from lbmpy.moments import is_shear_moment
+from lbmpy.moments import is_bulk_moment, is_shear_moment
 
 
 def relaxation_rate_from_lattice_viscosity(nu):
@@ -46,6 +46,18 @@ def get_shear_relaxation_rate(method):
                                       "Can not determine their relaxation rate automatically")
 
 
+def get_bulk_relaxation_rate(method):
+    """
+    The bulk moment is x^2 + y^2 + z^2, plus a constant for orthogonalization.
+    If this moment does not exist, the bulk relaxation is part of the shear relaxation.
+    The bulk relaxation rate determines the bulk viscosity in hydrodynamic LBM schemes.
+    """
+    for moment, relax_info in method.relaxation_info_dict.items():
+        if is_bulk_moment(moment, method.dim):
+            return relax_info.relaxation_rate
+    return get_shear_relaxation_rate(method)
+
+
 def relaxation_rate_scaling(omega, level_scale_factor):
     """Computes adapted omega for refinement.