From 4651c191bc366300c025d35e98427f57019fa3e8 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Thu, 26 Apr 2018 19:17:42 +0200
Subject: [PATCH] Documentation overhaul Part2
---
__init__.py | 5 +-
cpu/cpujit.py | 6 ++
display_utils.py | 2 +-
field.py | 147 +++++++++++++++++++++++++++++++-------------
gpucuda/__init__.py | 5 +-
gpucuda/indexing.py | 2 +-
kernelcreation.py | 76 ++++++++++++++++++-----
plot2d.py | 30 +++++++++
transformations.py | 1 +
9 files changed, 207 insertions(+), 67 deletions(-)
diff --git a/__init__.py b/__init__.py
index 20db8d6b9..5dd9cd572 100644
--- a/__init__.py
+++ b/__init__.py
@@ -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']
-
diff --git a/cpu/cpujit.py b/cpu/cpujit.py
index b793302e2..ca915d38b 100644
--- a/cpu/cpujit.py
+++ b/cpu/cpujit.py
@@ -1,5 +1,11 @@
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``
diff --git a/display_utils.py b/display_utils.py
index 612d69aa3..2d0fa87f7 100644
--- a/display_utils.py
+++ b/display_utils.py
@@ -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.
"""
diff --git a/field.py b/field.py
index 87b4785a2..69513a3ab 100644
--- a/field.py
+++ b/field.py
@@ -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([20, 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[10, 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
+ Examples 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"
diff --git a/gpucuda/__init__.py b/gpucuda/__init__.py
index f360794ee..9b178c470 100644
--- a/gpucuda/__init__.py
+++ b/gpucuda/__init__.py
@@ -1,4 +1,5 @@
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']
diff --git a/gpucuda/indexing.py b/gpucuda/indexing.py
index 862143abb..d39a54688 100644
--- a/gpucuda/indexing.py
+++ b/gpucuda/indexing.py
@@ -1,5 +1,5 @@
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
diff --git a/kernelcreation.py b/kernelcreation.py
index 17d804b2d..87bcf5d49 100644
--- a/kernelcreation.py
+++ b/kernelcreation.py
@@ -1,10 +1,10 @@
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 equations or a AssignmentCollection object
+ assignments: can be a single assignment, sequence of assignments or an `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/indexing.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.AbstractIndexing`
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
diff --git a/plot2d.py b/plot2d.py
index 6157f9314..68cc8d729 100644
--- a/plot2d.py
+++ b/plot2d.py
@@ -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
diff --git a/transformations.py b/transformations.py
index 62b78258d..727aae027 100644
--- a/transformations.py
+++ b/transformations.py
@@ -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
--
GitLab