pystencils_tests/test_datahandling.py → tests/test_datahandling.py

-import numpy as np
 import os
 from tempfile import TemporaryDirectory
+from pathlib import Path
+
+import numpy as np
+
 import pystencils as ps
-from pystencils import create_kernel, create_data_handling
+from pystencils import create_data_handling, create_kernel
+from pystencils.gpu.gpu_array_handler import GPUArrayHandler
+from pystencils.enums import Target
+
+try:
+    import pytest
+except ImportError:
+    import unittest.mock
+    pytest = unittest.mock.MagicMock()
+
+try:
+    import cupy.cuda.runtime
+    device_numbers = range(cupy.cuda.runtime.getDeviceCount())
+except ImportError:
+    device_numbers = []
+
+SCRIPT_FOLDER = Path(__file__).parent.absolute()
+INPUT_FOLDER = SCRIPT_FOLDER / "test_data"
 
 
 def basic_iteration(dh):
@@ -15,7 +35,7 @@ def basic_iteration(dh):
 
 
 def access_and_gather(dh, domain_size):
-    dh.add_array('f1', dtype=np.dtype(np.int32))
+    dh.add_array('f1', dtype=np.dtype(np.int8))
     dh.add_array_like('f2', 'f1')
     dh.add_array('v1', values_per_cell=3, dtype=np.int64, ghost_layers=2)
     dh.add_array_like('v2', 'v1')
@@ -26,7 +46,7 @@ def access_and_gather(dh, domain_size):
     # Check symbolic field properties
     assert dh.fields.f1.index_dimensions == 0
     assert dh.fields.f1.spatial_dimensions == len(domain_size)
-    assert dh.fields.f1.dtype.numpy_dtype == np.int32
+    assert dh.fields.f1.dtype.numpy_dtype == np.int8
 
     assert dh.fields.v1.index_dimensions == 1
     assert dh.fields.v1.spatial_dimensions == len(domain_size)
@@ -71,14 +91,10 @@ def access_and_gather(dh, domain_size):
 def synchronization(dh, test_gpu=False):
     field_name = 'comm_field_test'
     if test_gpu:
-        try:
-            from pycuda import driver
-            import pycuda.autoinit
-        except ImportError:
-            return
+        pytest.importorskip("cupy")
         field_name += 'Gpu'
 
-    dh.add_array(field_name, ghost_layers=1, dtype=np.int32, cpu=True, gpu=test_gpu)
+    dh.add_array(field_name, ghost_layers=1, dtype=np.int8, cpu=True, gpu=test_gpu)
 
     # initialize everything with 1
     for b in dh.iterate(ghost_layers=1):
@@ -88,8 +104,10 @@ def synchronization(dh, test_gpu=False):
 
     if test_gpu:
         dh.to_gpu(field_name)
+        dh.synchronization_function_gpu(field_name)()
+    else:
+        dh.synchronization_function_cpu(field_name)()
 
-    dh.synchronization_function(field_name, target='gpu' if test_gpu else 'cpu')()
 
     if test_gpu:
         dh.to_cpu(field_name)
@@ -98,17 +116,16 @@ def synchronization(dh, test_gpu=False):
         np.testing.assert_equal(42, b[field_name])
 
 
-def kernel_execution_jacobi(dh, test_gpu=False):
-    if test_gpu:
-        try:
-            from pycuda import driver
-            import pycuda.autoinit
-        except ImportError:
-            print("Skipping kernel_execution_jacobi GPU version, because pycuda not available")
-            return
+def kernel_execution_jacobi(dh, target):
 
+    test_gpu = target == Target.GPU
     dh.add_array('f', gpu=test_gpu)
     dh.add_array('tmp', gpu=test_gpu)
+
+    if test_gpu:
+        assert dh.is_on_gpu('f')
+        assert dh.is_on_gpu('tmp')
+
     stencil_2d = [(1, 0), (-1, 0), (0, 1), (0, -1)]
     stencil_3d = [(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1)]
     stencil = stencil_2d if dh.dim == 2 else stencil_3d
@@ -117,7 +134,7 @@ def kernel_execution_jacobi(dh, test_gpu=False):
     def jacobi():
         dh.fields.tmp.center @= sum(dh.fields.f.neighbors(stencil)) / len(stencil)
 
-    kernel = create_kernel(jacobi, target='gpu' if test_gpu else 'cpu').compile()
+    kernel = create_kernel(jacobi, config=ps.CreateKernelConfig(target=target)).compile()
     for b in dh.iterate(ghost_layers=1):
         b['f'].fill(42)
     dh.run_kernel(kernel)
@@ -126,6 +143,7 @@ def kernel_execution_jacobi(dh, test_gpu=False):
 
 
 def vtk_output(dh):
+    pytest.importorskip('pyevtk')
     dh.add_array('scalar_field')
     dh.add_array('vector_field', values_per_cell=dh.dim)
     dh.add_array('multiple_scalar_field', values_per_cell=9)
@@ -194,18 +212,200 @@ def test_access_and_gather():
 def test_kernel():
     for domain_shape in [(4, 5), (3, 4, 5)]:
         dh = create_data_handling(domain_size=domain_shape, periodicity=True)
-        kernel_execution_jacobi(dh, test_gpu=True)
+        assert all(dh.periodicity)
+        kernel_execution_jacobi(dh, Target.CPU)
         reduction(dh)
 
         try:
-            import pycuda
+            import cupy
             dh = create_data_handling(domain_size=domain_shape, periodicity=True)
-            kernel_execution_jacobi(dh, test_gpu=False)
+            kernel_execution_jacobi(dh, Target.GPU)
         except ImportError:
             pass
 
 
+@pytest.mark.parametrize('target', (Target.CPU, Target.GPU))
+def test_kernel_param(target):
+    for domain_shape in [(4, 5), (3, 4, 5)]:
+        if target == Target.GPU:
+            pytest.importorskip('cupy')
+
+        dh = create_data_handling(domain_size=domain_shape, periodicity=True, default_target=target)
+        kernel_execution_jacobi(dh, target)
+        reduction(dh)
+
+
 def test_vtk_output():
+    pytest.importorskip('pyevtk')
     for domain_shape in [(4, 5), (3, 4, 5)]:
         dh = create_data_handling(domain_size=domain_shape, periodicity=True)
         vtk_output(dh)
