boundaryconditions.py 11.8 KB
Newer Older
Martin Bauer's avatar
Martin Bauer committed
1
import sympy as sp
Martin Bauer's avatar
Martin Bauer committed
2
from pystencils import Assignment, Field
3
from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
Martin Bauer's avatar
Martin Bauer committed
4
from pystencils.data_types import create_type
Martin Bauer's avatar
Martin Bauer committed
5
from pystencils.sympyextensions import get_symmetric_part
6
from lbmpy.simplificationfactory import create_simplification_strategy
Frederik Hennig's avatar
Frederik Hennig committed
7
from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays
Martin Bauer's avatar
Martin Bauer committed
8
9


10
11
class LbBoundary:
    """Base class that all boundaries should derive from"""
Martin Bauer's avatar
Martin Bauer committed
12

Martin Bauer's avatar
Martin Bauer committed
13
14
15
    inner_or_boundary = True
    single_link = False

16
17
18
    def __init__(self, name=None):
        self._name = name

19
    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
20
21
        """
        This function defines the boundary behavior and must therefore be implemented by all boundaries.
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
        The boundary is defined through a list of sympy equations from which a boundary kernel is generated.


        :param f_out: a pystencils field acting as a proxy to access the populations streaming out of the current
                cell, i.e. the post-collision PDFs of the previous LBM step
        :param f_in: a pystencils field acting as a proxy to access the populations streaming into the current
                cell, i.e. the pre-collision PDFs for the next LBM step
        :param dir_symbol: a sympy symbol that can be used as an index to f_out and f_in. It describes the direction
                    pointing from the fluid to the boundary cell.
        :param inv_dir: an indexed sympy symbol which describes the inversion of a direction index. It can be used in
                the indices of f_out and f_in for retrieving the PDF of the inverse direction.
        :param lb_method: an instance of the LB method used. Use this to adapt the boundary to the method
                    (e.g. compressibility)
        :param index_field: the boundary index field that can be used to retrieve and update boundary data

        :return: list of pystencils assignments, or pystencils.AssignmentCollection
38
39
40
41
        """
        raise NotImplementedError("Boundary class has to overwrite __call__")

    @property
Martin Bauer's avatar
Martin Bauer committed
42
    def additional_data(self):
43
        """Return a list of (name, type) tuples for additional data items required in this boundary
Martin Bauer's avatar
Martin Bauer committed
44
45
        These data items can either be initialized in separate kernel see additional_data_kernel_init or by
        Python callbacks - see additional_data_callback """
46
47
        return []

48
    @property
Martin Bauer's avatar
Martin Bauer committed
49
    def additional_data_init_callback(self):
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"""
52
53
        return None

54
55
56
57
    def get_additional_code_nodes(self, lb_method):
        """Return a list of code nodes that will be added in the generated code before the index field loop."""
        return []

58
59
    @property
    def name(self):
60
61
62
63
64
65
        if self._name:
            return self._name
        else:
            return type(self).__name__

    @name.setter
Martin Bauer's avatar
Martin Bauer committed
66
67
    def name(self, new_value):
        self._name = new_value
68

69
70
# end class Boundary

71

72
class NoSlip(LbBoundary):
73
74
75

    def __init__(self, name=None):
        """Set an optional name here, to mark boundaries, for example for force evaluations"""
76
        super(NoSlip, self).__init__(name)
77

78
79
80
81
82
83
84
    """
        No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle.
        Extended for use with any streaming pattern.
    """

    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
        return Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol))
85
86

    def __hash__(self):
87
        return hash(self.name)
88
89

    def __eq__(self, other):
90
91
92
93
        if not isinstance(other, NoSlip):
            return False
        return self.name == other.name

94
95
# end class NoSlip

96

97
class UBB(LbBoundary):
Martin Bauer's avatar
Martin Bauer committed
98
99
    """Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""

Martin Bauer's avatar
Martin Bauer committed
100
    def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None):
101
        """
102
103
104
105
106
        Args:
            velocity: can either be a constant, an access into a field, or a callback function.
                      The callback functions gets a numpy record array with members, 'x','y','z', 'dir' (direction)
                      and 'velocity' which has to be set to the desired velocity of the corresponding link
            adapt_velocity_to_force:
107
        """
108
        super(UBB, self).__init__(name)
109
        self._velocity = velocity
Martin Bauer's avatar
Martin Bauer committed
110
        self._adaptVelocityToForce = adapt_velocity_to_force
111
112
113
114
115
        if callable(self._velocity) and not dim:
            raise ValueError("When using a velocity callback the dimension has to be specified with the dim parameter")
        elif not callable(self._velocity):
            dim = len(velocity)
        self.dim = dim
Martin Bauer's avatar
Martin Bauer committed
116

117
    @property
Martin Bauer's avatar
Martin Bauer committed
118
    def additional_data(self):
119
        if callable(self._velocity):
Martin Bauer's avatar
Martin Bauer committed
120
            return [('vel_%d' % (i,), create_type("double")) for i in range(self.dim)]
121
122
123
124
        else:
            return []

    @property
Martin Bauer's avatar
Martin Bauer committed
125
    def additional_data_init_callback(self):
126
127
128
        if callable(self._velocity):
            return self._velocity

129
    def get_additional_code_nodes(self, lb_method):
Frederik Hennig's avatar
Frederik Hennig committed
130
        return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)]
131
132

    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
Martin Bauer's avatar
Martin Bauer committed
133
        vel_from_idx_field = callable(self._velocity)
134
135
        vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity
        direction = dir_symbol
136

Frederik Hennig's avatar
Frederik Hennig committed
137
        assert self.dim == lb_method.dim, \
138
            f"Dimension of UBB ({self.dim}) does not match dimension of method ({lb_method.dim})"
139

Frederik Hennig's avatar
Frederik Hennig committed
140
        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(direction, lb_method.stencil)
141

142
143
144
        velocity = tuple(v_i.get_shifted(*neighbor_offset)
                         if isinstance(v_i, Field.Access) and not vel_from_idx_field
                         else v_i
145
                         for v_i in vel)
146
147

        if self._adaptVelocityToForce:
Martin Bauer's avatar
Martin Bauer committed
148
149
150
            cqc = lb_method.conserved_quantity_computation
            shifted_vel_eqs = cqc.equilibrium_input_equations_from_init_values(velocity=velocity)
            velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.first_order_moment_symbols).main_assignments]
151
152

        c_s_sq = sp.Rational(1, 3)
Martin Bauer's avatar
Martin Bauer committed
153
        weight_of_direction = LbmWeightInfo.weight_of_direction
154
        vel_term = 2 / c_s_sq \
Frederik Hennig's avatar
Frederik Hennig committed
155
            * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) \
Frederik Hennig's avatar
Frederik Hennig committed
156
            * weight_of_direction(direction, lb_method)
157
158

        # Better alternative: in conserved value computation
Martin Bauer's avatar
Martin Bauer committed
159
        # rename what is currently called density to "virtual_density"
160
        # provide a new quantity density, which is constant in case of incompressible models
Martin Bauer's avatar
Martin Bauer committed
161
162
163
        if not lb_method.conserved_quantity_computation.zero_centered_pdfs:
            cqc = lb_method.conserved_quantity_computation
            density_symbol = sp.Symbol("rho")
164
            pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))]
Martin Bauer's avatar
Martin Bauer committed
165
166
167
            density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol})
            density_symbol = lb_method.conserved_quantity_computation.defined_symbols()['density']
            result = density_equations.all_assignments
168
169
            result += [Assignment(f_in(inv_dir[direction]),
                                  f_out(direction) - vel_term * density_symbol)]
170
171
            return result
        else:
172
173
174
            return [Assignment(f_in(inv_dir[direction]),
                               f_out(direction) - vel_term)]
# end class UBB
175

Frederik Hennig's avatar
Frederik Hennig committed
176

177
class FixedDensity(LbBoundary):
178

179
    def __init__(self, density, name=None):
180
181
        if name is None:
            name = "Fixed Density " + str(density)
182
        super(FixedDensity, self).__init__(name)
183
        self._density = density
184

Frederik Hennig's avatar
Frederik Hennig committed
185
    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
186
187
        """Boundary condition that fixes the density/pressure at the obstacle"""

Martin Bauer's avatar
Martin Bauer committed
188
189
190
191
        def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom):
            new_main_assignments = [Assignment(a.lhs, get_symmetric_part(a.rhs, degrees_of_freedom))
                                    for a in assignment_collection.main_assignments]
            return assignment_collection.copy(new_main_assignments)
192

Martin Bauer's avatar
Martin Bauer committed
193
        cqc = lb_method.conserved_quantity_computation
Martin Bauer's avatar
Martin Bauer committed
194
        velocity = cqc.defined_symbols()['velocity']
Martin Bauer's avatar
Martin Bauer committed
195
196
        symmetric_eq = remove_asymmetric_part_of_main_assignments(lb_method.get_equilibrium(),
                                                                  degrees_of_freedom=velocity)
