lbmpy/phasefield/simplex_projection.pyx → src/lbmpy/phasefield/simplex_projection.pyx

-# Workaround for cython bug
-# see https://stackoverflow.com/questions/8024805/cython-compiled-c-extension-importerror-dynamic-module-does-not-define-init-fu
-WORKAROUND = "Something"
-
+# cython: language_level=3str
 import cython
 
 

src/lbmpy/phasefield_allen_cahn/analytical.py

+import numpy as np
+from math import sinh, cosh, cos, sin, pi
+
+
+def analytic_rising_speed(gravitational_acceleration, bubble_diameter, viscosity_gas):
+    r"""
+    Calculated the analytical rising speed of a bubble. This is the expected end rising speed.
+    Args:
+        gravitational_acceleration: the gravitational acceleration acting in the simulation scenario. Usually it gets
+                                    calculated based on dimensionless parameters which describe the scenario
+        bubble_diameter: the diameter of the bubble at the beginning of the simulation
+        viscosity_gas: the viscosity of the fluid inside the bubble
+    """
+    result = -(gravitational_acceleration * bubble_diameter * bubble_diameter) / (12.0 * viscosity_gas)
+    return result
+
+
+def analytical_solution_microchannel(reference_length, length_x, length_y,
+                                     kappa_top, kappa_bottom,
+                                     t_h, t_c, t_0,
+                                     reference_surface_tension, dynamic_viscosity_light_phase,
+                                     transpose=True):
+    """
+    https://www.sciencedirect.com/science/article/pii/S0021999113005986
+    """
+    l_ref = reference_length
+    sigma_t = reference_surface_tension
+
+    kstar = kappa_top / kappa_bottom
+    mp = (l_ref // 2) - 1
+
+    w = pi / l_ref
+    a = mp * w
+    b = mp * w
+
+    f = 1.0 / (kstar * sinh(b) * cosh(a) + sinh(a) * cosh(b))
+    g = sinh(a) * f
+
+    h = (sinh(a) ** 2 - a ** 2) * (sinh(b) ** 2 - b ** 2) / \
+        ((sinh(b) ** 2 - b ** 2) * (sinh(2.0 * a) - 2.0 * a)
+         + (sinh(a) ** 2 - a ** 2) * (sinh(2.0 * b) - 2.0 * b))
+
+    Ca1 = sinh(a) ** 2 / (sinh(a) ** 2 - a ** 2)
+    Ca2 = -1.0 * mp * a / (sinh(a) ** 2 - a ** 2)
+    Ca3 = (2 * a - sinh(2 * a)) / (2.0 * (sinh(a) ** 2 - a ** 2))
+
+    Cb1 = sinh(b) ** 2 / (sinh(b) ** 2 - b ** 2)
+    Cb2 = -1.0 * mp * b / (sinh(b) ** 2 - b ** 2)
+    Cb3 = (-2 * b + sinh(2 * b)) / (2.0 * (sinh(b) ** 2 - b ** 2))
+
+    umax = -1.0 * (t_0 * sigma_t / dynamic_viscosity_light_phase) * g * h
+    jj = 0
+    xx = np.linspace(-l_ref - 0.5, l_ref - 0.5, length_x)
+    yy = np.linspace(-mp, mp, length_y)
+    u_x = np.zeros([len(xx), len(yy)])
+    u_y = np.zeros([len(xx), len(yy)])
+    t_a = np.zeros([len(xx), len(yy)])
+    tt = t_c - t_h
+    nom = kstar * t_c * mp + t_h * mp
+    denom = mp + kstar * mp
+    for y in yy:
+        ii = 0
+        for x in xx:
+            swx = sin(w * x)
+            cwx = cos(w * x)
+
+            if y > 0:
+                tmp1 = ((Ca1 + w * (Ca2 + Ca3 * y)) * cosh(w * y) + (Ca3 + w * Ca1 * y) * sinh(w * y))
+                tmp2 = (Ca1 * y * cosh(w * y) + (Ca2 + Ca3 * y) * sinh(w * y))
+
+                t_a[ii, jj] = (tt * y + nom) / denom + t_0 * f * sinh(a - y * w) * cwx
+                u_x[ii, jj] = umax * tmp1 * swx
+                u_y[ii, jj] = -w * umax * tmp2 * cwx
+
+            elif y <= 0:
+                tmp3 = (sinh(a) * cosh(w * y) - kstar * sinh(w * y) * cosh(a))
+                tmp4 = ((Cb1 + w * (Cb2 + Cb3 * y)) * cosh(w * y) + (Cb3 + w * Cb1 * y) * sinh(w * y))
+
+                t_a[ii, jj] = (kstar * tt * y + nom) / denom + t_0 * f * tmp3 * cwx
+                u_x[ii, jj] = umax * tmp4 * swx
+                u_y[ii, jj] = -w * umax * (Cb1 * y * cosh(w * y) + (Cb2 + Cb3 * y) * sinh(w * y)) * cwx
+
+            ii += 1
+        jj += 1
+    x, y = np.meshgrid(xx, yy)
+    if transpose:
+        return x, y, u_x.T, u_y.T, t_a.T
+    else:
+        return x, y, u_x, u_y, t_a

src/lbmpy/phasefield_allen_cahn/contact_angle.py

+import math
+import sympy as sp
+
+from pystencils.astnodes import Block, Conditional, SympyAssignment
+
+from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
+from pystencils.boundaries.boundaryconditions import Boundary
+
+from pystencils.typing import TypedSymbol
+from pystencils.typing import CastFunc
+
+
+class ContactAngle(Boundary):
+    r"""
+    Wettability condition on solid boundaries according to equation 25 in :cite:`Fakhari2018`.
+
+    Args:
+        contact_angle: contact angle in degrees which is applied between the fluid and the solid boundary.
+        interface_width: interface width of the phase field model.
+        name: optional name of the boundary
+        data_type: data type for temporary variables which are used.
+    """
+
+    inner_or_boundary = False
+    single_link = True
+
+    def __init__(self, contact_angle, interface_width, name=None, data_type='double'):
+        self._contact_angle = contact_angle
+        self._interface_width = interface_width
+        self._data_type = data_type
+
+        super(ContactAngle, self).__init__(name)
+
+    def __call__(self, field, direction_symbol, **kwargs):
+
+        neighbor = BoundaryOffsetInfo.offset_from_dir(direction_symbol, field.spatial_dimensions)
+        dist = TypedSymbol("h", self._data_type)
+        angle = TypedSymbol("a", self._data_type)
+        d = CastFunc(sum([x * x for x in neighbor]), self._data_type)
+
+        var = - dist * (4.0 / self._interface_width) * angle
+        tmp = 1 + var
+        else_branch = (tmp - sp.sqrt(tmp * tmp - 4.0 * var * field[neighbor])) / var - field[neighbor]
+        if field.index_dimensions == 0:
+            if isinstance(self._contact_angle, (int, float)):
+                result = [SympyAssignment(angle, math.cos(math.radians(self._contact_angle))),
+                          SympyAssignment(dist, 0.5 * sp.sqrt(d)),
+                          Conditional(sp.LessThan(var * var, 0.000001),
+                                      Block([SympyAssignment(field.center, field[neighbor])]),
+                                      Block([SympyAssignment(field.center, else_branch)]))]
+                return result
+            else:
+                result = [SympyAssignment(angle, sp.cos(self._contact_angle * (sp.pi / sp.Number(180)))),
+                          SympyAssignment(dist, 0.5 * sp.sqrt(d)),
+                          Conditional(sp.LessThan(var * var, 0.000001),
+                                      Block([SympyAssignment(field.center, field[neighbor])]),
+                                      Block([SympyAssignment(field.center, else_branch)]))]
+                return result
+
+        else:
+            raise NotImplementedError("Contact angle only implemented for phase-fields which have a single "
+                                      "value for each cell")
+
+    def __hash__(self):
+        return hash("ContactAngle")
+
+    def __eq__(self, other):
+        if not isinstance(other, ContactAngle):
+            return False
+        return self.__dict__ == other.__dict__

src/lbmpy/phasefield_allen_cahn/derivatives.py

+from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
+
+import sympy as sp
+
+
+def laplacian_symbolic(field, stencil):
+    r"""
+    Get a symbolic expression for the laplacian of a field.
+    Args:
+        field: the field on which the laplacian is applied to
+        stencil: stencil to derive the finite difference for the laplacian (2nd order isotropic)
+    """
+    lap = sp.simplify(0)
+    for i in range(stencil.D):
+        deriv = FiniteDifferenceStencilDerivation((i, i), stencil)
+        for j in range(stencil.D):
+            # assume the stencil is symmetric
+            deriv.assume_symmetric(dim=j, anti_symmetric=False)
+
+        # set weights for missing degrees of freedom in the calculation and assume the stencil is isotropic
+        if stencil.D == 2 and stencil.Q == 9:
+            res = deriv.get_stencil(isotropic=True)
+            lap += res.apply(field.center)
+        elif stencil.D == 2 and stencil.Q == 25:
+            deriv.set_weight((2, 0), sp.Rational(1, 10))
+
+            res = deriv.get_stencil(isotropic=True)
+            lap += res.apply(field.center)
+        elif stencil.D == 3 and stencil.Q == 15:
+            deriv.set_weight((0, 0, 0), sp.Rational(-32, 27))
+            res = deriv.get_stencil(isotropic=True)
+            lap += res.apply(field.center)
+        elif stencil.D == 3 and stencil.Q == 19:
+            res = deriv.get_stencil(isotropic=True)
+            lap += res.apply(field.center)
+        elif stencil.D == 3 and stencil.Q == 27:
+            deriv.set_weight((0, 0, 0), sp.Rational(-38, 27))
+            res = deriv.get_stencil(isotropic=True)
+            lap += res.apply(field.center)
+        else:
+            raise ValueError(f"stencil with {stencil.D} dimensions and {stencil.Q} entries is not supported")
+
+    return lap
+
+
+def isotropic_gradient_symbolic(field, stencil):
+    r"""
+    Get a symbolic expression for the isotropic gradient of the phase-field
+    Args:
+        field: the field on which the isotropic gradient is applied
+        stencil: stencil to derive the finite difference for the gradient (2nd order isotropic)
+    """
+    deriv = FiniteDifferenceStencilDerivation((0,), stencil)
+
+    deriv.assume_symmetric(0, anti_symmetric=True)
+    deriv.assume_symmetric(1, anti_symmetric=False)
+    if stencil.D == 3:
+        deriv.assume_symmetric(2, anti_symmetric=False)
+
+    # set weights for missing degrees of freedom in the calculation and assume the stencil is isotropic
+    # furthermore the stencils gets rotated to get the y and z components
+    if stencil.D == 2 and stencil.Q == 9:
+        res = deriv.get_stencil(isotropic=True)
+        grad = [res.apply(field.center), res.rotate_weights_and_apply(field.center, (0, 1)), 0]
+    elif stencil.D == 2 and stencil.Q == 25:
+        deriv.set_weight((2, 0), sp.Rational(1, 10))
+
+        res = deriv.get_stencil(isotropic=True)
+        grad = [res.apply(field.center), res.rotate_weights_and_apply(field.center, (0, 1)), 0]
+    elif stencil.D == 3 and stencil.Q == 15:
+        res = deriv.get_stencil(isotropic=True)
+        grad = [res.apply(field.center),
+                res.rotate_weights_and_apply(field.center, (0, 1)),
+                res.rotate_weights_and_apply(field.center, (1, 2))]
+    elif stencil.D == 3 and stencil.Q == 19:
+        deriv.set_weight((0, 0, 0), sp.sympify(0))
+        deriv.set_weight((1, 0, 0), sp.Rational(1, 6))
+
+        res = deriv.get_stencil(isotropic=True)
+        grad = [res.apply(field.center),
+                res.rotate_weights_and_apply(field.center, (0, 1)),
+                res.rotate_weights_and_apply(field.center, (1, 2))]
+    elif stencil.D == 3 and stencil.Q == 27:
+        deriv.set_weight((0, 0, 0), sp.sympify(0))
+        deriv.set_weight((1, 0, 0), sp.Rational(2, 9))
+
+        res = deriv.get_stencil(isotropic=True)
+        grad = [res.apply(field.center),
+                res.rotate_weights_and_apply(field.center, (0, 1)),
+                res.rotate_weights_and_apply(field.center, (1, 2))]
+    else:
+        raise ValueError(f"stencil with {stencil.D} dimensions and {stencil.Q} entries is not supported")
+
+    return grad

src/lbmpy/phasefield_allen_cahn/kernel_equations.py

+from typing import Union
+
+from pystencils import Assignment, AssignmentCollection, Field
+from pystencils.sympyextensions import scalar_product
+
+from lbmpy import pdf_initialization_assignments
+from lbmpy.advanced_streaming.utility import get_accessor
+from lbmpy.creationfunctions import LBMConfig
+from lbmpy.fieldaccess import StreamPushTwoFieldsAccessor
+from lbmpy.methods.abstractlbmethod import LbmCollisionRule
+from lbmpy.utils import second_order_moment_tensor
+
+from lbmpy.phasefield_allen_cahn.derivatives import isotropic_gradient_symbolic, laplacian_symbolic
+from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters, ThermocapillaryParameters
+
+import sympy as sp
+
+
+def chemical_potential_symbolic(phi_field, stencil, beta, kappa):
+    r"""
+    Get a symbolic expression for the chemical potential according to equation (5) in PhysRevE.96.053301.
+    Args:
+        phi_field: the phase-field on which the chemical potential is applied
+        stencil: stencil to derive the finite difference for the laplacian (2nd order isotropic)
+        beta: coefficient related to surface tension and interface thickness
+        kappa: coefficient related to surface tension and interface thickness
+    """
+    lap = laplacian_symbolic(phi_field, stencil)
+    # get the chemical potential
+    four = sp.Rational(4, 1)
+    one = sp.Rational(1, 1)
+    half = sp.Rational(1, 2)
+    mu = four * beta * phi_field.center * (phi_field.center - one) * (phi_field.center - half) - kappa * lap
+    return mu
+
+
+def normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil=None):
+    r"""
+    Get a symbolic expression for the normalized isotropic gradient of the phase-field
+    Args:
+        phi_field: the phase-field on which the normalized isotropic gradient is applied
+        stencil: stencil of the lattice Boltzmann step
+        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+        field. If it is not given the stencil of the LB method will be applied.
+    """
+    if fd_stencil is None:
+        fd_stencil = stencil
+
+    tmp = (sum(map(lambda x: x * x, isotropic_gradient_symbolic(phi_field, fd_stencil))) + sp.Float(1e-32)) ** 0.5
+
+    result = [x / tmp for x in isotropic_gradient_symbolic(phi_field, fd_stencil)]
+    return result
+
+
+def pressure_force(phi_field, lb_method, stencil, density_heavy, density_light, fd_stencil=None):
+    r"""
+    Get a symbolic expression for the pressure force
+    Args:
+        phi_field: phase-field
+        lb_method: lattice boltzmann method used for hydrodynamics
+        stencil: stencil of the lattice Boltzmann step
+        density_heavy: density of the heavier fluid
+        density_light: density of the lighter fluid
+        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+        field. If it is not given the stencil of the LB method will be applied.
+    """
+    if fd_stencil is None:
+        fd_stencil = stencil
+
+    cqc = lb_method.conserved_quantity_computation
+    rho = cqc.density_deviation_symbol
+
+    iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
+    result = list(map(lambda x: sp.Rational(-1, 3) * rho * (density_heavy - density_light) * x, iso_grad))
+    return result
+
+
+def viscous_force(phi_field, lb_method, tau, density_heavy, density_light, fd_stencil=None):
+    r"""
+    Get a symbolic expression for the viscous force
+    Args:
+        phi_field: phase-field
+        lb_method: lattice boltzmann method used for hydrodynamics
+        tau: relaxation time of the hydrodynamic lattice boltzmann step
+        density_heavy: density of the heavier fluid
+        density_light: density of the lighter fluid
+        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+        field. If it is not given the stencil of the LB method will be applied.
+    """
+    stencil = lb_method.stencil
+
+    if fd_stencil is None:
+        fd_stencil = stencil
+
+    iso_grad = sp.Matrix(isotropic_gradient_symbolic(phi_field, fd_stencil)[:stencil.D])
+
+    f_neq = sp.Matrix(lb_method.pre_collision_pdf_symbols) - lb_method.get_equilibrium_terms()
+    stress_tensor = second_order_moment_tensor(f_neq, lb_method.stencil)
+    normal_stress_tensor = stress_tensor * iso_grad
+
+    density_difference = density_heavy - density_light
+
+    # Calculate Viscous Force MRT
+    half = sp.Rational(1, 2)
+    fmx = (half - tau) * normal_stress_tensor[0] * density_difference
+    fmy = (half - tau) * normal_stress_tensor[1] * density_difference
+    fmz = (half - tau) * normal_stress_tensor[2] * density_difference if stencil.D == 3 else 0
+
+    return [fmx, fmy, fmz]
+
+
+def surface_tension_force(phi_field, stencil, parameters, fd_stencil=None):
+    r"""
+    Get a symbolic expression for the surface tension force
+    Args:
+        phi_field: the phase-field on which the chemical potential is applied
+        stencil: stencil of the lattice Boltzmann step
+        parameters: AllenCahnParameters
+        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+        field. If it is not given the stencil of the LB method will be applied.
+    """
+    beta = parameters.beta
+    kappa = parameters.kappa
+
+    if fd_stencil is None:
+        fd_stencil = stencil
+
+    chemical_potential = chemical_potential_symbolic(phi_field, fd_stencil, beta, kappa)
+    iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
+    return [chemical_potential * x for x in iso_grad]
+
+
+def thermocapillary_surface_tension_force(phi_field, temperature_field,
+                                          stencil, parameters, fd_stencil=None):
+    r"""
+    Get a symbolic expression for the surface tension force
+    Args:
+        phi_field: the phase-field on which the chemical potential is applied
+        temperature_field: the temperature field which contains the temperature for each cell
+        stencil: stencil of the lattice Boltzmann step
+        parameters: AllenCahnParameters
+        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+        field. If it is not given the stencil of the LB method will be applied.
+    """
+    if fd_stencil is None:
+        fd_stencil = stencil
+
+    sigma_ref = parameters.symbolic_sigma_ref
+    sigma_t = parameters.symbolic_sigma_t
+    tmp_ref = parameters.symbolic_tmp_ref
+    interface_thickness = parameters.symbolic_interface_thickness
+
+    sigma = sigma_ref + sigma_t * (temperature_field.center - tmp_ref)
+    beta = sp.Rational(12, 1) * (sigma / interface_thickness)
+    kappa = sp.Rational(3, 2) * sigma * interface_thickness
+
+    chemical_potential = chemical_potential_symbolic(phi_field, fd_stencil, beta, kappa)
+    gradient_phi = isotropic_gradient_symbolic(phi_field, fd_stencil)
+    gradient_temp = isotropic_gradient_symbolic(temperature_field, fd_stencil)
+    magnitude_phi = sum([x * x for x in gradient_phi])
+
+    dot_temperature_phase = scalar_product(gradient_temp, gradient_phi)
+    delta_s = sp.Rational(3, 2) * interface_thickness * sigma_t
+
+    return [chemical_potential * gp + delta_s * (magnitude_phi * gt - dot_temperature_phase * gp) for gp, gt in
+            zip(gradient_phi, gradient_temp)]
+
+
+def hydrodynamic_force(phi_field, lb_method,
+                       parameters: Union[AllenCahnParameters, ThermocapillaryParameters],
+                       body_force, fd_stencil=None,
+                       temperature_field=None):
+    r"""
+    Get a symbolic expression for the hydrodynamic force. If a temperature field is provided the hydrodynamic force
+    for thermocapillary simulations is derived.
+    Args:
+        phi_field: phase-field
+        lb_method: Lattice boltzmann method used for hydrodynamics
+        parameters: AllenCahnParameters
+        body_force: force acting on the fluids. Usually the gravity
+        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+        field. If it is not given the stencil of the LB method will be applied.
+        temperature_field: the temperature field which contains the temperature for each cell
+    """
+    stencil = lb_method.stencil
+
+    if fd_stencil is None:
+        fd_stencil = stencil
+
+    density_heavy = parameters.symbolic_density_heavy
+    density_light = parameters.symbolic_density_light
+    tau_L = parameters.symbolic_tau_light
+    tau_H = parameters.symbolic_tau_heavy
+    tau = sp.Rational(1, 2) + tau_L + phi_field.center * (tau_H - tau_L)
+
+    fp = pressure_force(phi_field, lb_method, stencil, density_heavy, density_light, fd_stencil)
+    fm = viscous_force(phi_field, lb_method, tau, density_heavy, density_light, fd_stencil)
+
+    if temperature_field is None:
+        fs = surface_tension_force(phi_field, stencil, parameters, fd_stencil)
+    else:
+        assertion_string = "For thermocapillary ThermocapillaryParameters needs to be passed"
+        assert isinstance(parameters, ThermocapillaryParameters), assertion_string
+        fs = thermocapillary_surface_tension_force(phi_field, temperature_field, stencil, parameters, fd_stencil)
+
+    result = []
+    for i in range(stencil.D):
+        result.append(fs[i] + fp[i] + fm[i] + body_force[i])
+
+    return result
+
+
+def interface_tracking_force(phi_field, stencil, parameters: AllenCahnParameters, fd_stencil=None,
+                             phi_heavy=1, phi_light=0):
+    r"""
+    Get a symbolic expression for the hydrodynamic force
+    Args:
+        phi_field: phase-field
+        stencil: stencil of the phase-field distribution lattice Boltzmann step
+        parameters: AllenCahnParameters
+        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+        field. If it is not given the stencil of the LB method will be applied.
+        phi_heavy: phase field value in the bulk of the heavy fluid
+        phi_light: phase field value in the bulk of the light fluid
+
+    """
+    if fd_stencil is None:
+        fd_stencil = stencil
+
+    phi_zero = sp.Rational(1, 2) * (phi_light + phi_heavy)
+
+    normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil)
+    result = []
+    interface_thickness = parameters.symbolic_interface_thickness
+    for i in range(stencil.D):
+        fraction = (sp.Rational(1, 1) - sp.Rational(4, 1) * (phi_field.center - phi_zero) ** 2) / interface_thickness
+        result.append(sp.Rational(1, 3) * fraction * normal_fd[i])
+
+    return result
+
+
+def hydrodynamic_force_assignments(velocity_field, phi_field, lb_method,
+                                   parameters: AllenCahnParameters,
+                                   body_force, fd_stencil=None, sub_iterations=2,
+                                   temperature_field=None):
+
+    r"""
+    Get a symbolic expression for the hydrodynamic force
+    Args:
+        velocity_field: velocity
+        phi_field: phase-field
+        lb_method: Lattice boltzmann method used for hydrodynamics
+        parameters: AllenCahnParameters
+        body_force: force acting on the fluids. Usually the gravity
+        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+        field. If it is not given the stencil of the LB method will be applied.
+        sub_iterations: number of sub iterations for the hydrodynamic force
+        temperature_field: the temperature field which contains the temperature for each cell
+    """
+
+    rho_L = parameters.symbolic_density_light
+    rho_H = parameters.symbolic_density_heavy
+    density = rho_L + phi_field.center * (rho_H - rho_L)
+
+    stencil = lb_method.stencil
+    # method has to have a force model
+    symbolic_force = lb_method.force_model.symbolic_force_vector
+
+    force = hydrodynamic_force(phi_field, lb_method, parameters, body_force, fd_stencil=fd_stencil,
+                               temperature_field=temperature_field)
+
+    cqc = lb_method.conserved_quantity_computation
+
+    u_symp = cqc.velocity_symbols
+    cqe = cqc.equilibrium_input_equations_from_pdfs(lb_method.pre_collision_pdf_symbols)
+    cqe = cqe.new_without_subexpressions()
+
+    cqe_velocity = [eq.rhs for eq in cqe.main_assignments[1:]]
+    index = 0
+    aleph = sp.symbols(f"aleph_:{stencil.D * sub_iterations}")
+
+    force_Assignments = []
+
+    for i in range(stencil.D):
+        force_Assignments.append(Assignment(aleph[i], velocity_field.center_vector[i]))
+        index += 1
+
+    for k in range(sub_iterations - 1):
+        subs_dict = dict(zip(u_symp, aleph[k * stencil.D:index]))
+        for i in range(stencil.D):
+            new_force = force[i].subs(subs_dict) / density
+            force_Assignments.append(Assignment(aleph[index], cqe_velocity[i].subs({symbolic_force[i]: new_force})))
+            index += 1
+
+    subs_dict = dict(zip(u_symp, aleph[index - stencil.D:index]))
+
+    for i in range(stencil.D):
+        force_Assignments.append(Assignment(symbolic_force[i], force[i].subs(subs_dict)))
+
+    return force_Assignments
+
+
+def add_interface_tracking_force(update_rule: LbmCollisionRule, force):
+    r"""
+     Adds the interface tracking force to a lattice Boltzmann update rule
+     Args:
+         update_rule: lattice Boltzmann update rule
+         force: interface tracking force
+     """
+    method = update_rule.method
+    symbolic_force = method.force_model.symbolic_force_vector
+
+    for i in range(method.stencil.D):
+        update_rule.subexpressions += [Assignment(symbolic_force[i], force[i])]
+
+    update_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
+
+    return update_rule
+
+
+def add_hydrodynamic_force(update_rule: LbmCollisionRule, force, phi_field,
+                           hydro_pdfs, parameters: AllenCahnParameters, lbm_config: LBMConfig = None):
+    r"""
+     Adds the interface tracking force to a lattice Boltzmann update rule
+     Args:
+         update_rule: lattice Boltzmann update rule
+         force: interface tracking force
+         phi_field: phase-field
+         hydro_pdfs: source field of the hydrodynamic PDFs
+         parameters: AllenCahnParameters
+         lbm_config: LBMConfig to determine the streaming scheme
+     """
+    if lbm_config is None:
+        accessor = StreamPushTwoFieldsAccessor()
+    else:
+        accessor = get_accessor(lbm_config.streaming_pattern, lbm_config.timestep)
+    method = update_rule.method
+    reads = accessor.read(hydro_pdfs, method.stencil)
+
+    # First apply force according to Allen Cahn model
+    rho_L = parameters.symbolic_density_light
+    rho_H = parameters.symbolic_density_heavy
+    density = rho_L + phi_field.center * (rho_H - rho_L)
+
+    force_subs = {f: f / density for f in method.force_model.symbolic_force_vector}
+    update_rule = update_rule.new_with_substitutions(force_subs)
+    update_rule.subexpressions += force
+
+    # Then add missing conversed quantities that occur in the force terms
+    cqc = method.conserved_quantity_computation
+    density = cqc.density_symbol
+    density_deviation = cqc.density_deviation_symbol
+    free_symbols = update_rule.free_symbols
+
+    if density_deviation in free_symbols:
+        t = cqc.output_equations_from_pdfs(reads,
+                                           {"density_deviation": density_deviation}).new_without_subexpressions()
+        update_rule.add_subexpression(t.main_assignments[0].rhs, t.main_assignments[0].lhs)
+
+    if density in free_symbols:
+        t = cqc.output_equations_from_pdfs(reads, {"density": density}).new_without_subexpressions()
+        update_rule.add_subexpression(t.main_assignments[0].rhs, t.main_assignments[0].lhs)
+
+    update_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
+
+    subs_dict = {f: read for f, read in zip(method.pre_collision_pdf_symbols, reads)}
+
+    up = update_rule.new_with_substitutions(subs_dict)
+    result = LbmCollisionRule(lb_method=method, main_assignments=up.main_assignments,
+                              subexpressions=up.subexpressions, simplification_hints=up.simplification_hints,
+                              subexpression_symbol_generator=up.subexpression_symbol_generator)
+
+    return result
+
+
+def initializer_kernel_phase_field_lb(lb_method, phi, velocity, ac_pdfs, parameters: AllenCahnParameters,
+                                      fd_stencil=None, **kwargs):
+    r"""
+    Returns an assignment list for initializing the phase-field distribution functions
+    Args:
+        lb_method: lattice Boltzmann method of the phase-field lattice Boltzmann step
+        phi: order parameter of the Allen-Cahn LB step (phase field)
+        velocity: initial velocity
+        ac_pdfs: source field of the Allen-Cahn PDFs
+        parameters: AllenCahnParameters
+        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+        field. If it is not given the stencil of the LB method will be applied.
+        kwargs: keyword arguments for pdf_initialization_assignments
+    """
+
+    h_updates = pdf_initialization_assignments(lb_method, phi, velocity, ac_pdfs, **kwargs)
+    force_h = interface_tracking_force(phi, lb_method.stencil, parameters,
+                                       fd_stencil=fd_stencil)
+
+    cqc = lb_method.conserved_quantity_computation
+
+    rho = cqc.density_symbol
+    u_symp = cqc.velocity_symbols
+    symbolic_force = lb_method.force_model.symbolic_force_vector
+
+    macro_quantities = []
+
+    if isinstance(velocity, Field):
+        velocity = velocity.center_vector
+
+    if isinstance(phi, Field):
+        phi = phi.center
+
+    for i in range(lb_method.stencil.D):
+        macro_quantities.append(Assignment(symbolic_force[i], force_h[i]))
+
+    for i in range(lb_method.stencil.D):
+        macro_quantities.append(Assignment(u_symp[i],
+                                           velocity[i] - sp.Rational(1, 2) * symbolic_force[i]))
+
+    h_updates = AssignmentCollection(main_assignments=h_updates.main_assignments, subexpressions=macro_quantities)
+    h_updates = h_updates.new_with_substitutions({rho: phi})
+
+    return h_updates
+
+
+def initializer_kernel_hydro_lb(lb_method, pressure, velocity, hydro_pdfs, **kwargs):
+    r"""
+    Returns an assignment list for initializing the velocity distribution functions
+    Args:
+        lb_method: lattice Boltzmann method of the hydrodynamic lattice Boltzmann step
+        pressure: order parameter of the hydrodynamic LB step (pressure)
+        velocity: initial velocity
+        hydro_pdfs: source field of the hydrodynamic PDFs
+        kwargs: keyword arguments for pdf_initialization_assignments
+    """
+    symbolic_force = lb_method.force_model.symbolic_force_vector
+    force_subs = {f: 0 for f in symbolic_force}
+
+    g_updates = pdf_initialization_assignments(lb_method, pressure, velocity, hydro_pdfs, **kwargs)
+    g_updates = g_updates.new_with_substitutions(force_subs)
+
+    return g_updates

src/lbmpy/phasefield_allen_cahn/numerical_solver.py

+from pystencils import Assignment, AssignmentCollection
+from pystencils.sympyextensions import scalar_product
+from pystencils.simp.subexpression_insertion import insert_constants
+
+from lbmpy.phasefield_allen_cahn.derivatives import isotropic_gradient_symbolic, laplacian_symbolic
+
+import sympy as sp
+
+VELOCITY_SYMBOLS = sp.symbols(f"u_:{3}")
+GRAD_T_SYMBOLS = sp.symbols(f"gratT_:{3}")
+GRAD_K_SYMBOLS = sp.symbols(f"gratK_:{3}")
+LAPLACIAN_SYMBOL = sp.Symbol("lap")
+
+
+def get_runge_kutta_update_assignments(stencil, phase_field, temperature_field, velocity_field, runge_kutta_fields,
+                                       conduction_h, conduction_l, heat_capacity_h, heat_capacity_l,
+                                       density, stabiliser=1):
+    dimensions = len(stencil[0])
+
+    grad_temperature = isotropic_gradient_symbolic(temperature_field, stencil)
+    lap_temperature = laplacian_symbolic(temperature_field, stencil)
+    grad_conduction = _get_conduction_gradient(stencil, phase_field, conduction_h, conduction_l)
+
+    grad_rk = [isotropic_gradient_symbolic(rk, stencil) for rk in runge_kutta_fields]
+    lap_rk = [laplacian_symbolic(rk, stencil) for rk in runge_kutta_fields]
+
+    dot_u_grad_t = scalar_product(VELOCITY_SYMBOLS[:dimensions], GRAD_T_SYMBOLS[:dimensions])
+    dot_grad_k_grad_t = scalar_product(GRAD_K_SYMBOLS[:dimensions], GRAD_T_SYMBOLS[:dimensions])
+
+    conduction = conduction_l + phase_field.center * sp.nsimplify(conduction_h - conduction_l)
+    conduction_times_lap = conduction * LAPLACIAN_SYMBOL
+
+    heat_capacity = heat_capacity_l + phase_field.center * sp.nsimplify(heat_capacity_h - heat_capacity_l)
+
+    rho_cp = 1.0 / (density * heat_capacity)
+    end_term = dot_grad_k_grad_t + conduction_times_lap
+
+    update_stage_1 = temperature_field.center + stabiliser * 0.5 * (-1.0 * dot_u_grad_t + rho_cp * end_term)
+    subexpressions_1 = _get_stage(dimensions, velocity_field, grad_temperature, grad_conduction, lap_temperature)
+    stage_1 = AssignmentCollection(main_assignments=[Assignment(runge_kutta_fields[0].center, update_stage_1)],
+                                   subexpressions=subexpressions_1)
+
+    if len(runge_kutta_fields) == 1:
+        update_stage_2 = temperature_field.center + stabiliser * (-1.0 * dot_u_grad_t + rho_cp * end_term)
+        subexpressions_2 = _get_stage(dimensions, velocity_field, grad_rk[0], grad_conduction, lap_rk[0])
+        stage_2 = AssignmentCollection(main_assignments=[Assignment(temperature_field.center, update_stage_2)],
+                                       subexpressions=subexpressions_2)
+
+        return [insert_constants(ac) for ac in [stage_1, stage_2]]
+
+    update_stage_2 = temperature_field.center + stabiliser * 0.5 * (-1.0 * dot_u_grad_t + rho_cp * end_term)
+    subexpressions_2 = _get_stage(dimensions, velocity_field, grad_rk[0], grad_conduction, lap_rk[0])
+    stage_2 = AssignmentCollection(main_assignments=[Assignment(runge_kutta_fields[1].center, update_stage_2)],
+                                   subexpressions=subexpressions_2)
+
+    update_stage_3 = temperature_field.center + stabiliser * 1.0 * (-1.0 * dot_u_grad_t + rho_cp * end_term)
+    subexpressions_3 = _get_stage(dimensions, velocity_field, grad_rk[1], grad_conduction, lap_rk[1])
+    stage_3 = AssignmentCollection(main_assignments=[Assignment(runge_kutta_fields[2].center, update_stage_3)],
+                                   subexpressions=subexpressions_3)
+
+    update_stage_4 = stabiliser * 1.0 * (-1.0 * dot_u_grad_t + rho_cp * end_term)
+    rk_update = 2.0 * runge_kutta_fields[0].center + 4.0 * runge_kutta_fields[1].center + 2.0 * runge_kutta_fields[
+        2].center
+    update_stage_4 = (1.0 - 4.0 / 3.0) * temperature_field.center + (rk_update - update_stage_4) / 6.0
+    subexpressions_4 = _get_stage(dimensions, velocity_field, grad_rk[2], grad_conduction, lap_rk[2])
+    stage_4 = AssignmentCollection(main_assignments=[Assignment(temperature_field.center, update_stage_4)],
+                                   subexpressions=subexpressions_4)
+
+    return [insert_constants(ac) for ac in [stage_1, stage_2, stage_3, stage_4]]
+
+
+def get_initialiser_assignments(temperature_field, runge_kutta_fields):
+    result = []
+    for i in range(len(runge_kutta_fields)):
+        result.append(Assignment(runge_kutta_fields[i].center, temperature_field.center))
+
+    return result
+
+
+def _get_conduction_gradient(stencil, phase_field, conduction_h, conduction_l):
+    dimensions = len(stencil[0])
+    grad_phase = isotropic_gradient_symbolic(phase_field, stencil)
+
+    free_symbols = set()
+    for i in range(dimensions):
+        free_symbols.update(grad_phase[i].free_symbols)
+
+    subs_dict = {}
+
+    for f in free_symbols:
+        subs_dict[f] = interpolate_field_access(f, conduction_h, conduction_l)
+
+    result = list()
+    for i in range(dimensions):
+        eq = grad_phase[i].subs(subs_dict)
+        # replace very small numbers by zero
+        eq = eq.xreplace(dict([(n, 0) for n in eq.atoms(sp.Float) if abs(n) < 1e-16]))
+        result.append(eq)
+
+    return result
+
+
+def interpolate_field_access(field_access, upper, lower):
+    return lower + field_access * sp.nsimplify(upper - lower)
+
+
+def _get_stage(dimensions, velocity_field, gradient_t, gradient_k, laplacian):
+    result = list()
+
+    for i in range(dimensions):
+        result.append(Assignment(VELOCITY_SYMBOLS[i], velocity_field.center_vector[i]))
+
+    for i in range(dimensions):
+        result.append(Assignment(GRAD_T_SYMBOLS[i], gradient_t[i]))
+
+    for i in range(dimensions):
+        result.append(Assignment(GRAD_K_SYMBOLS[i], gradient_k[i]))
+
+    result.append(Assignment(LAPLACIAN_SYMBOL, laplacian))
+
+    return result

src/lbmpy/phasefield_allen_cahn/parameter_calculation.py

+import math
+import sympy as sp
+
+
+class AllenCahnParameters:
+    def __init__(self, density_heavy: float, density_light: float,
+                 dynamic_viscosity_heavy: float, dynamic_viscosity_light: float,
+                 surface_tension: float, mobility: float = 0.2,
+                 gravitational_acceleration: float = 0.0, interface_thickness: int = 5):
+
+        self.density_heavy = density_heavy
+        self.density_light = density_light
+        self.dynamic_viscosity_heavy = dynamic_viscosity_heavy
+        self.dynamic_viscosity_light = dynamic_viscosity_light
+        self.surface_tension = surface_tension
+        self.mobility = mobility
+        self.gravitational_acceleration = gravitational_acceleration
+        self.interface_thickness = interface_thickness
+
+    @property
+    def kinematic_viscosity_heavy(self):
+        return self.dynamic_viscosity_heavy / self.density_heavy
+
+    @property
+    def kinematic_viscosity_light(self):
+        return self.dynamic_viscosity_light / self.density_light
+
+    @property
+    def relaxation_time_heavy(self):
+        return 3.0 * self.kinematic_viscosity_heavy
+
+    @property
+    def relaxation_time_light(self):
+        return 3.0 * self.kinematic_viscosity_light
+
+    @property
+    def omega_phi(self):
+        return 1.0 / (0.5 + (3.0 * self.mobility))
+
+    @property
+    def symbolic_density_heavy(self):
+        return sp.Symbol("rho_H")
+
+    @property
+    def symbolic_density_light(self):
+        return sp.Symbol("rho_L")
+
+    @property
+    def symbolic_tau_heavy(self):
+        return sp.Symbol("tau_H")
+
+    @property
+    def symbolic_tau_light(self):
+        return sp.Symbol("tau_L")
+
+    @property
+    def symbolic_omega_phi(self):
+        return sp.Symbol("omega_phi")
+
+    @property
+    def symbolic_surface_tension(self):
+        return sp.Symbol("sigma")
+
+    @property
+    def symbolic_mobility(self):
+        return sp.Symbol("M_m")
+
+    @property
+    def symbolic_gravitational_acceleration(self):
+        return sp.Symbol("F_g")
+
+    @property
+    def symbolic_interface_thickness(self):
+        return sp.Symbol("W")
+
+    @property
+    def beta(self):
+        return sp.Rational(12, 1) * (self.symbolic_surface_tension / self.symbolic_interface_thickness)
+
+    @property
+    def kappa(self):
+        return sp.Rational(3, 2) * self.symbolic_surface_tension * self.symbolic_interface_thickness
+
+    def omega(self, phase_field):
+        tau_L = self.symbolic_tau_light
+        tau_H = self.symbolic_tau_heavy
+        tau = sp.Rational(1, 2) + tau_L + phase_field.center * (tau_H - tau_L)
+        return sp.simplify(1 / tau)
+
+    def parameter_map(self):
+        result = {self.symbolic_density_heavy: self.density_heavy,
+                  self.symbolic_density_light: self.density_light,
+                  self.symbolic_tau_heavy: self.relaxation_time_heavy,
+                  self.symbolic_tau_light: self.relaxation_time_light,
+                  self.symbolic_omega_phi: self.omega_phi,
+                  self.symbolic_gravitational_acceleration: self.gravitational_acceleration,
+                  self.symbolic_interface_thickness: self.interface_thickness,
+                  self.symbolic_mobility: self.mobility,
+                  self.symbolic_surface_tension: self.surface_tension}
+        return result
+
+    @property
+    def symbolic_to_numeric_map(self):
+        return {t.name: self.parameter_map()[t] for t in self.parameter_map()}
+
+    @staticmethod
+    def _parameter_strings():
+        names = ("Density heavy phase",
+                 "Density light phase",
+                 "Relaxation time heavy phase",
+                 "Relaxation time light phase",
+                 "Relaxation rate Allen Cahn LB",
+                 "Gravitational acceleration",
+                 "Interface thickness",
+                 "Mobility",
+                 "Surface tension")
+        return names
+
+    def _repr_html_(self):
+        names = self._parameter_strings()
+        table = """
+        <table style="border:none; width: 100%">
+            <tr {nb}>
+                <th {nb} >Name</th>
+                <th {nb} >SymPy Symbol </th>
+                <th {nb} >Value</th>
+            </tr>
+            {content}
+        </table>
+        """
+        content = ""
+        for name, (symbol, value) in zip(names, self.parameter_map().items()):
+            vals = {
+                'Name': name,
+                'Sympy Symbol': sp.latex(symbol),
+                'numeric value': sp.latex(value),
+                'nb': 'style="border:none"',
+            }
+            content += """<tr {nb}>
+                            <td {nb}>{Name}</td>
+                            <td {nb}>${Sympy Symbol}$</td>
+                            <td {nb}>${numeric value}$</td>
+                         </tr>\n""".format(**vals)
+        return table.format(content=content, nb='style="border:none"')
+
+
+class ThermocapillaryParameters(AllenCahnParameters):
+    def __init__(self, density_heavy: float, density_light: float,
+                 dynamic_viscosity_heavy: float, dynamic_viscosity_light: float,
+                 surface_tension: float,
+                 heat_conductivity_heavy: float, heat_conductivity_light: float,
+                 mobility: float = 0.2,
+                 gravitational_acceleration: float = 0.0, interface_thickness: int = 5,
+                 sigma_ref: float = 5e-3, sigma_t: float = 2e-4, tmp_ref: float = 0,
+                 velocity_wall: float = 0.0, reference_time: int = 10):
+
+        self.heat_conductivity_heavy = heat_conductivity_heavy
+        self.heat_conductivity_light = heat_conductivity_light
+        self.sigma_ref = sigma_ref
+        self.sigma_t = sigma_t
+        self.tmp_ref = tmp_ref
+        self.velocity_wall = velocity_wall
+        self.reference_time = reference_time
+
+        super(ThermocapillaryParameters, self).__init__(density_heavy, density_light,
+                                                        dynamic_viscosity_heavy, dynamic_viscosity_light,
+                                                        surface_tension, mobility,
+                                                        gravitational_acceleration, interface_thickness)
+
+    @property
+    def symbolic_heat_conductivity_heavy(self):
+        return sp.Symbol("kappa_H")
+
+    @property
+    def symbolic_heat_conductivity_light(self):
+        return sp.Symbol("kappa_L")
+
+    @property
+    def symbolic_sigma_ref(self):
+        return sp.Symbol("sigma_ref")
+
+    @property
+    def symbolic_sigma_t(self):
+        return sp.Symbol("sigma_T")
+
+    @property
+    def symbolic_tmp_ref(self):
+        return sp.Symbol("T_ref")
+
+    def parameter_map(self):
+        result = {self.symbolic_density_heavy: self.density_heavy,
+                  self.symbolic_density_light: self.density_light,
+                  self.symbolic_tau_heavy: self.relaxation_time_heavy,
+                  self.symbolic_tau_light: self.relaxation_time_light,
+                  self.symbolic_omega_phi: self.omega_phi,
+                  self.symbolic_gravitational_acceleration: self.gravitational_acceleration,
+                  self.symbolic_interface_thickness: self.interface_thickness,
+                  self.symbolic_mobility: self.mobility,
+                  self.symbolic_surface_tension: self.surface_tension,
+                  self.symbolic_heat_conductivity_heavy: self.heat_conductivity_heavy,
+                  self.symbolic_heat_conductivity_light: self.heat_conductivity_light,
+                  self.symbolic_sigma_ref: self.sigma_ref,
+                  self.symbolic_sigma_t: self.sigma_t,
+                  self.symbolic_tmp_ref: self.tmp_ref}
+        return result
+
+    @staticmethod
+    def _parameter_strings():
+        names = ("Density heavy phase",
+                 "Density light phase",
+                 "Relaxation time heavy phase",
+                 "Relaxation time light phase",
+                 "Relaxation rate Allen Cahn LB",
+                 "Gravitational acceleration",
+                 "Interface thickness",
+                 "Mobility",
+                 "Surface tension",
+                 "Heat Conductivity Heavy",
+                 "Heat Conductivity Light",
+                 "Sigma Ref",
+                 "Sigma T",
+                 "Temperature Ref")
+        return names
+
+
+def calculate_parameters_rti(reference_length=256,
+                             reference_time=16000,
+                             density_heavy=1.0,
+                             capillary_number=0.26,
+                             reynolds_number=3000,
+                             atwood_number=0.5,
+                             peclet_number=1000,
+                             density_ratio=3,
+                             viscosity_ratio=1):
+    r"""
+    Calculate the simulation parameters for the Rayleigh-Taylor instability.
+    The calculation is done according to the description in part B of PhysRevE.96.053301.
+
+    Args:
+        reference_length: reference length of the RTI
+        reference_time: chosen reference time
+        density_heavy: density of the heavier fluid
+        capillary_number: capillary number of the simulation represents the relative effect of viscous drag
+                          forces versus surface tension forces
+        reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
+                         and the inertial forces
+        atwood_number: atwood number of the simulation is a dimensionless density ratio
+        peclet_number: peclet number of the simulation is the ratio of the rate of advection
+                       of a physical quantity by the flow to the rate of diffusion of the same quantity
+                       driven by an appropriate gradient
+        density_ratio: density ratio of the heavier and the lighter fluid
+        viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
+    """
+
+    # calculate the gravitational acceleration and the reference velocity
+    gravity_lattice_units = reference_length / (reference_time ** 2 * atwood_number)
+    reference_velocity = math.sqrt(gravity_lattice_units * reference_length)
+
+    dynamic_viscosity_heavy = (density_heavy * reference_velocity * reference_length) / reynolds_number
+    dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
+
+    density_light = density_heavy / density_ratio
+
+    surface_tension = (dynamic_viscosity_heavy * reference_velocity) / capillary_number
+    mobility = (reference_velocity * reference_length) / peclet_number
+
+    parameters = AllenCahnParameters(density_heavy=density_heavy,
+                                     density_light=density_light,
+                                     dynamic_viscosity_heavy=dynamic_viscosity_heavy,
+                                     dynamic_viscosity_light=dynamic_viscosity_light,
+                                     surface_tension=surface_tension,
+                                     mobility=mobility,
+                                     gravitational_acceleration=-gravity_lattice_units)
+    return parameters
+
+
+def calculate_dimensionless_rising_bubble(reference_time=18000,
+                                          density_heavy=1.0,
+                                          bubble_radius=16,
+                                          bond_number=1,
+                                          reynolds_number=40,
+                                          density_ratio=1000,
+                                          viscosity_ratio=100):
+    r"""
+    Calculate the simulation parameters for a rising bubble. The parameter calculation follows the ideas of Weber and
+    is based on the Reynolds number. This means the rising velocity of the bubble is implicitly stated with the
+    Reynolds number
+
+    Args:
+        reference_time: chosen reference time
+        density_heavy: density of the heavier fluid
+        bubble_radius: initial radius of the rising bubble
+        bond_number: also called eötvös number is measuring the importance of gravitational forces compared to
+                     surface tension forces
+        reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
+                 and the inertial forces
+        density_ratio: density ratio of the heavier and the lighter fluid
+        viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
+    """
+
+    bubble_d = bubble_radius * 2
+    gravity_lattice_units = bubble_d / (reference_time ** 2)
+
+    density_light = density_heavy / density_ratio
+
+    dynamic_viscosity_heavy = (density_heavy * math.sqrt(gravity_lattice_units * bubble_d ** 3)) / reynolds_number
+    dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
+
+    surface_tension = (density_heavy - density_light) * gravity_lattice_units * bubble_d ** 2 / bond_number
+    # calculation of the Morton number
+    # Mo = gravitational_acceleration * dynamic_viscosity_heavy / (density_heavy * surface_tension ** 3)
+
+    parameters = AllenCahnParameters(density_heavy=density_heavy,
+                                     density_light=density_light,
+                                     dynamic_viscosity_heavy=dynamic_viscosity_heavy,
+                                     dynamic_viscosity_light=dynamic_viscosity_light,
+                                     surface_tension=surface_tension,
+                                     gravitational_acceleration=-gravity_lattice_units)
+    return parameters
+
+
+def calculate_parameters_taylor_bubble(reference_length=128,
+                                       reference_time=16000,
+                                       density_heavy=1.0,
+                                       diameter_outer=0.0254,
+                                       diameter_inner=0.0127):
+    r"""
+    Calculate the simulation parameters for a rising Taylor bubble in an annulus pipe. The calculation can be found in
+    10.1016/S0009-2509(97)00210-8 by G. Das
+
+    Args:
+        reference_length: chosen reference length in lattice cells
+        reference_time: chosen reference time in latte timesteps
+        density_heavy: density of water in lattice units
+        diameter_outer: diameter of the outer tube
+        diameter_inner: diameter of the inner cylinder
+    """
+
+    # physical parameters #
+    water_rho = 998  # kg/m3
+    air_rho = 1.2047  # kg/m3
+    surface_tension = 0.07286  # kg/s2
+    water_mu = 1.002e-3  # kg/ms
+    air_mu = 1.8205e-5  # kg/ms
+    gravity = 9.81  # m/s2
+
+    # water_nu = water_mu / water_rho  # m2/s
+    # air_nu = air_mu / air_rho  # m2/s
+
+    diameter_fluid = diameter_outer - diameter_inner
+    diameter_ratio = diameter_outer / diameter_inner
+
+    inverse_viscosity_number = math.sqrt((water_rho - air_rho) * water_rho * gravity * diameter_fluid ** 3) / water_mu
+    bond_number = (water_rho - air_rho) * gravity * diameter_fluid ** 2 / surface_tension
+    # morton_number = gravity * water_mu ** 4 * (water_rho - air_rho) / (water_rho ** 2 * surface_tension ** 3)
+
+    diameter_lattice_untis = reference_length / diameter_ratio
+
+    density_light = 1.0 / (water_rho / air_rho)
+    diameter_fluid = reference_length - diameter_lattice_untis
+    gravity_lattice_units = diameter_fluid / reference_time ** 2
+
+    density_diff = density_heavy - density_light
+
+    grav_df_cube = gravity_lattice_units * diameter_fluid ** 3
+    water_mu_lattice_units = math.sqrt(density_diff * density_heavy * grav_df_cube) / inverse_viscosity_number
+    air_mu_lattice_units = water_mu_lattice_units / (water_mu / air_mu)
+
+    dynamic_viscosity_heavy = water_mu_lattice_units
+    dynamic_viscosity_light = air_mu_lattice_units
+
+    surface_tension_lattice_units = density_diff * gravity_lattice_units * diameter_fluid ** 2 / bond_number
+
+    parameters = AllenCahnParameters(density_heavy=density_heavy,
+                                     density_light=density_light,
+                                     dynamic_viscosity_heavy=dynamic_viscosity_heavy,
+                                     dynamic_viscosity_light=dynamic_viscosity_light,
+                                     surface_tension=surface_tension_lattice_units,
+                                     gravitational_acceleration=-gravity_lattice_units)
+    return parameters
+
+
+def calculate_parameters_droplet_migration(radius=32,
+                                           reynolds_number=0.16,
+                                           capillary_number=0.01,
+                                           marangoni_number=0.08,
+                                           peclet_number=1,
+                                           viscosity_ratio=1,
+                                           heat_ratio=1,
+                                           interface_width=4,
+                                           reference_surface_tension=5e-3,
+                                           height=None):
+    r"""
+    Calculate the simulation parameters moving droplet with a laser heat source. The calculation can be found in
+    https://doi.org/10.1016/j.jcp.2013.08.054 by Liu et al.
+
+    Args:
+        radius: radius of the droplet which functions as the reference length
+        reynolds_number: Reynolds number of the simulation
+        capillary_number: Capillary number of the simulation
+        marangoni_number: Marangoni number of the simulation
+        peclet_number: Peclet number of the simulation
+        viscosity_ratio: viscosity ratio between the two fluids
+        heat_ratio: ratio of the heat conductivity in the two fluids
+        interface_width: Resolution of cells for the interface
+        reference_surface_tension: reference surface tension
+        height: height of the simulation domain. If not defined it is asumed to be 2 * radius of the droplet.
+
+    """
+
+    if height is None:
+        height = 2 * radius
+
+    u_char = math.sqrt(reynolds_number * capillary_number * reference_surface_tension / radius)
+    gamma = u_char / radius
+    u_wall = gamma * height
+
+    kinematic_viscosity_heavy = radius * u_char / reynolds_number
+    kinematic_viscosity_light = kinematic_viscosity_heavy / viscosity_ratio
+
+    heat_conductivity_heavy = radius * u_char / marangoni_number
+    heat_conductivity_light = heat_conductivity_heavy / heat_ratio
+
+    mobility_of_interface = gamma * radius * interface_width / peclet_number
+
+    parameters = ThermocapillaryParameters(density_heavy=1.0,
+                                           density_light=1.0,
+                                           dynamic_viscosity_heavy=kinematic_viscosity_heavy,
+                                           dynamic_viscosity_light=kinematic_viscosity_light,
+                                           surface_tension=0.0,
+                                           heat_conductivity_heavy=heat_conductivity_heavy,
+                                           heat_conductivity_light=heat_conductivity_light,
+                                           mobility=mobility_of_interface,
+                                           velocity_wall=u_wall, reference_time=int(1.0 / gamma))
+    return parameters

