Commit 63b38380 by Markus Holzer Committed by Michael Kuron

### Fix minor issues and remove depricated usage of cumulant LB method

parent 7df9ca58
 %% Cell type:code id: tags:  python from lbmpy.session import *  %% Cell type:markdown id: tags: # Tutorial 03: Defining LB methods in *lbmpy* ## A) General Form %% Cell type:markdown id: tags: The lattice Boltzmann equation in its most general form is: $$f_q(\mathbf{x} + \mathbf{c}_q \delta t, t+\delta t) = K\left( f_q(\mathbf{x}, t) \right)$$ with a discrete velocity set $\mathbf{c}_q$ (stencil) and a generic collision operator $K$. So a lattice Boltzmann method can be fully defined by picking a stencil and a collision operator. The collision operator $K$ has the following structure: - Transformation of particle distribution function $f$ into collision space. This transformation has to be invertible and may be nonlinear. - The collision operation is an convex combination of the pdf representation in collision space $c$ and some equilibrium vector $c^{(eq)}$. This equilibrium can also be defined in physical space, then $c^{(eq)} = C( f^{(eq)} )$. The convex combination is done elementwise using a diagonal matrix $S$ where the diagonal entries are the relaxation rates. - After collision, the collided state $c'$ is transformed back into physical space ![](../img/collision.svg) The full collision operator is: $$K(f) = C^{-1}\left( (I-S)C(f) + SC(f^{(eq}) \right)$$ or $$K(f) = C^{-1}\left( C(f) - S (C(f) - C(f^{(eq})) \right)$$ %% Cell type:markdown id: tags: ## B) Moment-based relaxation The most commonly used LBM collision operator is the multi relaxation time (MRT) collision. In MRT methods the collision space is spanned by moments of the distribution function. This is a very natural approach, since the pdf moments are the quantities that go into the Chapman Enskog analysis that is used to show that LB methods can solve the Navier Stokes equations. Also the lower order moments correspond to the macroscopic quantities of interest (density/pressure, velocity, shear rates, heat flux). Furthermore the transformation to collision space is linear in this case, simplifying the collision equations: $$K(f) = C^{-1}\left( C(f) - S (C(f) - C(f^{(eq})) \right)$$ $$K(f) = f - \underbrace{ C^{-1}SC}_{A}(f - f^{(eq)})$$ in *lbmpy* the following formulation is used, since it is more natural to define the equilibrium in moment-space instead of physical space: $$K(f) = f - C^{-1}S(Cf - c^{(eq)})$$ %% Cell type:markdown id: tags: ### Use a pre-defined method Lets create a moment-based method in *lbmpy* and see how the moment transformation $C$ and the relaxation rates that comprise the diagonal matrix $S$ can be defined. We start with a function that creates a basic MRT model. Don't use this for real simulations, there orthogonalized MRT methods should be used, as discussed in the next section. %% Cell type:code id: tags:  python from lbmpy.creationfunctions import create_lb_method method = create_lb_method(stencil='D2Q9', method='mrt_raw') # check also method='srt', 'trt', 'mrt' method  %% Output %% Cell type:markdown id: tags: The first column labeled "Moment" defines the collision space and thus the transformation matrix $C$. The remaining columns specify the equilibrium vector in moment space $c^{(eq)}$ and the corresponding relaxation rate. Each row of the "Moment" column defines one row of $C$. In the next cells this matrix and the discrete velocity set (stencil) of our method are shown. Check for example the second last row of the table $x^2 y$: In the corresponding second last row of the moment matrix $C$ where each column stands for a lattice velocity (for ordering visualized stencil below) and each entry is the expression $x^2 y$ where $x$ and $y$ are the components of the lattice velocity. In general the transformation matrix $C_{iq}$ is defined as; $$c_i = C_{iq} f_q = \sum_q m_i(c_q)$$ where $m_i(c_q)$ is the $i$'th moment polynomial where $x$ and $y$ are substituted with the components of the $q$'th lattice velocity %% Cell type:code id: tags:  python # Transformation matrix C method.moment_matrix  %% Output $\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\0 & 0 & 0 & -1 & 1 & -1 & 1 & -1 & 1\\0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1\\0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1\\0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1\\0 & 0 & 0 & 0 & 0 & -1 & 1 & 1 & -1\\0 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1\\0 & 0 & 0 & 0 & 0 & -1 & 1 & -1 & 1\\0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\end{matrix}\right]$ ⎡1 1 1 1 1 1 1 1 1 ⎤ ⎢ ⎥ ⎢0 0 0 -1 1 -1 1 -1 1 ⎥ ⎢ ⎥ ⎢0 1 -1 0 0 1 1 -1 -1⎥ ⎢ ⎥ ⎢0 0 0 1 1 1 1 1 1 ⎥ ⎢ ⎥ ⎢0 1 1 0 0 1 1 1 1 ⎥ ⎢ ⎥ ⎢0 0 0 0 0 -1 1 1 -1⎥ ⎢ ⎥ ⎢0 0 0 0 0 1 1 -1 -1⎥ ⎢ ⎥ ⎢0 0 0 0 0 -1 1 -1 1 ⎥ ⎢ ⎥ ⎣0 0 0 0 0 1 1 1 1 ⎦ %% Cell type:code id: tags:  python ps.stencil.plot(method.stencil) method.stencil  %% Output $\displaystyle \left( \left( 0, \ 0\right), \ \left( 0, \ 1\right), \ \left( 0, \ -1\right), \ \left( -1, \ 0\right), \ \left( 1, \ 0\right), \ \left( -1, \ 1\right), \ \left( 1, \ 1\right), \ \left( -1, \ -1\right), \ \left( 1, \ -1\right)\right)$ ((0, 0), (0, 1), (0, -1), (-1, 0), (1, 0), (-1, 1), (1, 1), (-1, -1), (1, -1)) %% Cell type:markdown id: tags: ### Orthogonal MRTs For a real MRT method, the moments should be orthogonalized. One can either orthogonalize using the standard scalar product or a scalar product that is weighted with the lattice weights. If unsure, use the weighted version. The next cell shows how to get both orthogonalized MRT versions in lbmpy. %% Cell type:code id: tags:  python weighted_ortho_mrt = create_lb_method(stencil="D2Q9", method="mrt", weighted=True) weighted_ortho_mrt  %% Output %% Cell type:code id: tags:  python ortho_mrt = create_lb_method(stencil="D2Q9", method="mrt", weighted=False) ortho_mrt  %% Output %% Cell type:markdown id: tags: One can check if a method is orthogonalized: %% Cell type:code id: tags:  python ortho_mrt.is_orthogonal, weighted_ortho_mrt.is_weighted_orthogonal  %% Output (True, True) %% Cell type:markdown id: tags: ### Define custom MRT method To choose custom values for the left moment column one can pass a nested list of moments. Moments that should be relaxed with the same paramter are grouped together. *lbmpy* also comes with a few templates for this list taken from literature: %% Cell type:code id: tags:  python from lbmpy.methods import mrt_orthogonal_modes_literature from lbmpy.stencils import get_stencil from lbmpy.moments import MOMENT_SYMBOLS x, y, z = MOMENT_SYMBOLS moments = mrt_orthogonal_modes_literature(get_stencil("D2Q9"), is_weighted=True, is_cumulant=False) moments = mrt_orthogonal_modes_literature(get_stencil("D2Q9"), is_weighted=True) moments  %% Output $\displaystyle \left[ \left[ 1\right], \ \left[ x, \ y\right], \ \left[ 3 x^{2} + 3 y^{2} - 2, \ x^{2} - y^{2}, \ x y\right], \ \left[ x \left(3 x^{2} + 3 y^{2} - 4\right), \ y \left(3 x^{2} + 3 y^{2} - 4\right)\right], \ \left[ - 15 x^{2} - 15 y^{2} + 9 \left(x^{2} + y^{2}\right)^{2} + 2\right]\right]$ ⎡ ⎢ ⎡ 2 2 2 2 ⎤ ⎡ ⎛ 2 2 ⎞ ⎛ 2 ⎣[1], [x, y], ⎣3⋅x + 3⋅y - 2, x - y , x⋅y⎦, ⎣x⋅⎝3⋅x + 3⋅y - 4⎠, y⋅⎝3⋅x + ⎡ 2 ⎤⎤ 2 ⎞⎤ ⎢ 2 2 ⎛ 2 2⎞ ⎥⎥ 3⋅y - 4⎠⎦, ⎣- 15⋅x - 15⋅y + 9⋅⎝x + y ⎠ + 2⎦⎦ %% Cell type:markdown id: tags: This nested moment list can be passed to create_lb_method: %% Cell type:code id: tags:  python create_lb_method(stencil="D2Q9", method="mrt", nested_moments=moments)  %% Output %% Cell type:markdown id: tags: If one needs to also specify custom equilibrium moments the following approach can be used %% Cell type:code id: tags:  python rho = sp.symbols("rho") u = sp.symbols("u_:3") omega = sp.symbols("omega_:4") method_table = [ # Conserved moments (1, rho, 0 ), (x, u[0], 0 ), (y, u[1], 0 ), # Shear moments (x*y, u[0]*u[1], omega[0]), (x**2-y**2, u[0]**2 - u[1]**2, omega[0]), (x**2+y**2, 2*rho/3 + u[0]**2 + u[1]**2, omega[1]), # Higher order (x * y**2, u[0]/3, omega[2]), (x**2 * y, u[1]/3, omega[2]), (x**2 * y**2, rho/9 + u[0]**2/3 + u[1]**2/3, omega[3]), ] method = create_generic_mrt(get_stencil("D2Q9"), method_table) method  %% Output %% Cell type:markdown id: tags: Instead of manually defining all entries in the method table, *lbmpy* has functions to fill the table according to a specific pattern. For example: - for a full stencil (D2Q9, D3Q27) there exist exactly 9 or 27 linearly independent moments. These can either be taken as they are, or orthogonalized using Gram-Schmidt, weighted Gram-Schmidt or a Hermite approach - equilibrium values can be computed from the standard discrete equilibrium of order 1,2 or 3. Alternatively they can also be computed as continuous moments of a Maxwellian distribution One option is to start with one of *lbmpy*'s built-in methods and modify it with create_lb_method_from_existing. In the next cell we fix the fourth order relaxation rate to a constant, by writing a function that defines how to alter each row of the collision table. This is for demonstration only, of course we could have done it right away when passing in the collision table. %% Cell type:code id: tags:  python def modification_func(moment, eq, rate): if rate == omega[3]: return moment, eq, 1.0 return moment, eq, rate method = create_lb_method_from_existing(method, modification_func) method  %% Output %% Cell type:markdown id: tags: Our customized method can be directly passed into one of the scenarios. We can for example set up a channel flow with it. Since we used symbols as relaxation rates, we have to pass them in as kernel_params. %% Cell type:code id: tags:  python ch = create_channel(domain_size=(100, 30), lb_method=method, u_max=0.05, kernel_params={'omega_0': 1.8, 'omega_1': 1.4, 'omega_2': 1.5}) ch.run(500) plt.figure(dpi=200) plt.vector_field(ch.velocity[:, :]);  %% Output %% Cell type:markdown id: tags: #### Bonus: Automatic analysis Above we have created a non-orthogonal MRT, where the shear viscosity and bulk viscosity can be independently controlled. For moment-based methods, *lbmpy* also offers an automatic Chapman Enskog analysis that can find the relation between viscosity and relaxation rate(s): %% Cell type:code id: tags:  python from lbmpy.chapman_enskog import ChapmanEnskogAnalysis analysis = ChapmanEnskogAnalysis(method) analysis.get_dynamic_viscosity()  %% Output $\displaystyle - \frac{\omega_{0} - 2}{6 \omega_{0}}$ -(ω₀ - 2) ────────── 6⋅ω₀ %% Cell type:code id: tags:  python analysis.get_bulk_viscosity()  %% Output $\displaystyle - \frac{1}{9} - \frac{1}{3 \omega_{1}} + \frac{5}{9 \omega_{0}}$ 1 1 5 - ─ - ──── + ──── 9 3⋅ω₁ 9⋅ω₀ ... ...
This diff is collapsed.
This diff is collapsed.
 ... ... @@ -53,7 +53,6 @@ General: - maxwellian_moments=True: way to compute equilibrium moments/cumulants, if False the standard discretized LBM equilibrium is used, otherwise the equilibrium moments are computed from the continuous Maxwellian. This makes only a difference if sparse stencils are used e.g. D2Q9 and D3Q27 are not affected, D319 and DQ15 are - cumulant=False: use cumulants instead of moments (deprecated: use method=cumulant directly) - initial_velocity=None: initial velocity in domain, can either be a tuple (x,y,z) velocity to set a constant velocity everywhere, or a numpy array with the same size of the domain, with a last coordinate of shape dim to set velocities on cell level ... ... @@ -87,8 +86,6 @@ Cumulant methods: - galilean_correction=False: Special correction for D3Q27 cumulant LBMs. For Details see :mod:lbmpy.methods.centeredcumulant.galilean_correction - pre_simplification=True: Simplifications applied during the derivaton of the collision rule for cumulant LBMs For details see :mod:lbmpy.methods.momentbased.moment_transforms. LES methods: ... ... @@ -112,7 +109,9 @@ Simplifications / Transformations: load/store streams and thus speeds up the kernel on most architectures - builtin_periodicity=(False,False,False): instead of handling periodicity by copying ghost layers, the periodicity is built into the kernel. This parameters specifies if the domain is periodic in (x,y,z) direction. Even if the periodicity is built into the kernel, the fields have one ghost layer to be consistent with other functions. periodicity is built into the kernel, the fields have one ghost layer to be consistent with other functions. - pre_simplification=True: Simplifications applied during the derivaton of the collision rule for cumulant LBMs For details see :mod:lbmpy.methods.momentbased.moment_transforms. Field size information: ... ... @@ -190,6 +189,7 @@ For example, to modify the AST one can run:: func = create_lb_function(ast=ast, ...) """ from collections import OrderedDict from copy import copy import sympy as sp ... ... @@ -199,8 +199,9 @@ import lbmpy.methods.centeredcumulant.force_model as cumulant_force_model from lbmpy.fieldaccess import CollideOnlyInplaceAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor from lbmpy.fluctuatinglb import add_fluctuations_to_collision_rule from lbmpy.methods import (create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc) from lbmpy.methods.abstractlbmethod import RelaxationInfo from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod from lbmpy.methods.momentbased.moment_transforms import PdfsToCentralMomentsByShiftMatrix from lbmpy.methods.momentbased.moment_transforms import FastCentralMomentTransform from lbmpy.methods.centeredcumulant.cumulant_transform import CentralMomentsToCumulantsByGeneratingFunc from lbmpy.methods.creationfunctions import ( create_with_monomial_cumulants, create_with_polynomial_cumulants, create_with_default_polynomial_cumulants) ... ... @@ -416,7 +417,6 @@ def create_lb_method(**params): 'equilibrium_order': params['equilibrium_order'], 'force_model': force_model, 'maxwellian_moments': params['maxwellian_moments'], 'cumulant': params['cumulant'], 'c_s_sq': params['c_s_sq'], } ... ... @@ -443,10 +443,6 @@ def create_lb_method(**params): if all(get_order(m) < 2 for m in moments): if params['entropic']: return relaxation_rates[0] elif params['cumulant']: result = relaxation_rates[next_relaxation_rate[0]] next_relaxation_rate[0] += 1 return result else: return 0 res = relaxation_rates[next_relaxation_rate[0]] ... ... @@ -494,11 +490,27 @@ def create_lb_method_from_existing(method, modification_function): modification_function: function receiving (moment, equilibrium_value, relaxation_rate) as arguments, i.e. one row of the relaxation table, returning a modified version """ relaxation_table = (modification_function(m, eq, rr) for m, eq, rr in zip(method.moments, method.moment_equilibrium_values, method.relaxation_rates)) compressible = method.conserved_quantity_computation.compressible cumulant = isinstance(method, CenteredCumulantBasedLbMethod) return create_generic_mrt(method.stencil, relaxation_table, compressible, method.force_model, cumulant) if isinstance(method, CenteredCumulantBasedLbMethod): rr_dict = OrderedDict() relaxation_table = (modification_function(m, eq, rr) for m, eq, rr in zip(method.cumulants, method.cumulant_equilibrium_values, method.relaxation_rates)) for cumulant, eq_value, rr in relaxation_table: cumulant = sp.sympify(cumulant) rr_dict[cumulant] = RelaxationInfo(eq_value, rr) return CenteredCumulantBasedLbMethod(method.stencil, rr_dict, method.conserved_quantity_computation, method.force_model, galilean_correction=method.galilean_correction, central_moment_transform_class=method.central_moment_transform_class, cumulant_transform_class=method.cumulant_transform_class) else: relaxation_table = (modification_function(m, eq, rr) for m, eq, rr in zip(method.moments, method.moment_equilibrium_values, method.relaxation_rates)) return create_generic_mrt(method.stencil, relaxation_table, compressible, method.force_model) # ---------------------------------------------------------------------------------------------------------------------- ... ... @@ -564,11 +576,10 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para 'force_model': 'none', 'force': (0, 0, 0), 'maxwellian_moments': True, 'cumulant': False, # Depricated usage. Cumulant is now an own method 'initial_velocity': None, 'galilean_correction': False, # only available for D3Q27 cumulant methods 'central_moment_transform_class': PdfsToCentralMomentsByShiftMatrix, 'central_moment_transform_class': FastCentralMomentTransform, 'cumulant_transform_class': CentralMomentsToCumulantsByGeneratingFunc, 'entropic': False, ... ...
 ... ... @@ -159,21 +159,45 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod): self.set_first_moment_relaxation_rate(self._force_model.override_momentum_relaxation_rate) self.force_model_rr_override = True @property def central_moment_transform_class(self): self._central_moment_transform_class @property def cumulants(self): return tuple(self._cumulant_to_relaxation_info_dict.keys()) @property def cumulant_equilibrium_values(self): return tuple([e.equilibrium_value for e in self._cumulant_to_relaxation_info_dict.values()]) @property def cumulant_transform_class(self): self._cumulant_transform_class @property def first_order_equilibrium_moment_symbols(self, ): return self._conserved_quantity_computation.first_order_moment_symbols @property def force_model(self): return self._force_model @property def galilean_correction(self): return self._galilean_correction @property def relaxation_info_dict(self): return self._cumulant_to_relaxation_info_dict @property def zeroth_order_equilibrium_moment_symbol(self, ): return self._conserved_quantity_computation.zeroth_order_moment_symbol def relaxation_rates(self): return tuple([e.relaxation_rate for e in self._cumulant_to_relaxation_info_dict.values()]) @property def first_order_equilibrium_moment_symbols(self, ): return self._conserved_quantity_computation.first_order_moment_symbols def zeroth_order_equilibrium_moment_symbol(self, ): return self._conserved_quantity_computation.zeroth_order_moment_symbol def set_zeroth_moment_relaxation_rate(self, relaxation_rate): e = sp.Rational(1, 1) ... ... @@ -200,18 +224,6 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod): def set_force_model(self, force_model): self._force_model = force_model @property def cumulants(self): return sorted(self._cumulant_to_relaxation_info_dict.keys(), key=moment_sort_key) @property def cumulant_equilibrium_values(self): return tuple([e.equilibrium_value for e in self._cumulant_to_relaxation_info_dict.values()]) @property def relaxation_rates(self): return tuple([e.relaxation_rate for e in self._cumulant_to_relaxation_info_dict.values()]) def _repr_html_(self): table = """ ... ...
This diff is collapsed.
 ... ... @@ -142,7 +142,7 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform): central_moments = self.forward_matrix * f_vec main_assignments = [Assignment(sq_sym(moment_symbol_base, e), eq) for e, eq in zip(self.moment_exponents, central_moments)] symbol_gen = SymbolGen(subexpression_base, dtype=float) symbol_gen = SymbolGen(subexpression_base) ac = AssignmentCollection(main_assignments, subexpression_symbol_generator=symbol_gen) if simplification: ... ... @@ -157,7 +157,7 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform): moment_vec = sp.Matrix(moments) pdfs_from_moments = self.backward_matrix * moment_vec main_assignments = [Assignment(f, eq) for f, eq in zip(pdf_symbols, pdfs_from_moments)] symbol_gen = SymbolGen(subexpression_base, dtype=float) symbol_gen = SymbolGen(subexpression_base) ac = AssignmentCollection(main_assignments, subexpression_symbol_generator=symbol_gen) if simplification: ... ... @@ -228,7 +228,7 @@ class FastCentralMomentTransform(AbstractMomentTransform): collect_partial_sums(e) subexpressions = [Assignment(lhs, rhs) for lhs, rhs in subexpressions_dict.items()] symbol_gen = SymbolGen(subexpression_base, dtype=float) symbol_gen = SymbolGen(subexpression_base) ac = AssignmentCollection(main_assignments, subexpressions=subexpressions, subexpression_symbol_generator=symbol_gen) if simplification: ... ... @@ -244,7 +244,7 @@ class FastCentralMomentTransform(AbstractMomentTransform): pdf_symbols, moment_symbol_base=POST_COLLISION_CENTRAL_MOMENT, simplification=False) raw_equations = raw_equations.new_without_subexpressions() symbol_gen = SymbolGen(subexpression_base, dtype=float)