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/analytical.py → src/lbmpy/phasefield/analytical.py

@@ -17,11 +17,11 @@ def symmetric_symbolic_surface_tension(i, j):
     if i == j:
         return 0
     index = (i, j) if i < j else (j, i)
-    return sp.Symbol("%s_%d_%d" % ((surface_tension_symbol_name,) + index))
+    return sp.Symbol(f"{surface_tension_symbol_name}_{index[0]}_{index[1]}")
 
 
 def symbolic_order_parameters(num_symbols):
-    return sp.symbols("%s_:%i" % (order_parameter_symbol_name, num_symbols))
+    return sp.symbols(f"{order_parameter_symbol_name}_:{num_symbols}")
 
 
 def free_energy_functional_3_phases(order_parameters=None, interface_width=interface_width_symbol, transformed=True,
@@ -65,7 +65,7 @@ def free_energy_functional_n_phases_penalty_term(order_parameters, interface_wid
                                                  penalty_term_factor=0.01):
     num_phases = len(order_parameters)
     if kappa is None:
-        kappa = sp.symbols("kappa_:%d" % (num_phases,))
+        kappa = sp.symbols(f"kappa_:{num_phases}")
     if not hasattr(kappa, "__len__"):
         kappa = [kappa] * num_phases
 
@@ -166,7 +166,8 @@ def free_energy_functional_n_phases(num_phases=None, surface_tensions=symmetric_
 
     def lambda_coeff(k, l):
         if symbolic_lambda:
-            return sp.Symbol("Lambda_%d%d" % ((k, l) if k < l else (l, k)))
+            symbol_names = (k, l) if k < l else (l, k)
+            return sp.Symbol(f"Lambda_{symbol_names[0]}{symbol_names[1]}")
         n = num_phases - 1
         if k == l:
             assert surface_tensions(l, l) == 0
@@ -241,7 +242,7 @@ def chemical_potentials_from_free_energy(free_energy, order_parameters=None):
 
 def force_from_phi_and_mu(order_parameters, dim, mu=None):
     if mu is None:
-        mu = sp.symbols("mu_:%d" % (len(order_parameters),))
+        mu = sp.symbols(f"mu_:{len(order_parameters)}")
 
     return sp.Matrix([sum(- c_i * Diff(mu_i, a) for c_i, mu_i in zip(order_parameters, mu))
                       for a in range(dim)])
@@ -298,7 +299,7 @@ def extract_gamma(free_energy, order_parameters):
             continue
 
         if len(diff_factors) != 2:
-            raise ValueError("Could not determine Λ because of term " + str(product))
+            raise ValueError(f"Could not determine Λ because of term {str(product)}")
 
         indices = sorted([order_parameters.index(d.args[0]) for d in diff_factors])
         increment = prod(e for e in product if e.func != Diff)

lbmpy/phasefield/cahn_hilliard_lbm.py → src/lbmpy/phasefield/cahn_hilliard_lbm.py

 import sympy as sp
 
 from lbmpy.maxwellian_equilibrium import get_weights
+from lbmpy.equilibrium import GenericDiscreteEquilibrium
+from lbmpy.methods import DensityVelocityComputation
 from lbmpy.methods.creationfunctions import create_from_equilibrium
-from lbmpy.stencils import get_stencil
 from pystencils.sympyextensions import kronecker_delta, multidimensional_sum
 
 
@@ -19,8 +20,6 @@ def cahn_hilliard_lb_method(stencil, mu, relaxation_rate=sp.Symbol("omega"), gam
         relaxation_rate: relaxation rate of method
         gamma: tunable parameter affecting the second order equilibrium moment
     """
-    if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
     weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
 
     kd = kronecker_delta
@@ -29,17 +28,21 @@ def cahn_hilliard_lb_method(stencil, mu, relaxation_rate=sp.Symbol("omega"), gam
         for r in multidimensional_sum(*args, dim=len(stencil[0])):
             yield r
 
-    op = sp.Symbol("rho")
-    v = sp.symbols("u_:%d" % (len(stencil[0]),))
+    compressible = True
+    zero_centered = False
+    cqc = DensityVelocityComputation(stencil, compressible, zero_centered)
+    rho = cqc.density_symbol
+    v = cqc.velocity_symbols
 
-    equilibrium = []
+    equilibrium_terms = []
     for d, w in zip(stencil, weights):
         c_s = sp.sqrt(sp.Rational(1, 3))
         result = gamma * mu / (c_s ** 2)
-        result += op * sum(d[i] * v[i] for i, in s(1)) / (c_s ** 2)
-        result += op * sum(v[i] * v[j] * (d[i] * d[j] - c_s ** 2 * kd(i, j)) for i, j in s(2)) / (2 * c_s ** 4)
-        equilibrium.append(w * result)
-
-    rho = sp.Symbol("rho")
-    equilibrium[0] = rho - sp.expand(sum(equilibrium[1:]))
-    return create_from_equilibrium(stencil, tuple(equilibrium), relaxation_rate, compressible=True)
+        result += rho * sum(d[i] * v[i] for i, in s(1)) / (c_s ** 2)
+        result += rho * sum(v[i] * v[j] * (d[i] * d[j] - c_s ** 2 * kd(i, j)) for i, j in s(2)) / (2 * c_s ** 4)
+        equilibrium_terms.append(w * result)
+
+    equilibrium_terms[0] = rho - sp.expand(sum(equilibrium_terms[1:]))
+    equilibrium = GenericDiscreteEquilibrium(stencil, equilibrium_terms, rho, v, deviation_only=False)
+    
+    return create_from_equilibrium(stencil, equilibrium, cqc, relaxation_rate, zero_centered=zero_centered)

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)
 
 

src/lbmpy/phasefield/fd_stencils.py

+import sympy as sp
+
+from pystencils.fd.spatial import fd_stencils_standard
+
+
+def fd_stencils_isotropic_high_density_code(indices, dx, fa):
+    dim = fa.field.spatial_dimensions
+    if dim == 1:
+        return fd_stencils_standard(indices, dx, fa)
+    elif dim == 2:
+        order = len(indices)
+
+        if order == 1:
+            idx = indices[0]
+            assert 0 <= idx < 2
+            other_idx = 1 if indices[0] == 0 else 0
+            weights = {-1: sp.Rational(1, 12) / dx,
+                       0: sp.Rational(1, 3) / dx,
+                       1: sp.Rational(1, 12) / dx}
+            upper_terms = sum(fa.neighbor(idx, +1).neighbor(other_idx, off) * w for off, w in weights.items())
+            lower_terms = sum(fa.neighbor(idx, -1).neighbor(other_idx, off) * w for off, w in weights.items())
+            return upper_terms - lower_terms
+        elif order == 2:
+            if indices[0] == indices[1]:
+                idx = indices[0]
+                diagonals = sp.Rational(1, 8) * sum(fa.neighbor(0, i).neighbor(1, j) for i in (-1, 1) for j in (-1, 1))
+                div_direction = sp.Rational(1, 2) * sum(fa.neighbor(idx, i) for i in (-1, 1))
+                center = - sp.Rational(3, 2) * fa
+                return (diagonals + div_direction + center) / (dx ** 2)
+            else:
+                return fd_stencils_standard(indices, dx, fa)
+    raise NotImplementedError("Supports only derivatives up to order 2 for 1D and 2D setups")

lbmpy/phasefield/kerneleqs.py → src/lbmpy/phasefield/kerneleqs.py

@@ -4,7 +4,7 @@ from lbmpy.phasefield.analytical import (
     chemical_potentials_from_free_energy, force_from_phi_and_mu, force_from_pressure_tensor,
     pressure_tensor_bulk_sqrt_term, pressure_tensor_from_free_energy, substitute_laplacian_by_sum,
     symmetric_tensor_linearization)
-from pystencils import Assignment
+from pystencils import Assignment, Target
 from pystencils.fd import Discretization2ndOrder, discretize_spatial
 
 # ---------------------------------- Kernels to compute force ----------------------------------------------------------
@@ -90,8 +90,8 @@ def cahn_hilliard_fd_eq(phase_idx, phi, mu, velocity, mobility, dx, dt):
 
 
 class CahnHilliardFDStep:
-    def __init__(self, data_handling, phi_field_name, mu_field_name, velocity_field_name, name='ch_fd', target='cpu',
-                 dx=1, dt=1, mobilities=1, equation_modifier=lambda eqs: eqs):
+    def __init__(self, data_handling, phi_field_name, mu_field_name, velocity_field_name, name='ch_fd',
+                 target=Target.CPU, dx=1, dt=1, mobilities=1, equation_modifier=lambda eqs: eqs):
         from pystencils import create_kernel
         self.data_handling = data_handling
 

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

-import pyximport
+try:
+    import pyximport
+
+    pyximport.install(language_level=3)
+    from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
+except ImportError:
+    raise ImportError("pyximport not available. Please install Cython to use simplex_projection_2d.")
+
 import sympy as sp
 
-from lbmpy.creationfunctions import create_lb_update_rule
+from lbmpy.creationfunctions import create_lb_update_rule, LBMConfig, LBMOptimisation
+from lbmpy.enums import Stencil
 from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
 from lbmpy.phasefield.analytical import chemical_potentials_from_free_energy, force_from_phi_and_mu
 from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
-from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
-from lbmpy.stencils import get_stencil
+from lbmpy.stencils import LBStencil
 from pystencils import Assignment, create_data_handling, create_kernel
 from pystencils.fd import Diff, discretize_spatial, expand_diff_full
 from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
 
-pyximport.install(language_level=3)
-
 
 def forth_order_isotropic_discretize(field):
     second_neighbor_stencil = [(i, j)
@@ -71,7 +76,7 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
     dh = create_data_handling(domain_size, periodicity=(True, True), default_ghost_layers=2)
 
     c = dh.add_array("c", values_per_cell=num_phases)
-    rho = dh.add_array("rho")
+    rho = dh.add_array("rho", values_per_cell=1)
     mu = dh.add_array("mu", values_per_cell=num_phases, latex_name="\\mu")
     force = dh.add_array("F", values_per_cell=dh.dim)
     u = dh.add_array("u", values_per_cell=dh.dim)
@@ -80,8 +85,8 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
     pdf_field = []
     pdf_dst_field = []
     for i in range(num_phases):
-        pdf_field_local = dh.add_array("pdf_ch_%d" % i, values_per_cell=9)  # 9 for D2Q9
-        pdf_dst_field_local = dh.add_array("pdfs_ch_%d_dst" % i, values_per_cell=9)
+        pdf_field_local = dh.add_array(f"pdf_ch_{i}", values_per_cell=9)  # 9 for D2Q9
+        pdf_dst_field_local = dh.add_array(f"pdfs_ch_{i}_dst", values_per_cell=9)
         pdf_field.append(pdf_field_local)
         pdf_dst_field.append(pdf_dst_field_local)
 
@@ -110,15 +115,16 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
     ch_collide_kernels = []
     ch_methods = []
     for i in range(num_phases):
-        ch_method = cahn_hilliard_lb_method(get_stencil("D2Q9"), mu(i),
+        ch_method = cahn_hilliard_lb_method(LBStencil(Stencil.D2Q9), mu(i),
                                             relaxation_rate=1.0, gamma=1.0)
         ch_methods.append(ch_method)
-        ch_update_rule = create_lb_update_rule(lb_method=ch_method,
-                                               kernel_type='collide_only',
-                                               density_input=c(i),
-                                               velocity_input=u.center_vector,
-                                               compressible=True,
-                                               optimization={"symbolic_field": pdf_field[i]})
+
+        lbm_config = LBMConfig(lb_method=ch_method, kernel_type='collide_only', density_input=c(i),
+                               velocity_input=u.center_vector, compressible=True, zero_centered=False)
+        lbm_opt = LBMOptimisation(symbolic_field=pdf_field[i])
+        ch_update_rule = create_lb_update_rule(lbm_config=lbm_config,
+                                               lbm_optimisation=lbm_opt)
+
         ch_assign = ch_update_rule.all_assignments
         ch_kernel = create_kernel(ch_assign).compile()
         ch_collide_kernels.append(ch_kernel)
@@ -127,10 +133,11 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
     for i in range(num_phases):
         ch_method = ch_methods[i]
 
-        ch_update_rule = create_lb_update_rule(lb_method=ch_method,
-                                               kernel_type='stream_pull_only',
-                                               temporary_field_name=pdf_dst_field[i].name,
-                                               optimization={"symbolic_field": pdf_field[i]})
+        lbm_config = LBMConfig(lb_method=ch_method, kernel_type='stream_pull_only',
+                               temporary_field_name=pdf_dst_field[i].name, zero_centered=False)
+        lbm_opt = LBMOptimisation(symbolic_field=pdf_field[i])
+        ch_update_rule = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+
         ch_assign = ch_update_rule.all_assignments
         ch_kernel = create_kernel(ch_assign).compile()
         ch_stream_kernels.append(ch_kernel)
@@ -140,7 +147,9 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
     for i in range(num_phases):
         ch_method = ch_methods[i]
         init_assign = pdf_initialization_assignments(lb_method=ch_method,
-                                                     density=c_vec[i], velocity=(0, 0), pdfs=pdf_field[i].center_vector)
+                                                     density=c_vec[i],
+                                                     velocity=(0, 0),
+                                                     pdfs=pdf_field[i].center_vector)
         init_kernel = create_kernel(init_assign).compile()
         init_kernels.append(init_kernel)
 
@@ -152,17 +161,17 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
         getter_kernel = create_kernel(output_assign).compile()
         getter_kernels.append(getter_kernel)
 
-    collide_assign = create_lb_update_rule(kernel_type='collide_only',
-                                           relaxation_rate=1.0,
-                                           force=force,
-                                           optimization={"symbolic_field": pdf_hydro_field},
-                                           compressible=True)
+    lbm_config = LBMConfig(kernel_type='collide_only', relaxation_rate=1.0, force=force,
+                           compressible=True, zero_centered=False)
+    lbm_opt = LBMOptimisation(symbolic_field=pdf_hydro_field)
+    collide_assign = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
     collide_kernel = create_kernel(collide_assign).compile()
 
-    stream_assign = create_lb_update_rule(kernel_type='stream_pull_only',
-                                          temporary_field_name=pdf_hydro_dst_field.name,
-                                          optimization={"symbolic_field": pdf_hydro_field},
-                                          output={"density": rho, "velocity": u})
+    lbm_config = LBMConfig(kernel_type='stream_pull_only', temporary_field_name=pdf_hydro_dst_field.name,
+                           zero_centered=False, output={"density": rho, "velocity": u})
+    lbm_opt = LBMOptimisation(symbolic_field=pdf_hydro_field)
+    stream_assign = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+
     stream_kernel = create_kernel(stream_assign).compile()
 
     method_collide = collide_assign.method

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)