diff --git a/doc/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb b/doc/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb
index 20e1653b9692f9488ad55a36113cf52d0e70be9c..6c2ff255e9826e28b1b5baebca6190c480a38996 100644
--- a/doc/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb
+++ b/doc/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb
@@ -419,7 +419,7 @@
    ],
    "source": [
     "method = create_lb_method(method='mrt', weighted=True, stencil='D2Q9', force=(1e-6, 0),\n",
-    "                          relaxation_rates=[0, 0, ω, 1.9, 1.9])\n",
+    "                          force_model='luo', relaxation_rates=[0, 0, ω, 1.9, 1.9])\n",
     "method"
    ]
   },
diff --git a/doc/sphinx/lbmpy.bib b/doc/sphinx/lbmpy.bib
index 9620dc345ce71d66b4fca4719f64ea83f654732b..bceb010e1ee47c31a8109031f90491bcd5b818d0 100644
--- a/doc/sphinx/lbmpy.bib
+++ b/doc/sphinx/lbmpy.bib
@@ -41,6 +41,13 @@
   publisher={APS}
 }
 
+@phdthesis{schiller2008thermal,
+  title={Thermal fluctuations and boundary conditions in the lattice Boltzmann method},
+  author={Schiller, Ulf Daniel},
+  year={2008},
+  school={Johannes Gutenberg Universit{\"a}t Mainz}
+}
+
 
 @article{Wohrwag2018,
 archivePrefix = {arXiv},
diff --git a/lbmpy/forcemodels.py b/lbmpy/forcemodels.py
index 0b96e844434556d409471b10c0f15f1998bfb8b4..5b4a0b19b0db656b6c8c070f8d4eec54ca24a116 100644
--- a/lbmpy/forcemodels.py
+++ b/lbmpy/forcemodels.py
@@ -201,6 +201,7 @@ class Buick:
         simple = Simple(self._force)
 
         shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
+        assert len(set(lb_method.relaxation_rates)) == 1, "Buick only works for SRT"
         correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
         return [correction_factor * t for t in simple(lb_method)]
 
diff --git a/lbmpy_tests/test_entropic.py b/lbmpy_tests/test_entropic.py
index 284570c1d4c874f8d11737a70346ade24ce94eba..47f6e841288e93f26037f5a38a29d2895f609de8 100644
--- a/lbmpy_tests/test_entropic.py
+++ b/lbmpy_tests/test_entropic.py
@@ -10,13 +10,13 @@ from lbmpy.stencils import get_stencil
 def test_entropic_methods():
     sc_kbc = create_lid_driven_cavity((20, 20), method='trt-kbc-n4', relaxation_rate=1.9999,
                                       entropic_newton_iterations=3, entropic=True, compressible=True,
-                                      force=(-1e-10, 0))
+                                      force=(-1e-10, 0), force_model="luo")
 
     sc_srt = create_lid_driven_cavity((40, 40), relaxation_rate=1.9999, lid_velocity=0.05, compressible=True,
-                                      force=(-1e-10, 0))
+                                      force=(-1e-10, 0), force_model="luo")
 
     sc_entropic = create_lid_driven_cavity((40, 40), method='entropic-srt', relaxation_rate=1.9999,
-                                           lid_velocity=0.05, compressible=True, force=(-1e-10, 0))
+                                           lid_velocity=0.05, compressible=True, force=(-1e-10, 0), force_model="luo")
 
     sc_srt.run(1000)
     sc_kbc.run(1000)
diff --git a/lbmpy_tests/test_fluctuation.py b/lbmpy_tests/test_fluctuation.py
index 63078aa12ef5373ae5d350f17c6143894c9a5f33..ded0db9149789bc95af08829a7e68b171752252e 100644
--- a/lbmpy_tests/test_fluctuation.py
+++ b/lbmpy_tests/test_fluctuation.py
@@ -5,7 +5,7 @@ from lbmpy.scenarios import create_channel
 
 def test_fluctuating_generation_pipeline():
     ch = create_channel((10, 10), stencil='D2Q9', method='mrt', weighted=True, relaxation_rates=[1.5] * 5, force=1e-5,
-                        fluctuating={'temperature': 1e-9}, kernel_params={'time_step': 1, 'seed': 312},
+                        force_model='luo', fluctuating={'temperature': 1e-9}, kernel_params={'time_step': 1, 'seed': 312},
                         optimization={'cse_global': True})
 
     ch.run(10)
diff --git a/lbmpy_tests/test_force.py b/lbmpy_tests/test_force.py
index 86d973f48b8d667a6f5dd3cc3dc194536f1a9db0..b42c03ebb214da089924a58073535b51060fa6d7 100644
--- a/lbmpy_tests/test_force.py
+++ b/lbmpy_tests/test_force.py
@@ -1,9 +1,11 @@
 from pystencils.session import *
 from lbmpy.session import *
 from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
-import pytest
 import lbmpy.forcemodels
 
+import pytest
+from contextlib import ExitStack as does_not_raise
+
 
 force_models = [fm.lower() for fm in dir(lbmpy.forcemodels) if fm[0].isupper()]
 
@@ -25,14 +27,22 @@ def test_total_momentum(method, force_model, omega):
     ρ = dh.add_array('rho')
     u = dh.add_array('u', values_per_cell=dh.dim)
 
-    collision = create_lb_update_rule(method=method,
-                                      stencil=stencil,
-                                      relaxation_rate=omega, 
-                                      compressible=True,
-                                      force_model=force_model, 
-                                      force=F,
-                                      kernel_type='collide_only',
-                                      optimization={'symbolic_field': src})
+    expectation = does_not_raise()
+    skip = False
+    if force_model in ['guo', 'buick'] and method != 'srt':
+        expectation = pytest.raises(AssertionError)
+        skip = True
+    with expectation:
+        collision = create_lb_update_rule(method=method,
+                                          stencil=stencil,
+                                          relaxation_rate=omega, 
+                                          compressible=True,
+                                          force_model=force_model, 
+                                          force=F,
+                                          kernel_type='collide_only',
+                                          optimization={'symbolic_field': src})
+    if skip:
+        return
 
     stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
                                                    {'density': ρ, 'velocity': u})
diff --git a/lbmpy_tests/test_split_optimization.py b/lbmpy_tests/test_split_optimization.py
index 821d6868f5dfbddc7e09aeb29a12643729c8133f..2e0b744191eb4c8c477f2ccb6617807d604ce1bf 100644
--- a/lbmpy_tests/test_split_optimization.py
+++ b/lbmpy_tests/test_split_optimization.py
@@ -16,6 +16,7 @@ def test_split_number_of_operations():
                 common_params = {'stencil': stencil,
                                  'method': method,
                                  'compressible': compressible,
+                                 'force_model': 'luo',
                                  'force': (1e-6, 1e-5, 1e-7)
                                  }
                 ast_with_splitting = create_lb_ast(optimization={'split': True}, **common_params)
diff --git a/lbmpy_tests/test_srt_trt_simplifications.py b/lbmpy_tests/test_srt_trt_simplifications.py
index 0d154edab4b016b2510cb04770d7adcd8bc8fdb6..eb799d4484ab538eaac3f6d2a26c67b2ca0539e6 100644
--- a/lbmpy_tests/test_srt_trt_simplifications.py
+++ b/lbmpy_tests/test_srt_trt_simplifications.py
@@ -4,7 +4,7 @@ known acceptable values.
 """
 import sympy as sp
 
-from lbmpy.forcemodels import Guo
+from lbmpy.forcemodels import Luo
 from lbmpy.methods import create_srt, create_trt, create_trt_with_magic_number
 from lbmpy.methods.momentbasedsimplifications import cse_in_opposing_directions
 from lbmpy.simplificationfactory import create_simplification_strategy
@@ -55,13 +55,13 @@ def test_simplifications_trt_d2q9_compressible():
 
 def test_simplifications_trt_d3q19_force_incompressible():
     o1, o2 = sp.symbols("omega_1 omega_2")
-    force_model = Guo([sp.Rational(1, 3), sp.Rational(1, 2), sp.Rational(1, 5)])
+    force_model = Luo([sp.Rational(1, 3), sp.Rational(1, 2), sp.Rational(1, 5)])
     method = create_trt(get_stencil("D3Q19"), o1, o2, compressible=False, force_model=force_model)
     check_method(method, [268, 281, 0], [241, 175, 1])
 
 
 def test_simplifications_trt_d3q19_force_compressible():
     o1, o2 = sp.symbols("omega_1 omega_2")
-    force_model = Guo([sp.Rational(1, 3), sp.Rational(1, 2), sp.Rational(1, 5)])
+    force_model = Luo([sp.Rational(1, 3), sp.Rational(1, 2), sp.Rational(1, 5)])
     method = create_trt_with_magic_number(get_stencil("D3Q19"), o1, compressible=False, force_model=force_model)
     check_method(method, [270, 283, 1], [243, 177, 1])