+
+
+def test_add_arrays():
+    domain_shape = (3, 4, 5)
+    field_description = 'x, y(9)'
+
+    dh = create_data_handling(domain_size=domain_shape, default_ghost_layers=0, default_layout='numpy')
+    x_, y_ = dh.add_arrays(field_description)
+
+    x, y = ps.fields(field_description + ': [3,4,5]')
+
+    assert x_ == x
+    assert y_ == y
+    assert x == dh.fields['x']
+    assert y == dh.fields['y']
+
+
+@pytest.mark.parametrize('shape', [(17, 12), (7, 11, 18)])
+@pytest.mark.parametrize('layout', ['zyxf', 'fzyx'])
+def test_add_arrays_with_layout(shape, layout):
+    pytest.importorskip('cupy')
+
+    dh = create_data_handling(domain_size=shape, default_layout=layout, default_target=ps.Target.GPU)
+    f1 = dh.add_array("f1", values_per_cell=19)
+    dh.fill(f1.name, 1.0)
+
+    assert dh.cpu_arrays[f1.name].shape == dh.gpu_arrays[f1.name].shape
+    assert dh.cpu_arrays[f1.name].strides == dh.gpu_arrays[f1.name].strides
+    assert dh.cpu_arrays[f1.name].dtype == dh.gpu_arrays[f1.name].dtype
+
+
+def test_get_kwarg():
+    domain_shape = (10, 10)
+    field_description = 'src, dst'
+
+    dh = create_data_handling(domain_size=domain_shape, default_ghost_layers=1)
+    src, dst = dh.add_arrays(field_description)
+    dh.fill("src", 1.0, ghost_layers=True)
+    dh.fill("dst", 0.0, ghost_layers=True)
+
+    with pytest.raises(ValueError):
+        dh.add_array('src', values_per_cell=1)
+
+    ur = ps.Assignment(src.center, dst.center)
+    kernel = ps.create_kernel(ur).compile()
+
+    kw = dh.get_kernel_kwargs(kernel)
+    assert np.all(kw[0]['src'] == dh.cpu_arrays['src'])
+    assert np.all(kw[0]['dst'] == dh.cpu_arrays['dst'])
+
+
+def test_add_custom_data():
+    pytest.importorskip('cupy')
+    import cupy as cp
+
+    def cpu_data_create_func():
+        return np.ones((2, 2), dtype=np.float64)
+
+    def gpu_data_create_func():
+        return cp.zeros((2, 2), dtype=np.float64)
+
+    def cpu_to_gpu_transfer_func(gpuarr, cpuarray):
+        gpuarr.set(cpuarray)
+
+    def gpu_to_cpu_transfer_func(gpuarr, cpuarray):
+        cpuarray[:] = gpuarr.get()
+
+    dh = create_data_handling(domain_size=(10, 10))
+    dh.add_custom_data('custom_data',
+                       cpu_data_create_func,
+                       gpu_data_create_func,
+                       cpu_to_gpu_transfer_func,
+                       gpu_to_cpu_transfer_func)
+
+    assert np.all(dh.custom_data_cpu['custom_data'] == 1)
+    assert np.all(dh.custom_data_gpu['custom_data'].get() == 0)
+
+    dh.to_cpu(name='custom_data')
+    dh.to_gpu(name='custom_data')
+
+    assert 'custom_data' in dh.custom_data_names
+
+
+def test_log():
+    dh = create_data_handling(domain_size=(10, 10))
+    dh.log_on_root()
+    assert dh.is_root
+    assert dh.world_rank == 0
+
+
+def test_save_data():
+    domain_shape = (2, 2)
+
+    dh = create_data_handling(domain_size=domain_shape, default_ghost_layers=1)
+    dh.add_array("src", values_per_cell=9)
+    dh.fill("src", 1.0, ghost_layers=True)
+    dh.add_array("dst", values_per_cell=9)
+    dh.fill("dst", 1.0, ghost_layers=True)
+
+    dh.save_all(str(INPUT_FOLDER) + '/datahandling_save_test')
+
+
+def test_load_data():
+    domain_shape = (2, 2)
+
+    dh = create_data_handling(domain_size=domain_shape, default_ghost_layers=1)
+    dh.add_array("src", values_per_cell=9)
+    dh.fill("src", 0.0, ghost_layers=True)
+    dh.add_array("dst", values_per_cell=9)
+    dh.fill("dst", 0.0, ghost_layers=True)
+
+    dh.load_all(str(INPUT_FOLDER) + '/datahandling_load_test')
+    assert np.all(dh.cpu_arrays['src']) == 1
+    assert np.all(dh.cpu_arrays['dst']) == 1
+
+    domain_shape = (3, 3)
+
+    dh = create_data_handling(domain_size=domain_shape, default_ghost_layers=1)
+    dh.add_array("src", values_per_cell=9)
+    dh.fill("src", 0.0, ghost_layers=True)
+    dh.add_array("dst", values_per_cell=9)
+    dh.fill("dst", 0.0, ghost_layers=True)
+    dh.add_array("dst2", values_per_cell=9)
+    dh.fill("dst2", 0.0, ghost_layers=True)
+
+    dh.load_all(str(INPUT_FOLDER) + '/datahandling_load_test')
+    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("device_number", device_numbers)
+def test_array_handler(device_number):
+    size = (2, 2)
+    pytest.importorskip('cupy')
+    array_handler = GPUArrayHandler(device_number)
+
+    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, order="F")
+    assert empty.strides == (8, 16)
+
+    random_array = array_handler.randn(size)
+
+    cpu_array = np.empty((20, 40), dtype=np.float64)
+    gpu_array = array_handler.to_gpu(cpu_array)
+
+    assert cpu_array.base is None
+    assert gpu_array.base is None
+    assert gpu_array.strides == cpu_array.strides
+
+    cpu_array2 = np.empty((20, 40), dtype=np.float64)
+    cpu_array2 = cpu_array2.swapaxes(0, 1)
+    gpu_array2 = array_handler.to_gpu(cpu_array2)
+
+    assert cpu_array2.base is not None
+    assert gpu_array2.base is not None
+    assert gpu_array2.strides == cpu_array2.strides
+

tests/test_datahandling_parallel.py

