lbmpy_tests/advanced_streaming/test_advanced_streaming_noslip.py → tests/advanced_streaming/test_advanced_streaming_noslip.py

@@ -11,7 +11,7 @@ from lbmpy.stencils import LBStencil
 import pystencils as ps
 
 from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
-from pystencils.typing import TypedSymbol, create_type
+from pystencils import TypedSymbol
 from pystencils.field import Field, FieldType
 
 import pytest
@@ -43,7 +43,7 @@ def test_advanced_streaming_noslip_single_cell(stencil, streaming_pattern, prev_
     index_struct_dtype = numpy_data_type_for_boundary_object(noslip, stencil.D)
 
     index_field = Field('indexVector', FieldType.INDEXED, index_struct_dtype, layout=[0],
-                        shape=(TypedSymbol("indexVectorSize", create_type(np.int64)), 1), strides=(1, 1))
+                        shape=(TypedSymbol("indexVectorSize", np.int32), 1), strides=(1, 1))
     index_vector = np.array([pos + (d,) for d in range(stencil.Q)], dtype=index_struct_dtype)
 
     ast = create_lattice_boltzmann_boundary_kernel(pdf_field,

lbmpy_tests/advanced_streaming/test_communication.py → tests/advanced_streaming/test_communication.py

@@ -5,36 +5,52 @@ import numpy as np
 from lbmpy.stencils import LBStencil
 from pystencils.slicing import get_slice_before_ghost_layer, get_ghost_region_slice
 from lbmpy.creationfunctions import create_lb_update_rule, LBMConfig, LBMOptimisation
-from lbmpy.advanced_streaming.communication import get_communication_slices, _fix_length_one_slices, \
-    LBMPeriodicityHandling
+from lbmpy.advanced_streaming.communication import (
+    get_communication_slices,
+    _fix_length_one_slices,
+    LBMPeriodicityHandling,
+    periodic_pdf_gpu_copy_kernel,
+)
 from lbmpy.advanced_streaming.utility import streaming_patterns, Timestep
 from lbmpy.enums import Stencil
 
 import pytest
 
 
-@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27])
-@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
-@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD])
+@pytest.mark.parametrize(
+    "stencil", [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]
+)
+@pytest.mark.parametrize("streaming_pattern", streaming_patterns)
+@pytest.mark.parametrize("timestep", [Timestep.EVEN, Timestep.ODD])
 def test_slices_not_empty(stencil, streaming_pattern, timestep):
     stencil = LBStencil(stencil)
     arr = np.zeros((4,) * stencil.D + (stencil.Q,))
-    slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep,
-                                      ghost_layers=1)
+    slices = get_communication_slices(
+        stencil,
+        streaming_pattern=streaming_pattern,
+        prev_timestep=timestep,
+        ghost_layers=1,
+    )
     for _, slices_list in slices.items():
         for src, dst in slices_list:
             assert all(s != 0 for s in arr[src].shape)
             assert all(s != 0 for s in arr[dst].shape)
 
 
-@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27])
-@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
-@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD])
+@pytest.mark.parametrize(
+    "stencil", [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]
+)
+@pytest.mark.parametrize("streaming_pattern", streaming_patterns)
+@pytest.mark.parametrize("timestep", [Timestep.EVEN, Timestep.ODD])
 def test_src_dst_same_shape(stencil, streaming_pattern, timestep):
     stencil = LBStencil(stencil)
     arr = np.zeros((4,) * stencil.D + (stencil.Q,))
-    slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep,
-                                      ghost_layers=1)
+    slices = get_communication_slices(
+        stencil,
+        streaming_pattern=streaming_pattern,
+        prev_timestep=timestep,
+        ghost_layers=1,
+    )
     for _, slices_list in slices.items():
         for src, dst in slices_list:
             src_shape = arr[src].shape
@@ -42,12 +58,15 @@ def test_src_dst_same_shape(stencil, streaming_pattern, timestep):
             assert src_shape == dst_shape
 
 
-@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.parametrize(
+    "stencil", [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]
+)
 def test_pull_communication_slices(stencil):
     stencil = LBStencil(stencil)
 
     slices = get_communication_slices(
-        stencil, streaming_pattern='pull', prev_timestep=Timestep.BOTH, ghost_layers=1)
+        stencil, streaming_pattern="pull", prev_timestep=Timestep.BOTH, ghost_layers=1
+    )
     for i, d in enumerate(stencil):
         if i == 0:
             continue
