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
Markus Holzer
lbmpy
Commits
482c8293
Commit
482c8293
authored
Apr 03, 2018
by
Martin Bauer
Browse files
PEP8 name refactoring
- test run again - notebooks not yet
parent
48999882
Changes
49
Expand all
Hide whitespace changes
Inline
Side-by-side
boundaries/boundaryconditions.py
View file @
482c8293
...
...
@@ -3,41 +3,41 @@ from pystencils import Field, Assignment
from
pystencils.astnodes
import
SympyAssignment
from
pystencils.sympyextensions
import
get_symmetric_part
from
pystencils.data_types
import
create_type
from
lbmpy.simplificationfactory
import
create
S
implification
S
trategy
from
lbmpy.simplificationfactory
import
create
_s
implification
_s
trategy
from
lbmpy.boundaries.boundaryhandling
import
BoundaryOffsetInfo
,
LbmWeightInfo
class
Boundary
(
object
)
:
class
Boundary
:
"""Base class all boundaries should derive from"""
def
__init__
(
self
,
name
=
None
):
self
.
_name
=
name
def
__call__
(
self
,
pdf
F
ield
,
direction
S
ymbol
,
lb
M
ethod
,
index
F
ield
):
def
__call__
(
self
,
pdf
_f
ield
,
direction
_s
ymbol
,
lb
_m
ethod
,
index
_f
ield
):
"""
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.
:param pdf
F
ield: pystencils field describing the pdf. The current cell is cell next to the boundary,
:param pdf
_f
ield: 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.
:param direction
S
ymbol: a sympy symbol that can be used as index to the pdfField. It describes
:param direction
_s
ymbol: a sympy symbol that can be used as index to the pdfField. It describes
the direction pointing from the fluid to the boundary cell
:param lb
M
ethod: an instance of the LB method used. Use this to adapt the boundary to the method
(e.g. compressiblity)
:param index
F
ield: the boundary index field that can be used to retrieve and update boundary data
:param lb
_m
ethod: an instance of the LB method used. Use this to adapt the boundary to the method
(e.g. compressib
i
lity)
:param index
_f
ield: the boundary index field that can be used to retrieve and update boundary data
:return: list of sympy equations
"""
raise
NotImplementedError
(
"Boundary class has to overwrite __call__"
)
@
property
def
additional
D
ata
(
self
):
def
additional
_d
ata
(
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 additionalDataKernelInit or by
Python callbacks - see additionalDataCallback """
return
[]
@
property
def
additional
D
ata
I
nit
C
allback
(
self
):
def
additional
_d
ata
_i
nit
_c
allback
(
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
...
...
@@ -50,8 +50,8 @@ class Boundary(object):
return
type
(
self
).
__name__
@
name
.
setter
def
name
(
self
,
new
V
alue
):
self
.
_name
=
new
V
alue
def
name
(
self
,
new
_v
alue
):
self
.
_name
=
new
_v
alue
class
NoSlip
(
Boundary
):
...
...
@@ -61,10 +61,10 @@ class NoSlip(Boundary):
super
(
NoSlip
,
self
).
__init__
(
name
)
"""No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle"""
def
__call__
(
self
,
pdf
F
ield
,
direction
S
ymbol
,
lb
M
ethod
,
**
kwargs
):
neighbor
=
BoundaryOffsetInfo
.
offset
F
rom
D
ir
(
direction
S
ymbol
,
lb
M
ethod
.
dim
)
inverse
D
ir
=
BoundaryOffsetInfo
.
inv
D
ir
(
direction
S
ymbol
)
return
[
Assignment
(
pdf
F
ield
[
neighbor
](
inverse
D
ir
),
pdf
F
ield
(
direction
S
ymbol
))]
def
__call__
(
self
,
pdf
_f
ield
,
direction
_s
ymbol
,
lb
_m
ethod
,
**
kwargs
):
neighbor
=
BoundaryOffsetInfo
.
offset
_f
rom
_d
ir
(
direction
_s
ymbol
,
lb
_m
ethod
.
dim
)
inverse
_d
ir
=
BoundaryOffsetInfo
.
inv
_d
ir
(
direction
_s
ymbol
)
return
[
Assignment
(
pdf
_f
ield
[
neighbor
](
inverse
_d
ir
),
pdf
_f
ield
(
direction
_s
ymbol
))]
def
__hash__
(
self
):
# All boundaries of these class behave equal -> should also be equal (as long as name is equal)
...
...
@@ -80,17 +80,17 @@ class UBB(Boundary):
"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""
def
__init__
(
self
,
velocity
,
adapt
V
elocity
ToF
orce
=
False
,
dim
=
None
,
name
=
None
):
def
__init__
(
self
,
velocity
,
adapt
_v
elocity
_to_f
orce
=
False
,
dim
=
None
,
name
=
None
):
"""
:param 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
:param adapt
V
elocity
ToF
orce:
:param adapt
_v
elocity
_to_f
orce:
"""
super
(
UBB
,
self
).
__init__
(
name
)
self
.
_velocity
=
velocity
self
.
_adaptVelocityToForce
=
adapt
V
elocity
ToF
orce
self
.
_adaptVelocityToForce
=
adapt
_v
elocity
_to_f
orce
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
):
...
...
@@ -98,56 +98,57 @@ class UBB(Boundary):
self
.
dim
=
dim
@
property
def
additional
D
ata
(
self
):
def
additional
_d
ata
(
self
):
if
callable
(
self
.
_velocity
):
return
[(
'vel_%d'
%
(
i
,),
create_type
(
"double"
))
for
i
in
range
(
self
.
dim
)]
else
:
return
[]
@
property
def
additional
D
ata
I
nit
C
allback
(
self
):
def
additional
_d
ata
_i
nit
_c
allback
(
self
):
if
callable
(
self
.
_velocity
):
return
self
.
_velocity
def
__call__
(
self
,
pdf
F
ield
,
direction
S
ymbol
,
lb
M
ethod
,
index
F
ield
,
**
kwargs
):
vel
F
rom
IdxF
ield
=
callable
(
self
.
_velocity
)
vel
=
[
index
F
ield
(
'vel_%d'
%
(
i
,))
for
i
in
range
(
self
.
dim
)]
if
vel
F
rom
IdxF
ield
else
self
.
_velocity
direction
=
direction
S
ymbol
def
__call__
(
self
,
pdf
_f
ield
,
direction
_s
ymbol
,
lb
_m
ethod
,
index
_f
ield
,
**
kwargs
):
vel
_f
rom
_idx_f
ield
=
callable
(
self
.
_velocity
)
vel
=
[
index
_f
ield
(
'vel_%d'
%
(
i
,))
for
i
in
range
(
self
.
dim
)]
if
vel
_f
rom
_idx_f
ield
else
self
.
_velocity
direction
=
direction
_s
ymbol
assert
self
.
dim
==
lb
M
ethod
.
dim
,
"Dimension of UBB (%d) does not match dimension of method (%d)"
\
%
(
self
.
dim
,
lb
M
ethod
.
dim
)
assert
self
.
dim
==
lb
_m
ethod
.
dim
,
"Dimension of UBB (%d) does not match dimension of method (%d)"
\
%
(
self
.
dim
,
lb
_m
ethod
.
dim
)
neighbor
=
BoundaryOffsetInfo
.
offset
F
rom
D
ir
(
direction
,
lb
M
ethod
.
dim
)
inverse
D
ir
=
BoundaryOffsetInfo
.
inv
D
ir
(
direction
)
neighbor
=
BoundaryOffsetInfo
.
offset
_f
rom
_d
ir
(
direction
,
lb
_m
ethod
.
dim
)
inverse
_d
ir
=
BoundaryOffsetInfo
.
inv
_d
ir
(
direction
)
velocity
=
tuple
(
v_i
.
get
S
hifted
(
*
neighbor
)
if
isinstance
(
v_i
,
Field
.
Access
)
and
not
vel
F
rom
IdxF
ield
else
v_i
velocity
=
tuple
(
v_i
.
get
_s
hifted
(
*
neighbor
)
if
isinstance
(
v_i
,
Field
.
Access
)
and
not
vel
_f
rom
_idx_f
ield
else
v_i
for
v_i
in
vel
)
if
self
.
_adaptVelocityToForce
:
cqc
=
lb
M
ethod
.
conserved
Q
uantity
C
omputation
shifted
VelE
qs
=
cqc
.
equilibrium
I
nput
E
quations
F
rom
I
nit
V
alues
(
velocity
=
velocity
)
velocity
=
[
eq
.
rhs
for
eq
in
shifted
VelE
qs
.
new_filtered
(
cqc
.
first
O
rder
M
oment
S
ymbols
).
main_assignments
]
cqc
=
lb
_m
ethod
.
conserved
_q
uantity
_c
omputation
shifted
_vel_e
qs
=
cqc
.
equilibrium
_i
nput
_e
quations
_f
rom
_i
nit
_v
alues
(
velocity
=
velocity
)
velocity
=
[
eq
.
rhs
for
eq
in
shifted
_vel_e
qs
.
new_filtered
(
cqc
.
first
_o
rder
_m
oment
_s
ymbols
).
main_assignments
]
c_s_sq
=
sp
.
Rational
(
1
,
3
)
weightOfDirection
=
LbmWeightInfo
.
weightOfDirection
velTerm
=
2
/
c_s_sq
*
sum
([
d_i
*
v_i
for
d_i
,
v_i
in
zip
(
neighbor
,
velocity
)])
*
weightOfDirection
(
direction
)
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
)
# Better alternative: in conserved value computation
# rename what is currently called density to "virtualDensity"
# provide a new quantity density, which is constant in case of incompressible models
if
not
lb
M
ethod
.
conserved
Q
uantity
C
omputation
.
zero
C
entered
P
dfs
:
cqc
=
lb
M
ethod
.
conserved
Q
uantity
C
omputation
density
S
ymbol
=
sp
.
Symbol
(
"rho"
)
pdf
F
ield
A
ccesses
=
[
pdf
F
ield
(
i
)
for
i
in
range
(
len
(
lb
M
ethod
.
stencil
))]
density
E
quations
=
cqc
.
output
E
quations
F
rom
P
dfs
(
pdf
F
ield
A
ccesses
,
{
'density'
:
density
S
ymbol
})
density
S
ymbol
=
lb
M
ethod
.
conserved
Q
uantity
C
omputation
.
defined_symbols
()[
'density'
]
result
=
density
E
quations
.
all_assignments
result
+=
[
Assignment
(
pdf
F
ield
[
neighbor
](
inverse
D
ir
),
pdf
F
ield
(
direction
)
-
vel
T
erm
*
density
S
ymbol
)]
if
not
lb
_m
ethod
.
conserved
_q
uantity
_c
omputation
.
zero
_c
entered
_p
dfs
:
cqc
=
lb
_m
ethod
.
conserved
_q
uantity
_c
omputation
density
_s
ymbol
=
sp
.
Symbol
(
"rho"
)
pdf
_f
ield
_a
ccesses
=
[
pdf
_f
ield
(
i
)
for
i
in
range
(
len
(
lb
_m
ethod
.
stencil
))]
density
_e
quations
=
cqc
.
output
_e
quations
_f
rom
_p
dfs
(
pdf
_f
ield
_a
ccesses
,
{
'density'
:
density
_s
ymbol
})
density
_s
ymbol
=
lb
_m
ethod
.
conserved
_q
uantity
_c
omputation
.
defined_symbols
()[
'density'
]
result
=
density
_e
quations
.
all_assignments
result
+=
[
Assignment
(
pdf
_f
ield
[
neighbor
](
inverse
_d
ir
),
pdf
_f
ield
(
direction
)
-
vel
_t
erm
*
density
_s
ymbol
)]
return
result
else
:
return
[
Assignment
(
pdf
F
ield
[
neighbor
](
inverse
D
ir
),
pdf
F
ield
(
direction
)
-
vel
T
erm
)]
return
[
Assignment
(
pdf
_f
ield
[
neighbor
](
inverse
_d
ir
),
pdf
_f
ield
(
direction
)
-
vel
_t
erm
)]
class
FixedDensity
(
Boundary
):
...
...
@@ -156,49 +157,52 @@ class FixedDensity(Boundary):
super
(
FixedDensity
,
self
).
__init__
(
name
)
self
.
_density
=
density
def
__call__
(
self
,
pdf
F
ield
,
direction
S
ymbol
,
lb
M
ethod
,
**
kwargs
):
def
__call__
(
self
,
pdf
_f
ield
,
direction
_s
ymbol
,
lb
_m
ethod
,
**
kwargs
):
"""Boundary condition that fixes the density/pressure at the obstacle"""
def
removeAsymmetricPartOfmain_assignments
(
eqColl
,
dofs
):
newmain_assignments
=
[
Assignment
(
e
.
lhs
,
get_symmetric_part
(
e
.
rhs
,
dofs
))
for
e
in
eqColl
.
main_assignments
]
return
eqColl
.
copy
(
newmain_assignments
)
def
remove_asymmetric_part_of_main_assignments
(
assignment_collection
,
degrees_of_freedom
):
new_main_assignments
=
[
Assignment
(
a
.
lhs
,
get_symmetric_part
(
a
.
rhs
,
degrees_of_freedom
))
for
a
in
assignment_collection
.
main_assignments
]
return
assignment_collection
.
copy
(
new_main_assignments
)
neighbor
=
BoundaryOffsetInfo
.
offset
F
rom
D
ir
(
direction
S
ymbol
,
lb
M
ethod
.
dim
)
inverse
D
ir
=
BoundaryOffsetInfo
.
inv
D
ir
(
direction
S
ymbol
)
neighbor
=
BoundaryOffsetInfo
.
offset
_f
rom
_d
ir
(
direction
_s
ymbol
,
lb
_m
ethod
.
dim
)
inverse
_d
ir
=
BoundaryOffsetInfo
.
inv
_d
ir
(
direction
_s
ymbol
)
cqc
=
lb
M
ethod
.
conserved
Q
uantity
C
omputation
cqc
=
lb
_m
ethod
.
conserved
_q
uantity
_c
omputation
velocity
=
cqc
.
defined_symbols
()[
'velocity'
]
symmetricEq
=
removeAsymmetricPartOfmain_assignments
(
lbMethod
.
getEquilibrium
(),
dofs
=
velocity
)
substitutions
=
{
sym
:
pdfField
(
i
)
for
i
,
sym
in
enumerate
(
lbMethod
.
preCollisionPdfSymbols
)}
symmetricEq
=
symmetricEq
.
new_with_substitutions
(
substitutions
)
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
)}
symmetric_eq
=
symmetric_eq
.
new_with_substitutions
(
substitutions
)
simplification
=
create
S
implification
S
trategy
(
lb
M
ethod
)
symmetric
E
q
=
simplification
(
symmetric
E
q
)
simplification
=
create
_s
implification
_s
trategy
(
lb
_m
ethod
)
symmetric
_e
q
=
simplification
(
symmetric
_e
q
)
density
S
ymbol
=
cqc
.
defined_symbols
()[
'density'
]
density
_s
ymbol
=
cqc
.
defined_symbols
()[
'density'
]
density
=
self
.
_density
densityEq
=
cqc
.
equilibriumInputEquationsFromInitValues
(
density
=
density
).
new_without_subexpressions
().
main_assignments
[
0
]
assert
densityEq
.
lhs
==
densitySymbol
transformedDensity
=
densityEq
.
rhs
equilibrium_input
=
cqc
.
equilibrium_input_equations_from_init_values
(
density
=
density
).
new_without_subexpressions
()
density_eq
=
equilibrium_input
.
main_assignments
[
0
]
assert
density_eq
.
lhs
==
density_symbol
transformed_density
=
density_eq
.
rhs
conditions
=
[(
eq_i
.
rhs
,
sp
.
Equality
(
direction
S
ymbol
,
i
))
for
i
,
eq_i
in
enumerate
(
symmetric
E
q
.
main_assignments
)]
+
[(
0
,
True
)]
conditions
=
[(
eq_i
.
rhs
,
sp
.
Equality
(
direction
_s
ymbol
,
i
))
for
i
,
eq_i
in
enumerate
(
symmetric
_e
q
.
main_assignments
)]
+
[(
0
,
True
)]
eq_component
=
sp
.
Piecewise
(
*
conditions
)
sub
E
xprs
=
[
Assignment
(
eq
.
lhs
,
transformed
D
ensity
if
eq
.
lhs
==
density
S
ymbol
else
eq
.
rhs
)
for
eq
in
symmetric
E
q
.
subexpressions
]
sub
e
xpr
ession
s
=
[
Assignment
(
eq
.
lhs
,
transformed
_d
ensity
if
eq
.
lhs
==
density
_s
ymbol
else
eq
.
rhs
)
for
eq
in
symmetric
_e
q
.
subexpressions
]
return
sub
E
xprs
+
[
SympyAssignment
(
pdf
F
ield
[
neighbor
](
inverse
D
ir
),
2
*
eq_component
-
pdf
F
ield
(
direction
S
ymbol
))]
return
sub
e
xpr
ession
s
+
[
SympyAssignment
(
pdf
_f
ield
[
neighbor
](
inverse
_d
ir
),
2
*
eq_component
-
pdf
_f
ield
(
direction
_s
ymbol
))]
class
NeumannByCopy
(
Boundary
):
def
__call__
(
self
,
pdf
F
ield
,
direction
S
ymbol
,
lb
M
ethod
,
**
kwargs
):
neighbor
=
BoundaryOffsetInfo
.
offset
F
rom
D
ir
(
direction
S
ymbol
,
lb
M
ethod
.
dim
)
inverse
D
ir
=
BoundaryOffsetInfo
.
inv
D
ir
(
direction
S
ymbol
)
return
[
Assignment
(
pdf
F
ield
[
neighbor
](
inverse
D
ir
),
pdf
F
ield
(
inverse
D
ir
)),
Assignment
(
pdf
F
ield
[
neighbor
](
direction
S
ymbol
),
pdf
F
ield
(
direction
S
ymbol
))]
def
__call__
(
self
,
pdf
_f
ield
,
direction
_s
ymbol
,
lb
_m
ethod
,
**
kwargs
):
neighbor
=
BoundaryOffsetInfo
.
offset
_f
rom
_d
ir
(
direction
_s
ymbol
,
lb
_m
ethod
.
dim
)
inverse
_d
ir
=
BoundaryOffsetInfo
.
inv
_d
ir
(
direction
_s
ymbol
)
return
[
Assignment
(
pdf
_f
ield
[
neighbor
](
inverse
_d
ir
),
pdf
_f
ield
(
inverse
_d
ir
)),
Assignment
(
pdf
_f
ield
[
neighbor
](
direction
_s
ymbol
),
pdf
_f
ield
(
direction
_s
ymbol
))]
def
__hash__
(
self
):
# All boundaries of these class behave equal -> should also be equal
...
...
@@ -209,11 +213,11 @@ class NeumannByCopy(Boundary):
class
StreamInZero
(
Boundary
):
def
__call__
(
self
,
pdf
F
ield
,
direction
S
ymbol
,
lb
M
ethod
,
**
kwargs
):
neighbor
=
BoundaryOffsetInfo
.
offset
F
rom
D
ir
(
direction
S
ymbol
,
lb
M
ethod
.
dim
)
inverse
D
ir
=
BoundaryOffsetInfo
.
inv
D
ir
(
direction
S
ymbol
)
return
[
Assignment
(
pdf
F
ield
[
neighbor
](
inverse
D
ir
),
0
),
Assignment
(
pdf
F
ield
[
neighbor
](
direction
S
ymbol
),
0
)]
def
__call__
(
self
,
pdf
_f
ield
,
direction
_s
ymbol
,
lb
_m
ethod
,
**
kwargs
):
neighbor
=
BoundaryOffsetInfo
.
offset
_f
rom
_d
ir
(
direction
_s
ymbol
,
lb
_m
ethod
.
dim
)
inverse
_d
ir
=
BoundaryOffsetInfo
.
inv
_d
ir
(
direction
_s
ymbol
)
return
[
Assignment
(
pdf
_f
ield
[
neighbor
](
inverse
_d
ir
),
0
),
Assignment
(
pdf
_f
ield
[
neighbor
](
direction
_s
ymbol
),
0
)]
def
__hash__
(
self
):
# All boundaries of these class behave equal -> should also be equal
...
...
@@ -226,12 +230,12 @@ class StreamInZero(Boundary):
class
AntiBounceBack
(
Boundary
):
"""No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle"""
def
__call__
(
self
,
pdf
F
ield
,
direction
S
ymbol
,
lb
M
ethod
,
**
kwargs
):
neighbor
=
BoundaryOffsetInfo
.
offset
F
rom
D
ir
(
direction
S
ymbol
,
lb
M
ethod
.
dim
)
inverse
D
ir
=
BoundaryOffsetInfo
.
inv
D
ir
(
direction
S
ymbol
)
rhs
=
pdf
F
ield
(
direction
S
ymbol
)
def
__call__
(
self
,
pdf
_f
ield
,
direction
_s
ymbol
,
lb
_m
ethod
,
**
kwargs
):
neighbor
=
BoundaryOffsetInfo
.
offset
_f
rom
_d
ir
(
direction
_s
ymbol
,
lb
_m
ethod
.
dim
)
inverse
_d
ir
=
BoundaryOffsetInfo
.
inv
_d
ir
(
direction
_s
ymbol
)
rhs
=
pdf
_f
ield
(
direction
_s
ymbol
)
t
=
-
rhs
return
[
SympyAssignment
(
pdf
F
ield
[
neighbor
](
inverse
D
ir
),
t
)]
return
[
SympyAssignment
(
pdf
_f
ield
[
neighbor
](
inverse
_d
ir
),
t
)]
def
__hash__
(
self
):
# All boundaries of these class behave equal -> should also be equal (as long as name is equal)
...
...
@@ -241,4 +245,3 @@ class AntiBounceBack(Boundary):
if
not
isinstance
(
other
,
AntiBounceBack
):
return
False
return
self
.
name
==
other
.
name
boundaries/boundaryhandling.py
View file @
482c8293
import
numpy
as
np
import
sympy
as
sp
from
pystencils
import
TypedSymbol
,
create
I
ndexed
K
ernel
,
Assignment
from
pystencils
import
TypedSymbol
,
create
_i
ndexed
_k
ernel
,
Assignment
from
pystencils.backends.cbackend
import
CustomCppCode
from
pystencils.boundaries
import
BoundaryHandling
from
pystencils.boundaries.boundaryhandling
import
BoundaryOffsetInfo
from
lbmpy.stencils
import
inverse
D
irection
from
lbmpy.stencils
import
inverse
_d
irection
class
LatticeBoltzmannBoundaryHandling
(
BoundaryHandling
):
def
__init__
(
self
,
lb
M
ethod
,
data
H
andling
,
pdf
F
ield
N
ame
,
name
=
"boundary
H
andling"
,
flag
I
nterface
=
None
,
target
=
'cpu'
,
open
MP
=
True
):
self
.
lb
M
ethod
=
lb
M
ethod
super
(
LatticeBoltzmannBoundaryHandling
,
self
).
__init__
(
data
H
andling
,
pdf
F
ield
N
ame
,
lb
M
ethod
.
stencil
,
name
,
flag
I
nterface
,
target
,
open
MP
)
def
__init__
(
self
,
lb
_m
ethod
,
data
_h
andling
,
pdf
_f
ield
_n
ame
,
name
=
"boundary
_h
andling"
,
flag
_i
nterface
=
None
,
target
=
'cpu'
,
open
mp
=
True
):
self
.
lb
_m
ethod
=
lb
_m
ethod
super
(
LatticeBoltzmannBoundaryHandling
,
self
).
__init__
(
data
_h
andling
,
pdf
_f
ield
_n
ame
,
lb
_m
ethod
.
stencil
,
name
,
flag
_i
nterface
,
target
,
open
mp
)
def
force
OnB
oundary
(
self
,
boundary
O
bj
):
def
force
_on_b
oundary
(
self
,
boundary
_o
bj
):
from
lbmpy.boundaries
import
NoSlip
if
isinstance
(
boundary
O
bj
,
NoSlip
):
return
self
.
_force
OnNoS
lip
(
boundary
O
bj
)
if
isinstance
(
boundary
_o
bj
,
NoSlip
):
return
self
.
_force
_on_no_s
lip
(
boundary
_o
bj
)
else
:
self
.
__call__
()
return
self
.
_force
OnB
oundary
(
boundary
O
bj
)
return
self
.
_force
_on_b
oundary
(
boundary
_o
bj
)
# ------------------------------ Implementation Details ------------------------------------------------------------
def
_force
OnNoS
lip
(
self
,
boundary
O
bj
):
dh
=
self
.
_data
H
andling
ff
G
host
L
ayers
=
dh
.
ghost
L
ayers
OfF
ield
(
self
.
flag
I
nterface
.
flag
F
ield
N
ame
)
method
=
self
.
lb
M
ethod
def
_force
_on_no_s
lip
(
self
,
boundary
_o
bj
):
dh
=
self
.
_data
_h
andling
ff
_g
host
_l
ayers
=
dh
.
ghost
_l
ayers
_of_f
ield
(
self
.
flag
_i
nterface
.
flag
_f
ield
_n
ame
)
method
=
self
.
lb
_m
ethod
stencil
=
np
.
array
(
method
.
stencil
)
result
=
np
.
zeros
(
self
.
dim
)
for
b
in
dh
.
iterate
(
ghost
L
ayers
=
ff
G
host
L
ayers
):
obj
ToIndL
ist
=
b
[
self
.
_index
A
rray
N
ame
].
boundary
O
bject
ToI
ndex
L
ist
pdf
A
rray
=
b
[
self
.
_field
N
ame
]
if
boundary
O
bj
in
obj
ToIndL
ist
:
ind
A
rr
=
obj
ToIndL
ist
[
boundary
O
bj
]
indices
=
[
ind
A
rr
[
name
]
for
name
in
(
'x'
,
'y'
,
'z'
)[:
method
.
dim
]]
indices
.
append
(
ind
A
rr
[
'dir'
])
values
=
2
*
pdf
A
rray
[
tuple
(
indices
)]
forces
=
stencil
[
ind
A
rr
[
'dir'
]]
*
values
[:,
np
.
newaxis
]
for
b
in
dh
.
iterate
(
ghost
_l
ayers
=
ff
_g
host
_l
ayers
):
obj
_to_ind_l
ist
=
b
[
self
.
_index
_a
rray
_n
ame
].
boundary
_o
bject
_to_i
ndex
_l
ist
pdf
_a
rray
=
b
[
self
.
_field
_n
ame
]
if
boundary
_o
bj
in
obj
_to_ind_l
ist
:
ind
_a
rr
=
obj
_to_ind_l
ist
[
boundary
_o
bj
]
indices
=
[
ind
_a
rr
[
name
]
for
name
in
(
'x'
,
'y'
,
'z'
)[:
method
.
dim
]]
indices
.
append
(
ind
_a
rr
[
'dir'
])
values
=
2
*
pdf
_a
rray
[
tuple
(
indices
)]
forces
=
stencil
[
ind
_a
rr
[
'dir'
]]
*
values
[:,
np
.
newaxis
]
result
+=
forces
.
sum
(
axis
=
0
)
return
dh
.
reduce
F
loat
S
equence
(
list
(
result
),
'sum'
)
return
dh
.
reduce
_f
loat
_s
equence
(
list
(
result
),
'sum'
)
def
_force
OnB
oundary
(
self
,
boundary
O
bj
):
dh
=
self
.
_data
H
andling
ff
G
host
L
ayers
=
dh
.
ghost
L
ayers
OfF
ield
(
self
.
flag
I
nterface
.
flag
F
ield
N
ame
)
method
=
self
.
lb
M
ethod
def
_force
_on_b
oundary
(
self
,
boundary
_o
bj
):
dh
=
self
.
_data
_h
andling
ff
_g
host
_l
ayers
=
dh
.
ghost
_l
ayers
_of_f
ield
(
self
.
flag
_i
nterface
.
flag
_f
ield
_n
ame
)
method
=
self
.
lb
_m
ethod
stencil
=
np
.
array
(
method
.
stencil
)
inv
D
irection
=
np
.
array
([
method
.
stencil
.
index
(
inverse
D
irection
(
d
))
inv
_d
irection
=
np
.
array
([
method
.
stencil
.
index
(
inverse
_d
irection
(
d
))
for
d
in
method
.
stencil
])
result
=
np
.
zeros
(
self
.
dim
)
for
b
in
dh
.
iterate
(
ghost
L
ayers
=
ff
G
host
L
ayers
):
obj
ToIndL
ist
=
b
[
self
.
_index
A
rray
N
ame
].
boundary
O
bject
ToI
ndex
L
ist
pdf
A
rray
=
b
[
self
.
_field
N
ame
]
if
boundary
O
bj
in
obj
ToIndL
ist
:
ind
A
rr
=
obj
ToIndL
ist
[
boundary
O
bj
]
indices
=
[
ind
A
rr
[
name
]
for
name
in
(
'x'
,
'y'
,
'z'
)[:
method
.
dim
]]
offsets
=
stencil
[
ind
A
rr
[
'dir'
]]
inv
D
ir
=
inv
D
irection
[
ind
A
rr
[
'dir'
]]
fluid
V
alues
=
pdf
A
rray
[
tuple
(
indices
)
+
(
ind
A
rr
[
'dir'
],)]
boundary
V
alues
=
pdf
A
rray
[
tuple
(
ind
+
offsets
[:,
i
]
for
i
,
ind
in
enumerate
(
indices
))
+
(
inv
D
ir
,)]
values
=
fluid
V
alues
+
boundary
V
alues
forces
=
stencil
[
ind
A
rr
[
'dir'
]]
*
values
[:,
np
.
newaxis
]
for
b
in
dh
.
iterate
(
ghost
_l
ayers
=
ff
_g
host
_l
ayers
):
obj
_to_ind_l
ist
=
b
[
self
.
_index
_a
rray
_n
ame
].
boundary
_o
bject
_to_i
ndex
_l
ist
pdf
_a
rray
=
b
[
self
.
_field
_n
ame
]
if
boundary
_o
bj
in
obj
_to_ind_l
ist
:
ind
_a
rr
=
obj
_to_ind_l
ist
[
boundary
_o
bj
]
indices
=
[
ind
_a
rr
[
name
]
for
name
in
(
'x'
,
'y'
,
'z'
)[:
method
.
dim
]]
offsets
=
stencil
[
ind
_a
rr
[
'dir'
]]
inv
_d
ir
=
inv
_d
irection
[
ind
_a
rr
[
'dir'
]]
fluid
_v
alues
=
pdf
_a
rray
[
tuple
(
indices
)
+
(
ind
_a
rr
[
'dir'
],)]
boundary
_v
alues
=
pdf
_a
rray
[
tuple
(
ind
+
offsets
[:,
i
]
for
i
,
ind
in
enumerate
(
indices
))
+
(
inv
_d
ir
,)]
values
=
fluid
_v
alues
+
boundary
_v
alues
forces
=
stencil
[
ind
_a
rr
[
'dir'
]]
*
values
[:,
np
.
newaxis
]
result
+=
forces
.
sum
(
axis
=
0
)
return
dh
.
reduce
F
loat
S
equence
(
list
(
result
),
'sum'
)
return
dh
.
reduce
_f
loat
_s
equence
(
list
(
result
),
'sum'
)
def
_create
B
oundary
K
ernel
(
self
,
symbolic
F
ield
,
symbolic
I
ndex
F
ield
,
boundary
Object
):
return
create
L
attice
B
oltzmann
B
oundary
K
ernel
(
symbolic
F
ield
,
symbolic
I
ndex
F
ield
,
self
.
lb
M
ethod
,
boundary
Object
,
target
=
self
.
_target
,
open
MP
=
self
.
_open
MP
)
def
_create
_b
oundary
_k
ernel
(
self
,
symbolic
_f
ield
,
symbolic
_i
ndex
_f
ield
,
boundary
_obj
):
return
create
_l
attice
_b
oltzmann
_b
oundary
_k
ernel
(
symbolic
_f
ield
,
symbolic
_i
ndex
_f
ield
,
self
.
lb
_m
ethod
,
boundary
_obj
,
target
=
self
.
_target
,
open
mp
=
self
.
_open
mp
)
class
LbmWeightInfo
(
CustomCppCode
):
...
...
@@ -81,24 +81,26 @@ class LbmWeightInfo(CustomCppCode):
# --------------------------- Functions to be used by boundaries --------------------------
@
staticmethod
def
weight
OfD
irection
(
dir
I
dx
):
return
sp
.
IndexedBase
(
LbmWeightInfo
.
WEIGHTS_SYMBOL
,
shape
=
(
1
,))[
dir
I
dx
]
def
weight
_of_d
irection
(
dir
_i
dx
):
return
sp
.
IndexedBase
(
LbmWeightInfo
.
WEIGHTS_SYMBOL
,
shape
=
(
1
,))[
dir
_i
dx
]
# ---------------------------------- Internal ---------------------------------------------
WEIGHTS_SYMBOL
=
TypedSymbol
(
"weights"
,
"double"
)
def
__init__
(
self
,
lb
M
ethod
):
weights
=
[
str
(
w
.
evalf
())
for
w
in
lb
M
ethod
.
weights
]
def
__init__
(
self
,
lb
_m
ethod
):
weights
=
[
str
(
w
.
evalf
())
for
w
in
lb
_m
ethod
.
weights
]
code
=
"const double %s [] = { %s };
\n
"
%
(
LbmWeightInfo
.
WEIGHTS_SYMBOL
.
name
,
","
.
join
(
weights
))
super
(
LbmWeightInfo
,
self
).
__init__
(
code
,
symbols_read
=
set
(),
symbols_defined
=
set
([
LbmWeightInfo
.
WEIGHTS_SYMBOL
]))
def
createLatticeBoltzmannBoundaryKernel
(
pdfField
,
indexField
,
lbMethod
,
boundaryFunctor
,
target
=
'cpu'
,
openMP
=
True
):
elements
=
[
BoundaryOffsetInfo
(
lbMethod
.
stencil
),
LbmWeightInfo
(
lbMethod
)]
indexArrDtype
=
indexField
.
dtype
.
numpy_dtype
dirSymbol
=
TypedSymbol
(
"dir"
,
indexArrDtype
.
fields
[
'dir'
][
0
])
elements
+=
[
Assignment
(
dirSymbol
,
indexField
[
0
](
'dir'
))]
elements
+=
boundaryFunctor
(
pdfField
=
pdfField
,
directionSymbol
=
dirSymbol
,
lbMethod
=
lbMethod
,
indexField
=
indexField
)
return
createIndexedKernel
(
elements
,
[
indexField
],
target
=
target
,
cpuOpenMP
=
openMP
)
symbols_defined
=
{
LbmWeightInfo
.
WEIGHTS_SYMBOL
})
def
create_lattice_boltzmann_boundary_kernel
(
pdf_field
,
index_field
,
lb_method
,
boundary_functor
,
target
=
'cpu'
,
openmp
=
True
):
elements
=
[
BoundaryOffsetInfo
(
lb_method
.
stencil
),
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
)
boundaries/createindexlistcython.pyx
View file @
482c8293
...
...
@@ -13,51 +13,51 @@ ctypedef fused IntegerType:
unsigned
int
unsigned
long
@
cython
.
boundscheck
(
False
)
# turn off bounds-checking for entire function
@
cython
.
wraparound
(
False
)
# turn off negative index wrapping for entire function
def
create
B
oundary
I
ndex
L
ist
2D
(
object
[
IntegerType
,
ndim
=
2
]
flag
F
ield
,
int
nr
OfG
host
L
ayers
,
IntegerType
boundary
M
ask
,
IntegerType
fluid
M
ask
,
object
[
int
,
ndim
=
2
]
stencil
):
@
cython
.
boundscheck
(
False
)
@
cython
.
wraparound
(
False
)
def
create
_b
oundary
_i
ndex
_l
ist
_2d
(
object
[
IntegerType
,
ndim
=
2
]
flag
_f
ield
,
int
nr
_of_g
host
_l
ayers
,
IntegerType
boundary
_m
ask
,
IntegerType
fluid
_m
ask
,
object
[
int
,
ndim
=
2
]
stencil
):
cdef
int
xs
,
ys
,
x
,
y
cdef
int
dirIdx
,
num
D
irections
,
dx
,
dy
cdef
int
dirIdx
,
num
_d
irections
,
dx
,
dy
xs
,
ys
=
flag
F
ield
.
shape
boundary
I
ndex
L
ist
=
[]
num
D
irections
=
stencil
.
shape
[
0
]
xs
,
ys
=
flag
_f
ield
.
shape
boundary
_i
ndex
_l
ist
=
[]
num
_d
irections
=
stencil
.
shape
[
0
]
for
y
in
range
(
nr
OfG
host
L
ayers
,
ys
-
nrOfG
host
L
ayers
):
for
x
in
range
(
nr
OfG
host
L
ayers
,
xs
-
nrOfG
host
L
ayers
):
if
flag
F
ield
[
x
,
y
]
&
fluid
M
ask
:
for
dirIdx
in
range
(
1
,
num
D
irections
):
for
y
in
range
(
nr
_of_g
host
_l
ayers
,
ys
-
nr_of_g
host
_l
ayers
):
for
x
in
range
(
nr
_of_g
host
_l
ayers
,
xs
-
nr_of_g
host
_l
ayers
):
if
flag
_f
ield
[
x
,
y
]
&
fluid
_m
ask
:
for
dirIdx
in
range
(
1
,
num
_d
irections
):
dx
=
stencil
[
dirIdx
,
0
]
dy
=
stencil
[
dirIdx
,
1
]
if
flag
F
ield
[
x
+
dx
,
y
+
dy
]
&
boundary
M
ask
:
boundary
I
ndex
L
ist
.
append
((
x
,
y
,
dirIdx
))
return
boundary
I
ndex
L
ist
if
flag
_f
ield
[
x
+
dx
,
y
+
dy
]
&
boundary
_m
ask
:
boundary
_i
ndex
_l
ist
.
append
((
x
,
y
,
dirIdx
))
return
boundary
_i
ndex
_l
ist
@
cython
.
boundscheck
(
False
)
# turn off bounds-checking for entire function
@
cython
.
wraparound
(
False
)
# turn off negative index wrapping for entire function
def
create
B
oundary
I
ndex
L
ist
3D
(
object
[
IntegerType
,
ndim
=
3
]
flag
F
ield
,
int
nr
OfG
host
L
ayers
,
IntegerType
boundary
M
ask
,
IntegerType
fluid
M
ask
,
object
[
int
,
ndim
=
2
]
stencil
):
@
cython
.
boundscheck
(
False
)
@
cython
.
wraparound
(
False
)
def
create
_b
oundary
_i
ndex
_l
ist
_3d
(
object
[
IntegerType
,
ndim
=
3
]
flag
_f
ield
,
int
nr
_of_g
host
_l
ayers
,
IntegerType
boundary
_m
ask
,
IntegerType
fluid
_m
ask
,
object
[
int
,
ndim
=
2
]
stencil
):
cdef
int
xs
,
ys
,
zs
,
x
,
y
,
z
cdef
int
dirIdx
,
num
D
irections
,
dx
,
dy
,
dz
cdef
int
dirIdx
,
num
_d
irections
,
dx
,
dy
,
dz
xs
,
ys
,
zs
=
flag
F
ield
.
shape
boundary
I
ndex
L
ist
=
[]
num
D
irections
=
stencil
.
shape
[
0
]
xs
,
ys
,
zs
=
flag
_f
ield
.
shape
boundary
_i
ndex
_l
ist
=
[]
num
_d
irections
=
stencil
.
shape
[
0
]
for
z
in
range
(
nr
OfG
host
L
ayers
,
zs
-
nrOfG
host
L
ayers
):
for
y
in
range
(
nr
OfG
host
L
ayers
,
ys
-
nrOfG
host
L
ayers
):
for
x
in
range
(
nr
OfG
host
L
ayers
,
xs
-
nrOfG
host
L
ayers
):
if
flag
F
ield
[
x
,
y
,
z
]
&
fluid
M
ask
:
for
dirIdx
in
range
(
1
,
num
D
irections
):
for
z
in
range
(
nr
_of_g
host
_l
ayers
,
zs
-
nr_of_g
host
_l
ayers
):
for
y
in
range
(
nr
_of_g
host
_l
ayers
,
ys
-
nr_of_g
host
_l
ayers
):
for
x
in
range
(
nr
_of_g
host
_l
ayers
,
xs
-
nr_of_g
host
_l
ayers
):
if
flag
_f
ield
[
x
,
y
,
z
]
&
fluid
_m
ask
:
for
dirIdx
in
range
(
1
,
num
_d
irections
):
dx
=
stencil
[
dirIdx
,
0
]
dy
=
stencil
[
dirIdx
,
1
]
dz
=
stencil
[
dirIdx
,
2
]
if
flag
F
ield
[
x
+
dx
,
y
+
dy
,
z
+
dz
]
&
boundary
M
ask
:
boundary
I
ndex
L
ist
.
append
((
x
,
y
,
z
,
dirIdx
))
return
boundary
I
ndex
L
ist
if
flag
_f
ield
[
x
+
dx
,
y
+
dy
,
z
+
dz
]
&
boundary
_m
ask
:
boundary
_i
ndex
_l
ist
.
append
((
x
,
y
,
z
,
dirIdx
))
return
boundary
_i
ndex
_l
ist
chapman_enskog/__init__.py
View file @
482c8293
from
lbmpy.chapman_enskog.derivative
import
DiffOperator
,
Diff
,
expand
U
sing
L
inearity
,
expand
U
sing
P
roduct
R
ule
,
\
normalize
D
iff
O
rder
,
chapman
E
nskog
D
erivative
E
xpansion
,
chapman
E
nskog
D
erivative
R
ecombination
from
lbmpy.chapman_enskog.chapman_enskog
import
LbMethodEqMoments
,
insert
M
oments
,
take
M
oments
,
\
CeMoment
,
chain
S
olve
AndS
ubstitute
,
time
D
iff
S
elector
,
moment
S
elector
,
ChapmanEnskogAnalysis