From 6654fa1727be20107e6702a5f915b8a22e78936a Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Tue, 13 Oct 2020 13:16:38 +0200
Subject: [PATCH] Extend testsuit

---
 .gitlab-ci.yml                                |   2 +
 pystencils/backends/cbackend.py               |   8 +-
 pystencils/backends/cuda_backend.py           |  13 +-
 pystencils/backends/json.py                   |  43 +---
 pystencils/backends/opencl_backend.py         |  13 +-
 pystencils/datahandling/pycuda.py             |  23 +-
 pystencils/datahandling/pyopencl.py           |  12 +-
 pystencils/fd/derivation.py                   |   3 +-
 pystencils/fd/finitedifferences.py            |  34 +--
 pystencils/fd/finitevolumes.py                |  53 ++---
 pystencils/gpucuda/cudajit.py                 |   5 +-
 pystencils/gpucuda/texture_utils.py           |  73 ------
 pystencils/jupyter.py                         |  50 +---
 pystencils/opencl/opencljit.py                |   5 +-
 pystencils_tests/test_datahandling.py         |  35 +++
 .../test_fd_derivation_via_rotation.ipynb     |   2 +-
 pystencils_tests/test_finite_differences.py   |  17 +-
 pystencils_tests/test_json_backend.py         |   8 +-
 .../test_jupyter_extensions.ipynb             | 224 ++++++++++++++++++
 .../{test_jacobi_llvm.py => test_llvm.py}     |  47 +++-
 .../test_simplification_strategy.py           |   5 +
 pystencils_tests/test_stencils.py             |   2 +-
 pytest.ini                                    |   4 +-
 23 files changed, 423 insertions(+), 258 deletions(-)
 create mode 100644 pystencils_tests/test_jupyter_extensions.ipynb
 rename pystencils_tests/{test_jacobi_llvm.py => test_llvm.py} (72%)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index f0481a590..39f6c8065 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -93,6 +93,7 @@ minimal-conda:
     - python setup.py quicktest
   tags:
     - docker
+    - cuda
 
 
 minimal-sympy-master:
@@ -107,6 +108,7 @@ minimal-sympy-master:
   allow_failure: true
   tags:
     - docker
+    - cuda
 
 
 pycodegen-integration:
