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

Flexible NoSlip Boundary

parent 235979c0
Pipeline #26778 failed with stage
in 16 minutes and 4 seconds
from .flexible_boundaries import AdvancedStreamingBoundaryOffsetInfo, \
FlexibleNoSlip, create_flexible_lbm_boundary_kernel
__all__ = [ 'AdvancedStreamingBoundaryOffsetInfo',
'FlexibleNoSlip',
'create_flexible_lbm_boundary_kernel' ]
\ No newline at end of file
import sympy as sp
import numpy as np
from lbmpy.boundaries.boundaryconditions import Boundary
from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor, \
StreamPushTwoFieldsAccessor, \
AAEvenTimeStepAccessor, \
AAOddTimeStepAccessor, \
EsoTwistEvenTimeStepAccessor, \
EsoTwistOddTimeStepAccessor
from pystencils import Assignment
from pystencils.data_types import TypedSymbol, create_type
from pystencils.backends.cbackend import CustomCodeNode
from pystencils import create_indexed_kernel
class AdvancedStreamingBoundaryOffsetInfo(CustomCodeNode):
# --------------------------- Functions to be used by boundaries --------------------------
@staticmethod
def outward_offset_from_dir(dir_idx, dim):
return tuple([sp.IndexedBase(symbol, shape=(1,))[dir_idx]
for symbol in AdvancedStreamingBoundaryOffsetInfo._outward_offset_symbols(dim)])
@staticmethod
def inverse_inward_offset_from_dir(dir_idx, dim):
return tuple([sp.IndexedBase(symbol, shape=(1,))[dir_idx]
for symbol in AdvancedStreamingBoundaryOffsetInfo._inverse_inward_offset_symbols(dim)])
@staticmethod
def outward_index_from_dir(dir_idx):
return sp.IndexedBase(AdvancedStreamingBoundaryOffsetInfo.OUTWARD_INDEX_SYMBOL, shape=(1,))[dir_idx]
@staticmethod
def inverse_inward_index_from_dir(dir_idx):
return sp.IndexedBase(AdvancedStreamingBoundaryOffsetInfo.INV_INWARD_INDEX_SYMBOL, shape=(1,))[dir_idx]
# ---------------------------------- Internal ---------------------------------------------
def __init__(self, pdf_field, stencil, between_timesteps='both', kernel_type='pull'):
if between_timesteps not in ['both', 'odd_to_even', 'even_to_odd']:
raise ValueError(
"Invalid value of parameter 'between_timesteps'. Must be one of 'both', 'odd_to_even' or 'even_to_odd'.",
between_timesteps)
if between_timesteps == 'both' and kernel_type in ['aa', 'esotwist']:
raise ValueError('Cannot create an offset info for both kinds of timesteps for kernel type ' + kernel_type)
odd_to_even = (between_timesteps == 'odd_to_even')
even_accessors = {
'pull': StreamPullTwoFieldsAccessor,
'push': StreamPushTwoFieldsAccessor,
'aa': AAEvenTimeStepAccessor,
'esotwist': EsoTwistEvenTimeStepAccessor
}
odd_accessors = {
'pull': StreamPullTwoFieldsAccessor,
'push': StreamPushTwoFieldsAccessor,
'aa': AAOddTimeStepAccessor,
'esotwist': EsoTwistOddTimeStepAccessor
}
try:
even_accessor = even_accessors[kernel_type]
odd_accessor = odd_accessors[kernel_type]
except KeyError:
raise ValueError(
"Invalid value of parameter 'kernel_type'.", kernel_type)
if between_timesteps == 'both':
assert even_accessor == odd_accessor
write_field_accesses = (
odd_accessor if odd_to_even else even_accessor).write(pdf_field, stencil)
read_field_accesses = (
even_accessor if odd_to_even else odd_accessor).read(pdf_field, stencil)
dim = len(stencil[0])
code = "\n"
# Offsets and Indices for the stream-out-populations
outward_offset_sym = self._outward_offset_symbols(dim)
for i in range(dim):
offset_str = ", ".join([str(a.offsets[i])
for a in write_field_accesses])
code += "const int64_t %s [] = { %s };\n" % (
outward_offset_sym[i].name, offset_str)
out_acc_indices = [str(a.index[0]) for a in write_field_accesses]
code += "const int %s [] = { %s };\n" % (
self.OUTWARD_INDEX_SYMBOL.name, ", ".join(out_acc_indices))
# Offsets and Indices for the stream-in-populations
# Permute the read_field_accesses for the inverse directions
inverse_indices = [self._inverse_dir_index(stencil, i)
for i in range(len(stencil))]
inv_read_field_accesses = [read_field_accesses[idx]
for idx in inverse_indices]
inv_inward_offset_sym = self._inverse_inward_offset_symbols(dim)
for i in range(dim):
offset_str = ", ".join([str(a.offsets[i])
for a in inv_read_field_accesses])
code += "const int64_t %s [] = { %s };\n" % (
inv_inward_offset_sym[i].name, offset_str)
inv_in_acc_indices = [str(a.index[0]) for a in inv_read_field_accesses]
code += "const int %s [] = { %s };\n" % (
self.INV_INWARD_INDEX_SYMBOL.name, ", ".join(inv_in_acc_indices))
defined_symbols = set(outward_offset_sym + inv_inward_offset_sym +
[self.OUTWARD_INDEX_SYMBOL, self.INV_INWARD_INDEX_SYMBOL])
super(AdvancedStreamingBoundaryOffsetInfo, self).__init__(
code, symbols_read=set(), symbols_defined=defined_symbols)
@staticmethod
def _inverse_dir_index(stencil, dir):
return stencil.index(tuple(-d for d in stencil[dir]))
@staticmethod
def _outward_offset_symbols(dim):
return [TypedSymbol(f"OutOffset_{d}", create_type(np.int64)) for d in ['x', 'y', 'z'][:dim]]
@staticmethod
def _inverse_inward_offset_symbols(dim):
return [TypedSymbol(f"InvInOffset_{d}", create_type(np.int64)) for d in ['x', 'y', 'z'][:dim]]
OUTWARD_INDEX_SYMBOL = TypedSymbol("outIdx", int)
INV_INWARD_INDEX_SYMBOL = TypedSymbol("invInIdx", int)
# end class AdvancedStreamingBoundaryOffsetInfo
class FlexibleNoSlip(Boundary):
def __init__(self, name=None):
"""Set an optional name here, to mark boundaries, for example for force evaluations"""
super(FlexibleNoSlip, self).__init__(name)
"""
No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle.
Extended for use with any streaming pattern.
"""
def __call__(self, pdf_field, direction_symbol, lb_method, **kwargs):
outward_offset = AdvancedStreamingBoundaryOffsetInfo.outward_offset_from_dir(direction_symbol, lb_method.dim)
outward_idx = AdvancedStreamingBoundaryOffsetInfo.outward_index_from_dir(direction_symbol)
inv_inward_offset = AdvancedStreamingBoundaryOffsetInfo.inverse_inward_offset_from_dir(direction_symbol, lb_method.dim)
inv_inward_idx = AdvancedStreamingBoundaryOffsetInfo.inverse_inward_index_from_dir(direction_symbol)
return [Assignment(pdf_field[inv_inward_offset](inv_inward_idx), pdf_field[outward_offset](outward_idx))]
def __hash__(self):
# All boundaries of these class behave equal -> should also be equal (as long as name is equal)
return hash(self.name)
def __eq__(self, other):
if not isinstance(other, FlexibleNoSlip):
return False
return self.name == other.name
# end class FlexibleNoSlip
def create_flexible_lbm_boundary_kernel(pdf_field, index_field, lb_method, boundary_functor,
between_timesteps='both', kernel_type='pull',
target='cpu', openmp=True):
elements = [AdvancedStreamingBoundaryOffsetInfo(pdf_field, lb_method.stencil, between_timesteps, kernel_type),
LbmWeightInfo(lb_method)]
index_arr_dtype = index_field.dtype.numpy_dtype
dir_symbol = TypedSymbol("dir", index_arr_dtype.fields['dir'][0])
elements += [Assignment(dir_symbol, index_field[0]('dir'))]
elements += boundary_functor(pdf_field=pdf_field, direction_symbol=dir_symbol,
lb_method=lb_method, index_field=index_field)
return create_indexed_kernel(elements, [index_field], target=target, cpu_openmp=openmp)
import numpy as np
import sympy as sp
from lbmpy.fieldaccess import *
from lbmpy.advanced_streaming import FlexibleNoSlip, create_flexible_lbm_boundary_kernel
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
\ No newline at end of file
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