lbmpy_tests/test_cpu_gpu_equivalence.py → tests/test_cpu_gpu_equivalence.py

@@ -5,8 +5,9 @@ import pytest
 from pystencils import Backend, Target, CreateKernelConfig
 
 from lbmpy.creationfunctions import LBMConfig, LBMOptimisation
+from lbmpy.enums import Method, Stencil
 from lbmpy.scenarios import create_channel
-from lbmpy.enums import Method
+from lbmpy.stencils import LBStencil
 
 
 def run_equivalence_test(domain_size, lbm_config, lbm_opt, config, time_steps=13):
@@ -31,14 +32,16 @@ def run_equivalence_test(domain_size, lbm_config, lbm_opt, config, time_steps=13
                                       ((18, 20), Method.MRT, True, (4, 2), 'zyxf'),
                                       ((7, 11, 18), Method.TRT, False, False, 'numpy')])
 def test_force_driven_channel_short(scenario):
-    pytest.importorskip("pycuda")
+    pytest.importorskip("cupy")
     ds = scenario[0]
+    stencil = LBStencil(Stencil.D2Q9) if len(ds) == 2 else LBStencil(Stencil.D3Q27)
     method = scenario[1]
     compressible = scenario[2]
     block_size = scenario[3]
     field_layout = scenario[4]
 
-    lbm_config = LBMConfig(method=method, compressible=compressible, relaxation_rates=[1.95, 1.9, 1.92, 1.92])
+    lbm_config = LBMConfig(stencil=stencil, method=method,
+                           compressible=compressible, relaxation_rates=[1.95, 1.9, 1.92, 1.92])
     lbm_opt = LBMOptimisation(field_layout=field_layout)
 
     # Different methods

lbmpy_tests/test_cumulant_transform.py → tests/test_cumulant_transform.py

@@ -6,11 +6,12 @@ from pystencils import AssignmentCollection
 
 from lbmpy.moment_transforms import (
     CentralMomentsToCumulantsByGeneratingFunc, PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT,
-    PRE_COLLISION_CUMULANT, PRE_COLLISION_MONOMIAL_CUMULANT
+    PRE_COLLISION_CUMULANT
 )
 from lbmpy.methods import cascaded_moment_sets_literature
 from lbmpy.stencils import Stencil, LBStencil
 
+
 @pytest.mark.parametrize("stencil", [Stencil.D2Q9, Stencil.D3Q19])
 def test_identity(stencil):
     stencil = LBStencil(stencil)
@@ -18,7 +19,7 @@ def test_identity(stencil):
     rho = sp.Symbol('rho')
     u = sp.symbols('u_:2')
     transform = CentralMomentsToCumulantsByGeneratingFunc(stencil, polys, rho, u,
-                    post_collision_symbol_base=PRE_COLLISION_CUMULANT)
+                                                          post_collision_symbol_base=PRE_COLLISION_CUMULANT)
 
     forward_eqs = transform.forward_transform()
     backward_eqs = transform.backward_transform(central_moment_base=PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT)

lbmpy_tests/test_diffusion.py → tests/test_diffusion.py

