boundaryconditions.py 23 KB
Newer Older
1
2
from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
from pystencils.simp.assignment_collection import AssignmentCollection
Martin Bauer's avatar
Martin Bauer committed
3
from pystencils import Assignment, Field
Frederik Hennig's avatar
Frederik Hennig committed
4
from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
Martin Bauer's avatar
Martin Bauer committed
5
from pystencils.data_types import create_type
Martin Bauer's avatar
Martin Bauer committed
6
from pystencils.sympyextensions import get_symmetric_part
Frederik Hennig's avatar
Frederik Hennig committed
7
8
from lbmpy.simplificationfactory import create_simplification_strategy
from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays
9
from pystencils.stencil import offset_to_direction_string, direction_string_to_offset, inverse_direction
Martin Bauer's avatar
Martin Bauer committed
10

Frederik Hennig's avatar
Frederik Hennig committed
11
12
import sympy as sp

Martin Bauer's avatar
Martin Bauer committed
13

Frederik Hennig's avatar
Frederik Hennig committed
14
class LbBoundary:
15
16
17
18
19
    """Base class that all boundaries should derive from.

    Args:
        name: optional name of the boundary.
    """
Martin Bauer's avatar
Martin Bauer committed
20

Martin Bauer's avatar
Martin Bauer committed
21
22
23
    inner_or_boundary = True
    single_link = False

24
25
26
    def __init__(self, name=None):
        self._name = name

Frederik Hennig's avatar
Frederik Hennig committed
27
    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
28
29
        """
        This function defines the boundary behavior and must therefore be implemented by all boundaries.
Frederik Hennig's avatar
Frederik Hennig committed
30
31
        The boundary is defined through a list of sympy equations from which a boundary kernel is generated.

32
33
34
35
36
37
38
39
40
41
42
43
        Args:
        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
        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
        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.
        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.
        lb_method:      an instance of the LB method used. Use this to adapt the boundary to the method
                        (e.g. compressibility)
        index_field:    the boundary index field that can be used to retrieve and update boundary data
Frederik Hennig's avatar
Frederik Hennig committed
44

45
46
        Returns:
            list of pystencils assignments, or pystencils.AssignmentCollection
47
48
49
50
        """
        raise NotImplementedError("Boundary class has to overwrite __call__")

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

57
    @property
Martin Bauer's avatar
Martin Bauer committed
58
    def additional_data_init_callback(self):
59
60
        """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"""
61
62
        return None

Frederik Hennig's avatar
Frederik Hennig committed
63
    def get_additional_code_nodes(self, lb_method):
64
65
66
67
68
        """Return a list of code nodes that will be added in the generated code before the index field loop.

        Args:
            lb_method: lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`
        """
Frederik Hennig's avatar
Frederik Hennig committed
69
70
        return []

71
72
    @property
    def name(self):
73
74
75
76
77
78
        if self._name:
            return self._name
        else:
            return type(self).__name__

    @name.setter
Martin Bauer's avatar
Martin Bauer committed
79
80
    def name(self, new_value):
        self._name = new_value
81

82

Frederik Hennig's avatar
Frederik Hennig committed
83
84
# end class Boundary

85

Frederik Hennig's avatar
Frederik Hennig committed
86
87
class NoSlip(LbBoundary):
    """
Markus Holzer's avatar
Markus Holzer committed
88
89
    No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle.
    Extended for use with any streaming pattern.
90
91
92

    Args:
        name: optional name of the boundary.
Frederik Hennig's avatar
Frederik Hennig committed
93
94
    """

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

Frederik Hennig's avatar
Frederik Hennig committed
99
100
    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))
101
102

    def __hash__(self):
103
        return hash(self.name)
104
105

    def __eq__(self, other):
106
107
108
109
        if not isinstance(other, NoSlip):
            return False
        return self.name == other.name

110

Frederik Hennig's avatar
Frederik Hennig committed
111
# end class NoSlip
112

Frederik Hennig's avatar
Frederik Hennig committed
113
114

class UBB(LbBoundary):
115
116
117
118
119
120
121
122
123
124
125
126
    """Velocity bounce back boundary condition, enforcing specified velocity at obstacle

    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: adapts the velocity to the correct equilibrium when the lattice Boltzmann method holds
                                 a forcing term. If no forcing term is set and adapt_velocity_to_force is set to True
                                 it has no effect.
        dim: number of spatial dimensions
        name: optional name of the boundary.
    """
Martin Bauer's avatar
Martin Bauer committed
127

128
    def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'):
129
        super(UBB, self).__init__(name)
130
        self._velocity = velocity
Martin Bauer's avatar
Martin Bauer committed
131
        self._adaptVelocityToForce = adapt_velocity_to_force
132
133
134
135
136
        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
137
        self.data_type = data_type
Martin Bauer's avatar
Martin Bauer committed
138

139
    @property
Martin Bauer's avatar
Martin Bauer committed
140
    def additional_data(self):
141
142
        """ In case of the UBB boundary additional data is a velocity vector. This vector is added to each cell to
            realize velocity profiles for the inlet."""
Frederik Hennig's avatar
Frederik Hennig committed
143
        if self.velocity_is_callable:
144
            return [('vel_%d' % (i,), create_type(self.data_type)) for i in range(self.dim)]
145
146
147
148
        else:
            return []

    @property
Martin Bauer's avatar
Martin Bauer committed
149
    def additional_data_init_callback(self):
150
151
152
        """Initialise additional data of the boundary. For an example see
            `tutorial 02 <https://pycodegen.pages.i10git.cs.fau.de/lbmpy/notebooks/02_tutorial_boundary_setup.html>`_
            or lbmpy.geometry.add_pipe_inflow_boundary"""
153
154
155
        if callable(self._velocity):
            return self._velocity

Frederik Hennig's avatar
Frederik Hennig committed
156
    def get_additional_code_nodes(self, lb_method):
157
158
159
160
161
162
163
164
        """Return a list of code nodes that will be added in the generated code before the index field loop.

        Args:
            lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`

        Returns:
            list containing LbmWeightInfo and NeighbourOffsetArrays
        """
Frederik Hennig's avatar
Frederik Hennig committed
165
166
        return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)]

Frederik Hennig's avatar
Frederik Hennig committed
167
168
169
170
171
172
    @property
    def velocity_is_callable(self):
        """Returns True is velocity is callable. This means the velocity should be initialised via a callback function.
        This is useful if the inflow velocity should have a certain profile for instance"""
        return callable(self._velocity)

Frederik Hennig's avatar
Frederik Hennig committed
173
    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
Martin Bauer's avatar
Martin Bauer committed
174
        vel_from_idx_field = callable(self._velocity)
Frederik Hennig's avatar
Frederik Hennig committed
175
176
        vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity
        direction = dir_symbol
177

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

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

Frederik Hennig's avatar
Frederik Hennig committed
183
184
185
        velocity = tuple(v_i.get_shifted(*neighbor_offset)
                         if isinstance(v_i, Field.Access) and not vel_from_idx_field
                         else v_i
186
                         for v_i in vel)
187
188

        if self._adaptVelocityToForce:
Martin Bauer's avatar
Martin Bauer committed
189
190
191
            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]
192
193

        c_s_sq = sp.Rational(1, 3)
Martin Bauer's avatar
Martin Bauer committed
194
        weight_of_direction = LbmWeightInfo.weight_of_direction
195
196
        vel_term = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) * weight_of_direction(
            direction, lb_method)
197
198

        # Better alternative: in conserved value computation
Martin Bauer's avatar
Martin Bauer committed
199
        # rename what is currently called density to "virtual_density"
200
        # provide a new quantity density, which is constant in case of incompressible models
Martin Bauer's avatar
Martin Bauer committed
201
202
203
        if not lb_method.conserved_quantity_computation.zero_centered_pdfs:
            cqc = lb_method.conserved_quantity_computation
            density_symbol = sp.Symbol("rho")
Frederik Hennig's avatar
Frederik Hennig committed
204
            pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))]
Martin Bauer's avatar
Martin Bauer committed
205
206
207
            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
Frederik Hennig's avatar
Frederik Hennig committed
208
209
            result += [Assignment(f_in(inv_dir[direction]),
                                  f_out(direction) - vel_term * density_symbol)]
210
211
            return result
        else:
Frederik Hennig's avatar
Frederik Hennig committed
212
213
            return [Assignment(f_in(inv_dir[direction]),
                               f_out(direction) - vel_term)]
214
215


Frederik Hennig's avatar
Frederik Hennig committed
216
# end class UBB
217
218


219
220
class SimpleExtrapolationOutflow(LbBoundary):
    r"""
Markus Holzer's avatar
Markus Holzer committed
221
222
223
    Simple Outflow boundary condition :cite:`geier2015`, equation F.1 (listed below).
    This boundary condition extrapolates missing populations from the last layer of
    fluid cells onto the boundary by copying them in the normal direction.
224
225
226
227
228
229
230

    .. math ::
        f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yzt}


    Args:
        normal_direction: direction vector normal to the outflow
231
        stencil: stencil used by the lattice boltzmann method
232
233
234
235
236
237
238
239
240
241
242
243
244
        name: optional name of the boundary.
    """

    def __init__(self, normal_direction, stencil, name=None):
        if isinstance(normal_direction, str):
            normal_direction = direction_string_to_offset(normal_direction, dim=len(stencil[0]))

        if name is None:
            name = f"Simple Outflow: {offset_to_direction_string(normal_direction)}"

        self.normal_direction = normal_direction
        super(SimpleExtrapolationOutflow, self).__init__(name)

245
246
247
248
249
250
251
252
253
254
255
256
    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.

        Args:
            lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`

        Returns:
            list containing NeighbourOffsetArrays

        """
        return [NeighbourOffsetArrays(lb_method.stencil)]

257
    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
258
259
260
261
        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
        tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction))

        return Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](inv_dir[dir_symbol]))
262
263
264
265
266
# end class SimpleExtrapolationOutflow


class ExtrapolationOutflow(LbBoundary):
    r"""
Markus Holzer's avatar
Markus Holzer committed
267
268
269
270
271
    Outflow boundary condition :cite:`geier2015`, equation F.2, with u neglected (listed below).
    This boundary condition interpolates populations missing on the boundary in normal direction.
    For this interpolation, the PDF values of the last time step are used. They are interpolated
    between fluid cell and boundary cell. To get the PDF values from the last time step an index
    array is used which stores them.
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294

    .. math ::
        f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yz(t - \Delta t)} c \theta^{\frac{1}{2}}
        \frac{\Delta t}{\Delta x} + \left(1 - c \theta^{\frac{1}{2}} \frac{\Delta t}{\Delta x} \right)
         f_{\overline{1}jk(x - \Delta x)yzt}


    Args:
        normal_direction: direction vector normal to the outflow
        lb_method: the lattice boltzman method to be used in the simulation
        dt: lattice time step size
        dx: lattice spacing distance
        name: optional name of the boundary.
        streaming_pattern: Streaming pattern to be used in the simulation
        zeroth_timestep: for in-place patterns, whether the initial setup corresponds to an even or odd time step
        initial_density: floating point constant or callback taking spatial coordinates (x, y [,z]) as
                         positional arguments, specifying the initial density on boundary nodes
        initial_velocity: tuple of floating point constants or callback taking spatial coordinates (x, y [,z]) as
                          positional arguments, specifying the initial velocity on boundary nodes
    """

    def __init__(self, normal_direction, lb_method, dt=1, dx=1, name=None,
                 streaming_pattern='pull', zeroth_timestep=Timestep.BOTH,
295
296
297
                 initial_density=None, initial_velocity=None, data_type='double'):

        self.data_type = data_type
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338

        self.lb_method = lb_method
        self.stencil = lb_method.stencil
        self.dim = len(self.stencil[0])
        if isinstance(normal_direction, str):
            normal_direction = direction_string_to_offset(normal_direction, dim=self.dim)

        if name is None:
            name = f"Outflow: {offset_to_direction_string(normal_direction)}"

        self.normal_direction = normal_direction
        self.streaming_pattern = streaming_pattern
        self.zeroth_timestep = zeroth_timestep
        self.dx = sp.Number(dx)
        self.dt = sp.Number(dt)
        self.c = sp.sqrt(sp.Rational(1, 3)) * (self.dx / self.dt)

        self.initial_density = initial_density
        self.initial_velocity = initial_velocity
        self.equilibrium_calculation = None

        if initial_density and initial_velocity:
            equilibrium = lb_method.get_equilibrium(conserved_quantity_equations=AssignmentCollection([]))
            rho = lb_method.zeroth_order_equilibrium_moment_symbol
            u_vec = lb_method.first_order_equilibrium_moment_symbols
            eq_lambda = equilibrium.lambdify((rho,) + u_vec)
            post_pdf_symbols = lb_method.post_collision_pdf_symbols

            def calc_eq_pdfs(density, velocity, j):
                return eq_lambda(density, *velocity)[post_pdf_symbols[j]]

            self.equilibrium_calculation = calc_eq_pdfs

        super(ExtrapolationOutflow, self).__init__(name)

    def init_callback(self, boundary_data, **_):
        dim = boundary_data.dim
        coord_names = ['x', 'y', 'z'][:dim]
        pdf_acc = AccessPdfValues(self.stencil, streaming_pattern=self.streaming_pattern,
                                  timestep=self.zeroth_timestep, streaming_dir='out')

Frederik Hennig's avatar
Frederik Hennig committed
339
        def get_boundary_cell_pdfs(f_cell, b_cell, direction):
340
341
            if self.equilibrium_calculation is not None:
                density = self.initial_density(
Frederik Hennig's avatar
Frederik Hennig committed
342
                    *b_cell) if callable(self.initial_density) else self.initial_density
343
                velocity = self.initial_velocity(
Frederik Hennig's avatar
Frederik Hennig committed
344
345
                    *b_cell) if callable(self.initial_velocity) else self.initial_velocity
                return self.equilibrium_calculation(density, velocity, direction)
346
            else:
Frederik Hennig's avatar
Frederik Hennig committed
347
                return pdf_acc.read_pdf(boundary_data.pdf_array, f_cell, direction)
348
349

        for entry in boundary_data.index_array:
350
351
352
353
354
355
            center = tuple(entry[c] for c in coord_names)
            direction = self.stencil[entry["dir"]]
            inv_dir = self.stencil.index(inverse_direction(direction))
            tangential_offset = tuple(offset - normal for offset, normal in zip(direction, self.normal_direction))
            domain_cell = tuple(f + o for f, o in zip(center, tangential_offset))
            outflow_cell = tuple(f + o for f, o in zip(center, direction))
356
357

            #   Initial fluid cell PDF values
358
359
            entry['pdf'] = pdf_acc.read_pdf(boundary_data.pdf_array, domain_cell, inv_dir)
            entry['pdf_nd'] = get_boundary_cell_pdfs(domain_cell, outflow_cell, inv_dir)
360
361
362

    @property
    def additional_data(self):
363
364
        """Used internally only. For the ExtrapolationOutflow information of the previous PDF values is needed. This
        information is stored in the index vector."""
365
        data = [('pdf', create_type(self.data_type)), ('pdf_nd', create_type(self.data_type))]
366
367
368
369
        return data

    @property
    def additional_data_init_callback(self):
370
371
372
373
374
375
376
377
378
379
380
381
382
        return self.init_callback

    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.

        Args:
            lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`

        Returns:
            list containing NeighbourOffsetArrays

        """
        return [NeighbourOffsetArrays(lb_method.stencil)]
383
384
385
386
387
388

    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
        subexpressions = []
        boundary_assignments = []
        dtdx = sp.Rational(self.dt, self.dx)

389
390
        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
        tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction))
391

392
393
394
395
        interpolated_pdf_sym = sp.Symbol('pdf_inter')
        interpolated_pdf_asm = Assignment(interpolated_pdf_sym, (index_field[0]('pdf') * (self.c * dtdx))
                                          + ((sp.Number(1) - self.c * dtdx) * index_field[0]('pdf_nd')))
        subexpressions.append(interpolated_pdf_asm)
396

397
398
        asm = Assignment(f_in.center(inv_dir[dir_symbol]), interpolated_pdf_sym)
        boundary_assignments.append(asm)
399

400
401
        asm = Assignment(index_field[0]('pdf'), f_out[tangential_offset](inv_dir[dir_symbol]))
        boundary_assignments.append(asm)
402

403
404
        asm = Assignment(index_field[0]('pdf_nd'), interpolated_pdf_sym)
        boundary_assignments.append(asm)
405
406

        return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
407
408


409
410
411
# end class ExtrapolationOutflow


Frederik Hennig's avatar
Frederik Hennig committed
412
class FixedDensity(LbBoundary):
413
414
415
416
417
418
    """Boundary condition that fixes the density/pressure at the obstacle.

    Args:
        density: value of the density which should be set.
        name: optional name of the boundary.
    """
419

420
    def __init__(self, density, name=None):
421
422
        if name is None:
            name = "Fixed Density " + str(density)
423
        super(FixedDensity, self).__init__(name)
424
        self._density = density
425

Frederik Hennig's avatar
Frederik Hennig committed
426
    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
Martin Bauer's avatar
Martin Bauer committed
427
428
429
430
        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)
431

Martin Bauer's avatar
Martin Bauer committed
432
        cqc = lb_method.conserved_quantity_computation
Martin Bauer's avatar
Martin Bauer committed
433
        velocity = cqc.defined_symbols()['velocity']
Martin Bauer's avatar
Martin Bauer committed
434
435
        symmetric_eq = remove_asymmetric_part_of_main_assignments(lb_method.get_equilibrium(),
                                                                  degrees_of_freedom=velocity)
Frederik Hennig's avatar
Frederik Hennig committed
436
        substitutions = {sym: f_out(i) for i, sym in enumerate(lb_method.pre_collision_pdf_symbols)}
Martin Bauer's avatar
Martin Bauer committed
437
        symmetric_eq = symmetric_eq.new_with_substitutions(substitutions)
438

Martin Bauer's avatar
Martin Bauer committed
439
440
        simplification = create_simplification_strategy(lb_method)
        symmetric_eq = simplification(symmetric_eq)
441

Martin Bauer's avatar
Martin Bauer committed
442
        density_symbol = cqc.defined_symbols()['density']
443
444

        density = self._density
Martin Bauer's avatar
Martin Bauer committed
445
446
        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
447
448
449
        density_eq = equilibrium_input.main_assignments[0]
        assert density_eq.lhs == density_symbol
        transformed_density = density_eq.rhs
450

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

Martin Bauer's avatar
Martin Bauer committed
455
456
        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
457

Frederik Hennig's avatar
Frederik Hennig committed
458
459
        return subexpressions + [Assignment(f_in(inv_dir[dir_symbol]),
                                            2 * eq_component - f_out(dir_symbol))]
460
461


Frederik Hennig's avatar
Frederik Hennig committed
462
# end class FixedDensity
Martin Bauer's avatar
Martin Bauer committed
463

464

Frederik Hennig's avatar
Frederik Hennig committed
465
class NeumannByCopy(LbBoundary):
466
467
    """Neumann boundary condition which is implemented by coping the PDF values to achieve similar values at the fluid
       and the boundary node"""
Frederik Hennig's avatar
Frederik Hennig committed
468
469

    def get_additional_code_nodes(self, lb_method):
470
471
472
473
474
475
476
477
        """Return a list of code nodes that will be added in the generated code before the index field loop.

        Args:
            lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`

        Returns:
            list containing NeighbourOffsetArrays
        """
Frederik Hennig's avatar
Frederik Hennig committed
478
479
480
481
482
483
        return [NeighbourOffsetArrays(lb_method.stencil)]

    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
        neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
        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))]
484
485
486
487
488
489
490

    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
491
492


Frederik Hennig's avatar
Frederik Hennig committed
493
# end class NeumannByCopy
Martin Bauer's avatar
Martin Bauer committed
494

Martin Bauer's avatar
Martin Bauer committed
495

Frederik Hennig's avatar
Frederik Hennig committed
496
class StreamInConstant(LbBoundary):
497
498
499
500
501
502
503
    """Boundary condition that takes a constant and overrides the boundary PDFs with this value. This is used for
    debugging mainly.

    Args:
        constant: value which should be set for the PDFs at the boundary cell.
        name: optional name of the boundary.
    """
504

505
506
507
    def __init__(self, constant, name=None):
        super(StreamInConstant, self).__init__(name)
        self._constant = constant
Martin Bauer's avatar
Martin Bauer committed
508

Frederik Hennig's avatar
Frederik Hennig committed
509
    def get_additional_code_nodes(self, lb_method):
510
511
512
513
514
515
516
517
        """Return a list of code nodes that will be added in the generated code before the index field loop.

        Args:
            lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`

        Returns:
            list containing NeighbourOffsetArrays
        """
Frederik Hennig's avatar
Frederik Hennig committed
518
519
520
521
522
523
        return [NeighbourOffsetArrays(lb_method.stencil)]

    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
        neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
        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
524
525

    def __hash__(self):
526
527
        # All boundaries of these class behave equal -> should also be equal
        return hash("StreamInConstant")
Martin Bauer's avatar
Martin Bauer committed
528
529

    def __eq__(self, other):
530
        return type(other) == StreamInConstant
Frederik Hennig's avatar
Frederik Hennig committed
531
# end class StreamInConstant