@@ -58,21 +77,115 @@ def test_pull_communication_slices(stencil):
                 dst = s[1][:-1]
                 break
 
-        inner_slice = _fix_length_one_slices(get_slice_before_ghost_layer(d, ghost_layers=1))
+        inner_slice = _fix_length_one_slices(
+            get_slice_before_ghost_layer(d, ghost_layers=1)
+        )
         inv_dir = (-e for e in d)
-        gl_slice = _fix_length_one_slices(get_ghost_region_slice(inv_dir, ghost_layers=1))
+        gl_slice = _fix_length_one_slices(
+            get_ghost_region_slice(inv_dir, ghost_layers=1)
+        )
         assert src == inner_slice
         assert dst == gl_slice
 
 
-@pytest.mark.parametrize('stencil_name', [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.parametrize("direction", LBStencil(Stencil.D3Q27).stencil_entries)
+@pytest.mark.parametrize("pull", [False, True])
+def test_gpu_comm_kernels(direction: tuple, pull: bool):
+    pytest.importorskip("cupy")
+
+    stencil = LBStencil(Stencil.D3Q27)
+    inv_dir = stencil[stencil.inverse_index(direction)]
+    target = ps.Target.GPU
+
+    domain_size = (4,) * stencil.D
+
+    dh: ps.datahandling.SerialDataHandling = ps.create_data_handling(
+        domain_size,
+        periodicity=(True,) * stencil.D,
+        parallel=False,
+        default_target=target,
+    )
+
+    field = dh.add_array("field", values_per_cell=2)
+
+    if pull:
+        dst_slice = get_ghost_region_slice(inv_dir)
+        src_slice = get_slice_before_ghost_layer(direction)
+    else:
+        dst_slice = get_slice_before_ghost_layer(direction)
+        src_slice = get_ghost_region_slice(inv_dir)
+
+    src_slice += (1,)
+    dst_slice += (1,)
+
+    kernel = periodic_pdf_gpu_copy_kernel(field, src_slice, dst_slice)
+
+    dh.cpu_arrays[field.name][src_slice] = 42.0
+    dh.all_to_gpu()
+
+    dh.run_kernel(kernel)
+
+    dh.all_to_cpu()
+    np.testing.assert_equal(dh.cpu_arrays[field.name][dst_slice], 42.0)
+
+
+@pytest.mark.parametrize("stencil", [Stencil.D2Q9, Stencil.D3Q19])
+@pytest.mark.parametrize("streaming_pattern", streaming_patterns)
+def test_direct_copy_and_kernels_equivalence(stencil: Stencil, streaming_pattern: str):
+    pytest.importorskip("cupy")
+
+    target = ps.Target.GPU
+    stencil = LBStencil(stencil)
+    domain_size = (4,) * stencil.D
+
+    dh: ps.datahandling.SerialDataHandling = ps.create_data_handling(
+        domain_size,
+        periodicity=(True,) * stencil.D,
+        parallel=False,
+        default_target=target,
+    )
+
+    pdfs_a = dh.add_array("pdfs_a", values_per_cell=stencil.Q)
+    pdfs_b = dh.add_array("pdfs_b", values_per_cell=stencil.Q)
+
+    dh.fill(pdfs_a.name, 0.0, ghost_layers=True)
+    dh.fill(pdfs_b.name, 0.0, ghost_layers=True)
+
+    for q in range(stencil.Q):
+        sl = ps.make_slice[:4, :4, q] if stencil.D == 2 else ps.make_slice[:4, :4, :4, q]
+        dh.cpu_arrays[pdfs_a.name][sl] = q
+        dh.cpu_arrays[pdfs_b.name][sl] = q
+
+    dh.all_to_gpu()
+
+    direct_copy = LBMPeriodicityHandling(stencil, dh, pdfs_a.name, streaming_pattern, cupy_direct_copy=True)
+    kernels_copy = LBMPeriodicityHandling(stencil, dh, pdfs_b.name, streaming_pattern, cupy_direct_copy=False)
+
+    direct_copy(Timestep.EVEN)
+    kernels_copy(Timestep.EVEN)
+
+    dh.all_to_cpu()
+
+    np.testing.assert_equal(
+        dh.cpu_arrays[pdfs_a.name],
+        dh.cpu_arrays[pdfs_b.name]
+    )
+
+
+@pytest.mark.parametrize(
+    "stencil_name", [Stencil.D2Q9, Stencil.D3Q15, Stencil.D3Q19, Stencil.D3Q27]
+)
 def test_optimised_and_full_communication_equivalence(stencil_name):
     target = ps.Target.CPU
     stencil = LBStencil(stencil_name)
-    domain_size = (4, ) * stencil.D
+    domain_size = (4,) * stencil.D
 
-    dh = ps.create_data_handling(domain_size, periodicity=(True, ) * stencil.D,
-                                 parallel=False, default_target=target)
+    dh = ps.create_data_handling(
+        domain_size,
+        periodicity=(True,) * stencil.D,
+        parallel=False,
+        default_target=target,
+    )
 
     pdf = dh.add_array("pdf", values_per_cell=len(stencil), dtype=np.int64)
     dh.fill("pdf", 0, ghost_layers=True)
@@ -82,34 +195,40 @@ def test_optimised_and_full_communication_equivalence(stencil_name):
     gl = dh.ghost_layers_of_field("pdf")
 
     num = 0
-    for idx, x in np.ndenumerate(dh.cpu_arrays['pdf']):
-        dh.cpu_arrays['pdf'][idx] = num
-        dh.cpu_arrays['pdf_tmp'][idx] = num
+    for idx, x in np.ndenumerate(dh.cpu_arrays["pdf"]):
+        dh.cpu_arrays["pdf"][idx] = num
+        dh.cpu_arrays["pdf_tmp"][idx] = num
         num += 1
 
     lbm_config = LBMConfig(stencil=stencil, kernel_type="stream_pull_only")
     lbm_opt = LBMOptimisation(symbolic_field=pdf, symbolic_temporary_field=pdf_tmp)
-    config = ps.CreateKernelConfig(target=dh.default_target, cpu_openmp=True)
+    config = ps.CreateKernelConfig(target=dh.default_target, data_type=np.int64, cpu_openmp=True)
 
     ac = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
     ast = ps.create_kernel(ac, config=config)
     stream = ast.compile()
 
-    full_communication = dh.synchronization_function(pdf.name, target=dh.default_target, optimization={"openmp": True})
+    full_communication = dh.synchronization_function(
+        pdf.name, target=dh.default_target, optimization={"openmp": True}
+    )
     full_communication()
 
     dh.run_kernel(stream)
     dh.swap("pdf", "pdf_tmp")
-    pdf_full_communication = np.copy(dh.cpu_arrays['pdf'])
+    pdf_full_communication = np.copy(dh.cpu_arrays["pdf"])
 
     num = 0
-    for idx, x in np.ndenumerate(dh.cpu_arrays['pdf']):
-        dh.cpu_arrays['pdf'][idx] = num
-        dh.cpu_arrays['pdf_tmp'][idx] = num
+    for idx, x in np.ndenumerate(dh.cpu_arrays["pdf"]):
+        dh.cpu_arrays["pdf"][idx] = num
+        dh.cpu_arrays["pdf_tmp"][idx] = num
         num += 1
 
-    optimised_communication = LBMPeriodicityHandling(stencil=stencil, data_handling=dh, pdf_field_name=pdf.name,
-                                                     streaming_pattern='pull')
+    optimised_communication = LBMPeriodicityHandling(
+        stencil=stencil,
+        data_handling=dh,
+        pdf_field_name=pdf.name,
+        streaming_pattern="pull",
+    )
     optimised_communication()
     dh.run_kernel(stream)
     dh.swap("pdf", "pdf_tmp")
@@ -119,9 +238,14 @@ def test_optimised_and_full_communication_equivalence(stencil_name):
             for j in range(gl, domain_size[1]):
                 for k in range(gl, domain_size[2]):
                     for f in range(len(stencil)):
-                        assert dh.cpu_arrays['pdf'][i, j, k, f] == pdf_full_communication[i, j, k, f], print(f)
+                        assert (
+                            dh.cpu_arrays["pdf"][i, j, k, f]
+                            == pdf_full_communication[i, j, k, f]
+                        ), print(f)
     else:
         for i in range(gl, domain_size[0]):
             for j in range(gl, domain_size[1]):
                 for f in range(len(stencil)):
-                    assert dh.cpu_arrays['pdf'][i, j, f] == pdf_full_communication[i, j, f]
+                    assert (
+                        dh.cpu_arrays["pdf"][i, j, f] == pdf_full_communication[i, j, f]
+                    )

lbmpy_tests/advanced_streaming/test_fully_periodic_flow.py → tests/advanced_streaming/test_fully_periodic_flow.py

@@ -20,7 +20,7 @@ all_results = dict()
 targets = [Target.CPU]
 
 try:
-    import pycuda.autoinit
+    import cupy
     targets += [Target.GPU]
 except Exception:
     pass

lbmpy_tests/advanced_streaming/test_periodic_pipe_with_force.py → tests/advanced_streaming/test_periodic_pipe_with_force.py

@@ -21,7 +21,7 @@ all_results = dict()
 targets = [Target.CPU]
 
 try:
-    import pycuda.autoinit
+    import cupy
     targets += [Target.GPU]
 except Exception:
     pass

tests/benchmark/db/config.json

+{"indexes": {"simulationresult": {"pk": {"key": "pk", "id": "76f7f2660bd4475b90697716ef655a46"}}}, "store_class": "transactional", "index_class": "transactional", "index_store_class": "basic", "serializer_class": "json", "autocommit": false, "version": "0.2.12"}
\ No newline at end of file

lbmpy_tests/benchmark/test_benchmark.py → tests/benchmark/test_benchmark.py

@@ -2,9 +2,16 @@ from statistics import median
 
 import numpy as np
 import sympy as sp
+import pytest
 
 from lbmpy.scenarios import create_lid_driven_cavity
-from pystencils.cpu.cpujit import add_or_change_compiler_flags
+from lbmpy._compat import IS_PYSTENCILS_2
+
+if IS_PYSTENCILS_2:
+    pytest.skip(reason="Benchmark test case needs to be rebuilt for pystencils 2.0", allow_module_level=True)
+else:
+    from pystencils.cpu.cpujit import add_or_change_compiler_flags
+
 from pystencils.runhelper import ParameterStudy
 
 

lbmpy_tests/centeredcumulant/test_equilibrium.py → tests/cumulantmethod/test_equilibrium.py

@@ -13,7 +13,8 @@ from lbmpy.methods import CollisionSpaceInfo
 from lbmpy.moment_transforms import (
     FastCentralMomentTransform,
     PdfsToCentralMomentsByMatrix,
-    PdfsToCentralMomentsByShiftMatrix)
+    PdfsToCentralMomentsByShiftMatrix,
+    BinomialChimeraTransform)
 
 sympy_numeric_version = [int(x, 10) for x in sp.__version__.split('.')]
 if len(sympy_numeric_version) < 3:
@@ -30,7 +31,8 @@ reference_equilibria = dict()
 @pytest.mark.parametrize('stencil_name', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
 @pytest.mark.parametrize('cm_transform', [PdfsToCentralMomentsByMatrix,
                                           FastCentralMomentTransform,
-                                          PdfsToCentralMomentsByShiftMatrix])
+                                          PdfsToCentralMomentsByShiftMatrix,
+                                          BinomialChimeraTransform])
 def test_equilibrium_pdfs(stencil_name, cm_transform):
     stencil = LBStencil(stencil_name)
     cspace = CollisionSpaceInfo(CollisionSpace.CUMULANTS, central_moment_transform_class=cm_transform)

