test_extrapolation_outflow.py 5.5 KB
Newer Older
1
2
from pystencils.stencil import inverse_direction

3
4
5
6
7
8
9
10
11
12
13
from lbmpy.stencils import get_stencil
from lbmpy.advanced_streaming.utility import AccessPdfValues, get_timesteps
import pytest
import numpy as np

from pystencils.datahandling import create_data_handling
from lbmpy.boundaries import LatticeBoltzmannBoundaryHandling, SimpleExtrapolationOutflow, ExtrapolationOutflow
from lbmpy.creationfunctions import create_lb_method
from lbmpy.advanced_streaming.utility import streaming_patterns
from pystencils.slicing import get_ghost_region_slice

14
15
from itertools import product

16
17
18
19
20
21
22
23
24

@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q27'])
@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
def test_pdf_simple_extrapolation(stencil, streaming_pattern):
    stencil = get_stencil(stencil)
    dim = len(stencil[0])
    values_per_cell = len(stencil)

    #   Field contains exactly one fluid cell
25
    domain_size = (3,) * dim
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
    for timestep in get_timesteps(streaming_pattern):
        dh = create_data_handling(domain_size, default_target='cpu')
        lb_method = create_lb_method(stencil=stencil)
        pdf_field = dh.add_array('f', values_per_cell=values_per_cell)
        dh.fill(pdf_field.name, np.nan, ghost_layers=True)
        bh = LatticeBoltzmannBoundaryHandling(lb_method, dh, pdf_field.name, streaming_pattern, target='cpu')

        #   Set up outflows in all directions
        for normal_dir in stencil[1:]:
            boundary_obj = SimpleExtrapolationOutflow(normal_dir, stencil)
            boundary_slice = get_ghost_region_slice(normal_dir)
            bh.set_boundary(boundary_obj, boundary_slice)

        pdf_arr = dh.cpu_arrays[pdf_field.name]

        #   Set up the domain with artificial PDF values
42
        # center = (1,) * dim
43
        out_access = AccessPdfValues(stencil, streaming_pattern, timestep, 'out')
44
45
46
        for cell in product(*(range(1, 4) for _ in range(dim))):
            for q in range(values_per_cell):
                out_access.write_pdf(pdf_arr, cell, q, q)
47
48
49
50
51
52
53
54

        #   Do boundary handling
        bh(prev_timestep=timestep)

        #   Check PDF values
        in_access = AccessPdfValues(stencil, streaming_pattern, timestep.next(), 'in')

        #   Inbound in center cell
55
56
57
58
        for cell in product(*(range(1, 4) for _ in range(dim))):
            for q in range(values_per_cell):
                f = in_access.read_pdf(pdf_arr, cell, q)
                assert f == q
59
60
61
62
63
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


def test_extrapolation_outflow_initialization_by_copy():
    stencil = get_stencil('D2Q9')
    values_per_cell = len(stencil)
    domain_size = (1, 5)

    streaming_pattern = 'esotwist'
    zeroth_timestep = 'even'

    pdf_acc = AccessPdfValues(stencil, streaming_pattern=streaming_pattern,
                              timestep=zeroth_timestep, streaming_dir='out')

    dh = create_data_handling(domain_size, default_target='cpu')
    lb_method = create_lb_method(stencil=stencil)
    pdf_field = dh.add_array('f', values_per_cell=values_per_cell)
    dh.fill(pdf_field.name, np.nan, ghost_layers=True)
    pdf_arr = dh.cpu_arrays[pdf_field.name]
    bh = LatticeBoltzmannBoundaryHandling(lb_method, dh, pdf_field.name,
                                          streaming_pattern=streaming_pattern, target='cpu')

    for y in range(1, 6):
        for j in range(values_per_cell):
            pdf_acc.write_pdf(pdf_arr, (1, y), j, j)

    normal_dir = (1, 0)
    outflow = ExtrapolationOutflow(normal_dir, lb_method, streaming_pattern=streaming_pattern,
                                   zeroth_timestep=zeroth_timestep)
    boundary_slice = get_ghost_region_slice(normal_dir)
    bh.set_boundary(outflow, boundary_slice)
    bh.prepare()

    blocks = list(dh.iterate())
    index_list = blocks[0][bh._index_array_name].boundary_object_to_index_list[outflow]
93
    assert len(index_list) == 13
94
    for entry in index_list:
95
96
97
98
        direction = stencil[entry["dir"]]
        inv_dir = stencil.index(inverse_direction(direction))
        assert entry[f'pdf'] == inv_dir
        assert entry[f'pdf_nd'] == inv_dir
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118


def test_extrapolation_outflow_initialization_by_callback():
    stencil = get_stencil('D2Q9')
    values_per_cell = len(stencil)
    domain_size = (1, 5)

    streaming_pattern = 'esotwist'
    zeroth_timestep = 'even'

    dh = create_data_handling(domain_size, default_target='cpu')
    lb_method = create_lb_method(stencil=stencil)
    pdf_field = dh.add_array('f', values_per_cell=values_per_cell)
    dh.fill(pdf_field.name, np.nan, ghost_layers=True)
    bh = LatticeBoltzmannBoundaryHandling(lb_method, dh, pdf_field.name,
                                          streaming_pattern=streaming_pattern, target='cpu')

    normal_dir = (1, 0)
    density_callback = lambda x, y: 1
    velocity_callback = lambda x, y: (0, 0)
119
    outflow = ExtrapolationOutflow(normal_dir, lb_method,
120
121
122
123
124
125
126
127
128
129
130
131
                                   streaming_pattern=streaming_pattern,
                                   zeroth_timestep=zeroth_timestep,
                                   initial_density=density_callback,
                                   initial_velocity=velocity_callback)
    boundary_slice = get_ghost_region_slice(normal_dir)
    bh.set_boundary(outflow, boundary_slice)
    bh.prepare()

    weights = [w.evalf() for w in lb_method.weights]

    blocks = list(dh.iterate())
    index_list = blocks[0][bh._index_array_name].boundary_object_to_index_list[outflow]
132
    assert len(index_list) == 13
133
    for entry in index_list:
134
135
136
        direction = stencil[entry["dir"]]
        inv_dir = stencil.index(inverse_direction(direction))
        assert entry[f'pdf_nd'] == weights[inv_dir]