Commit b5ac8c01 authored by Frederik Hennig's avatar Frederik Hennig
Browse files

Integrated exotic boundaries

parent dd1cbd09
%% Cell type:code id: tags:
```
%load_ext autoreload
%autoreload 2
```
%% Cell type:code id: tags:
```
import sympy as sp
import numpy as np
import pystencils as ps
from lbmpy.stencils import get_stencil
from lbmpy.boundaries import NoSlip
from lbmpy.boundaries.boundaryhandling import create_lattice_boltzmann_boundary_kernel
from lbmpy.creationfunctions import create_lb_method
stencil = get_stencil('D2Q9')
pdf_field = ps.fields('pdf_field(9): [2D]')
from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
from pystencils.data_types import TypedSymbol, create_type
from pystencils.field import FieldType
lb_method = create_lb_method(stencil=stencil)
noslip = NoSlip()
index_struct_dtype = numpy_data_type_for_boundary_object(noslip, lb_method.dim)
index_field = ps.Field('indexVector', FieldType.INDEXED, index_struct_dtype, layout=[0],
shape=(TypedSymbol("indexVectorSize", create_type(np.int64)), 1), strides=(1, 1))
```
%% Cell type:markdown id: tags:
## Defining Boundaries: Current Implementaiton
The boundary cell is accessed directly from the fluid cell through the `neighbour` offset with the inverse direction `inverse_dir`. Both of those are symbolic and are replaced by array accesses during code generation.
```
class NoSlip:
def __call__(self, pdf_field, direction_symbol, lb_method, **kwargs):
neighbor = BoundaryOffsetInfo.offset_from_dir(direction_symbol, lb_method.dim)
inverse_dir = BoundaryOffsetInfo.inv_dir(direction_symbol)
return [Assignment(pdf_field[neighbor](inverse_dir), pdf_field(direction_symbol))]
```
%% Cell type:code id: tags:
```
dir_symbol = TypedSymbol('dir', create_type(np.int64))
noslip(pdf_field, dir_symbol, lb_method)[0]
```
%%%% Output: execute_result
$\displaystyle {{pdf_field}_{({cx}_{dir},{cy}_{dir})}^{invdir[dir]}} \leftarrow {{pdf_field}_{(0,0)}^{dir}}$
Assignment(pdf_field_dd4358bcc3a3, pdf_field_8ffb2303630e)
%% Cell type:code id: tags:
```
noslip_ast = create_lattice_boltzmann_boundary_kernel(pdf_field, index_field, lb_method, noslip)
ps.show_code(noslip_ast)
```
%%%% Output: display_data
%%%% Output: display_data
%% Cell type:markdown id: tags:
## Defining Boundaries: Flexible Approach
The arrangment of the populations in memory is abstracted away through two proxy fields
$$
f_{out}, \quad f_{in}
$$
for the populations streaming *out* of the current cell and *into* the current cell. With those, boundaries can be defined in general notation, without knowledge about the streaming pattern. For example, the NoSlip boundary is defined as:
%% Cell type:code id: tags:
```
class FlexibleNoSlip:
def __call__(self, f_out, f_in, dir_symbol, inv_dir):
return ps.Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol))
```
%% Cell type:code id: tags:
```
f_out, f_in = ps.fields('f_out(9), f_in(9): [2D]')
dir_symbol = TypedSymbol('dir', create_type(np.int64))
inv_dir = sp.IndexedBase('invdir')
flex_noslip = FlexibleNoSlip()
noslip_assignments = flex_noslip(f_out, f_in, dir_symbol, inv_dir)
noslip_assignments
```
%%%% Output: execute_result
$\displaystyle {{f_in}_{(0,0)}^{invdir[dir]}} \leftarrow {{f_out}_{(0,0)}^{dir}}$
Assignment(f_in_1ba4ca851771, f_out_8ffb2303630e)
%% Cell type:markdown id: tags:
During Code Generation, all accesses to $f_{out}$ and $f_{in}$ are replaced by accesses to `pdf_field` according to the streaming pattern.
### Example: Boundary handling after an even time step with the AA-Pattern
For the current cell, the populations flowing out are the ones written in the previous step, and the populations flowing in are those begin read in the next step. The previous step was controlled by the `AAEvenTimeStepAccessor`, and the next step is controlled by `AAOddTimeStepAccessor`. The required memory locations are thus given by:
- Outbound populations: `AAEvenTimeStepAccessor.write`
- Inbound populations: `AAOddTimeStepAccessor.read`
%% Cell type:code id: tags:
```
from lbmpy.fieldaccess import AAEvenTimeStepAccessor, AAOddTimeStepAccessor
def substitute_proxies(assignments):
# Get the accessor lists for the streaming pattern (here AA)
outward_accesses = AAEvenTimeStepAccessor.write(pdf_field, stencil)
inward_accesses = AAOddTimeStepAccessor.read(pdf_field, stencil)
# Define the symbols which will correspond to the index translation arrays
# in the generated code
f_out_dir_idx = sp.IndexedBase('f_out_dir_idx')
f_out_offsets_x = sp.IndexedBase('f_out_offset_x')
f_out_offsets_y = sp.IndexedBase('f_out_offset_y')
f_in_inv_dir_idx = sp.IndexedBase('f_in_inv_dir_idx')
f_in_inv_offsets_x = sp.IndexedBase('f_in_inv_offsets_x')
f_in_inv_offsets_y = sp.IndexedBase('f_in_inv_offsets_y')
# Substitute all proxy field accesses
if not isinstance(assignments, ps.AssignmentCollection):
assignments = ps.AssignmentCollection([assignments])
accessor_subs = dict()
for fa in assignments.atoms(ps.Field.Access):
if fa.field == f_out:
if fa.offsets == (0,0):
if isinstance(fa.index[0], int):
accessor_subs[fa] = outward_accesses[fa.index[0]]
elif fa.index[0] == dir_symbol:
accessor_subs[fa] = pdf_field[
f_out_offsets_x[dir_symbol], f_out_offsets_y[dir_symbol]
]( f_out_dir_idx[dir_symbol] )
else:
# other cases like inverse direction, etc.
pass
else:
# non-trivial neighbour accesses -> add neighbour offset to streaming pattern offsets
pass
elif fa.field == f_in:
if fa.offsets == (0,0):
if isinstance(fa.index[0], int):
accessor_subs[fa] = inward_accesses[fa.index[0]]
elif fa.index[0] == inv_dir[dir_symbol]:
accessor_subs[fa] = pdf_field[
f_in_inv_offsets_x[dir_symbol], f_in_inv_offsets_y[dir_symbol]
]( f_in_inv_dir_idx[dir_symbol] )
else:
# other cases
pass
else:
# non-trivial neighbour accesses -> add neighbour offset to streaming pattern offsets
pass
else:
pass
return assignments.new_with_substitutions(accessor_subs)
substitute_proxies(noslip_assignments).main_assignments[0]
```
%%%% Output: execute_result
$\displaystyle {{pdf_field}_{({f_{in inv offsets x}}_{dir},{f_{in inv offsets y}}_{dir})}^{f_in_inv_dir_idx[dir]}} \leftarrow {{pdf_field}_{({f_{out offset x}}_{dir},{f_{out offset y}}_{dir})}^{f_out_dir_idx[dir]}}$
Assignment(pdf_field_5ed29f4dfe26, pdf_field_9520cd7bd10d)
%% Cell type:markdown id: tags:
Accesses with fixed indices are replaced directly:
%% Cell type:code id: tags:
```
as2 = ps.Assignment(f_in(3), f_out(4) + f_out(5) + f_out(dir_symbol))
as2
```
%%%% Output: execute_result
$\displaystyle {{f_in}_{(0,0)}^{3}} \leftarrow {{f_out}_{(0,0)}^{dir}} + {{f_out}_{(0,0)}^{4}} + {{f_out}_{(0,0)}^{5}}$
Assignment(f_in_C^3, f_out_8ffb2303630e + f_out_C^4 + f_out_C^5)
%% Cell type:code id: tags:
```
substitute_proxies(as2).main_assignments[0]
```
%%%% Output: execute_result
$\displaystyle {{pdf_field}_{(1,0)}^{4}} \leftarrow {{pdf_field}_{({f_{out offset x}}_{dir},{f_{out offset y}}_{dir})}^{f_out_dir_idx[dir]}} + {{pdf_field}_{(0,0)}^{3}} + {{pdf_field}_{(0,0)}^{8}}$
Assignment(pdf_field_E^4, pdf_field_9520cd7bd10d + pdf_field_C^3 + pdf_field_C^8)
%% Cell type:markdown id: tags:
The index translation arrays are prepended in the generated code:
```
const int64_t f_out_dir_idx [] = { 0, 2, 1, 4, 3, 8, 7, 6, 5 };
const int64_t f_out_offsets_x [] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
const int64_t f_out_offsets_y [] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
const int64_t f_in_inv_dir_idx [] = { 0, 1, 2, 3, 4, 5, 6, 7, 8 };
const int64_t f_in_inv_offsets_x [] = { 0, 0, 0, -1, 1, -1, 1, -1, 1 };
const int64_t f_in_inv_offsets_y [] = { 0, 1, -1, 0, 0, 1, 1, -1, -1 };
```
The generated boundary kernel is of course only fit for application after an even time step. For the AA-pattern's odd time steps, a second kernel needs to be generated. The definition of the NoSlip boundary itself is the same; only the index translation arrays are different.
## Integration with the pystencils Boundary Handling Infrastructure
The existing `pystencils.boundaries.BoundaryHandling` class is extended by a new subclass `FlexibleLBMBoundaryHandling` analogously to the implementation of the existing `LatticeBoltzmannBoundaryHandling`. Its behaviour depends on the streaming pattern which is passed to its constructor. It overrides `__call__` and `_create_boundary_kernel`:
```
class FlexibleLBMBoundaryHandling(BoundaryHandling):
def __init__(self, kernel_type): # kernel_type in ['pull', 'push', 'aa', 'esotwist']
(...)
def __call__(self, **kwargs):
# Handle boundaries like in base class, but call different kernels
# depending on the time step modulus
def _create_boundary_kernel(self, **kwargs):
if self.kernel_type in ['aa', 'esotwist']:
# in-place methods: create two boundary kernels, for even and odd timesteps
else:
# two-fields methods: create only one boundary kernel
```
## Integration with waLBerla
The waLBerla boundary code generation can also be extended to generate implementations which track the time step internally and selectively call one of two implementations:
```
class BoundarySweep{
public:
void operator() (IBlock * block){
if(even_timestep) evenSweep(block);
else oddSweep(block);
}
}
```
......@@ -9,9 +9,11 @@ from lbmpy.advanced_streaming.utility import get_accessor, inverse_dir_index, is
from itertools import product
def _array_pattern(dtype, name, content):
return f"const {str(dtype)} {name} [] = {{ {','.join(str(c) for c in content)} }}; \n"
class BetweenTimestepsIndexing:
# ==============================================
......@@ -187,12 +189,13 @@ class BetweenTimestepsIndexing:
# end class AdvancedStreamingIndexing
class NeighbourOffsetArraysForStencil(CustomCodeNode):
class NeighbourOffsetArrays(CustomCodeNode):
@staticmethod
def symbolic_neighbour_offset_from_dir(dir_idx, dim):
def symbolic_neighbour_offset(dir_idx, dim):
return tuple([sp.IndexedBase(symbol, shape=(1,))[dir_idx]
for symbol in NeighbourOffsetArraysForStencil._offset_symbols(dim)])
for symbol in NeighbourOffsetArrays._offset_symbols(dim)])
@staticmethod
def _offset_symbols(dim):
......@@ -202,11 +205,11 @@ class NeighbourOffsetArraysForStencil(CustomCodeNode):
offsets_dtype = create_type(offsets_dtype)
dim = len(stencil[0])
array_symbols = NeighbourOffsetArraysForStencil._offset_symbols(dim)
array_symbols = NeighbourOffsetArrays._offset_symbols(dim)
code = "\n"
for i, arrsymb in enumerate(array_symbols):
code += _array_pattern(offsets_dtype, arrsymb.name, (d[i] for d in stencil))
offset_symbols = NeighbourOffsetArraysForStencil._offset_symbols(dim)
super(NeighbourOffsetArraysForStencil, self).__init__(code, symbols_read=set(),
symbols_defined=set(offset_symbols))
\ No newline at end of file
offset_symbols = NeighbourOffsetArrays._offset_symbols(dim)
super(NeighbourOffsetArrays, self).__init__(code, symbols_read=set(),
symbols_defined=set(offset_symbols))
......@@ -114,5 +114,3 @@ class AccessPdfValues:
return tuple(idx[v] for v in ('x', 'y', 'z')[:len(idx) - 1] + ('dir',))
return self.read_multiple(pdf_arr, (to_index_tuple(idx) for idx in index_list))
......@@ -4,7 +4,7 @@ from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
from pystencils.data_types import create_type
from pystencils.sympyextensions import get_symmetric_part
from lbmpy.simplificationfactory import create_simplification_strategy
from lbmpy.advanced_streaming.indexing import NeighbourOffsetArraysForStencil
from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays
class Boundary:
......@@ -127,18 +127,17 @@ class UBB(Boundary):
return self._velocity
def get_additional_code_nodes(self, lb_method):
return [LbmWeightInfo(lb_method), NeighbourOffsetArraysForStencil(lb_method.stencil)]
return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
vel_from_idx_field = callable(self._velocity)
vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity
direction = dir_symbol
assert self.dim == lb_method.dim, f"Dimension of UBB ({self.dim}) does not match dimension of method ({lb_method.dim})"
assert self.dim == lb_method.dim, \
f"Dimension of UBB ({self.dim}) does not match dimension of method ({lb_method.dim})"
neighbor_offset = NeighbourOffsetArraysForStencil.symbolic_neighbour_offset_from_dir(
direction, lb_method.dim)
neighbor_offset = NeighbourOffsetArrays.symbolic_neighbour_offset(direction, lb_method.dim)
velocity = tuple(v_i.get_shifted(*neighbor_offset)
if isinstance(v_i, Field.Access) and not vel_from_idx_field
......@@ -153,8 +152,8 @@ class UBB(Boundary):
c_s_sq = sp.Rational(1, 3)
weight_of_direction = LbmWeightInfo.weight_of_direction
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)
* sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) \
* weight_of_direction(direction)
# Better alternative: in conserved value computation
# rename what is currently called density to "virtual_density"
......@@ -174,17 +173,16 @@ class UBB(Boundary):
f_out(direction) - vel_term)]
# end class UBB
# TODO
class FixedDensity(Boundary):
def __init__(self, density, name=None):
raise NotImplementedError() # not yet operable
if name is None:
name = "Fixed Density " + str(density)
super(FixedDensity, self).__init__(name)
self._density = density
def __call__(self, pdf_field, direction_symbol, lb_method, **kwargs):
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
"""Boundary condition that fixes the density/pressure at the obstacle"""
def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom):
......@@ -192,14 +190,11 @@ class FixedDensity(Boundary):
for a in assignment_collection.main_assignments]
return assignment_collection.copy(new_main_assignments)
neighbor = BoundaryOffsetInfo.offset_from_dir(direction_symbol, lb_method.dim)
inverse_dir = BoundaryOffsetInfo.inv_dir(direction_symbol)
cqc = lb_method.conserved_quantity_computation
velocity = cqc.defined_symbols()['velocity']
symmetric_eq = remove_asymmetric_part_of_main_assignments(lb_method.get_equilibrium(),
degrees_of_freedom=velocity)
substitutions = {sym: pdf_field(i) for i, sym in enumerate(lb_method.pre_collision_pdf_symbols)}
substitutions = {sym: f_out(i) for i, sym in enumerate(lb_method.pre_collision_pdf_symbols)}
symmetric_eq = symmetric_eq.new_with_substitutions(substitutions)
simplification = create_simplification_strategy(lb_method)
......@@ -214,24 +209,27 @@ class FixedDensity(Boundary):
assert density_eq.lhs == density_symbol
transformed_density = density_eq.rhs
conditions = [(eq_i.rhs, sp.Equality(direction_symbol, i))
conditions = [(eq_i.rhs, sp.Equality(dir_symbol, i))
for i, eq_i in enumerate(symmetric_eq.main_assignments)] + [(0, True)]
eq_component = sp.Piecewise(*conditions)
subexpressions = [Assignment(eq.lhs, transformed_density if eq.lhs == density_symbol else eq.rhs)
for eq in symmetric_eq.subexpressions]
return subexpressions + [Assignment(pdf_field[neighbor](inverse_dir),
2 * eq_component - pdf_field(direction_symbol))]
return subexpressions + [Assignment(f_in(inv_dir[dir_symbol]),
2 * eq_component - f_out(dir_symbol))]
# end class FixedDensity
class NeumannByCopy(Boundary):
def __call__(self, pdf_field, direction_symbol, lb_method, **kwargs):
raise NotImplementedError() # not yet operable
neighbor = BoundaryOffsetInfo.offset_from_dir(direction_symbol, lb_method.dim)
inverse_dir = BoundaryOffsetInfo.inv_dir(direction_symbol)
return [Assignment(pdf_field[neighbor](inverse_dir), pdf_field(inverse_dir)),
Assignment(pdf_field[neighbor](direction_symbol), pdf_field(direction_symbol))]
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):
neighbour_offset = NeighbourOffsetArrays.symbolic_neighbour_offset(dir_symbol, lb_method.dim)
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))]
def __hash__(self):
# All boundaries of these class behave equal -> should also be equal
......@@ -239,19 +237,21 @@ class NeumannByCopy(Boundary):
def __eq__(self, other):
return type(other) == NeumannByCopy
# end class NeumannByCopy
class StreamInConstant(Boundary):
def __init__(self, constant, name=None):
raise NotImplementedError() # not yet operable
super(StreamInConstant, self).__init__(name)
self._constant = constant
def __call__(self, pdf_field, direction_symbol, lb_method, **kwargs):
neighbor = BoundaryOffsetInfo.offset_from_dir(direction_symbol, lb_method.dim)
inverse_dir = BoundaryOffsetInfo.inv_dir(direction_symbol)
return [Assignment(pdf_field[neighbor](inverse_dir), self._constant),
Assignment(pdf_field[neighbor](direction_symbol), self._constant)]
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):
neighbour_offset = NeighbourOffsetArrays.symbolic_neighbour_offset(dir_symbol, lb_method.dim)
return [Assignment(f_in(inv_dir[dir_symbol]), self._constant),
Assignment(f_out[neighbour_offset](dir_symbol), self._constant)]
def __hash__(self):
# All boundaries of these class behave equal -> should also be equal
......@@ -259,3 +259,4 @@ class StreamInConstant(Boundary):
def __eq__(self, other):
return type(other) == StreamInConstant
# end class StreamInConstant
......@@ -23,7 +23,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
self._inplace = is_inplace(streaming_pattern)
self._prev_timestep = None
super(LatticeBoltzmannBoundaryHandling, self).__init__(data_handling, pdf_field_name, lb_method.stencil,
name, flag_interface, target, openmp)
name, flag_interface, target, openmp)
# ------------------------- Overridden methods of pystencils.BoundaryHandling -------------------------
......@@ -121,7 +121,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
method = self._lb_method
stencil = np.array(method.stencil)
inv_direction = np.array([method.stencil.index(inverse_direction(d))
for d in method.stencil])
for d in method.stencil])
result = np.zeros(self.dim)
for b in dh.iterate(ghost_layers=ff_ghost_layers):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment