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
Sebastian Bindgen
pystencils
Commits
d72cd721
Commit
d72cd721
authored
Apr 03, 2018
by
Martin Bauer
Browse files
PEP8 name refactoring
- test run again - notebooks not yet
parent
3bcfac93
Changes
52
Expand all
Hide whitespace changes
Inline
Side-by-side
__init__.py
View file @
d72cd721
from
pystencils.field
import
Field
,
FieldType
,
extractCommonSubexpressions
from
pystencils.field
import
Field
,
FieldType
from
pystencils.data_types
import
TypedSymbol
from
pystencils.slicing
import
makeSlice
from
pystencils.kernelcreation
import
create
K
ernel
,
create
I
ndexed
K
ernel
from
pystencils.kernelcreation
import
create
_k
ernel
,
create
_i
ndexed
_k
ernel
from
pystencils.display_utils
import
show_code
,
to_dot
from
pystencils.assignment_collection
import
AssignmentCollection
from
pystencils.assignment
import
Assignment
from
pystencils.sympyextensions
import
SymbolCreator
__all__
=
[
'Field'
,
'FieldType'
,
'extractCommonSubexpressions'
,
__all__
=
[
'Field'
,
'FieldType'
,
'TypedSymbol'
,
'makeSlice'
,
'create
K
ernel'
,
'create
I
ndexed
K
ernel'
,
'create
_k
ernel'
,
'create
_i
ndexed
_k
ernel'
,
'show_code'
,
'to_dot'
,
'AssignmentCollection'
,
'Assignment'
,
...
...
alignedarray.py
View file @
d72cd721
import
numpy
as
np
def
aligned_empty
(
shape
,
byte
A
lignment
=
32
,
dtype
=
np
.
float64
,
byte
O
ffset
=
0
,
order
=
'C'
,
align
I
nner
C
oordinate
=
True
):
def
aligned_empty
(
shape
,
byte
_a
lignment
=
32
,
dtype
=
np
.
float64
,
byte
_o
ffset
=
0
,
order
=
'C'
,
align
_i
nner
_c
oordinate
=
True
):
"""
Creates an aligned empty numpy array
:param shape: size of the array
:param byte
A
lignment: alignment in bytes, for the start address of the array holds (a % byteAlignment) == 0
:param byte
_a
lignment: alignment in bytes, for the start address of the array holds (a % byteAlignment) == 0
:param dtype: numpy data type
:param byte
O
ffset: offset in bytes for position that should be aligned i.e. (a+byte
O
ffset) % byteAlignment == 0
:param byte
_o
ffset: offset in bytes for position that should be aligned i.e. (a+byte
_o
ffset) % byteAlignment == 0
typically used to align first inner cell instead of ghost layer
:param order: storage linearization order
:param align
I
nner
C
oordinate: if True, the start of the innermost coordinate lines are aligned as well
:param align
_i
nner
_c
oordinate: if True, the start of the innermost coordinate lines are aligned as well
:return:
"""
if
(
not
align
I
nner
C
oordinate
)
or
(
not
hasattr
(
shape
,
'__len__'
)):
N
=
np
.
prod
(
shape
)
if
(
not
align
_i
nner
_c
oordinate
)
or
(
not
hasattr
(
shape
,
'__len__'
)):
size
=
np
.
prod
(
shape
)
d
=
np
.
dtype
(
dtype
)
tmp
=
np
.
empty
(
N
*
d
.
itemsize
+
byte
A
lignment
,
dtype
=
np
.
uint8
)
tmp
=
np
.
empty
(
size
*
d
.
itemsize
+
byte
_a
lignment
,
dtype
=
np
.
uint8
)
address
=
tmp
.
__array_interface__
[
'data'
][
0
]
offset
=
(
byteAlignment
-
(
address
+
byteOffset
)
%
byteAlignment
)
%
byteAlignment
t1
=
tmp
[
offset
:
offset
+
N
*
d
.
itemsize
]
return
tmp
[
offset
:
offset
+
N
*
d
.
itemsize
].
view
(
dtype
=
d
).
reshape
(
shape
,
order
=
order
)
offset
=
(
byte_alignment
-
(
address
+
byte_offset
)
%
byte_alignment
)
%
byte_alignment
return
tmp
[
offset
:
offset
+
size
*
d
.
itemsize
].
view
(
dtype
=
d
).
reshape
(
shape
,
order
=
order
)
else
:
if
order
==
'C'
:
n
dim0
=
shape
[
-
1
]
dim0
_size
=
shape
[
-
1
]
dim0
=
-
1
n
dim1
=
np
.
prod
(
shape
[:
-
1
])
dim1
_size
=
np
.
prod
(
shape
[:
-
1
])
else
:
n
dim0
=
shape
[
0
]
dim0
_size
=
shape
[
0
]
dim0
=
0
n
dim1
=
np
.
prod
(
shape
[
1
:])
dim1
_size
=
np
.
prod
(
shape
[
1
:])
d
=
np
.
dtype
(
dtype
)
assert
byte
A
lignment
>=
d
.
itemsize
and
byte
A
lignment
%
d
.
itemsize
==
0
padding
=
(
byte
A
lignment
-
((
n
dim0
*
d
.
itemsize
)
%
byte
A
lignment
))
%
byte
A
lignment
assert
byte
_a
lignment
>=
d
.
itemsize
and
byte
_a
lignment
%
d
.
itemsize
==
0
padding
=
(
byte
_a
lignment
-
((
dim0
_size
*
d
.
itemsize
)
%
byte
_a
lignment
))
%
byte
_a
lignment
N
=
ndim1
*
padding
+
np
.
prod
(
shape
)
*
d
.
itemsize
tmp
=
aligned_empty
(
N
,
byteAlignment
=
byteAlignment
,
dtype
=
np
.
uint8
,
byteOffset
=
byteOffset
).
view
(
dtype
=
dtype
)
bshape
=
[
i
for
i
in
shape
]
bshape
[
dim0
]
=
ndim0
+
padding
//
d
.
itemsize
tmp
=
tmp
.
reshape
(
bshape
,
order
=
order
)
size
=
dim1_size
*
padding
+
np
.
prod
(
shape
)
*
d
.
itemsize
tmp
=
aligned_empty
(
size
,
byte_alignment
=
byte_alignment
,
dtype
=
np
.
uint8
,
byte_offset
=
byte_offset
)
tmp
=
tmp
.
view
(
dtype
=
dtype
)
shape_in_bytes
=
[
i
for
i
in
shape
]
shape_in_bytes
[
dim0
]
=
dim0_size
+
padding
//
d
.
itemsize
tmp
=
tmp
.
reshape
(
shape_in_bytes
,
order
=
order
)
if
tmp
.
flags
[
'C_CONTIGUOUS'
]:
tmp
=
tmp
[...,
:
shape
[
-
1
]]
else
:
...
...
@@ -48,17 +48,17 @@ def aligned_empty(shape, byteAlignment=32, dtype=np.float64, byteOffset=0, order
return
tmp
def
aligned_zeros
(
shape
,
byte
A
lignment
=
16
,
dtype
=
float
,
byte
O
ffset
=
0
,
order
=
'C'
,
align
I
nner
C
oordinate
=
True
):
arr
=
aligned_empty
(
shape
,
dtype
=
dtype
,
byte
O
ffset
=
byte
O
ffset
,
order
=
order
,
byte
A
lignment
=
byte
A
lignment
,
align
I
nner
C
oordinate
=
align
I
nner
C
oordinate
)
def
aligned_zeros
(
shape
,
byte
_a
lignment
=
16
,
dtype
=
float
,
byte
_o
ffset
=
0
,
order
=
'C'
,
align
_i
nner
_c
oordinate
=
True
):
arr
=
aligned_empty
(
shape
,
dtype
=
dtype
,
byte
_o
ffset
=
byte
_o
ffset
,
order
=
order
,
byte
_a
lignment
=
byte
_a
lignment
,
align
_i
nner
_c
oordinate
=
align
_i
nner
_c
oordinate
)
x
=
np
.
zeros
((),
arr
.
dtype
)
arr
[...]
=
x
return
arr
def
aligned_ones
(
shape
,
byte
A
lignment
=
16
,
dtype
=
float
,
byte
O
ffset
=
0
,
order
=
'C'
,
align
I
nner
C
oordinate
=
True
):
arr
=
aligned_empty
(
shape
,
dtype
=
dtype
,
byte
O
ffset
=
byte
O
ffset
,
order
=
order
,
byte
A
lignment
=
byte
A
lignment
,
align
I
nner
C
oordinate
=
align
I
nner
C
oordinate
)
def
aligned_ones
(
shape
,
byte
_a
lignment
=
16
,
dtype
=
float
,
byte
_o
ffset
=
0
,
order
=
'C'
,
align
_i
nner
_c
oordinate
=
True
):
arr
=
aligned_empty
(
shape
,
dtype
=
dtype
,
byte
_o
ffset
=
byte
_o
ffset
,
order
=
order
,
byte
_a
lignment
=
byte
_a
lignment
,
align
_i
nner
_c
oordinate
=
align
_i
nner
_c
oordinate
)
x
=
np
.
ones
((),
arr
.
dtype
)
arr
[...]
=
x
return
arr
astnodes.py
View file @
d72cd721
import
sympy
as
sp
from
sympy.tensor
import
IndexedBase
from
pystencils.field
import
Field
from
pystencils.data_types
import
TypedSymbol
,
create_type
,
cast
F
unc
from
pystencils.data_types
import
TypedSymbol
,
create_type
,
cast
_f
unc
from
pystencils.sympyextensions
import
fast_subs
from
typing
import
List
,
Set
,
Optional
,
Union
,
Any
...
...
@@ -113,34 +113,34 @@ class KernelFunction(Node):
class
Argument
:
def
__init__
(
self
,
name
,
dtype
,
symbol
,
kernel_function_node
):
from
pystencils.transformations
import
symbol
N
ame
ToV
ariable
N
ame
from
pystencils.transformations
import
symbol
_n
ame
_to_v
ariable
_n
ame
self
.
name
=
name
self
.
dtype
=
dtype
self
.
isFieldPtrArgument
=
False
self
.
isFieldShapeArgument
=
False
self
.
isFieldStrideArgument
=
False
self
.
isFieldArgument
=
False
self
.
field
N
ame
=
""
self
.
field
_n
ame
=
""
self
.
coordinate
=
None
self
.
symbol
=
symbol
if
name
.
startswith
(
Field
.
DATA_PREFIX
):
self
.
isFieldPtrArgument
=
True
self
.
isFieldArgument
=
True
self
.
field
N
ame
=
name
[
len
(
Field
.
DATA_PREFIX
):]
self
.
field
_n
ame
=
name
[
len
(
Field
.
DATA_PREFIX
):]
elif
name
.
startswith
(
Field
.
SHAPE_PREFIX
):
self
.
isFieldShapeArgument
=
True
self
.
isFieldArgument
=
True
self
.
field
N
ame
=
name
[
len
(
Field
.
SHAPE_PREFIX
):]
self
.
field
_n
ame
=
name
[
len
(
Field
.
SHAPE_PREFIX
):]
elif
name
.
startswith
(
Field
.
STRIDE_PREFIX
):
self
.
isFieldStrideArgument
=
True
self
.
isFieldArgument
=
True
self
.
field
N
ame
=
name
[
len
(
Field
.
STRIDE_PREFIX
):]
self
.
field
_n
ame
=
name
[
len
(
Field
.
STRIDE_PREFIX
):]
self
.
field
=
None
if
self
.
isFieldArgument
:
field_map
=
{
symbol
N
ame
ToV
ariable
N
ame
(
f
.
name
):
f
for
f
in
kernel_function_node
.
fields_accessed
}
self
.
field
=
field_map
[
self
.
field
N
ame
]
field_map
=
{
symbol
_n
ame
_to_v
ariable
_n
ame
(
f
.
name
):
f
for
f
in
kernel_function_node
.
fields_accessed
}
self
.
field
=
field_map
[
self
.
field
_n
ame
]
def
__lt__
(
self
,
other
):
def
score
(
l
):
...
...
@@ -167,12 +167,12 @@ class KernelFunction(Node):
self
.
_body
=
body
body
.
parent
=
self
self
.
_parameters
=
None
self
.
function
N
ame
=
function_name
self
.
function
_n
ame
=
function_name
self
.
_body
.
parent
=
self
self
.
compile
=
None
self
.
ghost
L
ayers
=
ghost_layers
self
.
ghost
_l
ayers
=
ghost_layers
# these variables are assumed to be global, so no automatic parameter is generated for them
self
.
global
V
ariables
=
set
()
self
.
global
_v
ariables
=
set
()
self
.
backend
=
backend
@
property
...
...
@@ -202,19 +202,19 @@ class KernelFunction(Node):
return
set
(
o
.
field
for
o
in
self
.
atoms
(
ResolvedFieldAccess
))
def
_update_parameters
(
self
):
undefined_symbols
=
self
.
_body
.
undefined_symbols
-
self
.
global
V
ariables
undefined_symbols
=
self
.
_body
.
undefined_symbols
-
self
.
global
_v
ariables
self
.
_parameters
=
[
KernelFunction
.
Argument
(
s
.
name
,
s
.
dtype
,
s
,
self
)
for
s
in
undefined_symbols
]
self
.
_parameters
.
sort
()
def
__str__
(
self
):
self
.
_update_parameters
()
return
'{0} {1}({2})
\n
{3}'
.
format
(
type
(
self
).
__name__
,
self
.
function
N
ame
,
self
.
parameters
,
return
'{0} {1}({2})
\n
{3}'
.
format
(
type
(
self
).
__name__
,
self
.
function
_n
ame
,
self
.
parameters
,
(
"
\t
"
+
"
\t
"
.
join
(
str
(
self
.
body
).
splitlines
(
True
))))
def
__repr__
(
self
):
self
.
_update_parameters
()
return
'{0} {1}({2})'
.
format
(
type
(
self
).
__name__
,
self
.
function
N
ame
,
self
.
parameters
)
return
'{0} {1}({2})'
.
format
(
type
(
self
).
__name__
,
self
.
function
_n
ame
,
self
.
parameters
)
class
Block
(
Node
):
...
...
@@ -392,8 +392,8 @@ class LoopOverCoordinate(Node):
@
property
def
is_outermost_loop
(
self
):
from
pystencils.transformations
import
get
N
ext
P
arent
OfT
ype
return
get
N
ext
P
arent
OfT
ype
(
self
,
LoopOverCoordinate
)
is
None
from
pystencils.transformations
import
get
_n
ext
_p
arent
_of_t
ype
return
get
_n
ext
_p
arent
_of_t
ype
(
self
,
LoopOverCoordinate
)
is
None
@
property
def
is_innermost_loop
(
self
):
...
...
@@ -417,7 +417,7 @@ class SympyAssignment(Node):
self
.
_lhsSymbol
=
lhs_symbol
self
.
rhs
=
rhs_expr
self
.
_isDeclaration
=
True
is_cast
=
self
.
_lhsSymbol
.
func
==
cast
F
unc
is_cast
=
self
.
_lhsSymbol
.
func
==
cast
_f
unc
if
isinstance
(
self
.
_lhsSymbol
,
Field
.
Access
)
or
isinstance
(
self
.
_lhsSymbol
,
ResolvedFieldAccess
)
or
is_cast
:
self
.
_isDeclaration
=
False
self
.
_isConst
=
is_const
...
...
@@ -430,7 +430,7 @@ class SympyAssignment(Node):
def
lhs
(
self
,
new_value
):
self
.
_lhsSymbol
=
new_value
self
.
_isDeclaration
=
True
is_cast
=
self
.
_lhsSymbol
.
func
==
cast
F
unc
is_cast
=
self
.
_lhsSymbol
.
func
==
cast
_f
unc
if
isinstance
(
self
.
_lhsSymbol
,
Field
.
Access
)
or
isinstance
(
self
.
_lhsSymbol
,
sp
.
Indexed
)
or
is_cast
:
self
.
_isDeclaration
=
False
...
...
backends/__init__.py
View file @
d72cd721
from
.cbackend
import
print
_c
from
.cbackend
import
generate
_c
try
:
from
.dot
import
print_dot
...
...
backends/cbackend.py
View file @
d72cd721
import
sympy
as
sp
from
collections
import
namedtuple
from
sympy.core
import
S
from
typing
import
Optional
from
typing
import
Optional
,
Set
try
:
from
sympy.printing.ccode
import
C99CodePrinter
as
CCodePrinter
except
ImportError
:
from
sympy.printing.ccode
import
CCodePrinter
# for sympy versions < 1.1
from
pystencils.bitoperations
import
xor
,
rightShift
,
leftShi
ft
,
bitwise
A
nd
,
bitwise
O
r
from
pystencils.bitoperations
import
bitwise_xor
,
bit_shift_right
,
bit_shift_le
ft
,
bitwise
_a
nd
,
bitwise
_o
r
from
pystencils.astnodes
import
Node
,
ResolvedFieldAccess
,
SympyAssignment
from
pystencils.data_types
import
create_type
,
PointerType
,
get_type_of_expression
,
VectorType
,
cast
F
unc
from
pystencils.data_types
import
create_type
,
PointerType
,
get_type_of_expression
,
VectorType
,
cast
_f
unc
from
pystencils.backends.simd_instruction_sets
import
selectedInstructionSet
__all__
=
[
'
print_c
'
]
__all__
=
[
'
generate_c'
,
'CustomCppCode'
,
'get_headers
'
]
def
print
_c
(
ast_node
:
Node
,
signature_only
:
bool
=
False
,
use_float_constants
:
Optional
[
bool
]
=
None
)
->
str
:
def
generate
_c
(
ast_node
:
Node
,
signature_only
:
bool
=
False
,
use_float_constants
:
Optional
[
bool
]
=
None
)
->
str
:
"""Prints an abstract syntax tree node as C or CUDA code.
This function does not need to distinguish between C, C++ or CUDA code, it just prints 'C-like' code as encoded
...
...
@@ -42,7 +42,8 @@ def print_c(ast_node: Node, signature_only: bool = False, use_float_constants: O
return
printer
(
ast_node
)
def
get_headers
(
ast_node
):
def
get_headers
(
ast_node
:
Node
)
->
Set
[
str
]:
"""Return a set of header files, necessary to compile the printed C-like code."""
headers
=
set
()
if
hasattr
(
ast_node
,
'headers'
):
...
...
@@ -131,7 +132,7 @@ class CBackend:
def
_print_KernelFunction
(
self
,
node
):
function_arguments
=
[
"%s %s"
%
(
str
(
s
.
dtype
),
s
.
name
)
for
s
in
node
.
parameters
]
func_declaration
=
"FUNC_PREFIX void %s(%s)"
%
(
node
.
function
N
ame
,
", "
.
join
(
function_arguments
))
func_declaration
=
"FUNC_PREFIX void %s(%s)"
%
(
node
.
function
_n
ame
,
", "
.
join
(
function_arguments
))
if
self
.
_signatureOnly
:
return
func_declaration
...
...
@@ -163,7 +164,7 @@ class CBackend:
return
"%s %s = %s;"
%
(
data_type
,
self
.
sympyPrinter
.
doprint
(
node
.
lhs
),
self
.
sympyPrinter
.
doprint
(
node
.
rhs
))
else
:
lhs_type
=
get_type_of_expression
(
node
.
lhs
)
if
type
(
lhs_type
)
is
VectorType
and
node
.
lhs
.
func
==
cast
F
unc
:
if
type
(
lhs_type
)
is
VectorType
and
node
.
lhs
.
func
==
cast
_f
unc
:
return
self
.
_vectorInstructionSet
[
'storeU'
].
format
(
"&"
+
self
.
sympyPrinter
.
doprint
(
node
.
lhs
.
args
[
0
]),
self
.
sympyPrinter
.
doprint
(
node
.
rhs
))
+
';'
else
:
...
...
@@ -231,13 +232,13 @@ class CustomSympyPrinter(CCodePrinter):
def
_print_Function
(
self
,
expr
):
function_map
=
{
xor
:
'^'
,
rightShif
t
:
'>>'
,
leftShi
ft
:
'<<'
,
bitwise
O
r
:
'|'
,
bitwise
A
nd
:
'&'
,
bitwise_
xor
:
'^'
,
bit_shift_righ
t
:
'>>'
,
bit_shift_le
ft
:
'<<'
,
bitwise
_o
r
:
'|'
,
bitwise
_a
nd
:
'&'
,
}
if
expr
.
func
==
cast
F
unc
:
if
expr
.
func
==
cast
_f
unc
:
arg
,
data_type
=
expr
.
args
return
"*((%s)(& %s))"
%
(
PointerType
(
data_type
),
self
.
_print
(
arg
))
elif
expr
.
func
in
function_map
:
...
...
@@ -263,7 +264,7 @@ class VectorizedCustomSympyPrinter(CustomSympyPrinter):
return
None
def
_print_Function
(
self
,
expr
):
if
expr
.
func
==
cast
F
unc
:
if
expr
.
func
==
cast
_f
unc
:
arg
,
data_type
=
expr
.
args
if
type
(
data_type
)
is
VectorType
:
if
type
(
arg
)
is
ResolvedFieldAccess
:
...
...
backends/dot.py
View file @
d72cd721
...
...
@@ -72,7 +72,7 @@ def __shortened(node):
elif
isinstance
(
node
,
KernelFunction
):
params
=
[
f
.
name
for
f
in
node
.
fields_accessed
]
params
+=
[
p
.
name
for
p
in
node
.
parameters
if
not
p
.
isFieldArgument
]
return
"Func: %s (%s)"
%
(
node
.
function
N
ame
,
","
.
join
(
params
))
return
"Func: %s (%s)"
%
(
node
.
function
_n
ame
,
","
.
join
(
params
))
elif
isinstance
(
node
,
SympyAssignment
):
return
repr
(
node
.
lhs
)
elif
isinstance
(
node
,
Block
):
...
...
bitoperations.py
View file @
d72cd721
import
sympy
as
sp
xor
=
sp
.
Function
(
"⊻"
)
rightShif
t
=
sp
.
Function
(
"rshift"
)
leftShi
ft
=
sp
.
Function
(
"lshift"
)
bitwise
A
nd
=
sp
.
Function
(
"Bit&"
)
bitwise
O
r
=
sp
.
Function
(
"Bit|"
)
bitwise_
xor
=
sp
.
Function
(
"⊻"
)
bit_shift_righ
t
=
sp
.
Function
(
"rshift"
)
bit_shift_le
ft
=
sp
.
Function
(
"lshift"
)
bitwise
_a
nd
=
sp
.
Function
(
"Bit&"
)
bitwise
_o
r
=
sp
.
Function
(
"Bit|"
)
boundaries/boundaryconditions.py
View file @
d72cd721
from
pystencils
import
Assignment
from
pystencils.boundaries.boundaryhandling
import
BoundaryOffsetInfo
from
typing
import
List
,
Tuple
,
Any
class
Boundary
(
object
):
...
...
@@ -8,29 +9,30 @@ class Boundary(object):
def
__init__
(
self
,
name
=
None
):
self
.
_name
=
name
def
__call__
(
self
,
field
,
directionSymbol
,
indexField
):
"""
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 field: pystencils field where boundary condition should be applied.
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 directionSymbol: 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 indexField: the boundary index field that can be used to retrieve and update boundary data
:return: list of sympy equations
def
__call__
(
self
,
field
,
direction_symbol
,
index_field
)
->
List
[
Assignment
]:
"""Defines the boundary behavior and must therefore be implemented by all boundaries.
Here the boundary is defined as a list of sympy assignments, from which a boundary kernel is generated.
Args:
field: pystencils field where boundary condition should be applied.
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 pdfField. It describes
the direction pointing from the fluid to the boundary cell
index_field: the boundary index field that can be used to retrieve and update boundary data
"""
raise
NotImplementedError
(
"Boundary class has to overwrite __call__"
)
@
property
def
additional
D
ata
(
self
):
def
additional
_d
ata
(
self
)
->
Tuple
[
str
,
Any
]
:
"""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
...
...
@@ -43,22 +45,22 @@ 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
Neumann
(
Boundary
):
def
__call__
(
self
,
field
,
direction
S
ymbol
,
**
kwargs
):
def
__call__
(
self
,
field
,
direction
_s
ymbol
,
**
kwargs
):
neighbor
=
BoundaryOffsetInfo
.
offset
F
rom
D
ir
(
direction
S
ymbol
,
field
.
spatial
D
imensions
)
if
field
.
index
D
imensions
==
0
:
neighbor
=
BoundaryOffsetInfo
.
offset
_f
rom
_d
ir
(
direction
_s
ymbol
,
field
.
spatial
_d
imensions
)
if
field
.
index
_d
imensions
==
0
:
return
[
Assignment
(
field
[
neighbor
],
field
.
center
)]
else
:
from
itertools
import
product
if
not
field
.
has
F
ixed
I
ndex
S
hape
:
if
not
field
.
has
_f
ixed
_i
ndex
_s
hape
:
raise
NotImplementedError
(
"Neumann boundary works only for fields with fixed index shape"
)
index
I
ter
=
product
(
*
(
range
(
i
)
for
i
in
field
.
index
S
hape
))
return
[
Assignment
(
field
[
neighbor
](
*
idx
),
field
(
*
idx
))
for
idx
in
index
I
ter
]
index
_i
ter
=
product
(
*
(
range
(
i
)
for
i
in
field
.
index
_s
hape
))
return
[
Assignment
(
field
[
neighbor
](
*
idx
),
field
(
*
idx
))
for
idx
in
index
_i
ter
]
def
__hash__
(
self
):
# All boundaries of these class behave equal -> should also be equal
...
...
boundaries/boundaryhandling.py
View file @
d72cd721
This diff is collapsed.
Click to expand it.
boundaries/createindexlist.py
View file @
d72cd721
...
...
@@ -6,7 +6,7 @@ try:
import
pyximport
;
pyximport
.
install
()
from
pystencils.boundaries.createindexlistcython
import
create
B
oundary
I
ndex
L
ist
2D
,
create
B
oundary
I
ndex
L
ist
3D
from
pystencils.boundaries.createindexlistcython
import
create
_b
oundary
_i
ndex
_l
ist
_2d
,
create
_b
oundary
_i
ndex
_l
ist
_3d
cythonFuncsAvailable
=
True
except
Exception
:
...
...
@@ -22,7 +22,7 @@ def numpyDataTypeForBoundaryObject(boundaryObject, dim):
coordinateNames
=
boundaryIndexArrayCoordinateNames
[:
dim
]
return
np
.
dtype
([(
name
,
np
.
int32
)
for
name
in
coordinateNames
]
+
[(
directionMemberName
,
np
.
int32
)]
+
[(
i
[
0
],
i
[
1
].
numpy_dtype
)
for
i
in
boundaryObject
.
additional
D
ata
],
align
=
True
)
[(
i
[
0
],
i
[
1
].
numpy_dtype
)
for
i
in
boundaryObject
.
additional
_d
ata
],
align
=
True
)
def
_createBoundaryIndexListPython
(
flagFieldArr
,
nrOfGhostLayers
,
boundaryMask
,
fluidMask
,
stencil
):
...
...
@@ -51,9 +51,9 @@ def createBoundaryIndexList(flagField, stencil, boundaryMask, fluidMask, nrOfGho
if
cythonFuncsAvailable
:
stencil
=
np
.
array
(
stencil
,
dtype
=
np
.
int32
)
if
dim
==
2
:
idxList
=
create
B
oundary
I
ndex
L
ist
2D
(
flagField
,
nrOfGhostLayers
,
boundaryMask
,
fluidMask
,
stencil
)
idxList
=
create
_b
oundary
_i
ndex
_l
ist
_2d
(
flagField
,
nrOfGhostLayers
,
boundaryMask
,
fluidMask
,
stencil
)
elif
dim
==
3
:
idxList
=
create
B
oundary
I
ndex
L
ist
3D
(
flagField
,
nrOfGhostLayers
,
boundaryMask
,
fluidMask
,
stencil
)
idxList
=
create
_b
oundary
_i
ndex
_l
ist
_3d
(
flagField
,
nrOfGhostLayers
,
boundaryMask
,
fluidMask
,
stencil
)
else
:
raise
ValueError
(
"Flag field has to be a 2 or 3 dimensional numpy array"
)
return
np
.
array
(
idxList
,
dtype
=
indexArrDtype
)
...
...
@@ -67,7 +67,7 @@ def createBoundaryIndexArray(flagField, stencil, boundaryMask, fluidMask, bounda
idxArray
=
createBoundaryIndexList
(
flagField
,
stencil
,
boundaryMask
,
fluidMask
,
nrOfGhostLayers
)
dim
=
len
(
flagField
.
shape
)
if
boundaryObject
.
additional
D
ata
:
if
boundaryObject
.
additional
_d
ata
:
coordinateNames
=
boundaryIndexArrayCoordinateNames
[:
dim
]
indexArrDtype
=
numpyDataTypeForBoundaryObject
(
boundaryObject
,
dim
)
extendedIdxField
=
np
.
empty
(
len
(
idxArray
),
dtype
=
indexArrDtype
)
...
...
boundaries/createindexlistcython.pyx
View file @
d72cd721
...
...
@@ -15,49 +15,49 @@ ctypedef fused IntegerType:
@
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
):
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
]
flagField
,
int
nrOfGhostLayers
,
IntegerType
boundaryMask
,
IntegerType
fluidMask
,
object
[
int
,
ndim
=
2
]
stencil
):
def
create
_b
oundary
_i
ndex
_l
ist
_3d
(
object
[
IntegerType
,
ndim
=
3
]
flagField
,
int
nrOfGhostLayers
,
IntegerType
boundaryMask
,
IntegerType
fluidMask
,
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
=
flagField
.
shape
boundary
I
ndex
L
ist
=
[]
num
D
irections
=
stencil
.
shape
[
0
]
boundary
_i
ndex
_l
ist
=
[]
num
_d
irections
=
stencil
.
shape
[
0
]
for
z
in
range
(
nrOfGhostLayers
,
zs
-
nrOfGhostLayers
):
for
y
in
range
(
nrOfGhostLayers
,
ys
-
nrOfGhostLayers
):
for
x
in
range
(
nrOfGhostLayers
,
xs
-
nrOfGhostLayers
):
if
flagField
[
x
,
y
,
z
]
&
fluidMask
:
for
dirIdx
in
range
(
1
,
num
D
irections
):
for
dirIdx
in
range
(
1
,
num
_d
irections
):
dx
=
stencil
[
dirIdx
,
0
]
dy
=
stencil
[
dirIdx
,
1
]
dz
=
stencil
[
dirIdx
,
2
]
if
flagField
[
x
+
dx
,
y
+
dy
,
z
+
dz
]
&
boundaryMask
:
boundary
I
ndex
L
ist
.
append
((
x
,
y
,
z
,
dirIdx
))
return
boundary
I
ndex
L
ist
boundary
_i
ndex
_l
ist
.
append
((
x
,
y
,
z
,
dirIdx
))
return
boundary
_i
ndex
_l
ist
boundaries/inkernel.py
View file @
d72cd721
import
sympy
as
sp
from
pystencils
import
Field
,
TypedSymbol
from
pystencils.bitoperations
import
bitwise
A
nd
from
pystencils.bitoperations
import
bitwise
_a
nd
from
pystencils.boundaries.boundaryhandling
import
FlagInterface
from
pystencils.data_types
import
create_type
def
add
N
eumann
B
oundary
(
eqs
,
fields
,
flag
F
ield
,
boundary
F
lag
=
"neumannFlag"
,
inverse
F
lag
=
False
):
def
add
_n
eumann
_b
oundary
(
eqs
,
fields
,
flag
_f
ield
,
boundary
_f
lag
=
"neumannFlag"
,
inverse
_f
lag
=
False
):
"""
Replaces all neighbor accesses by flag field guarded accesses.
If flag in neighboring cell is set, the center value is used instead
:param eqs: list of equations containing field accesses to direct neighbors
:param fields: fields for which the Neumann boundary should be applied
:param flagField: integer field marking boundary cells
:param boundaryFlag: if flag field has value 'boundaryFlag' (no bitoperations yet) the cell is assumed to be boundary
:param inverseFlag: if true, boundary cells are where flagfield has not the value of boundaryFlag
:param flag_field: integer field marking boundary cells
:param boundary_flag: if flag field has value 'boundaryFlag' (no bit operations yet)
the cell is assumed to be boundary
:param inverse_flag: if true, boundary cells are where flag field has not the value of boundaryFlag
:return: list of equations with guarded field accesses
"""
if
not
hasattr
(
fields
,
"__len__"
):
fields
=
[
fields
]
fields
=
set
(
fields
)
if
type
(
boundary
F
lag
)
is
str
:
boundary
F
lag
=
TypedSymbol
(
boundary
F
lag
,
dtype
=
create_type
(
FlagInterface
.
FLAG_DTYPE
))
if
type
(
boundary
_f
lag
)
is
str
:
boundary
_f
lag
=
TypedSymbol
(
boundary
_f
lag
,
dtype
=
create_type
(
FlagInterface
.
FLAG_DTYPE
))
substitutions
=
{}
for
eq
in
eqs
:
...
...
@@ -33,10 +34,10 @@ def addNeumannBoundary(eqs, fields, flagField, boundaryFlag="neumannFlag", inver
if
all
(
offset
==
0
for
offset
in
fa
.
offsets
):
continue
if
inverse
F
lag
:
condition
=
sp
.
Eq
(
bitwise
A
nd
(
flag
F
ield
[
tuple
(
fa
.
offsets
)],
boundary
F
lag
),
0
)
if
inverse
_f
lag
:
condition
=
sp
.
Eq
(
bitwise
_a
nd
(
flag
_f
ield
[
tuple
(
fa
.
offsets
)],
boundary
_f
lag
),
0
)
else
:
condition
=
sp
.
Ne
(
bitwise
A
nd
(
flag
F
ield
[
tuple
(
fa
.
offsets
)],
boundary
F
lag
),
0
)
condition
=
sp
.
Ne
(
bitwise
_a
nd
(
flag
_f
ield
[
tuple
(
fa
.
offsets
)],
boundary
_f
lag
),
0
)