lbmpy_tests/centeredcumulant/test_flow_around_sphere.py → tests/cumulantmethod/test_flow_around_sphere.py

@@ -16,7 +16,7 @@ from pystencils import create_kernel, create_data_handling, Assignment, Target,
 from pystencils.slicing import slice_from_direction, get_slice_before_ghost_layer
 
 
-def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
+def flow_around_sphere(stencil, galilean_correction, fourth_order_correction, L_LU, total_steps):
     if galilean_correction and stencil.Q != 27:
         pytest.skip()
 
@@ -38,7 +38,7 @@ def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
 
     lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, relaxation_rate=omega_v,
                            compressible=True,
-                           galilean_correction=galilean_correction)
+                           galilean_correction=galilean_correction, fourth_order_correction=fourth_order_correction)
 
     lbm_opt = LBMOptimisation(pre_simplification=True)
     config = CreateKernelConfig(target=target)
@@ -58,7 +58,12 @@ def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
             boundary_assignments)
         iter_slice = get_slice_before_ghost_layer((1,) + (0,) * (stencil.D - 1))
         extrapolation_ast = create_kernel(
-            boundary_assignments, config=CreateKernelConfig(iteration_slice=iter_slice, ghost_layers=1, target=target))
+            boundary_assignments,
+            config=CreateKernelConfig(
+                iteration_slice=iter_slice,
+                target=target,
+                skip_independence_check=True)
+            )
         return extrapolation_ast.compile()
 
     dh = create_data_handling(channel_size, periodicity=False, default_layout='fzyx', default_target=target)
@@ -135,14 +140,22 @@ def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
 
 @pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
 @pytest.mark.parametrize('galilean_correction', [False, True])
-def test_flow_around_sphere_short(stencil, galilean_correction):
-    pytest.importorskip('pycuda')
-    flow_around_sphere(LBStencil(stencil), galilean_correction, 5, 200)
+@pytest.mark.parametrize('fourth_order_correction', [False, True])
+def test_flow_around_sphere_short(stencil, galilean_correction, fourth_order_correction):
+    pytest.importorskip('cupy')
+    stencil = LBStencil(stencil)
+    if fourth_order_correction and stencil.Q != 27:
+        pytest.skip("Fourth-order correction only defined for D3Q27 stencil.")
+    flow_around_sphere(stencil, galilean_correction, fourth_order_correction, 5, 200)
 
 
 @pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
 @pytest.mark.parametrize('galilean_correction', [False, True])
+@pytest.mark.parametrize('fourth_order_correction', [False, 0.01])
 @pytest.mark.longrun
-def test_flow_around_sphere_long(stencil, galilean_correction):
-    pytest.importorskip('pycuda')
-    flow_around_sphere(LBStencil(stencil), galilean_correction, 20, 3000)
+def test_flow_around_sphere_long(stencil, galilean_correction, fourth_order_correction):
+    pytest.importorskip('cupy')
+    stencil = LBStencil(stencil)
+    if fourth_order_correction and stencil.Q != 27:
+        pytest.skip("Fourth-order correction only defined for D3Q27 stencil.")
+    flow_around_sphere(stencil, galilean_correction, fourth_order_correction, 20, 3000)

tests/cumulantmethod/test_periodic_pipe_flow.ipynb

+%% Cell type:code id: tags:
+
+``` python
+%load_ext autoreload
+%autoreload 2
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from dataclasses import replace
+
+import sympy as sp
+import numpy as np
+from numpy.testing import assert_allclose
+
+from pystencils.session import *
+from pystencils import Target
+from lbmpy.session import *
+
+from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
+```
+
+%% Cell type:code id: tags:
+
+``` python
+try:
+    import cupy
+except ImportError:
+    cupy = None
+    target = ps.Target.CPU
+    print('No cupy installed')
+
+if cupy:
+    target = ps.Target.GPU
+```
+
+%% Cell type:markdown id: tags:
+
+## Pipe Flow Scenario
+
+%% Cell type:code id: tags:
+
+``` python
+class PeriodicPipeFlow:
+    def __init__(self, lbm_config, lbm_optimisation, config):
+
+        wall_boundary = NoSlip()
+
+        self.stencil = lbm_config.stencil
+        self.streaming_pattern = lbm_config.streaming_pattern
+        self.target = config.target
+
+        self.gpu = self.target == Target.GPU
+
+        #   Stencil
+        self.q = stencil.Q
+        self.dim = stencil.D
+
+        #   Streaming
+        self.inplace = is_inplace(self.streaming_pattern)
+        self.timesteps = get_timesteps(self.streaming_pattern)
+        self.zeroth_timestep = self.timesteps[0]
+
+        #   Domain, Data Handling and PDF fields
+        self.pipe_length = 60
+        self.pipe_radius = 15
+        self.domain_size = (self.pipe_length, ) + (2 * self.pipe_radius,) * (self.dim - 1)
+        self.periodicity = (True, ) + (False, ) * (self.dim - 1)
+        force = (0.0001, ) + (0.0,) * (self.dim - 1)
+        self.force = lbm_config.force
+
+        self.dh = create_data_handling(domain_size=self.domain_size,
+                                       periodicity=self.periodicity, default_target=self.target)
+
+        self.pdfs = self.dh.add_array('pdfs', self.q)
+        if not self.inplace:
+            self.pdfs_tmp = self.dh.add_array_like('pdfs_tmp', self.pdfs.name)
+
+        lbm_optimisation = replace(lbm_optimisation, symbolic_field=self.pdfs)
+
+        if not self.inplace:
+            lbm_optimisation = replace(lbm_optimisation, symbolic_temporary_field=self.pdfs_tmp)
+
+        self.lb_collision = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_optimisation)
+        self.lb_method = self.lb_collision.method
+
+        self.lb_kernels = []
+        for t in self.timesteps:
+            lbm_config = replace(lbm_config, timestep=t, collision_rule=self.lb_collision)
+            self.lb_kernels.append(create_lb_function(lbm_config=lbm_config,
+                                                      lbm_optimisation=lbm_optimisation,
+                                                      config=config))
+
+        #   Macroscopic Values
+        self.density = 1.0
+        self.density_field = self.dh.add_array('rho', 1)
+        u_x = 0.0
+        self.velocity = (u_x,) * self.dim
+        self.velocity_field = self.dh.add_array('u', self.dim)
+
+        single_gl_config = replace(config, ghost_layers=1)
+        setter = macroscopic_values_setter(
+            self.lb_method, self.density, self.velocity, self.pdfs,
+            streaming_pattern=self.streaming_pattern, previous_timestep=self.zeroth_timestep)
+        self.init_kernel = create_kernel(setter, config=single_gl_config).compile()
+
+        self.getter_kernels = []
+        for t in self.timesteps:
+            getter = macroscopic_values_getter(
+                self.lb_method, self.density_field, self.velocity_field, self.pdfs,
+                streaming_pattern=self.streaming_pattern, previous_timestep=t)
+            self.getter_kernels.append(create_kernel(getter, config=single_gl_config).compile())
+
+        #   Periodicity
+        self.periodicity_handler = LBMPeriodicityHandling(
+            self.stencil, self.dh, self.pdfs.name, streaming_pattern=self.streaming_pattern)
+
+        #   Boundary Handling
+        self.wall = wall_boundary
+        self.bh = LatticeBoltzmannBoundaryHandling(
+            self.lb_method, self.dh, self.pdfs.name,
+            streaming_pattern=self.streaming_pattern, target=self.target)
+
+        self.bh.set_boundary(boundary_obj=self.wall, mask_callback=self.mask_callback)
+
+        self.current_timestep = self.zeroth_timestep
+
+    def mask_callback(self, x, y, z=None):
+        y = y - self.pipe_radius
+        z = z - self.pipe_radius if z is not None else 0
+        return np.sqrt(y**2 + z**2) >= self.pipe_radius
+
+    def init(self):
+        self.current_timestep = self.zeroth_timestep
+        self.dh.run_kernel(self.init_kernel)
+
+    def step(self):
+        #   Order matters! First communicate, then boundaries, otherwise
+        #   periodicity handling overwrites reflected populations
+        # Periodicty
+        self.periodicity_handler(self.current_timestep)
+
+        # Boundaries
+        self.bh(prev_timestep=self.current_timestep)
+
+        # Here, the next time step begins
+        self.current_timestep = self.current_timestep.next()
+
+        # LBM Step
+        self.dh.run_kernel(self.lb_kernels[self.current_timestep.idx])
+
+        # Field Swaps
+        if not self.inplace:
+            self.dh.swap(self.pdfs.name, self.pdfs_tmp.name)
+
+        # Macroscopic Values
+        self.dh.run_kernel(self.getter_kernels[self.current_timestep.idx])
+
+    def run(self, iterations):
+        for _ in range(iterations):
+            self.step()
+
+    @property
+    def velocity_array(self):
+        if self.gpu:
+            self.dh.to_cpu(self.velocity_field.name)
+        return self.dh.gather_array(self.velocity_field.name)
+
+    def get_trimmed_velocity_array(self):
+        if self.gpu:
+            self.dh.to_cpu(self.velocity_field.name)
+        u = np.copy(self.dh.gather_array(self.velocity_field.name))
+        mask = self.bh.get_mask(None, self.wall)
+        for idx in np.ndindex(u.shape[:-1]):
+            if mask[idx] != 0:
+                u[idx] = np.full((self.dim, ), np.nan)
+        return u
+```
+
+%% Cell type:markdown id: tags:
+
+## General Setup
+
+%% Cell type:code id: tags:
+
+``` python
+stencil = LBStencil(Stencil.D3Q19)
+dim = stencil.D
+streaming_pattern = 'aa'
+
+config = ps.CreateKernelConfig(target=target)
+
+viscous_rr = 1.1
+force = (0.0001, ) + (0.0,) * (dim - 1)
+```
+
+%% Cell type:markdown id: tags:
+
+## 1. Reference: SRT Method
+
+%% Cell type:code id: tags:
+
+``` python
+srt_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=viscous_rr,
+                       force_model=ForceModel.SIMPLE, force=force, streaming_pattern=streaming_pattern)
+
+srt_flow = PeriodicPipeFlow(srt_config, LBMOptimisation(), config)
+srt_flow.init()
+srt_flow.run(400)
+```
+
+%% Output
+
+    ---------------------------------------------------------------------------
+    AttributeError                            Traceback (most recent call last)
+Cell     In[6], line 4
+          1 srt_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=viscous_rr,
+          2                        force_model=ForceModel.SIMPLE, force=force, streaming_pattern=streaming_pattern)
+    ----> 4 srt_flow = PeriodicPipeFlow(srt_config, LBMOptimisation(), config)
+          5 srt_flow.init()
+          6 srt_flow.run(400)
+Cell     In[4], line 47, in PeriodicPipeFlow.__init__(self, lbm_config, lbm_optimisation, config)
+         45 for t in self.timesteps:
+         46     lbm_config = replace(lbm_config, timestep=t, collision_rule=self.lb_collision)
+    ---> 47     self.lb_kernels.append(create_lb_function(lbm_config=lbm_config,
+         48                                               lbm_optimisation=lbm_optimisation,
+         49                                               config=config))
+         51 #   Macroscopic Values
+         52 self.density = 1.0
+File     ~/pystencils/lbmpy/lbmpy/creationfunctions.py:505, in create_lb_function(ast, lbm_config, lbm_optimisation, config, optimization, **kwargs)
+        502     ast = lbm_config.ast
+        504 if ast is None:
+    --> 505     ast = create_lb_ast(lbm_config.update_rule, lbm_config=lbm_config,
+        506                         lbm_optimisation=lbm_optimisation, config=config)
+        508 res = ast.compile()
+        510 res.method = ast.method
+File     ~/pystencils/lbmpy/lbmpy/creationfunctions.py:530, in create_lb_ast(update_rule, lbm_config, lbm_optimisation, config, optimization, **kwargs)
+        525     update_rule = create_lb_update_rule(lbm_config.collision_rule, lbm_config=lbm_config,
+        526                                         lbm_optimisation=lbm_optimisation, config=config)
+        528 field_types = set(fa.field.dtype for fa in update_rule.defined_symbols if isinstance(fa, Field.Access))
+    --> 530 config = replace(config, data_type=collate_types(field_types), ghost_layers=1)
+        531 ast = create_kernel(update_rule, config=config)
+        533 ast.method = update_rule.method
+File     /usr/lib/python3.11/dataclasses.py:1492, in replace(obj, **changes)
+       1485         changes[f.name] = getattr(obj, f.name)
+       1487 # Create the new object, which calls __init__() and
+       1488 # __post_init__() (if defined), using all of the init fields we've
+       1489 # added and/or left in 'changes'.  If there are values supplied in
+       1490 # changes that aren't fields, this will correctly raise a
+       1491 # TypeError.
+    -> 1492 return obj.__class__(**changes)
+File     <string>:24, in __init__(self, target, backend, function_name, data_type, default_number_float, default_number_int, iteration_slice, ghost_layers, cpu_openmp, cpu_vectorize_info, cpu_blocking, omp_single_loop, gpu_indexing, gpu_indexing_params, default_assignment_simplifications, cpu_prepend_optimizations, use_auto_for_assignments, index_fields, coordinate_names, allow_double_writes, skip_independence_check)
+File     ~/pystencils/pystencils/pystencils/config.py:177, in CreateKernelConfig.__post_init__(self)
+        174     self._check_type(dtype)
+        176 if not isinstance(self.data_type, dict):
+    --> 177     dt = copy(self.data_type)  # The copy is necessary because BasicType has sympy shinanigans
+        178     self.data_type = defaultdict(self.DataTypeFactory(dt))
+        180 if isinstance(self.data_type, dict) and not isinstance(self.data_type, defaultdict):
+File     /usr/lib/python3.11/copy.py:102, in copy(x)
+        100 if isinstance(rv, str):
+        101     return x
+    --> 102 return _reconstruct(x, None, *rv)
+File     /usr/lib/python3.11/copy.py:273, in _reconstruct(x, memo, func, args, state, listiter, dictiter, deepcopy)
+        271     state = deepcopy(state, memo)
+        272 if hasattr(y, '__setstate__'):
+    --> 273     y.__setstate__(state)
+        274 else:
+        275     if isinstance(state, tuple) and len(state) == 2:
+File     ~/.local/lib/python3.11/site-packages/sympy/core/basic.py:144, in Basic.__setstate__(self, state)
+        143 def __setstate__(self, state):
+    --> 144     for name, value in state.items():
+        145         setattr(self, name, value)
+    AttributeError: 'tuple' object has no attribute 'items'
+
+%% Cell type:code id: tags:
+
+``` python
+srt_u = srt_flow.get_trimmed_velocity_array()
+ps.plot.vector_field_magnitude(srt_u[30,:,:,:])
+ps.plot.colorbar()
+```
+
+%% Cell type:markdown id: tags:
+
+## 2. Centered Cumulant Method with Central Moment-Space Forcing
+
+%% Cell type:code id: tags:
+
+``` python
+cm_config = LBMConfig(method=Method.MONOMIAL_CUMULANT, stencil=stencil, compressible=True,
+                      relaxation_rate=viscous_rr,
+                      force=force, streaming_pattern=streaming_pattern)
+
+lbm_opt = LBMOptimisation(pre_simplification=True)
+
+cm_impl_f_flow = PeriodicPipeFlow(cm_config, lbm_opt, config)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+lb_method = cm_impl_f_flow.lb_method
+assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
+assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+cm_impl_f_flow.init()
+cm_impl_f_flow.run(400)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+cm_impl_f_u = cm_impl_f_flow.get_trimmed_velocity_array()
+ps.plot.vector_field_magnitude(cm_impl_f_u[30,:,:,:])
+ps.plot.colorbar()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+assert_allclose(cm_impl_f_u, srt_u, rtol=1, atol=1e-4)
+```
+
+%% Cell type:markdown id: tags:
+
+## 3. Centered Cumulant Method with Simple force model
+
+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.forcemodels import Simple
+
+cm_expl_force_config = LBMConfig(method=Method.CUMULANT, stencil=stencil, compressible=True,
+                      relaxation_rate=viscous_rr, force_model=Simple(force),
+                      force=force, streaming_pattern=streaming_pattern)
+
+lbm_opt = LBMOptimisation(pre_simplification=True)
+
+cm_expl_f_flow = PeriodicPipeFlow(cm_expl_force_config, lbm_opt, config)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+lb_method = cm_expl_f_flow.lb_method
+assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
+assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+cm_expl_f_flow.init()
+cm_expl_f_flow.run(400)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+cm_expl_f_u = cm_expl_f_flow.get_trimmed_velocity_array()
+ps.plot.vector_field_magnitude(cm_expl_f_u[30,:,:,:])
+ps.plot.colorbar()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+assert_allclose(cm_expl_f_u, srt_u, rtol=1, atol=1e-4)
+assert_allclose(cm_expl_f_u, cm_impl_f_u, rtol=1, atol=1e-4)
+```

lbmpy_tests/full_scenarios/schaefer_turek/scenario_schaefer_turek.py → tests/full_scenarios/schaefer_turek/scenario_schaefer_turek.py

@@ -10,7 +10,9 @@ import warnings
 import numpy as np
 import pytest
 
-from pystencils.backends.simd_instruction_sets import get_supported_instruction_sets
+from pystencils import Target
+
+from lbmpy._compat import get_supported_instruction_sets
 from lbmpy.boundaries.boundaryconditions import NoSlip
 from lbmpy.geometry import get_pipe_velocity_field
 from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity