Merge branch 'lees-edwards' into 'master'

Lees Edwards boundary conditions

Closes #7

See merge request pycodegen/lbmpy!49

.gitlab-ci.yml

@@ -13,6 +13,8 @@ tests-and-coverage:
       - $ENABLE_NIGHTLY_BUILDS
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
   script:
+    - env
+    - pip list
     - export NUM_CORES=$(nproc --all)
     - mkdir -p ~/.config/matplotlib
     - echo "backend:template" > ~/.config/matplotlib/matplotlibrc
@@ -38,6 +40,8 @@ test-longrun:
       - $ENABLE_NIGHTLY_BUILDS
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
   script:
+    - env 
+    - pip list
     - export NUM_CORES=$(nproc --all)
     - mkdir -p ~/.config/matplotlib
     - echo "backend:template" > ~/.config/matplotlib/matplotlibrc
@@ -103,6 +107,8 @@ pycodegen-integration:
   when: manual
   allow_failure: true
   script:
+    - env
+    - pip list
     - git clone https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pycodegen.git
     - cd pycodegen
     - git submodule sync --recursive

lbmpy_tests/test_lees_edwards.py

+from functools import partial
+
+import pystencils as ps
+from pystencils.astnodes import LoopOverCoordinate
+from pystencils.slicing import get_periodic_boundary_functor
+
+from lbmpy.creationfunctions import create_lb_update_rule
+from lbmpy.stencils import get_stencil
+from lbmpy.updatekernels import create_stream_pull_with_output_kernel
+from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
+from lbmpy.relaxationrates import lattice_viscosity_from_relaxation_rate
+
+import sympy as sp
+import numpy as np
+
+
+def get_le_boundary_functor(neighbor_stencil, shear_offset, ghost_layers=1, thickness=None, n=64):
+    functor_2 = get_periodic_boundary_functor(neighbor_stencil, ghost_layers, thickness)
+
+    def functor(pdfs, **_):
+
+        functor_2(pdfs)
+        weight = np.fmod(shear_offset[0] + n, 1.)
+
+        # First, we interpolate the upper pdfs
+        for i in range(len(pdfs[:, ghost_layers, :])):
+            ind1 = int(np.floor(i - shear_offset[0]) % n)
+            ind2 = int(np.ceil(i - shear_offset[0]) % n)
+
+            res = (1 - weight) * pdfs[ind1, ghost_layers, :] + weight * pdfs[ind2, ghost_layers, :]
+            pdfs[i, -ghost_layers, :] = res
+
+        # Second, we interpolate the lower pdfs
+        for i in range(len(pdfs[:, -ghost_layers, :])):
+            ind1 = int(np.floor(i + shear_offset[0]) % n)
+            ind2 = int(np.ceil(i + shear_offset[0]) % n)
+
+            res = (1 - weight) * pdfs[ind1, -ghost_layers - 1, :] + weight * pdfs[ind2, -ghost_layers - 1, :]
+            pdfs[i, ghost_layers - 1, :] = res
+
+    return functor
+
+
+def get_solution_navier_stokes(x, t, viscosity, velocity=1.0, height=1.0, max_iterations=10):
+    w = x / height - 0.5
+    for k in np.arange(1, max_iterations + 1):
+        w += 1.0 / (np.pi * k) * np.exp(-4 * np.pi ** 2 * viscosity * k ** 2 / height ** 2 * t) * np.sin(
+            2 * np.pi / height * k * x)
+    return velocity * w
+
+
+def test_lees_edwards():
+
+    domain_size = (64, 64)
+    omega = 1.0  # relaxation rate of first component
+    shear_velocity = 0.1  # shear velocity
+    shear_dir = 0  # direction of shear flow
+    shear_dir_normal = 1  # direction normal to shear plane, for interpolation
+
+    stencil = get_stencil("D2Q9")
+
+    dim = len(stencil[0])
+    dh = ps.create_data_handling(domain_size, periodicity=True, default_target='cpu')
+
+    src = dh.add_array('src', values_per_cell=len(stencil))
+    dh.fill('src', 1.0, ghost_layers=True)
+
+    dst = dh.add_array_like('dst', 'src')
+    dh.fill('dst', 0.0, ghost_layers=True)
+
+    force = dh.add_array('force', values_per_cell=dh.dim)
+    dh.fill('force', 0.0, ghost_layers=True)
+
+    rho = dh.add_array('rho', values_per_cell=1)
+    dh.fill('rho', 1.0, ghost_layers=True)
+    u = dh.add_array('u', values_per_cell=dh.dim)
+    dh.fill('u', 0.0, ghost_layers=True)
+
+    counters = [LoopOverCoordinate.get_loop_counter_symbol(i) for i in range(dim)]
+    points_up = sp.Symbol('points_up')
+    points_down = sp.Symbol('points_down')
+
+    u_p = sp.Piecewise((1, sp.And(ps.data_types.type_all_numbers(counters[1] <= 1, 'int'), points_down)),
+                       (-1, sp.And(ps.data_types.type_all_numbers(counters[1] >= src.shape[1] - 2, 'int'),
+                                   points_up)), (0, True)) * shear_velocity
+
+    collision = create_lb_update_rule(stencil=stencil,
+                                      relaxation_rate=omega,
+                                      compressible=True,
+                                      velocity_input=u.center_vector + sp.Matrix([u_p, 0]),
+                                      density_input=rho,
+                                      force_model='luo',
+                                      force=force.center_vector,
+                                      kernel_type='collide_only',
+                                      optimization={'symbolic_field': src})
+
+    to_insert = [s.lhs for s in collision.subexpressions
+                 if collision.method.first_order_equilibrium_moment_symbols[shear_dir]
+                 in s.free_symbols]
+    for s in to_insert:
+        collision = collision.new_with_inserted_subexpression(s)
+    ma = []
+    for a, c in zip(collision.main_assignments, collision.method.stencil):
+        if c[shear_dir_normal] == -1:
+            b = (True, False)
+        elif c[shear_dir_normal] == 1:
+            b = (False, True)
+        else:
+            b = (False, False)
+        a = ps.Assignment(a.lhs, a.rhs.replace(points_down, b[0]))
+        a = ps.Assignment(a.lhs, a.rhs.replace(points_up, b[1]))
+        ma.append(a)
+    collision.main_assignments = ma
+
+    stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
+                                                   {'density': rho, 'velocity': u})
+
+    stream_kernel = ps.create_kernel(stream, target=dh.default_target).compile()
+    collision_kernel = ps.create_kernel(collision, target=dh.default_target).compile()
+
+    init = macroscopic_values_setter(collision.method, velocity=(0, 0),
+                                     pdfs=src.center_vector, density=rho.center)
+    init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
+
+    offset = [0.0]
+
+    sync_pdfs = dh.synchronization_function([src.name],
+                                            functor=partial(get_le_boundary_functor, shear_offset=offset))
+
+    dh.run_kernel(init_kernel)
+
+    time = 500
+
+    dh.all_to_gpu()
+    for i in range(time):
+        dh.run_kernel(collision_kernel)
+
+        sync_pdfs()
+        dh.run_kernel(stream_kernel)
+
+        dh.swap(src.name, dst.name)
+        offset[0] += shear_velocity
+    dh.all_to_cpu()
+
+    nu = lattice_viscosity_from_relaxation_rate(omega)
+    h = domain_size[0]
+    k_max = 100
+
+    analytical_solution = get_solution_navier_stokes(np.linspace(0.5, h - 0.5, h), time, nu, shear_velocity, h, k_max)
+    np.testing.assert_array_almost_equal(analytical_solution, dh.gather_array(u.name)[0, :, 0], decimal=5)
+
+    dh.fill(rho.name, 1.0, ghost_layers=True)
+    dh.run_kernel(init_kernel)
+    dh.fill(u.name, 0.0, ghost_layers=True)
+    dh.fill('force', 0.0, ghost_layers=True)
+    dh.cpu_arrays[force.name][64 // 3, 1, :] = [1e-2, -1e-1]
+
+    offset[0] = 0
+
+    time = 20
+
+    dh.all_to_gpu()
+    for i in range(time):
+        dh.run_kernel(collision_kernel)
+
+        sync_pdfs()
+        dh.run_kernel(stream_kernel)
+
+        dh.swap(src.name, dst.name)
+    dh.all_to_cpu()
+
+    vel_unshifted = np.array(dh.gather_array(u.name)[:, -3:-1, :])
+
+    dh.fill(rho.name, 1.0, ghost_layers=True)
+    dh.run_kernel(init_kernel)
+    dh.fill(u.name, 0.0, ghost_layers=True)
+    dh.fill('force', 0.0, ghost_layers=True)
+    dh.cpu_arrays[force.name][64 // 3, 1, :] = [1e-2, -1e-1]
+
+    offset[0] = 10
+
+    time = 20
+
+    dh.all_to_gpu()
+    for i in range(time):
+        dh.run_kernel(collision_kernel)
+
+        sync_pdfs()
+        dh.run_kernel(stream_kernel)
+
+        dh.swap(src.name, dst.name)
+    dh.all_to_cpu()
+
+    vel_shifted = np.array(dh.gather_array(u.name)[:, -3:-1, :])
+
+    vel_rolled = np.roll(vel_shifted, -offset[0], axis=0)
+
+    np.testing.assert_array_almost_equal(vel_unshifted, vel_rolled)

lbmpy_tests/test_stokes_setup.ipynb

@@ -793,7 +793,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.8.2"
   }
  },
  "nbformat": 4,