+import numpy as np
+import waLBerla as wlb
+
+import pystencils
+from pystencils import make_slice
+
+from pathlib import Path
+
+from pystencils.boundaries import BoundaryHandling, Neumann
+from pystencils.slicing import slice_from_direction
+
+from pystencils.datahandling.parallel_datahandling import ParallelDataHandling
+from pystencils.datahandling import create_data_handling
+from tests.test_datahandling import (
+    access_and_gather, kernel_execution_jacobi, reduction, synchronization, vtk_output)
+
+SCRIPT_FOLDER = Path(__file__).parent.absolute()
+INPUT_FOLDER = SCRIPT_FOLDER / "test_data"
+
+try:
+    import pytest
+except ImportError:
+    import unittest.mock
+    pytest = unittest.mock.MagicMock()
+
+
+def test_access_and_gather():
+    block_size = (4, 7, 1)
+    num_blocks = (3, 2, 1)
+    cells = tuple(a * b for a, b in zip(block_size, num_blocks))
+    blocks = wlb.createUniformBlockGrid(blocks=num_blocks, cellsPerBlock=block_size, oneBlockPerProcess=False,
+                                        periodic=(1, 1, 1))
+    dh = ParallelDataHandling(blocks, default_ghost_layers=2)
+    access_and_gather(dh, cells)
+    synchronization(dh, test_gpu=False)
+    if hasattr(wlb, 'gpu'):
+        synchronization(dh, test_gpu=True)
+
+
+def test_gpu():
+    pytest.importorskip('waLBerla.gpu')
+
+    block_size = (4, 7, 1)
+    num_blocks = (3, 2, 1)
+    blocks = wlb.createUniformBlockGrid(blocks=num_blocks, cellsPerBlock=block_size, oneBlockPerProcess=False)
+    dh = ParallelDataHandling(blocks, default_ghost_layers=2)
+    dh.add_array('v', values_per_cell=3, dtype=np.int64, ghost_layers=2, gpu=True)
+
+    for b in dh.iterate():
+        b['v'].fill(42)
+    dh.all_to_gpu()
+    for b in dh.iterate():
+        b['v'].fill(0)
+    dh.to_cpu('v')
+    for b in dh.iterate():
+        np.testing.assert_equal(b['v'], 42)
+
+
+@pytest.mark.parametrize('target', (pystencils.Target.CPU, pystencils.Target.GPU))
+def test_kernel(target):
+    if target == pystencils.Target.GPU:
+        pytest.importorskip('waLBerla.gpu')
+
+    # 3D
+    blocks = wlb.createUniformBlockGrid(blocks=(3, 2, 4), cellsPerBlock=(3, 2, 5), oneBlockPerProcess=False)
+    dh = ParallelDataHandling(blocks, default_target=target)
+    kernel_execution_jacobi(dh, target)
+    reduction(dh)
+
+    # 2D
+    blocks = wlb.createUniformBlockGrid(blocks=(3, 2, 1), cellsPerBlock=(3, 2, 1), oneBlockPerProcess=False)
+    dh = ParallelDataHandling(blocks, dim=2, default_target=target)
+    kernel_execution_jacobi(dh, target)
+    reduction(dh)
+
+
+def test_vtk_output():
+    blocks = wlb.createUniformBlockGrid(blocks=(3, 2, 4), cellsPerBlock=(3, 2, 5), oneBlockPerProcess=False)
+    dh = ParallelDataHandling(blocks)
+    vtk_output(dh)
+
+
+def test_block_iteration():
+    block_size = (16, 16, 16)
+    num_blocks = (2, 2, 2)
+    blocks = wlb.createUniformBlockGrid(blocks=num_blocks, cellsPerBlock=block_size, oneBlockPerProcess=False)
+    dh = ParallelDataHandling(blocks, default_ghost_layers=2)
+    dh.add_array('v', values_per_cell=1, dtype=np.int64, ghost_layers=2)
+
+    for b in dh.iterate():
+        b['v'].fill(1)
+
+    s = 0
+    for b in dh.iterate():
+        s += np.sum(b['v'])
+
+    assert s == 40*40*40
+
+    sl = make_slice[0:18, 0:18, 0:18]
+    for b in dh.iterate(slice_obj=sl):
+        b['v'].fill(0)
+
+    s = 0
+    for b in dh.iterate():
+        s += np.sum(b['v'])
+
+    assert s == 40*40*40 - 20*20*20
+
+
+def test_getter_setter():
+    pytest.importorskip('waLBerla.gpu')
+
+    block_size = (2, 2, 2)
+    num_blocks = (2, 2, 2)
+    blocks = wlb.createUniformBlockGrid(blocks=num_blocks, cellsPerBlock=block_size, oneBlockPerProcess=False)
+    dh = ParallelDataHandling(blocks, default_ghost_layers=2, default_target=pystencils.Target.GPU)
+    dh.add_array('v', values_per_cell=1, dtype=np.int64, ghost_layers=2, gpu=True)
+
+    assert dh.shape == (4, 4, 4)
+    assert dh.periodicity == (False, False, False)
+    assert dh.values_per_cell('v') == 1
+    assert dh.has_data('v') is True
+    assert 'v' in dh.array_names
+    dh.log_on_root()
+    assert dh.is_root is True
+    assert dh.world_rank == 0
+
+    dh.to_gpu('v')
+    assert dh.is_on_gpu('v') is True
+    dh.all_to_cpu()
+
+
+def test_parallel_datahandling_boundary_conditions():
+    pytest.importorskip('waLBerla.gpu')
+
+    dh = create_data_handling(domain_size=(7, 7), periodicity=True, parallel=True,
+                              default_target=pystencils.Target.GPU)
+
+    src = dh.add_array('src', values_per_cell=1)
+    dh.fill(src.name, 0.0, ghost_layers=True)
+    dh.fill(src.name, 1.0, ghost_layers=False)
+
+    src2 = dh.add_array('src2', values_per_cell=1)
+
+    src_cpu = dh.add_array('src_cpu', values_per_cell=1, gpu=False)
+    dh.fill(src_cpu.name, 0.0, ghost_layers=True)
+    dh.fill(src_cpu.name, 1.0, ghost_layers=False)
+
+    boundary_stencil = [(1, 0), (-1, 0), (0, 1), (0, -1)]
+    boundary_handling_cpu = BoundaryHandling(dh, src_cpu.name, boundary_stencil,
+                                             name="boundary_handling_cpu", target=pystencils.Target.CPU)
+
+    boundary_handling = BoundaryHandling(dh, src.name, boundary_stencil,
+                                         name="boundary_handling_gpu", target=pystencils.Target.GPU)
+
+    neumann = Neumann()
+    for d in ('N', 'S', 'W', 'E'):
+        boundary_handling.set_boundary(neumann, slice_from_direction(d, dim=2))
+        boundary_handling_cpu.set_boundary(neumann, slice_from_direction(d, dim=2))
+
+    boundary_handling.prepare()
+    boundary_handling_cpu.prepare()
+
+    boundary_handling_cpu()
+
+    dh.all_to_gpu()
+    boundary_handling()
+    dh.all_to_cpu()
+    for block in dh.iterate():
+        np.testing.assert_almost_equal(block[src_cpu.name], block[src.name])
+
+    assert dh.custom_data_names == ('boundary_handling_cpuIndexArrays', 'boundary_handling_gpuIndexArrays')
+    dh.swap(src.name, src2.name, gpu=True)
+
+
+def test_save_data():
+    domain_shape = (2, 2)
+
+    dh = create_data_handling(domain_size=domain_shape, default_ghost_layers=1, parallel=True)
+    dh.add_array("src", values_per_cell=9)
+    dh.fill("src", 1.0, ghost_layers=True)
+    dh.add_array("dst", values_per_cell=9)
+    dh.fill("dst", 1.0, ghost_layers=True)
+
+    dh.save_all(str(INPUT_FOLDER) + '/datahandling_parallel_save_test')
+
+
+def test_load_data():
+    domain_shape = (2, 2)
+
+    dh = create_data_handling(domain_size=domain_shape, default_ghost_layers=1, parallel=True)
+    dh.add_array("src", values_per_cell=9)
+    dh.fill("src", 0.0, ghost_layers=True)
+    dh.add_array("dst", values_per_cell=9)
+    dh.fill("dst", 0.0, ghost_layers=True)
+
+    dh.load_all(str(INPUT_FOLDER) + '/datahandling_parallel_load_test')
+    assert np.all(dh.gather_array('src')) == 1
+    assert np.all(dh.gather_array('src')) == 1

pystencils_tests/test_derivative.py → tests/test_derivative.py

 import sympy as sp
-import pystencils as ps
+from sympy.abc import a, b, t, x, y, z
 from sympy.printing.latex import LatexPrinter
+
+import pystencils as ps
 from pystencils.fd import *
-from sympy.abc import a, b, t, x, y, z
 
 
 def test_derivative_basic():

tests/test_dot_printer.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+from pystencils.astnodes import Block, Conditional, SympyAssignment
+```
+
+%% Cell type:code id: tags:
+
+``` python
+src, dst = ps.fields("src, dst: double[2D]", layout='c')
+
+true_block = Block([SympyAssignment(dst[0, 0], src[-1, 0])])
+false_block = Block([SympyAssignment(dst[0, 0], src[1, 0])])
+ur = [true_block, Conditional(dst.center() > 0.0, true_block, false_block)]
+
+ast = ps.create_kernel(ur)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+ps.show_code(ast)
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+```

tests/test_dot_printer.py

+import pystencils as ps
+
+from pystencils.astnodes import Block, Conditional, SympyAssignment
+
+
+def test_dot_print():
+    src, dst = ps.fields("src, dst: double[2D]", layout='c')
+
+    true_block = Block([SympyAssignment(dst[0, 0], src[-1, 0])])
+    false_block = Block([SympyAssignment(dst[0, 0], src[1, 0])])
+    ur = [true_block, Conditional(dst.center() > 0.0, true_block, false_block)]
+
+    ast = ps.create_kernel(ur)

pystencils_tests/test_dtype_check.py → tests/test_dtype_check.py

-import pystencils
 import numpy as np
 import pytest
 
+import pystencils
+
 
 def test_dtype_check_wrong_type():
     array = np.ones((10, 20)).astype(np.float32)
@@ -15,7 +16,7 @@ def test_dtype_check_wrong_type():
 
     with pytest.raises(ValueError) as e:
         kernel(x=array, y=output)
-    assert 'Wrong data type' in str(e)
+    assert 'Wrong data type' in str(e.value)
 
 
 def test_dtype_check_correct_type():

pystencils_tests/test_fast_approximation.py → tests/test_fast_approximation.py

+import pytest
 import sympy as sp
+
 import pystencils as ps
-from pystencils.fast_approximation import insert_fast_divisions, insert_fast_sqrts, fast_sqrt, fast_inv_sqrt, \
-    fast_division
+from pystencils.fast_approximation import (
+    fast_division, fast_inv_sqrt, fast_sqrt, insert_fast_divisions, insert_fast_sqrts)
 
 
 def test_fast_sqrt():
+    pytest.importorskip('cupy')
     f, g = ps.fields("f, g: double[2D]")
     expr = sp.sqrt(f[0, 0] + f[1, 0])
 
     assert len(insert_fast_sqrts(expr).atoms(fast_sqrt)) == 1
     assert len(insert_fast_sqrts([expr])[0].atoms(fast_sqrt)) == 1
-    ast = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_sqrts(expr)), target='gpu')
-    code_str = str(ps.show_code(ast))
+    ast_gpu = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_sqrts(expr)), target=ps.Target.GPU)
+    ast_gpu.compile()
+    code_str = ps.get_code_str(ast_gpu)
     assert '__fsqrt_rn' in code_str
 
     expr = ps.Assignment(sp.Symbol("tmp"), 3 / sp.sqrt(f[0, 0] + f[1, 0]))
@@ -19,12 +23,14 @@ def test_fast_sqrt():
 
     ac = ps.AssignmentCollection([expr], [])
     assert len(insert_fast_sqrts(ac).main_assignments[0].atoms(fast_inv_sqrt)) == 1
-    ast = ps.create_kernel(insert_fast_sqrts(ac), target='gpu')
-    code_str = str(ps.show_code(ast))
+    ast_gpu = ps.create_kernel(insert_fast_sqrts(ac), target=ps.Target.GPU)
+    ast_gpu.compile()
+    code_str = ps.get_code_str(ast_gpu)
     assert '__frsqrt_rn' in code_str
 
 
 def test_fast_divisions():
+    pytest.importorskip('cupy')
     f, g = ps.fields("f, g: double[2D]")
     expr = f[0, 0] / f[1, 0]
     assert len(insert_fast_divisions(expr).atoms(fast_division)) == 1
@@ -32,6 +38,7 @@ def test_fast_divisions():
     expr = 1 / f[0, 0] * 2 / f[0, 1]
     assert len(insert_fast_divisions(expr).atoms(fast_division)) == 1
 
-    ast = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_divisions(expr)), target='gpu')
-    code_str = str(ps.show_code(ast))
+    ast = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_divisions(expr)), target=ps.Target.GPU)
+    ast.compile()
+    code_str = ps.get_code_str(ast)
     assert '__fdividef' in code_str

tests/test_fd_derivation.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+from pystencils.fd.derivation import *
+```
+
+%% Cell type:markdown id: tags:
+
+# 2D standard stencils
+
+%% Cell type:code id: tags:
+
+``` python
+stencil = [(-1, 0), (1, 0), (0, -1), (0, 1), (0, 0)]
+standard_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)
+f = ps.fields("f: [2D]")
+standard_2d_00_res = standard_2d_00.get_stencil()
+res = standard_2d_00_res.apply(f.center)
+expected = f[-1, 0] - 2 * f[0, 0] + f[1, 0]
+assert res == expected
+```
+
+%% Cell type:code id: tags:
+
+``` python
+assert standard_2d_00_res.accuracy == 2
+assert not standard_2d_00_res.is_isotropic
+standard_2d_00_res
+```
+
+%% Output
+
+    Finite difference stencil of accuracy 2, isotropic error: False
+
+%% Cell type:code id: tags:
+
+``` python
+standard_2d_00.get_stencil().as_array()
+```
+
+%% Output
+
+    $\displaystyle \left[\begin{matrix}0 & 0 & 0\\1 & -2 & 1\\0 & 0 & 0\end{matrix}\right]$
+    [[0, 0, 0], [1, -2, 1], [0, 0, 0]]
+
+%% Cell type:markdown id: tags:
+
+# 2D isotropic stencils
+
+## second x-derivative
+
+%% Cell type:code id: tags:
+
+``` python
+stencil = [(i, j) for i in (-1, 0, 1) for j in (-1, 0, 1)]
+isotropic_2d_00 = FiniteDifferenceStencilDerivation((0,0), stencil)
+isotropic_2d_00_res = isotropic_2d_00.get_stencil(isotropic=True)
+assert isotropic_2d_00_res.is_isotropic
+assert isotropic_2d_00_res.accuracy == 2
+isotropic_2d_00_res
+```
+
+%% Output
+
+    Finite difference stencil of accuracy 2, isotropic error: True
+
+%% Cell type:code id: tags:
+
+``` python
+isotropic_2d_00_res.as_array()
+```
+
+%% Output
+
+    $\displaystyle \left[\begin{matrix}\frac{1}{12} & - \frac{1}{6} & \frac{1}{12}\\\frac{5}{6} & - \frac{5}{3} & \frac{5}{6}\\\frac{1}{12} & - \frac{1}{6} & \frac{1}{12}\end{matrix}\right]$
+    [[1/12, -1/6, 1/12], [5/6, -5/3, 5/6], [1/12, -1/6, 1/12]]
+
+%% Cell type:code id: tags:
+
+``` python
+plt.figure(figsize=(2,2))
+isotropic_2d_00_res.visualize()
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+expected_result = sp.Array([[1, -2, 1], [10, -20, 10], [1, -2, 1]]) / 12
+assert expected_result == isotropic_2d_00_res.as_array()
+```
+
+%% Cell type:markdown id: tags:
+
+## Isotropic laplacian
+
+%% Cell type:code id: tags:
+
+``` python
+isotropic_2d_11 = FiniteDifferenceStencilDerivation((1,1), stencil)
+isotropic_2d_11_res = isotropic_2d_11.get_stencil(isotropic=True)
+iso_laplacian = isotropic_2d_00_res.as_array() + isotropic_2d_11_res.as_array()
+iso_laplacian
+```
+
+%% Output
+
+    $\displaystyle \left[\begin{matrix}\frac{1}{6} & \frac{2}{3} & \frac{1}{6}\\\frac{2}{3} & - \frac{10}{3} & \frac{2}{3}\\\frac{1}{6} & \frac{2}{3} & \frac{1}{6}\end{matrix}\right]$
+    [[1/6, 2/3, 1/6], [2/3, -10/3, 2/3], [1/6, 2/3, 1/6]]
+
+%% Cell type:code id: tags:
+
+``` python
+expected_result = sp.Array([[1, 4, 1], [4, -20, 4], [1, 4, 1]]) / 6
+assert iso_laplacian == expected_result
+```
+
+%% Cell type:markdown id: tags:
+
+# stencils for staggered fields
+
+%% Cell type:code id: tags:
+
+``` python
+half = sp.Rational(1, 2)
+
+fd_points_ex = (
+    (half, 0),
+    (-half, 0),
+    (half, 1),
+    (half, -1),
+    (-half, 1),
+    (-half, -1)
+)
+assert set(fd_points_ex) == set(FiniteDifferenceStaggeredStencilDerivation("E", 2).stencil)
+
+fd_points_ey = (
+    (0, half),
+    (0, -half),
+    (-1,-half),
+    (-1, half),
+    (1, -half),
+    (1, half)
+)
+assert set(fd_points_ey) == set(FiniteDifferenceStaggeredStencilDerivation("N",2).stencil)
+
+fd_points_c = (
+    (half, half),
+    (-half, -half),
+    (half, -half),
+    (-half, half)
+)
+assert set(fd_points_c) ==  set(FiniteDifferenceStaggeredStencilDerivation("NE",2).stencil)
+
+assert len(FiniteDifferenceStaggeredStencilDerivation("E",3).points) == 10
+assert len(FiniteDifferenceStaggeredStencilDerivation("NE",3).points) == 12
+assert len(FiniteDifferenceStaggeredStencilDerivation("TNE",3).points) == 8
+```
+
+%% Cell type:code id: tags:
+
+``` python
+c = ps.fields("c: [2D]")
+c3 = ps.fields("c3: [3D]")
+
+assert FiniteDifferenceStaggeredStencilDerivation("E", 2, (0,)).apply(c.center) == c[1, 0] - c[0, 0]
+assert FiniteDifferenceStaggeredStencilDerivation("W", 2, (0,)).apply(c.center) == c[0, 0] - c[-1, 0]
+assert FiniteDifferenceStaggeredStencilDerivation("N", 2, (1,)).apply(c.center) == c[0, 1] - c[0, 0]
+assert FiniteDifferenceStaggeredStencilDerivation("S", 2, (1,)).apply(c.center) == c[0, 0] - c[0, -1]
+
+assert FiniteDifferenceStaggeredStencilDerivation("E", 3, (0,)).apply(c3.center) == c3[1, 0, 0] - c3[0, 0, 0]
+assert FiniteDifferenceStaggeredStencilDerivation("W", 3, (0,)).apply(c3.center) == c3[0, 0, 0] - c3[-1, 0, 0]
+assert FiniteDifferenceStaggeredStencilDerivation("N", 3, (1,)).apply(c3.center) == c3[0, 1, 0] - c3[0, 0, 0]
+assert FiniteDifferenceStaggeredStencilDerivation("S", 3, (1,)).apply(c3.center) == c3[0, 0, 0] - c3[0, -1, 0]
+assert FiniteDifferenceStaggeredStencilDerivation("T", 3, (2,)).apply(c3.center) == c3[0, 0, 1] - c3[0, 0, 0]
+assert FiniteDifferenceStaggeredStencilDerivation("B", 3, (2,)).apply(c3.center) == c3[0, 0, 0] - c3[0, 0, -1]
+
+assert FiniteDifferenceStaggeredStencilDerivation("S", 2, (0,)).apply(c.center) == \
+       (c[1, 0] + c[1, -1] - c[-1, 0] - c[-1, -1])/4
+
+assert FiniteDifferenceStaggeredStencilDerivation("NE", 2, (0,)).apply(c.center) + \
+       FiniteDifferenceStaggeredStencilDerivation("NE", 2, (1,)).apply(c.center) == c[1, 1] - c[0, 0]
+assert FiniteDifferenceStaggeredStencilDerivation("NE", 3, (0,)).apply(c3.center) + \
+       FiniteDifferenceStaggeredStencilDerivation("NE", 3, (1,)).apply(c3.center) == c3[1, 1, 0] - c3[0, 0, 0]
+```
+
+%% Cell type:code id: tags:
+
+``` python
+d = FiniteDifferenceStaggeredStencilDerivation("NE", 2, (0, 1))
+assert d.apply(c.center) == c[0,0] + c[1,1] - c[1,0] - c[0,1]
+d.visualize()
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+v3 = ps.fields("v(3): [3D]")
+for i in range(*v3.index_shape):
+    assert FiniteDifferenceStaggeredStencilDerivation("E", 3, (0,)).apply(v3.center_vector[i]) == \
+           v3[1,0,0](i) - v3[0,0,0](i)
+```

