lbmpy/phasefield/phasefieldstep.py → src/lbmpy/phasefield/phasefieldstep.py

@@ -12,6 +12,7 @@ from lbmpy.phasefield.kerneleqs import (
 from pystencils import create_data_handling, create_kernel
 from pystencils.boundaries.boundaryhandling import FlagInterface
 from pystencils.boundaries.inkernel import add_neumann_boundary
+from pystencils.enums import Target
 from pystencils.simp import sympy_cse_on_assignment_list
 from pystencils.slicing import SlicedGetter, make_slice
 
@@ -34,9 +35,9 @@ class PhaseFieldStep:
                  discretization='standard'):
 
         if optimization is None:
-            optimization = {'openmp': False, 'target': 'cpu'}
+            optimization = {'openmp': False, 'target': Target.CPU}
         openmp = optimization.get('openmp', False)
-        target = optimization.get('target', 'cpu')
+        target = optimization.get('target', Target.CPU)
 
         if data_handling is None:
             data_handling = create_data_handling(domain_size, periodicity=True, parallel=False)
@@ -57,7 +58,7 @@ class PhaseFieldStep:
         self.chemical_potentials = chemical_potentials_from_free_energy(free_energy, order_parameters)
 
         # ------------------ Adding arrays ---------------------
-        gpu = target == 'gpu'
+        gpu = target == Target.GPU
         self.gpu = gpu
         self.num_order_parameters = len(order_parameters)
         self.order_parameters = order_parameters
@@ -134,7 +135,8 @@ class PhaseFieldStep:
 
         self.hydro_lbm_step = LatticeBoltzmannStep(data_handling=data_handling, name=name + '_hydroLBM',
                                                    relaxation_rate=hydro_dynamic_relaxation_rate,
-                                                   compute_velocity_in_every_step=True, force=self.force_field,
+                                                   compute_velocity_in_every_step=True,
+                                                   force=tuple([self.force_field(i) for i in range(dh.dim)]),
                                                    velocity_data_name=self.vel_field_name, kernel_params=kernel_params,
                                                    flag_interface=self.flag_interface,
                                                    time_step_order='collide_stream',
@@ -172,7 +174,7 @@ class PhaseFieldStep:
                                                compute_density_in_every_step=True,
                                                density_data_index=i,
                                                flag_interface=self.hydro_lbm_step.boundary_handling.flag_interface,
-                                               name=name + "_chLbm_%d" % (i,),
+                                               name=name + f"_chLbm_{i}",
                                                optimization=optimization)
                 self.cahn_hilliard_steps.append(ch_step)
 

lbmpy/phasefield/phasefieldstep_direct.py → src/lbmpy/phasefield/phasefieldstep_direct.py

 import numpy as np
 
 from lbmpy.creationfunctions import create_lb_function
+from lbmpy.enums import Stencil
 from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
 from lbmpy.phasefield.analytical import force_from_phi_and_mu
 from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
 from lbmpy.phasefield.kerneleqs import mu_kernel
-from lbmpy.stencils import get_stencil
-from pystencils import Assignment, create_data_handling, create_kernel
+from lbmpy.stencils import LBStencil
+from pystencils import Assignment, create_data_handling, create_kernel, Target
 from pystencils.fd import discretize_spatial
 from pystencils.fd.spatial import fd_stencils_forth_order_isotropic
 from pystencils.simp import sympy_cse_on_assignment_list
@@ -21,9 +22,9 @@ class PhaseFieldStepDirect:
                  cahn_hilliard_gammas=1,
                  optimization=None):
         if optimization is None:
-            optimization = {'openmp': False, 'target': 'cpu'}
+            optimization = {'openmp': False, 'target': Target.CPU}
         openmp = optimization.get('openmp', False)
-        target = optimization.get('target', 'cpu')
+        target = optimization.get('target', Target.CPU)
 
         if not hasattr(cahn_hilliard_relaxation_rates, '__len__'):
             cahn_hilliard_relaxation_rates = [cahn_hilliard_relaxation_rates] * len(order_parameters)
@@ -36,23 +37,23 @@ class PhaseFieldStepDirect:
         dh = data_handling
         self.data_handling = dh
 
-        stencil = get_stencil('D3Q19' if dh.dim == 3 else 'D2Q9')
+        stencil = LBStencil(Stencil.D3Q19) if dh.dim == 3 else LBStencil(Stencil.D2Q9)
 
         self.free_energy = free_energy
 
         # Data Handling
         kernel_parameters = {'cpu_openmp': openmp, 'target': target, 'ghost_layers': 2}
-        gpu = target == 'gpu'
+        gpu = target == Target.GPU
         phi_size = len(order_parameters)
         gl = kernel_parameters['ghost_layers']
-        self.phi_field = dh.add_array('{}_phi'.format(name), values_per_cell=phi_size, gpu=gpu, latex_name='φ',
+        self.phi_field = dh.add_array(f'{name}_phi', values_per_cell=phi_size, gpu=gpu, latex_name='φ',
                                       ghost_layers=gl)
-        self.mu_field = dh.add_array("{}_mu".format(name), values_per_cell=phi_size, gpu=gpu, latex_name="μ",
+        self.mu_field = dh.add_array(f"{name}_mu", values_per_cell=phi_size, gpu=gpu, latex_name="μ",
                                      ghost_layers=gl)
-        self.vel_field = dh.add_array("{}_u".format(name), values_per_cell=data_handling.dim, gpu=gpu, latex_name="u",
+        self.vel_field = dh.add_array(f"{name}_u", values_per_cell=data_handling.dim, gpu=gpu, latex_name="u",
                                       ghost_layers=gl)
 
-        self.force_field = dh.add_array("{}_force".format(name), values_per_cell=dh.dim, gpu=gpu, latex_name="F",
+        self.force_field = dh.add_array(f"{name}_force", values_per_cell=dh.dim, gpu=gpu, latex_name="F",
                                         ghost_layers=gl)
 
         self.phi = SlicedGetterDataHandling(self.data_handling, self.phi_field.name)
@@ -62,11 +63,11 @@ class PhaseFieldStepDirect:
 
         self.ch_pdfs = []
         for i in range(len(order_parameters)):
-            src = dh.add_array("{}_ch_src{}".format(name, i), values_per_cell=len(stencil), ghost_layers=gl)
-            dst = dh.add_array("{}_ch_dst{}".format(name, i), values_per_cell=len(stencil), ghost_layers=gl)
+            src = dh.add_array(f"{name}_ch_src{i}", values_per_cell=stencil.Q, ghost_layers=gl)
+            dst = dh.add_array(f"{name}_ch_dst{i}", values_per_cell=stencil.Q, ghost_layers=gl)
             self.ch_pdfs.append((src, dst))
-        self.hydro_pdfs = (dh.add_array("{}_hydro_src".format(name), values_per_cell=len(stencil), ghost_layers=gl),
-                           dh.add_array("{}_hydro_dst".format(name), values_per_cell=len(stencil), ghost_layers=gl))
+        self.hydro_pdfs = (dh.add_array(f"{name}_hydro_src", values_per_cell=stencil.Q, ghost_layers=gl),
+                           dh.add_array(f"{name}_hydro_dst", values_per_cell=stencil.Q, ghost_layers=gl))
 
         # Compute Kernels
         mu_assignments = mu_kernel(self.free_energy, order_parameters, self.phi_field, self.mu_field)

lbmpy/phasefield/post_processing.py → src/lbmpy/phasefield/post_processing.py

@@ -135,8 +135,7 @@ def get_triple_points(phase_arr, phase_indices, contour_line_eps=0.01, threshold
 def analytic_neumann_angles(kappas):
     """Computes analytic Neumann angles using surface tension parameters.
 
-    >>> analytic_neumann_angles([0.1, 0.1, 0.1])
-    [120.00000000000001, 120.00000000000001, 120.00000000000001]
+    >>> assert analytic_neumann_angles([0.1, 0.1, 0.1]) == [120.00000000000001, 120.00000000000001, 120.00000000000001]
     >>> r = analytic_neumann_angles([0.1, 0.2, 0.3])
     >>> assert np.allclose(sum(r), 360)
 

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 .._compat import IS_PYSTENCILS_2
+
+if IS_PYSTENCILS_2:
+    raise ImportError("`lbmpy.phasefield_allen_cahn.contact_angle` is only available when running with pystencils 1.x")
+
+from pystencils.astnodes import Block, Conditional, SympyAssignment
+
+from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
+from pystencils.boundaries.boundaryconditions import Boundary
+
+from pystencils 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 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/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 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.
 

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:
@@ -27,9 +27,11 @@ that defines the viscosity of the fluid (valid values being between 0 and 2).
 """
 import numpy as np
 
+from lbmpy._compat import IS_PYSTENCILS_2
 from lbmpy.boundaries import UBB, FixedDensity, NoSlip
 from lbmpy.geometry import add_pipe_inflow_boundary, add_pipe_walls
 from lbmpy.lbstep import LatticeBoltzmannStep
+from pystencils import Target
 from pystencils.datahandling import create_data_handling
 from pystencils.slicing import slice_from_direction
 
@@ -42,10 +44,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,9 +78,9 @@ 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`
@@ -86,7 +88,7 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No
     assert domain_size is not None or data_handling is not None
     if data_handling is None:
         optimization = kwargs.get('optimization', None)
-        target = optimization.get('target', None) if optimization else None
+        target = optimization.get('target', None) if optimization else Target.CPU
         data_handling = create_data_handling(domain_size,
                                              periodicity=False,
                                              default_ghost_layers=1,
@@ -94,7 +96,13 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No
                                              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])
+    if IS_PYSTENCILS_2:
+        my_ubb = UBB(
+            velocity=[lid_velocity, 0, 0][:step.method.dim],
+            data_type=step.config.get_option("default_dtype")
+        )
+    else:  # legacy version does not set data type of boundary correctly
+        my_ubb = UBB(velocity=[lid_velocity, 0, 0][:step.method.dim])
     step.boundary_handling.set_boundary(my_ubb, slice_from_direction('N', step.dim))
     for direction in ('W', 'E', 'S') if step.dim == 2 else ('W', 'E', 'S', 'T', 'B'):
         step.boundary_handling.set_boundary(NoSlip(), slice_from_direction(direction, step.dim))