boundaryconditions.py 4.21 KB
Newer Older
1
from pystencils import Assignment
2
from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
3
from typing import List, Tuple, Any
4
from pystencils.data_types import create_type
5
6


7
class Boundary:
8
9
    """Base class all boundaries should derive from"""

Martin Bauer's avatar
Martin Bauer committed
10
11
12
    inner_or_boundary = True
    single_link = False

13
14
15
    def __init__(self, name=None):
        self._name = name

Martin Bauer's avatar
Martin Bauer committed
16
17
18
19
20
21
22
23
24
    def __call__(self, field, direction_symbol, index_field) -> List[Assignment]:
        """Defines the boundary behavior and must therefore be implemented by all boundaries.

        Here the boundary is defined as a list of sympy assignments, from which a boundary kernel is generated.

        Args:
            field: pystencils field where boundary condition should be applied.
                   The current cell is cell next to the boundary, which is influenced by the boundary
                   cell i.e. has a link from the boundary cell to itself.
Martin Bauer's avatar
Martin Bauer committed
25
            direction_symbol: a sympy symbol that can be used as index to the pdf_field. It describes
Martin Bauer's avatar
Martin Bauer committed
26
27
                              the direction pointing from the fluid to the boundary cell
            index_field: the boundary index field that can be used to retrieve and update boundary data
28
29
30
31
        """
        raise NotImplementedError("Boundary class has to overwrite __call__")

    @property
Martin Bauer's avatar
Martin Bauer committed
32
    def additional_data(self) -> Tuple[str, Any]:
33
        """Return a list of (name, type) tuples for additional data items required in this boundary
Martin Bauer's avatar
Martin Bauer committed
34
35
        These data items can either be initialized in separate kernel see additional_data_kernel_init or by
        Python callbacks - see additional_data_callback """
36
        return ()
37
38

    @property
Martin Bauer's avatar
Martin Bauer committed
39
    def additional_data_init_callback(self):
40
41
42
43
44
45
46
47
48
49
50
51
        """Return a callback function called with a boundary data setter object and returning a dict of
        data-name to data for each element that should be initialized"""
        return None

    @property
    def name(self):
        if self._name:
            return self._name
        else:
            return type(self).__name__

    @name.setter
Martin Bauer's avatar
Martin Bauer committed
52
53
    def name(self, new_value):
        self._name = new_value
54
55
56


class Neumann(Boundary):
Martin Bauer's avatar
Martin Bauer committed
57
58
59
60

    inner_or_boundary = False
    single_link = True

Martin Bauer's avatar
Martin Bauer committed
61
    def __call__(self, field, direction_symbol, **kwargs):
62

Martin Bauer's avatar
Martin Bauer committed
63
64
        neighbor = BoundaryOffsetInfo.offset_from_dir(direction_symbol, field.spatial_dimensions)
        if field.index_dimensions == 0:
Martin Bauer's avatar
Martin Bauer committed
65
            return [Assignment(field.center, field[neighbor])]
66
67
        else:
            from itertools import product
Martin Bauer's avatar
Martin Bauer committed
68
            if not field.has_fixed_index_shape:
69
                raise NotImplementedError("Neumann boundary works only for fields with fixed index shape")
Martin Bauer's avatar
Martin Bauer committed
70
            index_iter = product(*(range(i) for i in field.index_shape))
Martin Bauer's avatar
Martin Bauer committed
71
            return [Assignment(field(*idx), field[neighbor](*idx)) for idx in index_iter]
72
73
74
75
76
77
78

    def __hash__(self):
        # All boundaries of these class behave equal -> should also be equal
        return hash("Neumann")

    def __eq__(self, other):
        return type(other) == Neumann
79
80
81


class Dirichlet(Boundary):
Martin Bauer's avatar
Martin Bauer committed
82
83
84
85

    inner_or_boundary = False
    single_link = True

86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
    def __init__(self, value, name="Dirchlet"):
        super().__init__(name)
        self._value = value

    @property
    def additional_data(self):
        if callable(self._value):
            return [('value', create_type("double"))]
        else:
            return []

    @property
    def additional_data_init_callback(self):
        if callable(self._value):
            return self._value

    def __call__(self, field, direction_symbol, index_field, **kwargs):
Martin Bauer's avatar
Martin Bauer committed
103

104
        if field.index_dimensions == 0:
Martin Bauer's avatar
Martin Bauer committed
105
106
107
108
109
110
111
112
            return [Assignment(field, index_field("value") if self.additional_data else self._value)]
        elif field.index_dimensions == 1:
            assert not self.additional_data
            if not field.has_fixed_index_shape:
                raise NotImplementedError("Field needs fixed index shape")
            assert len(self._value) == field.index_shape[0], "Dirichlet value does not match index shape of field"
            return [Assignment(field(i), self._value[i]) for i in range(field.index_shape[0])]
        raise NotImplementedError("Dirichlet boundary not implemented for fields with more than one index dimension")