lbmpy/phasefield/contact_angle_circle_fitting.py → src/lbmpy/phasefield/contact_angle_circle_fitting.py

@@ -39,8 +39,8 @@ def circle_intersections(midpoint0, midpoint1, radius0, radius1):
 def interface_region(concentration_arr, phase0, phase1, area=3):
     import scipy.ndimage as sc_image
 
-    area_phase0 = sc_image.morphology.binary_dilation(concentration_arr[..., phase0] > 0.5, iterations=area)
-    area_phase1 = sc_image.morphology.binary_dilation(concentration_arr[..., phase1] > 0.5, iterations=area)
+    area_phase0 = sc_image.binary_dilation(concentration_arr[..., phase0] > 0.5, iterations=area)
+    area_phase1 = sc_image.binary_dilation(concentration_arr[..., phase1] > 0.5, iterations=area)
     return np.logical_and(area_phase0, area_phase1)
 
 

lbmpy/phasefield/nphase_nestler.py → src/lbmpy/phasefield/nphase_nestler.py

@@ -4,10 +4,7 @@ try:
     pyximport.install(language_level=3)
     from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
 except ImportError:
-    try:
-        from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
-    except ImportError:
-        raise ImportError("neither pyximport nor binary module simplex_projection_2d available.")
+    raise ImportError("pyximport not available. Please install Cython to use simplex_projection_2d.")
 
 import sympy as sp
 

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

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 .._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.typing import TypedSymbol
+from pystencils import TypedSymbol
 from pystencils.typing import CastFunc
 
 
@@ -34,26 +39,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])]
-
-            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))
+            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
 
-            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")