Commit 40bd6745 authored by Frederik Hennig's avatar Frederik Hennig
Browse files

FlexibleNoSlip with test case

parent ab28a143
Pipeline #26901 passed with stage
in 16 minutes and 6 seconds
from .boundaryindexing import AdvancedStreamingBoundaryIndexing
from .boundaryconditions import FlexibleBoundary, FlexibleNoSlip
from .boundaryhandling import create_flexible_lbm_boundary_kernel
__all__ = ['AdvancedStreamingBoundaryIndexing', 'FlexibleBoundary', 'FlexibleNoSlip']
__all__ = ['AdvancedStreamingBoundaryIndexing', 'FlexibleBoundary',
'FlexibleNoSlip', 'create_flexible_lbm_boundary_kernel']
......@@ -17,8 +17,11 @@ def create_flexible_lbm_boundary_kernel(pdf_field, index_field, lb_method, bound
boundary_assignments = boundary_functor(f_out, f_in, dir_symbol, inv_dir, lb_method, index_field)
boundary_assignments = indexing.substitute_proxies(boundary_assignments)
elements = [indexing.code_node]
elements = []
# elements = [indexing.code_node]
# elements += boundary_functor.additional_code_nodes
elements += [Assignment(dir_symbol, index_field[0]('dir'))]
elements += boundary_functor.additional_code_nodes
elements += boundary_assignments.all_assignments
return create_indexed_kernel(elements, [index_field], target=target, cpu_openmp=openmp)
kernel = create_indexed_kernel(elements, [index_field], target=target, cpu_openmp=openmp)
kernel.body.insert_front(indexing.code_node)
return kernel
......@@ -93,30 +93,30 @@ class AdvancedStreamingBoundaryIndexing:
self._required_arrays = set()
self._trivial_translations = self._collect_trivial_translations()
def _index_array_name(self, f_dir, inverse):
def _index_array_symbol(self, f_dir, inverse):
assert f_dir in ['in', 'out']
inv = '_inv' if inverse else ''
name = f"f_{f_dir}{inv}_dir_idx"
return name
return TypedSymbol(name, self._index_dtype)
def _offset_array_names(self, f_dir, inverse):
def _offset_array_symbols(self, f_dir, inverse):
assert f_dir in ['in', 'out']
inv = '_inv' if inverse else ''
name_base = f"f_{f_dir}{inv}_offsets_"
names = [name_base + d for d in self._coordinate_names]
return names
symbols = [TypedSymbol(name_base + d, self._index_dtype) for d in self._coordinate_names]
return symbols
def _array_symbols(self, f_dir, inverse, index):
if (f_dir, inverse) in self._trivial_translations:
return {'index': index, 'offsets': (0, ) * self._dim}
index_array_name = self._index_array_name(f_dir, inverse)
offset_array_names = self._offset_array_names(f_dir, inverse)
index_array_symbol = self._index_array_symbol(f_dir, inverse)
offset_array_symbols = self._offset_array_symbols(f_dir, inverse)
self._required_arrays.add((f_dir, inverse))
return {
'index': sp.IndexedBase(TypedSymbol(index_array_name, self._index_dtype), shape=(1,))[index],
'offsets': tuple(sp.IndexedBase(TypedSymbol(n, self._index_dtype), shape=(1,))[index]
for n in offset_array_names)
'index': sp.IndexedBase(index_array_symbol, shape=(1,))[index],
'offsets': tuple(sp.IndexedBase(s, shape=(1,))[index]
for s in offset_array_symbols)
}
# =================================
......@@ -198,22 +198,33 @@ class AdvancedStreamingBoundaryIndexing:
class TranslationArraysNode(CustomCodeNode):
def __init__(self, indexing):
code = '\n'
code = ''
symbols_defined = set()
for f_dir, inv in indexing._required_arrays:
indices, offsets = indexing._get_translated_indices_and_offsets(f_dir, inv)
index_array_symbol = indexing._index_array_symbol(f_dir, inv)
symbols_defined.add(index_array_symbol)
acc_indices = ", ".join([str(i) for i in indices])
code += self._array_pattern(indexing._index_dtype, indexing._index_array_name(f_dir, inv), acc_indices)
code += self._array_pattern(indexing._index_dtype, index_array_symbol.name, acc_indices)
for d, arrname in enumerate(indexing._offset_array_names(f_dir, inv)):
offset_array_symbols = indexing._offset_array_symbols(f_dir, inv)
symbols_defined |= set(offset_array_symbols)
for d, arrsymb in enumerate(offset_array_symbols):
acc_offsets = ", ".join([str(o) for o in offsets[d]])
code += self._array_pattern(indexing._offsets_dtype, arrname, acc_offsets)
code += self._array_pattern(indexing._offsets_dtype, arrsymb.name, acc_offsets)
super(AdvancedStreamingBoundaryIndexing.TranslationArraysNode, self).__init__(
code, symbols_read=set(), symbols_defined=indexing._required_arrays)
code, symbols_read=set(), symbols_defined=symbols_defined)
def _array_pattern(self, dtype, name, content):
return f"const {str(dtype)} {name} [] = {{ {content} }}; \n"
def __str__(self):
return "Boundary Access Translation Arrays"
def __repr__(self):
return "Boundary Access Translation Arrays"
# end class AdvancedStreamingIndexing
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, AssignmentCollection
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]
@staticmethod
def stream_out_pdfs(stencil):
return sp.symbols('f_out_:{}'.format(len(stencil)))
@staticmethod
def stream_in_pdfs(stencil):
return sp.symbols('f_in_:{}'.format(len(stencil)))
@staticmethod
def stream_out_pdf_symbolic_dir_indexed(dir_index):
return sp.IndexedBase('f_out', shape=(1,))[dir_index]
@staticmethod
def stream_in_pdf_symbolic_inverse_dir_indexed(dir_index):
return sp.IndexedBase('f_in', shape=(1,))[dir_index]
@staticmethod
def inverse_dir_symbol(dir_symbol):
symbol_type = dir_symbol.dtype if isinstance(dir_symbol, TypedSymbol) else create_type(np.int64)
return TypedSymbol("{}_inv".format(dir_symbol.name), symbol_type)
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
self.pdf_field = pdf_field
self.stencil = stencil
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)
self.write_field_accesses = write_field_accesses
self.read_field_accesses = read_field_accesses
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 int64_t %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 int64_t %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)
def substitute_pdf_proxies(self, boundary_assignments, dir_symbol):
if not isinstance(boundary_assignments, AssignmentCollection):
boundary_assignments = AssignmentCollection(boundary_assignments)
out_pdfs = self.stream_out_pdfs(self.stencil)
in_pdfs = self.stream_in_pdfs(self.stencil)
dir_indexed_out_pdf = self.stream_out_pdf_symbolic_dir_indexed(dir_symbol)
inv_dir_symbol = self.inverse_dir_symbol(dir_symbol)
inv_dir_indexed_in_pdf = self.stream_in_pdf_symbolic_inverse_dir_indexed(inv_dir_symbol)
q = len(self.stencil)
dim = len(self.stencil[0])
pdf_subs = dict()
# Direct accesses
for i in range(q):
pdf_subs[out_pdfs[i]] = self.write_field_accesses[i]
pdf_subs[in_pdfs[i]] = self.read_field_accesses[i]
# Symbolically indexed accesses
outward_offset = self.outward_offset_from_dir(dir_symbol, dim)
outward_idx = self.outward_index_from_dir(dir_symbol)
inv_inward_offset = self.inverse_inward_offset_from_dir(dir_symbol, dim)
inv_inward_idx = self.inverse_inward_index_from_dir(dir_symbol)
pdf_out_access_by_dir = self.pdf_field[outward_offset](outward_idx)
pdf_subs[dir_indexed_out_pdf] = pdf_out_access_by_dir
pdf_inv_in_access_by_dir = self.pdf_field[inv_inward_offset](inv_inward_idx)
pdf_subs[inv_dir_indexed_in_pdf] = pdf_inv_in_access_by_dir
return boundary_assignments.new_with_substitutions(pdf_subs)
@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.flexible_boundaries_prototype import FlexibleNoSlip, create_flexible_lbm_boundary_kernel
from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor, \
StreamPushTwoFieldsAccessor, \
AAEvenTimeStepAccessor, \
AAOddTimeStepAccessor, \
EsoTwistEvenTimeStepAccessor, \
EsoTwistOddTimeStepAccessor
from lbmpy.advanced_streaming import FlexibleNoSlip, create_flexible_lbm_boundary_kernel
from lbmpy.creationfunctions import create_lb_method
from lbmpy.stencils import get_stencil
......@@ -19,7 +25,7 @@ from itertools import product
# INFRASTRUCTURE
# ===========================
stencils = [ 'D2Q9', 'D3Q19', 'D3Q27']
timesteps = ['odd_to_even', 'even_to_odd']
kernel_types = ['push', 'pull', 'aa', 'esotwist']
......@@ -50,21 +56,17 @@ class AccessPdfValues:
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):
def write_pdf(self, pdf_arr, pos, d, value):
offsets = self.write_accs[d].offsets
x += offsets[0]
y += offsets[1]
pos = tuple( p + o for p, o in zip(pos, offsets) )
i = self.write_accs[d].index[0]
pdf_arr[pos + (i,)] = value
pdf_arr[x,y,i] = value
def read_pdf(self, pdf_arr, x, y, d):
def read_pdf(self, pdf_arr, pos, d):
offsets = self.read_accs[d].offsets
x += offsets[0]
y += offsets[1]
pos = tuple( p + o for p, o in zip(pos, offsets) )
i = self.read_accs[d].index[0]
return pdf_arr[x,y,i]
return pdf_arr[pos + (i,)]
......@@ -74,17 +76,18 @@ class AccessPdfValues:
import pytest
arguments = product(timesteps, kernel_types)
arguments = product(stencils, timesteps, kernel_types)
@pytest.mark.parametrize("timestep_type, kernel_type", arguments)
def test_flexible_noslip(timestep_type, kernel_type):
@pytest.mark.parametrize("stencil, timestep_type, kernel_type", arguments)
def test_flexible_noslip_single_cell(stencil, timestep_type, kernel_type):
"""
Flexible NoSlip Test
"""
stencil = get_stencil('D2Q9')
pdf_field = ps.fields('pdfs(9): [2D]')
stencil = get_stencil(stencil)
q = len(stencil)
dim = len(stencil[0])
pdf_field = ps.fields(f'pdfs({q}): [{dim}D]')
even_accessor = even_accessors[kernel_type]
odd_accessor = odd_accessors[kernel_type]
......@@ -95,18 +98,19 @@ def test_flexible_noslip(timestep_type, kernel_type):
prev_pdf_access = AccessPdfValues(pdf_field, stencil, prev_accessor)
next_pdf_access = AccessPdfValues(pdf_field, stencil, next_accessor)
pdfs = np.zeros((3,3,9))
pdfs = np.zeros((3,) * dim + (q,))
pos = (1,) * dim
for d in range(q):
prev_pdf_access.write_pdf(pdfs, 1, 1, d, d)
prev_pdf_access.write_pdf(pdfs, pos, d, d)
lb_method = create_lb_method()
lb_method = create_lb_method(stencil=stencil, method='srt')
flex_noslip = FlexibleNoSlip()
index_struct_dtype = numpy_data_type_for_boundary_object(flex_noslip, lb_method.dim)
index_struct_dtype = numpy_data_type_for_boundary_object(flex_noslip, 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)
index_vector = np.array([ pos + (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,
......@@ -117,10 +121,6 @@ def test_flexible_noslip(timestep_type, kernel_type):
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
reflected_pdfs = [ next_pdf_access.read_pdf(pdfs, pos, d) for d in range(q)]
inverse_pdfs = [ inverse_dir_index(stencil, d) for d in range(q) ]
assert reflected_pdfs == inverse_pdfs
\ 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