@@ -77,7 +77,7 @@ def test_diffusion():
 
       The hydrodynamic field is not simulated, instead a constant velocity is assumed.
     """
-    pytest.importorskip("pycuda")
+    pytest.importorskip("cupy")
     # Parameters
     domain_size = (1600, 160)
     omega = 1.38

lbmpy_tests/test_float_kernel.py → tests/test_float_kernel.py

@@ -3,8 +3,9 @@ import pytest
 import pystencils as ps
 
 from lbmpy.creationfunctions import create_lb_function, LBMConfig
-from lbmpy.enums import Method
+from lbmpy.enums import Method, Stencil
 from lbmpy.scenarios import create_lid_driven_cavity
+from lbmpy.stencils import LBStencil
 
 
 @pytest.mark.parametrize('double_precision', [False, True])
@@ -28,7 +29,8 @@ def test_creation(double_precision, method_enum):
 @pytest.mark.parametrize('double_precision', [False, True])
 @pytest.mark.parametrize('method_enum', [Method.SRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
 def test_scenario(method_enum, double_precision):
-    lbm_config = LBMConfig(method=method_enum, relaxation_rate=1.45, compressible=True)
+    lbm_config = LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=method_enum,
+                           relaxation_rate=1.45, compressible=True)
     config = ps.CreateKernelConfig(data_type="float64" if double_precision else "float32",
                                    default_number_float="float64" if double_precision else "float32")
     sc = create_lid_driven_cavity((16, 16, 8), lbm_config=lbm_config, config=config)

lbmpy_tests/test_fluctuating_lb.py → tests/test_fluctuating_lb.py

 """Tests velocity and stress fluctuations for thermalized LB"""
 
+import warnings
+import pytest
 import pystencils as ps
+
+from pystencils import get_code_str
+from pystencils.backends.simd_instruction_sets import (
+    get_supported_instruction_sets,
+    get_vector_instruction_set,
+)
+from pystencils.cpu.cpujit import get_compiler_config
+from pystencils.enums import Target
+from pystencils.rng import PhiloxTwoDoubles
+
 from lbmpy.creationfunctions import *
-from lbmpy.forcemodels import Guo
+from lbmpy.forcemodels import Guo 
 from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
 import numpy as np
 from lbmpy.enums import Stencil
 from lbmpy.moments import is_bulk_moment, is_shear_moment, get_order
 from lbmpy.stencils import LBStencil
-from pystencils.rng import PhiloxTwoDoubles
-
-import pytest
-from pystencils import get_code_str
-from pystencils.backends.simd_instruction_sets import get_supported_instruction_sets, get_vector_instruction_set
-from pystencils.cpu.cpujit import get_compiler_config
-from pystencils.enums import Target
 
 
 def _skip_instruction_sets_windows(instruction_sets):
-    if get_compiler_config()['os'] == 'windows':
+    if get_compiler_config()["os"] == "windows":
         # skip instruction sets supported by the CPU but not by the compiler
-        if 'avx' in instruction_sets and ('/arch:avx2' not in get_compiler_config()['flags'].lower()
-                                          and '/arch:avx512' not in get_compiler_config()['flags'].lower()):
-            instruction_sets.remove('avx')
-        if 'avx512' in instruction_sets and '/arch:avx512' not in get_compiler_config()['flags'].lower():
-            instruction_sets.remove('avx512')
+        if "avx" in instruction_sets and (
+            "/arch:avx2" not in get_compiler_config()["flags"].lower()
+            and "/arch:avx512" not in get_compiler_config()["flags"].lower()
+        ):
+            instruction_sets.remove("avx")
+        if (
+            "avx512" in instruction_sets
+            and "/arch:avx512" not in get_compiler_config()["flags"].lower()
+        ):
+            instruction_sets.remove("avx512")
     return instruction_sets
 
 
 def single_component_maxwell(x1, x2, kT, mass):
     """Integrate the probability density from x1 to x2 using the trapezoidal rule"""
     x = np.linspace(x1, x2, 1000)
-    return np.trapz(np.exp(-mass * x ** 2 / (2. * kT)), x) / np.sqrt(2. * np.pi * kT / mass)
+
+    try:
+        trapezoid = np.trapezoid  # since numpy 2.0
+    except AttributeError:
+        trapezoid = np.trapz
+
+    return trapezoid(np.exp(-mass * x**2 / (2.0 * kT)), x) / np.sqrt(
+        2.0 * np.pi * kT / mass
+    )
 
 
 def rr_getter(moment_group):
@@ -63,53 +81,86 @@ def second_order_moment_tensor_assignments(function_values, stencil, output_fiel
     """Assignments for calculating the pressure tensor"""
     assert len(function_values) == len(stencil)
     dim = len(stencil[0])
-    return [ps.Assignment(output_field(i, j),
-                          sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)))
-            for i in range(dim) for j in range(dim)]
+    return [
+        ps.Assignment(
+            output_field(i, j),
+            sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)),
+        )
+        for i in range(dim)
+        for j in range(dim)
+    ]
 
 
 def add_pressure_output_to_collision_rule(collision_rule, pressure_field):
-    pressure_ouput = second_order_moment_tensor_assignments(collision_rule.method.pre_collision_pdf_symbols,
-                                                            collision_rule.method.stencil, pressure_field)
+    pressure_ouput = second_order_moment_tensor_assignments(
+        collision_rule.method.pre_collision_pdf_symbols,
+        collision_rule.method.stencil,
+        pressure_field,
+    )
     collision_rule.main_assignments = collision_rule.main_assignments + pressure_ouput
 
 
-def get_fluctuating_lb(size=None, kT=None,
-                       omega_shear=None, omega_bulk=None, omega_odd=None, omega_even=None,
-                       rho_0=None, target=None):
+def get_fluctuating_lb(
+    size=None,
+    kT=None,
+    omega_shear=None,
+    omega_bulk=None,
+    omega_odd=None,
+    omega_even=None,
+    rho_0=None,
+    target=None,
+    zero_centered: bool = False,
+):
     # Parameters
     stencil = LBStencil(Stencil.D3Q19)
 
     # Setup data handling
-    dh = ps.create_data_handling((size,) * stencil.D, periodicity=True, default_target=target)
-    src = dh.add_array('src', values_per_cell=stencil.Q, layout='f')
-    dst = dh.add_array_like('dst', 'src')
-    rho = dh.add_array('rho', layout='f', latex_name='\\rho', values_per_cell=1)
-    u = dh.add_array('u', values_per_cell=dh.dim, layout='f')
-    pressure_field = dh.add_array('pressure', values_per_cell=(
-        3, 3), layout='f', gpu=target == Target.GPU)
+    dh = ps.create_data_handling(
+        (size,) * stencil.D, periodicity=True, default_target=target
+    )
+    src = dh.add_array("src", values_per_cell=stencil.Q, layout="f")
+    dst = dh.add_array_like("dst", "src")
+    rho = dh.add_array("rho", layout="f", latex_name="\\rho", values_per_cell=1)
+    u = dh.add_array("u", values_per_cell=dh.dim, layout="f")
+    pressure_field = dh.add_array(
+        "pressure", values_per_cell=(3, 3), layout="f", gpu=target == Target.GPU
+    )
     force_field = dh.add_array(
-        'force', values_per_cell=stencil.D, layout='f', gpu=target == Target.GPU)
+        "force", values_per_cell=stencil.D, layout="f", gpu=target == Target.GPU
+    )
 
     # Method setup
-    lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, compressible=True,
-                           weighted=True, zero_centered=False, relaxation_rates=rr_getter,
-                           force_model=Guo(force=force_field.center_vector),
-                           fluctuating={'temperature': kT},
-                           kernel_type='collide_only')
+    lbm_config = LBMConfig(
+        stencil=stencil,
+        method=Method.MRT,
+        compressible=True,
+        weighted=True,
+        zero_centered=zero_centered,
+        relaxation_rates=rr_getter,
+        force_model=Guo(force=force_field.center_vector),
+        fluctuating={"temperature": kT},
+        kernel_type="collide_only",
+    )
 
     lb_method = create_lb_method(lbm_config=lbm_config)
     lbm_config.lb_method = lb_method
-    
-    lbm_opt = LBMOptimisation(symbolic_field=src, cse_global=True)
-    collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
-
-    add_pressure_output_to_collision_rule(collision_rule, pressure_field)
 
-    collision = create_lb_update_rule(collision_rule=collision_rule,
-                                      lbm_config=lbm_config, lbm_optimisation=lbm_opt)
-    stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
-                                                   {'density': rho, 'velocity': u})
+    lbm_opt = LBMOptimisation(symbolic_field=src, cse_global=True)
+    collision_rule = create_lb_collision_rule(
+        lbm_config=lbm_config, lbm_optimisation=lbm_opt
+    )
+
+    # add_pressure_output_to_collision_rule(collision_rule, pressure_field)
+
+    collision = create_lb_update_rule(
+        collision_rule=collision_rule, lbm_config=lbm_config, lbm_optimisation=lbm_opt
+    )
+    stream = create_stream_pull_with_output_kernel(
+        collision.method,
+        src,
+        dst,
+        {"density": rho, "velocity": u, "moment2": pressure_field},
+    )
 
     config = ps.CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
 
@@ -120,15 +171,18 @@ def get_fluctuating_lb(size=None, kT=None,
     sync_pdfs = dh.synchronization_function([src.name])
 
     # Initialization
-    init = macroscopic_values_setter(collision.method, velocity=(0,) * dh.dim,
-                                     pdfs=src.center_vector, density=rho.center)
+    init = macroscopic_values_setter(
+        collision.method,
+        velocity=(0,) * dh.dim,
+        pdfs=src.center_vector,
+        density=rho_0
+    )
     init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
 
     dh.fill(rho.name, rho_0)
     dh.fill(u.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
     dh.fill(u.name, 0)
-    dh.fill(force_field.name, np.nan,
-            ghost_layers=True, inner_ghost_layers=True)
+    dh.fill(force_field.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
     dh.fill(force_field.name, 0)
     dh.run_kernel(init_kernel)
 
@@ -136,8 +190,15 @@ def get_fluctuating_lb(size=None, kT=None,
     def time_loop(start, steps):
         dh.all_to_gpu()
         for i in range(start, start + steps):
-            dh.run_kernel(collision_kernel, omega_shear=omega_shear, omega_bulk=omega_bulk,
-                          omega_odd=omega_odd, omega_even=omega_even, seed=42, time_step=i)
+            dh.run_kernel(
+                collision_kernel,
+                omega_shear=omega_shear,
+                omega_bulk=omega_bulk,
+                omega_odd=omega_odd,
+                omega_even=omega_even,
+                seed=42,
+                time_step=i,
+            )
 
             sync_pdfs()
             dh.run_kernel(stream_kernel)
@@ -148,13 +209,27 @@ def get_fluctuating_lb(size=None, kT=None,
     return dh, time_loop
 
 
-def test_resting_fluid(target=Target.CPU):
-    rho_0 = 0.86
-    kT = 4E-4
-    L = [60] * 3
-    dh, time_loop = get_fluctuating_lb(size=L[0], target=target,
-                                       rho_0=rho_0, kT=kT,
-                                       omega_shear=0.8, omega_bulk=0.5, omega_even=.04, omega_odd=0.3)
+@pytest.mark.parametrize(
+    "zero_centered", [False, True], ids=["regular-storage", "zero-centered"]
+)
+@pytest.mark.parametrize(
+    "domain_size", [8, 60]
+)
+def test_resting_fluid(zero_centered: bool, domain_size: int, target=Target.CPU):
+    rho_0 = 0.86 
+    kT = 4e-4
+    L = [domain_size] * 3 
+    dh, time_loop = get_fluctuating_lb(
+        size=L[0],
+        target=target,
+        rho_0=rho_0,
+        kT=kT,
+        omega_shear=0.8,
+        omega_bulk=0.5,
+        omega_even=0.04,
+        omega_odd=0.3,
+        zero_centered=zero_centered,
+    )
 
     # Test
     t = 0
@@ -168,38 +243,43 @@ def test_resting_fluid(target=Target.CPU):
         res_u = dh.gather_array("u").reshape((-1, 3))
         res_rho = dh.gather_array("rho").reshape((-1,))
 
-        # mass conservation
+        # mass conservationo
+        # density per cell fluctuates, but toal mass is conserved
         np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3E-12)
 
         # momentum conservation
         momentum = np.dot(res_rho, res_u)
-        np.testing.assert_allclose(momentum, [0, 0, 0], atol=1E-10)
+        np.testing.assert_allclose(momentum, [0, 0, 0], atol=1e-10)
 
-        # temperature
-        kinetic_energy = 1 / 2 * np.dot(res_rho, res_u * res_u) / np.product(L)
-        np.testing.assert_allclose(
-            kinetic_energy, [kT / 2] * 3, atol=kT * 0.01)
+        # temperature (fluctuates around pre-set kT)
+        kinetic_energy = 1 / 2 * np.dot(res_rho, res_u * res_u) / np.prod(L)
+        kT_tol = 0.075 *(16/domain_size)**(3/2) 
+        np.testing.assert_allclose(kinetic_energy, [kT / 2] * 3, rtol=kT_tol)
 
         # velocity distribution
         v_hist, v_bins = np.histogram(
-            res_u, bins=11, range=(-.075, .075), density=True)
+            res_u, bins=11, range=(-0.075, 0.075), density=True
+        )
 
         # Calculate expected values from single
         v_expected = []
         for j in range(len(v_hist)):
             # Maxwell distribution
-            res = 1. / (v_bins[j + 1] - v_bins[j]) * \
-                  single_component_maxwell(
-                      v_bins[j], v_bins[j + 1], kT, rho_0)
+            res = (
+                1.0
+                / (v_bins[j + 1] - v_bins[j])
+                * single_component_maxwell(v_bins[j], v_bins[j + 1], kT, rho_0)
+            )
             v_expected.append(res)
         v_expected = np.array(v_expected)
 
-        # 10% accuracy on the entire histogram
-        np.testing.assert_allclose(v_hist, v_expected, rtol=0.1)
-        # 1% accuracy on the middle part
+        hist_tol_all = 0.75 *(16/domain_size)**(3/2)
+        np.testing.assert_allclose(v_hist, v_expected, rtol=hist_tol_all)
+        hist_tol_center = hist_tol_all/10
         remove = 3
         np.testing.assert_allclose(
-            v_hist[remove:-remove], v_expected[remove:-remove], rtol=0.01)
+            v_hist[remove:-remove], v_expected[remove:-remove], rtol=hist_tol_center
+        )
 
         # pressure tensor against expressions from
         # Duenweg, Schiller, Ladd, https://arxiv.org/abs/0707.1581
@@ -212,19 +292,35 @@ def test_resting_fluid(target=Target.CPU):
         # Diagonal elements are rho c_s^22 +<u,u>. When the fluid is
         # thermalized, the expectation value of <u,u> = kT due to the
         # equi-partition theorem.
-        p_av_expected = np.diag([rho_0 * c_s ** 2 + kT] * 3)
+        p_av_expected = np.diag([rho_0 * c_s**2 + kT] * 3)
+        pressure_atol = c_s**2 / 200 *(16/domain_size)**(3/2)
         np.testing.assert_allclose(
-            np.mean(res_pressure, axis=0), p_av_expected, atol=c_s ** 2 / 2000)
+            np.mean(res_pressure, axis=0), p_av_expected, atol=pressure_atol)
+
 
 
-def test_point_force(target=Target.CPU):
+@pytest.mark.parametrize(
+    "zero_centered", [False, True], ids=["regular-storage", "zero-centered"]
+)
+@pytest.mark.parametrize(
+    "domain_size", [8, 60]
+)
+def test_point_force(zero_centered: bool, domain_size: int, target=Target.CPU):
     """Test momentum balance for thermalized fluid with applied poitn forces"""
     rho_0 = 0.86
-    kT = 4E-4
-    L = [8] * 3
-    dh, time_loop = get_fluctuating_lb(size=L[0], target=target,
-                                       rho_0=rho_0, kT=kT,
-                                       omega_shear=0.8, omega_bulk=0.5, omega_even=.04, omega_odd=0.3)
+    kT = 4e-4
+    L = [domain_size] * 3
+    dh, time_loop = get_fluctuating_lb(
+        size=L[0],
+        target=target,
+        rho_0=rho_0,
+        kT=kT,
+        omega_shear=0.8,
+        omega_bulk=0.5,
+        omega_even=0.8,
+        omega_odd=0.8,
+        zero_centered=zero_centered
+    )
 
     # Test
     t = 0
@@ -233,17 +329,17 @@ def test_point_force(target=Target.CPU):
 
     introduced_momentum = np.zeros(3)
     for i in range(100):
-        point_force = 1E-5 * (np.random.random(3) - .5)
+        point_force = 1e-2/domain_size**(3/2) * (np.random.random(3) - 0.5)
         introduced_momentum += point_force
 
         # Note that ghost layers are included in the indexing
         force_pos = np.random.randint(1, L[0] - 2, size=3)
 
-        dh.cpu_arrays["force"][force_pos[0],
-                               force_pos[1], force_pos[2]] = point_force
+        dh.cpu_arrays["force"][force_pos[0], force_pos[1], force_pos[2]] = point_force
         t = time_loop(t, 1)
         res_u = dh.gather_array("u").reshape((-1, 3))
         res_rho = dh.gather_array("rho").reshape((-1,))
+        
 
         # mass conservation
         np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3E-12)
@@ -251,86 +347,127 @@ def test_point_force(target=Target.CPU):
         # momentum conservation
         momentum = np.dot(res_rho, res_u)
         np.testing.assert_allclose(
-            momentum, introduced_momentum + 0.5 * point_force, atol=1E-10)
-        dh.cpu_arrays["force"][force_pos[0],
-                               force_pos[1], force_pos[2]] = np.zeros(3)
-
-
-@pytest.mark.skipif(not get_supported_instruction_sets(), reason="No vector instruction sets supported")
-@pytest.mark.parametrize('data_type', ("float32", "float64"))
-@pytest.mark.parametrize('assume_aligned', (True, False))
-@pytest.mark.parametrize('assume_inner_stride_one', (True, False))
-@pytest.mark.parametrize('assume_sufficient_line_padding', (True, False))
-def test_vectorization(data_type, assume_aligned, assume_inner_stride_one, assume_sufficient_line_padding):
+            momentum, introduced_momentum + 0.5 * point_force, atol=1e-10
+        )
+        dh.cpu_arrays["force"][force_pos[0], force_pos[1], force_pos[2]] = np.zeros(3)
+
+
+@pytest.mark.skipif(
+    not get_supported_instruction_sets(), reason="No vector instruction sets supported"
+)
+@pytest.mark.parametrize("data_type", ("float32", "float64"))
+@pytest.mark.parametrize("assume_aligned", (True, False))
+@pytest.mark.parametrize("assume_inner_stride_one", (True, False))
+@pytest.mark.parametrize("assume_sufficient_line_padding", (True, False))
+def test_vectorization(
+    data_type, assume_aligned, assume_inner_stride_one, assume_sufficient_line_padding
+):
     stencil = LBStencil(Stencil.D3Q19)
-    pdfs, pdfs_tmp = ps.fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[3D]", layout='fzyx')
+    pdfs, pdfs_tmp = ps.fields(
+        f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[3D]", layout="fzyx"
+    )
 
     method = create_mrt_orthogonal(
-        stencil=stencil,
+        stencil=stencil, compressible=True, weighted=True, relaxation_rates=rr_getter
+    )
+
+    rng_node = (
+        ps.rng.PhiloxTwoDoubles if data_type == "float64" else ps.rng.PhiloxFourFloats
+    )
+    lbm_config = LBMConfig(
+        lb_method=method,
+        fluctuating={
+            "temperature": sp.Symbol("kT"),
+            "rng_node": rng_node,
+            "block_offsets": tuple([0] * stencil.D),
+        },
         compressible=True,
-        weighted=True,
-        relaxation_rates=rr_getter)
-
-    rng_node = ps.rng.PhiloxTwoDoubles if data_type == "float64" else ps.rng.PhiloxFourFloats
-    lbm_config = LBMConfig(lb_method=method, fluctuating={'temperature': sp.Symbol("kT"),
-                                                          'rng_node': rng_node,
-                                                          'block_offsets': tuple([0] * stencil.D)},
-                           compressible=True, zero_centered=False, 
-                           stencil=method.stencil, kernel_type='collide_only')
-    lbm_opt = LBMOptimisation(cse_global=True, symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp)
+        zero_centered=False,
+        stencil=method.stencil,
+        kernel_type="collide_only",
+    )
+    lbm_opt = LBMOptimisation(
+        cse_global=True, symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp
+    )
 
     collision = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
 
     instruction_sets = _skip_instruction_sets_windows(get_supported_instruction_sets())
     instruction_set = instruction_sets[-1]
 
-    config = ps.CreateKernelConfig(target=Target.CPU,
-                                   data_type=data_type, default_number_float=data_type,
-                                   cpu_vectorize_info={'instruction_set': instruction_set,
-                                                       'assume_aligned': assume_aligned,
-                                                       'assume_inner_stride_one': assume_inner_stride_one,
-                                                       'assume_sufficient_line_padding': assume_sufficient_line_padding,
-                                                       }
-                                   )
-
-    if not assume_inner_stride_one and 'storeS' not in get_vector_instruction_set(data_type, instruction_set):
-        with pytest.warns(UserWarning) as warn:
+    config = ps.CreateKernelConfig(
+        target=Target.CPU,
+        data_type=data_type,
+        default_number_float=data_type,
+        cpu_vectorize_info={
+            "instruction_set": instruction_set,
+            "assume_aligned": assume_aligned,
+            "assume_inner_stride_one": assume_inner_stride_one,
+            "assume_sufficient_line_padding": assume_sufficient_line_padding,
+        },
+    )
+
+    if not assume_inner_stride_one and "storeS" not in get_vector_instruction_set(
+        data_type, instruction_set
+    ):
+        with pytest.warns(UserWarning) as pytest_warnings:
             ast = ps.create_kernel(collision, config=config)
-            assert 'Could not vectorize loop' in warn[0].message.args[0]
+            assert "Could not vectorize loop" in pytest_warnings[0].message.args[0]
     else:
-        with pytest.warns(None) as warn:
-            ast = ps.create_kernel(collision, config=config)
-            assert len(warn) == 0
+        ast = ps.create_kernel(collision, config=config)
     ast.compile()
     code = get_code_str(ast)
     print(code)
 
 
-@pytest.mark.parametrize('data_type', ("float32", "float64"))
-@pytest.mark.parametrize('assume_aligned', (True, False))
-@pytest.mark.parametrize('assume_inner_stride_one', (True, False))
-@pytest.mark.parametrize('assume_sufficient_line_padding', (True, False))
-def test_fluctuating_lb_issue_188_wlb(data_type, assume_aligned,
-                                      assume_inner_stride_one, assume_sufficient_line_padding):
+@pytest.mark.parametrize("data_type", ("float32", "float64"))
+@pytest.mark.parametrize("assume_aligned", (True, False))
+@pytest.mark.parametrize("assume_inner_stride_one", (True, False))
+@pytest.mark.parametrize("assume_sufficient_line_padding", (True, False))
+def test_fluctuating_lb_issue_188_wlb(
+    data_type, assume_aligned, assume_inner_stride_one, assume_sufficient_line_padding
+):
     stencil = LBStencil(Stencil.D3Q19)
     temperature = sp.symbols("temperature")
-    pdfs, pdfs_tmp = ps.fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[3D]", layout='fzyx')
-
-    rng_node = ps.rng.PhiloxTwoDoubles if data_type == "float64" else ps.rng.PhiloxFourFloats
-    fluctuating = {'temperature': temperature,
-                   'block_offsets': 'walberla',
-                   'rng_node': rng_node}
-
-    lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, compressible=True,
-                           weighted=True, zero_centered=False, relaxation_rate=1.4,
-                           fluctuating=fluctuating)
-    lbm_opt = LBMOptimisation(symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp, cse_global=True)
+    pdfs, pdfs_tmp = ps.fields(
+        f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[3D]", layout="fzyx"
+    )
+
+    rng_node = (
+        ps.rng.PhiloxTwoDoubles if data_type == "float64" else ps.rng.PhiloxFourFloats
+    )
+    fluctuating = {
+        "temperature": temperature,
+        "block_offsets": "walberla",
+        "rng_node": rng_node,
+    }
+
+    lbm_config = LBMConfig(
+        stencil=stencil,
+        method=Method.MRT,
+        compressible=True,
+        weighted=True,
+        zero_centered=False,
+        relaxation_rate=1.4,
+        fluctuating=fluctuating,
+    )
+    lbm_opt = LBMOptimisation(
+        symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp, cse_global=True
+    )
 
     up = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
 
-    cpu_vectorize_info = {'instruction_set': 'avx', 'assume_inner_stride_one': True, 'assume_aligned': True}
-    config = ps.CreateKernelConfig(target=ps.Target.CPU, data_type=data_type, default_number_float=data_type,
-                                   cpu_vectorize_info=cpu_vectorize_info)
+    cpu_vectorize_info = {
+        "instruction_set": "avx",
+        "assume_inner_stride_one": True,
+        "assume_aligned": True,
+    }
+    config = ps.CreateKernelConfig(
+        target=ps.Target.CPU,
+        data_type=data_type,
+        default_number_float=data_type,
+        cpu_vectorize_info=cpu_vectorize_info,
+    )
 
     ast = create_kernel(up, config=config)
     code = ps.get_code_str(ast)

lbmpy_tests/test_force.py → tests/test_force.py

@@ -12,10 +12,9 @@ from lbmpy.macroscopic_value_kernels import macroscopic_values_setter, macroscop
 from lbmpy.moments import (is_bulk_moment, moments_up_to_component_order,
                            exponents_to_polynomial_representations, exponent_tuple_sort_key)
 from lbmpy.stencils import LBStencil
-from lbmpy.updatekernels import create_stream_pull_with_output_kernel
 
 # all force models available are defined in the ForceModel enum, but Cumulant is not a "real" force model
-force_models = [f for f in ForceModel]
+force_models = [f for f in ForceModel if f is not ForceModel.CENTRALMOMENT]
 
 
 @pytest.mark.parametrize("method_enum", [Method.SRT, Method.TRT, Method.MRT, Method.CUMULANT])
@@ -25,7 +24,7 @@ force_models = [f for f in ForceModel]
 def test_total_momentum(method_enum, zero_centered, force_model, omega):
     if method_enum == Method.CUMULANT and \
             force_model not in (ForceModel.SIMPLE, ForceModel.LUO, ForceModel.GUO, ForceModel.HE):
-        return True
+        return
 
     L = (16, 16)
     stencil = LBStencil(Stencil.D2Q9)
@@ -198,7 +197,7 @@ def test_literature(force_model, stencil, method):
     assert len(omega_momentum) == 1
     omega_momentum = omega_momentum[0]
 
-    subs_dict = lb_method.subs_dict_relxation_rate
+    subs_dict = lb_method.subs_dict_relaxation_rate
     force_term = sp.simplify(lb_method.force_model(lb_method).subs(subs_dict))
     u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)
     rho = lb_method.conserved_quantity_computation.density_symbol
@@ -283,7 +282,7 @@ def test_modes_central_moment(force_model, compressible):
                            compressible=compressible, force_model=force_model, force=tuple(F))
     method = create_lb_method(lbm_config=lbm_config)
 
-    subs_dict = method.subs_dict_relxation_rate
+    subs_dict = method.subs_dict_relaxation_rate
     force_moments = method.force_model.central_moment_space_forcing(method)
     force_moments = force_moments.subs(subs_dict)
 
@@ -309,9 +308,9 @@ def test_symmetric_forcing_equivalence(force_model, compressible):
                            nested_moments=moment_polys, compressible=True, force_model=force_model, force=tuple(F))
     method = create_lb_method(lbm_config=lbm_config)
     if not method.force_model.has_symmetric_central_moment_forcing:
-        return True
+        return
 
-    subs_dict = method.subs_dict_relxation_rate
+    subs_dict = method.subs_dict_relaxation_rate
     force_moments = method.force_model.central_moment_space_forcing(method)
     force_moments = force_moments.subs(subs_dict)
 
@@ -336,7 +335,7 @@ def test_modes_central_moment_longrun(stencil, force_model, compressible):
                            compressible=compressible, force_model=force_model, force=tuple(F))
     method = create_lb_method(lbm_config=lbm_config)
 
-    subs_dict = method.subs_dict_relxation_rate
+    subs_dict = method.subs_dict_relaxation_rate
     force_moments = method.force_model.moment_space_forcing(method)
     force_moments = force_moments.subs(subs_dict)
 
@@ -360,7 +359,7 @@ def _check_modes(stencil, force_model, compressible):
                            compressible=compressible, force_model=force_model, force=tuple(F))
     method = create_lb_method(lbm_config=lbm_config)
 
-    subs_dict = method.subs_dict_relxation_rate
+    subs_dict = method.subs_dict_relaxation_rate
     force_moments = method.force_model.moment_space_forcing(method)
     force_moments = force_moments.subs(subs_dict)
 

lbmpy_tests/test_gpu_block_size_limiting.py → tests/test_gpu_block_size_limiting.py

@@ -6,7 +6,7 @@ from pystencils import Target, CreateKernelConfig
 
 
 def test_gpu_block_size_limiting():
-    pytest.importorskip("pycuda")
+    pytest.importorskip("cupy")
     too_large = 2048*2048
     lbm_config = LBMConfig(method=Method.CUMULANT, stencil=LBStencil(Stencil.D3Q19),
                            relaxation_rate=1.8, compressible=True)
@@ -17,5 +17,5 @@ def test_gpu_block_size_limiting():
     kernel = ast.compile()
     assert all(b < too_large for b in limited_block_size['block'])
     bs = [too_large, too_large, too_large]
-    ast.indexing.limit_block_size_by_register_restriction(bs, kernel.num_regs)
+    bs = ast.indexing.limit_block_size_by_register_restriction(bs, kernel.num_regs)
     assert all(b < too_large for b in bs)