Write eddy-viscosity directly in SGS calculation.

src/lbmpy/creationfunctions.py

@@ -241,6 +241,10 @@ class LBMConfig:
     A pystencils Field can be passed here, where the calculated free relaxation rate of
     an entropic or subgrid-scale method is written to
     """
+    eddy_viscosity_field: Field = None
+    """
+    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
@@ -727,7 +731,8 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
             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)
+                                       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"))

src/lbmpy/turbulence_models.py

@@ -8,20 +8,21 @@ 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):
+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)
+                                     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)
+                            omega_output_field=omega_output_field, eddy_viscosity_field=eddy_viscosity_field)
 
 
-def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_field=None):
+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
@@ -71,13 +72,19 @@ def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_fie
     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):
+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
     """
@@ -134,6 +141,9 @@ def add_qr_model(collision_rule, qr_constant, omega_output_field=None):
     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))