centered_cumulants.py 3.85 KB
Newer Older
Frederik Hennig's avatar
Frederik Hennig committed
1
2
3
4
5
from lbmpy.stencils import get_stencil
import sympy as sp

from pystencils.stencil import have_same_entries

Frederik Hennig's avatar
Frederik Hennig committed
6
from lbmpy.moments import MOMENT_SYMBOLS
Frederik Hennig's avatar
Frederik Hennig committed
7
8
9
10
11
12


def get_default_polynomial_cumulants_for_stencil(stencil):
    """
    Returns default groups of cumulants to be relaxed with common relaxation rates as stated in literature.
    Groups are ordered like this:
Markus Holzer's avatar
Markus Holzer committed
13
14
15
16
17
18
19
20
21

    - First group is density
    - Second group are the momentum modes
    - Third group are the shear modes
    - Fourth group is the bulk mode
    - Remaining groups do not govern hydrodynamic properties

    Args:
        stencil: can be D2Q9, D2Q19 or D3Q27
Frederik Hennig's avatar
Frederik Hennig committed
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
    """
    x, y, z = MOMENT_SYMBOLS
    if have_same_entries(stencil, get_stencil("D2Q9")):
        #   Cumulants of the D2Q9 stencil up to third order are equal to
        #   the central moments; only the fourth-order cumulant x**2 * y**2
        #   has a more complicated form. They can be arranged into groups
        #   for the preservation of rotational invariance as described by
        #   Martin Geier in his dissertation.
        #
        #   Reference: Martin Geier. Ab inito derivation of the cascaded Lattice Boltzmann
        #   Automaton. Dissertation. University of Freiburg. 2006.
        return [
            [sp.sympify(1)],        # density is conserved
            [x, y],                 # momentum is relaxed for cumulant forcing

            [x * y, x**2 - y**2],   # shear

            [x**2 + y**2],          # bulk

            [x**2 * y, x * y**2],
            [x**2 * y**2]
        ]

    elif have_same_entries(stencil, get_stencil("D3Q19")):
        #   D3Q19 cumulants are obtained by pruning the D3Q27 cumulant set as
        #   described by Coreixas, 2019.
        return [
            [sp.sympify(1)],                # density is conserved
            [x, y, z],                      # momentum might be affected by forcing

            [x * y,
             x * z,
             y * z,
             x ** 2 - y ** 2,
             x ** 2 - z ** 2],              # shear

            [x ** 2 + y ** 2 + z ** 2],     # bulk

            [x * y ** 2 + x * z ** 2,
             x ** 2 * y + y * z ** 2,
             x ** 2 * z + y ** 2 * z],

            [x * y ** 2 - x * z ** 2,
             x ** 2 * y - y * z ** 2,
             x ** 2 * z - y ** 2 * z],

            [x ** 2 * y ** 2 - 2 * x ** 2 * z ** 2 + y ** 2 * z ** 2,
             x ** 2 * y ** 2 + x ** 2 * z ** 2 - 2 * y ** 2 * z ** 2],

            [x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2]
        ]

    elif have_same_entries(stencil, get_stencil("D3Q27")):
        #   Cumulants grouped to preserve rotational invariance as described by Geier et al, 2015
        return [
            [sp.sympify(1)],                # density is conserved
            [x, y, z],                      # momentum might be affected by forcing

            [x * y,
             x * z,
             y * z,
             x ** 2 - y ** 2,
             x ** 2 - z ** 2],              # shear

            [x ** 2 + y ** 2 + z ** 2],     # bulk

            [x * y ** 2 + x * z ** 2,
             x ** 2 * y + y * z ** 2,
             x ** 2 * z + y ** 2 * z],

            [x * y ** 2 - x * z ** 2,
             x ** 2 * y - y * z ** 2,
             x ** 2 * z - y ** 2 * z],

            [x * y * z],

            [x ** 2 * y ** 2 - 2 * x ** 2 * z ** 2 + y ** 2 * z ** 2,
             x ** 2 * y ** 2 + x ** 2 * z ** 2 - 2 * y ** 2 * z ** 2],

            [x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2],

            [x ** 2 * y * z,
             x * y ** 2 * z,
             x * y * z ** 2],

            [x ** 2 * y ** 2 * z,
             x ** 2 * y * z ** 2,
             x * y ** 2 * z ** 2],

            [x ** 2 * y ** 2 * z ** 2]
        ]
    else:
        raise ValueError("No default set of cumulants is available for this stencil. "
                         "Please specify your own set of polynomial cumulants.")