Pressure tensor test

lbmpy_tests/test_shear_flow.py

 """
-Test Poiseuille flow against analytical solution
+Test shear flow velocity and pressureagainst analytical solutions
 """
 
 
@@ -49,13 +49,13 @@ def shear_flow(x, t, nu, v, h, k_max):
     return -v * u
 
 
-rho_0 = 1.0  # density
+rho_0 = 2.2  # density
 eta = 1.6  # kinematic viscosity
-width = 41  # of box
+width = 40  # 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
+actual_width = width - wall_thickness  # subtract boundary layer from box width
+shear_velocity = 0.2  # scale by width to keep stable
+t_max = 20000
 
 
 @pytest.mark.parametrize('target', ('cpu', 'gpu', 'opencl'))
@@ -118,7 +118,7 @@ def test_poiseuille_channel(target, stencil_name):
     getter_eqs = cqc.output_equations_from_pdfs(src.center_vector,
                                                 {'moment2': p})
 
-    kernel_compute_p = ps.create_kernel(getter_eqs, ghost_layers=0 ).compile()
+    kernel_compute_p = ps.create_kernel(getter_eqs, ghost_layers=0).compile()
 
     # ## Set up the simulation
 
@@ -126,13 +126,13 @@ def test_poiseuille_channel(target, stencil_name):
                                      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))
+    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:])
+        lbbh.set_boundary(UBB(velocity=-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:, :])
+        lbbh.set_boundary(UBB(velocity=-vel_vec), ps.make_slice[:, -wall_thickness:, :])
     else:
         raise Exception()
 
@@ -146,7 +146,6 @@ def test_poiseuille_channel(target, stencil_name):
 
         dh.run_kernel(init_kernel)
 
-
     sync_pdfs = dh.synchronization_function([src.name])
 
     # Time loop
@@ -177,41 +176,44 @@ def test_poiseuille_channel(target, stencil_name):
             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))
+        if(dh.dim == 2):
+            pp = pp.reshape((len(pp), 2, 2))
+        if(dh.dim == 3):
+            pp = pp.reshape((len(pp), 3, 3))
+        return uu, pp
 
     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)
+    expected = shear_flow(x=(np.arange(0, actual_width) + .5),
+                          t=t_max,
+                          nu=eta / rho_0,
+                          v=shear_velocity,
+                          h=actual_width,
+                          k_max=100)
 
-
-    shear_direction = np.array((1, 0, 0), dtype=float)
-    shear_plane_normal = np.array((0, 1, 0), dtype=float)
+    if(dh.dim == 2):
+        shear_direction = np.array((1, 0), dtype=float)
+        shear_plane_normal = np.array((0, 1), dtype=float)
+    if(dh.dim == 3):
+        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
+    correction_factor = eta / (eta - 1. / 6.)
+
+    p_expected = rho_0 * np.identity(dh.dim) / 3.0 + dynamic_viscosity * shear_rate / correction_factor * (
+        np.outer(shear_plane_normal, shear_direction) + np.transpose(np.outer(shear_plane_normal, shear_direction)))
 
-    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
+    # Sustract the tensorproduct of the velosity to get the pressure
+    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)
+    for i in range(actual_width - 2):
+        np.testing.assert_allclose(pressure_profile[i], p_expected, atol=1E-9, rtol=1E-3)