Frederik Hennig's avatar
Frederik Hennig committed
197
        substitutions = {sym: f_out(i) for i, sym in enumerate(lb_method.pre_collision_pdf_symbols)}
Martin Bauer's avatar
Martin Bauer committed
198
        symmetric_eq = symmetric_eq.new_with_substitutions(substitutions)
199

Martin Bauer's avatar
Martin Bauer committed
200
201
        simplification = create_simplification_strategy(lb_method)
        symmetric_eq = simplification(symmetric_eq)
202

Martin Bauer's avatar
Martin Bauer committed
203
        density_symbol = cqc.defined_symbols()['density']
204
205

        density = self._density
Martin Bauer's avatar
Martin Bauer committed
206
207
        equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density)
        equilibrium_input = equilibrium_input.new_without_subexpressions()
Martin Bauer's avatar
Martin Bauer committed
208
209
210
        density_eq = equilibrium_input.main_assignments[0]
        assert density_eq.lhs == density_symbol
        transformed_density = density_eq.rhs
211

Frederik Hennig's avatar
Frederik Hennig committed
212
        conditions = [(eq_i.rhs, sp.Equality(dir_symbol, i))
Martin Bauer's avatar
Martin Bauer committed
213
                      for i, eq_i in enumerate(symmetric_eq.main_assignments)] + [(0, True)]
214
215
        eq_component = sp.Piecewise(*conditions)

Martin Bauer's avatar
Martin Bauer committed
216
217
        subexpressions = [Assignment(eq.lhs, transformed_density if eq.lhs == density_symbol else eq.rhs)
                          for eq in symmetric_eq.subexpressions]
Martin Bauer's avatar
Martin Bauer committed
218

Frederik Hennig's avatar
Frederik Hennig committed
219
220
221
        return subexpressions + [Assignment(f_in(inv_dir[dir_symbol]),
                                            2 * eq_component - f_out(dir_symbol))]
# end class FixedDensity
Martin Bauer's avatar
Martin Bauer committed
222

223

224
class NeumannByCopy(LbBoundary):
Frederik Hennig's avatar
Frederik Hennig committed
225
226
227
228
229

    def get_additional_code_nodes(self, lb_method):
        return [NeighbourOffsetArrays(lb_method.stencil)]

    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
Frederik Hennig's avatar
Frederik Hennig committed
230
        neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
Frederik Hennig's avatar
Frederik Hennig committed
231
232
        return [Assignment(f_in(inv_dir[dir_symbol]), f_out(inv_dir[dir_symbol])),
                Assignment(f_out[neighbour_offset](dir_symbol), f_out(dir_symbol))]
233
234
235
236
237
238
239

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

    def __eq__(self, other):
        return type(other) == NeumannByCopy
Frederik Hennig's avatar
Frederik Hennig committed
240
# end class NeumannByCopy
Martin Bauer's avatar
Martin Bauer committed
241

Martin Bauer's avatar
Martin Bauer committed
242

243
class StreamInConstant(LbBoundary):
244
245
246
    def __init__(self, constant, name=None):
        super(StreamInConstant, self).__init__(name)
        self._constant = constant
Martin Bauer's avatar
Martin Bauer committed
247

Frederik Hennig's avatar
Frederik Hennig committed
248
249
250
251
    def get_additional_code_nodes(self, lb_method):
        return [NeighbourOffsetArrays(lb_method.stencil)]

    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
Frederik Hennig's avatar
Frederik Hennig committed
252
        neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
Frederik Hennig's avatar
Frederik Hennig committed
253
254
        return [Assignment(f_in(inv_dir[dir_symbol]), self._constant),
                Assignment(f_out[neighbour_offset](dir_symbol), self._constant)]
Martin Bauer's avatar
Martin Bauer committed
255
256

    def __hash__(self):
257
258
        # All boundaries of these class behave equal -> should also be equal
        return hash("StreamInConstant")
Martin Bauer's avatar
Martin Bauer committed
259
260

    def __eq__(self, other):
261
        return type(other) == StreamInConstant
Frederik Hennig's avatar
Frederik Hennig committed
262
# end class StreamInConstant
263
264
265
266


#   ------------------------- Old, Deprecated Implementation -------------------------

Frederik Hennig's avatar
Frederik Hennig committed
267
class Boundary(LbBoundary):
268
269
270
271
272
273
274
275

    def __init__(self, name=None):
        from lbmpy.boundaries.boundaryhandling import deprecation_message
        deprecation_message()
        self._name = name

    def __call__(self, pdf_field, direction_symbol, lb_method, index_field):
        raise NotImplementedError("Boundary class has to overwrite __call__")