test_shear_flow.py 7.13 KB
Newer Older
1
"""
itischler's avatar
itischler committed
2
Test shear flow velocity and pressureagainst analytical solutions
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
"""


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


itischler's avatar
itischler committed
52
rho_0 = 2.2  # density
53
eta = 1.6  # kinematic viscosity
itischler's avatar
itischler committed
54
width = 40  # of box
55
wall_thickness = 2
itischler's avatar
itischler committed
56
actual_width = width - wall_thickness  # subtract boundary layer from box width
itischler's avatar
itischler committed
57
shear_velocity = 0.2  # scale by width to keep stable
itischler's avatar
itischler committed
58
t_max = 2000
59
60
61
62


@pytest.mark.parametrize('target', ('cpu', 'gpu', 'opencl'))
@pytest.mark.parametrize('stencil_name', ("D2Q9", "D3Q19",))
itischler's avatar
itischler committed
63
def test_shear_flow(target, stencil_name):
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
    # 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})

itischler's avatar
itischler committed
121
    kernel_compute_p = ps.create_kernel(getter_eqs, **opts).compile()
122
123
124
125
126
127
128

    # ## 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()

itischler's avatar
itischler committed
129
    vel_vec = sp.Matrix([0.5 * shear_velocity] + [0] * (dim - 1))
130
131
    if dim == 2:
        lbbh.set_boundary(UBB(velocity=vel_vec), ps.make_slice[:, :wall_thickness])
itischler's avatar
itischler committed
132
        lbbh.set_boundary(UBB(velocity=-vel_vec), ps.make_slice[:, -wall_thickness:])
133
134
    elif dim == 3:
        lbbh.set_boundary(UBB(velocity=vel_vec), ps.make_slice[:, :wall_thickness, :])
itischler's avatar
itischler committed
135
        lbbh.set_boundary(UBB(velocity=-vel_vec), ps.make_slice[:, -wall_thickness:, :])
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
    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]

itischler's avatar
itischler committed
183
        if(dh.dim == 2):
itischler's avatar
itischler committed
184
            pp = pp.reshape((len(pp), 2, 2))
itischler's avatar
itischler committed
185
        if(dh.dim == 3):
itischler's avatar
itischler committed
186
            pp = pp.reshape((len(pp), 3, 3))
itischler's avatar
itischler committed
187
        return uu, pp
188
189
190
191
192

    init()
    # Simulation
    profile, pressure_profile = time_loop(t_max)

itischler's avatar
itischler committed
193
194
195
196
197
198
    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)
199

itischler's avatar
itischler committed
200
201
202
203
204
205
    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)
206
207
208

    shear_rate = shear_velocity / actual_width
    dynamic_viscosity = eta * rho_0
itischler's avatar
itischler committed
209
    correction_factor = eta / (eta - 1. / 6.)
210

itischler's avatar
itischler committed
211
212
    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)))
itischler's avatar
itischler committed
213
214

    # Sustract the tensorproduct of the velosity to get the pressure
itischler's avatar
itischler committed
215
    pressure_profile[:, 0, 0] -= rho_0 * profile[:, 0]**2
216
217
    
    np.testing.assert_allclose(profile[:, 0], expected[1:-1], atol=1E-9)
itischler's avatar
itischler committed
218
219
    for i in range(actual_width - 2):
        np.testing.assert_allclose(pressure_profile[i], p_expected, atol=1E-9, rtol=1E-3)