Commit a04d5150 authored by itischler's avatar itischler
Browse files

formatting

parent 15698e80
Pipeline #33325 failed with stages
in 12 minutes and 31 seconds
"""
Test Poiseuille flow against analytical solution
Test shear flow velocity and pressureagainst analytical solutions
"""
......@@ -54,7 +54,7 @@ eta = 1.6 # kinematic viscosity
width = 40 # of box
wall_thickness = 2
actual_width = width - wall_thickness # subtract boundary layer from box width
shear_velocity = 0.001 # scale by width to keep stable
shear_velocity = 0.2 # scale by width to keep stable
t_max = 20000
......@@ -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,28 +176,26 @@ 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]
if(dh.dim == 2):
pp = pp.reshape((len(pp),2,2))
pp = pp.reshape((len(pp), 2, 2))
if(dh.dim == 3):
pp = pp.reshape((len(pp),3,3))
pp = pp.reshape((len(pp), 3, 3))
return uu, pp
init()
# Simulation
profile, pressure_profile = time_loop(t_max)
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)
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)
if(dh.dim == 2):
shear_direction = np.array((1, 0), dtype=float)
......@@ -211,12 +208,12 @@ def test_poiseuille_channel(target, stencil_name):
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(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)))
# Sustract the tensorproduct of the velosity to get the pressure
pressure_profile[:,0,0] -= rho_0*profile[:,0]**2
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-2):
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