tests/test_interpolation_boundaries.py

+import pytest
+import numpy as np
+from pystencils.slicing import slice_from_direction
+
+from lbmpy import LBMConfig, LBStencil, Stencil, Method
+from lbmpy.boundaries.boundaryconditions import NoSlip, NoSlipLinearBouzidi, QuadraticBounceBack, UBB
+from lbmpy.lbstep import LatticeBoltzmannStep
+
+
+def check_velocity(noslip_velocity, interpolated_velocity, wall_distance):
+    # First we check if the flow is fully developed
+    np.testing.assert_almost_equal(np.gradient(np.gradient(noslip_velocity)), 0, decimal=8)
+    np.testing.assert_almost_equal(np.gradient(np.gradient(interpolated_velocity)), 0, decimal=8)
+
+    # If the wall is closer to the first fluid cell we expect a lower velocity at the first fluid cell
+    if wall_distance < 0.5:
+        assert noslip_velocity[0] > interpolated_velocity[0]
+    # If the wall is further away from the first fluid cell we expect a higher velocity at the first fluid cell
+    if wall_distance > 0.5:
+        assert noslip_velocity[0] < interpolated_velocity[0]
+    # If the wall cuts the cell halfway the interpolation BC should approximately fall back to a noslip boundary
+    if wall_distance == 0.5:
+        np.testing.assert_almost_equal(noslip_velocity[0], interpolated_velocity[0], decimal=2)
+
+
+def couette_flow(stencil, method_enum, zero_centered, wall_distance, compressible):
+    dim = stencil.D
+    if dim == 2:
+        domain_size = (50, 25)
+        wall_velocity = (0.01, 0)
+        periodicity = (True, False)
+    else:
+        domain_size = (50, 25, 25)
+        wall_velocity = (0.01, 0, 0)
+        periodicity = (True, False, True)
+
+    timesteps = 10000
+    omega = 1.1
+
+    lbm_config = LBMConfig(stencil=stencil,
+                           method=method_enum,
+                           relaxation_rate=omega,
+                           zero_centered=zero_centered,
+                           compressible=compressible,
+                           weighted=True)
+
+    lb_step_noslip = LatticeBoltzmannStep(domain_size=domain_size, periodicity=periodicity,
+                                          lbm_config=lbm_config, compute_velocity_in_every_step=True)
+    lb_step_bouzidi = LatticeBoltzmannStep(domain_size=domain_size, periodicity=periodicity,
+                                           lbm_config=lbm_config, compute_velocity_in_every_step=True)
+
+    lb_step_quadratic_bb = LatticeBoltzmannStep(domain_size=domain_size, periodicity=periodicity,
+                                                lbm_config=lbm_config, compute_velocity_in_every_step=True)
+
+    def init_wall_distance(boundary_data, **_):
+        for cell in boundary_data.index_array:
+            cell['q'] = wall_distance
+
+    moving_wall = UBB(wall_velocity)
+    noslip = NoSlip("wall")
+    bouzidi = NoSlipLinearBouzidi("wall", init_wall_distance=init_wall_distance)
+    quadratic_bb = QuadraticBounceBack(omega, "wall", init_wall_distance=init_wall_distance)
+
+    lb_step_noslip.boundary_handling.set_boundary(noslip, slice_from_direction('S', dim))
+    lb_step_noslip.boundary_handling.set_boundary(moving_wall, slice_from_direction('N', dim))
+
+    lb_step_bouzidi.boundary_handling.set_boundary(bouzidi, slice_from_direction('S', dim))
+    lb_step_bouzidi.boundary_handling.set_boundary(moving_wall, slice_from_direction('N', dim))
+
+    lb_step_quadratic_bb.boundary_handling.set_boundary(quadratic_bb, slice_from_direction('S', dim))
+    lb_step_quadratic_bb.boundary_handling.set_boundary(moving_wall, slice_from_direction('N', dim))
+
+    lb_step_noslip.run(timesteps)
+    lb_step_bouzidi.run(timesteps)
+    lb_step_quadratic_bb.run(timesteps)
+
+    if dim == 2:
+        noslip_velocity = lb_step_noslip.velocity[domain_size[0] // 2, :, 0]
+        bouzidi_velocity = lb_step_bouzidi.velocity[domain_size[0] // 2, :, 0]
+        quadratic_bb_velocity = lb_step_quadratic_bb.velocity[domain_size[0] // 2, :, 0]
+    else:
+        noslip_velocity = lb_step_noslip.velocity[domain_size[0] // 2, :, domain_size[2] // 2, 0]
+        bouzidi_velocity = lb_step_bouzidi.velocity[domain_size[0] // 2, :, domain_size[2] // 2, 0]
+        quadratic_bb_velocity = lb_step_quadratic_bb.velocity[domain_size[0] // 2, :, domain_size[2] // 2, 0]
+
+    check_velocity(noslip_velocity, bouzidi_velocity, wall_distance)
+    check_velocity(noslip_velocity, quadratic_bb_velocity, wall_distance)
+
+
+@pytest.mark.parametrize("zero_centered", [False, True])
+@pytest.mark.parametrize("wall_distance", [0.1, 0.5, 0.9])
+@pytest.mark.parametrize("compressible", [True, False])
+def test_couette_flow_short(zero_centered, wall_distance, compressible):
+    stencil = LBStencil(Stencil.D2Q9)
+    couette_flow(stencil, Method.SRT, zero_centered, wall_distance, compressible)
+
+
+@pytest.mark.parametrize("method_enum", [Method.MRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
+@pytest.mark.parametrize("zero_centered", [False, True])
+@pytest.mark.parametrize("wall_distance", [0.1, 0.5, 0.9])
+@pytest.mark.parametrize("compressible", [True, False])
+@pytest.mark.longrun
+def test_couette_flow_long(method_enum, zero_centered, wall_distance, compressible):
+    if method_enum is Method.CUMULANT and compressible is False:
+        pytest.skip("incompressible cumulant is not supported")
+
+    stencil = LBStencil(Stencil.D2Q9)
+    couette_flow(stencil, method_enum, zero_centered, wall_distance, compressible)
+
+
+@pytest.mark.parametrize("method_enum", [Method.SRT, Method.MRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
+@pytest.mark.parametrize("wall_distance", [0.1, 0.5, 0.9])
+@pytest.mark.parametrize("stencil", [Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.longrun
+def test_couette_flow_d3d(method_enum, wall_distance, stencil):
+    stencil = LBStencil(stencil)
+    couette_flow(stencil, method_enum, True, wall_distance, True)
+
+

tests/test_json_serializer.py

+"""
+Test the lbmpy-specific JSON encoder and serializer as used in the Database class.
+"""
+
+import tempfile
+import pystencils as ps
+
+from lbmpy import Stencil, Method, ForceModel, SubgridScaleModel
+from lbmpy.advanced_streaming import Timestep
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, LBStencil
+from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor
+
+from pystencils.runhelper import Database
+from lbmpy.db import LbmpyJsonSerializer
+
+
+def test_json_serializer():
+
+    stencil = LBStencil(Stencil.D3Q27)
+    q = stencil.Q
+    pdfs, pdfs_tmp = ps.fields(f"pdfs({q}), pdfs_tmp({q}): double[3D]", layout='fzyx')
+    density = ps.fields(f"rho: double[3D]", layout='fzyx')
+
+    from lbmpy.non_newtonian_models import CassonsParameters
+    cassons_params = CassonsParameters(0.2)
+
+    # create dummy lbmpy config
+    lbm_config = LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, force_model=ForceModel.GUO,
+                           compressible=True, relaxation_rate=1.999, subgrid_scale_model=SubgridScaleModel.SMAGORINSKY,
+                           galilean_correction=True, cassons=cassons_params, density_input=density,
+                           kernel_type=StreamPullTwoFieldsAccessor, timestep=Timestep.BOTH)
+
+    lbm_optimization = LBMOptimisation(cse_pdfs=False, cse_global=False, builtin_periodicity=(True, False, False),
+                                       symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp)
+
+    # create dummy database
+    temp_dir = tempfile.TemporaryDirectory()
+    db = Database(file=temp_dir.name, serializer_info=('lbmpy_serializer', LbmpyJsonSerializer))
+
+    db.save(params={'lbm_config': lbm_config, 'lbm_optimization': lbm_optimization}, result={'test': 'dummy'})

lbmpy_tests/test_lbstep.py → tests/test_lbstep.py

@@ -4,9 +4,10 @@ from types import MappingProxyType
 from pystencils import Target, CreateKernelConfig
 
 from lbmpy.scenarios import create_fully_periodic_flow, create_lid_driven_cavity
+from lbmpy import SubgridScaleModel
 
 try:
-    import pycuda.driver
+    import cupy
     gpu_available = True
 except ImportError:
     gpu_available = False
@@ -77,7 +78,13 @@ def test_data_handling_2d():
 
 
 def test_smagorinsky_setup():
-    step = create_lid_driven_cavity((30, 30), smagorinsky=0.16, relaxation_rate=1.99)
+    step = create_lid_driven_cavity((30, 30), subgrid_scale_model=(SubgridScaleModel.SMAGORINSKY, 0.16),
+                                    relaxation_rate=1.99)
+    step.run(10)
+
+def test_qr_setup():
+    step = create_lid_driven_cavity((30, 30), subgrid_scale_model=(SubgridScaleModel.QR, 0.33),
+                                    relaxation_rate=1.99)
     step.run(10)
 
 

lbmpy_tests/test_macroscopic_value_kernels.py → tests/test_macroscopic_value_kernels.py

 import pytest
 import numpy as np
+import sympy as sp
 
 from lbmpy.creationfunctions import create_lb_method, LBMConfig
 from lbmpy.enums import Stencil, ForceModel, Method
 from lbmpy.macroscopic_value_kernels import (
-    compile_macroscopic_values_getter, compile_macroscopic_values_setter)
+    compile_macroscopic_values_getter, compile_macroscopic_values_setter, strain_rate_tensor_getter)
 from lbmpy.advanced_streaming.utility import streaming_patterns, Timestep
 from lbmpy.stencils import LBStencil
 
@@ -83,3 +84,18 @@ def test_set_get_constant_velocity(stencil, force_model, compressible):
         computed_velocity = velocity_output_field[0, 0, 0, :]
 
     np.testing.assert_almost_equal(np.array(ref_velocity[:method.dim]), computed_velocity)
+
+
+@pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19])
+@pytest.mark.parametrize("method_enum", [Method.SRT, Method.TRT, Method.MRT, Method.CUMULANT])
+def test_get_strain_rate_tensor(stencil, method_enum):
+    stencil = LBStencil(stencil)
+    lbm_config = LBMConfig(stencil=stencil, method=method_enum, compressible=True)
+    method = create_lb_method(lbm_config=lbm_config)
+
+    pdfs = sp.symbols(f"f_:{stencil.Q}")
+    strain_rate_tensor = sp.symbols(f"strain_rate_tensor_:{method.dim * method.dim}")
+
+    getter_assignments = strain_rate_tensor_getter(method, strain_rate_tensor, pdfs)
+
+    assert len(getter_assignments) == method.dim * method.dim

tests/test_mrt_simplifications.py

+import pytest
+
+import sympy as sp
+
+from pystencils.sympyextensions import is_constant
+
+from lbmpy import Stencil, LBStencil, Method, create_lb_collision_rule, LBMConfig, LBMOptimisation
+
+@pytest.mark.parametrize('method', [Method.MRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
+def test_mrt_simplifications(method: Method):
+    stencil = Stencil.D3Q19
+    lbm_config = LBMConfig(stencil=stencil, method=method, compressible=True)
+    lbm_opt = LBMOptimisation(simplification='auto')
+
+    cr = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+
+    
+    for subexp in cr.subexpressions:
+        rhs = subexp.rhs
+        #   Check for aliases
+        assert not isinstance(rhs, sp.Symbol)
+
+        #   Check for logarithms
+        assert not rhs.atoms(sp.log)
+
+        #   Check for nonextracted constant summands or factors
+        exprs = rhs.atoms(sp.Add, sp.Mul)
+        for expr in exprs:
+            for arg in expr.args:
+                if isinstance(arg, sp.Number):
+                    if arg not in {sp.Number(1), sp.Number(-1), sp.Float(1), sp.Float(-1)}:
+                        breakpoint()
+                    
+        #   Check for divisions
+        if not (isinstance(rhs, sp.Pow) and rhs.args[1] < 0):
+            powers = rhs.atoms(sp.Pow)
+            for p in powers:
+                assert p.args[1] > 0

lbmpy_tests/test_poisuille_channel.py → tests/test_poisuille_channel.py

@@ -7,7 +7,7 @@ import pytest
 from lbmpy.enums import Stencil, CollisionSpace
 
 import pystencils as ps
-from poiseuille import poiseuille_channel
+from . import poiseuille
 
 
 @pytest.mark.parametrize('target', (ps.Target.CPU, ps.Target.GPU))
@@ -18,7 +18,10 @@ def test_poiseuille_channel(target, stencil_name, zero_centered, moment_space_co
     # Cuda
     if target == ps.Target.GPU:
         import pytest
-        pytest.importorskip("pycuda")
+        pytest.importorskip("cupy")
 
     cspace_info = CollisionSpace.RAW_MOMENTS if moment_space_collision else CollisionSpace.POPULATIONS
-    poiseuille_channel(target=target, stencil_name=stencil_name, zero_centered=zero_centered, collision_space_info=cspace_info)
+    poiseuille.poiseuille_channel(target=target,
+                                  stencil_name=stencil_name,
+                                  zero_centered=zero_centered,
+                                  collision_space_info=cspace_info)

tests/test_psm.py

+import pytest
+import numpy as np
+
+from pystencils import fields, CreateKernelConfig, Target, create_kernel, get_code_str, create_data_handling
+from lbmpy import LBMConfig, Stencil, Method, LBStencil, create_lb_method, create_lb_collision_rule, LBMOptimisation, \
+    create_lb_update_rule
+from lbmpy.partially_saturated_cells import PSMConfig
+
+
+@pytest.mark.parametrize("stencil", [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
+@pytest.mark.parametrize("method_enum", [Method.SRT, Method.MRT, Method.CUMULANT])
+def test_psm_integration(stencil, method_enum):
+    stencil = LBStencil(stencil)
+
+    fraction_field = fields("fraction_field: double[10, 10, 5]", layout=(2, 1, 0))
+    object_vel = fields("object_vel(3): double[10, 10, 5]", layout=(3, 2, 1, 0))
+    psm_config = PSMConfig(fraction_field=fraction_field, object_velocity_field=object_vel)
+
+    lbm_config = LBMConfig(stencil=stencil, method=method_enum, relaxation_rate=1.5, compressible=True,
+                           psm_config=psm_config)
+    config = CreateKernelConfig(target=Target.CPU)#
+
+    collision_rule = create_lb_collision_rule(lbm_config=lbm_config)
+
+    ast = create_kernel(collision_rule, config=config)
+    code_str = get_code_str(ast)
+
+
+def get_data_handling_for_psm(use_psm):
+    domain_size = (10, 10)
+    stencil = LBStencil(Stencil.D2Q9)
+    dh = create_data_handling(domain_size=domain_size, periodicity=(True, False))
+
+    f = dh.add_array('f', values_per_cell=len(stencil))
+    dh.fill(f.name, 0.0, ghost_layers=True)
+    f_tmp = dh.add_array('f_tmp', values_per_cell=len(stencil))
+    dh.fill(f_tmp.name, 0.0, ghost_layers=True)
+
+    psm_config = None
+    if use_psm:
+        fraction_field = dh.add_array('fraction_field', values_per_cell=1)
+        dh.fill(fraction_field.name, 0.0, ghost_layers=True)
+        object_vel = dh.add_array('object_vel', values_per_cell=dh.dim)
+        dh.fill(object_vel.name, 0.0, ghost_layers=True)
+        psm_config = PSMConfig(fraction_field=fraction_field, object_velocity_field=object_vel)
+
+    lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.5, psm_config=psm_config)
+
+    lbm_optimisation = LBMOptimisation(symbolic_field=f, symbolic_temporary_field=f_tmp)
+    update_rule = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_optimisation)
+    kernel_lb_step_with_psm = create_kernel(update_rule).compile()
+    return (dh, kernel_lb_step_with_psm)
+
+
+def test_lbm_vs_psm():
+
+    psm_dh, psm_kernel = get_data_handling_for_psm(True)
+    lbm_dh, lbm_kernel = get_data_handling_for_psm(False)
+
+    for i in range(20):
+        psm_dh.run_kernel(psm_kernel)
+        psm_dh.swap('f', 'f_tmp')
+
+        lbm_dh.run_kernel(lbm_kernel)
+        lbm_dh.swap('f', 'f_tmp')
+
+    max_vel_error = np.max(np.abs(psm_dh.gather_array('f') - lbm_dh.gather_array('f')))
+    np.testing.assert_allclose(max_vel_error, 0, atol=1e-14)

lbmpy_tests/test_quicktests.py → tests/test_quicktests.py

@@ -5,11 +5,11 @@ import pystencils as ps
 from lbmpy.enums import ForceModel, Method, Stencil
 from lbmpy.scenarios import create_lid_driven_cavity
 
-from lbmpy_tests.poiseuille import poiseuille_channel
+from . import poiseuille
 
 
 def test_poiseuille_channel_quicktest():
-    poiseuille_channel(target=ps.Target.CPU, stencil_name=Stencil.D2Q9)
+    poiseuille.poiseuille_channel(target=ps.Target.CPU, stencil_name=Stencil.D2Q9)
 
 
 def test_entropic_methods():

lbmpy_tests/test_relaxation_rate.py → tests/test_relaxation_rate.py

 import pytest
 import sympy as sp
 from lbmpy.creationfunctions import create_lb_method, LBMConfig
+from lbmpy.moments import is_shear_moment, get_order
 from lbmpy.enums import Method, Stencil
 from lbmpy.relaxationrates import get_shear_relaxation_rate
 from lbmpy.stencils import LBStencil
@@ -19,3 +20,17 @@ def test_relaxation_rate():
                            relaxation_rates=omegas)
     method = create_lb_method(lbm_config=lbm_config)
     assert get_shear_relaxation_rate(method) == omegas[0]
+
+
+@pytest.mark.parametrize("method", [Method.MRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
+def test_default_mrt_behaviour(method):
+    lbm_config = LBMConfig(stencil=LBStencil(Stencil.D3Q19), method=method, compressible=True)
+    method = create_lb_method(lbm_config=lbm_config)
+
+    for moment, relax_info in method.relaxation_info_dict.items():
+        if get_order(moment) <= 1:
+            assert relax_info.relaxation_rate == 0
+        elif is_shear_moment(moment, method.dim):
+            assert relax_info.relaxation_rate == sp.Symbol('omega')
+        else:
+            assert relax_info.relaxation_rate == 1