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
Tom Harke
pystencils
Commits
4651c191
Commit
4651c191
authored
Apr 26, 2018
by
Martin Bauer
Browse files
Documentation overhaul Part2
parent
01ab38e8
Changes
9
Hide whitespace changes
Inline
Side-by-side
__init__.py
View file @
4651c191
...
...
@@ -3,7 +3,7 @@ from . import sympy_gmpy_bug_workaround # NOQA
from
.field
import
Field
,
FieldType
,
fields
from
.data_types
import
TypedSymbol
from
.slicing
import
make_slice
from
.kernelcreation
import
create_kernel
,
create_indexed_kernel
from
.kernelcreation
import
create_kernel
,
create_indexed_kernel
,
create_staggered_kernel
from
.display_utils
import
show_code
,
to_dot
from
.simp
import
AssignmentCollection
from
.assignment
import
Assignment
...
...
@@ -15,7 +15,7 @@ from . import fd
__all__
=
[
'Field'
,
'FieldType'
,
'fields'
,
'TypedSymbol'
,
'make_slice'
,
'create_kernel'
,
'create_indexed_kernel'
,
'create_kernel'
,
'create_indexed_kernel'
,
'create_staggered_kernel'
,
'show_code'
,
'to_dot'
,
'AssignmentCollection'
,
'Assignment'
,
...
...
@@ -23,4 +23,3 @@ __all__ = ['Field', 'FieldType', 'fields',
'create_data_handling'
,
'kernel'
,
'fd'
]
cpu/cpujit.py
View file @
4651c191
r
"""
*pystencils* automatically searches for a compiler, so in most cases no explicit configuration is required.
On Linux make sure that 'gcc' and 'g++' are installed and in your path.
On Windows a recent Visual Studio installation is required.
In case anything does not work as expected or a special compiler should be used, changes can be specified
in a configuration file.
*pystencils* looks for a configuration file in JSON format at the following locations in the listed order.
1. at the path specified in the environment variable ``PYSTENCILS_CONFIG``
...
...
display_utils.py
View file @
4651c191
...
...
@@ -33,7 +33,7 @@ def highlight_cpp(code: str):
def
show_code
(
ast
:
KernelFunction
):
"""Returns an object to display
C code.
"""Returns an object to display
generated code (C/C++ or CUDA)
Can either be displayed as HTML in Jupyter notebooks or printed as normal string.
"""
...
...
field.py
View file @
4651c191
...
...
@@ -9,38 +9,40 @@ from pystencils.alignedarray import aligned_empty
from
pystencils.data_types
import
TypedSymbol
,
create_type
,
create_composite_type_from_string
,
StructType
from
pystencils.sympyextensions
import
is_integer_sequence
__all__
=
[
'Field'
,
'fields'
,
'FieldType'
]
def
fields
(
description
=
None
,
index_dimensions
=
0
,
layout
=
None
,
**
kwargs
):
"""Creates pystencils fields from a string description.
Examples:
Create a 2D scalar and vector field
>>> s, v = fields("s, v(2): double[2D]")
>>> assert s.spatial_dimensions == 2 and s.index_dimensions == 0
>>> assert v.spatial_dimensions == 2 and v.index_dimensions == 1 and v.index_shape == (2,)
Create an integer field of shape (10, 20)
>>> f = fields("f : int32[10, 20]")
>>> f.has_fixed_shape, f.shape
(True, (10, 20))
Create a 2D scalar and vector field:
>>> s, v = fields("s, v(2): double[2D]")
>>> assert s.spatial_dimensions == 2 and s.index_dimensions == 0
>>> assert (v.spatial_dimensions, v.index_dimensions, v.index_shape) == (2, 1, (2,))
Numpy arrays can be used as template for shape and data type of field
>>> arr_s, arr_v = np.zeros([20, 20]), np.zeros([2
0, 20
, 2
])
>>> s, v = fields("s, v(2)", s=arr_s, v=arr_v)
>>> assert s.index_dimensions == 0 and v.index_shape == (2,) and s.dtype.numpy_dtype == arr_s.dtype
Create an integer field of shape (10, 20):
>>> f = fields("f : int32[1
0, 20]
"
)
>>> f.has_fixed_shape, f.shape
(True, (10, 20))
Format string can be left out, field names are taken from keyword arguments.
>>> fields(f1=arr_s, f2=arr_s)
[f1, f2]
Numpy arrays can be used as template for shape and data type of field:
>>> arr_s, arr_v = np.zeros([20, 20]), np.zeros([20, 20, 2])
>>> s, v = fields("s, v(2)", s=arr_s, v=arr_v)
>>> assert s.index_dimensions == 0 and s.dtype.numpy_dtype == arr_s.dtype
>>> assert v.index_shape == (2,)
The keyword names 'index_dimension' and 'layout' have special meaning and thus can not be used to pass
numpy arrays:
>>> f = fields(f=arr_v, index_dimensions=1)
>>> assert f.index_dimensions == 1
>>> f = fields("pdfs(19) : float32[3D]", layout='fzyx')
>>> f.layout
(2, 1, 0)
Format string can be left out, field names are taken from keyword arguments.
>>> fields(f1=arr_s, f2=arr_s)
[f1, f2]
The keyword names ``index_dimension`` and ``layout`` have special meaning, don't use them for field names
>>> f = fields(f=arr_v, index_dimensions=1)
>>> assert f.index_dimensions == 1
>>> f = fields("pdfs(19) : float32[3D]", layout='fzyx')
>>> f.layout
(2, 1, 0)
"""
result
=
[]
if
description
:
...
...
@@ -108,36 +110,41 @@ class Field:
This Field class knows about the dimension, memory layout (strides) and optionally about the size of an array.
Creating Fields:
To create a field use one of the static create* members. There are two options:
1. create a kernel with fixed loop sizes i.e. the shape of the array is already known. This is usually the
case if just-in-time compilation directly from Python is done. (see :func:`Field.create_from_numpy_array`)
The preferred method to create fields is the `fields` function.
Alternatively one can use one of the static functions `Field.create_generic`, `Field.create_from_numpy_array`
and `Field.create_fixed_size`. Don't instantiate the Field directly!
Fields can be created with known or unknown shapes:
1. If you want to create a kernel with fixed loop sizes i.e. the shape of the array is already known.
This is usually the case if just-in-time compilation directly from Python is done.
(see `Field.create_from_numpy_array`
2. create a more general kernel that works for variable array sizes. This can be used to create kernels
beforehand for a library. (see
:func:
`Field.create_generic`)
beforehand for a library. (see `Field.create_generic`)
Dimensions:
Dimensions
and Indexing
:
A field has spatial and index dimensions, where the spatial dimensions come first.
The interpretation is that the field has multiple cells in (usually) two or three dimensional space which are
looped over. Additionally
N values are stored per cell. In this case spatial_dimensions is two or three,
looped over. Additionally N values are stored per cell. In this case spatial_dimensions is two or three,
and index_dimensions equals N. If you want to store a matrix on each point in a two dimensional grid, there
are four dimensions, two spatial and two index dimensions: ``len(arr.shape) == spatial_dims + index_dims``
Indexing:
When accessing (indexing) a field the result is a FieldAccess which is derived from sympy Symbol.
The shape of the index dimension does not have to be specified. Just use the 'index_dimensions' parameter.
However, it is good practice to define the shape, since out of bounds accesses can be directly detected in this
case. The shape can be passed with the 'index_shape' parameter of the field creation functions.
When accessing (indexing) a field the result is a `Field.Access` which is derived from sympy Symbol.
First specify the spatial offsets in [], then in case index_dimension>0 the indices in ()
e.g. ``f[-1,0,0](7)``
Example
without
index dimensions:
Example
using no
index dimensions:
>>> a = np.zeros([10, 10])
>>> f = Field.create_from_numpy_array("f", a, index_dimensions=0)
>>> jacobi = (
f[-1,0] + f[1,0] + f[0,-1] + f[0,1]
) / 4
>>> jacobi = (f[-1,0] + f[1,0] + f[0,-1] + f[0,1]) / 4
Example
with
index dimensions
: LBM D2Q9
stream pull
Example
s for
index dimensions
to create LB field and implement
stream pull
:
>>> from pystencils import Assignment
>>> stencil = np.array([[0,0], [0,1], [0,-1]])
>>> src = Field.create_generic("src", spatial_dimensions=2, index_dimensions=1)
>>> dst = Field.create_generic("dst", spatial_dimensions=2, index_dimensions=1)
>>> src, dst = fields("src(3), dst(3) : double[2D]")
>>> for i, offset in enumerate(stencil):
... Assignment(dst[0,0](i), src[-offset](i))
Assignment(dst_C^0, src_C^0)
...
...
@@ -380,6 +387,24 @@ class Field:
# noinspection PyAttributeOutsideInit,PyUnresolvedReferences
class
Access
(
sp
.
Symbol
):
"""Class representing a relative access into a `Field`.
This class behaves like a normal sympy Symbol, it is actually derived from it. One can built up
sympy expressions using field accesses, solve for them, etc.
Examples:
>>> vector_field_2d = fields("v(2): double[2D]") # create a 2D vector field
>>> northern_neighbor_y_component = vector_field_2d[0, 1](1)
>>> northern_neighbor_y_component
v_N^1
>>> central_y_component = vector_field_2d(1)
>>> central_y_component
v_C^1
>>> central_y_component.get_shifted(1, 0) # move the existing access
v_E^1
>>> central_y_component.at_index(0) # change component
v_C^0
"""
def
__new__
(
cls
,
name
,
*
args
,
**
kwargs
):
obj
=
Field
.
Access
.
__xnew_cached_
(
cls
,
name
,
*
args
,
**
kwargs
)
return
obj
...
...
@@ -458,15 +483,18 @@ class Field:
raise
TypeError
(
"Field access is not iterable"
)
@
property
def
field
(
self
):
def
field
(
self
)
->
'Field'
:
"""Field that the Access points to"""
return
self
.
_field
@
property
def
offsets
(
self
):
def
offsets
(
self
)
->
Tuple
:
"""Spatial offset as tuple"""
return
tuple
(
self
.
_offsets
)
@
property
def
required_ghost_layers
(
self
):
def
required_ghost_layers
(
self
)
->
int
:
"""Largest spatial distance that is accessed."""
return
int
(
np
.
max
(
np
.
abs
(
self
.
_offsets
)))
@
property
...
...
@@ -475,21 +503,54 @@ class Field:
@
property
def
offset_name
(
self
)
->
str
:
"""Spatial offset as string, East-West for x, North-South for y and Top-Bottom for z coordinate.
Example:
>>> f = fields("f: double[2D]")
>>> f[1, 1].offset_name # north-east
'NE'
"""
return
self
.
_offsetName
@
property
def
index
(
self
):
"""Value of index coordinates as tuple."""
return
self
.
_index
def
neighbor
(
self
,
coord_id
:
int
,
offset
:
Sequence
[
int
])
->
'Field.Access'
:
"""Returns a new Access with changed spatial coordinates.
Args:
coord_id: index of the coordinate to change (0 for x, 1 for y,...)
offset: incremental change of this coordinate
Example:
>>> f = fields('f: [2D]')
>>> f[0,0].neighbor(coord_id=1, offset=-1)
f_S
"""
offset_list
=
list
(
self
.
offsets
)
offset_list
[
coord_id
]
+=
offset
return
Field
.
Access
(
self
.
field
,
tuple
(
offset_list
),
self
.
index
)
def
get_shifted
(
self
,
*
shift
)
->
'Field.Access'
:
"""Returns a new Access with changed spatial coordinates
Example:
>>> f = fields("f: [2D]")
>>> f[0,0].get_shifted(1, 1)
f_NE
"""
return
Field
.
Access
(
self
.
field
,
tuple
(
a
+
b
for
a
,
b
in
zip
(
shift
,
self
.
offsets
)),
self
.
index
)
def
at_index
(
self
,
*
idx_tuple
):
def
at_index
(
self
,
*
idx_tuple
)
->
'Field.Access'
:
"""Returns new Access with changed index.
Example:
>>> f = fields("f(9): [2D]")
>>> f(0).at_index(8)
f_C^8
"""
return
Field
.
Access
(
self
.
field
,
self
.
offsets
,
idx_tuple
)
def
_hashable_content
(
self
):
...
...
@@ -729,12 +790,12 @@ def direction_string_to_offset(direction: str, dim: int = 3):
direction
=
direction
[
1
:]
return
offset
[:
dim
]
def
_parse_type_description
(
type_description
):
if
not
type_description
:
return
np
.
float64
,
None
elif
'['
in
type_description
:
assert
type_description
[
-
1
]
==
']'
splitted
=
type_description
[:
-
1
].
split
(
"["
,
)
type_part
,
size_part
=
type_description
[:
-
1
].
split
(
"["
,
)
if
not
type_part
:
type_part
=
"float64"
...
...
gpucuda/__init__.py
View file @
4651c191
from
pystencils.gpucuda.kernelcreation
import
create_cuda_kernel
,
created_indexed_cuda_kernel
from
pystencils.gpucuda.cudajit
import
make_python_function
__all__
=
[
'create_cuda_kernel'
,
'created_indexed_cuda_kernel'
,
'make_python_function'
]
from
.indexing
import
AbstractIndexing
,
BlockIndexing
,
LineIndexing
__all__
=
[
'create_cuda_kernel'
,
'created_indexed_cuda_kernel'
,
'make_python_function'
,
'AbstractIndexing'
,
'BlockIndexing'
,
'LineIndexing'
]
gpucuda/indexing.py
View file @
4651c191
import
abc
from
typing
import
Tuple
from
typing
import
Tuple
# noqa
import
sympy
as
sp
from
pystencils.astnodes
import
Conditional
,
Block
from
pystencils.slicing
import
normalize_slice
...
...
kernelcreation.py
View file @
4651c191
from
types
import
MappingProxyType
import
sympy
as
sp
from
.assignment
import
Assignment
from
.astnodes
import
LoopOverCoordinate
,
Conditional
,
Block
,
SympyAssignment
from
.simp.assignment_collection
import
AssignmentCollection
from
.gpucuda.indexing
import
indexing_creator_from_params
from
.transformations
import
remove_conditionals_in_staggered_kernel
from
pystencils
.assignment
import
Assignment
from
pystencils
.astnodes
import
LoopOverCoordinate
,
Conditional
,
Block
,
SympyAssignment
from
pystencils
.simp.assignment_collection
import
AssignmentCollection
from
pystencils
.gpucuda.indexing
import
indexing_creator_from_params
from
pystencils
.transformations
import
remove_conditionals_in_staggered_kernel
def
create_kernel
(
assignments
,
target
=
'cpu'
,
data_type
=
"double"
,
iteration_slice
=
None
,
ghost_layers
=
None
,
...
...
@@ -14,7 +14,7 @@ def create_kernel(assignments, target='cpu', data_type="double", iteration_slice
Creates abstract syntax tree (AST) of kernel, using a list of update equations.
Args:
assignments:
either be a plain list of equation
s or a
AssignmentCollection
object
assignments:
can be a single assignment, sequence of assignment
s or a
n `
AssignmentCollection
`
target: 'cpu', 'llvm' or 'gpu'
data_type: data type used for all untyped symbols (i.e. non-fields), can also be a dict from symbol name
to type
...
...
@@ -22,18 +22,34 @@ def create_kernel(assignments, target='cpu', data_type="double", iteration_slice
part of the field is iterated over
ghost_layers: if left to default, the number of necessary ghost layers is determined automatically
a single integer specifies the ghost layer count at all borders, can also be a sequence of
pairs [(x_lower_gl, x_upper_gl), .... ]
pairs
``
[(x_lower_gl, x_upper_gl), .... ]
``
cpu_openmp: True or number of threads for OpenMP parallelization, False for no OpenMP
cpu_vectorize_info: pair of instruction set name
(
'sse, 'avx'
,
'avx512'
) and data type ('float', 'double')
gpu_indexing: either 'block' or 'line' , or custom indexing class
(
see
gpucuda/i
ndexing
.py)
cpu_vectorize_info: pair of instruction set name
, i.e. one of
'sse, 'avx'
or
'avx512'
and data type 'float' or 'double'. For example ``('avx', 'double')``
gpu_indexing: either 'block' or 'line' , or custom indexing class
,
see
`pystencils.gpucuda.AbstractI
ndexing
`
gpu_indexing_params: dict with indexing parameters (constructor parameters of indexing class)
e.g. for 'block' one can specify {'block_size': (20, 20, 10) }
e.g. for 'block' one can specify
'
{'block_size': (20, 20, 10) }
'
Returns:
abstract syntax tree object, that can either be printed as source code with `show_code` or can be compiled with
through its `compile()` member
abstract syntax tree (AST) object, that can either be printed as source code with `show_code` or
can be compiled with through its `compile()` member
Example:
>>> import pystencils as ps
>>> import numpy as np
>>> s, d = ps.fields('s, d: [2D]')
>>> assignment = ps.Assignment(d[0,0], s[0, 1] + s[0, -1] + s[1, 0] + s[-1, 0])
>>> ast = ps.create_kernel(assignment, target='cpu', cpu_openmp=True)
>>> kernel = ast.compile()
>>> d_arr = np.zeros([5, 5])
>>> kernel(d=d_arr, s=np.ones([5, 5]))
>>> d_arr
array([[0., 0., 0., 0., 0.],
[0., 4., 4., 4., 0.],
[0., 4., 4., 4., 0.],
[0., 4., 4., 4., 0.],
[0., 0., 0., 0., 0.]])
"""
# ---- Normalizing parameters
...
...
@@ -89,9 +105,33 @@ def create_indexed_kernel(assignments, index_fields, target='cpu', data_type="do
index_fields: list of index fields, i.e. 1D fields with struct data type
coordinate_names: name of the coordinate fields in the struct data type
"""
if
isinstance
(
assignments
,
AssignmentCollection
):
Example:
>>> import pystencils as ps
>>> import numpy as np
>>>
>>> # Index field stores the indices of the cell to visit together with optional values
>>> index_arr_dtype = np.dtype([('x', np.int32), ('y', np.int32), ('val', np.double)])
>>> index_arr = np.array([(1, 1, 0.1), (2, 2, 0.2), (3, 3, 0.3)], dtype=index_arr_dtype)
>>> idx_field = ps.fields(idx=index_arr)
>>>
>>> # Additional values stored in index field can be accessed in the kernel as well
>>> s, d = ps.fields('s, d: [2D]')
>>> assignment = ps.Assignment(d[0,0], 2 * s[0, 1] + 2 * s[1, 0] + idx_field('val'))
>>> ast = create_indexed_kernel(assignment, [idx_field], coordinate_names=('x', 'y'))
>>> kernel = ast.compile()
>>> d_arr = np.zeros([5, 5])
>>> kernel(s=np.ones([5, 5]), d=d_arr, idx=index_arr)
>>> d_arr
array([[0. , 0. , 0. , 0. , 0. ],
[0. , 4.1, 0. , 0. , 0. ],
[0. , 0. , 4.2, 0. , 0. ],
[0. , 0. , 0. , 4.3, 0. ],
[0. , 0. , 0. , 0. , 0. ]])
"""
if
isinstance
(
assignments
,
Assignment
):
assignments
=
[
assignments
]
elif
isinstance
(
assignments
,
AssignmentCollection
):
assignments
=
assignments
.
all_assignments
if
target
==
'cpu'
:
from
pystencils.cpu
import
create_indexed_kernel
...
...
@@ -116,10 +156,12 @@ def create_indexed_kernel(assignments, index_fields, target='cpu', data_type="do
def
create_staggered_kernel
(
staggered_field
,
expressions
,
subexpressions
=
(),
target
=
'cpu'
,
**
kwargs
):
"""Kernel that updates a staggered field.
.. image:: /img/staggered_grid.svg
Args:
staggered_field: field that has one index coordinate and
where e.g. f[0,0](0) is interpreted as value at the left cell boundary, f[1,0](0) the right cell
boundary and f[0,0](1) the southern cell boundary etc.
where e.g.
``
f[0,0](0)
``
is interpreted as value at the left cell boundary,
``
f[1,0](0)
``
the right cell
boundary and
``
f[0,0](1)
``
the southern cell boundary etc.
expressions: sequence of expressions of length dim, defining how the east, southern, (bottom) cell boundary
should be update
subexpressions: optional sequence of Assignments, that define subexpressions used in the main expressions
...
...
plot2d.py
View file @
4651c191
...
...
@@ -253,6 +253,36 @@ def vector_field_magnitude_animation(run_function, plot_setup_function=lambda *_
return
animation
.
FuncAnimation
(
fig
,
update_figure
,
interval
=
interval
,
frames
=
frames
)
def
scalar_field_animation
(
run_function
,
plot_setup_function
=
lambda
*
_
:
None
,
rescale
=
True
,
plot_update_function
=
lambda
*
_
:
None
,
interval
=
30
,
frames
=
180
,
**
kwargs
):
import
matplotlib.animation
as
animation
fig
=
gcf
()
im
=
None
field
=
run_function
()
if
rescale
:
f_min
,
f_max
=
np
.
min
(
field
),
np
.
max
(
field
)
field
=
(
field
-
f_min
)
/
(
f_max
-
f_min
)
im
=
scalar_field
(
field
,
vmin
=
0.0
,
vmax
=
1.0
,
**
kwargs
)
else
:
im
=
scalar_field
(
field
,
**
kwargs
)
plot_setup_function
(
im
)
def
update_figure
(
*
_
):
f
=
run_function
()
if
rescale
:
f_min
,
f_max
=
np
.
min
(
f
),
np
.
max
(
f
)
f
=
(
f
-
f_min
)
/
(
f_max
-
f_min
)
if
hasattr
(
f
,
'mask'
):
f
=
np
.
ma
.
masked_array
(
f
,
mask
=
f
.
mask
[:,
:,
0
])
f
=
np
.
swapaxes
(
f
,
0
,
1
)
im
.
set_array
(
f
)
plot_update_function
(
im
)
return
im
,
return
animation
.
FuncAnimation
(
fig
,
update_figure
,
interval
=
interval
,
frames
=
frames
)
def
surface_plot_animation
(
run_function
,
frames
=
90
,
interval
=
30
):
from
mpl_toolkits.mplot3d
import
Axes3D
import
matplotlib.animation
as
animation
...
...
transformations.py
View file @
4651c191
...
...
@@ -687,6 +687,7 @@ class KernelConstraintsCheck:
def
process_expression
(
self
,
rhs
):
self
.
_update_accesses_rhs
(
rhs
)
if
isinstance
(
rhs
,
Field
.
Access
):
self
.
fields_read
.
add
(
rhs
.
field
)
return
rhs
elif
isinstance
(
rhs
,
TypedSymbol
):
return
rhs
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment