shear flow test added to test the second moment

lbmpy_tests/test_shear_flow.py

+"""
+Test Poiseuille flow against analytical solution
+"""
+
+
+import numpy as np
+import pytest
+import sympy as sp
+
+import lbmpy
+from lbmpy.boundaries import UBB
+from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
+from lbmpy.creationfunctions import create_lb_update_rule, create_stream_pull_with_output_kernel
+from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
+from lbmpy.stencils import get_stencil
+
+import pystencils as ps
+
+
+def shear_flow(x, t, nu, v, h, k_max):
+    """
+    Analytical solution for driven shear flow between two plates.
+
+    Parameters
+    ----------
+    x : :obj:`float`
+        Position from the left plane.
+    t : :obj:`float`
+        Time since start of the shearing.
+    nu : :obj:`float`
+        Kinematic viscosity.
+    v : :obj:`float`
+        Shear rate.
+    h : :obj:`float`
+        Distance between the plates.
+    k_max : :obj:`int`
+        Maximum considered wave number.
+
+    Returns
+    -------
+    :obj:`double` : Analytical velocity
+
+    """
+
+    u = x / h - 0.5
+    for k in np.arange(1, k_max + 1):
+        u += 1.0 / (np.pi * k) * np.exp(
+            -4 * np.pi ** 2 * nu * k ** 2 / h ** 2 * t) * np.sin(2 * np.pi / h * k * x)
+    return -v * u
+
+
+rho_0 = 1.0  # density
+eta = 1.6  # kinematic viscosity
+width = 41  # of box
+wall_thickness = 2
+actual_width = width - 2*wall_thickness  # subtract boundary layer from box width
+shear_velocity = 0.2 # scale by width to keep stable
+t_max = 2000
+
+
+@pytest.mark.parametrize('target', ('cpu', 'gpu', 'opencl'))
+@pytest.mark.parametrize('stencil_name', ("D2Q9", "D3Q19",))
+def test_poiseuille_channel(target, stencil_name):
+    # OpenCL and Cuda
+    if target == 'opencl':
+        import pytest
+        pytest.importorskip("pyopencl")
+        import pystencils.opencl.autoinit
+    elif target == 'gpu':
+        import pytest
+        pytest.importorskip("pycuda")
+
+    # LB parameters
+    lb_stencil = get_stencil(stencil_name)
+    dim = len(lb_stencil[0])
+
+    if dim == 2:
+        L = [4, width]
+    elif dim == 3:
+        L = [4, width, 4]
+    else:
+        raise Exception()
+    periodicity = [True, False] + [True] * (dim - 2)
+
+    omega = lbmpy.relaxationrates.relaxation_rate_from_lattice_viscosity(eta)
+
+    # ## Data structures
+    dh = ps.create_data_handling(L, periodicity=periodicity, default_target=target)
+
+    src = dh.add_array('src', values_per_cell=len(lb_stencil))
+    dst = dh.add_array_like('dst', 'src')
+    ρ = dh.add_array('rho', latex_name='\\rho')
+    u = dh.add_array('u', values_per_cell=dh.dim)
+    p = dh.add_array('p', values_per_cell=dh.dim**2)
+
+    # LB Setup
+    collision = create_lb_update_rule(stencil=lb_stencil,
+                                      relaxation_rate=omega,
+                                      method="trt",
+                                      compressible=True,
+                                      kernel_type='collide_only',
+                                      optimization={'symbolic_field': src})
+
+    stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'velocity': u})
+
+    opts = {'cpu_openmp': False,
+            'cpu_vectorize_info': None,
+            'target': dh.default_target}
+
+    stream_kernel = ps.create_kernel(stream, **opts).compile()
+    collision_kernel = ps.create_kernel(collision, **opts).compile()
+
+    # Boundaries
+    lbbh = LatticeBoltzmannBoundaryHandling(collision.method, dh, src.name, target=dh.default_target)
+
+    # Second moment test setup
+    cqc = collision.method.conserved_quantity_computation
+    getter_eqs = cqc.output_equations_from_pdfs(src.center_vector,
+                                                {'moment2': p})
+
+    kernel_compute_p = ps.create_kernel(getter_eqs, ghost_layers=0 ).compile()
+
+    # ## Set up the simulation
+
+    init = macroscopic_values_setter(collision.method, velocity=(0,) * dh.dim,
+                                     pdfs=src.center_vector, density=ρ.center)
+    init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
+
+    vel_vec = sp.Matrix([0.5*shear_velocity] + [0] * (dim-1))
+    if dim == 2:
+        lbbh.set_boundary(UBB(velocity=vel_vec), ps.make_slice[:, :wall_thickness])
+        lbbh.set_boundary(UBB(velocity=-1*vel_vec), ps.make_slice[:, -wall_thickness:])
+    elif dim == 3:
+        lbbh.set_boundary(UBB(velocity=vel_vec), ps.make_slice[:, :wall_thickness, :])
+        lbbh.set_boundary(UBB(velocity=-1*vel_vec), ps.make_slice[:, -wall_thickness:, :])
+    else:
+        raise Exception()
+
+    for bh in lbbh, :
+        assert len(bh._boundary_object_to_boundary_info) == 2, "Restart kernel to clear boundaries"
+
+    def init():
+        dh.fill(ρ.name, rho_0)
+        dh.fill(u.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
+        dh.fill(u.name, 0)
+
+        dh.run_kernel(init_kernel)
+
+
+    sync_pdfs = dh.synchronization_function([src.name])
+
+    # Time loop
+    def time_loop(steps):
+        dh.all_to_gpu()
+        for i in range(steps):
+            dh.run_kernel(collision_kernel)
+            sync_pdfs()
+            lbbh()
+            dh.run_kernel(stream_kernel)
+            dh.run_kernel(kernel_compute_p)
+
+            dh.swap(src.name, dst.name)
+
+        if u.name in dh.gpu_arrays:
+            dh.to_cpu(u.name)
+        uu = dh.gather_array(u.name)
+        # average periodic directions
+        if dim == 3:  # dont' swap order
+            uu = np.average(uu, axis=2)
+        uu = np.average(uu, axis=0)
+
+        if p.name in dh.gpu_arrays:
+            dh.to_cpu(p.name)
+        pp = dh.gather_array(p.name)
+        # average periodic directions
+        if dim == 3:  # dont' swap order
+            pp = np.average(pp, axis=2)
+        pp = np.average(pp, axis=0)
+
+
+        # cut off wall regions
+        uu = uu[wall_thickness:-wall_thickness]
+        pp = pp[wall_thickness:-wall_thickness]
+
+        return uu, pp.reshape((len(pp),3,3))
+
+    init()
+    # Simulation
+    profile, pressure_profile = time_loop(t_max)
+
+    # compare against analytical solution
+    # The profile is of shape (n,3). Force is in x-direction
+    y = np.arange(len(profile[:, 0]))
+    mid = (y[-1] - y[0]) / 2  # Mid point of channel
+
+    expected = shear_flow(
+                x=(np.arange(0, actual_width+2) + .5),
+                t=t_max,
+                nu=eta/rho_0,
+                v=shear_velocity,
+                h=actual_width+2,
+                k_max=100)
+
+
+    shear_direction = np.array((1, 0, 0), dtype=float)
+    shear_plane_normal = np.array((0, 1, 0), dtype=float)
+
+    shear_rate = shear_velocity / actual_width
+    dynamic_viscosity = eta * rho_0
+
+    p_expected = rho_0 * np.identity(3) /3.0 - dynamic_viscosity * shear_rate * (
+            np.outer(shear_plane_normal, shear_direction) + np.transpose(np.outer(shear_plane_normal, shear_direction)))
+    pressure_profile[:,0,0] -= rho_0*profile[:,0]**2
+    
+    np.testing.assert_allclose(profile[:, 0], expected[1:-1], atol=1E-9)
+    for i in range(actual_width):
+      np.testing.assert_allclose(pressure_profile[i], p_expected, atol=1E-9, rtol=1E-3)