Skip to content
Snippets Groups Projects
Commit 4dd3683d authored by Markus Holzer's avatar Markus Holzer
Browse files

Merge branch 'fix-freeslip' into 'master'

Fix FreeSlip boundary condition

See merge request !115
parents 41fbad6f dfe9776b
Branches
No related tags found
1 merge request!115Fix FreeSlip boundary condition
Pipeline #38110 failed
...@@ -120,9 +120,10 @@ class FreeSlip(LbBoundary): ...@@ -120,9 +120,10 @@ class FreeSlip(LbBoundary):
Args: Args:
stencil: LBM stencil which is used for the simulation stencil: LBM stencil which is used for the simulation
normal_direction: optional normal direction. If the Free slip boundary is applied to a certain side in the normal_direction: optional normal direction pointing from wall to fluid.
domain it is not necessary to calculate the normal direction since it can be stated for all If the Free slip boundary is applied to a certain side in the domain it is not necessary
boundary cells. This reduces the memory space for the index array significantly. to calculate the normal direction since it can be stated for all boundary cells.
This reduces the memory space for the index array significantly.
name: optional name of the boundary. name: optional name of the boundary.
""" """
...@@ -182,7 +183,12 @@ class FreeSlip(LbBoundary): ...@@ -182,7 +183,12 @@ class FreeSlip(LbBoundary):
normal_direction[i] = direction[i] normal_direction[i] = direction[i]
ref_direction = MirroredStencilDirections.mirror_stencil(ref_direction, i) ref_direction = MirroredStencilDirections.mirror_stencil(ref_direction, i)
ref_direction = inverse_direction(ref_direction) # convex corner special case:
if all(n == 0 for n in normal_direction):
normal_direction = direction
else:
ref_direction = inverse_direction(ref_direction)
for i, cell_name in zip(range(dim), self.additional_data): for i, cell_name in zip(range(dim), self.additional_data):
cell[cell_name[0]] = -normal_direction[i] cell[cell_name[0]] = -normal_direction[i]
cell['ref_dir'] = self.stencil.index(ref_direction) cell['ref_dir'] = self.stencil.index(ref_direction)
...@@ -208,13 +214,14 @@ class FreeSlip(LbBoundary): ...@@ -208,13 +214,14 @@ class FreeSlip(LbBoundary):
def get_additional_code_nodes(self, lb_method): def get_additional_code_nodes(self, lb_method):
if self.normal_direction: if self.normal_direction:
return [MirroredStencilDirections(self.stencil, self.mirror_axis)] return [MirroredStencilDirections(self.stencil, self.mirror_axis), NeighbourOffsetArrays(lb_method.stencil)]
else: else:
return [] return [NeighbourOffsetArrays(lb_method.stencil)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
if self.normal_direction: if self.normal_direction:
normal_direction = self.normal_direction tangential_offset = tuple(offset + normal for offset, normal in zip(neighbor_offset, self.normal_direction))
mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(self.mirror_axis) mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(self.mirror_axis)
mirrored_direction = inv_dir[sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]] mirrored_direction = inv_dir[sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]]
else: else:
...@@ -222,10 +229,11 @@ class FreeSlip(LbBoundary): ...@@ -222,10 +229,11 @@ class FreeSlip(LbBoundary):
for i, cell_name in zip(range(self.dim), self.additional_data): for i, cell_name in zip(range(self.dim), self.additional_data):
normal_direction.append(index_field[0](cell_name[0])) normal_direction.append(index_field[0](cell_name[0]))
normal_direction = tuple(normal_direction) normal_direction = tuple(normal_direction)
tangential_offset = tuple(offset + normal for offset, normal in zip(neighbor_offset, normal_direction))
mirrored_direction = index_field[0]('ref_dir') mirrored_direction = index_field[0]('ref_dir')
return Assignment(f_in(inv_dir[dir_symbol]), f_in[normal_direction](mirrored_direction)) return Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](mirrored_direction))
# end class FreeSlip # end class FreeSlip
......
...@@ -112,6 +112,122 @@ def test_simple(target): ...@@ -112,6 +112,122 @@ def test_simple(target):
assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 8] == 5)) assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 8] == 5))
@pytest.mark.parametrize("given_normal", [True, False])
def test_free_slip(given_normal):
# check if Free slip BC is applied correctly
stencil = LBStencil(Stencil.D2Q9)
dh = create_data_handling(domain_size=(4, 4),)
src1 = dh.add_array('src1', values_per_cell=stencil.Q)
dh.fill('src1', 0.0, ghost_layers=True)
shape = dh.gather_array('src1', ghost_layers=True).shape
num = 0
for x in range(shape[0]):
for y in range(shape[1]):
for direction in range(shape[2]):
dh.cpu_arrays[src1.name][x, y, direction] = num
num += 1
method = create_lb_method(lbm_config=LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.8))
bh = LatticeBoltzmannBoundaryHandling(method, dh, 'src1', name="bh1")
if given_normal:
free_slipN = FreeSlip(stencil=stencil, normal_direction=(0, -1))
free_slipS = FreeSlip(stencil=stencil, normal_direction=(0, 1))
free_slipE = FreeSlip(stencil=stencil, normal_direction=(-1, 0))
free_slipW = FreeSlip(stencil=stencil, normal_direction=(1, 0))
bh.set_boundary(free_slipN, slice_from_direction('N', dh.dim))
bh.set_boundary(free_slipS, slice_from_direction('S', dh.dim))
bh.set_boundary(free_slipE, slice_from_direction('E', dh.dim))
bh.set_boundary(free_slipW, slice_from_direction('W', dh.dim))
else:
free_slip = FreeSlip(stencil=stencil)
bh.set_boundary(free_slip, slice_from_direction('N', dh.dim))
bh.set_boundary(free_slip, slice_from_direction('S', dh.dim))
bh.set_boundary(free_slip, slice_from_direction('E', dh.dim))
bh.set_boundary(free_slip, slice_from_direction('W', dh.dim))
bh()
mirrored_dirN = {6: 8, 1: 2, 5: 7}
mirrored_dirS = {7: 5, 2: 1, 8: 6}
mirrored_dirE = {6: 5, 4: 3, 8: 7}
mirrored_dirW = {5: 6, 3: 4, 7: 8}
# check North
assert dh.cpu_arrays[src1.name][1, -1, mirrored_dirN[6]] == dh.cpu_arrays[src1.name][1, -2, 6]
assert dh.cpu_arrays[src1.name][1, -1, mirrored_dirN[1]] == dh.cpu_arrays[src1.name][1, -2, 1]
for i in range(2, 4):
assert dh.cpu_arrays[src1.name][i, -1, mirrored_dirN[6]] == dh.cpu_arrays[src1.name][i, -2, 6]
assert dh.cpu_arrays[src1.name][i, -1, mirrored_dirN[1]] == dh.cpu_arrays[src1.name][i, -2, 1]
assert dh.cpu_arrays[src1.name][i, -1, mirrored_dirN[5]] == dh.cpu_arrays[src1.name][i, -2, 5]
assert dh.cpu_arrays[src1.name][4, -1, mirrored_dirN[1]] == dh.cpu_arrays[src1.name][4, -2, 1]
assert dh.cpu_arrays[src1.name][4, -1, mirrored_dirN[5]] == dh.cpu_arrays[src1.name][4, -2, 5]
# check East
assert dh.cpu_arrays[src1.name][-1, 1, mirrored_dirE[6]] == dh.cpu_arrays[src1.name][-2, 1, 6]
assert dh.cpu_arrays[src1.name][-1, 1, mirrored_dirE[4]] == dh.cpu_arrays[src1.name][-2, 1, 4]
for i in range(2, 4):
assert dh.cpu_arrays[src1.name][-1, i, mirrored_dirE[6]] == dh.cpu_arrays[src1.name][-2, i, 6]
assert dh.cpu_arrays[src1.name][-1, i, mirrored_dirE[4]] == dh.cpu_arrays[src1.name][-2, i, 4]
assert dh.cpu_arrays[src1.name][-1, i, mirrored_dirE[8]] == dh.cpu_arrays[src1.name][-2, i, 8]
assert dh.cpu_arrays[src1.name][-1, 4, mirrored_dirE[4]] == dh.cpu_arrays[src1.name][-2, 4, 4]
assert dh.cpu_arrays[src1.name][-1, 4, mirrored_dirE[8]] == dh.cpu_arrays[src1.name][-2, 4, 8]
# check South
assert dh.cpu_arrays[src1.name][1, 0, mirrored_dirS[8]] == dh.cpu_arrays[src1.name][1, 1, 8]
assert dh.cpu_arrays[src1.name][1, 0, mirrored_dirS[2]] == dh.cpu_arrays[src1.name][1, 1, 2]
for i in range(2, 4):
assert dh.cpu_arrays[src1.name][i, 0, mirrored_dirS[7]] == dh.cpu_arrays[src1.name][i, 1, 7]
assert dh.cpu_arrays[src1.name][i, 0, mirrored_dirS[2]] == dh.cpu_arrays[src1.name][i, 1, 2]
assert dh.cpu_arrays[src1.name][i, 0, mirrored_dirS[8]] == dh.cpu_arrays[src1.name][i, 1, 8]
assert dh.cpu_arrays[src1.name][4, 0, mirrored_dirS[2]] == dh.cpu_arrays[src1.name][4, 1, 2]
assert dh.cpu_arrays[src1.name][4, 0, mirrored_dirS[7]] == dh.cpu_arrays[src1.name][4, 1, 7]
# check West
assert dh.cpu_arrays[src1.name][0, 1, mirrored_dirW[5]] == dh.cpu_arrays[src1.name][1, 1, 5]
assert dh.cpu_arrays[src1.name][0, 1, mirrored_dirW[3]] == dh.cpu_arrays[src1.name][1, 1, 3]
for i in range(2, 4):
assert dh.cpu_arrays[src1.name][0, i, mirrored_dirW[5]] == dh.cpu_arrays[src1.name][1, i, 5]
assert dh.cpu_arrays[src1.name][0, i, mirrored_dirW[3]] == dh.cpu_arrays[src1.name][1, i, 3]
assert dh.cpu_arrays[src1.name][0, i, mirrored_dirW[7]] == dh.cpu_arrays[src1.name][1, i, 7]
assert dh.cpu_arrays[src1.name][0, 4, mirrored_dirW[3]] == dh.cpu_arrays[src1.name][1, 4, 3]
assert dh.cpu_arrays[src1.name][0, 4, mirrored_dirW[7]] == dh.cpu_arrays[src1.name][1, 4, 7]
if given_normal:
# check corners --> determined by the last boundary applied there.
# SouthWest --> West
assert dh.cpu_arrays[src1.name][0, 0, mirrored_dirW[5]] == dh.cpu_arrays[src1.name][1, 0, 5]
# NorthWest --> West
assert dh.cpu_arrays[src1.name][0, -1, mirrored_dirW[7]] == dh.cpu_arrays[src1.name][1, -1, 7]
# NorthEast --> East
assert dh.cpu_arrays[src1.name][-1, -1, mirrored_dirE[8]] == dh.cpu_arrays[src1.name][-2, -1, 8]
# SouthEast --> East
assert dh.cpu_arrays[src1.name][-1, 0, mirrored_dirE[6]] == dh.cpu_arrays[src1.name][-2, 0, 6]
else:
# check corners --> this time the normals are calculated correctly in the corners
# SouthWest --> Normal = (1, 1); dir 7 --> 6
assert dh.cpu_arrays[src1.name][0, 0, 6] == dh.cpu_arrays[src1.name][1, 1, 7]
# NorthWest --> Normal = (1, -1); dir 8 --> 5
assert dh.cpu_arrays[src1.name][0, -1, 8] == dh.cpu_arrays[src1.name][1, -2, 5]
# NorthEast --> Normal = (-1, -1); dir 7 --> 6
assert dh.cpu_arrays[src1.name][-1, -1, 7] == dh.cpu_arrays[src1.name][-2, -2, 6]
# SouthEast --> Normal = (-1, 1); dir 5 --> 8
assert dh.cpu_arrays[src1.name][-1, 0, 5] == dh.cpu_arrays[src1.name][-2, 1, 8]
def test_free_slip_index_list(): def test_free_slip_index_list():
stencil = LBStencil(Stencil.D2Q9) stencil = LBStencil(Stencil.D2Q9)
dh = create_data_handling(domain_size=(4, 4), periodicity=(False, False)) dh = create_data_handling(domain_size=(4, 4), periodicity=(False, False))
...@@ -181,6 +297,50 @@ def test_free_slip_index_list(): ...@@ -181,6 +297,50 @@ def test_free_slip_index_list():
assert normal == normal_north_east assert normal == normal_north_east
def test_free_slip_index_list_convex_corner():
stencil = LBStencil(Stencil.D2Q9)
dh = create_data_handling(domain_size=(4, 4))
src = dh.add_array('src', values_per_cell=len(stencil))
dh.fill('src', 0.0, ghost_layers=True)
lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.8)
method = create_lb_method(lbm_config=lbm_config)
def bh_callback(x, y):
radius = 2
x_mid = 2
y_mid = 2
return (x - x_mid) ** 2 + (y - y_mid) ** 2 > radius ** 2
bh = LatticeBoltzmannBoundaryHandling(method, dh, 'src', name="bh")
free_slip = FreeSlip(stencil=stencil)
bh.set_boundary(free_slip, mask_callback=bh_callback)
bh.prepare()
for b in dh.iterate():
for b_obj, idx_arr in b[bh._index_array_name].boundary_object_to_index_list.items():
index_array = idx_arr
# correct index array for this case with convex corners
test = [(2, 1, 2, 0, 1, 2), (2, 1, 3, 1, 0, 3), (2, 1, 7, 1, 1, 7),
(2, 1, 8, 0, 1, 7), (3, 1, 2, 0, 1, 2), (3, 1, 4, -1, 0, 4),
(3, 1, 7, 0, 1, 8), (3, 1, 8, -1, 1, 8), (1, 2, 2, 0, 1, 2),
(1, 2, 3, 1, 0, 3), (1, 2, 5, 1, 0, 7), (1, 2, 7, 1, 1, 7),
(2, 2, 7, 1, 1, 7), (3, 2, 8, -1, 1, 8), (4, 2, 2, 0, 1, 2),
(4, 2, 4, -1, 0, 4), (4, 2, 6, -1, 0, 8), (4, 2, 8, -1, 1, 8),
(1, 3, 1, 0, -1, 1), (1, 3, 3, 1, 0, 3), (1, 3, 5, 1, -1, 5),
(1, 3, 7, 1, 0, 5), (2, 3, 5, 1, -1, 5), (3, 3, 6, -1, -1, 6),
(4, 3, 1, 0, -1, 1), (4, 3, 4, -1, 0, 4), (4, 3, 6, -1, -1, 6),
(4, 3, 8, -1, 0, 6), (2, 4, 1, 0, -1, 1), (2, 4, 3, 1, 0, 3),
(2, 4, 5, 1, -1, 5), (2, 4, 6, 0, -1, 5), (3, 4, 1, 0, -1, 1),
(3, 4, 4, -1, 0, 4), (3, 4, 5, 0, -1, 6), (3, 4, 6, -1, -1, 6)]
for i, cell in enumerate(index_array):
for j in range(len(cell)):
assert cell[j] == test[i][j]
def test_free_slip_equivalence(): def test_free_slip_equivalence():
# check if Free slip BC does the same if the normal direction is specified or not # check if Free slip BC does the same if the normal direction is specified or not
...@@ -214,7 +374,7 @@ def test_free_slip_equivalence(): ...@@ -214,7 +374,7 @@ def test_free_slip_equivalence():
bh1() bh1()
bh2() bh2()
assert np.array_equal(dh.cpu_arrays['src1'], dh.cpu_arrays['src2']) assert np.array_equal(dh.gather_array('src1'), dh.gather_array('src2'))
def test_exotic_boundaries(): def test_exotic_boundaries():
......
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment