Merge branch 'FixIBC' into 'master'

[Fix] Interpolation Bounce back

See merge request pycodegen/lbmpy!160

src/lbmpy/boundaries/boundaryconditions.py

@@ -244,12 +244,12 @@ class QuadraticBounceBack(LbBoundary):
         return [LbmWeightInfo(lb_method, self.data_type), inverse_dir_node, NeighbourOffsetArrays(lb_method.stencil)]
 
     @staticmethod
-    def get_equilibrium(v, u, rho, drho, weight, compressible, deviation_only):
+    def get_equilibrium(v, u, rho, drho, weight, compressible, zero_centered):
         rho_background = sp.Integer(1)
 
         result = discrete_equilibrium(v, u, rho, weight,
                                       order=2, c_s_sq=sp.Rational(1, 3), compressible=compressible)
-        if deviation_only:
+        if zero_centered:
             shift = discrete_equilibrium(v, [0] * len(u), rho_background, weight,
                                          order=0, c_s_sq=sp.Rational(1, 3), compressible=False)
             result = simplify_by_equality(result - shift, rho, drho, rho_background)
@@ -257,6 +257,7 @@ class QuadraticBounceBack(LbBoundary):
         return result
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        omega = self.relaxation_rate
         inv = sp.IndexedBase(self.inv_dir_symbol, shape=(1,))[dir_symbol]
         weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
         weight_of_direction = weight_info.weight_of_direction
@@ -272,7 +273,7 @@ class QuadraticBounceBack(LbBoundary):
         v = [TypedSymbol(f"c_{i}", self.data_type) for i in range(lb_method.stencil.D)]
         v_inv = [TypedSymbol(f"c_inv_{i}", self.data_type) for i in range(lb_method.stencil.D)]
         one = sp.Float(1.0)
-        two = sp.Float(2.0)
+        half = sp.Rational(1, 2)
 
         subexpressions = [Assignment(pdf_symbols[i], pdf) for i, pdf in enumerate(pdf_field_accesses)]
         subexpressions.append(Assignment(f_xf, f_out(dir_symbol)))
@@ -294,24 +295,19 @@ class QuadraticBounceBack(LbBoundary):
         drho = cqc.density_deviation_symbol
         u = sp.Matrix(cqc.velocity_symbols)
         compressible = cqc.compressible
-        deviation_only = lb_method.equilibrium_distribution.deviation_only
+        zero_centered = cqc.zero_centered_pdfs
 
         cqe = cqc.equilibrium_input_equations_from_pdfs(pdf_symbols, False)
         subexpressions.append(cqe.all_assignments)
 
-        eq_dir = self.get_equilibrium(v, u, rho, drho, weight, compressible, deviation_only)
-        eq_inv = self.get_equilibrium(v_inv, u, rho, drho, weight_inv, compressible, deviation_only)
+        eq_dir = self.get_equilibrium(v, u, rho, drho, weight, compressible, zero_centered)
+        eq_inv = self.get_equilibrium(v_inv, u, rho, drho, weight_inv, compressible, zero_centered)
 
         subexpressions.append(Assignment(feq, eq_dir + eq_inv))
 
-        # equation E.4 first summand
-        e41 = (f_xf_inv - f_xf) / two
-        # equation E.4 second summand
-        e42 = (f_xf_inv + f_xf - self.relaxation_rate * feq) / (two - two * self.relaxation_rate)
-        # equation E.3
-        fw = ((one - q) * (e41 + e42)) + q * f_xf_inv
-        # equation E.1
-        result = (one / (q + one)) * fw + (q / (q + one)) * f_xf
+        t1 = (f_xf - f_xf_inv + (f_xf + f_xf_inv - feq * omega) / (one - omega))
+        t2 = (q * (f_xf + f_xf_inv)) / (one + q)
+        result = (one - q) / (one + q) * t1 * half + t2
 
         boundary_assignments = [Assignment(f_in(inv_dir[dir_symbol]), result)]
         return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)

src/lbmpy/maxwellian_equilibrium.py

@@ -85,7 +85,7 @@ def discrete_maxwellian_equilibrium(stencil, rho=sp.Symbol("rho"), u=sp.symbols(
 def discrete_equilibrium(v=sp.symbols("v_:3"), u=sp.symbols("u_:3"), rho=sp.Symbol("rho"), weight=sp.Symbol("w"),
                          order=2, c_s_sq=sp.Symbol("c_s") ** 2, compressible=True):
     """
-    Returns the common discrete LBM equilibrium as a list of sympy expressions
+    Returns the common discrete LBM equilibrium depending on the mesoscopic velocity and the directional lattice weight
 
     Args:
         v: symbols for mesoscopic velocity

src/lbmpy/oldroydb.py

@@ -177,9 +177,10 @@ 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)
 

tests/test_cumulant_transform.py

@@ -6,11 +6,12 @@ from pystencils import AssignmentCollection
 
 from lbmpy.moment_transforms import (
     CentralMomentsToCumulantsByGeneratingFunc, PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT,
-    PRE_COLLISION_CUMULANT, PRE_COLLISION_MONOMIAL_CUMULANT
+    PRE_COLLISION_CUMULANT
 )
 from lbmpy.methods import cascaded_moment_sets_literature
 from lbmpy.stencils import Stencil, LBStencil
 
+
 @pytest.mark.parametrize("stencil", [Stencil.D2Q9, Stencil.D3Q19])
 def test_identity(stencil):
     stencil = LBStencil(stencil)

tests/test_interpolation_boundaries.py

+import pytest
+import numpy as np
+from pystencils.slicing import slice_from_direction
+
+from lbmpy import LBMConfig, LBStencil, Stencil, Method
+from lbmpy.boundaries.boundaryconditions import NoSlip, NoSlipLinearBouzidi, QuadraticBounceBack, UBB
+from lbmpy.lbstep import LatticeBoltzmannStep
+
+
+def check_velocity(noslip_velocity, interpolated_velocity, wall_distance):
+    # First we check if the flow is fully developed
+    np.testing.assert_almost_equal(np.gradient(np.gradient(noslip_velocity)), 0, decimal=8)
+    np.testing.assert_almost_equal(np.gradient(np.gradient(interpolated_velocity)), 0, decimal=8)
+
+    # If the wall is closer to the first fluid cell we expect a lower velocity at the first fluid cell
+    if wall_distance < 0.5:
+        assert noslip_velocity[0] > interpolated_velocity[0]
+    # If the wall is further away from the first fluid cell we expect a higher velocity at the first fluid cell
+    if wall_distance > 0.5:
+        assert noslip_velocity[0] < interpolated_velocity[0]
+    # If the wall cuts the cell halfway the interpolation BC should approximately fall back to a noslip boundary
+    if wall_distance == 0.5:
+        np.testing.assert_almost_equal(noslip_velocity[0], interpolated_velocity[0], decimal=2)
+
+
+def couette_flow(stencil, method_enum, zero_centered, wall_distance, compressible):
+    dim = stencil.D
+    if dim == 2:
+        domain_size = (50, 25)
+        wall_velocity = (0.01, 0)
+        periodicity = (True, False)
+    else:
+        domain_size = (50, 25, 25)
+        wall_velocity = (0.01, 0, 0)
+        periodicity = (True, False, True)
+
+    timesteps = 10000
+    omega = 1.1
+
+    lbm_config = LBMConfig(stencil=stencil,
+                           method=method_enum,
+                           relaxation_rate=omega,
+                           zero_centered=zero_centered,
+                           compressible=compressible,
+                           weighted=True)
+
+    lb_step_noslip = LatticeBoltzmannStep(domain_size=domain_size, periodicity=periodicity,
+                                          lbm_config=lbm_config, compute_velocity_in_every_step=True)
+    lb_step_bouzidi = LatticeBoltzmannStep(domain_size=domain_size, periodicity=periodicity,
+                                           lbm_config=lbm_config, compute_velocity_in_every_step=True)
+
+    lb_step_quadratic_bb = LatticeBoltzmannStep(domain_size=domain_size, periodicity=periodicity,
+                                                lbm_config=lbm_config, compute_velocity_in_every_step=True)
+
+    def init_wall_distance(boundary_data, **_):
+        for cell in boundary_data.index_array:
+            cell['q'] = wall_distance
+
+    moving_wall = UBB(wall_velocity)
+    noslip = NoSlip("wall")
+    bouzidi = NoSlipLinearBouzidi("wall", init_wall_distance=init_wall_distance)
+    quadratic_bb = QuadraticBounceBack(omega, "wall", init_wall_distance=init_wall_distance)
+
+    lb_step_noslip.boundary_handling.set_boundary(noslip, slice_from_direction('S', dim))
+    lb_step_noslip.boundary_handling.set_boundary(moving_wall, slice_from_direction('N', dim))
+
+    lb_step_bouzidi.boundary_handling.set_boundary(bouzidi, slice_from_direction('S', dim))
+    lb_step_bouzidi.boundary_handling.set_boundary(moving_wall, slice_from_direction('N', dim))
+
+    lb_step_quadratic_bb.boundary_handling.set_boundary(quadratic_bb, slice_from_direction('S', dim))
+    lb_step_quadratic_bb.boundary_handling.set_boundary(moving_wall, slice_from_direction('N', dim))
+
+    lb_step_noslip.run(timesteps)
+    lb_step_bouzidi.run(timesteps)
+    lb_step_quadratic_bb.run(timesteps)
+
+    if dim == 2:
+        noslip_velocity = lb_step_noslip.velocity[domain_size[0] // 2, :, 0]
+        bouzidi_velocity = lb_step_bouzidi.velocity[domain_size[0] // 2, :, 0]
+        quadratic_bb_velocity = lb_step_quadratic_bb.velocity[domain_size[0] // 2, :, 0]
+    else:
+        noslip_velocity = lb_step_noslip.velocity[domain_size[0] // 2, :, domain_size[2] // 2, 0]
+        bouzidi_velocity = lb_step_bouzidi.velocity[domain_size[0] // 2, :, domain_size[2] // 2, 0]
+        quadratic_bb_velocity = lb_step_quadratic_bb.velocity[domain_size[0] // 2, :, domain_size[2] // 2, 0]
+
+    check_velocity(noslip_velocity, bouzidi_velocity, wall_distance)
+    check_velocity(noslip_velocity, quadratic_bb_velocity, wall_distance)
+
+
+@pytest.mark.parametrize("zero_centered", [False, True])
+@pytest.mark.parametrize("wall_distance", [0.1, 0.5, 0.9])
+@pytest.mark.parametrize("compressible", [True, False])
+def test_couette_flow_short(zero_centered, wall_distance, compressible):
+    stencil = LBStencil(Stencil.D2Q9)
+    couette_flow(stencil, Method.SRT, zero_centered, wall_distance, compressible)
+
+
+@pytest.mark.parametrize("method_enum", [Method.MRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
+@pytest.mark.parametrize("zero_centered", [False, True])
+@pytest.mark.parametrize("wall_distance", [0.1, 0.5, 0.9])
+@pytest.mark.parametrize("compressible", [True, False])
+@pytest.mark.longrun
+def test_couette_flow_long(method_enum, zero_centered, wall_distance, compressible):
+    if method_enum is Method.CUMULANT and compressible is False:
+        pytest.skip("incompressible cumulant is not supported")
+
+    stencil = LBStencil(Stencil.D2Q9)
+    couette_flow(stencil, method_enum, zero_centered, wall_distance, compressible)
+
+
+@pytest.mark.parametrize("method_enum", [Method.SRT, Method.MRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
+@pytest.mark.parametrize("wall_distance", [0.1, 0.5, 0.9])
+@pytest.mark.parametrize("stencil", [Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.longrun
+def test_couette_flow_d3d(method_enum, wall_distance, stencil):
+    stencil = LBStencil(stencil)
+    couette_flow(stencil, method_enum, True, wall_distance, True)
+
+