From c5b73cbb755183492ef7cf07c32808a391b0432a Mon Sep 17 00:00:00 2001
From: Helen Schottenhamml <helen.schottenhamml@fau.de>
Date: Fri, 26 Jan 2024 15:10:17 +0100
Subject: [PATCH] Rework SGS models

---
 ...utorial_modifying_method_smagorinsky.ipynb |   4 +-
 doc/sphinx/lbmpy.bib                          |  12 ++
 src/lbmpy/__init__.py                         |   3 +-
 src/lbmpy/creationfunctions.py                |  54 +++++----
 src/lbmpy/db.py                               |   4 +-
 src/lbmpy/enums.py                            |  17 ++-
 src/lbmpy/turbulence_models.py                | 104 ++++++++++++++++--
 tests/test_json_serializer.py                 |   8 +-
 tests/test_lbstep.py                          |   9 +-
 tests/test_smagorinsky.py                     |  10 +-
 10 files changed, 183 insertions(+), 42 deletions(-)

diff --git a/doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb b/doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb
index 574932fb..f3f3fe52 100644
--- a/doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb
+++ b/doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb
@@ -445,10 +445,10 @@
     }
    ],
    "source": [
-    "lbm_conifg = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, force=(1e-6, 0),\n",
+    "lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, force=(1e-6, 0),\n",
     "                       force_model=ForceModel.LUO, relaxation_rates=[ω, 1.9, 1.9, 1.9])\n",
     "\n",
-    "method = create_lb_method(lbm_config=lbm_conifg)\n",
+    "method = create_lb_method(lbm_config=lbm_config)\n",
     "method"
    ]
   },
diff --git a/doc/sphinx/lbmpy.bib b/doc/sphinx/lbmpy.bib
index 551f3f5f..efce2c1a 100644
--- a/doc/sphinx/lbmpy.bib
+++ b/doc/sphinx/lbmpy.bib
@@ -276,4 +276,16 @@ journal = {Communications in Computational Physics}
     url = {https://doi.org/10.1063/1.1399290},
 }
 
+@Article{rozema15,
+  doi       = {10.1063/1.4928700},
+  year      = {2015},
+  month     = {aug},
+  publisher = {AIP Publishing},
+  volume    = {27},
+  number    = {8},
+  author    = {Wybe Rozema and Hyun J. Bae and Parviz Moin and Roel Verstappen},
+  title     = {Minimum-dissipation models for large-eddy simulation},
+  journal   = {Physics of Fluids}
+}
+
 @Comment{jabref-meta: databaseType:bibtex;}
diff --git a/src/lbmpy/__init__.py b/src/lbmpy/__init__.py
index 542e736a..2c41a0df 100644
--- a/src/lbmpy/__init__.py
+++ b/src/lbmpy/__init__.py
@@ -7,7 +7,7 @@ from .creationfunctions import (
     LBMConfig,
     LBMOptimisation,
 )
-from .enums import Stencil, Method, ForceModel, CollisionSpace
+from .enums import Stencil, Method, ForceModel, CollisionSpace, SubgridScaleModel
 from .lbstep import LatticeBoltzmannStep
 from .macroscopic_value_kernels import (
     pdf_initialization_assignments,
@@ -38,6 +38,7 @@ __all__ = [
     "Method",
     "ForceModel",
     "CollisionSpace",
+    "SubgridScaleModel",
     "LatticeBoltzmannStep",
     "pdf_initialization_assignments",
     "macroscopic_values_getter",
diff --git a/src/lbmpy/creationfunctions.py b/src/lbmpy/creationfunctions.py
index b0497e6b..85cef388 100644
--- a/src/lbmpy/creationfunctions.py
+++ b/src/lbmpy/creationfunctions.py
@@ -20,21 +20,21 @@ and belongs to pystencils, not lbmpy. This can be found in the pystencils module
 of the generated code is specified.
 
 1. *Method*:
-         the method defines the collision process. Currently there are two big categories:
+         the method defines the collision process. Currently, there are two big categories:
          moment and cumulant based methods. A method defines how each moment or cumulant is relaxed by
          storing the equilibrium value and the relaxation rate for each moment/cumulant.
 2. *Collision/Update Rule*:
          Methods can generate a "collision rule" which is an equation collection that define the
          post collision values as a function of the pre-collision values. On these equation collection
          simplifications are applied to reduce the number of floating point operations.
-         At this stage an entropic optimization step can also be added to determine one relaxation rate by an
+         At this stage an entropic optimisation step can also be added to determine one relaxation rate by an
          entropy condition.
          Then a streaming rule is added which transforms the collision rule into an update rule.
          The streaming step depends on the pdf storage (source/destination, AABB pattern, EsoTwist).
          Currently only the simple source/destination  pattern is supported.
 3. *AST*:
         The abstract syntax tree describes the structure of the kernel, including loops and conditionals.
-        The ast can be modified e.g. to add OpenMP pragmas, reorder loops or apply other optimizations.
+        The ast can be modified, e.g., to add OpenMP pragmas, reorder loops or apply other optimisations.
 4. *Function*:
         This step compiles the AST into an executable function, either for CPU or GPUs. This function
         behaves like a normal Python function and runs one LBM time step.
@@ -63,7 +63,7 @@ import pystencils.astnodes
 import sympy as sp
 import sympy.core.numbers
 
-from lbmpy.enums import Stencil, Method, ForceModel, CollisionSpace
+from lbmpy.enums import Stencil, Method, ForceModel, CollisionSpace, SubgridScaleModel
 import lbmpy.forcemodels as forcemodels
 from lbmpy.fieldaccess import CollideOnlyInplaceAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor
 from lbmpy.fluctuatinglb import add_fluctuations_to_collision_rule
@@ -78,7 +78,7 @@ from lbmpy.methods.momentbased.entropic import add_entropy_condition, add_iterat
 from lbmpy.relaxationrates import relaxation_rate_from_magic_number
 from lbmpy.simplificationfactory import create_simplification_strategy
 from lbmpy.stencils import LBStencil
-from lbmpy.turbulence_models import add_smagorinsky_model
+from lbmpy.turbulence_models import add_sgs_model
 from lbmpy.updatekernels import create_lbm_kernel, create_stream_pull_with_output_kernel
 from lbmpy.advanced_streaming.utility import Timestep, get_accessor
 from pystencils import CreateKernelConfig, create_kernel
@@ -145,7 +145,7 @@ class LBMConfig:
     """
     delta_equilibrium: bool = None
     """
-    Determines whether or not the (continuous or discrete, see `continuous_equilibrium`) maxwellian equilibrium is
+    Determines whether or not the (continuous or discrete, see `continuous_equilibrium`) Maxwellian equilibrium is
     expressed in its absolute form, or only by its deviation from the rest state (typically given by the reference
     density and zero velocity). This parameter is only effective if `zero_centered` is set to `True`. Then, if
     `delta_equilibrium` is `False`, the rest state must be reintroduced to the populations during collision. Otherwise,
@@ -175,7 +175,7 @@ class LBMConfig:
     """
     A list of lists of modes, grouped by common relaxation times. This is usually used in
     conjunction with `lbmpy.methods.default_moment_sets.mrt_orthogonal_modes_literature`.
-    If this argument is not provided, Gram-Schmidt orthogonalization of the default modes is performed.
+    If this argument is not provided, Gram-Schmidt orthogonalisation of the default modes is performed.
     """
 
     force_model: Union[lbmpy.forcemodels.AbstractForceModel, ForceModel] = None
@@ -239,12 +239,17 @@ class LBMConfig:
     omega_output_field: Field = None
     """
     A pystencils Field can be passed here, where the calculated free relaxation rate of
-    an entropic or Smagorinsky method is written to
+    an entropic or subgrid-scale method is written to
     """
-    smagorinsky: Union[float, bool] = False
+    eddy_viscosity_field: Field = None
     """
-    set to Smagorinsky constant to activate turbulence model, ``omega_output_field`` can be set to
-    write out adapted relaxation rates. If set to `True`, 0.12 is used as default smagorinsky constant.
+    A pystencils Field can be passed here, where the eddy-viscosity of a subgrid-scale model is written.
+    """
+    subgrid_scale_model: Union[SubgridScaleModel, tuple[SubgridScaleModel, float], tuple[SubgridScaleModel, int]] = None
+    """
+    Choose a subgrid-scale model (SGS) for large-eddy simulations. ``omega_output_field`` can be set to
+    write out adapted relaxation rates. Either provide just the SGS and use the default model constants or provide a
+    tuple of the SGS and its corresponding model constant.
     """
     cassons: CassonsParameters = False
     """
@@ -702,8 +707,8 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
                                                      limiter=cumulant_limiter)
 
     if lbm_config.entropic:
-        if lbm_config.smagorinsky or lbm_config.cassons:
-            raise ValueError("Choose either entropic, smagorinsky or cassons")
+        if lbm_config.subgrid_scale_model or lbm_config.cassons:
+            raise ValueError("Choose either entropic, subgrid-scale or cassons")
         if lbm_config.entropic_newton_iterations:
             if isinstance(lbm_config.entropic_newton_iterations, bool):
                 iterations = 3
@@ -714,14 +719,23 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
         else:
             collision_rule = add_entropy_condition(collision_rule, omega_output_field=lbm_config.omega_output_field)
 
-    elif lbm_config.smagorinsky:
+    elif lbm_config.subgrid_scale_model:
         if lbm_config.cassons:
-            raise ValueError("Cassons model can not be combined with Smagorinsky model")
-        smagorinsky_constant = 0.12 if lbm_config.smagorinsky is True else lbm_config.smagorinsky
-        collision_rule = add_smagorinsky_model(collision_rule, smagorinsky_constant,
-                                               omega_output_field=lbm_config.omega_output_field)
-        if 'split_groups' in collision_rule.simplification_hints:
-            collision_rule.simplification_hints['split_groups'][0].append(sp.Symbol("smagorinsky_omega"))
+            raise ValueError("Cassons model can not be combined with a subgrid-scale model")
+
+        model_constant = None
+        sgs_model = lbm_config.subgrid_scale_model
+
+        if isinstance(lbm_config.subgrid_scale_model, tuple):
+            sgs_model = lbm_config.subgrid_scale_model[0]
+            model_constant = lbm_config.subgrid_scale_model[1]
+
+        collision_rule = add_sgs_model(collision_rule=collision_rule, subgrid_scale_model=sgs_model,
+                                       model_constant=model_constant, omega_output_field=lbm_config.omega_output_field,
+                                       eddy_viscosity_field=lbm_config.eddy_viscosity_field)
+
+    if 'split_groups' in collision_rule.simplification_hints:
+        collision_rule.simplification_hints['split_groups'][0].append(sp.Symbol("sgs_omega"))
 
     elif lbm_config.cassons:
         collision_rule = add_cassons_model(collision_rule, parameter=lbm_config.cassons,
diff --git a/src/lbmpy/db.py b/src/lbmpy/db.py
index bca65fb1..cd7983cb 100644
--- a/src/lbmpy/db.py
+++ b/src/lbmpy/db.py
@@ -4,7 +4,7 @@ import inspect
 
 from pystencils.runhelper.db import PystencilsJsonEncoder
 from pystencils.simp import SimplificationStrategy
-from lbmpy import LBStencil, Method, CollisionSpace
+from lbmpy import LBStencil, Method, CollisionSpace, SubgridScaleModel
 from lbmpy.creationfunctions import LBMConfig, LBMOptimisation
 from lbmpy.methods import CollisionSpaceInfo
 from lbmpy.forcemodels import AbstractForceModel
@@ -16,7 +16,7 @@ class LbmpyJsonEncoder(PystencilsJsonEncoder):
     def default(self, obj):
         if isinstance(obj, (LBMConfig, LBMOptimisation, CollisionSpaceInfo, CassonsParameters)):
             return obj.__dict__
-        if isinstance(obj, (LBStencil, Method, CollisionSpace)):
+        if isinstance(obj, (LBStencil, Method, CollisionSpace, SubgridScaleModel)):
             return obj.name
         if isinstance(obj, AbstractForceModel):
             return obj.__class__.__name__
diff --git a/src/lbmpy/enums.py b/src/lbmpy/enums.py
index 3eee11f0..1f017050 100644
--- a/src/lbmpy/enums.py
+++ b/src/lbmpy/enums.py
@@ -219,5 +219,20 @@ class ForceModel(Enum):
     """
     CENTRALMOMENT = auto()
     """
-    See :class:`lbmpy.forcemodels`
+    See :class:`lbmpy.forcemodels.CentralMoment`
+    """
+
+
+class SubgridScaleModel(Enum):
+    """
+    The SubgridScaleModel enumeration defines which subgrid-scale model (SGS) is used to perform
+    Large-Eddy-Simulations (LES).
+    """
+    SMAGORINSKY = auto()
+    """
+    See :func:`lbmpy.turbulence_models.add_smagorinsky_model`
+    """
+    QR = auto()
+    """
+    See :func:`lbmpy.turbulence_models.add_qr_model`
     """
diff --git a/src/lbmpy/turbulence_models.py b/src/lbmpy/turbulence_models.py
index 8f98385a..84cd6bef 100644
--- a/src/lbmpy/turbulence_models.py
+++ b/src/lbmpy/turbulence_models.py
@@ -1,12 +1,29 @@
 import sympy as sp
 
-from lbmpy.relaxationrates import get_shear_relaxation_rate
+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_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_field=None):
-    r""" Adds a smagorinsky model to a lattice Boltzmann collision rule. To add the Smagorinsky model to a LB scheme
+
+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`.
@@ -34,7 +51,7 @@ def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_fie
     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 "
+        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
@@ -42,18 +59,91 @@ def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_fie
     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("smagorinsky_omega")
+    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, 1 / omega_s),
+    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,
-                      2 / (tau_0 + sp.sqrt(18 * smagorinsky_constant ** 2 * second_order_neq_moments + tau_0 ** 2)))]
+                      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))
 
diff --git a/tests/test_json_serializer.py b/tests/test_json_serializer.py
index 3429f12b..aee893b5 100644
--- a/tests/test_json_serializer.py
+++ b/tests/test_json_serializer.py
@@ -5,7 +5,7 @@ Test the lbmpy-specific JSON encoder and serializer as used in the Database clas
 import tempfile
 import pystencils as ps
 
-from lbmpy import Stencil, Method, ForceModel
+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
@@ -26,9 +26,9 @@ def test_json_serializer():
 
     # create dummy lbmpy config
     lbm_config = LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, force_model=ForceModel.GUO,
-                           compressible=True, relaxation_rate=1.999, smagorinsky=True, galilean_correction=True,
-                           cassons=cassons_params, density_input=density, kernel_type=StreamPullTwoFieldsAccessor,
-                           timestep=Timestep.BOTH)
+                           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)
diff --git a/tests/test_lbstep.py b/tests/test_lbstep.py
index 629d106e..73d36a25 100644
--- a/tests/test_lbstep.py
+++ b/tests/test_lbstep.py
@@ -4,6 +4,7 @@ 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 cupy
@@ -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)
 
 
diff --git a/tests/test_smagorinsky.py b/tests/test_smagorinsky.py
index bb5fc99b..de71f39e 100644
--- a/tests/test_smagorinsky.py
+++ b/tests/test_smagorinsky.py
@@ -1,19 +1,21 @@
 import sympy as sp
 from pystencils.simp.subexpression_insertion import insert_constants
 
-from lbmpy import create_lb_collision_rule, LBMConfig, LBStencil, Stencil, Method
+from lbmpy import create_lb_collision_rule, LBMConfig, LBStencil, Stencil, Method, SubgridScaleModel
 
 
 def test_smagorinsky_with_constant_omega():
     stencil = LBStencil(Stencil.D2Q9)
 
-    config = LBMConfig(stencil=stencil, method=Method.SRT, smagorinsky=True, relaxation_rate=sp.Symbol("omega"))
+    config = LBMConfig(stencil=stencil, method=Method.SRT, subgrid_scale_model=SubgridScaleModel.SMAGORINSKY,
+                       relaxation_rate=sp.Symbol("omega"))
     collision_rule = create_lb_collision_rule(lbm_config=config)
 
-    config = LBMConfig(stencil=stencil, method=Method.SRT, smagorinsky=True, relaxation_rate=1.5)
+    config = LBMConfig(stencil=stencil, method=Method.SRT, subgrid_scale_model=SubgridScaleModel.SMAGORINSKY,
+                       relaxation_rate=1.5)
     collision_rule2 = create_lb_collision_rule(lbm_config=config)
 
     collision_rule = collision_rule.subs({sp.Symbol("omega"): 1.5})
     collision_rule = insert_constants(collision_rule)
 
-    assert collision_rule == collision_rule2
\ No newline at end of file
+    assert collision_rule == collision_rule2
-- 
GitLab