tests/test_fd_derivation_via_rotation.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from pystencils.session import *
+from pystencils.fd.derivation import *
+```
+
+%% Cell type:markdown id: tags:
+
+# 2D isotropic stencils
+
+%% Cell type:code id: tags:
+
+``` python
+stencil_2D = ((0, 0),
+              (0, 1), (0, -1), (-1, 0), (1, 0),
+              (-1, 1), (1, 1), (-1, -1), (1, -1))
+f = ps.fields("f: [2D]")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+Isotropic_Gradient = list()
+deriv_x = FiniteDifferenceStencilDerivation((0, ), stencil_2D)
+
+deriv_x.assume_symmetric(0, anti_symmetric=True)
+deriv_x.assume_symmetric(1, anti_symmetric=False)
+
+deriv_x.set_weight((0, 0), sp.Integer(0))
+deriv_x.set_weight((1, 0), sp.Rational(1, 3))
+deriv_x.set_weight((1, 1), sp.Rational(1, 12))
+
+res_x = deriv_x.get_stencil(isotropic=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+Isotropic_Gradient = list()
+deriv_y = FiniteDifferenceStencilDerivation((1, ), stencil_2D)
+
+deriv_y.assume_symmetric(0, anti_symmetric=False)
+deriv_y.assume_symmetric(1, anti_symmetric=True)
+
+deriv_y.set_weight((0, 0), sp.Integer(0))
+deriv_y.set_weight((0, 1), sp.Rational(1, 3))
+deriv_y.set_weight((1, 1), sp.Rational(1, 12))
+
+res_y = deriv_y.get_stencil(isotropic=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+assert res_x.apply(f.center) - res_y.rotate_weights_and_apply(f.center, (1, 0)) == 0
+assert res_y.apply(f.center) - res_x.rotate_weights_and_apply(f.center, (0, 1)) == 0
+```
+
+%% Cell type:markdown id: tags:
+
+# 3D isotropic stencils
+
+%% Cell type:code id: tags:
+
+``` python
+stencil_3D = ((0, 0, 0),
+             (0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
+             (-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
+             (0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
+             (0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1),
+             (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1),
+             (-1, -1, -1))
+f = ps.fields("f: [3D]")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+deriv_x = FiniteDifferenceStencilDerivation((0, ), stencil_3D)
+
+deriv_x.assume_symmetric(0, anti_symmetric=True)
+deriv_x.assume_symmetric(1, anti_symmetric=False)
+deriv_x.assume_symmetric(2, anti_symmetric=False)
+
+deriv_x.set_weight((0, 0, 0), sp.Integer(0))
+deriv_x.set_weight((1, 1, 1), sp.Rational(1, 3360))
+
+res_x = deriv_x.get_stencil(isotropic=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+deriv_y = FiniteDifferenceStencilDerivation((1, ), stencil_3D)
+
+deriv_y.assume_symmetric(0, anti_symmetric=False)
+deriv_y.assume_symmetric(1, anti_symmetric=True)
+deriv_y.assume_symmetric(2, anti_symmetric=False)
+
+deriv_y.set_weight((0, 0, 0), sp.Integer(0))
+deriv_y.set_weight((1, 1, 1), sp.Rational(1, 3360))
+
+res_y = deriv_y.get_stencil(isotropic=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+deriv_z = FiniteDifferenceStencilDerivation((2, ), stencil_3D)
+
+deriv_z.assume_symmetric(0, anti_symmetric=False)
+deriv_z.assume_symmetric(1, anti_symmetric=False)
+deriv_z.assume_symmetric(2, anti_symmetric=True)
+
+deriv_z.set_weight((0, 0, 0), sp.Integer(0))
+deriv_z.set_weight((1, 1, 1), sp.Rational(1, 3360))
+
+res_z = deriv_z.get_stencil(isotropic=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+assert res_x.apply(f.center) - res_y.rotate_weights_and_apply(f.center, (1, 0)) == 0
+assert res_x.apply(f.center) - res_z.rotate_weights_and_apply(f.center, (2, 1)) == 0
+
+assert res_y.apply(f.center) - res_x.rotate_weights_and_apply(f.center, (0, 1)) == 0
+assert res_y.apply(f.center) - res_z.rotate_weights_and_apply(f.center, (0, 2)) == 0
+
+assert res_z.apply(f.center) - res_x.rotate_weights_and_apply(f.center, (1, 2)) == 0
+assert res_z.apply(f.center) - res_y.rotate_weights_and_apply(f.center, (2, 0)) == 0
+```

tests/test_fd_derivative.py

+import sympy as sp
+from pystencils import fields
+from pystencils.fd import Diff, diff, collect_diffs
+from pystencils.fd.derivative import replace_generic_laplacian
+
+
+def test_fs():
+    f = sp.Symbol("f", commutative=False)
+
+    a = Diff(Diff(Diff(f, 1), 0), 0)
+    assert a.is_commutative is False
+    print(str(a))
+
+    assert diff(f) == f
+
+    x, y = sp.symbols("x, y")
+    collected_terms = collect_diffs(diff(x, 0, 0))
+    assert collected_terms == Diff(Diff(x, 0, -1), 0, -1)
+
+    src = fields("src : double[2D]")
+    expr = sp.Add(Diff(Diff(src[0, 0])), 10)
+    expected = Diff(Diff(src[0, 0], 0, -1), 0, -1) + Diff(Diff(src[0, 0], 1, -1), 1, -1) + 10
+    result = replace_generic_laplacian(expr, 3)
+    assert result == expected
\ No newline at end of file

tests/test_field.py

+import numpy as np
+import pytest
+import sympy as sp
+
+import pystencils as ps
+from pystencils import TypedSymbol
+from pystencils.typing import create_type
+from pystencils.field import Field, FieldType, layout_string_to_tuple, spatial_layout_string_to_tuple
+
+
+def test_field_basic():
+    f = Field.create_generic("f", spatial_dimensions=2)
+    assert FieldType.is_generic(f)
+    assert f["E"] == f[1, 0]
+    assert f["N"] == f[0, 1]
+    assert "_" in f.center._latex("dummy")
+
+    assert (
+        f.index_to_physical(index_coordinates=sp.Matrix([0, 0]), staggered=False)[0]
+        == 0
+    )
+    assert (
+        f.index_to_physical(index_coordinates=sp.Matrix([0, 0]), staggered=False)[1]
+        == 0
+    )
+
+    assert (
+        f.physical_to_index(physical_coordinates=sp.Matrix([0, 0]), staggered=False)[0]
+        == 0
+    )
+    assert (
+        f.physical_to_index(physical_coordinates=sp.Matrix([0, 0]), staggered=False)[1]
+        == 0
+    )
+
+    f1 = f.new_field_with_different_name("f1")
+    assert f1.ndim == f.ndim
+    assert f1.values_per_cell() == f.values_per_cell()
+
+    fixed = ps.fields("f(5, 5) : double[20, 20]")
+    assert fixed.neighbor_vector((1, 1)).shape == (5, 5)
+
+    f = Field.create_fixed_size("f", (10, 10), strides=(80, 8), dtype=np.float64)
+    assert f.spatial_strides == (10, 1)
+    assert f.index_strides == ()
+    assert f.center_vector == sp.Matrix([f.center])
+
+    f1 = f.new_field_with_different_name("f1")
+    assert f1.ndim == f.ndim
+    assert f1.values_per_cell() == f.values_per_cell()
+
+    f = Field.create_fixed_size("f", (8, 8, 2, 2), index_dimensions=2)
+    assert f.center_vector == sp.Matrix([[f(0, 0), f(0, 1)], [f(1, 0), f(1, 1)]])
+    field_access = f[1, 1]
+    assert field_access.nr_of_coordinates == 2
+    assert field_access.offset_name == "NE"
+    neighbor = field_access.neighbor(coord_id=0, offset=-2)
+    assert neighbor.offsets == (-1, 1)
+    assert "_" in neighbor._latex("dummy")
+
+    f = Field.create_fixed_size("f", (8, 8, 2, 2, 2), index_dimensions=3)
+    assert f.center_vector == sp.Array(
+        [[[f(i, j, k) for k in range(2)] for j in range(2)] for i in range(2)]
+    )
+
+    f = Field.create_generic("f", spatial_dimensions=5, index_dimensions=2)
+    field_access = f[1, -1, 2, -3, 0](1, 0)
+    assert field_access.offsets == (1, -1, 2, -3, 0)
+    assert field_access.index == (1, 0)
+
+
+def test_error_handling():
+    struct_dtype = np.dtype([('a', np.int32), ('b', np.float64), ('c', np.uint32)])
+    Field.create_generic('f', spatial_dimensions=2, index_dimensions=0, dtype=struct_dtype)
+    with pytest.raises(ValueError) as e:
+        Field.create_generic(
+            "f", spatial_dimensions=2, index_dimensions=1, dtype=struct_dtype
+        )
+    assert "index dimension" in str(e.value)
+
+    arr = np.array([[[(1,) * 3, (2,) * 3, (3,) * 3]] * 2], dtype=struct_dtype)
+    Field.create_from_numpy_array("f", arr, index_dimensions=0)
+    with pytest.raises(ValueError) as e:
+        Field.create_from_numpy_array("f", arr, index_dimensions=1)
+    assert "Structured arrays" in str(e.value)
+
+    arr = np.zeros([3, 3, 3])
+    Field.create_from_numpy_array("f", arr, index_dimensions=2)
+    with pytest.raises(ValueError) as e:
+        Field.create_from_numpy_array("f", arr, index_dimensions=3)
+    assert "Too many" in str(e.value)
+
+    Field.create_fixed_size(
+        "f", (3, 2, 4), index_dimensions=0, dtype=struct_dtype, layout="reverse_numpy"
+    )
+    with pytest.raises(ValueError) as e:
+        Field.create_fixed_size(
+            "f",
+            (3, 2, 4),
+            index_dimensions=1,
+            dtype=struct_dtype,
+            layout="reverse_numpy",
+        )
+    assert "Structured arrays" in str(e.value)
+
+    f = Field.create_fixed_size("f", (10, 10))
+    with pytest.raises(ValueError) as e:
+        f[1]
+    assert "Wrong number of spatial indices" in str(e.value)
+
+    f = Field.create_generic("f", spatial_dimensions=2, index_shape=(3,))
+    with pytest.raises(ValueError) as e:
+        f(3)
+    assert "out of bounds" in str(e.value)
+
+    f = Field.create_fixed_size("f", (10, 10, 3, 4), index_dimensions=2)
+    with pytest.raises(ValueError) as e:
+        f(3, 0)
+    assert "out of bounds" in str(e.value)
+
+    with pytest.raises(ValueError) as e:
+        f(1, 0)(1, 0)
+    assert "Indexing an already indexed" in str(e.value)
+
+    with pytest.raises(ValueError) as e:
+        f(1)
+    assert "Wrong number of indices" in str(e.value)
+
+    with pytest.raises(ValueError) as e:
+        Field.create_generic("f", spatial_dimensions=2, layout="wrong")
+    assert "Unknown layout descriptor" in str(e.value)
+
+    assert layout_string_to_tuple("fzyx", dim=4) == (3, 2, 1, 0)
+    with pytest.raises(ValueError) as e:
+        layout_string_to_tuple("wrong", dim=4)
+    assert "Unknown layout descriptor" in str(e.value)
+
+
+def test_decorator_scoping():
+    dst = ps.fields("dst : double[2D]")
+
+    def f1():
+        a = sp.Symbol("a")
+
+        def f2():
+            b = sp.Symbol("b")
+
+            @ps.kernel
+            def decorated_func():
+                dst[0, 0] @= a + b
+
+            return decorated_func
+
+        return f2
+
+    assert f1()(), ps.Assignment(dst[0, 0], sp.Symbol("a") + sp.Symbol("b"))
+
+
+def test_string_creation():
+    x, y, z = ps.fields("  x(4),    y(3,5) z : double[  3,  47]")
+    assert x.index_shape == (4,)
+    assert y.index_shape == (3, 5)
+    assert z.spatial_shape == (3, 47)
+
+
+def test_itemsize():
+
+    x = ps.fields("x: float32[1d]")
+    y = ps.fields("y:  float64[2d]")
+    i = ps.fields("i:  int16[1d]")
+
+    assert x.itemsize == 4
+    assert y.itemsize == 8
+    assert i.itemsize == 2
+
+
+def test_spatial_memory_layout_descriptors():
+    assert (
+        spatial_layout_string_to_tuple("AoS", 3)
+        == spatial_layout_string_to_tuple("aos", 3)
+        == spatial_layout_string_to_tuple("ZYXF", 3)
+        == spatial_layout_string_to_tuple("zyxf", 3)
+        == (2, 1, 0)
+    )
+    assert (
+        spatial_layout_string_to_tuple("SoA", 3)
+        == spatial_layout_string_to_tuple("soa", 3)
+        == spatial_layout_string_to_tuple("FZYX", 3)
+        == spatial_layout_string_to_tuple("fzyx", 3)
+        == spatial_layout_string_to_tuple("f", 3)
+        == spatial_layout_string_to_tuple("F", 3)
+        == (2, 1, 0)
+    )
+    assert (
+        spatial_layout_string_to_tuple("c", 3)
+        == spatial_layout_string_to_tuple("C", 3)
+        == (0, 1, 2)
+    )
+
+    assert spatial_layout_string_to_tuple("C", 5) == (0, 1, 2, 3, 4)
+
+    with pytest.raises(ValueError):
+        spatial_layout_string_to_tuple("aos", -1)
+
+    with pytest.raises(ValueError):
+        spatial_layout_string_to_tuple("aos", 4)
+
+
+def test_memory_layout_descriptors():
+    assert (
+        layout_string_to_tuple("AoS", 4)
+        == layout_string_to_tuple("aos", 4)
+        == layout_string_to_tuple("ZYXF", 4)
+        == layout_string_to_tuple("zyxf", 4)
+        == (2, 1, 0, 3)
+    )
+    assert (
+        layout_string_to_tuple("SoA", 4)
+        == layout_string_to_tuple("soa", 4)
+        == layout_string_to_tuple("FZYX", 4)
+        == layout_string_to_tuple("fzyx", 4)
+        == layout_string_to_tuple("f", 4)
+        == layout_string_to_tuple("F", 4)
+        == (3, 2, 1, 0)
+    )
+    assert (
+        layout_string_to_tuple("c", 4)
+        == layout_string_to_tuple("C", 4)
+        == (0, 1, 2, 3)
+    )
+
+    assert layout_string_to_tuple("C", 5) == (0, 1, 2, 3, 4)
+
+    with pytest.raises(ValueError):
+        layout_string_to_tuple("aos", -1)
+
+    with pytest.raises(ValueError):
+        layout_string_to_tuple("aos", 5)
+
+
+def test_staggered():
+
+    # D2Q5
+    j1, j2, j3 = ps.fields(
+        "j1(2), j2(2,2), j3(2,2,2) : double[2D]", field_type=FieldType.STAGGERED
+    )
+
+    assert j1[0, 1](1) == j1.staggered_access((0, sp.Rational(1, 2)))
+    assert j1[0, 1](1) == j1.staggered_access(np.array((0, sp.Rational(1, 2))))
+    assert j1[1, 1](1) == j1.staggered_access((1, sp.Rational(1, 2)))
+    assert j1[0, 2](1) == j1.staggered_access((0, sp.Rational(3, 2)))
+    assert j1[0, 1](1) == j1.staggered_access("N")
+    assert j1[0, 0](1) == j1.staggered_access("S")
+    assert j1.staggered_vector_access("N") == sp.Matrix([j1.staggered_access("N")])
+    assert j1.staggered_stencil_name == "D2Q5"
+
+    assert j1.physical_coordinates[0] == TypedSymbol("ctr_0", create_type("int"), nonnegative=True)
+    assert j1.physical_coordinates[1] == TypedSymbol("ctr_1", create_type("int"), nonnegative=True)
+    assert j1.physical_coordinates_staggered[0] == TypedSymbol("ctr_0", create_type("int"), nonnegative=True) + 0.5
+    assert j1.physical_coordinates_staggered[1] == TypedSymbol("ctr_1", create_type("int"), nonnegative=True) + 0.5
+    assert j1.index_to_physical(index_coordinates=sp.Matrix([0, 0]), staggered=True)[0] == 0.5
+    assert j1.index_to_physical(index_coordinates=sp.Matrix([0, 0]), staggered=True)[1] == 0.5
+    assert j1.physical_to_index(physical_coordinates=sp.Matrix([0, 0]), staggered=True)[0] == -0.5
+    assert j1.physical_to_index(physical_coordinates=sp.Matrix([0, 0]), staggered=True)[1] == -0.5
+
+    assert j2[0, 1](1, 1) == j2.staggered_access((0, sp.Rational(1, 2)), 1)
+    assert j2[0, 1](1, 1) == j2.staggered_access("N", 1)
+    assert j2.staggered_vector_access("N") == sp.Matrix(
+        [j2.staggered_access("N", 0), j2.staggered_access("N", 1)]
+    )
+
+    assert j3[0, 1](1, 1, 1) == j3.staggered_access((0, sp.Rational(1, 2)), (1, 1))
+    assert j3[0, 1](1, 1, 1) == j3.staggered_access("N", (1, 1))
+    assert j3.staggered_vector_access("N") == sp.Matrix(
+        [[j3.staggered_access("N", (i, j)) for j in range(2)] for i in range(2)]
+    )
+
+    # D2Q9
+    k1, k2 = ps.fields("k1(4), k2(2) : double[2D]", field_type=FieldType.STAGGERED)
+
+    assert k1[1, 1](2) == k1.staggered_access("NE")
+    assert k1[0, 0](2) == k1.staggered_access("SW")
+    assert k1[0, 0](3) == k1.staggered_access("NW")
+
+    a = k1.staggered_access("NE")
+    assert a._staggered_offset(a.offsets, a.index[0]) == [
+        sp.Rational(1, 2),
+        sp.Rational(1, 2),
+    ]
+    a = k1.staggered_access("SW")
+    assert a._staggered_offset(a.offsets, a.index[0]) == [
+        sp.Rational(-1, 2),
+        sp.Rational(-1, 2),
+    ]
+    a = k1.staggered_access("NW")
+    assert a._staggered_offset(a.offsets, a.index[0]) == [
+        sp.Rational(-1, 2),
+        sp.Rational(1, 2),
+    ]
+
+    # sign reversed when using as flux field
+    r = ps.fields("r(2) : double[2D]", field_type=FieldType.STAGGERED_FLUX)
+    assert r[0, 0](0) == r.staggered_access("W")
+    assert -r[1, 0](0) == r.staggered_access("E")

tests/test_field_access_poly.py

+#
+# Copyright © 2020 Stephan Seitz <stephan.seitz@fau.de>
+#
+# Distributed under terms of the GPLv3 license.
+
+"""
+
+"""
+
+import pytest
+from pystencils.session import *
+
+from sympy import poly
+
+
+def test_field_access_poly():
+    dh = ps.create_data_handling((20, 20))
+    ρ = dh.add_array('rho')
+    rho = ρ.center
+    a = poly(rho+0.5, rho)
+    print(a)
+
+
+def test_field_access_piecewise():
+    try:
+        a = sp.Piecewise((0, 1 < sp.Max(-0.5, sp.Symbol("test") + 0.5)), (1, True))
+        a.simplify()
+    except Exception as e:
+        pytest.skip(f"Bug in SymPy 1.10: {e}")
+    else:
+        dh = ps.create_data_handling((20, 20))
+        ρ = dh.add_array('rho')
+        pw = sp.Piecewise((0, 1 < sp.Max(-0.5, ρ.center+0.5)), (1, True))
+        a = sp.simplify(pw)
+        print(a)

pystencils_tests/test_field_equality.ipynb → tests/test_field_equality.ipynb

 %% Cell type:code id: tags:
 
 ``` python
 from pystencils.session import *
-from pystencils.data_types import cast_func
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Test field equality behaviour
 
 %% Cell type:markdown id: tags:
 
 Fields create with same parameters are equal
 
 %% Cell type:code id: tags:
 
 ``` python
 f1 = ps.Field.create_generic('f', spatial_dimensions=2, index_dimensions=0)
 f2 = ps.Field.create_generic('f', spatial_dimensions=2, index_dimensions=0)
 
 assert f1 == f2
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print("Field ids equal in accesses: ", id(f1.center._field) == id(f2.center._field))
 print("Field accesses equal: ", f1.center == f2.center)
 ```
 
 %% Output
 
     Field ids equal in accesses:  True
     Field accesses equal:  True
 
 %% Cell type:code id: tags:
 
 ``` python
 f1 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
 f2 = ps.Field.create_generic('f', spatial_dimensions=2, index_dimensions=0)
 assert f1 != f2
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 f1 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
 f2 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0, dtype=np.float32)
 assert f1 != f2
 ```
 
 %% Cell type:markdown id: tags:
 
 Properties of fields:
 - `field_type`: enum distinguishing normal from index or buffer fields
 - `_dtype`: data type of field elements
 - `_layout`: tuple indicating the memory linearization order
 - `shape`: size of field for each dimension
 - `strides`: number of elements to jump over to increase coordinate of this dimension by one
 - `latex_name`: optional display name when field is printed as latex
 
 Equality compare of fields:
 - field has `__eq__` and ``__hash__`` overridden
 - all parameter but `latex_name` are considered for equality
 
 %% Cell type:markdown id: tags:
 
 ## Test field access equality behaviour
 
 %% Cell type:code id: tags:
 
 ``` python
 f1 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
 f2 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
 assert f1.center == f2.center
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 f1 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0)
 f2 = ps.Field.create_generic('f', spatial_dimensions=1, index_dimensions=0, dtype=np.float32)
 assert f1.center != f2.center
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def print_field_accesses_debug(expr):
     from pystencils import Field
     fas = list(expr.atoms(Field.Access))
     fields = list({e.field for e in fas})
     print("Field Accesses:")
     for fa in fas:
         print(f"   - {fa}, hash {hash(fa)}, offsets {fa._offsets}, index {fa._index}, {fa._hashable_content()}")
     print("")
     for i in range(len(fas)):
         for j in range(len(fas)):
             if not i < j:
                 continue
             print( f"   -> {i},{j}  {fas[i]} == {fas[j]}: {fas[i] == {fas[j]}}")
 
     print("Fields")
     for f in fields:
         print(f"  - {f}, {id(f)}, shape {f.shape}, strides {f.strides}, {f._dtype}, {f.field_type}, layout {f.layout}")
     print("")
     for i in range(len(fields)):
         for j in range(len(fields)):
             if not i < j:
                 continue
             print(f"  - {fields[i]} == {fields[j]}: {fields[i] == fields[j]}, ids equal {id(fields[i])==id(fields[j])}, hash equal {hash(fields[i])==hash(fields[j])}")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print_field_accesses_debug(f1.center * f2.center)
 ```
 
 %% Output
 
     Field Accesses:
-       - f[0], hash -3276894289571194847, offsets (0,), index (), (('f_C', ('commutative', True)), ((0,), (_size_f_0,), (_stride_f_0,), 3146377891102027609, <FieldType.GENERIC: 0>, 'f', None), 0)
-       - f[0], hash -1516451775709390846, offsets (0,), index (), (('f_C', ('commutative', True)), ((0,), (_size_f_0,), (_stride_f_0,), -1421177580377734245, <FieldType.GENERIC: 0>, 'f', None), 0)
+       - f[0], hash -8859424145258271267, offsets (0,), index (), ((('f_C', ('commutative', True), ('complex', True), ('extended_real', True), ('finite', True), ('hermitian', True), ('imaginary', False), ('infinite', False), ('real', True)), 2305067722319023373), ((0,), (_size_f_0,), (_stride_f_0,), <FieldType.GENERIC: 0>, 'f', None, double), 0)
+       - f[0], hash -6454673863007224785, offsets (0,), index (), ((('f_C', ('commutative', True), ('complex', True), ('extended_real', True), ('finite', True), ('hermitian', True), ('imaginary', False), ('infinite', False), ('real', True)), 4093629613697528859), ((0,), (_size_f_0,), (_stride_f_0,), <FieldType.GENERIC: 0>, 'f', None, float), 0)
     
        -> 0,1  f[0] == f[0]: False
     Fields
-      - f, 140548694371968, shape (_size_f_0,), strides (_stride_f_0,), double, FieldType.GENERIC, layout (0,)
-      - f, 140548693963104, shape (_size_f_0,), strides (_stride_f_0,), float, FieldType.GENERIC, layout (0,)
+      - f, 4881406800, shape (_size_f_0,), strides (_stride_f_0,), double, FieldType.GENERIC, layout (0,)
+      - f, 4881445024, shape (_size_f_0,), strides (_stride_f_0,), float, FieldType.GENERIC, layout (0,)
     
       - f == f: False, ids equal False, hash equal False

pystencils_tests/test_finite_differences.py → tests/test_finite_differences.py

 import sympy as sp
+import pytest
+
 import pystencils as ps
 from pystencils.astnodes import LoopOverCoordinate
-from pystencils.fd.spatial import fd_stencils_standard, fd_stencils_isotropic, discretize_spatial
-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, \
+    fd_stencils_forth_order_isotropic
 
 
 def test_spatial_2d_unit_sum():
@@ -16,7 +19,7 @@ def test_spatial_2d_unit_sum():
              diff(f, 1, 1),
              diff(f, 0, 0) + diff(f, 1, 1)]
 
-    schemes = [fd_stencils_standard, fd_stencils_isotropic]
+    schemes = [fd_stencils_standard, fd_stencils_isotropic, 'standard', 'isotropic']
 
     for term in terms:
         for scheme in schemes:
@@ -32,7 +35,7 @@ def test_spatial_1d_unit_sum():
     terms = [diff(f, 0),
              diff(f, 0, 0)]
 
-    schemes = [fd_stencils_standard, fd_stencils_isotropic]
+    schemes = [fd_stencils_standard, fd_stencils_isotropic, 'standard', 'isotropic']
 
     for term in terms:
         for scheme in schemes:
@@ -41,6 +44,17 @@ def test_spatial_1d_unit_sum():
             assert sum(coefficients) == 0
 
 
+def test_fd_stencils_forth_order_isotropic():
+    f = ps.fields("f: double[2D]")
+    a = fd_stencils_forth_order_isotropic([0], 1, f[0, 0](0))
+    sten, coefficients = ps.stencil.coefficients(a)
+    assert sum(coefficients) == 0
+
+    for i, direction in enumerate(sten):
+        counterpart = sten.index((direction[0] * -1, direction[1] * -1))
+        assert coefficients[i] + coefficients[counterpart] == 0
+
+
 def test_staggered_laplacian():
     f = ps.fields("f : double[2D]")
     a, dx = sp.symbols("a, dx")
@@ -68,3 +82,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)

tests/test_floor_ceil_int_optimization.py

+# -*- coding: utf-8 -*-
+#
+# Copyright © 2019 Stephan Seitz <stephan.seitz@fau.de>
+#
+# Distributed under terms of the GPLv3 license.
+
+"""
+
+"""
+
+import sympy as sp
+
+import pystencils
+from pystencils.typing import create_type
+
+
+def test_floor_ceil_int_optimization():
+    x, y = pystencils.fields('x,y: int32[2d]')
+    a, b, c = sp.symbols('a, b, c')
+    int_symbol = sp.Symbol('int_symbol', integer=True)
+    typed_symbol = pystencils.TypedSymbol('typed_symbol', create_type('int64'))
+
+    assignments = pystencils.AssignmentCollection({
+        a:  sp.floor(1),
+        b:  sp.ceiling(typed_symbol),
+        c:  sp.floor(int_symbol),
+        y.center():  sp.ceiling(x.center()) + sp.floor(x.center())
+    })
+
+    assert(typed_symbol.is_integer)
+    print(sp.simplify(sp.ceiling(typed_symbol)))
+
+    print(assignments)
+
+    wild_floor = sp.floor(sp.Wild('w1'))
+
+    assert not sp.floor(int_symbol).match(wild_floor)
+    assert sp.floor(a).match(wild_floor)
+
+    assert not assignments.find(wild_floor)
+
+
+def test_floor_ceil_float_no_optimization():
+    x, y = pystencils.fields('x,y: float32[2d]')
+    a, b, c = sp.symbols('a, b, c')
+    int_symbol = sp.Symbol('int_symbol', integer=True)
+    typed_symbol = pystencils.TypedSymbol('typed_symbol', create_type('float32'))
+
+    assignments = pystencils.AssignmentCollection({
+        a:  sp.floor(1),
+        b:  sp.ceiling(typed_symbol),
+        c:  sp.floor(int_symbol),
+        y.center():  sp.ceiling(x.center()) + sp.floor(x.center())
+    })
+
+    assert not typed_symbol.is_integer
+    print(sp.simplify(sp.ceiling(typed_symbol)))
+
+    print(assignments)
+
+    wild_floor = sp.floor(sp.Wild('w1'))
+
+    assert not sp.floor(int_symbol).match(wild_floor)
+    assert sp.floor(a).match(wild_floor)
+
+    assert assignments.find(wild_floor)