Cleaned wall model implementation

lbmpy/boundaries/__init__.py

@@ -2,9 +2,10 @@ from lbmpy.boundaries.boundaryconditions import (
     UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow,
     ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant, FreeSlip)
 from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
-from lbmpy.boundaries.wall_models import WallFunctionBounce
+from lbmpy.boundaries.wall_treatment.wall_models import WallFunctionBounce
+from lbmpy.boundaries.wall_treatment.spaldings_law import spaldings_law
 
 __all__ = ['NoSlip', 'FreeSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
            'FixedDensity', 'DiffusionDirichlet', 'NeumannByCopy',
            'LatticeBoltzmannBoundaryHandling', 'StreamInConstant',
-           'WallFunctionBounce']
+           'WallFunctionBounce', 'spaldings_law']

lbmpy/boundaries/wall_treatment/spaldings_law.py

+import sympy as sp
+
+
+def spaldings_law(u_plus, y_plus, kappa=0.41, B=5.5):
+    """
+    Returns a symbolic expression for spaldings law
+
+    Args:
+        u_plus: velocity nondimensionalized by the friction velocity u_tau
+        y_plus: distances from the wall nondimensionalized by the friction velocity u_tau
+        kappa: free parameter
+        B: free parameter
+    """
+    k_times_u = kappa * u_plus
+    fraction_1 = (k_times_u ** 2) / 2
+    fraction_2 = (k_times_u ** 3) / 6
+    return u_plus + sp.exp(-kappa * B) * (sp.exp(k_times_u) - 1 - k_times_u - fraction_1 - fraction_2) - y_plus

lbmpy/boundaries/wall_models.py → lbmpy/boundaries/wall_treatment/wall_models.py

 import sympy as sp
-from pystencils import Assignment
+import numpy as np
+
+from pystencils import Assignment, TypedSymbol
 from pystencils.stencil import offset_to_direction_string
 
 from lbmpy.advanced_streaming.indexing import MirroredStencilDirections, NeighbourOffsetArrays
 from lbmpy.boundaries.boundaryconditions import LbBoundary
+from lbmpy.boundaries.wall_treatment.spaldings_law import spaldings_law
 from lbmpy.relaxationrates import lattice_viscosity_from_relaxation_rate
 
 
@@ -13,17 +16,33 @@ class WallFunctionBounce(LbBoundary):
 
     Args:
         stencil: LBM stencil which is used for the simulation
+        pdfs: symbolic representation of the particle distribution functions.
         normal_direction: optional normal direction. If the Free slip boundary is applied to a certain side in the
                           domain it is not necessary to calculate the normal direction since it can be stated for all
                           boundary cells. This reduces the memory space for the index array significantly.
