test_flexible_noslip.py 3.75 KB
Newer Older
Frederik Hennig's avatar
Frederik Hennig committed
1
2
3
4
import numpy as np
import sympy as sp

from lbmpy.fieldaccess import *
5
from lbmpy.advanced_streaming.flexible_boundaries_prototype import FlexibleNoSlip, create_flexible_lbm_boundary_kernel
Frederik Hennig's avatar
Frederik Hennig committed
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
52
53
54
55
56
57
58
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
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
121
122
123
124
125
126
from lbmpy.creationfunctions import create_lb_method
from lbmpy.stencils import get_stencil

import pystencils as ps

from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
from pystencils.data_types import TypedSymbol, create_type
from pystencils.field import Field, FieldType

from itertools import product


#   ===========================
#         INFRASTRUCTURE
#   ===========================


timesteps = ['odd_to_even', 'even_to_odd']
kernel_types = ['push', 'pull', 'aa', 'esotwist']

even_accessors = {
    'pull': StreamPullTwoFieldsAccessor,
    'push': StreamPushTwoFieldsAccessor,
    'aa': AAEvenTimeStepAccessor,
    'esotwist': EsoTwistEvenTimeStepAccessor
}

odd_accessors = {
    'pull': StreamPullTwoFieldsAccessor,
    'push': StreamPushTwoFieldsAccessor,
    'aa': AAOddTimeStepAccessor,
    'esotwist': EsoTwistOddTimeStepAccessor
}


def inverse_dir_index(stencil, dir):
    return stencil.index(tuple(-d for d in stencil[dir]))

"""
    Allows to access values from a PDF array correctly depending on 
    the streaming pattern.
"""
class AccessPdfValues:
    def __init__(self, pdf_field, stencil, accessor):
        self.read_accs = accessor.read(pdf_field, stencil)
        self.write_accs = accessor.write(pdf_field, stencil)

    def write_pdf(self, pdf_arr, x, y, d, value):
        offsets = self.write_accs[d].offsets
        x += offsets[0]
        y += offsets[1]
        i = self.write_accs[d].index[0]

        pdf_arr[x,y,i] = value

    def read_pdf(self, pdf_arr, x, y, d):
        offsets = self.read_accs[d].offsets
        x += offsets[0]
        y += offsets[1]
        i = self.read_accs[d].index[0]

        return pdf_arr[x,y,i]



#   ===========================
#              TESTS
#   ===========================

import pytest

arguments = product(timesteps, kernel_types)

@pytest.mark.parametrize("timestep_type, kernel_type", arguments)
def test_flexible_noslip(timestep_type, kernel_type):
    """
    Flexible NoSlip Test
    """

    stencil = get_stencil('D2Q9')
    pdf_field = ps.fields('pdfs(9): [2D]')
    q = len(stencil)

    even_accessor = even_accessors[kernel_type]
    odd_accessor = odd_accessors[kernel_type]

    prev_accessor = (odd_accessor if timestep_type == 'odd_to_even' else even_accessor)
    next_accessor = (even_accessor if timestep_type == 'odd_to_even' else odd_accessor)

    prev_pdf_access = AccessPdfValues(pdf_field, stencil, prev_accessor)
    next_pdf_access = AccessPdfValues(pdf_field, stencil, next_accessor)

    pdfs = np.zeros((3,3,9))
    for d in range(q):
        prev_pdf_access.write_pdf(pdfs, 1, 1, d, d)

    lb_method = create_lb_method()
    flex_noslip = FlexibleNoSlip()

    index_struct_dtype = numpy_data_type_for_boundary_object(flex_noslip, lb_method.dim)

    index_field = Field('indexVector', FieldType.INDEXED, index_struct_dtype, layout=[0],
                    shape=(TypedSymbol("indexVectorSize", create_type(np.int64)), 1), strides=(1, 1))
    index_vector = np.array([ (1,1,d) for d in range(q) ], dtype=index_struct_dtype)

    flex_ast = create_flexible_lbm_boundary_kernel(pdf_field, 
                    index_field, lb_method, flex_noslip,
                    between_timesteps=timestep_type,
                    kernel_type=kernel_type)

    flex_kernel = flex_ast.compile()
    
    flex_kernel(pdfs=pdfs, indexVector=index_vector, indexVectorSize=len(index_vector))

    reflected_pdfs = [ next_pdf_access.read_pdf(pdfs, 1, 1, d) for d in range(q)] 

    for d, f in enumerate(reflected_pdfs):
        if int(f) != inverse_dir_index(stencil, d):
            return False

    return True