Merge branch 'poiseuille' into 'master'

Test Poiseuille channel against analytical solution for several stencils

See merge request !73

lbmpy_tests/test_poisuille_channel.py

 """
-This test revealed a problem with OpenCL (https://i10git.cs.fau.de/pycodegen/lbmpy/issues/9#note_9521)
+Test Poiseuille flow against analytical solution
 """
 
 
 import numpy as np
 import pytest
+import sympy as sp
 
 import lbmpy
 from lbmpy.boundaries import NoSlip
 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.session import *
-from pystencils.session import *
+from lbmpy.stencils import get_stencil
 
-previous_vmids = {}
-previous_endresults = {}
+import pystencils as ps
+
+
+def poiseuille_flow(z, H, ext_force_density, dyn_visc):
+    """
+    Analytical solution for plane Poiseuille flow.
+
+    Parameters
+    ----------
+    z : :obj:`float`
+        Distance to the mid plane of the channel.
+    H : :obj:`float`
+        Distance between the boundaries.
+    ext_force_density : :obj:`float`
+        Force density on the fluid normal to the boundaries.
+    dyn_visc : :obj:`float`
+        Dynamic viscosity of the LB fluid.
+
+    """
+    return ext_force_density * 1. / (2 * dyn_visc) * (H**2.0 / 4.0 - z**2.0)
+
+
+rho_0 = 1.2  # density
+eta = 0.2  # kinematic viscosity
+width = 41  # of box
+actual_width = width - 2  # subtract boundary layer from box width
+ext_force_density = 0.2 / actual_width**2  # scale by width to keep stable
 
 
 @pytest.mark.parametrize('target', ('cpu', 'gpu', 'opencl'))
-def test_poiseuille_channel(target):
-    # # Lattice Boltzmann
-    #
-    # ## Definitions
+@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")
@@ -30,65 +55,72 @@ def test_poiseuille_channel(target):
         import pytest
         pytest.importorskip("pycuda")
 
-    # In[2]:
+    # LB parameters
+    lb_stencil = get_stencil(stencil_name)
+    dim = len(lb_stencil[0])
 
-    L = (34, 34)
+    if dim == 2:
+        L = [4, width]
+    elif dim == 3:
+        L = [4, width, 4]
+    else:
+        raise Exception()
+    periodicity = [True, False] + [True] * (dim - 2)
 
-    lb_stencil = get_stencil("D2Q9")
-    eta = 1
     omega = lbmpy.relaxationrates.relaxation_rate_from_lattice_viscosity(eta)
 
-    f_pre = 0.00001
-
     # ## Data structures
-
-    # In[3]:
-
-    dh = ps.create_data_handling(L, periodicity=(True, False), default_target=target)
-
-    opts = {'cpu_openmp': True,
-            'cpu_vectorize_info': None,
-            'target': dh.default_target}
+    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)
 
-    # In[4]:
-
+    # LB Setup
     collision = create_lb_update_rule(stencil=lb_stencil,
                                       relaxation_rate=omega,
+                                      method="trt",
                                       compressible=True,
-                                      force_model='guo',
-                                      force=sp.Matrix([f_pre, 0]),
+                                      force_model="schiller",
+                                      force=sp.Matrix([ext_force_density] + [0] * (dim - 1)),
                                       kernel_type='collide_only',
                                       optimization={'symbolic_field': src})
 
     stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'velocity': u})
 
-    lbbh = LatticeBoltzmannBoundaryHandling(collision.method, dh, src.name, target=dh.default_target)
+    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()
 
-    # ## Set up the simulation
+    # Boundaries
+    lbbh = LatticeBoltzmannBoundaryHandling(collision.method, dh, src.name, target=dh.default_target)
 
-    # In[5]:
+    # ## Set up the simulation
 
-    init = macroscopic_values_setter(collision.method, velocity=(0,)*dh.dim,
+    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()
 
     noslip = NoSlip()
-    lbbh.set_boundary(noslip, make_slice[:, :4])
-    lbbh.set_boundary(noslip, make_slice[:, -4:])
+    wall_thickness = 2
+    if dim == 2:
+        lbbh.set_boundary(noslip, ps.make_slice[:, :wall_thickness])
+        lbbh.set_boundary(noslip, ps.make_slice[:, -wall_thickness:])
+    elif dim == 3:
+        lbbh.set_boundary(noslip, ps.make_slice[:, :wall_thickness, :])
+        lbbh.set_boundary(noslip, ps.make_slice[:, -wall_thickness:, :])
+    else:
+        raise Exception()
 
     for bh in lbbh, :
         assert len(bh._boundary_object_to_boundary_info) == 1, "Restart kernel to clear boundaries"
 
     def init():
-        dh.fill(ρ.name, 1)
+        dh.fill(ρ.name, rho_0)
         dh.fill(u.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
         dh.fill(u.name, 0)
 
@@ -98,10 +130,11 @@ def test_poiseuille_channel(target):
 
     sync_pdfs = dh.synchronization_function([src.name])
 
+    # Time loop
     def time_loop(steps):
         dh.all_to_gpu()
-        vmid = np.empty((2, steps//10+1))
         i = -1
+        last_max_vel = -1
         for i in range(steps):
             dh.run_kernel(collision_kernel)
             sync_pdfs()
@@ -110,71 +143,41 @@ def test_poiseuille_channel(target):
 
             dh.swap(src.name, dst.name)
 
-            if i % 10 == 0:
+            # Consider early termination
+            if i % 100 == 0:
                 if u.name in dh.gpu_arrays:
                     dh.to_cpu(u.name)
                 uu = dh.gather_array(u.name)
-                uu = uu[L[0]//2-1:L[0]//2+1, L[1]//2-1:L[1]//2+1, 0].mean()
-                vmid[:, i//10] = [i, uu]
-                if 1/np.sqrt(3) < uu:
-                    break
-        if dh.is_on_gpu(u.name):
-            dh.to_cpu(u.name)
-
-        return vmid[:, :i//10+1]
+                # average periodic directions
+                if dim == 3:  # dont' swap order
+                    uu = np.average(uu, axis=2)
+                uu = np.average(uu, axis=0)
 
-    # In[7]:
-
-    # def plot():
-
-        # plt.subplot(2, 2, 1)
-        # plt.title("$u$")
-        # plt.xlabel("$x$")
-        # plt.ylabel("$y$")
-        # plt.vector_field_magnitude(uu)
-        # plt.colorbar()
-
-        # plt.subplot(2, 2, 2)
-        # plt.title("$u$")
-        # plt.xlabel("$x/2$")
-        # plt.ylabel("$y/2$")
-        # plt.vector_field(uu, step=2)
-
-        # actualwidth = np.sum(1-np.isnan(uu[0, :, 0]))
-        # uref = f_pre*actualwidth**2/(8*(eta))
-
-        # plt.subplot(2, 2, 3)
-        # plt.title("flow profile")
-        # plt.xlabel("$y$")
-        # plt.ylabel(r"$u_x$")
-        # plt.plot((uu[L[0]//2-1, :, 0]+uu[L[0]//2, :, 0])/2)
+                max_vel = np.nanmax(uu)
+                if np.abs(max_vel / last_max_vel - 1) < 5E-6:
+                    break
+                last_max_vel = max_vel
 
-        # plt.subplot(2, 2, 4)
-        # plt.title("convergence")
-        # plt.xlabel("$t$")
-        # plt.plot()
+        # cut off wall regions
+        uu = uu[wall_thickness:-wall_thickness]
 
-    # ## Run the simulation
+        # correct for f/2 term
+        uu -= np.array([ext_force_density / 2 / rho_0] + [0] * (dim - 1))
 
-    # In[8]:
+        return uu
 
     init()
-    vmid = time_loop(1000)
-    # plot()
-
-    uu = dh.gather_array(u.name)
+    # Simulation
+    profile = time_loop(5000)
 
-    for target, prev_endresult in previous_endresults.items():
-        assert np.allclose(uu[4:-4, 4:-4, 0], prev_endresult[4:-4, 4:-4, 0],
-                           atol=1e-5), f'uu does not agree with result from {target}'
+    # 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
 
-    for target, prev_vmid in previous_vmids.items():
-        assert np.allclose(vmid, prev_vmid, atol=1e-5), f'vmid does not agree with result from {target}'
+    expected = poiseuille_flow((y - mid), actual_width, ext_force_density, rho_0 * eta)
 
-    # uref = f_pre*actualwidth**2/(8*(eta))
-    # actualwidth = np.sum(1-np.isnan(uu[0, :, 0]))
-    # assert np.allclose(vmid[1, -1]/uref, 1, atol=1e-3)
-    # print(vmid[1, -1])
+    np.testing.assert_allclose(profile[:, 0], expected, rtol=0.006)
 
-    previous_vmids[target] = vmid
-    previous_endresults[target] = uu
+    # Test zero vel in other directions
+    np.testing.assert_allclose(profile[:, 1:], np.zeros_like(profile[:, 1:]), atol=1E-9)