creationfunctions.py 35.2 KB
Newer Older
Martin Bauer's avatar
Martin Bauer committed
1
2
3
4
import itertools
import operator
from collections import OrderedDict
from functools import reduce
Martin Bauer's avatar
Martin Bauer committed
5

Martin Bauer's avatar
Martin Bauer committed
6
import sympy as sp
Martin Bauer's avatar
Martin Bauer committed
7
8
9
10

from lbmpy.maxwellian_equilibrium import (
    compressible_to_incompressible_moment_value, get_cumulants_of_continuous_maxwellian_equilibrium,
    get_moments_of_continuous_maxwellian_equilibrium,
11
    get_moments_of_discrete_maxwellian_equilibrium, get_weights)
Frederik Hennig's avatar
Frederik Hennig committed
12

Markus Holzer's avatar
Markus Holzer committed
13
14
15
from lbmpy.methods.abstractlbmethod import RelaxationInfo

from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod
Frederik Hennig's avatar
Frederik Hennig committed
16
17
18
from lbmpy.methods.centeredcumulant.centered_cumulants import get_default_polynomial_cumulants_for_stencil
from lbmpy.methods.centeredcumulant.cumulant_transform import CentralMomentsToCumulantsByGeneratingFunc

Martin Bauer's avatar
Martin Bauer committed
19
from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation
Markus Holzer's avatar
Markus Holzer committed
20

Frederik Hennig's avatar
Frederik Hennig committed
21
from lbmpy.methods.momentbased.momentbasedmethod import MomentBasedLbMethod
Markus Holzer's avatar
Markus Holzer committed
22
from lbmpy.methods.momentbased.moment_transforms import PdfsToCentralMomentsByShiftMatrix
Frederik Hennig's avatar
Frederik Hennig committed
23

Martin Bauer's avatar
Martin Bauer committed
24
25
from lbmpy.moments import (
    MOMENT_SYMBOLS, discrete_moment, exponents_to_polynomial_representations,
Markus Holzer's avatar
Markus Holzer committed
26
    get_default_moment_set_for_stencil, gram_schmidt, is_even, moments_of_order,
Frederik Hennig's avatar
Frederik Hennig committed
27
28
29
    moments_up_to_component_order, sort_moments_into_groups_of_same_order,
    is_bulk_moment, is_shear_moment, get_order)

Martin Bauer's avatar
Martin Bauer committed
30
from lbmpy.relaxationrates import relaxation_rate_from_magic_number
31
from lbmpy.stencils import get_stencil
32
from pystencils.stencil import have_same_entries
Martin Bauer's avatar
Martin Bauer committed
33
from pystencils.sympyextensions import common_denominator
Martin Bauer's avatar
Martin Bauer committed
34
35


Martin Bauer's avatar
Martin Bauer committed
36
def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, compressible=False,
37
                                               force_model=None, equilibrium_order=2, c_s_sq=sp.Rational(1, 3)):
Martin Bauer's avatar
Martin Bauer committed
38
39
40
    r"""Creates a moment-based LBM by taking a list of moments with corresponding relaxation rate.

    These moments are relaxed against the moments of the discrete Maxwellian distribution.
Martin Bauer's avatar
Martin Bauer committed
41

Martin Bauer's avatar
Martin Bauer committed
42
    Args:
Martin Bauer's avatar
Martin Bauer committed
43
        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
Martin Bauer's avatar
Martin Bauer committed
44
        moment_to_relaxation_rate_dict: dict that has as many entries as the stencil. Each moment, which can be
Martin Bauer's avatar
Martin Bauer committed
45
46
                                        represented by an exponent tuple or in polynomial form
                                        (see `lbmpy.moments`), is mapped to a relaxation rate.
Martin Bauer's avatar
Martin Bauer committed
47
        compressible: incompressible LBM methods split the density into :math:`\rho = \rho_0 + \Delta \rho`
Martin Bauer's avatar
Martin Bauer committed
48
49
50
                      where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
                      This approximates the incompressible Navier-Stokes equations better than the standard
                      compressible model.
Martin Bauer's avatar
Martin Bauer committed
51
52
53
54
55
        force_model: force model instance, or None if no external forces
        equilibrium_order: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
        c_s_sq: Speed of sound squared

    Returns:
Frederik Hennig's avatar
Frederik Hennig committed
56
        `lbmpy.methods.momentbased.MomentBasedLbMethod` instance
Martin Bauer's avatar
Martin Bauer committed
57
    """
58
    if isinstance(stencil, str):
Martin Bauer's avatar
Martin Bauer committed
59
60
61
        stencil = get_stencil(stencil)
    mom_to_rr_dict = OrderedDict(moment_to_relaxation_rate_dict)
    assert len(mom_to_rr_dict) == len(stencil), \
Martin Bauer's avatar
Martin Bauer committed
62
63
        "The number of moments has to be the same as the number of stencil entries"

Martin Bauer's avatar
Martin Bauer committed
64
    density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
65
66
67
    eq_values = get_moments_of_discrete_maxwellian_equilibrium(stencil, tuple(mom_to_rr_dict.keys()),
                                                               c_s_sq=c_s_sq, compressible=compressible,
                                                               order=equilibrium_order)
Martin Bauer's avatar
Martin Bauer committed
68

Martin Bauer's avatar
Martin Bauer committed
69
70
    rr_dict = OrderedDict([(mom, RelaxationInfo(eq_mom, rr))
                           for mom, rr, eq_mom in zip(mom_to_rr_dict.keys(), mom_to_rr_dict.values(), eq_values)])
Frederik Hennig's avatar
Frederik Hennig committed
71
72

    return MomentBasedLbMethod(stencil, rr_dict, density_velocity_computation, force_model)
Martin Bauer's avatar
Martin Bauer committed
73
74


Martin Bauer's avatar
Martin Bauer committed
75
def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, compressible=False,
76
                                                 force_model=None, equilibrium_order=2, c_s_sq=sp.Rational(1, 3)):
Martin Bauer's avatar
Martin Bauer committed
77
78
79
    r"""
    Creates a moment-based LBM by taking a list of moments with corresponding relaxation rate. These moments are
    relaxed against the moments of the continuous Maxwellian distribution.
Martin Bauer's avatar
Martin Bauer committed
80
    For parameter description see :func:`lbmpy.methods.create_with_discrete_maxwellian_eq_moments`.
Martin Bauer's avatar
Martin Bauer committed
81
    By using the continuous Maxwellian we automatically get a compressible model.
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

    Args:
        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
        moment_to_relaxation_rate_dict: dict that has as many entries as the stencil. Each moment, which can be
                                        represented by an exponent tuple or in polynomial form
                                        (see `lbmpy.moments`), is mapped to a relaxation rate.
        compressible: incompressible LBM methods split the density into :math:`\rho = \rho_0 + \Delta \rho`
                      where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
                      This approximates the incompressible Navier-Stokes equations better than the standard
                      compressible model.
        force_model: force model instance, or None if no external forces
        equilibrium_order: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
        c_s_sq: Speed of sound squared

    Returns:
        `lbmpy.methods.momentbased.MomentBasedLbMethod` instance
Martin Bauer's avatar
Martin Bauer committed
98
    """
99
    if isinstance(stencil, str):
Martin Bauer's avatar
Martin Bauer committed
100
101
102
        stencil = get_stencil(stencil)
    mom_to_rr_dict = OrderedDict(moment_to_relaxation_rate_dict)
    assert len(mom_to_rr_dict) == len(stencil), "The number of moments has to be equal to the number of stencil entries"
Martin Bauer's avatar
Martin Bauer committed
103
    dim = len(stencil[0])
Martin Bauer's avatar
Martin Bauer committed
104
    density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
105
106
    eq_values = get_moments_of_continuous_maxwellian_equilibrium(tuple(mom_to_rr_dict.keys()), dim, c_s_sq=c_s_sq,
                                                                 order=equilibrium_order)
Martin Bauer's avatar
Martin Bauer committed
107
108

    if not compressible:
Martin Bauer's avatar
Martin Bauer committed
109
110
111
        rho = density_velocity_computation.defined_symbols(order=0)[1]
        u = density_velocity_computation.defined_symbols(order=1)[1]
        eq_values = [compressible_to_incompressible_moment_value(em, rho, u) for em in eq_values]
Martin Bauer's avatar
Martin Bauer committed
112

Martin Bauer's avatar
Martin Bauer committed
113
    rr_dict = OrderedDict([(mom, RelaxationInfo(eq_mom, rr))
Frederik Hennig's avatar
Frederik Hennig committed
114
                           for mom, rr, eq_mom in zip(mom_to_rr_dict.keys(), mom_to_rr_dict.values(), eq_values)])
115
116

    return MomentBasedLbMethod(stencil, rr_dict, density_velocity_computation, force_model)
Martin Bauer's avatar
Martin Bauer committed
117
118


Martin Bauer's avatar
Martin Bauer committed
119
def create_generic_mrt(stencil, moment_eq_value_relaxation_rate_tuples, compressible=False,
120
                       force_model=None):
121
    r"""
Martin Bauer's avatar
Martin Bauer committed
122
123
124
125
126
127
128
    Creates a generic moment-based LB method.

    Args:
        stencil: sequence of lattice velocities
        moment_eq_value_relaxation_rate_tuples: sequence of tuples containing (moment, equilibrium value, relax. rate)
        compressible: compressibility, determines calculation of velocity for force models
        force_model: see create_with_discrete_maxwellian_eq_moments
129
    """
Martin Bauer's avatar
Martin Bauer committed
130
    density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
131

Martin Bauer's avatar
Martin Bauer committed
132
    rr_dict = OrderedDict()
Martin Bauer's avatar
Martin Bauer committed
133
    for moment, eq_value, rr in moment_eq_value_relaxation_rate_tuples:
134
        moment = sp.sympify(moment)
Martin Bauer's avatar
Martin Bauer committed
135
        rr_dict[moment] = RelaxationInfo(eq_value, rr)
136
    return MomentBasedLbMethod(stencil, rr_dict, density_velocity_computation, force_model)
137
138


139
140
def create_from_equilibrium(stencil, equilibrium, moment_to_relaxation_rate_dict,
                            compressible=False, force_model=None):
141
142
    r"""
    Creates a moment-based LB method using a given equilibrium distribution function
Martin Bauer's avatar
Martin Bauer committed
143
144
145
146
147
148
149
150

    Args:
        stencil: see create_with_discrete_maxwellian_eq_moments
        equilibrium: list of equilibrium terms, dependent on rho and u, one for each stencil direction
        moment_to_relaxation_rate_dict: relaxation rate for each moment, or a symbol/float if all should relaxed with
                                        the same rate
        compressible: see create_with_discrete_maxwellian_eq_moments
        force_model: see create_with_discrete_maxwellian_eq_moments
151
    """
152
    if isinstance(stencil, str):
Martin Bauer's avatar
Martin Bauer committed
153
154
155
156
        stencil = get_stencil(stencil)
    if any(isinstance(moment_to_relaxation_rate_dict, t) for t in (sp.Symbol, float, int)):
        moment_to_relaxation_rate_dict = {m: moment_to_relaxation_rate_dict
                                          for m in get_default_moment_set_for_stencil(stencil)}
157

Martin Bauer's avatar
Martin Bauer committed
158
159
160
    mom_to_rr_dict = OrderedDict(moment_to_relaxation_rate_dict)
    assert len(mom_to_rr_dict) == len(stencil), "The number of moments has to be equal to the number of stencil entries"
    density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
161

Martin Bauer's avatar
Martin Bauer committed
162
    rr_dict = OrderedDict([(mom, RelaxationInfo(discrete_moment(equilibrium, mom, stencil).expand(), rr))
Frederik Hennig's avatar
Frederik Hennig committed
163
                           for mom, rr in zip(mom_to_rr_dict.keys(), mom_to_rr_dict.values())])
Martin Bauer's avatar
Martin Bauer committed
164
    return MomentBasedLbMethod(stencil, rr_dict, density_velocity_computation, force_model)
165
166


Martin Bauer's avatar
Martin Bauer committed
167
168
169
# ------------------------------------ SRT / TRT/ MRT Creators ---------------------------------------------------------


Martin Bauer's avatar
Martin Bauer committed
170
171
def create_srt(stencil, relaxation_rate, maxwellian_moments=False, **kwargs):
    r"""Creates a single relaxation time (SRT) lattice Boltzmann model also known as BGK model.
Martin Bauer's avatar
Martin Bauer committed
172

Martin Bauer's avatar
Martin Bauer committed
173
174
175
176
177
178
179
180
    Args:
        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
        relaxation_rate: relaxation rate (inverse of the relaxation time)
                        usually called :math:`\omega` in LBM literature
        maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
                        used to compute the equilibrium moments

    Returns:
Frederik Hennig's avatar
Frederik Hennig committed
181
        :class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
Martin Bauer's avatar
Martin Bauer committed
182
    """
183
    if isinstance(stencil, str):
Martin Bauer's avatar
Martin Bauer committed
184
185
186
187
188
        stencil = get_stencil(stencil)
    moments = get_default_moment_set_for_stencil(stencil)
    rr_dict = OrderedDict([(m, relaxation_rate) for m in moments])
    if maxwellian_moments:
        return create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, **kwargs)
Martin Bauer's avatar
Martin Bauer committed
189
    else:
Martin Bauer's avatar
Martin Bauer committed
190
        return create_with_discrete_maxwellian_eq_moments(stencil, rr_dict, **kwargs)
Martin Bauer's avatar
Martin Bauer committed
191
192


Martin Bauer's avatar
Martin Bauer committed
193
194
def create_trt(stencil, relaxation_rate_even_moments, relaxation_rate_odd_moments,
               maxwellian_moments=False, **kwargs):
Martin Bauer's avatar
Martin Bauer committed
195
196
197
198
199
    """
    Creates a two relaxation time (TRT) lattice Boltzmann model, where even and odd moments are relaxed differently.
    In the SRT model the exact wall position of no-slip boundaries depends on the viscosity, the TRT method does not
    have this problem.

Martin Bauer's avatar
Martin Bauer committed
200
    Parameters are similar to :func:`lbmpy.methods.create_srt`, but instead of one relaxation rate there are
Martin Bauer's avatar
Martin Bauer committed
201
    two relaxation rates: one for even moments (determines viscosity) and one for odd moments.
Martin Bauer's avatar
Martin Bauer committed
202
    If unsure how to choose the odd relaxation rate, use the function :func:`lbmpy.methods.create_trt_with_magic_number`
Martin Bauer's avatar
Martin Bauer committed
203
    """
204
    if isinstance(stencil, str):
Martin Bauer's avatar
Martin Bauer committed
205
206
207
208
209
210
        stencil = get_stencil(stencil)
    moments = get_default_moment_set_for_stencil(stencil)
    rr_dict = OrderedDict([(m, relaxation_rate_even_moments if is_even(m) else relaxation_rate_odd_moments)
                           for m in moments])
    if maxwellian_moments:
        return create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, **kwargs)
Martin Bauer's avatar
Martin Bauer committed
211
    else:
Martin Bauer's avatar
Martin Bauer committed
212
        return create_with_discrete_maxwellian_eq_moments(stencil, rr_dict, **kwargs)
Martin Bauer's avatar
Martin Bauer committed
213
214


Martin Bauer's avatar
Martin Bauer committed
215
def create_trt_with_magic_number(stencil, relaxation_rate, magic_number=sp.Rational(3, 16), **kwargs):
216
    r"""
Martin Bauer's avatar
Martin Bauer committed
217
218
    Creates a two relaxation time (TRT) lattice Boltzmann method, where the relaxation time for odd moments is
    determines from the even moment relaxation time and a "magic number".
Martin Bauer's avatar
Martin Bauer committed
219
    For possible parameters see :func:`lbmpy.methods.create_trt`
Markus Holzer's avatar
Markus Holzer committed
220

221
222
223
224
225
226
227
228
    Args:
        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
        relaxation_rate: relaxation rate (inverse of the relaxation time)
                        usually called :math:`\omega` in LBM literature
        magic_number: magic number which is used to calculate the relxation rate for the odd moments.

    Returns:
        :class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
Martin Bauer's avatar
Martin Bauer committed
229
    """
Martin Bauer's avatar
Martin Bauer committed
230
231
232
    rr_odd = relaxation_rate_from_magic_number(relaxation_rate, magic_number)
    return create_trt(stencil, relaxation_rate_even_moments=relaxation_rate,
                      relaxation_rate_odd_moments=rr_odd, **kwargs)
Martin Bauer's avatar
Martin Bauer committed
233
234


Martin Bauer's avatar
Martin Bauer committed
235
def create_mrt_raw(stencil, relaxation_rates, maxwellian_moments=False, **kwargs):
236
237
    r"""
    Creates a MRT method using non-orthogonalized moments.
Markus Holzer's avatar
Markus Holzer committed
238

239
240
241
242
243
244
245
246
    Args:
        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
        relaxation_rates: relaxation rates (inverse of the relaxation times) for each moment
        maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
                        used to compute the equilibrium moments.
    Returns:
        :class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
    """
247
    if isinstance(stencil, str):
Martin Bauer's avatar
Martin Bauer committed
248
249
250
251
252
        stencil = get_stencil(stencil)
    moments = get_default_moment_set_for_stencil(stencil)
    rr_dict = OrderedDict(zip(moments, relaxation_rates))
    if maxwellian_moments:
        return create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, **kwargs)
253
    else:
Martin Bauer's avatar
Martin Bauer committed
254
        return create_with_discrete_maxwellian_eq_moments(stencil, rr_dict, **kwargs)
255
256


Martin Bauer's avatar
Martin Bauer committed
257
258
def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, method_name='KBC-N4',
                   maxwellian_moments=False, **kwargs):
259
260
261
    """
    Creates a method with two relaxation rates, one for lower order moments which determines the viscosity and
    one for higher order moments. In entropic models this second relaxation rate is chosen subject to an entropy
Martin Bauer's avatar
Martin Bauer committed
262
    condition. Which moments are relaxed by which rate is determined by the method_name
263

264
265
266
267
268
269
270
    Args:
        dim: 2 or 3, leads to stencil D2Q9 or D3Q27
        shear_relaxation_rate: relaxation rate that determines viscosity
        higher_order_relaxation_rate: relaxation rate for higher order moments
        method_name: string 'KBC-Nx' where x can be an number from 1 to 4, for details see
                    "Karlin 2015: Entropic multi relaxation lattice Boltzmann models for turbulent flows"
        maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
271
                        used to compute the equilibrium moments.
272
273
274
275
    """
    def product(iterable):
        return reduce(operator.mul, iterable, 1)

Martin Bauer's avatar
Martin Bauer committed
276
    the_moment = MOMENT_SYMBOLS[:dim]
277
278

    rho = [sp.Rational(1, 1)]
Martin Bauer's avatar
Martin Bauer committed
279
    velocity = list(the_moment)
280

Martin Bauer's avatar
Martin Bauer committed
281
282
283
284
    shear_tensor_off_diagonal = [product(t) for t in itertools.combinations(the_moment, 2)]
    shear_tensor_diagonal = [m_i * m_i for m_i in the_moment]
    shear_tensor_trace = sum(shear_tensor_diagonal)
    shear_tensor_trace_free_diagonal = [dim * d - shear_tensor_trace for d in shear_tensor_diagonal]
285

Martin Bauer's avatar
Martin Bauer committed
286
287
288
    poly_repr = exponents_to_polynomial_representations
    energy_transport_tensor = list(poly_repr([a for a in moments_of_order(3, dim, True)
                                              if 3 not in a]))
289

Martin Bauer's avatar
Martin Bauer committed
290
291
    explicitly_defined = set(rho + velocity + shear_tensor_off_diagonal
                             + shear_tensor_diagonal + energy_transport_tensor)
Martin Bauer's avatar
Martin Bauer committed
292
    rest = list(set(poly_repr(moments_up_to_component_order(2, dim))) - explicitly_defined)
Martin Bauer's avatar
Martin Bauer committed
293
    assert len(rest) + len(explicitly_defined) == 3**dim
294

Martin Bauer's avatar
Martin Bauer committed
295
    # naming according to paper Karlin 2015: Entropic multi relaxation lattice Boltzmann models for turbulent flows
Martin Bauer's avatar
Martin Bauer committed
296
297
298
299
    d = shear_tensor_off_diagonal + shear_tensor_trace_free_diagonal[:-1]
    t = [shear_tensor_trace]
    q = energy_transport_tensor
    if method_name == 'KBC-N1':
Martin Bauer's avatar
Martin Bauer committed
300
        decomposition = [d, t + q + rest]
Martin Bauer's avatar
Martin Bauer committed
301
    elif method_name == 'KBC-N2':
Martin Bauer's avatar
Martin Bauer committed
302
        decomposition = [d + t, q + rest]
Martin Bauer's avatar
Martin Bauer committed
303
    elif method_name == 'KBC-N3':
Martin Bauer's avatar
Martin Bauer committed
304
        decomposition = [d + q, t + rest]
Martin Bauer's avatar
Martin Bauer committed
305
    elif method_name == 'KBC-N4':
Martin Bauer's avatar
Martin Bauer committed
306
        decomposition = [d + t + q, rest]
307
308
309
    else:
        raise ValueError("Unknown model. Supported models KBC-Nx where x in (1,2,3,4)")

Martin Bauer's avatar
Martin Bauer committed
310
311
    omega_s, omega_h = shear_relaxation_rate, higher_order_relaxation_rate
    shear_part, rest_part = decomposition
312

Martin Bauer's avatar
Martin Bauer committed
313
314
315
316
    relaxation_rates = [omega_s] + \
                       [omega_s] * len(velocity) + \
                       [omega_s] * len(shear_part) + \
                       [omega_h] * len(rest_part)
317

Martin Bauer's avatar
Martin Bauer committed
318
319
320
    stencil = get_stencil("D2Q9") if dim == 2 else get_stencil("D3Q27")
    all_moments = rho + velocity + shear_part + rest_part
    moment_to_rr = OrderedDict((m, rr) for m, rr in zip(all_moments, relaxation_rates))
321

Martin Bauer's avatar
Martin Bauer committed
322
323
    if maxwellian_moments:
        return create_with_continuous_maxwellian_eq_moments(stencil, moment_to_rr, **kwargs)
324
    else:
Martin Bauer's avatar
Martin Bauer committed
325
        return create_with_discrete_maxwellian_eq_moments(stencil, moment_to_rr, **kwargs)
326
327


Martin Bauer's avatar
Martin Bauer committed
328
def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=False, weighted=None,
329
                          nested_moments=None, **kwargs):
Martin Bauer's avatar
Martin Bauer committed
330
    r"""
331
    Returns an orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
Martin Bauer's avatar
Martin Bauer committed
332
333
    These MRT methods are just one specific version - there are many MRT methods possible for all these stencils
    which differ by the linear combination of moments and the grouping into equal relaxation times.
Martin Bauer's avatar
Martin Bauer committed
334
    To create a generic MRT method use `create_with_discrete_maxwellian_eq_moments`
Martin Bauer's avatar
Martin Bauer committed
335

Martin Bauer's avatar
Martin Bauer committed
336
337
338
339
    Args:
        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
        relaxation_rate_getter: function getting a list of moments as argument, returning the associated relaxation
        maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
340
                                               used to compute the equilibrium moments
341
        weighted: whether to use weighted or unweighted orthogonality
Martin Bauer's avatar
Martin Bauer committed
342
343
344
        nested_moments: a list of lists of modes, grouped by common relaxation times. This is usually used in
                        conjunction with `mrt_orthogonal_modes_literature`. If this argument is not provided,
                        Gram-Schmidt orthogonalization of the default modes is performed.
Martin Bauer's avatar
Martin Bauer committed
345
    """
346
    dim = len(stencil[0])
347
    if isinstance(stencil, str):
Martin Bauer's avatar
Martin Bauer committed
348
        stencil = get_stencil(stencil)
Martin Bauer's avatar
Martin Bauer committed
349

350
    if weighted is None and not nested_moments:
Martin Bauer's avatar
Martin Bauer committed
351
        raise ValueError("Please specify whether you want weighted or unweighted orthogonality using 'weighted='")
352
353
354
355
    elif weighted:
        weights = get_weights(stencil, sp.Rational(1, 3))
    else:
        weights = None
Martin Bauer's avatar
Martin Bauer committed
356

Martin Bauer's avatar
Martin Bauer committed
357
    moment_to_relaxation_rate_dict = OrderedDict()
358
    if not nested_moments:
Martin Bauer's avatar
Martin Bauer committed
359
        moments = get_default_moment_set_for_stencil(stencil)
360
361
362
363
364
365
366
        x, y, z = MOMENT_SYMBOLS
        if dim == 2:
            diagonal_viscous_moments = [x ** 2 + y ** 2, x ** 2]
        else:
            diagonal_viscous_moments = [x ** 2 + y ** 2 + z ** 2, x ** 2, y ** 2 - z ** 2]
        for i, d in enumerate(MOMENT_SYMBOLS[:dim]):
            moments[moments.index(d**2)] = diagonal_viscous_moments[i]
367
        orthogonal_moments = gram_schmidt(moments, stencil, weights)
Martin Bauer's avatar
Martin Bauer committed
368
369
        orthogonal_moments_scaled = [e * common_denominator(e) for e in orthogonal_moments]
        nested_moments = list(sort_moments_into_groups_of_same_order(orthogonal_moments_scaled).values())
Martin Bauer's avatar
Martin Bauer committed
370
371
372
373
374
375
376
        # second order moments: separate bulk from shear
        second_order_moments = nested_moments[2]
        bulk_moment = [m for m in second_order_moments if is_bulk_moment(m, dim)]
        shear_moments = [m for m in second_order_moments if is_shear_moment(m, dim)]
        assert len(shear_moments) + len(bulk_moment) == len(second_order_moments)
        nested_moments[2] = shear_moments
        nested_moments.insert(3, bulk_moment)
377
378
379
380
381
382
383
384
385
386
387
    for moment_list in nested_moments:
        rr = relaxation_rate_getter(moment_list)
        for m in moment_list:
            moment_to_relaxation_rate_dict[m] = rr

    if maxwellian_moments:
        return create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)
    else:
        return create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)


388
def mrt_orthogonal_modes_literature(stencil, is_weighted):
389
390
391
392
393
394
395
    """
    Returns a list of lists of modes, grouped by common relaxation times.
    This is for commonly used MRT models found in literature.

    Args:
        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
        is_weighted: whether to use weighted or unweighted orthogonality
Martin Bauer's avatar
Martin Bauer committed
396
397

    MRT schemes as described in the following references are used
398
399
400
401
    """
    x, y, z = MOMENT_SYMBOLS
    one = sp.Rational(1, 1)

402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
    if have_same_entries(stencil, get_stencil("D2Q9")) and is_weighted:
        # Reference:
        # Duenweg, B., Schiller, U. D., & Ladd, A. J. (2007). Statistical mechanics of the fluctuating
        # lattice Boltzmann equation. Physical Review E, 76(3)
        sq = x ** 2 + y ** 2
        all_moments = [one, x, y, 3 * sq - 2, 2 * x ** 2 - sq, x * y,
                       (3 * sq - 4) * x, (3 * sq - 4) * y, 9 * sq ** 2 - 15 * sq + 2]
        nested_moments = list(sort_moments_into_groups_of_same_order(all_moments).values())
        return nested_moments
    elif have_same_entries(stencil, get_stencil("D3Q15")) and is_weighted:
        sq = x ** 2 + y ** 2 + z ** 2
        nested_moments = [
            [one, x, y, z],  # [0, 3, 5, 7]
            [sq - 1],  # [1]
            [3 * sq ** 2 - 9 * sq + 4],  # [2]
            [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z],  # [4, 6, 8]
            [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z],  # [9, 10, 11, 12, 13]
            [x * y * z]
        ]
    elif have_same_entries(stencil, get_stencil("D3Q19")) and is_weighted:
        # This MRT variant mentioned in the dissertation of Ulf Schiller
        # "Thermal fluctuations and boundary conditions in the lattice Boltzmann method" (2008), p. 24ff
        # There are some typos in the moment matrix on p.27
        # The here implemented ordering of the moments is however different from that reference (Eq. 2.61-2.63)
        # The moments are weighted-orthogonal (Eq. 2.58)

        # Further references:
        # Duenweg, B., Schiller, U. D., & Ladd, A. J. (2007). Statistical mechanics of the fluctuating
        # lattice Boltzmann equation. Physical Review E, 76(3)
        # Chun, B., & Ladd, A. J. (2007). Interpolated boundary condition for lattice Boltzmann simulations of
        # flows in narrow gaps. Physical review E, 75(6)
        sq = x ** 2 + y ** 2 + z ** 2
        nested_moments = [
            [one, x, y, z],  # [0, 3, 5, 7]
            [sq - 1],  # [1]
            [3 * sq ** 2 - 6 * sq + 1],  # [2]
            [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z],  # [4, 6, 8]
            [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z],  # [9, 11, 13, 14, 15]
            [(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)],  # [10, 12]
            [(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z]  # [16, 17, 18]
        ]
    elif have_same_entries(stencil, get_stencil("D3Q27")) and not is_weighted:
        xsq, ysq, zsq = x ** 2, y ** 2, z ** 2
        all_moments = [
            sp.Rational(1, 1),  # 0
            x, y, z,  # 1, 2, 3
            x * y, x * z, y * z,  # 4, 5, 6
            xsq - ysq,  # 7
            (xsq + ysq + zsq) - 3 * zsq,  # 8
            (xsq + ysq + zsq) - 2,  # 9
            3 * (x * ysq + x * zsq) - 4 * x,  # 10
            3 * (xsq * y + y * zsq) - 4 * y,  # 11
            3 * (xsq * z + ysq * z) - 4 * z,  # 12
            x * ysq - x * zsq,  # 13
            xsq * y - y * zsq,  # 14
            xsq * z - ysq * z,  # 15
            x * y * z,  # 16
            3 * (xsq * ysq + xsq * zsq + ysq * zsq) - 4 * (xsq + ysq + zsq) + 4,  # 17
            3 * (xsq * ysq + xsq * zsq - 2 * ysq * zsq) - 2 * (2 * xsq - ysq - zsq),  # 18
            3 * (xsq * ysq - xsq * zsq) - 2 * (ysq - zsq),  # 19
            3 * (xsq * y * z) - 2 * (y * z),  # 20
            3 * (x * ysq * z) - 2 * (x * z),  # 21
            3 * (x * y * zsq) - 2 * (x * y),  # 22
            9 * (x * ysq * zsq) - 6 * (x * ysq + x * zsq) + 4 * x,  # 23
            9 * (xsq * y * zsq) - 6 * (xsq * y + y * zsq) + 4 * y,  # 24
            9 * (xsq * ysq * z) - 6 * (xsq * z + ysq * z) + 4 * z,  # 25
            27 * (xsq * ysq * zsq) - 18 * (xsq * ysq + xsq * zsq + ysq * zsq) + 12 * (xsq + ysq + zsq) - 8,  # 26
        ]
        nested_moments = list(sort_moments_into_groups_of_same_order(all_moments).values())
Martin Bauer's avatar
Martin Bauer committed
471
    else:
472
473
        raise NotImplementedError("No MRT model is available (yet) for this stencil. "
                                  "Create a custom MRT using 'create_with_discrete_maxwellian_eq_moments'")
Martin Bauer's avatar
Martin Bauer committed
474

475
    return nested_moments
Martin Bauer's avatar
Martin Bauer committed
476

Frederik Hennig's avatar
Frederik Hennig committed
477
478
479
480
481
482
483
484
485
486
487
488
489
# ----------------------------------------- Cumulant method creators ---------------------------------------------------


def create_centered_cumulant_model(stencil, cumulant_to_rr_dict, force_model=None,
                                   equilibrium_order=None, c_s_sq=sp.Rational(1, 3),
                                   galilean_correction=False,
                                   central_moment_transform_class=PdfsToCentralMomentsByShiftMatrix,
                                   cumulant_transform_class=CentralMomentsToCumulantsByGeneratingFunc):
    r"""Creates a cumulant lattice Boltzmann model.

    Args:
        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
        cumulant_to_rr_dict: dict that has as many entries as the stencil. Each cumulant, which can be
Markus Holzer's avatar
Markus Holzer committed
490
491
                             represented by an exponent tuple or in polynomial form is mapped to a relaxation rate.
                             See `get_default_polynomial_cumulants_for_stencil`
Frederik Hennig's avatar
Frederik Hennig committed
492
        force_model: force model used for the collision. For cumulant LB method a good choice is
Markus Holzer's avatar
Markus Holzer committed
493
                     `lbmpy.methods.centeredcumulant.CenteredCumulantForceModel`
Frederik Hennig's avatar
Frederik Hennig committed
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
        equilibrium_order: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
        c_s_sq: Speed of sound squared
        galilean_correction: special correction for D3Q27 cumulant collisions. See Appendix H in
                             :cite:`geier2015`. Implemented in :mod:`lbmpy.methods.centeredcumulant.galilean_correction`
        central_moment_transform_class: Class which defines the transformation to the central moment space
                                        (see :mod:`lbmpy.methods.momentbased.moment_transforms`)
        cumulant_transform_class: Class which defines the transformation from the central moment space to the
                                  cumulant space (see :mod:`lbmpy.methods.centeredcumulant.cumulant_transform`)

    Returns:
        :class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` instance
    """

    one = sp.Integer(1)
    if isinstance(stencil, str):
        stencil = get_stencil(stencil)

    assert len(cumulant_to_rr_dict) == len(
        stencil), "The number of moments has to be equal to the number of stencil entries"
    dim = len(stencil[0])

    # CQC
    cqc = DensityVelocityComputation(stencil, True, force_model=force_model)
    density_symbol = cqc.zeroth_order_moment_symbol
    velocity_symbols = cqc.first_order_moment_symbols

    #   Equilibrium Values
    higher_order_polynomials = list(filter(lambda x: get_order(x) > 1, cumulant_to_rr_dict.keys()))

    #   Lower Order Equilibria
    cumulants_to_relaxation_info_dict = {one: RelaxationInfo(density_symbol, cumulant_to_rr_dict[one])}
    for d in MOMENT_SYMBOLS[:dim]:
        cumulants_to_relaxation_info_dict[d] = RelaxationInfo(0, cumulant_to_rr_dict[d])

    #   Polynomial Cumulant Equilibria
    polynomial_equilibria = get_cumulants_of_continuous_maxwellian_equilibrium(
        higher_order_polynomials, dim, rho=density_symbol, u=velocity_symbols, c_s_sq=c_s_sq, order=equilibrium_order)
    polynomial_equilibria = [density_symbol * v for v in polynomial_equilibria]

    for i, c in enumerate(higher_order_polynomials):
        cumulants_to_relaxation_info_dict[c] = RelaxationInfo(polynomial_equilibria[i], cumulant_to_rr_dict[c])

    return CenteredCumulantBasedLbMethod(stencil, cumulants_to_relaxation_info_dict, cqc, force_model,
                                         galilean_correction=galilean_correction,
                                         central_moment_transform_class=central_moment_transform_class,
                                         cumulant_transform_class=cumulant_transform_class)


def create_with_polynomial_cumulants(stencil, relaxation_rates, cumulant_groups, **kwargs):
    r"""Creates a cumulant lattice Boltzmann model based on a default polynomial set.

    Args:
        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
        relaxation_rates: relaxation rates for each cumulant group. If None are provided a list of symbolic relaxation
                          rates is created and used. If only a list with one entry is provided this relaxation rate is
                          used for determine the viscosity of the simulation. All other cumulants are relaxed with one.
                          If a cumulant force model is provided the first order cumulants are relaxed with two to ensure
                          that the force is applied correctly to the momentum groups
        cumulant_groups: Nested sequence of polynomial expressions defining the cumulants to be relaxed. All cumulants 
                         within one group are relaxed with the same relaxation rate.
        kwargs: See :func:`create_centered_cumulant_model`

    Returns:
        :class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` instance
    """
    if isinstance(stencil, str):
        stencil = get_stencil(stencil)

    cumulant_to_rr_dict = OrderedDict()
    rr_iter = iter(relaxation_rates)
    for group in cumulant_groups:
        if all(get_order(c) <= 1 for c in group):
            for cumulant in group:
                cumulant_to_rr_dict[cumulant] = 0
        else:
            try:
                rr = next(rr_iter)
                for cumulant in group:
                    cumulant_to_rr_dict[cumulant] = rr
            except StopIteration:
                raise ValueError('Not enough relaxation rates specified.')

    return create_centered_cumulant_model(stencil, cumulant_to_rr_dict, **kwargs)


def create_with_monomial_cumulants(stencil, relaxation_rates, **kwargs):
    r"""Creates a cumulant lattice Boltzmann model based on a default polinomial set.

    Args:
        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
        relaxation_rates: relaxation rates for each cumulant group. If None are provided a list of symbolic relaxation
                          rates is created and used. If only a list with one entry is provided this relaxation rate is
                          used for determine the viscosity of the simulation. All other cumulants are relaxed with one.
                          If a cumulant force model is provided the first order cumulants are relaxed with two to ensure
                          that the force is applied correctly to the momentum groups
        kwargs: See :func:`create_centered_cumulant_model`

    Returns:
        :class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` instance
    """

    if isinstance(stencil, str):
        stencil = get_stencil(stencil)

    dim = len(stencil[0])

    # Get monomial moments
    cumulants = get_default_moment_set_for_stencil(stencil)

    if len(relaxation_rates) == 1:
        r_rates = []
        for c in cumulants:
            order = get_order(c)
            if order <= 1:
                #   Conserved moments
                continue
            elif is_shear_moment(c, dim):
                #   Viscosity-governing moments
                r_rates.append(relaxation_rates[0])
            else:
                #   The rest
                r_rates.append(1)
    else:
        r_rates = relaxation_rates

    cumulant_groups = [(c,) for c in cumulants]

    return create_with_polynomial_cumulants(stencil, r_rates, cumulant_groups, **kwargs)


def create_with_default_polynomial_cumulants(stencil, relaxation_rates, **kwargs):
    r"""Creates a cumulant lattice Boltzmann model based on a default polynomial set.

    Args: See :func:`create_with_polynomial_cumulants`.

    Returns:
        :class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` instance
    """

    if isinstance(stencil, str):
        stencil = get_stencil(stencil)

    # Get polynomial groups
    cumulant_groups = get_default_polynomial_cumulants_for_stencil(stencil)

    if len(relaxation_rates) == 1:
        r_rates = [relaxation_rates[0]]     # For correct viscosity
        r_rates += [1] * (len(cumulant_groups) - 3)
    else:
        assert len(relaxation_rates) >= len(cumulant_groups) - 2, \
            f"Number of relaxation rates must at least match the number of non-conserved cumulant groups. " \
            f"For this stencil we have {len(cumulant_groups) - 2} such cumulant groups"
        r_rates = relaxation_rates

    return create_with_polynomial_cumulants(stencil, r_rates, cumulant_groups, **kwargs)

Martin Bauer's avatar
Martin Bauer committed
650
651
652
653

# ----------------------------------------- Comparison view for notebooks ----------------------------------------------


Martin Bauer's avatar
Martin Bauer committed
654
def compare_moment_based_lb_methods(reference, other, show_deviations_only=False):
Martin Bauer's avatar
Martin Bauer committed
655
656
    import ipy_table
    table = []
Martin Bauer's avatar
Martin Bauer committed
657
    caption_rows = [len(table)]
Martin Bauer's avatar
Martin Bauer committed
658
659
    table.append(['Shared Moment', 'ref', 'other', 'difference'])

Martin Bauer's avatar
Martin Bauer committed
660
661
662
    reference_moments = set(reference.moments)
    other_moments = set(other.moments)
    for moment in reference_moments.intersection(other_moments):
Martin Bauer's avatar
Martin Bauer committed
663
664
        reference_value = reference.relaxation_info_dict[moment].equilibrium_value
        other_value = other.relaxation_info_dict[moment].equilibrium_value
Martin Bauer's avatar
Martin Bauer committed
665
666
        diff = sp.simplify(reference_value - other_value)
        if show_deviations_only and diff == 0:
Martin Bauer's avatar
Martin Bauer committed
667
668
669
            pass
        else:
            table.append(["$%s$" % (sp.latex(moment), ),
Martin Bauer's avatar
Martin Bauer committed
670
671
                          "$%s$" % (sp.latex(reference_value), ),
                          "$%s$" % (sp.latex(other_value), ),
Martin Bauer's avatar
Martin Bauer committed
672
                          "$%s$" % (sp.latex(diff),)])
Martin Bauer's avatar
Martin Bauer committed
673

Martin Bauer's avatar
Martin Bauer committed
674
675
676
    only_in_ref = reference_moments - other_moments
    if only_in_ref:
        caption_rows.append(len(table))
Martin Bauer's avatar
Martin Bauer committed
677
        table.append(['Only in Ref', 'value', '', ''])
Martin Bauer's avatar
Martin Bauer committed
678
        for moment in only_in_ref:
Martin Bauer's avatar
Martin Bauer committed
679
            val = reference.relaxation_info_dict[moment].equilibrium_value
Martin Bauer's avatar
Martin Bauer committed
680
681
682
683
            table.append(["$%s$" % (sp.latex(moment),),
                          "$%s$" % (sp.latex(val),),
                          " ", " "])

Martin Bauer's avatar
Martin Bauer committed
684
685
686
    only_in_other = other_moments - reference_moments
    if only_in_other:
        caption_rows.append(len(table))
Martin Bauer's avatar
Martin Bauer committed
687
        table.append(['Only in Other', '', 'value', ''])
Martin Bauer's avatar
Martin Bauer committed
688
        for moment in only_in_other:
Martin Bauer's avatar
Martin Bauer committed
689
            val = other.relaxation_info_dict[moment].equilibrium_value
Martin Bauer's avatar
Martin Bauer committed
690
691
692
693
694
            table.append(["$%s$" % (sp.latex(moment),),
                          " ",
                          "$%s$" % (sp.latex(val),),
                          " "])

Martin Bauer's avatar
Martin Bauer committed
695
696
    table_display = ipy_table.make_table(table)
    for row_idx in caption_rows:
Martin Bauer's avatar
Martin Bauer committed
697
        for col in range(4):
Martin Bauer's avatar
Martin Bauer committed
698
699
            ipy_table.set_cell_style(row_idx, col, color='#bbbbbb')
    return table_display