lbmpy_tests/test_boundary_handling.py → tests/test_boundary_handling.py

 import numpy as np
 import pytest
 
-from lbmpy.boundaries import NoSlip, UBB, SimpleExtrapolationOutflow, ExtrapolationOutflow, \
-    FixedDensity, DiffusionDirichlet, NeumannByCopy, StreamInConstant, FreeSlip
-from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
+from lbmpy.boundaries import (NoSlip, NoSlipLinearBouzidi, QuadraticBounceBack,
+                              UBB, SimpleExtrapolationOutflow, ExtrapolationOutflow, FixedDensity, DiffusionDirichlet,
+                              NeumannByCopy, StreamInConstant, FreeSlip)
+from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling, create_lattice_boltzmann_boundary_kernel
 from lbmpy.creationfunctions import create_lb_function, create_lb_method, LBMConfig
 from lbmpy.enums import Stencil, Method
 from lbmpy.geometry import add_box_boundary
 from lbmpy.lbstep import LatticeBoltzmannStep
 from lbmpy.stencils import LBStencil
+import pystencils as ps
 from pystencils import create_data_handling, make_slice, Target, CreateKernelConfig
 from pystencils.slicing import slice_from_direction
 from pystencils.stencil import inverse_direction
@@ -26,7 +28,7 @@ def mirror_stencil(direction, mirror_axis):
 def test_simple(target):
     if target == Target.GPU:
         import pytest
-        pytest.importorskip('pycuda')
+        pytest.importorskip('cupy')
 
     dh = create_data_handling((4, 4), parallel=False, default_target=target)
     dh.add_array('pdfs', values_per_cell=9, cpu=True, gpu=target != Target.CPU)
@@ -45,7 +47,7 @@ def test_simple(target):
     bh = LatticeBoltzmannBoundaryHandling(lb_func.method, dh, 'pdfs', target=target)
 
     wall = NoSlip()
-    moving_wall = UBB((1, 0))
+    moving_wall = UBB((1, 0), density=1.0)
     bh.set_boundary(wall, make_slice[0, :])
     bh.set_boundary(wall, make_slice[-1, :])
     bh.set_boundary(wall, make_slice[:, 0])
@@ -435,3 +437,45 @@ def test_boundary_utility_functions():
     assert stream == StreamInConstant(constant=1.0, name="stream")
     assert not stream == StreamInConstant(constant=1.0, name="test")
     assert not stream == noslip
+
+
+@pytest.mark.parametrize("given_force_vector", [True, False])
+@pytest.mark.parametrize("dtype", ["float32", "float64"])
+def test_force_on_boundary(given_force_vector, dtype):
+    stencil = LBStencil(Stencil.D2Q9)
+    pdfs = ps.fields(f"pdfs_src({stencil.Q}): {dtype}[{stencil.D}D]", layout='fzyx')
+
+    method = create_lb_method(lbm_config=LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.8))
+
+    noslip = NoSlip(name="noslip", calculate_force_on_boundary=True)
+    bouzidi = NoSlipLinearBouzidi(name="bouzidi", calculate_force_on_boundary=True)
+    qq_bounce_Back = QuadraticBounceBack(name="qqBB", relaxation_rate=1.8, calculate_force_on_boundary=True)
+
+    boundary_objects = [noslip, bouzidi, qq_bounce_Back]
+    for boundary in boundary_objects:
+        if given_force_vector:
+            force_vector_type = np.dtype([(f"F_{i}", dtype) for i in range(stencil.D)], align=True)
+            force_vector = ps.Field('forceVector', ps.FieldType.INDEXED, force_vector_type, layout=[0],
+                                    shape=(ps.TypedSymbol("forceVectorSize", "int32"), 1), strides=(1, 1))
+        else:
+            force_vector = None
+
+        index_struct_dtype = _numpy_data_type_for_boundary_object(boundary, stencil.D)
+
+        index_field = ps.Field('indexVector', ps.FieldType.INDEXED, index_struct_dtype, layout=[0],
+                               shape=(ps.TypedSymbol("indexVectorSize", "int32"), 1), strides=(1, 1))
+
+        create_lattice_boltzmann_boundary_kernel(pdfs, index_field, method, boundary, force_vector=force_vector)
+
+
+def _numpy_data_type_for_boundary_object(boundary_object, dim):
+    boundary_index_array_coordinate_names = ["x", "y", "z"]
+    direction_member_name = "dir"
+    default_index_array_dtype = np.int32
+
+    coordinate_names = boundary_index_array_coordinate_names[:dim]
+    return np.dtype(
+        [(name, default_index_array_dtype) for name in coordinate_names]
+        + [(direction_member_name, default_index_array_dtype)]
+        + [(i[0], i[1].numpy_dtype) for i in boundary_object.additional_data],
+        align=True,)

tests/test_compiled_in_boundaries.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import pystencils as ps
+from lbmpy.session import *
+from lbmpy.boundaries.boundaries_in_kernel import update_rule_with_push_boundaries
+from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
+from collections import OrderedDict
+from time import perf_counter
+```
+
+%% Cell type:markdown id: tags:
+
+# Version 1: compile-in boundaries
+
+%% Cell type:code id: tags:
+
+``` python
+domain_size = (32, 32, 32)
+relaxation_rate = 1.8
+time_steps = 100
+lid_velocity = 0.05
+stencil = LBStencil(Stencil.D3Q19)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+dh = create_data_handling(domain_size, default_target=ps.Target.CPU)
+pdfs = dh.add_array('pdfs', values_per_cell=stencil.Q)
+u = dh.add_array('u', values_per_cell=stencil.D)
+streaming_pattern = 'aa'
+```
+
+%% Cell type:code id: tags:
+
+``` python
+boundaries = OrderedDict((
+    ((0, 1, 0), UBB([lid_velocity, 0, 0])),
+    ((1, 0, 0), NoSlip()),
+    ((-1, 0, 0), NoSlip()),
+    ((0, -1, 0), NoSlip()),
+    ((0, 0, 1), NoSlip()),
+    ((0, 0, -1), NoSlip()),
+))
+
+lbm_opt = LBMOptimisation(symbolic_field=pdfs, cse_global=False, cse_pdfs=False)
+lbm_config = LBMConfig(stencil=stencil, relaxation_rate=relaxation_rate, compressible=False)
+
+cr_even = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+cr_odd = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+update_rule_aa_even = update_rule_with_push_boundaries(cr_even, pdfs, boundaries, streaming_pattern, Timestep.EVEN)
+update_rule_aa_odd = update_rule_with_push_boundaries(cr_odd, pdfs, boundaries, streaming_pattern, Timestep.ODD)
+
+getter_assignments = macroscopic_values_getter(update_rule_aa_even.method, velocity=u.center_vector,
+                                               pdfs=pdfs, density=None,
+                                               streaming_pattern=streaming_pattern,
+                                               previous_timestep=Timestep.EVEN)
+
+config = ps.CreateKernelConfig(target=dh.default_target)
+getter_kernel = ps.create_kernel(getter_assignments, config=config).compile()
+even_kernel = ps.create_kernel(update_rule_aa_even, config=config).compile()
+odd_kernel = ps.create_kernel(update_rule_aa_odd, config=config).compile()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+def init():
+    dh.fill(pdfs.name, 0, ghost_layers=True)
+
+def aa_time_loop(steps=100):
+    assert steps % 2 == 0, "Works only for an even number of time steps"
+    dh.all_to_gpu()
+    for i in range(steps // 2):
+        dh.run_kernel(odd_kernel)
+        dh.run_kernel(even_kernel)
+    dh.run_kernel(getter_kernel)
+    dh.all_to_cpu()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+init()
+aa_time_loop(time_steps)
+vel_version1 = dh.gather_array(u.name, ghost_layers=False).copy()
+plt.vector_field_magnitude(vel_version1[:, :, domain_size[2]//2, :])
+plt.colorbar()
+```
+
+%% Output
+
+    <matplotlib.colorbar.Colorbar at 0x11b578710>
+
+
+
+%% Cell type:markdown id: tags:
+
+# Version 2: Normal boundary handling
+
+%% Cell type:code id: tags:
+
+``` python
+ldc = create_lid_driven_cavity(domain_size, relaxation_rate=relaxation_rate, lid_velocity=lid_velocity)
+ldc.run(time_steps)
+vel_version2 = ldc.velocity[:, :, :, :]
+
+plt.vector_field_magnitude(vel_version2[:, :, domain_size[2]//2, :])
+plt.colorbar()
+```
+
+%% Output
+
+    <matplotlib.colorbar.Colorbar at 0x12e912a90>
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+np.testing.assert_almost_equal(vel_version1[1:-1, 1:-1, :], vel_version2[1:-1, 1:-1, :], decimal=2)
+```