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

lbmpy/phasefield_allen_cahn/contact_angle.py → src/lbmpy/phasefield_allen_cahn/contact_angle.py

 import math
 import sympy as sp
 
-from pystencils.astnodes import SympyAssignment
+from pystencils.astnodes import Block, Conditional, SympyAssignment
 
 from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
 from pystencils.boundaries.boundaryconditions import Boundary
@@ -34,26 +34,29 @@ class ContactAngle(Boundary):
     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 math.isclose(90, self._contact_angle, abs_tol=1e-5):
-                return [SympyAssignment(field.center, field[neighbor])]
+            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
 
-            dist = TypedSymbol("h", self._data_type)
-            angle = TypedSymbol("a", self._data_type)
-            tmp = TypedSymbol("tmp", self._data_type)
-
-            result = [SympyAssignment(tmp, CastFunc(sum([x * x for x in neighbor]), self._data_type)),
-                      SympyAssignment(dist, 0.5 * sp.sqrt(tmp)),
-                      SympyAssignment(angle, math.cos(math.radians(self._contact_angle)))]
-
-            var = - dist * (4.0 / self._interface_width) * angle
-            tmp = 1 + var
-            else_branch = (tmp - sp.sqrt(tmp * tmp - 4 * var * field[neighbor])) / var - field[neighbor]
-            update = sp.Piecewise((field[neighbor], dist < 0.001), (else_branch, True))
-
-            result.append(SympyAssignment(field.center, update))
-            return result
         else:
             raise NotImplementedError("Contact angle only implemented for phase-fields which have a single "
                                       "value for each cell")

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

lbmpy/phasefield_allen_cahn/kernel_equations.py → src/lbmpy/phasefield_allen_cahn/kernel_equations.py

-from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
+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.parameter_calculation import AllenCahnParameters
+
+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
 
@@ -18,29 +25,7 @@ def chemical_potential_symbolic(phi_field, stencil, beta, kappa):
         beta: coefficient related to surface tension and interface thickness
         kappa: coefficient related to surface tension and interface thickness
     """
-    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.Q == 9:
-            res = deriv.get_stencil(isotropic=True)
-            lap += res.apply(phi_field.center)
-        elif stencil.Q == 15:
-            deriv.set_weight((0, 0, 0), sp.Rational(-32, 27))
-            res = deriv.get_stencil(isotropic=True)
-            lap += res.apply(phi_field.center)
-        elif stencil.Q == 19:
-            res = deriv.get_stencil(isotropic=True)
-            lap += res.apply(phi_field.center)
-        else:
-            deriv.set_weight((0, 0, 0), sp.Rational(-38, 27))
-            res = deriv.get_stencil(isotropic=True)
-            lap += res.apply(phi_field.center)
-
+    lap = laplacian_symbolic(phi_field, stencil)
     # get the chemical potential
     four = sp.Rational(4, 1)
     one = sp.Rational(1, 1)
@@ -49,50 +34,6 @@ def chemical_potential_symbolic(phi_field, stencil, beta, kappa):
     return mu
 
 
-def isotropic_gradient_symbolic(phi_field, stencil):
-    r"""
-    Get a symbolic expression for the isotropic gradient of the phase-field
-    Args:
-        phi_field: the phase-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.Q == 9:
-        res = deriv.get_stencil(isotropic=True)
-        grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), 0]
-    elif stencil.Q == 15:
-        res = deriv.get_stencil(isotropic=True)
-        grad = [res.apply(phi_field.center),
-                res.rotate_weights_and_apply(phi_field.center, (0, 1)),
-                res.rotate_weights_and_apply(phi_field.center, (1, 2))]
-    elif 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(phi_field.center),
-                res.rotate_weights_and_apply(phi_field.center, (0, 1)),
-                res.rotate_weights_and_apply(phi_field.center, (1, 2))]
-    else:
-        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(phi_field.center),
-                res.rotate_weights_and_apply(phi_field.center, (0, 1)),
-                res.rotate_weights_and_apply(phi_field.center, (1, 2))]
-
-    return grad
-
-
 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
@@ -134,11 +75,10 @@ def pressure_force(phi_field, lb_method, stencil, density_heavy, density_light,
     return result
 
 
-def viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, density_light, fd_stencil=None):
+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:
-        lb_velocity_field: hydrodynamic distribution function
         phi_field: phase-field
         lb_method: lattice boltzmann method used for hydrodynamics
         tau: relaxation time of the hydrodynamic lattice boltzmann step
@@ -154,7 +94,7 @@ def viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, d
 
     iso_grad = sp.Matrix(isotropic_gradient_symbolic(phi_field, fd_stencil)[:stencil.D])
 
-    f_neq = lb_velocity_field.center_vector - lb_method.get_equilibrium_terms()
+    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
 
@@ -169,17 +109,19 @@ def viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, d
     return [fmx, fmy, fmz]
 
 
-def surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil=None):
+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
-        beta: coefficient related to surface tension and interface thickness
-        kappa: coefficient related to surface tension and interface thickness
+        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
 
@@ -188,18 +130,57 @@ def surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil=None):
     return [chemical_potential * x for x in iso_grad]
 
 
-def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, parameters: AllenCahnParameters,
-                       body_force, fd_stencil=None):
+def thermocapillary_surface_tension_force(phi_field, temperature_field,
+                                          stencil, parameters, fd_stencil=None):
     r"""
-    Get a symbolic expression for the hydrodynamic force
+    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:
-        lb_velocity_field: hydrodynamic distribution function
         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
 
@@ -211,12 +192,16 @@ def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, parameters: Alle
     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)
-    beta = parameters.beta
-    kappa = parameters.kappa
 
     fp = pressure_force(phi_field, lb_method, stencil, density_heavy, density_light, fd_stencil)
-    fm = viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, density_light, fd_stencil)
-    fs = surface_tension_force(phi_field, stencil, beta, kappa, 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):
@@ -254,14 +239,14 @@ def interface_tracking_force(phi_field, stencil, parameters: AllenCahnParameters
     return result
 
 
-def hydrodynamic_force_assignments(lb_velocity_field, velocity_field, phi_field, lb_method,
+def hydrodynamic_force_assignments(velocity_field, phi_field, lb_method,
                                    parameters: AllenCahnParameters,
-                                   body_force, fd_stencil=None, sub_iterations=2):
+                                   body_force, fd_stencil=None, sub_iterations=2,
+                                   temperature_field=None):
 
     r"""
     Get a symbolic expression for the hydrodynamic force
     Args:
-        lb_velocity_field: hydrodynamic distribution function
         velocity_field: velocity
         phi_field: phase-field
         lb_method: Lattice boltzmann method used for hydrodynamics
@@ -270,6 +255,7 @@ def hydrodynamic_force_assignments(lb_velocity_field, velocity_field, phi_field,
         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
@@ -280,12 +266,13 @@ def hydrodynamic_force_assignments(lb_velocity_field, velocity_field, phi_field,
     # method has to have a force model
     symbolic_force = lb_method.force_model.symbolic_force_vector
 
-    force = hydrodynamic_force(lb_velocity_field, phi_field, lb_method, parameters, body_force, fd_stencil=fd_stencil)
+    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_velocity_field.center_vector)
+    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:]]
@@ -332,7 +319,7 @@ def add_interface_tracking_force(update_rule: LbmCollisionRule, force):
 
 
 def add_hydrodynamic_force(update_rule: LbmCollisionRule, force, phi_field,
-                           hydro_pdfs, parameters: AllenCahnParameters):
+                           hydro_pdfs, parameters: AllenCahnParameters, lbm_config: LBMConfig = None):
     r"""
      Adds the interface tracking force to a lattice Boltzmann update rule
      Args:
@@ -341,29 +328,53 @@ def add_hydrodynamic_force(update_rule: LbmCollisionRule, force, phi_field,
          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)
 
-    method = update_rule.method
-    symbolic_force = method.force_model.symbolic_force_vector
+    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
-    rho = cqc.density_deviation_symbol
+    density = cqc.density_symbol
+    density_deviation = cqc.density_deviation_symbol
+    free_symbols = update_rule.free_symbols
 
-    force_subs = {f: f / density for f in symbolic_force}
+    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)
 
-    update_rule = update_rule.subs(force_subs)
+    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.subexpressions += [Assignment(rho, sum(hydro_pdfs.center_vector))]
-    update_rule.subexpressions += force
     update_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
 
-    return update_rule
+    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):
+                                      fd_stencil=None, **kwargs):
     r"""
     Returns an assignment list for initializing the phase-field distribution functions
     Args:
@@ -374,9 +385,10 @@ def initializer_kernel_phase_field_lb(lb_method, phi, velocity, ac_pdfs, paramet
         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)
+    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)
 
@@ -407,7 +419,7 @@ def initializer_kernel_phase_field_lb(lb_method, phi, velocity, ac_pdfs, paramet
     return h_updates
 
 
-def initializer_kernel_hydro_lb(lb_method, pressure, velocity, hydro_pdfs):
+def initializer_kernel_hydro_lb(lb_method, pressure, velocity, hydro_pdfs, **kwargs):
     r"""
     Returns an assignment list for initializing the velocity distribution functions
     Args:
@@ -415,11 +427,12 @@ def initializer_kernel_hydro_lb(lb_method, pressure, velocity, hydro_pdfs):
         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)
+    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

lbmpy/phasefield_allen_cahn/parameter_calculation.py → src/lbmpy/phasefield_allen_cahn/parameter_calculation.py

@@ -103,7 +103,8 @@ class AllenCahnParameters:
     def symbolic_to_numeric_map(self):
         return {t.name: self.parameter_map()[t] for t in self.parameter_map()}
 
-    def _repr_html_(self):
+    @staticmethod
+    def _parameter_strings():
         names = ("Density heavy phase",
                  "Density light phase",
                  "Relaxation time heavy phase",
@@ -113,7 +114,10 @@ class AllenCahnParameters:
                  "Interface thickness",
                  "Mobility",
                  "Surface tension")
+        return names
 
+    def _repr_html_(self):
+        names = self._parameter_strings()
         table = """
         <table style="border:none; width: 100%">
             <tr {nb}>
@@ -140,6 +144,85 @@ class AllenCahnParameters:
         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,
@@ -170,8 +253,8 @@ def calculate_parameters_rti(reference_length=256,
     """
 
     # calculate the gravitational acceleration and the reference velocity
-    g = reference_length / (reference_time ** 2 * atwood_number)
-    reference_velocity = math.sqrt(g * reference_length)
+    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
@@ -187,7 +270,7 @@ def calculate_parameters_rti(reference_length=256,
                                      dynamic_viscosity_light=dynamic_viscosity_light,
                                      surface_tension=surface_tension,
                                      mobility=mobility,
-                                     gravitational_acceleration=-g)
+                                     gravitational_acceleration=-gravity_lattice_units)
     return parameters
 
 
@@ -215,15 +298,15 @@ def calculate_dimensionless_rising_bubble(reference_time=18000,
         viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
     """
 
-    bubble_diameter = bubble_radius * 2
-    g = bubble_diameter / (reference_time ** 2)
+    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(g * bubble_diameter ** 3)) / reynolds_number
+    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) * g * bubble_diameter ** 2 / bond_number
+    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)
 
@@ -232,5 +315,121 @@ def calculate_dimensionless_rising_bubble(reference_time=18000,
                                      dynamic_viscosity_heavy=dynamic_viscosity_heavy,
                                      dynamic_viscosity_light=dynamic_viscosity_light,
                                      surface_tension=surface_tension,
-                                     gravitational_acceleration=-g)
+                                     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/scenarios.py → src/lbmpy/scenarios.py

@@ -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,9 +76,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`

lbmpy/simplificationfactory.py → src/lbmpy/simplificationfactory.py

@@ -92,6 +92,7 @@ def _cumulant_space_simplification(split_inner_loop):
     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:

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))