lbmpy/plot.py → src/lbmpy/plot.py

@@ -139,14 +139,17 @@ class LbGrid:
                 self.annotations[(cell, arrow_position)] = Text(arrow_end[0], arrow_end[1], annotation,
                                                                 **annotation_kwargs)
 
-    def fill_with_default_arrows(self, **kwargs):
+    def fill_with_default_arrows(self, inverted=False, **kwargs):
         """Fills the complete grid with the default pdf arrows"""
         for x in range(self._xCells):
             for y in range(self._yCells):
                 for dx in [-1, 0, 1]:
                     for dy in [-1, 0, 1]:
                         kwargs.setdefault('color', '#bbbbbb')
-                        self.add_arrow((x, y), (dx, dy), (dx, dy), **kwargs)
+                        if not inverted:
+                            self.add_arrow((x, y), (dx, dy), (dx, dy), **kwargs)
+                        else:
+                            self.add_arrow((x, y), (dx, dy), (-dx, -dy), **kwargs)
 
     def draw(self, ax):
         """Draw the grid into a given matplotlib axes object"""

lbmpy/quadratic_equilibrium_construction.py → src/lbmpy/quadratic_equilibrium_construction.py

@@ -46,9 +46,9 @@ def match_generic_equilibrium_ansatz(stencil, equilibrium, u=sp.symbols("u_:3"))
         generic quadratic form, a ValueError is raised
 
     Example:
-          >>> from lbmpy.stencils import get_stencil
+          >>> from lbmpy import LBStencil, Stencil
           >>> from lbmpy.maxwellian_equilibrium import discrete_maxwellian_equilibrium
-          >>> stencil = get_stencil("D2Q9")
+          >>> stencil = LBStencil(Stencil.D2Q9)
           >>> eq = discrete_maxwellian_equilibrium(stencil)
           >>> result = match_generic_equilibrium_ansatz(stencil, eq)
           >>> result[sp.Symbol('A_0')]
@@ -94,12 +94,13 @@ def moment_constraint_equations(stencil, equilibrium, moment_to_value_dict, u=sp
 
 def hydrodynamic_moment_values(up_to_order=3, dim=3, compressible=True):
     """Returns the values of moments that are required to approximate Navier Stokes (if up_to_order=3)"""
-    from lbmpy.maxwellian_equilibrium import get_moments_of_continuous_maxwellian_equilibrium
+    from lbmpy.maxwellian_equilibrium import get_equilibrium_values_of_maxwell_boltzmann_function
     from lbmpy.moments import moments_up_to_order
 
     moms = moments_up_to_order(up_to_order, dim)
     c_s_sq = sp.Symbol("p") / sp.Symbol("rho")
-    moment_values = get_moments_of_continuous_maxwellian_equilibrium(moms, dim=dim, c_s_sq=c_s_sq, order=2)
+    moment_values = get_equilibrium_values_of_maxwell_boltzmann_function(moms, dim=dim, c_s_sq=c_s_sq, order=2,
+                                                                         space="moment")
     if not compressible:
         moment_values = [compressible_to_incompressible_moment_value(m, sp.Symbol("rho"), sp.symbols("u_:3")[:dim])
                          for m in moment_values]

lbmpy/relaxationrates.py → src/lbmpy/relaxationrates.py

 import sympy as sp
 
-from lbmpy.moments import get_order, is_shear_moment
+from lbmpy.moments import is_bulk_moment, is_shear_moment
 
 
 def relaxation_rate_from_lattice_viscosity(nu):
@@ -31,7 +31,7 @@ def get_shear_relaxation_rate(method):
 
     relaxation_rates = set()
     for moment, relax_info in method.relaxation_info_dict.items():
-        if is_shear_moment(moment):
+        if is_shear_moment(moment, method.dim):
             relaxation_rates.add(relax_info.relaxation_rate)
     if len(relaxation_rates) == 1:
         return relaxation_rates.pop()
@@ -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.
 
@@ -57,26 +69,3 @@ def relaxation_rate_scaling(omega, level_scale_factor):
         relaxation rate on refined grid
     """
     return omega / (omega / 2 + level_scale_factor * (1 - omega / 2))
-
-
-def default_relaxation_rate_names():
-    next_index = [0]
-
-    def result(moment_list):
-        shear_moment_inside = False
-        all_conserved_moments = True
-        for m in moment_list:
-            if is_shear_moment(m):
-                shear_moment_inside = True
-            if not (get_order(m) == 0 or get_order(m) == 1):
-                all_conserved_moments = False
-
-        if shear_moment_inside:
-            return sp.Symbol("omega")
-        elif all_conserved_moments:
-            return 0
-        else:
-            next_index[0] += 1
-            return sp.Symbol("omega_%d" % (next_index[0],))
-
-    return result

lbmpy/scenarios.py → src/lbmpy/scenarios.py

@@ -4,9 +4,9 @@ Scenario setup
 
 This module contains functions to set up pre-defined scenarios like a lid-driven cavity or channel flows.
 It is a good starting point if you are new to lbmpy.
-
+>>> from lbmpy.enums import Method
 >>> scenario = create_channel(domain_size=(20, 10), force=1e-5,
-...                          method='srt', relaxation_rate=1.9)
+...                          method=Method.SRT, relaxation_rate=1.9)
 >>> scenario.run(100)
 
 All scenarios can be modified, for example you can create a simple channel first, then place an object in it:
@@ -42,10 +42,10 @@ def create_fully_periodic_flow(initial_velocity, periodicity_in_kernel=False, lb
         initial_velocity: numpy array that defines an initial velocity for each cell. The shape of this
                          array determines the domain size.
         periodicity_in_kernel: don't use boundary handling for periodicity, but directly generate the kernel periodic
-        lbm_kernel: a LBM function, which would otherwise automatically created
+        lbm_kernel: an LBM function, which would otherwise automatically created
         data_handling: data handling instance that is used to create the necessary fields. If a data handling is
                        passed the periodicity and parallel arguments are ignored.
-        parallel: True for distributed memory parallelization with walberla
+        parallel: True for distributed memory parallelization with waLBerla
         kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
 
     Returns:
@@ -76,16 +76,22 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No
     Args:
         domain_size: tuple specifying the number of cells in each dimension
         lid_velocity: x velocity of lid in lattice coordinates.
-        lbm_kernel: a LBM function, which would otherwise automatically created
+        lbm_kernel: an LBM function, which would otherwise automatically created
         kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
-        parallel: True for distributed memory parallelization with walberla
+        parallel: True for distributed memory parallelization with waLBerla
         data_handling: see documentation of :func:`create_fully_periodic_flow`
     Returns:
         instance of :class:`Scenario`
     """
     assert domain_size is not None or data_handling is not None
     if data_handling is None:
-        data_handling = create_data_handling(domain_size, periodicity=False, default_ghost_layers=1, parallel=parallel)
+        optimization = kwargs.get('optimization', None)
+        target = optimization.get('target', None) if optimization else None
+        data_handling = create_data_handling(domain_size,
+                                             periodicity=False,
+                                             default_ghost_layers=1,
+                                             parallel=parallel,
+                                             default_target=target)
     step = LatticeBoltzmannStep(data_handling=data_handling, lbm_kernel=lbm_kernel, name="ldc", **kwargs)
 
     my_ubb = UBB(velocity=[lid_velocity, 0, 0][:step.method.dim])
@@ -99,7 +105,7 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No
 def create_channel(domain_size=None, force=None, pressure_difference=None, u_max=None, diameter_callback=None,
                    duct=False, wall_boundary=NoSlip(), parallel=False, data_handling=None, **kwargs):
     """Create a channel scenario (2D or 3D).
-    
+
     The channel can be driven either by force, velocity inflow or pressure difference. Choose one and pass
     exactly one of the parameters 'force', 'pressure_difference' or 'u_max'.
 

lbmpy/session.py → src/lbmpy/session.py

 import numpy as np
 import sympy as sp
 
-import lbmpy.plot as plt
 import pystencils as ps
+from pystencils import make_slice, show_code
+from pystencils.jupyter import *
+from lbmpy.advanced_streaming import *
 from lbmpy.boundaries import *
 from lbmpy.creationfunctions import *
+from lbmpy.enums import *
 from lbmpy.geometry import *
 from lbmpy.lbstep import LatticeBoltzmannStep
+from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
 from lbmpy.parameterization import ScalingWidget
+import lbmpy.plot as plt
 from lbmpy.postprocessing import *
 from lbmpy.scenarios import *
-from pystencils import make_slice, show_code
-from pystencils.jupyter import *
-from pystencils.sympy_gmpy_bug_workaround import *
+from lbmpy.stencils import LBStencil

src/lbmpy/simplificationfactory.py

+import sympy as sp
+
+from lbmpy.innerloopsplit import create_lbm_split_groups
+from lbmpy.methods.momentbased.momentbasedmethod import MomentBasedLbMethod
+from lbmpy.methods.momentbased.centralmomentbasedmethod import CentralMomentBasedLbMethod
+from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
+
+from lbmpy.methods.cumulantbased.cumulant_simplifications import (
+    insert_log_products, expand_post_collision_central_moments)
+from lbmpy.methods.momentbased.momentbasedsimplifications import (
+    factor_density_after_factoring_relaxation_times, factor_relaxation_rates,
+    replace_common_quadratic_and_constant_term, replace_density_and_velocity, replace_second_order_velocity_products,
+    insert_half_force, insert_conserved_quantity_products)
+from pystencils.simp import (
+    SimplificationStrategy, add_subexpressions_for_divisions, apply_to_all_assignments,
+    subexpression_substitution_in_main_assignments, insert_aliases, insert_constants,
+    add_subexpressions_for_constants)
+# add_subexpressions_for_constants)
+
+
+def create_simplification_strategy(lb_method, split_inner_loop=False):
+    if isinstance(lb_method, MomentBasedLbMethod):
+        if lb_method.moment_space_collision:
+            return _moment_space_simplification(split_inner_loop)
+        else:
+            if len(set(lb_method.relaxation_rates)) <= 2:
+                # SRT and TRT methods with population-space collision
+                return _srt_trt_population_space_simplification(split_inner_loop)
+            else:
+                # General MRT methods with population-space collision
+                return _mrt_population_space_simplification(split_inner_loop)
+    elif isinstance(lb_method, CentralMomentBasedLbMethod):
+        return _moment_space_simplification(split_inner_loop)
+    elif isinstance(lb_method, CumulantBasedLbMethod):
+        return _cumulant_space_simplification(split_inner_loop)
+    else:
+        return SimplificationStrategy()
+
+#   --------------- Internal ----------------------------------------------------------------------------
+
+
+def _srt_trt_population_space_simplification(split_inner_loop):
+    s = SimplificationStrategy()
+    expand = apply_to_all_assignments(sp.expand)
+    s.add(expand)
+    s.add(replace_second_order_velocity_products)
+    s.add(expand)
+    s.add(factor_relaxation_rates)
+    s.add(replace_density_and_velocity)
+    s.add(replace_common_quadratic_and_constant_term)
+    s.add(factor_density_after_factoring_relaxation_times)
+    s.add(subexpression_substitution_in_main_assignments)
+    if split_inner_loop:
+        s.add(create_lbm_split_groups)
+    s.add(add_subexpressions_for_divisions)
+    s.add(insert_constants)
+    s.add(insert_aliases)
+    s.add(lambda ac: ac.new_without_unused_subexpressions())
+    return s
+
+
+def _mrt_population_space_simplification(split_inner_loop):
+    s = SimplificationStrategy()
+    s.add(subexpression_substitution_in_main_assignments)
+    s.add(add_subexpressions_for_divisions)
+    if split_inner_loop:
+        s.add(create_lbm_split_groups)
+    s.add(lambda ac: ac.new_without_unused_subexpressions())
+    return s
+
+
+def _moment_space_simplification(split_inner_loop):
+    s = SimplificationStrategy()
+    s.add(insert_constants)
+    s.add(insert_half_force)
+    s.add(insert_aliases)
+    s.add(add_subexpressions_for_divisions)
+    s.add(add_subexpressions_for_constants)
+    if split_inner_loop:
+        s.add(create_lbm_split_groups)
+    s.add(lambda ac: ac.new_without_unused_subexpressions())
+    return s
+
+
+def _cumulant_space_simplification(split_inner_loop):
+    s = SimplificationStrategy()
+    s.add(insert_constants)
+    s.add(insert_aliases)
+    s.add(insert_log_products)
+    s.add(insert_conserved_quantity_products)
+    s.add(insert_half_force)
+    s.add(expand_post_collision_central_moments)
+    s.add(insert_aliases)
+    s.add(insert_constants)
+    s.add(insert_aliases)
+    s.add(add_subexpressions_for_divisions)
+    s.add(add_subexpressions_for_constants)
+    if split_inner_loop:
+        s.add(create_lbm_split_groups)
+    s.add(lambda ac: ac.new_without_unused_subexpressions())
+    return s

lbmpy/sparse/mapping.py → src/lbmpy/sparse/mapping.py

@@ -3,7 +3,7 @@ from typing import Tuple
 import numpy as np
 import sympy as sp
 
-from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
+from lbmpy.custom_code_nodes import LbmWeightInfo
 from pystencils import Assignment, Field, TypedSymbol
 from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
 from pystencils.boundaries.createindexlist import (
@@ -123,7 +123,10 @@ class SparseLbMapper:
                 assert direction_idx == 0
                 continue
             for own_cell_idx, cell in enumerate(self.fluid_coordinates):
-                inv_neighbor_cell = np.array([cell_i - dir_i for cell_i, dir_i in zip(cell, direction)])
+                inv_neighbor_cell = np.array(
+                    [np.int32(cell_i) - dir_i for cell_i, dir_i in zip(cell, direction)],
+                    dtype=np.uint32
+                )
                 if flag_arr[tuple(inv_neighbor_cell)] & fluid_boundary_mask:
                     neighbor_cell_idx = self.cell_idx(tuple(inv_neighbor_cell))
                     result.append(pdf_index(neighbor_cell_idx, direction_idx))

lbmpy/sparse/update_rule_sparse.py → src/lbmpy/sparse/update_rule_sparse.py

@@ -74,7 +74,7 @@ def create_macroscopic_value_setter_sparse(method, pdfs, density, velocity) -> A
         velocity: similar to density parameter
     """
     cqc = method.conserved_quantity_computation
-    inp_eqs = cqc.equilibrium_input_equations_from_init_values(density, velocity)
+    inp_eqs = cqc.equilibrium_input_equations_from_init_values(density, velocity, force_substitution=False)
     result = method.get_equilibrium(conserved_quantity_equations=inp_eqs)
     substitutions = {a: b for a, b in zip(method.post_collision_pdf_symbols, pdfs.center_vector)}
     return result.new_with_substitutions(substitutions).new_without_subexpressions()

src/lbmpy/stencils.py

+import pystencils as ps
+from pystencils.stencil import have_same_entries
+from lbmpy.enums import Stencil
+
+import sympy as sp
+
+
+class LBStencil:
+    r"""
+    Class representing a lattice Boltzmann stencil in DxQy notation, where d is the dimension
+    (length of the velocity tuples) and y is number of discrete velocities. For every dimension many different version
+    of a certain stencil is available. The reason for that is to ensure comparability with the literature.
+    Internally the stencil is represented as a tuple of tuples, where the ordering of the tuples plays no role.
+
+    Args:
+        stencil: Can be tuple of tuples which represents a DxQy stencil, a string like 'D2Q9' or an enum of
+                 lbmpy.enums.Stencil
+        ordering: the LBM literature does not use a common order of the discrete velocities, therefore here
+                  different common orderings are available. All orderings lead to the same method, it just has
+                  to be used consistently. Here more orderings are available to compare intermediate results with
+                  the literature.
+    """
+
+    def __init__(self, stencil, ordering='walberla'):
+        if isinstance(stencil, tuple):
+            ordering = None
+            self._stencil_entries = stencil
+        elif isinstance(stencil, str):
+            self._stencil_entries = _predefined_stencils(stencil, ordering)
+        elif isinstance(stencil, Stencil):
+            self._stencil_entries = _predefined_stencils(stencil.name, ordering)
+        else:
+            raise ValueError("The LBStencil can only be created with either a tuple of tuples which defines the "
+                             "stencil, a string or an Enum of type lbmpy.enums.Stencil")
+
+        valid_stencil = ps.stencil.is_valid(self._stencil_entries)
+        if valid_stencil is False:
+            raise ValueError("The stencil you have created is not valid. "
+                             "It probably contains elements with different lengths")
+
+        if len(set(self._stencil_entries)) < len(self._stencil_entries):
+            raise ValueError("The stencil you have created is not valid. "
+                             "It contains duplicated elements")
+
+        self._ordering = ordering
+        self._dim = len(self._stencil_entries[0])
+        self._q = len(self._stencil_entries)
+
+    @property
+    def D(self):
+        return self._dim
+
+    @property
+    def Q(self):
+        return self._q
+
+    @property
+    def ordering(self):
+        return self._ordering
+
+    @property
+    def name(self):
+        return f"D{self.D}Q{self.Q}"
+
+    @property
+    def stencil_entries(self):
+        return self._stencil_entries
+
+    @property
+    def inverse_stencil_entries(self):
+        return tuple([ps.stencil.inverse_direction(d) for d in self._stencil_entries])
+
+    def plot(self, slice=False, **kwargs):
+        ps.stencil.plot(stencil=self._stencil_entries, slice=slice, **kwargs)
+
+    def index(self, direction):
+        assert len(direction) == self.D, "direction must match stencil.D"
+        return self._stencil_entries.index(direction)
+
+    def inverse_index(self, direction):
+        assert len(direction) == self.D, "direction must match stencil.D"
+        direction = ps.stencil.inverse_direction(direction)
+        return self._stencil_entries.index(direction)
+
+    def __getitem__(self, index):
+        return self._stencil_entries[index]
+
+    def __iter__(self):
+        yield from self._stencil_entries
+
+    def __eq__(self, other):
+        return self.ordering == other.ordering and have_same_entries(self._stencil_entries, other.stencil_entries)
+
+    def __len__(self):
+        return len(self._stencil_entries)
+
+    def __str__(self):
+        return str(self._stencil_entries)
+
+    def __hash__(self):
+        return hash(self._stencil_entries)
+
+    def _repr_html_(self):
+        table = """
+        <table style="border:none; width: 100%">
+            <tr {nb}>
+                <th {nb} >Nr.</th>
+                <th {nb} >Direction Name</th>
+                <th {nb} >Direction </th>
+            </tr>
+            {content}
+        </table>
+        """
+        content = ""
+        for i, direction in enumerate(self._stencil_entries):
+            vals = {
+                'nr': sp.latex(i),
+                'name': sp.latex(ps.stencil.offset_to_direction_string(direction)),
+                'entry': sp.latex(direction),
+                'nb': 'style="border:none"'
+            }
+            content += """<tr {nb}>
+                            <td {nb}>${nr}$</td>
+                            <td {nb}>${name}$</td>
+                            <td {nb}>${entry}$</td>
+                         </tr>\n""".format(**vals)
+        return table.format(content=content, nb='style="border:none"')
+
+
+def _predefined_stencils(stencil: str, ordering: str):
+    predefined_stencils = {
+        'D2Q9': {
+            'walberla': ((0, 0),
+                         (0, 1), (0, -1), (-1, 0), (1, 0),
+                         (-1, 1), (1, 1), (-1, -1), (1, -1),),
+            'counterclockwise': ((0, 0),
+                                 (1, 0), (0, 1), (-1, 0), (0, -1),
+                                 (1, 1), (-1, 1), (-1, -1), (1, -1)),
+            'braunschweig': ((0, 0),
+                             (-1, 1), (-1, 0), (-1, -1), (0, -1),
+                             (1, -1), (1, 0), (1, 1), (0, 1)),
+            'uk': ((0, 0),
+                   (1, 0), (-1, 0), (0, 1), (0, -1),
+                   (1, 1), (-1, -1), (-1, 1), (1, -1),
+                   ),
+            'lehmann': ((0, 0),
+                        (1, 0), (-1, 0), (0, 1), (0, -1),
+                        (1, 1), (-1, -1), (1, -1), (-1, 1),
+                        )
+        },
+        'D2V17': {
+            'walberla': (
+                (0, 0), (0, -1), (-1, 0), (1, 0), (0, 1), (-1, -1), (1, -1), (-1, 1), (1, 1), (-2, -2), (2, -2),
+                (-2, 2),
+                (2, 2), (0, -3), (-3, 0), (3, 0), (0, 3)),
+        },
+        'D2V37': {
+            'walberla': (
+                (0, 0), (0, -1), (-1, 0), (1, 0), (0, 1), (-1, -1), (1, -1), (-1, 1), (1, 1), (0, -2), (-2, 0),
+                (2, 0),
+                (0, 2), (-1, -2), (1, -2), (-2, -1), (2, -1), (-2, 1), (2, 1), (-1, 2), (1, 2), (-2, -2), (2, -2),
+                (-2, 2),
+                (2, 2), (0, -3), (-3, 0), (3, 0), (0, 3), (-1, -3), (1, -3), (-3, -1), (3, -1), (-3, 1), (3, 1),
+                (-1, 3),
+                (1, 3))
+        },
+        'D3Q7': {
+            'walberla': ((0, 0, 0),
+                         (0, 1, 0), (0, -1, 0),
+                         (-1, 0, 0), (1, 0, 0),
+                         (0, 0, 1), (0, 0, -1))
+
+        },
+        'D3Q15': {
+            'walberla':
+                ((0, 0, 0),
+                 (0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
+                 (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1),
+                 (-1, -1, -1)),
+            'premnath': ((0, 0, 0),
+                         (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                         (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
+                         (1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
+            'lehmann': ((0, 0, 0),
+                        (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                        (1, 1, 1), (-1, -1, -1), (1, 1, -1), (-1, -1, 1),
+                        (1, -1, 1), (-1, 1, -1), (-1, 1, 1), (1, -1, -1)),
+            'fakhari': ((0, 0, 0),
+                        (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                        (1, 1, 1), (-1, -1, -1), (-1, 1, 1), (1, -1, -1),
+                        (1, -1, 1), (-1, 1, -1), (1, 1, -1), (-1, -1, 1)),
+        },
+        'D3Q19': {
+            'walberla': ((0, 0, 0),
+                         (0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
+                         (-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
+                         (0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
+                         (0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1)),
+            'counterclockwise': ((0, 0, 0),
+                                 (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                                 (1, 1, 0), (-1, -1, 0), (1, 0, 1), (-1, 0, -1),
+                                 (0, 1, 1), (0, -1, -1), (1, -1, 0), (-1, 1, 0),
+                                 (1, 0, -1), (-1, 0, 1), (0, 1, -1), (0, -1, 1)),
+            'braunschweig': ((0, 0, 0),
+                             (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                             (1, 1, 0), (-1, -1, 0), (1, -1, 0), (-1, 1, 0),
+                             (1, 0, 1), (-1, 0, -1), (1, 0, -1), (-1, 0, 1),
+                             (0, 1, 1), (0, -1, -1), (0, 1, -1), (0, -1, 1)),
+            'premnath': ((0, 0, 0),
+                         (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                         (1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
+                         (1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1),
+                         (0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1)),
+            'lehmann': ((0, 0, 0),
+                        (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                        (1, 1, 0), (-1, -1, 0), (1, 0, 1), (-1, 0, -1),
+                        (0, 1, 1), (0, -1, -1), (1, -1, 0), (-1, 1, 0),
+                        (1, 0, -1), (-1, 0, 1), (0, 1, -1), (0, -1, 1)),
+        },
+        'D3Q27': {
+            'walberla': ((0, 0, 0),
+                         (0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
+                         (-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
+                         (0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
+                         (0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1),
+                         (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1),
+                         (-1, -1, -1)),
+            'premnath': ((0, 0, 0),
+                         (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                         (1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
+                         (1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1),
+                         (0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1),
+                         (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
+                         (1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
+            'fakhari': ((0, 0, 0),
+                        (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                        (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
+                        (1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1),
+                        (1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
+                        (1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1), (0, 1, 1), (0, -1, 1), (0, 1, -1),
+                        (0, -1, -1)),
+            'lehmann': ((0, 0, 0),
+                        (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                        (1, 1, 0), (-1, -1, 0), (1, 0, 1), (-1, 0, -1),
+                        (0, 1, 1), (0, -1, -1), (1, -1, 0), (-1, 1, 0),
+                        (1, 0, -1), (-1, 0, 1), (0, 1, -1), (0, -1, 1),
+                        (1, 1, 1), (-1, -1, -1), (1, 1, -1), (-1, -1, 1), (1, -1, 1), (-1, 1, -1), (-1, 1, 1),
+                        (1, -1, -1)),
+            'braunschweig': ((0, 0, 0),
+                             (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
+                             (1, 1, 0), (-1, -1, 0), (1, -1, 0), (-1, 1, 0),
+                             (1, 0, 1), (-1, 0, -1), (1, 0, -1), (-1, 0, 1),
+                             (0, 1, 1), (0, -1, -1), (0, 1, -1), (0, -1, 1),
+                             (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1),
+                             (-1, -1, -1)),
+        }
+    }
+    try:
+        return predefined_stencils[stencil][ordering]
+    except KeyError:
+        err_msg = ""
+        for stencil, ordering_names in predefined_stencils.items():
+            err_msg += "  %s: %s\n" % (stencil, ", ".join(ordering_names.keys()))
+
+        raise ValueError("No such stencil available. "
+                         "Available stencils: <stencil_name>( <ordering_names> )\n" + err_msg)