lbmpy/oldroydb.py → src/lbmpy/oldroydb.py

@@ -13,7 +13,7 @@ class OldroydB:
         assert not ps.FieldType.is_staggered(F)
         assert ps.FieldType.is_staggered(tauflux)
         assert ps.FieldType.is_staggered(tauface)
-        
+
         self.dim = dim
         self.u = u
         self.tau = tau
@@ -22,7 +22,7 @@ class OldroydB:
         self.tauface_field = tauface
         self.lambda_p = lambda_p
         self.eta_p = eta_p
-        
+
         full_stencil = ["C"] + self.tauflux.staggered_stencil + \
             list(map(inverse_direction_string, self.tauflux.staggered_stencil))
         self.stencil = tuple(map(lambda d: tuple(ps.stencil.direction_string_to_offset(d, self.dim)), full_stencil))
@@ -30,27 +30,27 @@ class OldroydB:
             list(map(inverse_direction_string, self.tauface_field.staggered_stencil))
         self.force_stencil = tuple(map(lambda d: tuple(ps.stencil.direction_string_to_offset(d, self.dim)),
                                    full_stencil))
-        
+
         self.disc = ps.fd.FVM1stOrder(self.tau, self._flux(), self._source())
         if vof:
             self.vof = ps.fd.VOF(self.tauflux, self.u, self.tau)
         else:
             self.vof = None
-    
+
     def _flux(self):
         return [self.tau.center_vector.applyfunc(lambda t: t * self.u.center_vector[i]) for i in range(self.dim)]
-    
+
     def _source(self):
         gradu = sp.Matrix([[ps.fd.diff(self.u.center_vector[j], i) for j in range(self.dim)] for i in range(self.dim)])
         gamma = gradu + gradu.transpose()
         return self.tau.center_vector * gradu + gradu.transpose() * self.tau.center_vector + \
             (self.eta_p * gamma - self.tau.center_vector) / self.lambda_p
-    
+
     def tauface(self):
         return ps.AssignmentCollection([ps.Assignment(self.tauface_field.staggered_vector_access(d),
                                         (self.tau.center_vector + self.tau.neighbor_vector(d)) / 2)
                                         for d in self.tauface_field.staggered_stencil])
-    
+
     def force(self):
         full_stencil = self.tauface_field.staggered_stencil + \
             list(map(inverse_direction_string, self.tauface_field.staggered_stencil))
@@ -60,13 +60,13 @@ class OldroydB:
                          for d in full_stencil]) for j in range(self.dim)])
         A0 = sum([sp.Matrix(direction_string_to_offset(d)).norm() for d in full_stencil])
         return ps.AssignmentCollection(ps.Assignment(self.F.center_vector, dtau / A0 * 2 * self.dim))
-    
+
     def flux(self):
         if self.vof is not None:
             return self.vof
         else:
             return self.disc.discrete_flux(self.tauflux)
-    
+
     def continuity(self):
         cont = self.disc.discrete_continuity(self.tauflux)
         tau_copy = sp.Matrix(self.dim, self.dim, lambda i, j: sp.Symbol("tau_old_%d_%d" % (i, j)))
@@ -79,24 +79,24 @@ class OldroydB:
 class Flux(Boundary):
     inner_or_boundary = True  # call the boundary condition with the fluid cell
     single_link = False  # needs to be called for all directional fluxes
-    
+
     def __init__(self, stencil, value=None):
         self.stencil = stencil
         self.value = value
 
     def __call__(self, field, direction_symbol, **kwargs):
         assert ps.FieldType.is_staggered(field)
-        
+
         assert all([s == 0 for s in self.stencil[0]])
         accesses = [field.staggered_vector_access(ps.stencil.offset_to_direction_string(d))
                     for d in self.stencil[1:]]
         conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))]
-        
+
         if self.value is None:
             val = sp.Matrix(np.zeros(accesses[0].shape, dtype=int))
         else:
             val = self.value
-        
+
         # use conditional
         conditional = None
         for a, c, d in zip(accesses, conds, self.stencil[1:]):
@@ -118,7 +118,7 @@ class Flux(Boundary):
         return hash((Flux, self.stencil, self.value))
 
     def __eq__(self, other):
-        return type(other) == Flux and other.stencil == self.stencil and self.value == other.value
+        return type(other) is Flux and other.stencil == self.stencil and self.value == other.value
 
 
 class Extrapolation(Boundary):
@@ -139,12 +139,12 @@ class Extrapolation(Boundary):
 
     def __call__(self, field, direction_symbol, **kwargs):
         assert ps.FieldType.is_staggered(field)
-        
+
         assert all([s == 0 for s in self.stencil[0]])
         accesses = [field.staggered_vector_access(ps.stencil.offset_to_direction_string(d))
                     for d in self.stencil[1:]]
         conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))]
-        
+
         # use conditional
         conditional = None
         for a, c, o in zip(accesses, conds, self.stencil[1:]):
@@ -169,28 +169,29 @@ class Extrapolation(Boundary):
         return hash((Extrapolation, self.stencil, self.src, self.weights))
 
     def __eq__(self, other):
-        return type(other) == Extrapolation and other.stencil == self.stencil and \
+        return type(other) is Extrapolation and other.stencil == self.stencil and \
             other.src == self.src and other.weights == self.weights
 
 
 class ForceOnBoundary(Boundary):
     inner_or_boundary = False  # call the boundary condition with the boundary cell
     single_link = False  # needs to be called for all directional fluxes
-    
-    def __init__(self, stencil, force_field):
+
+    def __init__(self, stencil, force_field, name=None):
         self.stencil = stencil
         self.force_field = force_field
-        
+        super(ForceOnBoundary).__init__(name)
+
         assert not ps.FieldType.is_staggered(force_field)
 
     def __call__(self, face_stress_field, direction_symbol, **kwargs):
         assert ps.FieldType.is_staggered(face_stress_field)
-        
+
         assert all([s == 0 for s in self.stencil[0]])
         accesses = [face_stress_field.staggered_vector_access(ps.stencil.offset_to_direction_string(d))
                     for d in self.stencil[1:]]
         conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))]
-        
+
         # use conditional
         conditional = None
         for a, c, o in zip(accesses, conds, self.stencil[1:]):
@@ -205,5 +206,5 @@ class ForceOnBoundary(Boundary):
         return hash((ForceOnBoundary, self.stencil, self.force_field))
 
     def __eq__(self, other):
-        return type(other) == ForceOnBoundary and other.stencil == self.stencil and \
+        return type(other) is ForceOnBoundary and other.stencil == self.stencil and \
             other.force_field == self.force_field

src/lbmpy/partially_saturated_cells.py

+import sympy as sp
+from dataclasses import dataclass
+
+from lbmpy.methods.abstractlbmethod import LbmCollisionRule
+from pystencils import Assignment, AssignmentCollection
+from pystencils.field import Field
+
+
+@dataclass
+class PSMConfig:
+    fraction_field: Field = None
+    """
+    Fraction field for PSM 
+    """
+
+    object_velocity_field: Field = None
+    """
+    Object velocity field for PSM 
+    """
+
+    SC: int = 1
+    """
+    Solid collision option for PSM
+    """
+
+    MaxParticlesPerCell: int = 1
+    """
+    Maximum number of particles overlapping with a cell 
+    """
+
+    individual_fraction_field: Field = None
+    """
+    Fraction field for each overlapping particle in PSM 
+    """
+
+    particle_force_field: Field = None
+    """
+    Force field for each overlapping particle in PSM 
+    """
+
+
+def add_psm_to_collision_rule(collision_rule, psm_config):
+    if psm_config.individual_fraction_field is None:
+        psm_config.individual_fraction_field = psm_config.fraction_field
+
+    method = collision_rule.method
+    pre_collision_pdf_symbols = method.pre_collision_pdf_symbols
+    stencil = method.stencil
+
+    # Get equilibrium from object velocity for solid collision
+    forces_rhs = [0] * psm_config.MaxParticlesPerCell * stencil.D
+    solid_collisions = [0] * stencil.Q
+    for p in range(psm_config.MaxParticlesPerCell):
+        equilibrium_fluid = method.get_equilibrium_terms()
+        equilibrium_solid = []
+        for eq in equilibrium_fluid:
+            eq_sol = eq
+            for i in range(stencil.D):
+                eq_sol = eq_sol.subs(sp.Symbol("u_" + str(i)),
+                                     psm_config.object_velocity_field.center(p * stencil.D + i), )
+            equilibrium_solid.append(eq_sol)
+
+        # Build solid collision
+        for i, (eqFluid, eqSolid, f, offset) in enumerate(
+                zip(equilibrium_fluid, equilibrium_solid, pre_collision_pdf_symbols, stencil)):
+            inverse_direction_index = stencil.stencil_entries.index(stencil.inverse_stencil_entries[i])
+            if psm_config.SC == 1:
+                solid_collision = psm_config.individual_fraction_field.center(p) * (
+                    (
+                        pre_collision_pdf_symbols[inverse_direction_index]
+                        - equilibrium_fluid[inverse_direction_index]
+                    )
+                    - (f - eqSolid)
+                )
+            elif psm_config.SC == 2:
+                # TODO get relaxation rate vector from method and use the right relaxation rate [i]
+                solid_collision = psm_config.individual_fraction_field.center(p) * (
+                    (eqSolid - f) + (1.0 - method.relaxation_rates[0]) * (f - eqFluid)
+                )
+            elif psm_config.SC == 3:
+                solid_collision = psm_config.individual_fraction_field.center(p) * (
+                    (
+                        pre_collision_pdf_symbols[inverse_direction_index]
+                        - equilibrium_solid[inverse_direction_index]
+                    )
+                    - (f - eqSolid)
+                )
+            else:
+                raise ValueError("Only SC=1, SC=2 and SC=3 are supported.")
+            solid_collisions[i] += solid_collision
+            for j in range(stencil.D):
+                forces_rhs[p * stencil.D + j] -= solid_collision * int(offset[j])
+
+    # Add solid collision to main assignments of collision rule
+    collision_assignments = []
+    for main, sc in zip(collision_rule.main_assignments, solid_collisions):
+        collision_assignments.append(Assignment(main.lhs, main.rhs + sc))
+
+    # Add hydrodynamic force calculations to collision assignments if two-way coupling is used
+    # (i.e., force field is not None)
+    if psm_config.particle_force_field is not None:
+        for p in range(psm_config.MaxParticlesPerCell):
+            for i in range(stencil.D):
+                collision_assignments.append(
+                    Assignment(
+                        psm_config.particle_force_field.center(p * stencil.D + i),
+                        forces_rhs[p * stencil.D + i],
+                    )
+                )
+
+    collision_assignments = AssignmentCollection(collision_assignments)
+    ac = LbmCollisionRule(method, collision_assignments, collision_rule.subexpressions,
+                          collision_rule.simplification_hints)
+    ac.topological_sort()
+    return ac

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)