diff --git a/doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb b/doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb index 574932fbe8b2a9b1a4a592e82a0de36772d4cee3..f3f3fe52ed603056e73c7b4df3bbbd2d1b7d3e85 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 551f3f5f9c2c23b4565e13f8bb514410fd1d298c..efce2c1ac21aaf48d1c7b2298e508186703e1284 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 542e736af303e868cc2f1c131789f4bc010e62ff..2c41a0dfa569c5d08a2b8edad816f3cff57c3ab8 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 b0497e6bbfe91b17504cddd1fa7531acc7e5e0d4..85cef3886dd4f32825ebb8810e6dc8c699be8e8b 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 bca65fb182a125600b2e309dd2d5784f6d4d8c47..cd7983cba7d95e8fff2274954f782e9a518668ca 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 3eee11f027f1d01db1f796a7390daf2162a89bcf..1f017050d37b9ed33789d08eb2da8584a4fd12ac 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 8f98385a7e95244232247ed964b98cbe51f51097..84cd6befd1e766425b816059f26f13d4d6403281 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 3429f12bb6ba0eff62e76e127707f08b290d8547..aee893b529e385e3455e889f9d80cf9a7e1374b3 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 629d106e3379a9c0996374f47bf46d5d817d569e..73d36a2565b129893b98caaa7120d5e334aa340e 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 bb5fc99bbd8e9a7f970c776c9c64d7400a9e73a2..de71f39ec1105dbe774d5f4f81e75e9f9e79ab68 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