Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Frederik Hennig
lbmpy
Commits
dd1cbd09
Commit
dd1cbd09
authored
Oct 24, 2020
by
Frederik Hennig
Browse files
replaced boundary handling by new implementation
parent
c91bf2ac
Changes
11
Expand all
Hide whitespace changes
Inline
Side-by-side
lbmpy/advanced_streaming/__init__.py
View file @
dd1cbd09
from
.indexing
import
BetweenTimestepsIndexing
from
.boundaryconditions
import
FlexibleBoundary
,
FlexibleNoSlip
from
.boundaryhandling
import
FlexibleLBMBoundaryHandling
,
\
create_advanced_streaming_boundary_kernel
from
.communication
import
get_communication_slices
,
PeriodicityHandling
from
.utility
import
Timestep
,
get_accessor
__all__
=
[
'BetweenTimestepsIndexing'
,
'FlexibleBoundary'
,
'FlexibleNoSlip'
,
'create_advanced_streaming_boundary_kernel'
,
'FlexibleLBMBoundaryHandling'
,
__all__
=
[
'BetweenTimestepsIndexing'
,
'get_communication_slices'
,
'PeriodicityHandling'
,
'Timestep'
,
'get_accessor'
]
lbmpy/advanced_streaming/boundaryconditions.py
deleted
100644 → 0
View file @
c91bf2ac
import
sympy
as
sp
from
pystencils
import
Assignment
,
Field
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
class
FlexibleBoundary
:
"""Base class for all boundary classes that use the AdvancedStreamingBoundaryIndexing"""
inner_or_boundary
=
True
single_link
=
False
def
__init__
(
self
,
name
=
None
):
self
.
_name
=
name
def
__call__
(
self
,
f_out
,
f_in
,
dir_symbol
,
inv_dir
,
lb_method
,
index_field
):
"""
This function defines the boundary behavior and must therefore be implemented by all boundaries.
The boundary is defined through a list of sympy equations from which a boundary kernel is generated.
:param 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
:param 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
:param 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.
:param 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.
:param lb_method: an instance of the LB method used. Use this to adapt the boundary to the method
(e.g. compressibility)
:param index_field: the boundary index field that can be used to retrieve and update boundary data
:return: list of pystencils assignments, or pystencils.AssignmentCollection
"""
raise
NotImplementedError
(
"Boundary class has to overwrite __call__"
)
@
property
def
additional_data
(
self
):
"""Return a list of (name, type) tuples for additional data items required in this boundary
These data items can either be initialized in separate kernel see additional_data_kernel_init or by
Python callbacks - see additional_data_callback """
return
[]
@
property
def
additional_data_init_callback
(
self
):
"""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"""
return
None
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."""
return
[]
@
property
def
name
(
self
):
if
self
.
_name
:
return
self
.
_name
else
:
return
type
(
self
).
__name__
@
name
.
setter
def
name
(
self
,
new_value
):
self
.
_name
=
new_value
# end class FlexibleBoundary
class
FlexibleNoSlip
(
FlexibleBoundary
):
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
,
f_out
,
f_in
,
dir_symbol
,
inv_dir
,
lb_method
,
index_field
):
return
Assignment
(
f_in
(
inv_dir
[
dir_symbol
]),
f_out
(
dir_symbol
))
def
__hash__
(
self
):
return
hash
(
self
.
name
)
def
__eq__
(
self
,
other
):
if
not
isinstance
(
other
,
FlexibleNoSlip
):
return
False
return
self
.
name
==
other
.
name
# end class FlexibleNoSlip
class
FlexibleUBB
(
FlexibleBoundary
):
"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""
def
__init__
(
self
,
velocity
,
adapt_velocity_to_force
=
False
,
dim
=
None
,
name
=
None
):
"""
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:
"""
super
(
FlexibleUBB
,
self
).
__init__
(
name
)
self
.
_velocity
=
velocity
self
.
_adaptVelocityToForce
=
adapt_velocity_to_force
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
@
property
def
additional_data
(
self
):
if
callable
(
self
.
_velocity
):
return
[(
'vel_%d'
%
(
i
,),
create_type
(
"double"
))
for
i
in
range
(
self
.
dim
)]
else
:
return
[]
@
property
def
additional_data_init_callback
(
self
):
if
callable
(
self
.
_velocity
):
return
self
.
_velocity
@
property
def
get_additional_code_nodes
(
self
,
lb_method
):
return
[
LbmWeightInfo
(
lb_method
),
NeighbourOffsetArraysForStencil
(
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
}
)"
neighbor_offset
=
NeighbourOffsetArraysForStencil
.
symbolic_neighbour_offset_from_dir
(
direction
,
lb_method
.
dim
)
velocity
=
tuple
(
v_i
.
get_shifted
(
*
neighbor_offset
)
if
isinstance
(
v_i
,
Field
.
Access
)
and
not
vel_from_idx_field
else
v_i
for
v_i
in
vel
)
if
self
.
_adaptVelocityToForce
:
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
]
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
)
# Better alternative: in conserved value computation
# rename what is currently called density to "virtual_density"
# provide a new quantity density, which is constant in case of incompressible models
if
not
lb_method
.
conserved_quantity_computation
.
zero_centered_pdfs
:
cqc
=
lb_method
.
conserved_quantity_computation
density_symbol
=
sp
.
Symbol
(
"rho"
)
pdf_field_accesses
=
[
f_out
(
i
)
for
i
in
range
(
len
(
lb_method
.
stencil
))]
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
result
+=
[
Assignment
(
f_in
(
inv_dir
[
direction
]),
f_out
(
direction
)
-
vel_term
*
density_symbol
)]
return
result
else
:
return
[
Assignment
(
f_in
(
inv_dir
[
direction
]),
f_out
(
direction
)
-
vel_term
)]
lbmpy/advanced_streaming/boundaryhandling.py
deleted
100644 → 0
View file @
c91bf2ac
import
numpy
as
np
from
lbmpy.advanced_streaming.indexing
import
BetweenTimestepsIndexing
from
lbmpy.advanced_streaming.utility
import
is_inplace
,
Timestep
,
AccessPdfValues
from
pystencils
import
Field
,
Assignment
,
create_indexed_kernel
from
pystencils.stencil
import
inverse_direction
from
pystencils.boundaries
import
BoundaryHandling
from
pystencils.boundaries.createindexlist
import
numpy_data_type_for_boundary_object
class
FlexibleLBMBoundaryHandling
(
BoundaryHandling
):
"""
Enables boundary handling for LBM simulations with advanced streaming patterns.
For the in-place patterns AA and EsoTwist, two kernels are generated for a boundary
object and the right one selected depending on the time step.
"""
def
__init__
(
self
,
lb_method
,
data_handling
,
pdf_field_name
,
streaming_pattern
=
'pull'
,
name
=
"boundary_handling"
,
flag_interface
=
None
,
target
=
'cpu'
,
openmp
=
True
):
self
.
_lb_method
=
lb_method
self
.
_streaming_pattern
=
streaming_pattern
self
.
_inplace
=
is_inplace
(
streaming_pattern
)
self
.
_prev_timestep
=
None
super
(
FlexibleLBMBoundaryHandling
,
self
).
__init__
(
data_handling
,
pdf_field_name
,
lb_method
.
stencil
,
name
,
flag_interface
,
target
,
openmp
)
# ------------------------- Overridden methods of pystencils.BoundaryHandling -------------------------
@
property
def
prev_timestep
(
self
):
return
self
.
_prev_timestep
def
__call__
(
self
,
prev_timestep
=
Timestep
.
BOTH
,
**
kwargs
):
self
.
_prev_timestep
=
prev_timestep
super
(
FlexibleLBMBoundaryHandling
,
self
).
__call__
(
**
kwargs
)
self
.
_prev_timestep
=
None
def
add_fixed_steps
(
self
,
fixed_loop
,
**
kwargs
):
raise
NotImplementedError
(
"Adding to fixed loop is not supported by FlexibleLBMBoundaryHandling"
)
def
_add_boundary
(
self
,
boundary_obj
,
flag
=
None
):
if
self
.
_inplace
:
return
self
.
_add_flexible_boundary
(
boundary_obj
,
flag
)
else
:
return
super
(
FlexibleLBMBoundaryHandling
,
self
).
_add_boundary
(
boundary_obj
,
flag
)
def
_add_flexible_boundary
(
self
,
boundary_obj
,
flag
=
None
):
if
boundary_obj
not
in
self
.
_boundary_object_to_boundary_info
:
sym_index_field
=
Field
.
create_generic
(
'indexField'
,
spatial_dimensions
=
1
,
dtype
=
numpy_data_type_for_boundary_object
(
boundary_obj
,
self
.
dim
))
kernels
=
[
self
.
_create_boundary_kernel
(
self
.
_data_handling
.
fields
[
self
.
_field_name
],
sym_index_field
,
boundary_obj
,
Timestep
.
EVEN
).
compile
(),
self
.
_create_boundary_kernel
(
self
.
_data_handling
.
fields
[
self
.
_field_name
],
sym_index_field
,
boundary_obj
,
Timestep
.
ODD
).
compile
()]
if
flag
is
None
:
flag
=
self
.
flag_interface
.
reserve_next_flag
()
boundary_info
=
self
.
InplaceStreamingBoundaryInfo
(
self
,
boundary_obj
,
flag
,
kernels
)
self
.
_boundary_object_to_boundary_info
[
boundary_obj
]
=
boundary_info
return
self
.
_boundary_object_to_boundary_info
[
boundary_obj
].
flag
def
_create_boundary_kernel
(
self
,
symbolic_field
,
symbolic_index_field
,
boundary_obj
,
prev_timestep
=
Timestep
.
BOTH
):
return
create_advanced_streaming_boundary_kernel
(
symbolic_field
,
symbolic_index_field
,
self
.
_lb_method
,
boundary_obj
,
prev_timestep
=
prev_timestep
,
streaming_pattern
=
self
.
_streaming_pattern
,
target
=
self
.
_target
,
openmp
=
self
.
_openmp
)
class
InplaceStreamingBoundaryInfo
(
object
):
@
property
def
kernel
(
self
):
prev_timestep
=
self
.
_boundary_handling
.
prev_timestep
if
prev_timestep
is
None
:
raise
Exception
(
"The boundary kernel property was accessed while "
+
"there was no boundary handling in progress."
)
return
self
.
_kernels
[
prev_timestep
]
def
__init__
(
self
,
boundary_handling
,
boundary_obj
,
flag
,
kernels
):
self
.
_boundary_handling
=
boundary_handling
self
.
boundary_object
=
boundary_obj
self
.
flag
=
flag
self
.
_kernels
=
kernels
# end class InplaceStreamingBoundaryInfo
# ------------------------------ Force On Boundary ------------------------------------------------------------
def
force_on_boundary
(
self
,
boundary_obj
,
prev_timestep
=
Timestep
.
BOTH
):
from
lbmpy.advanced_streaming
import
FlexibleNoSlip
if
isinstance
(
boundary_obj
,
FlexibleNoSlip
):
return
self
.
_force_on_no_slip
(
boundary_obj
,
prev_timestep
)
else
:
self
.
__call__
()
return
self
.
_force_on_boundary
(
boundary_obj
,
prev_timestep
)
def
_force_on_no_slip
(
self
,
boundary_obj
,
prev_timestep
):
dh
=
self
.
_data_handling
ff_ghost_layers
=
dh
.
ghost_layers_of_field
(
self
.
flag_interface
.
flag_field_name
)
method
=
self
.
_lb_method
stencil
=
np
.
array
(
method
.
stencil
)
result
=
np
.
zeros
(
self
.
dim
)
for
b
in
dh
.
iterate
(
ghost_layers
=
ff_ghost_layers
):
obj_to_ind_list
=
b
[
self
.
_index_array_name
].
boundary_object_to_index_list
pdf_array
=
b
[
self
.
_field_name
]
if
boundary_obj
in
obj_to_ind_list
:
ind_arr
=
obj_to_ind_list
[
boundary_obj
]
acc
=
AccessPdfValues
(
dh
.
fields
[
self
.
_field_name
],
self
.
_lb_method
.
stencil
,
streaming_pattern
=
self
.
_streaming_pattern
,
timestep
=
prev_timestep
,
streaming_dir
=
'out'
)
values
=
2
*
acc
.
collect_from_index_list
(
pdf_array
,
ind_arr
)
forces
=
stencil
[
ind_arr
[
'dir'
]]
*
values
[:,
np
.
newaxis
]
result
+=
forces
.
sum
(
axis
=
0
)
return
dh
.
reduce_float_sequence
(
list
(
result
),
'sum'
)
def
_force_on_boundary
(
self
,
boundary_obj
,
prev_timestep
):
dh
=
self
.
_data_handling
ff_ghost_layers
=
dh
.
ghost_layers_of_field
(
self
.
flag_interface
.
flag_field_name
)
method
=
self
.
_lb_method
stencil
=
np
.
array
(
method
.
stencil
)
result
=
np
.
zeros
(
self
.
dim
)
for
b
in
dh
.
iterate
(
ghost_layers
=
ff_ghost_layers
):
obj_to_ind_list
=
b
[
self
.
_index_array_name
].
boundary_object_to_index_list
pdf_array
=
b
[
self
.
_field_name
]
if
boundary_obj
in
obj_to_ind_list
:
ind_arr
=
obj_to_ind_list
[
boundary_obj
]
acc_out
=
AccessPdfValues
(
dh
.
fields
[
self
.
_field_name
],
self
.
_lb_method
.
stencil
,
streaming_pattern
=
self
.
_streaming_pattern
,
timestep
=
prev_timestep
,
streaming_dir
=
'out'
)
acc_in
=
AccessPdfValues
(
dh
.
fields
[
self
.
_field_name
],
self
.
_lb_method
.
stencil
,
streaming_pattern
=
self
.
_streaming_pattern
,
timestep
=
prev_timestep
.
next
(),
streaming_dir
=
'in'
)
acc_fluid
=
acc_out
if
boundary_obj
.
inner_or_boundary
else
acc_in
acc_boundary
=
acc_in
if
boundary_obj
.
inner_or_boundary
else
acc_out
fluid_values
=
acc_fluid
.
collect_from_index_list
(
pdf_array
,
ind_arr
)
boundary_values
=
acc_boundary
.
collect_from_index_list
(
pdf_array
,
ind_arr
)
values
=
fluid_values
+
boundary_values
forces
=
stencil
[
ind_arr
[
'dir'
]]
*
values
[:,
np
.
newaxis
]
result
+=
forces
.
sum
(
axis
=
0
)
return
dh
.
reduce_float_sequence
(
list
(
result
),
'sum'
)
# end class FlexibleLBMBoundaryHandling
def
create_advanced_streaming_boundary_kernel
(
pdf_field
,
index_field
,
lb_method
,
boundary_functor
,
prev_timestep
=
Timestep
.
BOTH
,
streaming_pattern
=
'pull'
,
target
=
'cpu'
,
openmp
=
True
):
index_dtype
=
index_field
.
dtype
.
numpy_dtype
.
fields
[
'dir'
][
0
]
offsets_dtype
=
index_field
.
dtype
.
numpy_dtype
.
fields
[
'x'
][
0
]
indexing
=
BetweenTimestepsIndexing
(
pdf_field
,
lb_method
.
stencil
,
prev_timestep
,
streaming_pattern
,
index_dtype
,
offsets_dtype
)
f_out
,
f_in
=
indexing
.
proxy_fields
dir_symbol
=
indexing
.
dir_symbol
inv_dir
=
indexing
.
inverse_dir_symbol
boundary_assignments
=
boundary_functor
(
f_out
,
f_in
,
dir_symbol
,
inv_dir
,
lb_method
,
index_field
)
boundary_assignments
=
indexing
.
substitute_proxies
(
boundary_assignments
)
# Code Elements inside the loop
elements
=
[
Assignment
(
dir_symbol
,
index_field
[
0
](
'dir'
))]
elements
+=
boundary_assignments
.
all_assignments
kernel
=
create_indexed_kernel
(
elements
,
[
index_field
],
target
=
target
,
cpu_openmp
=
openmp
)
# Code Elements ahead of the loop
index_arrs_node
=
indexing
.
create_code_node
()
for
node
in
boundary_functor
.
get_additional_code_nodes
(
lb_method
)[::
-
1
]:
kernel
.
body
.
insert_front
(
node
)
kernel
.
body
.
insert_front
(
index_arrs_node
)
return
kernel
lbmpy/boundaries/boundaryconditions.py
View file @
dd1cbd09
import
sympy
as
sp
from
lbmpy.boundaries.boundaryhandling
import
BoundaryOffsetInfo
,
LbmWeightInfo
from
lbmpy.simplificationfactory
import
create_simplification_strategy
from
pystencils
import
Assignment
,
Field
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
class
Boundary
:
"""Base class all boundar
ies should derive from
"""
"""Base class
for
all boundar
y classes that use the AdvancedStreamingBoundaryIndexing
"""
inner_or_boundary
=
True
single_link
=
False
...
...
@@ -16,23 +16,25 @@ class Boundary:
def
__init__
(
self
,
name
=
None
):
self
.
_name
=
name
def
__call__
(
self
,
pdf_field
,
direction_symbol
,
lb_method
,
index_field
):
def
__call__
(
self
,
f_out
,
f_in
,
dir_symbol
,
inv_dir
,
lb_method
,
index_field
):
"""
This function defines the boundary behavior and must therefore be implemented by all boundaries.
Here the boundary is defined as a list of sympy equations, from which a boundary kernel is generated.
Args:
pdf_field: pystencils field describing the pdf. The current cell is cell next to the boundary,
which is influenced by the boundary cell i.e. has a link from the boundary cell to
itself.
direction_symbol: a sympy symbol that can be used as index to the pdf_field. It describes
the direction pointing from the fluid to the boundary cell
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
Returns:
:return: list of sympy equations
The boundary is defined through a list of sympy equations from which a boundary kernel is generated.
:param 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
:param 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
:param 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.
:param 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.
:param lb_method: an instance of the LB method used. Use this to adapt the boundary to the method
(e.g. compressibility)
:param index_field: the boundary index field that can be used to retrieve and update boundary data
:return: list of pystencils assignments, or pystencils.AssignmentCollection
"""
raise
NotImplementedError
(
"Boundary class has to overwrite __call__"
)
...
...
@@ -49,6 +51,10 @@ class Boundary:
data-name to data for each element that should be initialized"""
return
None
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."""
return
[]
@
property
def
name
(
self
):
if
self
.
_name
:
...
...
@@ -60,6 +66,8 @@ class Boundary:
def
name
(
self
,
new_value
):
self
.
_name
=
new_value
# end class Boundary
class
NoSlip
(
Boundary
):
...
...
@@ -67,14 +75,15 @@ class NoSlip(Boundary):
"""Set an optional name here, to mark boundaries, for example for force evaluations"""
super
(
NoSlip
,
self
).
__init__
(
name
)
"""No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle"""
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
))]
"""
No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle.
Extended for use with any streaming pattern.
"""
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
))
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
):
...
...
@@ -82,6 +91,8 @@ class NoSlip(Boundary):
return
False
return
self
.
name
==
other
.
name
# end class NoSlip
class
UBB
(
Boundary
):
"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""
...
...
@@ -115,18 +126,23 @@ class UBB(Boundary):
if
callable
(
self
.
_velocity
):
return
self
.
_velocity
def
__call__
(
self
,
pdf_field
,
direction_symbol
,
lb_method
,
index_field
,
**
kwargs
):
def
get_additional_code_nodes
(
self
,
lb_method
):
return
[
LbmWeightInfo
(
lb_method
),
NeighbourOffsetArraysForStencil
(
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
(
'vel_
%d'
%
(
i
,)
)
for
i
in
range
(
self
.
dim
)]
if
vel_from_idx_field
else
self
.
_velocity
direction
=
dir
ection
_symbol
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
,
"Dimension of UBB (%d) does not match dimension of method (%d)"
\
%
(
self
.
dim
,
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
=
BoundaryOffsetInfo
.
offset_from_dir
(
direction
,
lb_method
.
dim
)
inverse_dir
=
BoundaryOffsetInfo
.
inv_dir
(
direction
)
neighbor
_offset
=
NeighbourOffsetArraysForStencil
.
symbolic_neighbour_offset_from_dir
(
direction
,
lb_method
.
dim
)
velocity
=
tuple
(
v_i
.
get_shifted
(
*
neighbor
)
if
isinstance
(
v_i
,
Field
.
Access
)
and
not
vel_from_idx_field
else
v_i
velocity
=
tuple
(
v_i
.
get_shifted
(
*
neighbor_offset
)
if
isinstance
(
v_i
,
Field
.
Access
)
and
not
vel_from_idx_field
else
v_i
for
v_i
in
vel
)
if
self
.
_adaptVelocityToForce
:
...
...
@@ -136,8 +152,9 @@ 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
,
velocity
)])
*
weight_of_direction
(
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
)
# Better alternative: in conserved value computation
# rename what is currently called density to "virtual_density"
...
...
@@ -145,21 +162,23 @@ class UBB(Boundary):
if
not
lb_method
.
conserved_quantity_computation
.
zero_centered_pdfs
:
cqc
=
lb_method
.
conserved_quantity_computation
density_symbol
=
sp
.
Symbol
(
"rho"
)
pdf_field_accesses
=
[
pdf_field
(
i
)
for
i
in
range
(
len
(
lb_method
.
stencil
))]
pdf_field_accesses
=
[
f_out
(
i
)
for
i
in
range
(
len
(
lb_method
.
stencil
))]
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
result
+=
[
Assignment
(
pdf_field
[
neighbor
](
inverse_dir
),
pdf_field
(
direction
)
-
vel_term
*
density_symbol
)]
result
+=
[
Assignment
(
f_in
(
inv_dir
[
direction
]
),
f_out
(
direction
)
-
vel_term
*
density_symbol
)]
return
result
else
:
return
[
Assignment
(
pdf_field
[
neighbor
](
inverse_dir
),
pdf_field
(
direction
)
-
vel_term
)]
return
[
Assignment
(
f_in
(
inv_dir
[
direction
]
),
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
)
...
...
@@ -208,6 +227,7 @@ class FixedDensity(Boundary):
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
)),
...
...
@@ -223,6 +243,7 @@ class NeumannByCopy(Boundary):
class
StreamInConstant
(
Boundary
):
def
__init__
(
self
,
constant
,
name
=
None
):
raise
NotImplementedError
()
# not yet operable
super
(
StreamInConstant
,
self
).
__init__
(
name
)
self
.
_constant
=
constant
...
...
lbmpy/boundaries/boundaryhandling.py
View file @
dd1cbd09
import
numpy
as
np
import
sympy
as
sp
from
pystencils
import
Assignment
,
TypedSymbol
,
create_indexed_kernel
from
pystencils.backends.cbackend
import
CustomCodeNode
from
pystencils.boundaries
import
BoundaryHandling
from
pystencils.boundaries.boundaryhandling
import
BoundaryOffsetInfo
from
lbmpy.advanced_streaming.indexing
import
BetweenTimestepsIndexing
from
lbmpy.advanced_streaming.utility
import
is_inplace
,
Timestep
,
AccessPdfValues
from
pystencils
import
Field
,
Assignment
,
TypedSymbol
,
create_indexed_kernel
from
pystencils.stencil
import
inverse_direction
from
pystencils.boundaries
import
BoundaryHandling
from
pystencils.boundaries.createindexlist
import
numpy_data_type_for_boundary_object
from
pystencils.backends.cbackend
import
CustomCodeNode
class
LatticeBoltzmannBoundaryHandling
(
BoundaryHandling
):
def
__init__
(
self
,
lb_method
,
data_handling
,
pdf_field_name
,
name
=
"boundary_handling"
,
flag_interface
=
None
,
target
=
'cpu'
,
openmp
=
True
):
self
.
lb_method
=
lb_method
"""
Enables boundary handling for LBM simulations with advanced streaming patterns.
For the in-place patterns AA and EsoTwist, two kernels are generated for a boundary
object and the right one selected depending on the time step.
"""