diff --git a/pystencils/backends/cbackend.py b/pystencils/backends/cbackend.py
index d4f431ffb..0227b791f 100644
--- a/pystencils/backends/cbackend.py
+++ b/pystencils/backends/cbackend.py
@@ -45,9 +45,11 @@ def generate_c(ast_node: Node,
     create_kernel functions.
 
     Args:
-        ast_node:
-        signature_only:
-        dialect: 'c' or 'cuda'
+        ast_node: ast representation of kernel
+        signature_only: generate signature without function body
+        dialect: 'c', 'cuda' or opencl
+        custom_backend: use own custom printer for code generation
+        with_globals: enable usage of global variables
     Returns:
         C-like code for the ast node and its descendants
     """
diff --git a/pystencils/backends/cuda_backend.py b/pystencils/backends/cuda_backend.py
index 2d7dc579e..b699c9ce4 100644
--- a/pystencils/backends/cuda_backend.py
+++ b/pystencils/backends/cuda_backend.py
@@ -10,17 +10,20 @@ with open(join(dirname(__file__), 'cuda_known_functions.txt')) as f:
     CUDA_KNOWN_FUNCTIONS = {l.strip(): l.strip() for l in lines if l}
 
 
-def generate_cuda(astnode: Node, signature_only: bool = False) -> str:
+def generate_cuda(ast_node: Node, signature_only: bool = False, custom_backend=None, with_globals=True) -> str:
     """Prints an abstract syntax tree node as CUDA code.
 
     Args:
-        astnode: KernelFunction node to generate code for
-        signature_only: if True only the signature is printed
+        ast_node: ast representation of kernel
+        signature_only: generate signature without function body
+        custom_backend: use own custom printer for code generation
+        with_globals: enable usage of global variables
 
     Returns:
-        C-like code for the ast node and its descendants
+        CUDA code for the ast node and its descendants
     """
-    return generate_c(astnode, signature_only, dialect='cuda')
+    return generate_c(ast_node, signature_only, dialect='cuda',
+                      custom_backend=custom_backend, with_globals=with_globals)
 
 
 class CudaBackend(CBackend):
diff --git a/pystencils/backends/json.py b/pystencils/backends/json.py
index a0fcc9aa5..954965c6e 100644
--- a/pystencils/backends/json.py
+++ b/pystencils/backends/json.py
@@ -12,22 +12,10 @@ import json
 from pystencils.astnodes import NodeOrExpr
 from pystencils.backends.cbackend import CustomSympyPrinter, generate_c
 
-try:
-    import toml
-except Exception:
-    class toml:
-        def dumps(self, *args):
-            raise ImportError('toml not installed')
-
-        def dump(self, *args):
-            raise ImportError('toml not installed')
-
 try:
     import yaml
-except Exception:
-    class yaml:
-        def dumps(self, *args):
-            raise ImportError('pyyaml not installed')
+except ImportError:
+    raise ImportError('yaml not installed')
 
 
 def expr_to_dict(expr_or_node: NodeOrExpr, with_c_code=True, full_class_names=False):
@@ -67,8 +55,8 @@ def print_json(expr_or_node: NodeOrExpr):
     Returns:
         str: JSON representation of AST
     """
-    dict = expr_to_dict(expr_or_node)
-    return json.dumps(dict, indent=4)
+    expr_or_node_dict = expr_to_dict(expr_or_node)
+    return json.dumps(expr_or_node_dict, indent=4)
 
 
 def write_json(filename: str, expr_or_node: NodeOrExpr):
@@ -78,28 +66,17 @@ def write_json(filename: str, expr_or_node: NodeOrExpr):
         filename (str): Path for the file to write
         expr_or_node (NodeOrExpr): a SymPy expression or a :class:`pystencils.astnodes.Node`
     """
-    dict = expr_to_dict(expr_or_node)
-    with open(filename, 'w') as f:
-        json.dump(dict, f, indent=4)
-
-
-def print_toml(expr_or_node):
-    dict = expr_to_dict(expr_or_node, full_class_names=False)
-    return toml.dumps(dict)
-
-
-def write_toml(filename, expr_or_node):
-    dict = expr_to_dict(expr_or_node)
+    expr_or_node_dict = expr_to_dict(expr_or_node)
     with open(filename, 'w') as f:
-        toml.dump(dict, f)
+        json.dump(expr_or_node_dict, f, indent=4)
 
 
 def print_yaml(expr_or_node):
-    dict = expr_to_dict(expr_or_node, full_class_names=False)
-    return yaml.dump(dict)
+    expr_or_node_dict = expr_to_dict(expr_or_node, full_class_names=False)
+    return yaml.dump(expr_or_node_dict)
 
 
 def write_yaml(filename, expr_or_node):
-    dict = expr_to_dict(expr_or_node)
+    expr_or_node_dict = expr_to_dict(expr_or_node)
     with open(filename, 'w') as f:
-        yaml.dump(dict, f)
+        yaml.dump(expr_or_node_dict, f)
diff --git a/pystencils/backends/opencl_backend.py b/pystencils/backends/opencl_backend.py
index 1813d63ef..9fe7b9788 100644
--- a/pystencils/backends/opencl_backend.py
+++ b/pystencils/backends/opencl_backend.py
@@ -11,17 +11,20 @@ with open(join(dirname(__file__), 'opencl1.1_known_functions.txt')) as f:
     OPENCL_KNOWN_FUNCTIONS = {l.strip(): l.strip() for l in lines if l}
 
 
-def generate_opencl(astnode: Node, signature_only: bool = False) -> str:
+def generate_opencl(ast_node: Node, signature_only: bool = False, custom_backend=None, with_globals=True) -> str:
     """Prints an abstract syntax tree node (made for target 'gpu') as OpenCL code.
 
     Args:
-        astnode: KernelFunction node to generate code for
-        signature_only: if True only the signature is printed
+        ast_node: ast representation of kernel
+        signature_only: generate signature without function body
+        custom_backend: use own custom printer for code generation
+        with_globals: enable usage of global variables
 
     Returns:
-        C-like code for the ast node and its descendants
+        OpenCL code for the ast node and its descendants
     """
-    return generate_c(astnode, signature_only, dialect='opencl')
+    return generate_c(ast_node, signature_only, dialect='opencl',
+                      custom_backend=custom_backend, with_globals=with_globals)
 
 
 class OpenClBackend(CudaBackend):
diff --git a/pystencils/datahandling/pycuda.py b/pystencils/datahandling/pycuda.py
index 14e255c85..1c65a1e9b 100644
--- a/pystencils/datahandling/pycuda.py
+++ b/pystencils/datahandling/pycuda.py
@@ -13,26 +13,31 @@ class PyCudaArrayHandler:
         import pycuda.autoinit  # NOQA
 
     def zeros(self, shape, dtype=np.float64, order='C'):
-        return gpuarray.zeros(shape, dtype, order)
+        cpu_array = np.zeros(shape=shape, dtype=dtype, order=order)
+        return self.to_gpu(cpu_array)
 
-    def ones(self, shape, dtype, order='C'):
-        return gpuarray.ones(shape, dtype, order)
+    def ones(self, shape, dtype=np.float64, order='C'):
+        cpu_array = np.ones(shape=shape, dtype=dtype, order=order)
+        return self.to_gpu(cpu_array)
 
     def empty(self, shape, dtype=np.float64, layout=None):
         if layout:
-            cpu_array = pystencils.field.create_numpy_array_with_layout(shape, dtype, layout)
+            cpu_array = pystencils.field.create_numpy_array_with_layout(shape=shape, dtype=dtype, layout=layout)
             return self.to_gpu(cpu_array)
         else:
             return gpuarray.empty(shape, dtype)
 
-    def to_gpu(self, array):
+    @staticmethod
+    def to_gpu(array):
         return gpuarray.to_gpu(array)
 
-    def upload(self, gpuarray, numpy_array):
-        gpuarray.set(numpy_array)
+    @staticmethod
+    def upload(array, numpy_array):
+        array.set(numpy_array)
 
-    def download(self, gpuarray, numpy_array):
-        gpuarray.get(numpy_array)
+    @staticmethod
+    def download(array, numpy_array):
+        array.get(numpy_array)
 
     def randn(self, shape, dtype=np.float64):
         cpu_array = np.random.randn(*shape).astype(dtype)
diff --git a/pystencils/datahandling/pyopencl.py b/pystencils/datahandling/pyopencl.py
index b4e53150f..2466c80cd 100644
--- a/pystencils/datahandling/pyopencl.py
+++ b/pystencils/datahandling/pyopencl.py
@@ -19,15 +19,17 @@ class PyOpenClArrayHandler:
         self.queue = queue
 
     def zeros(self, shape, dtype=np.float64, order='C'):
-        return gpuarray.zeros(shape, dtype, order)
+        cpu_array = np.zeros(shape=shape, dtype=dtype, order=order)
+        return self.to_gpu(cpu_array)
 
-    def ones(self, shape, dtype, order='C'):
-        return gpuarray.ones(self.queue, shape, dtype, order)
+    def ones(self, shape, dtype=np.float64, order='C'):
+        cpu_array = np.ones(shape=shape, dtype=dtype, order=order)
+        return self.to_gpu(cpu_array)
 
     def empty(self, shape, dtype=np.float64, layout=None):
         if layout:
-            cpu_array = pystencils.field.create_numpy_array_with_layout(shape, dtype, layout)
-            return self.from_numpy(cpu_array)
+            cpu_array = pystencils.field.create_numpy_array_with_layout(shape=shape, dtype=dtype, layout=layout)
+            return self.to_gpu(cpu_array)
         else:
             return gpuarray.empty(self.queue, shape, dtype)
 
diff --git a/pystencils/fd/derivation.py b/pystencils/fd/derivation.py
index 23579f488..868a6872e 100644
--- a/pystencils/fd/derivation.py
+++ b/pystencils/fd/derivation.py
@@ -126,7 +126,6 @@ class FiniteDifferenceStencilDerivation:
 
     def isotropy_equations(self, order):
         def cycle_int_sequence(sequence, modulus):
-            import numpy as np
             result = []
             arr = np.array(sequence, dtype=int)
             while True:
@@ -200,7 +199,7 @@ class FiniteDifferenceStencilDerivation:
                via rotation and apply them to a field."""
             dim = len(self.stencil[0])
             assert (dim == 2 or dim == 3), "This function is only for 2D or 3D stencils available"
-            rotated_weights = np.rot90(np.array(self), 1, axes)
+            rotated_weights = np.rot90(np.array(self.__array__()), 1, axes)
 
             result = []
             max_offset = max(max(abs(e) for e in direction) for direction in self.stencil)
diff --git a/pystencils/fd/finitedifferences.py b/pystencils/fd/finitedifferences.py
index 1a119f71b..20d071793 100644
--- a/pystencils/fd/finitedifferences.py
+++ b/pystencils/fd/finitedifferences.py
@@ -23,9 +23,9 @@ def diffusion(scalar, diffusion_coeff, idx=None):
         >>> f = Field.create_generic('f', spatial_dimensions=2)
         >>> d = sp.Symbol("d")
         >>> dx = sp.Symbol("dx")
-        >>> diffusion_term = diffusion(scalar=f, diffusion_coeff=sp.Symbol("d"))
+        >>> diffusion_term = diffusion(scalar=f, diffusion_coeff=d)
         >>> discretization = Discretization2ndOrder()
-        >>> expected_output = ((f[-1, 0] + f[0, -1] - 4*f[0, 0] + f[0, 1] + f[1, 0]) * d) / dx**2
+        >>> expected_output = ((f[-1, 0] + f[0, -1] - 4 * f[0, 0] + f[0, 1] + f[1, 0]) * d) / dx**2
         >>> sp.simplify(discretization(diffusion_term) - expected_output)
         0
     """
@@ -79,13 +79,6 @@ class Discretization2ndOrder:
         self.dt = dt
         self.spatial_stencil = discretization_stencil_func
 
-    @staticmethod
-    def _diff_order(e):
-        if not isinstance(e, Diff):
-            return 0
-        else:
-            return 1 + Discretization2ndOrder._diff_order(e.args[0])
-
     def _discretize_diffusion(self, e):
         result = 0
         for c in range(e.dim):
@@ -121,29 +114,6 @@ class Discretization2ndOrder:
             new_args = [self._discretize_spatial(a) for a in e.args]
             return e.func(*new_args) if new_args else e
 
-    def _discretize_diff(self, e):
-        order = self._diff_order(e)
-        if order == 1:
-            fa = e.args[0]
-            index = e.target
-            return (fa.neighbor(index, 1) - fa.neighbor(index, -1)) / (2 * self.dx)
-        elif order == 2:
-            indices = sorted([e.target, e.args[0].target])
-            fa = e.args[0].args[0]
-            if indices[0] == indices[1] and all(i >= 0 for i in indices):
-                result = (-2 * fa + fa.neighbor(indices[0], -1) + fa.neighbor(indices[0], +1))
-            elif indices[0] == indices[1]:
-                result = 0
-                for d in range(fa.field.spatial_dimensions):
-                    result += (-2 * fa + fa.neighbor(d, -1) + fa.neighbor(d, +1))
-            else:
-                assert all(i >= 0 for i in indices)
-                offsets = [(1, 1), [-1, 1], [1, -1], [-1, -1]]
-                result = sum(o1 * o2 * fa.neighbor(indices[0], o1).neighbor(indices[1], o2) for o1, o2 in offsets) / 4
-            return result / (self.dx ** 2)
-        else:
-            raise NotImplementedError("Term contains derivatives of order > 2")
-
     def __call__(self, expr):
         if isinstance(expr, list):
             return [self(e) for e in expr]
diff --git a/pystencils/fd/finitevolumes.py b/pystencils/fd/finitevolumes.py
index d2ddc3c79..d8f46f394 100644
--- a/pystencils/fd/finitevolumes.py
+++ b/pystencils/fd/finitevolumes.py
@@ -7,6 +7,25 @@ from collections import defaultdict
 from collections.abc import Iterable
 
 
+def get_access_and_direction(term):
+    direction1 = term.args[1]
+    if isinstance(term.args[0], ps.Field.Access):  # first derivative
+        access = term.args[0]
+        direction = (direction1,)
+    elif isinstance(term.args[0], ps.fd.Diff):  # nested derivative
+        if isinstance(term.args[0].args[0], ps.fd.Diff):  # third or higher derivative
+            raise ValueError("can only handle first and second derivatives")
+        elif not isinstance(term.args[0].args[0], ps.Field.Access):
+            raise ValueError("can only handle derivatives of field accesses")
+
+        access, direction2 = term.args[0].args[:2]
+        direction = (direction1, direction2)
+    else:
+        raise NotImplementedError(f"can only deal with derivatives of field accesses, "
+                                  f"but not {type(term.args[0])}; expansion of derivatives probably failed")
+    return access, direction
+
+
 class FVM1stOrder:
     """Finite-volume discretization
 
@@ -48,22 +67,7 @@ class FVM1stOrder:
                 avg = (term.get_shifted(*neighbor) + term) * sp.Rational(1, 2)
                 return avg
             elif isinstance(term, ps.fd.Diff):
-                direction1 = term.args[1]
-                if isinstance(term.args[0], ps.Field.Access):  # first derivative
-                    access = term.args[0]
-                    direction = (direction1,)
-                elif isinstance(term.args[0], ps.fd.Diff):  # nested derivative
-                    if isinstance(term.args[0].args[0], ps.fd.Diff):  # third or higher derivative
-                        raise ValueError("can only handle first and second derivatives")
-                    elif not isinstance(term.args[0].args[0], ps.Field.Access):
-                        raise ValueError("can only handle derivatives of field accesses")
-
-                    access, direction2 = term.args[0].args[:2]
-                    direction = (direction1, direction2)
-                else:
-                    raise NotImplementedError("can only deal with derivatives of field accesses, "
-                                              "but not {}; expansion of derivatives probably failed"
-                                              .format(type(term.args[0])))
+                access, direction = get_access_and_direction(term)
 
                 fds = FDS(neighbor, access.field.spatial_dimensions, direction)
                 return fds.apply(access)
@@ -103,22 +107,7 @@ class FVM1stOrder:
 
         def discretize(term):
             if isinstance(term, ps.fd.Diff):
-                direction1 = term.args[1]
-                if isinstance(term.args[0], ps.Field.Access):  # first derivative
-                    access = term.args[0]
-                    direction = (direction1,)
-                elif isinstance(term.args[0], ps.fd.Diff):  # nested derivative
-                    if isinstance(term.args[0].args[0], ps.fd.Diff):  # third or higher derivative
-                        raise ValueError("can only handle first and second derivatives")
-                    elif not isinstance(term.args[0].args[0], ps.Field.Access):
-                        raise ValueError("can only handle derivatives of field accesses")
-
-                    access, direction2 = term.args[0].args[:2]
-                    direction = (direction1, direction2)
-                else:
-                    raise NotImplementedError("can only deal with derivatives of field accesses, "
-                                              "but not {}; expansion of derivatives probably failed"
-                                              .format(type(term.args[0])))
+                access, direction = get_access_and_direction(term)
 
                 if self.dim == 2:
                     stencil = ["".join(a).replace(" ", "") for a in itertools.product("NS ", "EW ")
diff --git a/pystencils/gpucuda/cudajit.py b/pystencils/gpucuda/cudajit.py
index 86f251c97..03dfc8245 100644
--- a/pystencils/gpucuda/cudajit.py
+++ b/pystencils/gpucuda/cudajit.py
@@ -1,6 +1,7 @@
 import numpy as np
 
-from pystencils.backends.cbackend import generate_c, get_headers
+from pystencils.backends.cbackend import get_headers
+from pystencils.backends.cuda_backend import generate_cuda
 from pystencils.data_types import StructType
 from pystencils.field import FieldType
 from pystencils.gpucuda.texture_utils import ndarray_to_tex
@@ -45,7 +46,7 @@ def make_python_function(kernel_function_node, argument_dict=None, custom_backen
     code = includes + "\n"
     code += "#define FUNC_PREFIX __global__\n"
     code += "#define RESTRICT __restrict__\n\n"
-    code += str(generate_c(kernel_function_node, dialect='cuda', custom_backend=custom_backend))
+    code += str(generate_cuda(kernel_function_node, custom_backend=custom_backend))
     textures = set(d.interpolator for d in kernel_function_node.atoms(
         InterpolatorAccess) if isinstance(d.interpolator, TextureCachedField))
 
diff --git a/pystencils/gpucuda/texture_utils.py b/pystencils/gpucuda/texture_utils.py
index 8d5ba56c8..0b383507f 100644
--- a/pystencils/gpucuda/texture_utils.py
+++ b/pystencils/gpucuda/texture_utils.py
@@ -7,8 +7,6 @@
 """
 
 """
-
-from os.path import dirname, isdir, join
 from typing import Union
 
 import numpy as np
@@ -21,15 +19,6 @@ except Exception:
     pass
 
 
-def pow_two_divider(n):
-    if n == 0:
-        return 0
-    divider = 1
-    while (n & divider) == 0:
-        divider <<= 1
-    return divider
-
-
 def ndarray_to_tex(tex_ref,  # type: Union[cuda.TextureReference, cuda.SurfaceReference]
                    ndarray,
                    address_mode=None,
@@ -66,65 +55,3 @@ def ndarray_to_tex(tex_ref,  # type: Union[cuda.TextureReference, cuda.SurfaceRe
 
     if not read_as_integer:
         tex_ref.set_flags(tex_ref.get_flags() & ~cuda.TRSF_READ_AS_INTEGER)
-
-
-def prefilter_for_cubic_bspline(gpuarray):
-    import pycuda.autoinit  # NOQA
-    from pycuda.compiler import SourceModule
-
-    ndim = gpuarray.ndim
-    assert ndim == 2 or ndim == 3, "Only 2d or 3d supported"
-    assert isdir(join(dirname(__file__), "CubicInterpolationCUDA", "code")), \
-        "Submodule CubicInterpolationCUDA does not exist"
-    nvcc_options = ["-w", "-std=c++11", "-Wno-deprecated-gpu-targets"]
-    nvcc_options += ["-I" + join(dirname(__file__), "CubicInterpolationCUDA", "code")]
-    nvcc_options += ["-I" + join(dirname(__file__), "CubicInterpolationCUDA", "code", "internal")]
-
-    file_name = join(dirname(__file__), "CubicInterpolationCUDA", "code", "cubicPrefilter%iD.cu" % ndim)
-    with open(file_name) as file:
-        code = file.read()
-
-    mod = SourceModule(code, options=nvcc_options)
-
-    if ndim == 2:
-        height, width = gpuarray.shape
-        block = min(pow_two_divider(height), 64)
-        grid = height // block
-        func = mod.get_function('SamplesToCoefficients2DXf')
-        func(gpuarray, np.uint32(gpuarray.strides[-2]), *(np.uint32(r)
-                                                          for r in reversed(gpuarray.shape)),
-             block=(block, 1, 1),
-             grid=(grid, 1, 1))
-
-        block = min(pow_two_divider(width), 64)
-        grid = width // block
-        func = mod.get_function('SamplesToCoefficients2DYf')
-        func(gpuarray, np.uint32(gpuarray.strides[-2]), *(np.uint32(r)
-                                                          for r in reversed(gpuarray.shape)),
-             block=(block, 1, 1),
-             grid=(grid, 1, 1))
-    elif ndim == 3:
-        depth, height, width = gpuarray.shape
-        dimX = min(min(pow_two_divider(width), pow_two_divider(height)), 64)
-        dimY = min(min(pow_two_divider(depth), pow_two_divider(height)), 512 / dimX)
-        block = (dimX, dimY, 1)
-
-        dimGridX = (height // block[0], depth // block[1], 1)
-        dimGridY = (width // block[0], depth // block[1], 1)
-        dimGridZ = (width // block[0], height // block[1], 1)
-
-        func = mod.get_function("SamplesToCoefficients3DXf")
-        func(gpuarray, np.uint32(gpuarray.strides[-2]), *(np.uint32(r)
-                                                          for r in reversed(gpuarray.shape)),
-             block=block,
-             grid=dimGridX)
-        func = mod.get_function("SamplesToCoefficients3DYf")
-        func(gpuarray, np.uint32(gpuarray.strides[-2]), *(np.uint32(r)
-                                                          for r in reversed(gpuarray.shape)),
-             block=block,
-             grid=dimGridY)
-        func = mod.get_function("SamplesToCoefficients3DZf")
-        func(gpuarray, np.uint32(gpuarray.strides[-2]), *(np.uint32(r)
-                                                          for r in reversed(gpuarray.shape)),
-             block=block,
-             grid=dimGridZ)
diff --git a/pystencils/jupyter.py b/pystencils/jupyter.py
index 5e86e0c2c..50bbb3a34 100644
--- a/pystencils/jupyter.py
+++ b/pystencils/jupyter.py
@@ -7,55 +7,7 @@ from IPython.display import HTML
 
 import pystencils.plot as plt
 
-__all__ = ['log_progress', 'make_imshow_animation', 'display_animation', 'set_display_mode']
-
-
-def log_progress(sequence, every=None, size=None, name='Items'):
-    """Copied from https://github.com/alexanderkuk/log-progress"""
-    from ipywidgets import IntProgress, HTML, VBox
-    from IPython.display import display
-
-    is_iterator = False
-    if size is None:
-        try:
-            size = len(sequence)
-        except TypeError:
-            is_iterator = True
-    if size is not None:
-        if every is None:
-            if size <= 200:
-                every = 1
-            else:
-                every = int(size / 200)     # every 0.5%
-    else:
-        assert every is not None, 'sequence is iterator, set every'
-
-    if is_iterator:
-        progress = IntProgress(min=0, max=1, value=1)
-        progress.bar_style = 'info'
-    else:
-        progress = IntProgress(min=0, max=size, value=0)
-    label = HTML()
-    box = VBox(children=[label, progress])
-    display(box)
-
-    index = 0
-    try:
-        for index, record in enumerate(sequence, 1):
-            if index == 1 or index % every == 0:
-                if is_iterator:
-                    label.value = f'{name}: {index} / ?'
-                else:
-                    progress.value = index
-                    label.value = f'{name}: {index} / {size}'
-            yield record
-    except:
-        progress.bar_style = 'danger'
-        raise
-    else:
-        progress.bar_style = 'success'
-        progress.value = index
-        label.value = f"{name}: {str(index or '?')}"
+__all__ = ['make_imshow_animation', 'display_animation', 'set_display_mode']
 
 
 VIDEO_TAG = """<video controls width="80%">
diff --git a/pystencils/opencl/opencljit.py b/pystencils/opencl/opencljit.py
index 6b5893865..a33f51c07 100644
--- a/pystencils/opencl/opencljit.py
+++ b/pystencils/opencl/opencljit.py
@@ -1,6 +1,7 @@
 import numpy as np
 
-from pystencils.backends.cbackend import generate_c, get_headers
+from pystencils.backends.cbackend import get_headers
+from pystencils.backends.opencl_backend import generate_opencl
 from pystencils.gpucuda.cudajit import _build_numpy_argument_list, _check_arguments
 from pystencils.include import get_pystencils_include_path
 from pystencils.kernel_wrapper import KernelWrapper
@@ -91,7 +92,7 @@ def make_python_function(kernel_function_node, opencl_queue, opencl_ctx, argumen
     code = includes + "\n"
     code += "#define FUNC_PREFIX __kernel\n"
     code += "#define RESTRICT restrict\n\n"
-    code += str(generate_c(kernel_function_node, dialect='opencl', custom_backend=custom_backend))
+    code += str(generate_opencl(kernel_function_node, custom_backend=custom_backend))
     options = []
     if USE_FAST_MATH:
         options.append("-cl-unsafe-math-optimizations")
diff --git a/pystencils_tests/test_datahandling.py b/pystencils_tests/test_datahandling.py
index 949f06789..a9b6878c8 100644
--- a/pystencils_tests/test_datahandling.py
+++ b/pystencils_tests/test_datahandling.py
@@ -6,6 +6,8 @@ import numpy as np
 
 import pystencils as ps
 from pystencils import create_data_handling, create_kernel
+from pystencils.datahandling.pycuda import PyCudaArrayHandler
+from pystencils.datahandling.pyopencl import PyOpenClArrayHandler
 
 try:
     import pytest
@@ -355,3 +357,36 @@ def test_load_data():
     assert np.all(dh.cpu_arrays['src']) == 0
     assert np.all(dh.cpu_arrays['dst']) == 0
     assert np.all(dh.cpu_arrays['dst2']) == 0
+
+
+@pytest.mark.parametrize('target', ('gpu', 'opencl'))
+def test_array_handler(target):
+    size = (2, 2)
+    if target == 'gpu':
+        array_handler = PyCudaArrayHandler()
+    if target == 'opencl':
+        pytest.importorskip('pyopencl')
+        import pyopencl as cl
+        from pystencils.opencl.opencljit import init_globally
+        init_globally()
+        ctx = cl.create_some_context(0)
+        queue = cl.CommandQueue(ctx)
+        array_handler = PyOpenClArrayHandler(queue)
+
+    zero_array = array_handler.zeros(size)
+    cpu_array = np.empty(size)
+    array_handler.download(zero_array, cpu_array)
+    assert np.all(cpu_array) == 0
+
+    ones_array = array_handler.ones(size)
+    cpu_array = np.empty(size)
+    array_handler.download(ones_array, cpu_array)
+    assert np.all(cpu_array) == 1
+
+    empty = array_handler.empty(size)
+    assert empty.strides == (16, 8)
+    empty = array_handler.empty(shape=size, layout=(1, 0))
+    assert empty.strides == (8, 16)
+
+    random_array = array_handler.randn(size)
+
diff --git a/pystencils_tests/test_fd_derivation_via_rotation.ipynb b/pystencils_tests/test_fd_derivation_via_rotation.ipynb
index 485a5452e..bebab9462 100644
--- a/pystencils_tests/test_fd_derivation_via_rotation.ipynb
+++ b/pystencils_tests/test_fd_derivation_via_rotation.ipynb
@@ -187,7 +187,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.8.2"
   }
  },
  "nbformat": 4,
diff --git a/pystencils_tests/test_finite_differences.py b/pystencils_tests/test_finite_differences.py
index 2491247c2..b6d528d17 100644
--- a/pystencils_tests/test_finite_differences.py
+++ b/pystencils_tests/test_finite_differences.py
@@ -1,8 +1,9 @@
 import sympy as sp
+import pytest
 
 import pystencils as ps
 from pystencils.astnodes import LoopOverCoordinate
-from pystencils.fd import diff
+from pystencils.fd import diff, diffusion, Discretization2ndOrder
 from pystencils.fd.spatial import discretize_spatial, fd_stencils_isotropic, fd_stencils_standard
 
 
@@ -69,3 +70,17 @@ def test_staggered_combined():
 
     to_test = ps.fd.discretize_spatial_staggered(expr, dx)
     assert sp.expand(reference - to_test) == 0
+
+
+def test_diffusion():
+    f = ps.fields("f(3): [2D]")
+    d = sp.Symbol("d")
+    dx = sp.Symbol("dx")
+    idx = 2
+    diffusion_term = diffusion(scalar=f, diffusion_coeff=sp.Symbol("d"), idx=idx)
+    discretization = Discretization2ndOrder()
+    expected_output = ((f[-1, 0](idx) + f[0, -1](idx) - 4 * f[0, 0](idx) + f[0, 1](idx) + f[1, 0](idx)) * d) / dx ** 2
+    assert sp.simplify(discretization(diffusion_term) - expected_output) == 0
+
+    with pytest.raises(ValueError):
+        diffusion(scalar=d, diffusion_coeff=sp.Symbol("d"), idx=idx)
diff --git a/pystencils_tests/test_json_backend.py b/pystencils_tests/test_json_backend.py
index bc9765c47..a3fb2420c 100644
--- a/pystencils_tests/test_json_backend.py
+++ b/pystencils_tests/test_json_backend.py
@@ -11,7 +11,8 @@
 import sympy
 
 import pystencils
-from pystencils.backends.json import print_json
+from pystencils.backends.json import print_json, print_yaml, write_json, write_yaml
+import tempfile
 
 
 def test_json_backend():
@@ -26,3 +27,8 @@ def test_json_backend():
     ast = pystencils.create_kernel(assignments)
 
     print(print_json(ast))
+    print(print_yaml(ast))
+
+    temp_dir = tempfile.TemporaryDirectory()
+    write_json(temp_dir.name + '/test.json', ast)
+    write_yaml(temp_dir.name + '/test.yaml', ast)
diff --git a/pystencils_tests/test_jupyter_extensions.ipynb b/pystencils_tests/test_jupyter_extensions.ipynb
new file mode 100644
index 000000000..da8d17890
--- /dev/null
+++ b/pystencils_tests/test_jupyter_extensions.ipynb
@@ -0,0 +1,224 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pystencils.session import *"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dh = ps.create_data_handling(domain_size=(256, 256), periodicity=True)\n",
+    "c_field = dh.add_array('c')\n",
+    "dh.fill(\"c\", 0.0, ghost_layers=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for x in range(129):\n",
+    "    for y in range(258):\n",
+    "        dh.cpu_arrays['c'][x, y] = 1.0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.image.AxesImage at 0x7ff5b00707c0>"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAFlCAYAAAATVk7bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAUcklEQVR4nO3db4xl9X3f8c834JCoTlQsFkoX1KXuRgpEDVZWyBJPaB0V4kbBrkQFUl2kWiKpsGRLrlpwHjh9gGSpjVO1jS2RGpmqrulWtmWUOk0IdWVFio0Xl5p/pl4Zx2ygsKlbmbYSKfDtgzmUKcwyy86Mvzv3vl7S6N77O+fc+0O/ucu+dc89W90dAAAAmPAj0xMAAABgfYlSAAAAxohSAAAAxohSAAAAxohSAAAAxohSAAAAxpw7PYEkueCCC/rQoUPT0wDWwf95ZHoGkCR5+PsHct6J/zU9DUiS/NTP/cXpKQAr7sEHH/yT7j6w1bazIkoPHTqUY8eOTU8DWAMv/9fD01OAJMlf+uyv5O0f/ur0NCBJct+xfzs9BWDFVdUfnWqb03cBAAAYI0oBAAAYI0oBAAAYI0oBAAAYs22UVtWlVfXlqnq8qh6tqg8u479WVX9cVQ8tP+/edMztVXW8qp6oqmv38j8AAACA/et0rr77YpIPd/c3quonkjxYVfct236ju//x5p2r6vIkNya5IsmfT/L7VfVT3f3Sbk4cAACA/W/bT0q7+5nu/sZy//kkjyc5+AaHXJ/knu5+obufTHI8yVW7MVkAAABWy5v6TmlVHUryjiRfW4Y+UFXfrKq7qur8Zexgkqc2HXYiW0RsVd1SVceq6tjJkyff9MQBAADY/047SqvqrUk+l+RD3f2DJJ9M8vYkVyZ5Jsmvv7LrFof36wa67+zuI9195MCBA2964gAAAOx/pxWlVfWWbATpZ7r780nS3c9290vd/XKS38qrp+ieSHLppsMvSfL07k0ZAACAVXE6V9+tJJ9K8nh3f3zT+MWbdntvkkeW+/cmubGqzquqy5IcTvLA7k0ZAACAVXE6V9+9Osn7kjxcVQ8tYx9JclNVXZmNU3O/m+SXk6S7H62qo0key8aVe2915V0AAAC2sm2UdvcfZOvviX7pDY65I8kdO5gXAAAAa+BNXX0XAAAAdpMoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYIwoBQAAYMy2UVpVl1bVl6vq8ap6tKo+uIy/raruq6pvL7fnbzrm9qo6XlVPVNW1e/kfAAAAwP51Op+Uvpjkw93900nemeTWqro8yW1J7u/uw0nuXx5n2XZjkiuSXJfkE1V1zl5MHgAAgP1t2yjt7me6+xvL/eeTPJ7kYJLrk9y97HZ3kvcs969Pck93v9DdTyY5nuSq3Z44AAAA+9+b+k5pVR1K8o4kX0tyUXc/k2yEa5ILl90OJnlq02EnlrHXPtctVXWsqo6dPHnyzc8cAACAfe+0o7Sq3prkc0k+1N0/eKNdtxjr1w1039ndR7r7yIEDB053GgAAAKyQ04rSqnpLNoL0M939+WX42aq6eNl+cZLnlvETSS7ddPglSZ7enekCAACwSk7n6ruV5FNJHu/uj2/adG+Sm5f7Nyf54qbxG6vqvKq6LMnhJA/s3pQBAABYFeeexj5XJ3lfkoer6qFl7CNJPpbkaFW9P8n3ktyQJN39aFUdTfJYNq7ce2t3v7TrMwcAAGDf2zZKu/sPsvX3RJPkXac45o4kd+xgXgAAAKyBN3X1XQAAANhNohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAxohQAAIAx20ZpVd1VVc9V1SObxn6tqv64qh5aft69advtVXW8qp6oqmv3auIAAADsf6fzSemnk1y3xfhvdPeVy8+XkqSqLk9yY5IrlmM+UVXn7NZkAQAAWC3bRml3fyXJ90/z+a5Pck93v9DdTyY5nuSqHcwPAACAFbaT75R+oKq+uZzee/4ydjDJU5v2ObGMAQAAwOucaZR+Msnbk1yZ5Jkkv76M1xb79lZPUFW3VNWxqjp28uTJM5wGAAAA+9kZRWl3P9vdL3X3y0l+K6+eonsiyaWbdr0kydOneI47u/tIdx85cODAmUwDAACAfe6MorSqLt708L1JXrky771Jbqyq86rqsiSHkzywsykCAACwqs7dboeq+mySa5JcUFUnknw0yTVVdWU2Ts39bpJfTpLufrSqjiZ5LMmLSW7t7pf2ZuoAAADsd9tGaXfftMXwp95g/zuS3LGTSQEAALAednL1XQAAANgRUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMAYUQoAAMCYbaO0qu6qqueq6pFNY2+rqvuq6tvL7fmbtt1eVcer6omqunavJg4AAMD+dzqflH46yXWvGbstyf3dfTjJ/cvjVNXlSW5McsVyzCeq6pxdmy0AAAArZdso7e6vJPn+a4avT3L3cv/uJO/ZNH5Pd7/Q3U8mOZ7kql2aKwAAACvmTL9TelF3P5Mky+2Fy/jBJE9t2u/EMgYAAACvs9sXOqotxnrLHatuqapjVXXs5MmTuzwNAAAA9oMzjdJnq+riJFlun1vGTyS5dNN+lyR5eqsn6O47u/tIdx85cODAGU4DAACA/exMo/TeJDcv929O8sVN4zdW1XlVdVmSw0ke2NkUAQAAWFXnbrdDVX02yTVJLqiqE0k+muRjSY5W1fuTfC/JDUnS3Y9W1dEkjyV5Mcmt3f3SHs0dAACAfW7bKO3um06x6V2n2P+OJHfsZFIAAACsh92+0BEAAACcNlEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAGFEKAADAmHN3cnBVfTfJ80leSvJidx+pqrcl+TdJDiX5bpK/2d3/fWfTBAAAYBXtxielf6W7r+zuI8vj25Lc392Hk9y/PAYAAIDX2YvTd69Pcvdy/+4k79mD1wAAAGAF7DRKO8nvVdWDVXXLMnZRdz+TJMvthVsdWFW3VNWxqjp28uTJHU4DAACA/WhH3ylNcnV3P11VFya5r6q+dboHdvedSe5MkiNHjvQO5wEAAMA+tKNPSrv76eX2uSRfSHJVkmer6uIkWW6f2+kkAQAAWE1nHKVV9Weq6ideuZ/kryV5JMm9SW5edrs5yRd3OkkAAABW005O370oyReq6pXn+dfd/e+r6utJjlbV+5N8L8kNO58mAAAAq+iMo7S7v5PkZ7cY/29J3rWTSQEAALAe9uKfhAEAAIDTIkoBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYs2dRWlXXVdUTVXW8qm7bq9cBAABg/9qTKK2qc5L8ZpJfSHJ5kpuq6vK9eC0AAAD2r736pPSqJMe7+zvd/adJ7kly/R69FgAAAPvUXkXpwSRPbXp8Yhn7f6rqlqo6VlXHTp48uUfTAAAA4Gx27h49b20x1v/fg+47k9yZJFV1sqr+aI/msp0LkvzJ0Gszy9qvJ+u+vs6ytf97eXJ6CuvhLFv3s1PVVn912/es/Xqy7mevv3CqDXsVpSeSXLrp8SVJnj7Vzt19YI/msa2qOtbdR6ZenznWfj1Z9/Vl7deTdV9f1n49Wff9aa9O3/16ksNVdVlV/WiSG5Pcu0evBQAAwD61J5+UdveLVfWBJL+b5Jwkd3X3o3vxWgAAAOxfe3X6brr7S0m+tFfPv4vunJ4AY6z9erLu68varyfrvr6s/Xqy7vtQdff2ewEAAMAe2KvvlAIAAMC21ipKq+ofVdW3quqbVfWFqvqzm7bdXlXHq+qJqrp20/jPVdXDy7Z/Wit6zfRVVlU3VNWjVfVyVR15zTbrvkaq6rplrY9X1W3T82H3VNVdVfVcVT2yaextVXVfVX17uT1/07Yt3/vsL1V1aVV9uaoeX/6c/+Aybu1XXFX9WFU9UFX/eVn7f7iMW/s1UFXnVNV/qqrfXh5b931uraI0yX1Jfqa7/3KS/5Lk9iSpqsuzcYXgK5Jcl+QTVXXOcswnk9yS5PDyc90Pe9Ls2CNJ/kaSr2wetO7rZVnb30zyC0kuT3LT8jvAavh0Xv8+vS3J/d19OMn9y+Pt3vvsLy8m+XB3/3SSdya5dVlfa7/6XkjyV7v7Z5NcmeS6qnpnrP26+GCSxzc9tu773FpFaXf/Xne/uDz8ajb+/dQkuT7JPd39Qnc/meR4kquq6uIkP9ndf9gbX779l0ne80OfODvS3Y939xNbbLLu6+WqJMe7+zvd/adJ7snG7wAroLu/kuT7rxm+Psndy/278+r7eMv3/g9louyq7n6mu7+x3H8+G39JPRhrv/J6w/9cHr5l+elY+5VXVZck+etJ/sWmYeu+z61VlL7G30nyO8v9g0me2rTtxDJ2cLn/2nFWg3VfL6dab1bXRd39TLIRL0kuXMb9LqygqjqU5B1JvhZrvxaWUzgfSvJckvu629qvh3+S5O8neXnTmHXf5/bsn4SZUlW/n+TPbbHpV7v7i8s+v5qNU34+88phW+zfbzDOWeZ01n2rw7YYs+6ry7ryCr8LK6aq3prkc0k+1N0/eIPLAFj7FdLdLyW5crlGyBeq6mfeYHdrvwKq6heTPNfdD1bVNadzyBZj1v0stHJR2t0//0bbq+rmJL+Y5F396r+HcyLJpZt2uyTJ08v4JVuMc5bZbt1Pwbqvl1OtN6vr2aq6uLufWU7Lf24Z97uwQqrqLdkI0s909+eXYWu/Rrr7f1TVf8zGdwat/Wq7OskvVdW7k/xYkp+sqn8V677vrdXpu1V1XZJ/kOSXuvt/b9p0b5Ibq+q8qrosGxe2eWD5+P/5qnrncvXVv53kVJ+6sf9Y9/Xy9SSHq+qyqvrRbFz44N7hObG37k1y83L/5rz6Pt7yvT8wP3Zo+TP6U0ke7+6Pb9pk7VdcVR1YPiFNVf14kp9P8q1Y+5XW3bd39yXdfSgb/x//D939t2Ld972V+6R0G/88yXlJ7ltO7flqd/9Kdz9aVUeTPJaN03pvXU4JSZK/m42rOv54Nr6D+juve1bOalX13iT/LMmBJP+uqh7q7mut+3rp7her6gNJfjfJOUnu6u5Hh6fFLqmqzya5JskFVXUiyUeTfCzJ0ap6f5LvJbkhSbZ577O/XJ3kfUkeXr5bmCQfibVfBxcnuXu5kuqPJDna3b9dVX8Ya7+OvOf3uXr1DFYAAAD44Vqr03cBAAA4u4hSAAAAxohSAAAAxohSAAAAxohSAAAAxohSAAAAxohSAAAAxohSAAAAxvxfiSrKybwcBzgAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 1152x432 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.scalar_field(dh.cpu_arrays[\"c\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ur = ps.Assignment(c_field[0, 0], c_field[1, 0])\n",
+    "ast = ps.create_kernel(ur, target=dh.default_target, cpu_openmp=True)\n",
+    "kernel = ast.compile()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c_sync = dh.synchronization_function(['c'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def timeloop(steps=10):\n",
+    "    for i in range(steps):\n",
+    "        c_sync()\n",
+    "        dh.run_kernel(kernel)\n",
+    "    return dh.gather_array('c')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ps.jupyter.set_display_mode('video')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<video controls width=\"80%\">\n",
+       " <source src=\"data:video/x-m4v;base64,AAAAIGZ0eXBpc29tAAACAGlzb21pc28yYXZjMW1wNDEAAAAIZnJlZQAAG7ptZGF0AAACcgYF//9u3EXpvebZSLeWLNgg2SPu73gyNjQgLSBjb3JlIDE2MCByMzAxMSBjZGU5YTkzIC0gSC4yNjQvTVBFRy00IEFWQyBjb2RlYyAtIENvcHlsZWZ0IDIwMDMtMjAyMCAtIGh0dHA6Ly93d3cudmlkZW9sYW4ub3JnL3gyNjQuaHRtbCAtIG9wdGlvbnM6IGNhYmFjPTAgcmVmPTMgZGVibG9jaz0xOjA6MCBhbmFseXNlPTB4MToweDExMSBtZT1oZXggc3VibWU9NyBwc3k9MSBwc3lfcmQ9MS4wMDowLjAwIG1peGVkX3JlZj0xIG1lX3JhbmdlPTE2IGNocm9tYV9tZT0xIHRyZWxsaXM9MSA4eDhkY3Q9MCBjcW09MCBkZWFkem9uZT0yMSwxMSBmYXN0X3Bza2lwPTEgY2hyb21hX3FwX29mZnNldD0tMiB0aHJlYWRzPTEzIGxvb2thaGVhZF90aHJlYWRzPTIgc2xpY2VkX3RocmVhZHM9MCBucj0wIGRlY2ltYXRlPTEgaW50ZXJsYWNlZD0wIGJsdXJheV9jb21wYXQ9MCBjb25zdHJhaW5lZF9pbnRyYT0wIGJmcmFtZXM9MCB3ZWlnaHRwPTAga2V5aW50PTI1MCBrZXlpbnRfbWluPTI1IHNjZW5lY3V0PTQwIGludHJhX3JlZnJlc2g9MCByY19sb29rYWhlYWQ9NDAgcmM9Y3JmIG1idHJlZT0xIGNyZj0yMy4wIHFjb21wPTAuNjAgcXBtaW49MCBxcG1heD02OSBxcHN0ZXA9NCBpcF9yYXRpbz0xLjQwIGFxPTE6MS4wMACAAAAOrGWIhAvyYoAAqcScnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnJycnXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX/jpjbGB7AAzGAv7sCipqrccEuOuIt9wAQv19XnkMpGlotPm4GczKBOwyaRQEfbjkIzAXIjSREc//2ziZbN/UzdeETcuEXwGTTXk02oxGpWgc6S67htGL9CBekqo/X9B3HgDmoP/f+ceQBarBL+ZXP9oApPT4HbR3wS/wYCzUhqYDZhE+OOUh99X0hW1mMGaeH/8J8ABkkSBohw/FnngsMbX9xCvvvvvvvvvvvvvvv8f+HgsE4Cjp/jaBACYDCwrG0AVtB8IgdLCBIQUcc8YNbW1tbW18f4fQbH5CVedRSS22/+4+HhwMYQEAAQD8y9cICFNYNwdAWiHcQ6Ku50Ixi1tbW1tbW18f4fQbDOONHfzltvSS/8fDoOBLBZkwHEM8lucISAQ5cCEAeRtKWCCMcj4QDAB/oNgnwcTD0t3xAPQgQMiAcPLDg4eX4MuWUuX7Xx+AfsFYawIHG1W5BCoGcBoAX8mYqjmeQAjl3MJWNspK5uIV99999999999999/j7/DhQTwAI4v2wAECXNli4J6h+uuuuuuuuuuuuu+8NQTKoJxQDfffffffffffffffTBIp4IRQNHbW1tbW1rp5iRQACc89ra2tra2trXT0y10wRC7gjvvvvvvvvvvvvvvvrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrp6euuuuuunp666666666enrp666666666666666666666666666666666666666666666666666enrrrrrrp6euuuuuuuunp66euuuuuuuuuuuuuuuuuuuuuuuuuuuuuuv/HTEWSMD2ABmMBf3YFFTVVR01FmVSu1wAQv19XnkMpGlotPm4Gc2mBzagwiE0SX77TQmZgPkwsRXv+/0GQ0YX86mXeP/dRwobpwi+BRLtbXf7/BeISU9BJ0k11HsbmAJ/Q5/r3Vh2FgDmyh7vrTEEydhW/crm+0AKTnpeLRySFvBgLNSGpgNmET445SH31fSFbWYwZoO1/+/9Q/XXXXXXXXXXXXXXXT09ddddddPT1111111109PXT111111111111111111111111111111133+uuuuuuuuuuuuuuuuunp6666666enrrrrrrrrp6eunrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrp6euuuuuunp666666666enrp666666666666666666666666666666666666666666666666666enrrrrrrp6euuuuuuuunp66euuuuuuuuuuuuuuuuuuuuuuuuuuuuuuv/HHAw2HuAAz3E2SM4hL8aBofMpMNdk05TYgdM1HOOsbocABMnENF+XJ0UYZ0HZwlzlpclJL2iwm0IhC8pS+3Do4KeXMLTXH6A5fF/Tm8gsc9ERLgBWNkHaxOojivmneEHiFcJ7criCogc8z8a0eGsdL/95wwywRiM3px9qsAFzBB0zWuQh98H47//f9Q/XXXXXXXXXXXXXXXT09ddddddPT1111111109PXT111111111111111111111111111111133+uuuuuuuuuuuuuuuuunp6666666enrrrrrrrrp6eunrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrp6euuuuuunp666666666enrp666666666666666666666666666666666666666666666666666enrrrrrrp6euuuuuuuunp66euuuuuuuuuuuuuuuuuuuuuuuuuuuuuuv/HHCB2HuAAz3E2SM4hL8aBofMpMNdk00FFKQVMZHCJOv/Q4ACZOIaL8uToowy5oAAy0FYg6tFhNoRCF5Sl9tbxC+RKcqjZ/93X+5fF8IKNBIW90H/3gArGyDtYnURxX2x4ABmB/HGxXEFRA55f41Fhbocf/e+GGWCMRm9OPtVgAuYIOma1yEPvg/Hf/7/qH6666666666666666enrrrrrrp6euuuuuuuunp66euuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu+/1111111111111111109PXXXXXXT09ddddddddPT109ddddddddddddddddddddddddddddddddddddddddddddddddddPT11111109PXXXXXXXXT09dPXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXT09ddddddPT1111111109PXT1111111111111111111111111111111//9OFA5wAH2BII3/bz8KKNeaHBWfj3dyZvdKwGEXFxzEqfbq+ELjwi/xENxXu7t8i4ncCqOcORXxdQ9TgiL8yGCuq2IDf/842AAIAQsnEIOrSQjYjEhJnL/7tA65zBFT8cNnN/7wAXMCBUzWuQh99/rqHa666666666666666enrrrrrrp6euuuuuuuunp66euuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu+/1111111111111111109PXXXXXXT09ddddddddPT109ddddddddddddddddddddddddddddddddddddddddddddddddddPT11111109PXXXXXXXXT09dPXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXT09ddddddPT1111111109PXT11111111111111111111111111111111/fSEFyxYewCR+vYEAO+moM2jrG+9xhFKdbvreLYU6+zwYMGicVtpb01Wc153BGwYZMVjkoffB/D/AIBwQQEUOYDwvlNES50MP/wDALQQ4ACbINJO7InbrkESoGMqeu84fEf8a1r7TW+++++//xh7QqJS8ABNkGk3dkRW65BEqBjKnpvOfffffff+HGMAtBdAQ+udwAzI0Ri6kVHf7mIxnXcVLgU2i4OECQgBolQnOgPRIwhMfh0ISYtf/GH0Co9/6CQR1+i0eRfYTCz6fWsej/8mk3JpN/hx6AgQXBYACb8yp4Tv4xBDslr46NWRuBotWB0ZZZ3cmk06vla2totrC6OHct/pJaSW3KZKZyGSGTrePRgAIZum09J//gAZBfIpmyIOTZkG0UxnnQ234ceAcFwSBD7sHEEIy9KoOIAmCaBKctrtSPJrN/4AEA0C0EoOORy3xqhAKLLoyD0ZwcOMg4cZa4fwDANgtCwCG01Wrn4QmRAC4AGdftiMv4weHB9vv/8foMEMABJmEl2d0VSZzzloCNUDmVPn1v77777774/AAhhujMjyI//wAHmYSWZ3RVbB0QP0gd2UlBz33333/wx/CCHAA1YI5zGOcS7n376h+uuuuuuuuuuuuuu+v+rF4vxQewAEx+KUG+qzCF8T67ZNWnpwAPHQRFbYitH8sMaJR1SXVur33GQKrY5Uw1bv/91De2oXv+DMjAAECj7y1mdoa7/e+/3CfA9/cBTxlEYz1se6p3cSQHRVJShzGa3//KVj1Xu19TewNj/2k6UaQ3/9d8QiWk1D//dvxAg5BZXfddvE4kzHN+b/3/Xx38t/62/7h8UAAzfffff9wMMIAfYegA1uc5J+7wDgAHgI8TLWOMQboiTOve/AfOmYRbHTHEk3hsmxWFfqqRzBJeQQhBXWrA7FeYDRtO1taHUDslHLl9dY0xKlMUrDe1/1fTuuXieAkp8hCmRlFf/7sCAGsXmUEXx1WF/86F2Kxbt64+kpsUtl13251fHO9//vvt+bbN8VQcSJRWykMNw+KAAZvvvvv/Zg+ZaCgcAAcZIJn4160IQQJ/+4EjqQ0OxJinnnEVuixfjbLIHPI3ojqwDuZMtG7fvLNjUtiv7b9/kQdAAG1pa6WlpfgC8MH0wQAANnE0O5YlZekcNzFG6FsNkGNMeQ3juYolSwrJH5JtoFZZNKh/+8xoHLN9H6N/+mg7ds87v77AAhFrI5i6p+2W2WjbbBh+qiTIohK6LnFC8pI1++D/DDh/QTKAp4TXOclrCSMm8CxZWOav5+fzEdPJB1z/8/7QdtbW1vteLM6hQdX0wQAANuRoOz5qR+OC4lo1aHpZF38fkrek+DU99rGVs52/95DeKpiiRaiJ1eSbaBWWTSof/vNaByzfP0lv9bQdu2rzm/v9XMjF6zui99/VZNaMRO3h7glDRIiDp//74oXlI1/95/4fh/hMoCnhNc5yWsJIybwcWVjmr+fn6xHTyDr//9v/Mg6KANIktLS119YF0fVewQAAuWBiJvFahZu0D6AZCUJtNHF3IBGk1ifciqy4GtWnHWv5/X5//MbV2MQS7aVlN8Yk0dji6ikYqcTBST+xujf/VCFloG8ZdoZLDHH5LdV78zNmPv+3aB75xdWe5/uyUS5lX5hP+/6h+X+Ez6BJx85TXIV5cf3vlQ7PLWf+//1/mZjvX//8/6huuuuuvqwKOgbr2CAABZiJYV7VQgYvBMCNEFa2kJ7+/sWRABwwpUzH62EQbu6Jk3d785Q24ioTavQ+P/9gmYziLRznULb75mQ7EoG8pm49H9g/531NuCGpy9+bU37jZ0Nyq7fWqNEkGfPff8/5S/DxwRxeBOP0Jk3z70EXaB/h1Dddddddddddddddddddddddd9ddddddd9dddddd9dddddddddddddddd9dddddd99dddddd9ddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddeAAAAT5BmjgX4D5gLigomOBmnCi/MJQPu5MEAwAP6BWCPAdSctjgSIABGQHscsPY5fgwJ5ZQJ5fiICS8DRFAk5l5l74COwb4a8RIeCfFsFTvnYK8QuIXELiFxC4hcQtG+H/sPgmmX6+DfzHfOwTxJ318G/mO+d4k76+DfzHfO8Sd9fBv5jvneJO+vg38x3zvEnfXwb+Y753iTvr4N/Md87xJ318G/mO+d4k76+DfzHfO8Sd9fBv5jvnZYk76+DfzHfO8Sd9fBv5jvneJO+vg38x3zvEnfXwb+Y753iTvr4N/Md87xJ318G/mO+d4k76+DfzHfO8Sd9fBvuY753iTvr4PYGaCQFkAMe1PSyF8t6MUGZfwKFQM9Dc0S8E8gy98MYKIQxEIZ8MeEIjDHhDwh4QijwWxB+CrxPjhHECOCnxPgkgAAADQQZpUBXgPjAQ3gGZDgJr3QEfh76UVc4CAxGAgOEEGXaU0U1l0u/AzThhQ9b/HJP+Dw8Eud9fxQtgo9B6d+AIhijvB6d+4o7wenfuKO8Hp37ijvB6d+4o7wenfuKO8Hp37ijvB6d+4o7wenfuKO8Hp37ijvB6d+4o7wenfuKO8Hp37ijvB6d+4o7wenfuKO8Hp37ijvB6d+4o7wenfuKO8HXoJ+t2BFhoIXugE/Birnu+f0EXvAzTjlByT5Fvg4EoFPxEt3xnxnxnxkCL8CR8EkAAAAOFBmmAnwHvgIDWBmnBM0q0Dr7n4jASHGIMnJePAzThiT/KJQdT8YuoPTB//2CuKi7PzvnfFD2SpXyvrXOwXweneKO+d4PTvFHfO8Hp3ijvneD07xR3zvB6d4o753g9O8Ud87weneKO+d4PTvFHfO8Hp3ijvneD07xR3zvB6d4o753g9O8Ud87weneKO+d4PTvFHfO8Hp3ijvneD07xR3zvB6d4o753g9O8Ud87wdhBBTH7CYWfT6sDNOFkiA6+5/xRfAzfw5AAwvJTekwcyoaTkfPk/wOqufjF1BwJQK5r4C3gAAAEdQZqAL8B7YCA1QDMoFjdc+AQD4CAoN3NrgZuIWDw8fugM/QKm65/pr+eC+D30C5mdgrn+jvnlg9wM1AuZ0A0M/0d88F8HvoFzOuf6O+eC+D30C5nXP9HfPBfB76Bczrn+jvngvg99AuZ1z/R3zwXwe+gXM65/o754L4PfQLmdc/0d88F8HvoFzOuf6O+eC+D30C5nXP9HfPBfB76Bczrn+jvngvg99AuZ1z/R3zwXwe+gXM65/o754L4PfQLmdc/0d88F8HvoFzOuf6O+eC+D30C5nXP9HfPBfB76Bczrn+jvngvg99AuZnYL5/o75+DvzhSL+fT6i//ghC0i8i6Qb/n+ieqeBm4BMO9cwvDMvPB7XECECGq4PvgSPgpgAAAArkGaoC/AeuAgNOsBneufAR2rAJl8DhRiJB4fxbBVfEIFsTfB8d8/E3wfHfPxN8Hx3z8TfB8d8/E3wfHfPxN8Hx3z8TfB8d8/E3wfHfPxN8Hx3z8TfB8d88sTfB8d8/E3wfHfPxN8Hx3z8TfB8d8/E3wfHfPxN8Hx3z8TfB8d8/E3wfHfPxN8HgQQYkX2Ewsblq5adYDO9c9gFw0784UWZicP/oDxQJQQ0H4hCYC4gAAAAQpBmsAvwHpgIDWBmnBUl/MJQObuaAZlAjQYTnwEBov4GjgZooEXHUFy8y/AQGD0Wya34GjEhNl/e/3vX8Hp3zsE8Sd9fwenfO8Sd87LB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB6d87xJ3zvB5YBdPgZpwVX/A5u55g/gAQDQLQQg07R8RA0IDlgczlhzOXwcoyDlGT0AQjWBmigRTLzLw7J/+AgMHh4LfkivgSPgSPghgAAAJFBmuArwHrQDMoEz+B+z4CI3gZp6eWcHp8wlA+7mD2vxnEQSxJ3zyweV0fiTvB7XFneD2uLO8HtcWd4Pa4s7we1xZ3g9rizvB7XFneD2uLO8HtcWd4Pa4s7we1xZ3g9rizvB7XFneD2uLO8HtcWd4Pa4s7we1xZ3g8r6/LnwENvAzTgoTwen8Hp/B78CR8CR8FkAAAAsUGbACfAetgEyziFuwCZUChP3we3ntO+d/xQ98pXyvg6nwdT++D0e9nzvnfFHfvg9O8Ud++D07xR374PTvFHfvg9O8Ud++D07xR374PTvFHfvg9O8Ud++D07xR374PTvFHfvg9O8Ud++D07xR374PTvFHfvg9O8Ud++D07xR374PTvFHfvg9O8Ud++D07xR3zsFsHuBmnBcoHX3P+fAQGgggzY6af7AJlOFnfB1P++A+YAAAAOFBmyArwHpQDM9cS7wCZe+D2vriXfoFz9cHteeC2Jd+gWP1we154LYl36BY+dgvg9rz8S79AufOwXwe15+Jd+gXPnYL4Pa8/Eu/QLnzsF8HtefiXfoFz52C+D2vPxLv0C587BfB7Xn4l36Bc+dgvg9rz8S79AufOwXwe15+Jd+gXPnYMYPa8/Eu/QMHzsGMHtefiXfoGD52DGD2vPxLv0DB87BjB7Xn4l36Bg+dgxg9rz8S79AwfOwYwe15+Jd+gYPnYMYPa8/Eu/QMHzsGMHtefiXfoGD98HteIRYl3/XAfUAAAAGfQZtAJ8B5vwKULgowDb5JTBy0SZg5cSZhWNvAUZrkeFD+BAxOBohomX1lEoOKY/gXMSbH/w0HC8AYr6Jk2lM4EqLhI5tKzuRzpLuTf/+w3wzURPAPLteC8knqLn2QgzHsln54Twr+JN8of7JHsGvrq6+JFsFkvzsgfpgMzvnYKYk753iTvneDM753iTvneJO+d4MzvneJO+d4k753gzO+d4k753iTvneDM753iTvneJO+d4MzvneJO+d4k753gzO+d4k753iTvneDM753iTvneJO+d4MzvneJO+d4k753gzO+d4k753iTvneDM753iTvneJO+d4MzvneJO+d4k753gzO+d4k753iTvneDM753iTvneJO+d4MzvneJO+d4k753gzO+d4k753iTvneDM753iTvneJO+d4MzvneJO+d4k753gzN4FjwCEFwKAAcuNGJWXV4Bk016ye+aZHMrOxXJt6fAgYnA0TkWDjMf8CxQf6JN8MfDQLg1ARK0498BZcGzuTX/UK7nm/wDANgth9SBw+rV/173+8HvwJHwJHwSQAAAANpBm2AnwH5gZqBY/gfMUbH4cOw3wCUegnHNL98LjwbigFfk9eHrIRnuZWWxm5g99HS5f/4oWyPKD0Nsv73+91+aE0J2CeKO8Hp3zvFHeD0753ijvB6d87xR3g9O+d4o7wenfO8Ud4PTvneKO8Hp3zvFHeD0753ijvB6d87xR3g9O+d4o7wenfO8Ud4PTvneKO8Hp3zvFHeD0753ijvB6d87xR3g9O+d4o7wenfO8Ud4PaAZ0gLJl5l9LFGD+AYBsFsBKkZG+b/hMsS1gJ1ExPvNNycmg1xmFK4CpgAAAy9tb292AAAAbG12aGQAAAAAAAAAAAAAAAAAAAPoAAABkAABAAABAAAAAAAAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAAAAAAAAAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAACWXRyYWsAAABcdGtoZAAAAAMAAAAAAAAAAAAAAAEAAAAAAAABkAAAAAAAAAAAAAAAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAAAAAAAAAEAAAAAEgAAAAbAAAAAAACRlZHRzAAAAHGVsc3QAAAAAAAAAAQAAAZAAAAAAAAEAAAAAAdFtZGlhAAAAIG1kaGQAAAAAAAAAAAAAAAAAADwAAAAYAFXEAAAAAAAtaGRscgAAAAAAAAAAdmlkZQAAAAAAAAAAAAAAAFZpZGVvSGFuZGxlcgAAAAF8bWluZgAAABR2bWhkAAAAAQAAAAAAAAAAAAAAJGRpbmYAAAAcZHJlZgAAAAAAAAABAAAADHVybCAAAAABAAABPHN0YmwAAACUc3RzZAAAAAAAAAABAAAAhGF2YzEAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAAEgAGwAEgAAABIAAAAAAAAAAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAY//8AAAAuYXZjQwFCwB7/4QAWZ0LAHtkASA3oQAAAAwBAAAAPA8WLkgEABWjLg8sgAAAAGHN0dHMAAAAAAAAAAQAAAAwAAAIAAAAAFHN0c3MAAAAAAAAAAQAAAAEAAAAcc3RzYwAAAAAAAAABAAAAAQAAAAwAAAABAAAARHN0c3oAAAAAAAAAAAAAAAwAABEmAAABQgAAANQAAADlAAABIQAAALIAAAEOAAAAlQAAALUAAADlAAABowAAAN4AAAAUc3RjbwAAAAAAAAABAAAAMAAAAGJ1ZHRhAAAAWm1ldGEAAAAAAAAAIWhkbHIAAAAAAAAAAG1kaXJhcHBsAAAAAAAAAAAAAAAALWlsc3QAAAAlqXRvbwAAAB1kYXRhAAAAAQAAAABMYXZmNTguNDUuMTAw\" type=\"video/mp4\">\n",
+       " Your browser does not support the video tag.\n",
+       "</video>"
+      ],
+      "text/plain": [
+       "<IPython.core.display.HTML object>"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ani = ps.plot.scalar_field_animation(timeloop, rescale=True, frames=12)\n",
+    "ps.jupyter.display_animation(ani)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ps.jupyter.set_display_mode('image_update')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAFlCAYAAAATVk7bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAUlUlEQVR4nO3db4xl9X3f8c/X4JCoTlosFkoX1CXRRgpEDVZWyBVSReu0YDcKdiWiRaqDVEukFZZsyVILzgOnD5AstXGqtrElUiNT1TUlsi2j1GmCiSsrUmK8uNQG1sQrg82GFWzqRqaqhLvw7YM5lCnM/mFnhu/O3NdLGt17f+ece3/oNzPsW/fcM9XdAQAAgAlvmp4AAAAAq0uUAgAAMEaUAgAAMEaUAgAAMEaUAgAAMEaUAgAAMOb86QkkyUUXXdT79u2bngawAv704e9MTyF/+aoTueS8H05PA/LkD9+S//Otl6anwTngp3/+J6enAOxyDz/88J93956Ntp0TUbpv374cOnRoehrACvi7b7ppegp55+/8RT544VPT04D8ynf/Vp79mz+YngbngAcO/c70FIBdrqq+e7JtTt8FAABgzGmjtKour6ovV9Xhqnqsqj6wjP96Vf1ZVT2yfL1r3TF3VNWRqnqiqq7fzv8AAAAAdq4zOX33RJIPdffXq+rHkzxcVQ8s236zu//l+p2r6sokB5NcleSvJflSVf10d7+4lRMHAABg5zvtO6Xdfay7v77cfz7J4SR7T3HIjUnu7e4XuvvJJEeSXLMVkwUAAGB3eV2fKa2qfUneluSry9D7q+obVXV3VV24jO1N8vS6w47m1BELAADAijrjKK2qtyT5bJIPdvcPknwiyU8luTrJsSS/8fKuGxzeGzzfrVV1qKoOHT9+/HVPHAAAgJ3vjKK0qt6ctSD9dHd/Lkm6+9nufrG7X0ry23nlFN2jSS5fd/hlSZ559XN2913dfaC7D+zZs+GfqwEAAGCXO5Or71aSTyY53N0fWzd+6brd3pPk0eX+/UkOVtUFVXVFkv1JHtq6KQMAALBbnMnVd69N8t4k36yqR5axDye5uaquztqpuU8l+dUk6e7Hquq+JI9n7cq9t7nyLgAAABs5bZR29x9l48+JfvEUx9yZ5M5NzAsAAIAV8LquvgsAAABbSZQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAwRpQCAAAw5rRRWlWXV9WXq+pwVT1WVR9Yxt9aVQ9U1beX2wvXHXNHVR2pqieq6vrt/A8AAABg5zqTd0pPJPlQd/9Mkrcnua2qrkxye5IHu3t/kgeXx1m2HUxyVZIbkny8qs7bjskDAACws502Srv7WHd/fbn/fJLDSfYmuTHJPctu9yR593L/xiT3dvcL3f1kkiNJrtnqiQMAALDzva7PlFbVviRvS/LVJJd097FkLVyTXLzstjfJ0+sOO7qMvfq5bq2qQ1V16Pjx469/5gAAAOx4ZxylVfWWJJ9N8sHu/sGpdt1grF8z0H1Xdx/o7gN79uw502kAAACwi5xRlFbVm7MWpJ/u7s8tw89W1aXL9kuTPLeMH01y+brDL0vyzNZMFwAAgN3kTK6+W0k+meRwd39s3ab7k9yy3L8lyRfWjR+sqguq6ook+5M8tHVTBgAAYLc4/wz2uTbJe5N8s6oeWcY+nOSjSe6rqvcl+V6Sm5Kkux+rqvuSPJ61K/fe1t0vbvnMAQAA2PFOG6Xd/UfZ+HOiSfKOkxxzZ5I7NzEvAAAAVsDruvouAAAAbCVRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwJjTRmlV3V1Vz1XVo+vGfr2q/qyqHlm+3rVu2x1VdaSqnqiq67dr4gAAAOx8Z/JO6aeS3LDB+G9299XL1xeTpKquTHIwyVXLMR+vqvO2arIAAADsLqeN0u7+SpLvn+Hz3Zjk3u5+obufTHIkyTWbmB8AAAC72GY+U/r+qvrGcnrvhcvY3iRPr9vn6DL2GlV1a1UdqqpDx48f38Q0AAAA2KnONko/keSnklyd5FiS31jGa4N9e6Mn6O67uvtAdx/Ys2fPWU4DAACAneysorS7n+3uF7v7pSS/nVdO0T2a5PJ1u16W5JnNTREAAIDd6qyitKouXffwPUlevjLv/UkOVtUFVXVFkv1JHtrcFAEAANitzj/dDlX1mSTXJbmoqo4m+UiS66rq6qydmvtUkl9Nku5+rKruS/J4khNJbuvuF7dn6gAAAOx0p43S7r55g+FPnmL/O5PcuZlJAQAAsBo2c/VdAAAA2BRRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwBhRCgAAwJjTRmlV3V1Vz1XVo+vG3lpVD1TVt5fbC9dtu6OqjlTVE1V1/XZNHAAAgJ3vTN4p/VSSG141dnuSB7t7f5IHl8epqiuTHExy1XLMx6vqvC2bLQAAALvKaaO0u7+S5PuvGr4xyT3L/XuSvHvd+L3d/UJ3P5nkSJJrtmiuAAAA7DJn+5nSS7r7WJIstxcv43uTPL1uv6PLGAAAALzGVl/oqDYY6w13rLq1qg5V1aHjx49v8TQAAADYCc42Sp+tqkuTZLl9bhk/muTydftdluSZjZ6gu+/q7gPdfWDPnj1nOQ0AAAB2srON0vuT3LLcvyXJF9aNH6yqC6rqiiT7kzy0uSkCAACwW51/uh2q6jNJrktyUVUdTfKRJB9Ncl9VvS/J95LclCTd/VhV3Zfk8SQnktzW3S9u09wBAADY4U4bpd1980k2veMk+9+Z5M7NTAoAAIDVsNUXOgIAAIAzJkoBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYI0oBAAAYc/5mDq6qp5I8n+TFJCe6+0BVvTXJf0qyL8lTSX65u//n5qYJAADAbrQV75T+7e6+ursPLI9vT/Jgd+9P8uDyGAAAAF5jO07fvTHJPcv9e5K8exteAwAAgF1gs1HaSf6gqh6uqluXsUu6+1iSLLcXb3RgVd1aVYeq6tDx48c3OQ0AAAB2ok19pjTJtd39TFVdnOSBqvrWmR7Y3XcluStJDhw40JucBwAAADvQpt4p7e5nltvnknw+yTVJnq2qS5NkuX1us5MEAABgdzrrKK2qv1RVP/7y/SR/L8mjSe5Pcsuy2y1JvrDZSQIAALA7beb03UuSfL6qXn6e/9jd/6Wqvpbkvqp6X5LvJblp89MEAABgNzrrKO3u7yT5uQ3G/0eSd2xmUgAAAKyG7fiTMAAAAHBGRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjRCkAAABjti1Kq+qGqnqiqo5U1e3b9ToAAADsXNsSpVV1XpLfSvLOJFcmubmqrtyO1wIAAGDn2q53Sq9JcqS7v9PdP0xyb5Ibt+m1AAAA2KG2K0r3Jnl63eOjy9j/U1W3VtWhqjp0/PjxbZoGAAAA57Lzt+l5a4Ox/v8edN+V5K4kqarjVfXdbZrLyVyU5M/f4Nfk3GH9V9f42n/pyuRDkxNYbePrf245Mj2BN5K1P4Wqjf7ptqtY/9Vl7c8df/1kG7YrSo8muXzd48uSPHOynbt7zzbN46Sq6lB3H3ijX5dzg/VfXdZ+tVn/1WXtV5v1X13WfmfYrtN3v5Zkf1VdUVU/kuRgkvu36bUAAADYobblndLuPlFV70/y+0nOS3J3dz+2Ha8FAADAzrVdp++mu7+Y5Ivb9fxb4K7pCTDK+q8ua7/arP/qsvarzfqvLmu/A1R3n34vAAAA2Abb9ZlSAAAAOK2ViNKq+hdV9a2q+kZVfb6q/sq6bXdU1ZGqeqKqrl83/vNV9c1l27+uFbhW+m5UVTdV1WNV9VJVHXjVNmu/YqrqhmW9j1TV7dPzYWtV1d1V9VxVPbpu7K1V9UBVfXu5vXDdtg1/B7DzVNXlVfXlqjq8/M7/wDJu/VdAVf1oVT1UVf99Wf9/voxb/xVRVedV1X+rqt9dHlv7HWYlojTJA0l+trv/RpI/TXJHklTVlVm7MvBVSW5I8vGqOm855hNJbk2yf/m64Y2eNFvi0ST/IMlX1g9a+9WzrO9vJXlnkiuT3Lx8H7B7fCqv/Xm9PcmD3b0/yYPL49P9DmDnOZHkQ939M0nenuS2ZY2t/2p4Icnf6e6fS3J1khuq6u2x/qvkA0kOr3ts7XeYlYjS7v6D7j6xPPyTrP3d1CS5Mcm93f1Cdz+Ztb8ifk1VXZrkJ7r7j3vtQ7f/Psm73/CJs2ndfbi7n9hgk7VfPdckOdLd3+nuHya5N2vfB+wS3f2VJN9/1fCNSe5Z7t+TV36eN/wd8IZMlC3X3ce6++vL/eez9o/TvbH+K6HX/K/l4ZuXr471XwlVdVmSv5/k360btvY7zEpE6av8oyS/t9zfm+TpdduOLmN7l/uvHmf3sPar52Rrzu52SXcfS9bCJcnFy7jvh12qqvYleVuSr8b6r4zl9M1HkjyX5IHutv6r418l+adJXlo3Zu13mG37kzBvtKr6UpK/usGmX+vuLyz7/FrWTvH59MuHbbB/n2Kcc9CZrP1Gh20wZu13N2vLer4fdqGqekuSzyb5YHf/4BSXBLD+u0x3v5jk6uW6IZ+vqp89xe7Wf5eoql9M8lx3P1xV153JIRuMWftzwK6J0u7+hVNtr6pbkvxiknf0K38H52iSy9ftdlmSZ5bxyzYY5xx0urU/CWu/ek625uxuz1bVpd19bDk9/7ll3PfDLlNVb85akH66uz+3DFv/FdPdf1FV/zVrnxe0/rvftUl+qareleRHk/xEVf2HWPsdZyVO362qG5L8syS/1N3/e92m+5McrKoLquqKrF3U5qHlbf7nq+rty5VXfyXJyd5xY2ey9qvna0n2V9UVVfUjWbvQwf3Dc2L73Z/kluX+LXnl53nD3wED82MLLL+vP5nkcHd/bN0m678CqmrP8g5pqurHkvxCkm/F+u963X1Hd1/W3fuy9v/1P+zufxhrv+PsmndKT+PfJrkgyQPLqTx/0t3/uLsfq6r7kjyetdN6b1tO/0iSf5K1Kzn+WNY+g/p7r3lWznlV9Z4k/ybJniT/uaoe6e7rrf3q6e4TVfX+JL+f5Lwkd3f3Y8PTYgtV1WeSXJfkoqo6muQjST6a5L6qel+S7yW5KUlO8zuAnefaJO9N8s3lc4VJ8uFY/1VxaZJ7lquovinJfd39u1X1x7H+q8rP/g5Tr5zJCgAAAG+slTh9FwAAgHOTKAUAAGCMKAUAAGCMKAUAAGCMKAUAAGCMKAUAAGCMKAUAAGCMKAUAAGDM/wVMdtKBi/7dVQAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 1152x432 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "ani = ps.plot.scalar_field_animation(timeloop, rescale=True, frames=12)\n",
+    "ps.jupyter.display_animation(ani)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def grid_update_function(image):\n",
+    "    for i in range(40):\n",
+    "        c_sync()\n",
+    "        dh.run_kernel(kernel)\n",
+    "    return dh.gather_array('c')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAFoCAYAAAB3+xGSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAPSklEQVR4nO3dX4zl5V3H8c9XaGmsmEBbCAKxqGuUJhbrZmuCMTWNgtxQL2q2F4aLJtsLmrSJXlC9sDckamx71yZrJBKjRRJt4IKIuDFpvIFuG6T8EVhbUrYQtoJJG5tgWb9ezG/DCDO7w87Mzn73vF7J5Jx5zu/MeZ6c4Z2zvzmHp7o7AMzxY3s9AQDeGuEGGEa4AYYRboBhhBtgGOEGGGbXwl1VN1fV01V1rKru2K3HAVg1tRvv466qi5I8k+Q3kxxP8rUkH+vuJ3f8wQBWzG694j6Q5Fh3f6u7/yfJPUlu3aXHAlgpF+/Sz706yfPrvj+e5IPrD6iqQ0kOJck7f7x+5Rd+7u27NBWAeZ57/kf5z1dO1ka37Va4N3qw/3dOprsPJzmcJPvf/45+5MFrd2kqAPMcuOn5TW/brVMlx5OsL/E1SV7YpccCWCm7Fe6vJdlXVddV1duTHExy/y49FsBK2ZVTJd39WlV9MsmDSS5Kcld3P7EbjwWwanbrHHe6+4EkD+zWzwdYVT45CTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMc/F27lxVzyX5QZKTSV7r7v1VdXmSv0vy3iTPJfnd7v6v7U0TgFN24hX3b3T3Dd29f/n+jiRHuntfkiPL9wDskN04VXJrkruX63cn+cguPAbAytpuuDvJP1XV16vq0DJ2ZXe/mCTL5RUb3bGqDlXV0ao6+r2XT25zGgCrY1vnuJPc2N0vVNUVSR6qqn/f6h27+3CSw0my//3v6G3OA2BlbOsVd3e/sFyeSPKVJAeSvFRVVyXJcnliu5ME4HVnHe6qemdVXXrqepLfSvJ4kvuT3LYcdluS+7Y7SQBet51TJVcm+UpVnfo5f9vd/1hVX0tyb1V9PMl3knx0+9ME4JSzDnd3fyvJ+zcYfznJh7czKQA255OTAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wzBnDXVV3VdWJqnp83djlVfVQVT27XF627rbPVNWxqnq6qm7arYkDrKqtvOL+qyQ3v2HsjiRHuntfkiPL96mq65McTPK+5T5frKqLdmy2AJw53N391SSvvGH41iR3L9fvTvKRdeP3dPer3f3tJMeSHNihuQKQsz/HfWV3v5gky+UVy/jVSZ5fd9zxZexNqupQVR2tqqPfe/nkWU4DYPXs9B8na4Ox3ujA7j7c3fu7e/973uVsCsBWnW24X6qqq5JkuTyxjB9Pcu26465J8sLZTw+ANzrbcN+f5Lbl+m1J7ls3frCqLqmq65LsS/LI9qYIwHoXn+mAqvpykg8leXdVHU/yx0n+JMm9VfXxJN9J8tEk6e4nqureJE8meS3J7d3tBDbADjpjuLv7Y5vc9OFNjr8zyZ3bmRQAm/PJSYBhhBtgGOEGGEa4AYYRboBhhBtgGOEGGEa4AYYRboBhhBtgGOEGGEa4AYY54/9k6lx45rEfz00/dcNeTwPgvPFMv7zpbV5xAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwxzxnBX1V1VdaKqHl839tmq+m5VPbp83bLuts9U1bGqerqqbtqtiQOsqq284v6rJDdvMP6F7r5h+XogSarq+iQHk7xvuc8Xq+qinZosAFsId3d/NckrW/x5tya5p7tf7e5vJzmW5MA25gfAG2znHPcnq+qx5VTKZcvY1UmeX3fM8WXsTarqUFUdraqjP8qr25gGwGo523B/KcnPJrkhyYtJPreM1wbH9kY/oLsPd/f+7t7/tlxyltMAWD1nFe7ufqm7T3b3/yb5i7x+OuR4kmvXHXpNkhe2N0UA1jurcFfVVeu+/Z0kp95xcn+Sg1V1SVVdl2Rfkke2N0UA1rv4TAdU1ZeTfCjJu6vqeJI/TvKhqroha6dBnkvyiSTp7ieq6t4kTyZ5Lcnt3X1yd6YOsJqqe8NT0OfUT9bl/cH68F5PA+C88XAfyff7lY3+buiTkwDTCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMGcMd1VdW1X/UlVPVdUTVfWpZfzyqnqoqp5dLi9bd5/PVNWxqnq6qm7azQUArJqtvOJ+Lcnvd/cvJvnVJLdX1fVJ7khypLv3JTmyfJ/ltoNJ3pfk5iRfrKqLdmPyAKvojOHu7he7+xvL9R8keSrJ1UluTXL3ctjdST6yXL81yT3d/Wp3fzvJsSQHdnriAKvqLZ3jrqr3JvnlJA8nubK7X0zW4p7kiuWwq5M8v+5ux5exN/6sQ1V1tKqO/iivvvWZA6yoLYe7qn4iyd8n+XR3f/90h24w1m8a6D7c3fu7e//bcslWpwGw8rYU7qp6W9ai/Tfd/Q/L8EtVddVy+1VJTizjx5Ncu+7u1yR5YWemC8BW3lVSSf4yyVPd/fl1N92f5Lbl+m1J7ls3frCqLqmq65LsS/LIzk0ZYLVdvIVjbkzye0m+WVWPLmN/mORPktxbVR9P8p0kH02S7n6iqu5N8mTW3pFye3ef3PGZA6yo6n7T6edz7ifr8v5gfXivpwFw3ni4j+T7/cpGfzP0yUmAaYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYS7e6wkkyc//0g/z4IOP7vU0AM4bB2764aa3ecUNMIxwAwwj3ADDCDfAMGcMd1VdW1X/UlVPVdUTVfWpZfyzVfXdqnp0+bpl3X0+U1XHqurpqrppNxcAsGq28q6S15L8fnd/o6ouTfL1qnpoue0L3f3n6w+uquuTHEzyviQ/leSfq+rnu/vkTk4cYFWd8RV3d7/Y3d9Yrv8gyVNJrj7NXW5Nck93v9rd305yLMmBnZgsAG/xHHdVvTfJLyd5eBn6ZFU9VlV3VdVly9jVSZ5fd7fj2SD0VXWoqo5W1dHvvezFOMBWbTncVfUTSf4+yae7+/tJvpTkZ5PckOTFJJ87degGd+83DXQf7u793b3/Pe+66C1PHGBVbSncVfW2rEX7b7r7H5Kku1/q7pPd/b9J/iKvnw45nuTadXe/JskLOzdlgNW2lXeVVJK/TPJUd39+3fhV6w77nSSPL9fvT3Kwqi6pquuS7EvyyM5NGWC1beVdJTcm+b0k36yqU/9DkT9M8rGquiFrp0GeS/KJJOnuJ6rq3iRPZu0dKbd7RwnAzjljuLv7X7PxeesHTnOfO5PcuY15AbAJn5wEGEa4AYYRboBhhBtgmOp+02djzv0kqr6X5L+T/Odez2WPvDvWvqpWef3Wfno/3d3v2eiG8yLcSVJVR7t7/17PYy9Y+2quPVnt9Vv72a/dqRKAYYQbYJjzKdyH93oCe8jaV9cqr9/az9J5c44bgK05n15xA7AFwg0wzJ6Hu6puXjYVPlZVd+z1fM6Fqnquqr65bLJ8dBm7vKoeqqpnl8vLzvRzJlh2RzpRVY+vG9t0rRfSRtObrH0lNtk+zSbjF/xzf042WO/uPftKclGS/0jyM0nenuTfkly/l3M6R+t+Lsm73zD2Z0nuWK7fkeRP93qeO7TWX0/ygSSPn2mtSa5ffgcuSXLd8rtx0V6vYYfX/tkkf7DBsRfa2q9K8oHl+qVJnlnWeME/96dZ+44993v9ivtAkmPd/a3u/p8k92Rts+FVdGuSu5frdyf5yB7OZcd091eTvPKG4c3WekFtNL3J2jdzoa19s03GL/jn/jRr38xbXvteh3tLGwtfgDrJP1XV16vq0DJ2ZXe/mKw98Umu2LPZ7b7N1roqvw9nvcn2RG/YZHylnvud3GB9vb0O95Y2Fr4A3djdH0jy20lur6pf3+sJnSdW4fdhW5tsT7PBJuObHrrB2Oj17/QG6+vtdbhXcmPh7n5huTyR5CtZ+2fRS6f28VwuT+zdDHfdZmu94H8feoU22d5ok/GsyHO/2xus73W4v5ZkX1VdV1VvT3Iwa5sNX7Cq6p1Vdemp60l+K2sbLd+f5LblsNuS3Lc3MzwnNlvrBb/R9Kpssr3ZJuNZgef+nGywfh78BfaWrP3V9T+S/NFez+ccrPdnsvYX5H9L8sSpNSd5V5IjSZ5dLi/f67nu0Hq/nLV/Fv4oa68sPn66tSb5o+V34ekkv73X89+Ftf91km8meWz5D/aqC3Ttv5a1f+4/luTR5euWVXjuT7P2HXvufeQdYJi9PlUCwFsk3ADDCDfAMMINMIxwAwwj3ADDCDfAMP8HxfXmalFAB3UAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 1152x432 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "animation = ps.jupyter.make_imshow_animation(dh.cpu_arrays[\"c\"], grid_update_function, frames=300)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/pystencils_tests/test_jacobi_llvm.py b/pystencils_tests/test_llvm.py
similarity index 72%
rename from pystencils_tests/test_jacobi_llvm.py
rename to pystencils_tests/test_llvm.py
index 953432c11..54c60ff45 100644
--- a/pystencils_tests/test_jacobi_llvm.py
+++ b/pystencils_tests/test_llvm.py
@@ -6,11 +6,11 @@ try:
     from pystencils.cpu.cpujit import get_llc_command
     from pystencils import Assignment, Field, show_code
     import numpy as np
+    import sympy as sp
 except ModuleNotFoundError:
     pytest.importorskip("llvmlite")
 
 
-
 def test_jacobi_fixed_field_size():
     size = (30, 20)
 
@@ -93,5 +93,50 @@ def test_jacobi_variable_field_size():
     np.testing.assert_almost_equal(error, 0.0)
 
 
+def test_pow_llvm():
+    size = (30, 20)
+
+    src_field_llvm = 4 * np.ones(size)
+    dst_field_llvm = np.zeros(size)
+
+    f = Field.create_from_numpy_array("f", src_field_llvm)
+    d = Field.create_from_numpy_array("d", dst_field_llvm)
+
+    ur = Assignment(d[0, 0], sp.Pow(f[0, 0], -1.0))
+    ast = create_kernel([ur])
+
+    jit = generate_and_jit(ast)
+    jit('kernel', dst_field_llvm, src_field_llvm)
+    assert np.all(0.25 == dst_field_llvm)
+
+    ur = Assignment(d[0, 0], sp.Pow(f[0, 0], 0.5))
+    ast = create_kernel([ur])
+
+    jit = generate_and_jit(ast)
+    jit('kernel', dst_field_llvm, src_field_llvm)
+    assert np.all(2.0 == dst_field_llvm)
+
+    ur = Assignment(d[0, 0], sp.Pow(f[0, 0], 2.0))
+    ast = create_kernel([ur])
+
+    jit = generate_and_jit(ast)
+    jit('kernel', dst_field_llvm, src_field_llvm)
+    assert np.all(16.0 == dst_field_llvm)
+
+    ur = Assignment(d[0, 0], sp.Pow(f[0, 0], 3.0))
+    ast = create_kernel([ur])
+
+    jit = generate_and_jit(ast)
+    jit('kernel', dst_field_llvm, src_field_llvm)
+    assert np.all(64.0 == dst_field_llvm)
+
+    ur = Assignment(d[0, 0], sp.Pow(f[0, 0], 4.0))
+    ast = create_kernel([ur])
+
+    jit = generate_and_jit(ast)
+    jit('kernel', dst_field_llvm, src_field_llvm)
+    assert np.all(256.0 == dst_field_llvm)
+
+
 if __name__ == "__main__":
     test_jacobi_fixed_field_size_gpu()
diff --git a/pystencils_tests/test_simplification_strategy.py b/pystencils_tests/test_simplification_strategy.py
index b160f4320..19ade3ddc 100644
--- a/pystencils_tests/test_simplification_strategy.py
+++ b/pystencils_tests/test_simplification_strategy.py
@@ -73,6 +73,11 @@ def test_split_inner_loop():
     code = ps.get_code_str(ast)
     # we have four inner loops as indicated in split groups (4 elements) plus one outer loop
     assert code.count('for') == 5
+    ast = ps.create_kernel(ac, target='gpu')
+
+    code = ps.get_code_str(ast)
+    # on GPUs is wouldn't be good to use loop splitting
+    assert code.count('for') == 0
 
     ac = AssignmentCollection(main, subexpressions)
     ast = ps.create_kernel(ac)
diff --git a/pystencils_tests/test_stencils.py b/pystencils_tests/test_stencils.py
index fa205ad7f..e66971316 100644
--- a/pystencils_tests/test_stencils.py
+++ b/pystencils_tests/test_stencils.py
@@ -27,4 +27,4 @@ def test_coefficient_list():
 
 def test_plot_expression():
     f = ps.fields("f: double[2D]")
-    plot_expression(2 * f[1, 0] + 3 * f[0, -1], matrix_form=True)
\ No newline at end of file
+    plot_expression(2 * f[1, 0] + 3 * f[0, -1], matrix_form=True)
diff --git a/pytest.ini b/pytest.ini
index 5cdf16be1..0cf6fbee5 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -27,10 +27,12 @@ exclude_lines =
        pragma: no cover
 
        def __repr__
+       def _repr_html_
 
        # Don't complain if tests don't hit defensive assertion code:
        raise AssertionError
        raise NotImplementedError
+       NotImplementedError()
        #raise ValueError
 
        # Don't complain if non-runnable code isn't run:
@@ -39,7 +41,7 @@ exclude_lines =
        if __name__ == .__main__.:
 
 skip_covered = True
-fail_under = 83
+fail_under = 85
 
 [html]
 directory = coverage_report
-- 
GitLab