+        omega: relaxation rate used in the simulation
+        kappa: free parameter used for spaldings law
+        B: free parameter used for spaldings law
+        dt: time discretisation. Usually one in LB units
+        dy: space discretisation. Usually one in LB units
+        y: distance from the wall
+        newton_steps: number of newton steps to evaluate spaldings law.
         name: optional name of the boundary.
     """
 
-    def __init__(self, stencil, normal_direction, omega, vel_debug, name=None):
+    def __init__(self, stencil, pdfs, normal_direction, omega, kappa=0.41, B=5.5,
+                 dt=1, dy=1, y=0.5, newton_steps=10, name=None):
         """Set an optional name here, to mark boundaries, for example for force evaluations"""
         self.stencil = stencil
+        self.pdfs = pdfs
         self.omega = omega
-        self.vel_debug = vel_debug
+
+        self.kappa = kappa
+        self.B = B
+        self.dt = dt
+        self.dy = dy
+        self.y = y
+        self.newton_steps = newton_steps
 
         if len(normal_direction) - normal_direction.count(0) != 1:
             raise ValueError("Only normal directions for straight walls are supported for example (0, 1, 0) for "
@@ -32,7 +51,7 @@ class WallFunctionBounce(LbBoundary):
         self.mirror_axis = normal_direction.index(*[dir for dir in normal_direction if dir != 0])
 
         self.normal_direction = normal_direction
-        self.dim = len(stencil[0])
+        self.dim = stencil.D
 
         if name is None:
             name = f"WFB : {offset_to_direction_string([-x for x in normal_direction])}"
@@ -43,71 +62,71 @@ class WallFunctionBounce(LbBoundary):
         return [MirroredStencilDirections(self.stencil, self.mirror_axis), NeighbourOffsetArrays(lb_method.stencil)]
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
-        normal_direction = self.normal_direction
+        # needed symbols for offsets and indices
         mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(self.mirror_axis)
         mirrored_direction = inv_dir[sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]]
 
-        dt = 1
-        dy = 1
-        factor = dt / (2 * dy)
+        name_base = "f_in_inv_offsets_"
+        offset_array_symbols = [TypedSymbol(name_base + d, np.int64) for d in ['x', 'y', 'z']]
+        mirrored_offset = sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]
+        offsets = tuple(sp.IndexedBase(s, shape=(1,))[mirrored_offset] for s in offset_array_symbols)
 
-        stencil = self.stencil
-        cqc = lb_method.conserved_quantity_computation
-        rho = cqc.zeroth_order_moment_symbol
-        u = cqc.first_order_moment_symbols
-
-        newton_steps = 10
-
-        u_t = sp.symbols(f"u_tau_:{newton_steps + 1}")
+        # needed symbols in the Assignments
+        u_t = sp.symbols(f"u_tau_:{self.newton_steps + 1}")
         u_m = sp.Symbol("u_m")
-        B = 5.5
-        k = 0.41
-        delta = sp.symbols(f"delta_:{newton_steps}")
+        delta = sp.symbols(f"delta_:{self.newton_steps}")
         tau_w = sp.Symbol("tau_w")
         tau_w_x, tau_w_z = sp.symbols("tau_w_x tau_w_z")
 
-        pdf_center_vector = sp.Matrix([0] * stencil.Q)
+        normal_direction = self.normal_direction
+
+        # get velocity and density
+        cqc = lb_method.conserved_quantity_computation
+        rho = cqc.zeroth_order_moment_symbol
+        u = cqc.first_order_moment_symbols
+
+        pdf_center_vector = sp.Matrix([0] * self.stencil.Q)
 
-        for i in range(stencil.Q):
-            m = sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[i]
-            pdf_center_vector[i] = f_in[normal_direction](i)
+        for i in range(self.stencil.Q):
+            pdf_center_vector[i] = self.pdfs[offsets[0] + normal_direction[0],
+                                             offsets[1] + normal_direction[1],
+                                             offsets[2] + normal_direction[2]](i)
 
         eq_equations = cqc.equilibrium_input_equations_from_pdfs(pdf_center_vector)
+        result = eq_equations.all_assignments
 
         u_mag = Assignment(u_m, sp.sqrt(sum([x**2 for x in u])))
+        result.append(u_mag)
+
+        # using spaldings law
+
         nu = lattice_viscosity_from_relaxation_rate(self.omega)
         up = u_m / u_t[0]
 
-        y = 0.5
-        y_plus = (y * u_t[0]) / nu
-
-        spaldings_law = up + sp.exp(-k * B) * (
-                sp.exp(k * up) - 1 - (k * up) - ((k * up) ** 2) / 2 - ((k * up) ** 3) / 6) - y_plus
+        y_plus = (self.y * u_t[0]) / nu
 
-        m = -spaldings_law / spaldings_law.diff(u_t[0])
+        wall_law = spaldings_law(up, y_plus, kappa=self.kappa, B=self.B)
+        m = -wall_law / wall_law.diff(u_t[0])
 
         init_guess = Assignment(u_t[0], u_m)
 
-        newton_assignmets = []
-        for i in range(newton_steps):
-            newton_assignmets.append(Assignment(delta[i], m.subs({u_t[0]: u_t[i]})))
-            newton_assignmets.append(Assignment(u_t[i + 1], u_t[i] + delta[i]))
+        newton_assignments = []
+        for i in range(self.newton_steps):
+            newton_assignments.append(Assignment(delta[i], m.subs({u_t[0]: u_t[i]})))
+            newton_assignments.append(Assignment(u_t[i + 1], u_t[i] + delta[i]))
 
-        shear_stress = Assignment(tau_w, u_t[newton_steps]**2 * rho)
+        shear_stress = Assignment(tau_w, u_t[self.newton_steps]**2 * rho)
 
-        result = eq_equations.all_assignments
-        # result.append(Assignment(f_in(1), f_out(1)))
-        result.append(Assignment(self.vel_debug[0, 0, 0](0), u[0]))
-        result.append(Assignment(self.vel_debug[0, 0, 0](1), u[1]))
-        result.append(Assignment(self.vel_debug[0, 0, 0](2), u[2]))
-        result.append(u_mag)
         result.append(init_guess)
-        result.extend(newton_assignmets)
+        result.extend(newton_assignments)
         result.append(shear_stress)
 
+        # calculate tau_wx and tau_wz and use it to calculate the drag
+
         result.append(Assignment(tau_w_x, - u[0] / (sp.sqrt(u[0]**2 + u[2]**2)) * tau_w))
         result.append(Assignment(tau_w_z, - u[2] / (sp.sqrt(u[0] ** 2 + u[2] ** 2)) * tau_w))
 
+        factor = self.dt / (2 * self.dy)
         neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
         drag = neighbor_offset[0] * factor * tau_w_x + neighbor_offset[2] * factor * tau_w_z