Commit 3a71611a authored by itischler's avatar itischler
Browse files

Merge branch 'pressure_tensor_test' into 'pressure_tensor'

Pressure tensor test

See merge request !2
parents d77a9064 a04d5150
Pipeline #33326 failed with stages
in 13 minutes and 11 seconds
"""
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)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment