src/lbmpy/turbulence_models.py

+import sympy as sp
+
+from lbmpy.relaxationrates import get_shear_relaxation_rate, relaxation_rate_from_lattice_viscosity, \
+    lattice_viscosity_from_relaxation_rate
+from lbmpy.utils import extract_shear_relaxation_rate, frobenius_norm, second_order_moment_tensor
+from pystencils import Assignment
+
+from lbmpy.enums import SubgridScaleModel
+
+
+def add_sgs_model(collision_rule, subgrid_scale_model: SubgridScaleModel, model_constant=None, omega_output_field=None,
+                  eddy_viscosity_field=None):
+    r""" Wrapper for SGS models to provide default constants and outsource SGS model handling from creation routines."""
+
+    if subgrid_scale_model == SubgridScaleModel.SMAGORINSKY:
+        model_constant = model_constant if model_constant else sp.Float(0.12)
+        return add_smagorinsky_model(collision_rule=collision_rule, smagorinsky_constant=model_constant,
+                                     omega_output_field=omega_output_field, eddy_viscosity_field=eddy_viscosity_field)
+    if subgrid_scale_model == SubgridScaleModel.QR:
+        model_constant = model_constant if model_constant else sp.Rational(1, 3)
+        return add_qr_model(collision_rule=collision_rule, qr_constant=model_constant,
+                            omega_output_field=omega_output_field, eddy_viscosity_field=eddy_viscosity_field)
+
+
+def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_field=None, eddy_viscosity_field=None):
+    r""" Adds a Smagorinsky model to a lattice Boltzmann collision rule. To add the Smagorinsky model to a LB scheme
+        one has to first compute the strain rate tensor $S_{ij}$ in each cell, and compute the turbulent
+        viscosity :math:`nu_t` from it. Then the local relaxation rate has to be adapted to match the total viscosity
+        :math `\nu_{total}` instead of the standard viscosity :math `\nu_0`.
+
+        A fortunate property of LB methods is, that the strain rate tensor can be computed locally from the
+        non-equilibrium part of the distribution function. This is somewhat surprising, since the strain rate tensor
+        contains first order derivatives. The strain rate tensor can be obtained by
+
+        .. math ::
+            S_{ij} = - \frac{3 \omega_s}{2 \rho_{(0)}} \Pi_{ij}^{(neq)}
+
+        where :math `\omega_s` is the relaxation rate that determines the viscosity, :math `\rho_{(0)}` is :math `\rho`
+        in compressible models and :math `1` for incompressible schemes.
+        :math `\Pi_{ij}^{(neq)}` is the second order moment tensor of the non-equilibrium part of
+        the distribution functions
+        :math `f^{(neq)} = f - f^{(eq)}` and can be computed as
+
+        .. math ::
+            \Pi_{ij}^{(neq)} = \sum_q c_{qi} c_{qj} \; f_q^{(neq)}
+
+
+    """
+    method = collision_rule.method
+    omega_s = get_shear_relaxation_rate(method)
+    omega_s, found_symbolic_shear_relaxation = extract_shear_relaxation_rate(collision_rule, omega_s)
+
+    if not found_symbolic_shear_relaxation:
+        raise ValueError("For the Smagorinsky model the shear relaxation rate has to be a symbol or it has to be "
+                         "assigned to a single equation in the assignment list")
+    f_neq = sp.Matrix(method.pre_collision_pdf_symbols) - method.get_equilibrium_terms()
+    equilibrium = method.equilibrium_distribution
+
+    tau_0 = sp.Symbol("tau_0_")
+    second_order_neq_moments = sp.Symbol("Pi")
+    rho = equilibrium.density if equilibrium.compressible else equilibrium.background_density
+    adapted_omega = sp.Symbol("sgs_omega")
+
+    collision_rule = collision_rule.new_with_substitutions({omega_s: adapted_omega}, substitute_on_lhs=False)
+    # for derivation see notebook demo_custom_LES_model.ipynb
+    eqs = [Assignment(tau_0, sp.Float(1) / omega_s),
+           Assignment(second_order_neq_moments,
+                      frobenius_norm(second_order_moment_tensor(f_neq, method.stencil), factor=2) / rho),
+           Assignment(adapted_omega,
+                      sp.Float(2) / (tau_0 + sp.sqrt(
+                          sp.Float(18) * smagorinsky_constant ** 2 * second_order_neq_moments + tau_0 ** 2)))]
+    collision_rule.subexpressions += eqs
+    collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
+
+    if eddy_viscosity_field:
+        collision_rule.main_assignments.append(Assignment(
+            eddy_viscosity_field.center,
+            smagorinsky_constant ** 2 * sp.Rational(3, 2) * adapted_omega * second_order_neq_moments
+        ))
+
+    if omega_output_field:
+        collision_rule.main_assignments.append(Assignment(omega_output_field.center, adapted_omega))
+
+    return collision_rule
+
+
+def add_qr_model(collision_rule, qr_constant, omega_output_field=None, eddy_viscosity_field=None):
+    r""" Adds a QR model to a lattice Boltzmann collision rule, see :cite:`rozema15`.
+    WARNING : this subgrid-scale model is only defined for isotropic grids
+    """
+    method = collision_rule.method
+    omega_s = get_shear_relaxation_rate(method)
+    omega_s, found_symbolic_shear_relaxation = extract_shear_relaxation_rate(collision_rule, omega_s)
+
+    if not found_symbolic_shear_relaxation:
+        raise ValueError("For the QR model the shear relaxation rate has to be a symbol or it has to be "
+                         "assigned to a single equation in the assignment list")
+    f_neq = sp.Matrix(method.pre_collision_pdf_symbols) - method.get_equilibrium_terms()
+    equilibrium = method.equilibrium_distribution
+
+    nu_0, nu_e = sp.symbols("qr_nu_0 qr_nu_e")
+    c_pi_s = sp.Symbol("qr_c_pi")
+    rho = equilibrium.density if equilibrium.compressible else equilibrium.background_density
+    adapted_omega = sp.Symbol("sgs_omega")
+
+    stencil = method.stencil
+
+    pi = second_order_moment_tensor(f_neq, stencil)
+
+    r_prime = sp.Symbol("qr_r_prime")
+    q_prime = sp.Symbol("qr_q_prime")
+
+    c_pi = qr_constant * sp.Piecewise(
+        (r_prime, sp.StrictGreaterThan(r_prime, sp.Float(0))),
+        (sp.Float(0), True)
+    ) / q_prime
+
+    collision_rule = collision_rule.new_with_substitutions({omega_s: adapted_omega}, substitute_on_lhs=False)
+
+    if stencil.D == 2:
+        nu_e_assignments = [Assignment(nu_e, c_pi_s)]
+    elif stencil.D == 3:
+        base_viscosity = sp.Symbol("qr_base_viscosity")
+        nu_e_assignments = [
+            Assignment(base_viscosity, sp.Float(6) * nu_0 + sp.Float(1)),
+            Assignment(nu_e, (-base_viscosity + sp.sqrt(base_viscosity ** 2 + sp.Float(72) * c_pi_s / rho))
+                       / sp.Float(12))
+        ]
+    else:
+        raise ValueError("QR-model is only defined for 2- or 3-dimensional flows")
+
+    matrix_entries = sp.Matrix(stencil.D, stencil.D, sp.symbols(f"sgs_qr_pi:{stencil.D ** 2}"))
+
+    eqs = [Assignment(nu_0, lattice_viscosity_from_relaxation_rate(omega_s)),
+           *[Assignment(matrix_entries[i], pi[i]) for i in range(stencil.D ** 2)],
+           Assignment(r_prime, sp.Float(-1) ** (stencil.D + 1) * matrix_entries.det()),
+           Assignment(q_prime, sp.Rational(1, 2) * (matrix_entries * matrix_entries).trace()),
+           Assignment(c_pi_s, c_pi),
+           *nu_e_assignments,
+           Assignment(adapted_omega, relaxation_rate_from_lattice_viscosity(nu_0 + nu_e))]
+    collision_rule.subexpressions += eqs
+    collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
+
+    if eddy_viscosity_field:
+        collision_rule.main_assignments.append(Assignment(eddy_viscosity_field.center, nu_e))
+
+    if omega_output_field:
+        collision_rule.main_assignments.append(Assignment(omega_output_field.center, adapted_omega))
+
+    return collision_rule

lbmpy/updatekernels.py → src/lbmpy/updatekernels.py

-import numpy as np
 import sympy as sp
 
 from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor
 from lbmpy.methods.abstractlbmethod import LbmCollisionRule
-from pystencils import Assignment, AssignmentCollection, Field
+from pystencils import Assignment, AssignmentCollection
 from pystencils.field import create_numpy_array_with_layout, layout_string_to_tuple
 from pystencils.simp import add_subexpressions_for_field_reads
 from pystencils.sympyextensions import fast_subs
@@ -11,30 +10,34 @@ from pystencils.sympyextensions import fast_subs
 # -------------------------------------------- LBM Kernel Creation -----------------------------------------------------
 
 
-def create_lbm_kernel(collision_rule, input_field, output_field, accessor):
+def create_lbm_kernel(collision_rule, src_field, dst_field=None, accessor=StreamPullTwoFieldsAccessor(),
+                      data_type=None):
     """Replaces the pre- and post collision symbols in the collision rule by field accesses.
 
     Args:
         collision_rule:  instance of LbmCollisionRule, defining the collision step
-        input_field: field used for reading pdf values
-        output_field: field used for writing pdf values
-                      if accessor.is_inplace this parameter is ignored
+        src_field: field used for reading pdf values
+        dst_field: field used for writing pdf values if accessor.is_inplace this parameter is ignored
         accessor: instance of PdfFieldAccessor, defining where to read and write values
-                  to create e.g. a fused stream-collide kernel
+                  to create e.g. a fused stream-collide kernel See 'fieldaccess.PdfFieldAccessor'
+        data_type: If a datatype is given the field reads are isolated and written to TypedSymbols of type data_type
 
     Returns:
         LbmCollisionRule where pre- and post collision symbols have been replaced
     """
     if accessor.is_inplace:
-        output_field = input_field
+        dst_field = src_field
+
+    if not accessor.is_inplace and dst_field is None:
+        raise ValueError("For two field accessors a destination field has to be provided")
 
     method = collision_rule.method
     pre_collision_symbols = method.pre_collision_pdf_symbols
     post_collision_symbols = method.post_collision_pdf_symbols
     substitutions = {}
 
-    input_accesses = accessor.read(input_field, method.stencil)
-    output_accesses = accessor.write(output_field, method.stencil)
+    input_accesses = accessor.read(src_field, method.stencil)
+    output_accesses = accessor.write(dst_field, method.stencil)
 
     for (idx, offset), input_access, output_access in zip(enumerate(method.stencil), input_accesses, output_accesses):
         substitutions[pre_collision_symbols[idx]] = input_access
@@ -48,41 +51,67 @@ def create_lbm_kernel(collision_rule, input_field, output_field, accessor):
             new_split_groups.append([fast_subs(e, substitutions) for e in split_group])
         result.simplification_hints['split_groups'] = new_split_groups
 
-    if accessor.is_inplace:
-        result = add_subexpressions_for_field_reads(result, subexpressions=True, main_assignments=True)
+    if accessor.is_inplace or (data_type is not None and src_field.dtype != data_type):
+        result = add_subexpressions_for_field_reads(result, subexpressions=True, main_assignments=True,
+                                                    data_type=data_type)
 
     return result
 
 
-def create_stream_pull_only_kernel(stencil, numpy_arr=None, src_field_name="src", dst_field_name="dst",
-                                   generic_layout='numpy', generic_field_type=np.float64):
-    """Creates a stream-pull kernel, without collision.
+def create_stream_only_kernel(stencil, src_field, dst_field=None, accessor=StreamPullTwoFieldsAccessor()):
+    """Creates a stream kernel, without collision.
 
-    For parameters see function ``create_stream_pull_collide_kernel``
+    Args:
+        stencil: lattice Boltzmann stencil which is used in the form of a tuple of tuples
+        src_field: field used for reading pdf values
+        dst_field: field used for writing pdf values if accessor.is_inplace this parameter is not necessary but it
+                   is used if provided
+        accessor: instance of PdfFieldAccessor, defining where to read and write values
+                  to create e.g. a fused stream-collide kernel See 'fieldaccess.PdfFieldAccessor'
+
+    Returns:
+        AssignmentCollection of the stream only update rule
     """
+    if accessor.is_inplace and dst_field is None:
+        dst_field = src_field
+
+    if not accessor.is_inplace and dst_field is None:
+        raise ValueError("For two field accessors a destination field has to be provided")
 
-    dim = len(stencil[0])
-    if numpy_arr is None:
-        src = Field.create_generic(src_field_name, dim, index_shape=(len(stencil),),
-                                   layout=generic_layout, dtype=generic_field_type)
-        dst = Field.create_generic(dst_field_name, dim, index_shape=(len(stencil),),
-                                   layout=generic_layout, dtype=generic_field_type)
-    else:
-        src = Field.create_from_numpy_array(src_field_name, numpy_arr, index_dimensions=1)
-        dst = Field.create_from_numpy_array(dst_field_name, numpy_arr, index_dimensions=1)
+    temporary_symbols = sp.symbols(f'streamed_:{stencil.Q}')
+    subexpressions = [Assignment(tmp, acc) for tmp, acc in zip(temporary_symbols, accessor.read(src_field, stencil))]
+    main_assignments = [Assignment(acc, tmp) for acc, tmp in zip(accessor.write(dst_field, stencil), temporary_symbols)]
+    return AssignmentCollection(main_assignments, subexpressions=subexpressions)
 
-    accessor = StreamPullTwoFieldsAccessor()
-    eqs = [Assignment(a, b) for a, b in zip(accessor.write(dst, stencil), accessor.read(src, stencil))]
-    return AssignmentCollection(eqs, [])
 
+def create_stream_pull_with_output_kernel(lb_method, src_field, dst_field=None, output=None,
+                                          accessor=StreamPullTwoFieldsAccessor()):
+    """Creates a stream kernel, without collision but macroscopic quantaties like density or velocity can be calculated.
+
+    Args:
+        lb_method: lattice Boltzmann method see 'creationfunctions.create_lb_method'
+        src_field: field used for reading pdf values
+        dst_field: field used for writing pdf values if accessor.is_inplace this parameter is ignored
+        output: dictonary which containes macroscopic quantities as keys which should be calculated and fields as
+                values which should be used to write the data e.g.: {'density': density_field}
+        accessor: instance of PdfFieldAccessor, defining where to read and write values
+                  to create e.g. a fused stream-collide kernel See 'fieldaccess.PdfFieldAccessor'
+
+    Returns:
+        AssignmentCollection of the stream only update rule
+    """
+    if accessor.is_inplace:
+        dst_field = src_field
+
+    if not accessor.is_inplace and dst_field is None:
+        raise ValueError("For two field accessors a destination field has to be provided")
 
-def create_stream_pull_with_output_kernel(lb_method, src_field, dst_field, output):
     stencil = lb_method.stencil
     cqc = lb_method.conserved_quantity_computation
-    streamed = sp.symbols("streamed_:%d" % (len(stencil),))
-    accessor = StreamPullTwoFieldsAccessor()
+    streamed = sp.symbols(f"streamed_:{stencil.Q}")
     stream_assignments = [Assignment(a, b) for a, b in zip(streamed, accessor.read(src_field, stencil))]
-    output_eq_collection = cqc.output_equations_from_pdfs(streamed, output)
+    output_eq_collection = cqc.output_equations_from_pdfs(streamed, output) if output\
+        else AssignmentCollection(main_assignments=[])
     write_eqs = [Assignment(a, b) for a, b in zip(accessor.write(dst_field, stencil), streamed)]
 
     subexpressions = stream_assignments + output_eq_collection.subexpressions
@@ -90,6 +119,27 @@ def create_stream_pull_with_output_kernel(lb_method, src_field, dst_field, outpu
     return LbmCollisionRule(lb_method, main_eqs, subexpressions,
                             simplification_hints=output_eq_collection.simplification_hints)
 
+
+def create_copy_kernel(stencil, src_field, dst_field, accessor=StreamPullTwoFieldsAccessor()):
+    """Creates a copy kernel, which can be used to transfer information from src to dst field.
+
+    Args:
+        stencil: lattice Boltzmann stencil which is used in the form of a tuple of tuples
+        src_field: field used for reading pdf values
+        dst_field: field used for writing pdf values
+        accessor: instance of PdfFieldAccessor, defining where to read and write values
+                  to create e.g. a fused stream-collide kernel See 'fieldaccess.PdfFieldAccessor'
+
+    Returns:
+        AssignmentCollection of a copy update rule
+    """
+
+    temporary_symbols = sp.symbols(f'copied:{stencil.Q}')
+    subexpressions = [Assignment(tmp, acc) for tmp, acc in zip(temporary_symbols, accessor.write(src_field, stencil))]
+    main_assignments = [Assignment(acc, tmp) for acc, tmp in zip(accessor.write(dst_field, stencil), temporary_symbols)]
+    return AssignmentCollection(main_assignments, subexpressions=subexpressions)
+
+
 # ---------------------------------- Pdf array creation for various layouts --------------------------------------------
 
 

src/lbmpy/utils.py

+import sympy as sp
+
+
+def second_order_moment_tensor(function_values, stencil):
+    """Returns (D x D) Matrix of second order moments of the given function where D is the dimension"""
+    assert len(function_values) == stencil.Q
+    return sp.Matrix(stencil.D, stencil.D, lambda i, j: sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)))
+
+
+def frobenius_norm(matrix, factor=1):
+    """Computes the Frobenius norm of a matrix defined as the square root of the sum of squared matrix elements
+    The optional factor is added inside the square root"""
+    return sp.sqrt(sum(i * i for i in matrix) * factor)
+
+
+def extract_shear_relaxation_rate(collision_rule, shear_relaxation_rate):
+    """Searches for the shear relaxation rate in the collision equations.
+       If the shear relaxation rate is assigned to a single assignment its lhs is returned.
+       Otherwise, the shear relaxation rate has to be a sympy symbol, or it can not be extracted"""
+    found_symbolic_shear_relaxation = True
+    if isinstance(shear_relaxation_rate, (float, int)):
+        found_symbolic_shear_relaxation = False
+        for eq in collision_rule.all_assignments:
+            if eq.rhs == shear_relaxation_rate:
+                found_symbolic_shear_relaxation = True
+                shear_relaxation_rate = eq.lhs
+
+    return shear_relaxation_rate, found_symbolic_shear_relaxation

tests/advanced_streaming/test_advanced_streaming_noslip.py

+import numpy as np
+
+from lbmpy.advanced_streaming import Timestep
+from lbmpy.boundaries import NoSlip
+from lbmpy.boundaries.boundaryhandling import create_lattice_boltzmann_boundary_kernel
+from lbmpy.advanced_streaming.utility import streaming_patterns, inverse_dir_index, AccessPdfValues
+from lbmpy.enums import Method, Stencil
+from lbmpy.creationfunctions import create_lb_method, LBMConfig
+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.field import Field, FieldType
+
+import pytest
+
+
+@pytest.mark.parametrize("stencil", [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.parametrize("streaming_pattern", streaming_patterns)
+@pytest.mark.parametrize("prev_timestep", [Timestep.EVEN, Timestep.ODD])
+def test_advanced_streaming_noslip_single_cell(stencil, streaming_pattern, prev_timestep):
+    """
+    Advanced Streaming NoSlip Test
+    """
+
+    stencil = LBStencil(stencil)
+    pdf_field = ps.fields(f'pdfs({stencil.Q}): [{stencil.D}D]')
+
+    prev_pdf_access = AccessPdfValues(stencil, streaming_pattern, prev_timestep, 'out')
+    next_pdf_access = AccessPdfValues(stencil, streaming_pattern, prev_timestep.next(), 'in')
+
+    pdfs = np.zeros((3,) * stencil.D + (stencil.Q,))
+    pos = (1,) * stencil.D
+    for d in range(stencil.Q):
+        prev_pdf_access.write_pdf(pdfs, pos, d, d)
+
+    lbm_config = LBMConfig(stencil=stencil, method=Method.SRT)
+    lb_method = create_lb_method(lbm_config=lbm_config)
+    noslip = NoSlip()
+
+    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))
+    index_vector = np.array([pos + (d,) for d in range(stencil.Q)], dtype=index_struct_dtype)
+
+    ast = create_lattice_boltzmann_boundary_kernel(pdf_field,
+                                                   index_field, lb_method, noslip,
+                                                   prev_timestep=prev_timestep,
+                                                   streaming_pattern=streaming_pattern)
+
+    flex_kernel = ast.compile()
+
+    flex_kernel(pdfs=pdfs, indexVector=index_vector, indexVectorSize=len(index_vector))
+
+    reflected_pdfs = [next_pdf_access.read_pdf(pdfs, pos, d) for d in range(stencil.Q)]
+    inverse_pdfs = [inverse_dir_index(stencil, d) for d in range(stencil.Q)]
+    assert reflected_pdfs == inverse_pdfs

tests/advanced_streaming/test_communication.py

+import pystencils as ps
+
+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,
+    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])
+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,
+    )
+    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])
+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,
+    )
+    for _, slices_list in slices.items():
+        for src, dst in slices_list:
+            src_shape = arr[src].shape
+            dst_shape = arr[dst].shape
+            assert src_shape == dst_shape
+
+
+@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
+    )
+    for i, d in enumerate(stencil):
+        if i == 0:
+            continue
+
+        for s in slices[d]:
+            if s[0][-1] == i:
+                src = s[0][:-1]
+                dst = s[1][:-1]
+                break
+
+        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)
+        )
+        assert src == inner_slice
+        assert dst == gl_slice
+
+
+@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
+
+    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)
+    pdf_tmp = dh.add_array("pdf_tmp", values_per_cell=len(stencil), dtype=np.int64)
+    dh.fill("pdf_tmp", 0, ghost_layers=True)
+
+    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
+        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)
+
+    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.run_kernel(stream)
+    dh.swap("pdf", "pdf_tmp")
+    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
+        num += 1
+
+    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")
+
+    if stencil.D == 3:
+        for i in range(gl, domain_size[0]):
+            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)
+    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]
+                    )

tests/advanced_streaming/test_fully_periodic_flow.py

+import numpy as np
+from dataclasses import replace
+
+from pystencils.datahandling import create_data_handling
+from pystencils import create_kernel, Target, CreateKernelConfig
+
+from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_function, LBMConfig, LBMOptimisation
+from lbmpy.enums import Method, Stencil
+from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
+from lbmpy.stencils import LBStencil
+
+from lbmpy.advanced_streaming import LBMPeriodicityHandling
+from lbmpy.advanced_streaming.utility import is_inplace, streaming_patterns, get_timesteps
+
+import pytest
+from numpy.testing import assert_allclose
+
+all_results = dict()
+
+targets = [Target.CPU]
+
+try:
+    import cupy
+    targets += [Target.GPU]
+except Exception:
+    pass
+
+
+@pytest.mark.parametrize('target', targets)
+@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
+@pytest.mark.longrun
+def test_fully_periodic_flow(target, stencil, streaming_pattern):
+    gpu = False
+    if target == Target.GPU:
+        gpu = True
+
+    #   Stencil
+    stencil = LBStencil(stencil)
+
+    #   Streaming
+    inplace = is_inplace(streaming_pattern)
+    timesteps = get_timesteps(streaming_pattern)
+    zeroth_timestep = timesteps[0]
+
+    #   Data Handling and PDF fields
+    domain_size = (30,) * stencil.D
+    periodicity = (True,) * stencil.D
+
+    dh = create_data_handling(domain_size=domain_size, periodicity=periodicity, default_target=target)
+
+    pdfs = dh.add_array('pdfs', stencil.Q)
+    if not inplace:
+        pdfs_tmp = dh.add_array_like('pdfs_tmp', pdfs.name)
+
+    #   LBM Streaming and Collision
+    lbm_config = LBMConfig(stencil=stencil, method=Method.SRT,
+                           relaxation_rate=1.0, streaming_pattern=streaming_pattern)
+
+    lbm_opt = LBMOptimisation(symbolic_field=pdfs)
+    config = CreateKernelConfig(target=target)
+
+    if not inplace:
+        lbm_opt = replace(lbm_opt, symbolic_temporary_field=pdfs_tmp)
+
+    lb_collision = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt, config=config)
+    lb_method = lb_collision.method
+
+    lb_kernels = []
+    for t in timesteps:
+        lbm_config = replace(lbm_config, timestep=t)
+        lb_kernels.append(create_lb_function(collision_rule=lb_collision,
+                                             lbm_config=lbm_config, lbm_optimisation=lbm_opt))
+
+    #   Macroscopic Values
+    density = 1.0
+    density_field = dh.add_array('rho', 1)
+    u_x = 0.01
+    velocity = (u_x,) * stencil.D
+    velocity_field = dh.add_array('u', stencil.D)
+
+    u_ref = np.full(domain_size + (stencil.D,), u_x)
+
+    setter = macroscopic_values_setter(
+        lb_method, density, velocity, pdfs,
+        streaming_pattern=streaming_pattern, previous_timestep=zeroth_timestep)
+    setter_kernel = create_kernel(setter,
+                                  config=CreateKernelConfig(target=target, ghost_layers=1)).compile()
+
+    getter_kernels = []
+    for t in timesteps:
+        getter = macroscopic_values_getter(
+            lb_method, density_field, velocity_field, pdfs,
+            streaming_pattern=streaming_pattern, previous_timestep=t)
+        getter_kernels.append(create_kernel(getter,
+                                            config=CreateKernelConfig(target=target, ghost_layers=1)).compile())
+
+    #   Periodicity
+    periodicity_handler = LBMPeriodicityHandling(stencil, dh, pdfs.name, streaming_pattern=streaming_pattern)
+
+    # Initialization and Timestep
+    current_timestep = zeroth_timestep
+
+    def init():
+        global current_timestep
+        current_timestep = zeroth_timestep
+        dh.run_kernel(setter_kernel)
+
+    def one_step():
+        global current_timestep
+
+        # Periodicty
+        periodicity_handler(current_timestep)
+
+        # Here, the next time step begins
+        current_timestep = current_timestep.next()
+
+        # LBM Step
+        dh.run_kernel(lb_kernels[current_timestep.idx])
+
+        # Field Swaps
+        if not inplace:
+            dh.swap(pdfs.name, pdfs_tmp.name)
+
+        # Macroscopic Values
+        dh.run_kernel(getter_kernels[current_timestep.idx])
+
+    #   Run the simulation
+    init()
+
+    for _ in range(100):
+        one_step()
+
+    #   Evaluation
+    if gpu:
+        dh.to_cpu(velocity_field.name)
+    u = dh.gather_array(velocity_field.name)
+
+    #   Equal to the steady-state velocity field up to numerical errors
+    assert_allclose(u, u_ref)
+
+    #   Flow must be equal up to numerical error for all streaming patterns
+    global all_results
+    for key, prev_u in all_results.items():
+        if key[0] == stencil:
+            prev_pattern = key[1]
+            assert_allclose(
+                u, prev_u, err_msg=f'Velocity field for {streaming_pattern} differed from {prev_pattern}!')
+    all_results[(stencil, streaming_pattern)] = u

tests/advanced_streaming/test_periodic_pipe_with_force.py

+import numpy as np
+from dataclasses import replace
+
+from pystencils.datahandling import create_data_handling
+from pystencils import create_kernel, Target, CreateKernelConfig
+
+from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_function, LBMConfig, LBMOptimisation
+from lbmpy.enums import ForceModel, Method, Stencil
+from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
+from lbmpy.stencils import LBStencil
+
+from lbmpy.advanced_streaming import LBMPeriodicityHandling
+from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling
+from lbmpy.advanced_streaming.utility import is_inplace, streaming_patterns, get_timesteps
+
+import pytest
+from numpy.testing import assert_allclose
+
+all_results = dict()
+
+targets = [Target.CPU]
+
+try:
+    import cupy
+    targets += [Target.GPU]
+except Exception:
+    pass
+
+
+class PeriodicPipeFlow:
+    def __init__(self, stencil, streaming_pattern, wall_boundary=None, target=Target.CPU):
+
+        if wall_boundary is None:
+            wall_boundary = NoSlip()
+
+        self.target = target
+        self.gpu = target in [Target.GPU]
+
+        #   Stencil
+        self.stencil = stencil
+        self.q = stencil.Q
+        self.dim = stencil.D
+
+        #   Streaming
+        self.streaming_pattern = streaming_pattern
+        self.inplace = is_inplace(self.streaming_pattern)
+        self.timesteps = get_timesteps(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)
+        self.force = (0.0001, ) + (0.0,) * (self.dim - 1)
+
+        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 Streaming and Collision
+        lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.0,
+                               force_model=ForceModel.GUO, force=self.force, streaming_pattern=streaming_pattern)
+
+        lbm_opt = LBMOptimisation(symbolic_field=self.pdfs)
+        config = CreateKernelConfig(target=self.target)
+
+        if not self.inplace:
+            lbm_opt = replace(lbm_opt, symbolic_temporary_field=self.pdfs_tmp)
+
+        self.lb_collision = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+        self.lb_method = self.lb_collision.method
+
+        self.lb_kernels = []
+        for t in self.timesteps:
+            lbm_config = replace(lbm_config, timestep=t)
+            self.lb_kernels.append(create_lb_function(collision_rule=self.lb_collision,
+                                                      lbm_config=lbm_config,
+                                                      lbm_optimisation=lbm_opt,
+                                                      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)
+
+        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=CreateKernelConfig(target=target, ghost_layers=1)).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=CreateKernelConfig(target=target, ghost_layers=1)).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
+
+
+@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
+@pytest.mark.parametrize('target', targets)
+@pytest.mark.longrun
+def test_periodic_pipe(stencil, streaming_pattern, target):
+    stencil = LBStencil(stencil)
+    pipeflow = PeriodicPipeFlow(stencil, streaming_pattern, target=target)
+    pipeflow.init()
+    pipeflow.run(100)
+    u = pipeflow.get_trimmed_velocity_array()
+
+    #   Flow must be equal up to numerical error for all streaming patterns
+    global all_results
+    for key, prev_u in all_results.items():
+        if key[0] == stencil:
+            prev_pattern = key[1]
+            assert_allclose(
+                u, prev_u,
+                rtol=1, atol=1e-16,
+                err_msg=f'Velocity field for {streaming_pattern} differed from {prev_pattern}!')
+    all_results[(stencil, streaming_pattern)] = u

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

 from statistics import median
 
 import numpy as np
-import pytest
 import sympy as sp
 
 from lbmpy.scenarios import create_lid_driven_cavity
@@ -86,14 +85,16 @@ def method_options(dim=2, with_srt=False, with_mrt=False,
     rr_free = "rr_free"
     methods = [{'method': 'trt'}]
     if with_mrt:
-        methods += [{'method': 'mrt3'}, {'method': 'mrt_raw'}]
+        methods += [{'method': 'mrt'}, {'method': 'mrt_raw'}]
 
     if with_srt:
         methods += [{'method': 'srt'}]
 
     if with_entropic:
-        methods += [{'entropic': True, 'method': 'mrt3', 'relaxation_rates': [1.5, 1.5, rr_free]},
-                    {'entropic': True, 'method': 'mrt3', 'relaxation_rates': [1.5, rr_free, rr_free]},
+        methods += [{'entropic': True, 'method': 'mrt', 'relaxation_rates': [1.5, 1.5, rr_free, rr_free,
+                                                                             rr_free, rr_free]},
+                    {'entropic': True, 'method': 'mrt', 'relaxation_rates': [1.5, rr_free, rr_free, rr_free,
+                                                                             rr_free, rr_free]},
                     {'entropic': True, 'method': 'trt-kbc-n1'},
                     {'entropic': True, 'method': 'trt-kbc-n2'},
                     {'entropic': True, 'method': 'trt-kbc-n3'},
@@ -226,11 +227,11 @@ def study_block_sizes_trt(study):
                 study.add_run(params, weight=ds[0] // domain_sizes[0][0])
 
 
-@pytest.mark.longrun
-def test_run():
-    """Called by test suite - ensures that benchmark can be run"""
-    s = ParameterStudy(run)
-    study_optimization_options(s, domain_sizes=((17, 23), (19, 17, 18))).run()
+# @pytest.mark.longrun
+# def test_run():
+#     """Called by test suite - ensures that benchmark can be run"""
+#     s = ParameterStudy(run)
+#     study_optimization_options(s, domain_sizes=((17, 23), (19, 17, 18))).run()
 
 
 def study_compiler_flags(study):

tests/cumulantmethod/test_equilibrium.py

+from lbmpy.creationfunctions import LBMConfig
+from lbmpy.enums import Method, Stencil, CollisionSpace
+from lbmpy.maxwellian_equilibrium import generate_equilibrium_by_matching_moments
+from lbmpy.moments import extract_monomials
+from lbmpy.stencils import LBStencil
+import pytest
+import sympy as sp
+
+from pystencils.simp import AssignmentCollection
+from pystencils import Assignment
+from lbmpy.creationfunctions import create_lb_method
+from lbmpy.methods import CollisionSpaceInfo
+from lbmpy.moment_transforms import (
+    FastCentralMomentTransform,
+    PdfsToCentralMomentsByMatrix,
+    PdfsToCentralMomentsByShiftMatrix,
+    BinomialChimeraTransform)
+
+sympy_numeric_version = [int(x, 10) for x in sp.__version__.split('.')]
+if len(sympy_numeric_version) < 3:
+    sympy_numeric_version.append(0)
+sympy_numeric_version.reverse()
+sympy_version = sum(x * (100 ** i) for i, x in enumerate(sympy_numeric_version))
+
+
+reference_equilibria = dict()    
+
+
+@pytest.mark.skipif(sympy_version < 10200,
+                    reason="Old Sympy Versions take too long to form the inverse moment matrix")
+@pytest.mark.parametrize('stencil_name', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.parametrize('cm_transform', [PdfsToCentralMomentsByMatrix,
+                                          FastCentralMomentTransform,
+                                          PdfsToCentralMomentsByShiftMatrix,
+                                          BinomialChimeraTransform])
+def test_equilibrium_pdfs(stencil_name, cm_transform):
+    stencil = LBStencil(stencil_name)
+    cspace = CollisionSpaceInfo(CollisionSpace.CUMULANTS, central_moment_transform_class=cm_transform)
+    lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, compressible=True, zero_centered=False,
+                           collision_space_info=cspace)
+
+    c_lb_method = create_lb_method(lbm_config=lbm_config)
+    rho = c_lb_method.zeroth_order_equilibrium_moment_symbol
+    u = c_lb_method.first_order_equilibrium_moment_symbols
+    pdfs = c_lb_method.post_collision_pdf_symbols
+    cqe_assignments = [Assignment(sym, sym) for sym in u] + [Assignment(rho, rho)]
+    cqe = AssignmentCollection(main_assignments=cqe_assignments)
+    method_equilibrium_eqs = c_lb_method.get_equilibrium(cqe, False, False).main_assignments_dict
+
+    #   Reference Equations
+    ref_equilibrium = reference_equilibria.get(stencil_name, None)
+    if ref_equilibrium is None:
+        raw_moments = list(extract_monomials(c_lb_method.cumulants, dim=stencil.D))
+        ref_equilibrium = generate_equilibrium_by_matching_moments(
+            stencil, tuple(raw_moments), rho=rho, u=u, c_s_sq=sp.Rational(1, 3), order=2*stencil.D)
+        reference_equilibria[stencil_name] = ref_equilibrium
+
+    for i in range(stencil.Q):
+        method_eq = method_equilibrium_eqs[pdfs[i]]
+        ref_eq = ref_equilibrium[i]
+        assert method_eq.expand() - ref_eq.expand() == sp.sympify(0), \
+            f"Equilibrium equation for pdf index {i} did not match." 

tests/cumulantmethod/test_flow_around_sphere.py

+import pytest
+import numpy as np
+from dataclasses import replace
+
+from lbmpy.advanced_streaming import BetweenTimestepsIndexing, Timestep, get_timesteps
+from lbmpy.creationfunctions import create_lb_function, create_lb_collision_rule,\
+    create_lb_method, LBMConfig, LBMOptimisation
+from lbmpy.boundaries import LatticeBoltzmannBoundaryHandling, NoSlip, UBB
+from lbmpy.enums import Method, Stencil
+
+from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
+from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
+from lbmpy.stencils import LBStencil
+
+from pystencils import create_kernel, create_data_handling, Assignment, Target, CreateKernelConfig
+from pystencils.slicing import slice_from_direction, get_slice_before_ghost_layer
+
+
+def flow_around_sphere(stencil, galilean_correction, fourth_order_correction, L_LU, total_steps):
+    if galilean_correction and stencil.Q != 27:
+        pytest.skip()
+
+    target = Target.GPU
+    streaming_pattern = 'aa'
+    timesteps = get_timesteps(streaming_pattern)
+
+    u_max = 0.05
+    Re = 500000
+
+    kinematic_viscosity = (L_LU * u_max) / Re
+    initial_velocity = (u_max,) + (0,) * (stencil.D - 1)
+
+    omega_v = relaxation_rate_from_lattice_viscosity(kinematic_viscosity)
+
+    channel_size = (10 * L_LU,) + (5 * L_LU,) * (stencil.D - 1)
+    sphere_position = (channel_size[0] // 3,) + (channel_size[1] // 2,) * (stencil.D - 1)
+    sphere_radius = L_LU // 2
+
+    lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, relaxation_rate=omega_v,
+                           compressible=True,
+                           galilean_correction=galilean_correction, fourth_order_correction=fourth_order_correction)
+
+    lbm_opt = LBMOptimisation(pre_simplification=True)
+    config = CreateKernelConfig(target=target)
+
+    lb_method = create_lb_method(lbm_config=lbm_config)
+
+    def get_extrapolation_kernel(timestep):
+        boundary_assignments = []
+        indexing = BetweenTimestepsIndexing(
+            pdf_field, stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep)
+        f_out, _ = indexing.proxy_fields
+        for i, d in enumerate(stencil):
+            if d[0] == -1:
+                asm = Assignment(f_out.neighbor(0, 1)(i), f_out.center(i))
+                boundary_assignments.append(asm)
+        boundary_assignments = indexing.substitute_proxies(
+            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,
+                skip_independence_check=True)
+            )
+        return extrapolation_ast.compile()
+
+    dh = create_data_handling(channel_size, periodicity=False, default_layout='fzyx', default_target=target)
+
+    u_field = dh.add_array('u', stencil.D)
+    rho_field = dh.add_array('rho', 1)
+    pdf_field = dh.add_array('pdfs', stencil.Q)
+
+    dh.fill(u_field.name, 0.0, ghost_layers=True)
+    dh.fill(rho_field.name, 0.0, ghost_layers=True)
+
+    dh.to_gpu(u_field.name)
+    dh.to_gpu(rho_field.name)
+
+    lbm_opt = replace(lbm_opt, symbolic_field=pdf_field)
+
+    bh = LatticeBoltzmannBoundaryHandling(lb_method, dh, pdf_field.name,
+                                          streaming_pattern=streaming_pattern, target=target)
+    wall = NoSlip()
+    inflow = UBB(initial_velocity)
+
+    bh.set_boundary(inflow, slice_from_direction('W', stencil.D))
+
+    directions = ('N', 'S', 'T', 'B') if stencil.D == 3 else ('N', 'S')
+    for direction in directions:
+        bh.set_boundary(wall, slice_from_direction(direction, stencil.D))
+
+    outflow_kernels = [get_extrapolation_kernel(Timestep.EVEN), get_extrapolation_kernel(Timestep.ODD)]
+
+    def sphere_boundary_callback(x, y, z=None):
+        x = x - sphere_position[0]
+        y = y - sphere_position[1]
+        z = z - sphere_position[2] if z is not None else 0
+        return np.sqrt(x ** 2 + y ** 2 + z ** 2) <= sphere_radius
+
+    bh.set_boundary(wall, mask_callback=sphere_boundary_callback)
+
+    init_eqs = pdf_initialization_assignments(lb_method, 1.0, initial_velocity, pdf_field,
+                                              streaming_pattern=streaming_pattern,
+                                              previous_timestep=timesteps[0])
+    init_kernel = create_kernel(init_eqs, config=config).compile()
+
+    output = {
+        'density': rho_field,
+        'velocity': u_field
+    }
+    lbm_config = replace(lbm_config, output=output)
+
+    lb_collision_rule = create_lb_collision_rule(lb_method=lb_method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+
+    lb_kernels = []
+    for t in timesteps:
+        lbm_config = replace(lbm_config, timestep=t)
+        lbm_config = replace(lbm_config, streaming_pattern=streaming_pattern)
+        lb_kernels.append(create_lb_function(collision_rule=lb_collision_rule,
+                                             lbm_config=lbm_config, lbm_optimisation=lbm_opt, config=config))
+
+    timestep = timesteps[0]
+
+    dh.run_kernel(init_kernel)
+
+    stability_check_frequency = 1000
+
+    for i in range(total_steps):
+        bh(prev_timestep=timestep)
+        dh.run_kernel(outflow_kernels[timestep.idx])
+        timestep = timestep.next()
+        dh.run_kernel(lb_kernels[timestep.idx])
+
+        if i % stability_check_frequency == 0:
+            dh.to_cpu(u_field.name)
+            assert np.isfinite(dh.cpu_arrays[u_field.name]).all()
+
+
+@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.parametrize('galilean_correction', [False, True])
+@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, 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/kida_vortex_flow/scenario_kida_vortex_flow.py → tests/full_scenarios/kida_vortex_flow/scenario_kida_vortex_flow.py

@@ -250,8 +250,8 @@ def run(re=6000, eval_interval=0.05, total_time=3.0, domain_size=100, u_0=0.05,
 
 def create_full_parameter_study(gpu=False):
     """Creates a parameter study that can run the Kida vortex flow with entropic, KBC, Smagorinsky and MRT methods."""
-    opt_cpu = {'target': 'cpu', 'openmp': 4}
-    opt_gpu = {'target': 'gpu'}
+    opt_cpu = {'target': ps.Target.CPU, 'openmp': 4}
+    opt_gpu = {'target': ps.Target.GPU}
 
     mrt_one = [{'method': 'mrt3', 'relaxation_rates': ['viscosity', 1, 1], 'stencil': stencil}
                for stencil in ('D3Q19', 'D3Q27')]