forcemodels.py 9.37 KB
Newer Older
1
r"""
Martin Bauer's avatar
Martin Bauer committed
2
3
4
.. module:: forcemodels
    :synopsis: Collection of forcing terms for hydrodynamic LBM simulations

5
6

Get started:
7
------------
8

Martin Bauer's avatar
Martin Bauer committed
9
This module offers different models to introduce a body force in the lattice Boltzmann scheme.
10
If you don't know which model to choose, use :class:`lbmpy.forcemodels.Schiller`.
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
For incompressible collision models the :class:`lbmpy.forcemodels.Buick` model can be better.


Detailed information:
---------------------

Force models add a term :math:`C_F` to the collision equation:

.. math ::

    f(\pmb{x} + c_q \Delta t, t + \Delta t) - f(\pmb{x},t) = \Omega(f, f^{(eq)})
                                                            + \underbrace{F_q}_{\mbox{forcing term}}

The form of this term depends on the concrete force model: the first moment of this forcing term is equal
to the acceleration :math:`\pmb{a}` for all force models.

.. math ::

Martin Bauer's avatar
Martin Bauer committed
29
    \sum_q \pmb{c}_q F_q = \pmb{a}
30
31


Martin Bauer's avatar
Martin Bauer committed
32
The second order moment is different for the forcing models - if it is zero the model is suited for
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
incompressible flows. For weakly compressible collision operators a force model with a corrected second order moment
should be chosen.

.. math ::

    \sum_q c_{qi} c_{qj} f_q = F_i u_j + F_j u_i  \hspace{1cm} \mbox{for Guo, Luo models}
    
    \sum_q c_{qi} c_{qj} f_q = 0  \hspace{1cm} \mbox{for Simple, Buick}
    
Models with zero second order moment have:

.. math ::
    
    F_q = \frac{w_q}{c_s^2} c_{qi} \; a_i

Models with nonzero second moment have:

.. math ::
    
    F_q = \frac{w_q}{c_s^2} c_{qi} \; a_i + \frac{w_q}{c_s^4} (c_{qi} c_{qj} - c_s^2 \delta_{ij} ) u_j \, a_i


Martin Bauer's avatar
Martin Bauer committed
55
For all force models the computation of the macroscopic velocity has to be adapted (shifted) by adding a term
56
57
58
59
60
61
62
63
64
:math:`S_{macro}` that we call "macroscopic velocity shift"
    
    .. math ::
    
        \pmb{u} = \sum_q \pmb{c}_q f_q + S_{macro}
        
        S_{macro} = \frac{\Delta t}{2} \sum_q F_q


Martin Bauer's avatar
Martin Bauer committed
65
Some models also shift the velocity entering the equilibrium distribution.
66
67
68
69
70
71
72

Comparison
----------
 
Force models can be distinguished by 2 options:

**Option 1**:
Martin Bauer's avatar
Martin Bauer committed
73
74
    :math:`C_F = 1` and equilibrium velocity is not shifted, or :math:`C_F=(1 - \frac{\omega}{2})`
    and equilibrum is shifted.
75
76
77
78
79
80
81
82
83
84
85
86
 
**Option 2**: 
    second velocity moment is zero or :math:`F_i u_j + F_j u_i`


=====================  ====================  =================
 Option2 \\ Option1    no equilibrium shift  equilibrium shift
=====================  ====================  =================
second moment zero      :class:`Simple`      :class:`Buick`
second moment nonzero   :class:`Luo`         :class:`Guo`         
=====================  ====================  =================
  
Martin Bauer's avatar
Martin Bauer committed
87
"""
88

Martin Bauer's avatar
Martin Bauer committed
89
import sympy as sp
Martin Bauer's avatar
Martin Bauer committed
90

Michael Kuron's avatar
Michael Kuron committed
91
from lbmpy.relaxationrates import get_bulk_relaxation_rate, get_shear_relaxation_rate
Martin Bauer's avatar
Martin Bauer committed
92
93


94
class Simple:
95
    r"""
Martin Bauer's avatar
Martin Bauer committed
96
    A simple force model which introduces the following additional force term in the
97
    collision process :math:`\frac{w_q}{c_s^2} c_{qi} \; a_i` (often: force = rho * acceleration)
Martin Bauer's avatar
Martin Bauer committed
98
99
100
    Should only be used with constant forces!
    Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
    """
Frederik Hennig's avatar
Frederik Hennig committed
101

Martin Bauer's avatar
Martin Bauer committed
102
103
104
    def __init__(self, force):
        self._force = force

Martin Bauer's avatar
Martin Bauer committed
105
106
    def __call__(self, lb_method, **kwargs):
        assert len(self._force) == lb_method.dim
107

Martin Bauer's avatar
Martin Bauer committed
108
        def scalar_product(a, b):
109
110
            return sum(a_i * b_i for a_i, b_i in zip(a, b))

Martin Bauer's avatar
Martin Bauer committed
111
        return [3 * w_i * scalar_product(self._force, direction)
Martin Bauer's avatar
Martin Bauer committed
112
                for direction, w_i in zip(lb_method.stencil, lb_method.weights)]
Martin Bauer's avatar
Martin Bauer committed
113

Frederik Hennig's avatar
Frederik Hennig committed
114
115
116
    def moment_space_forcing(self, lb_method):
        return (lb_method.moment_matrix * sp.Matrix(self(lb_method))).expand()

Martin Bauer's avatar
Martin Bauer committed
117
    def macroscopic_velocity_shift(self, density):
Martin Bauer's avatar
Martin Bauer committed
118
        return default_velocity_shift(density, self._force)
Martin Bauer's avatar
Martin Bauer committed
119

120
121
122
    def macroscopic_momentum_density_shift(self, density):
        return default_momentum_density_shift(self._force)

Martin Bauer's avatar
Martin Bauer committed
123

124
class Luo:
Martin Bauer's avatar
Martin Bauer committed
125
    r"""Force model by Luo :cite:`luo1993lattice`.
Martin Bauer's avatar
Martin Bauer committed
126
127
128

    Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity.
    """
Frederik Hennig's avatar
Frederik Hennig committed
129

Martin Bauer's avatar
Martin Bauer committed
130
131
132
    def __init__(self, force):
        self._force = force

Martin Bauer's avatar
Martin Bauer committed
133
134
    def __call__(self, lb_method):
        u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)
Martin Bauer's avatar
Martin Bauer committed
135
136
137
        force = sp.Matrix(self._force)

        result = []
Martin Bauer's avatar
Martin Bauer committed
138
        for direction, w_i in zip(lb_method.stencil, lb_method.weights):
Martin Bauer's avatar
Martin Bauer committed
139
140
141
142
            direction = sp.Matrix(direction)
            result.append(3 * w_i * force.dot(direction - u + 3 * direction * direction.dot(u)))
        return result

Frederik Hennig's avatar
Frederik Hennig committed
143
144
145
    def moment_space_forcing(self, lb_method):
        return (lb_method.moment_matrix * sp.Matrix(self(lb_method))).expand()

Martin Bauer's avatar
Martin Bauer committed
146
    def macroscopic_velocity_shift(self, density):
Martin Bauer's avatar
Martin Bauer committed
147
        return default_velocity_shift(density, self._force)
Martin Bauer's avatar
Martin Bauer committed
148

149
150
151
    def macroscopic_momentum_density_shift(self, density):
        return default_momentum_density_shift(self._force)

Martin Bauer's avatar
Martin Bauer committed
152

153
class Guo:
154
    r"""
155
    Force model by Guo  :cite:`guo2002discrete`
Martin Bauer's avatar
Martin Bauer committed
156
157
    Adapts the calculation of the macroscopic velocity as well as the equilibrium velocity (both shifted by F/2)!
    """
Frederik Hennig's avatar
Frederik Hennig committed
158

159
    def __init__(self, force):
Martin Bauer's avatar
Martin Bauer committed
160
161
        self._force = force

Martin Bauer's avatar
Martin Bauer committed
162
    def __call__(self, lb_method):
Martin Bauer's avatar
Martin Bauer committed
163
        luo = Luo(self._force)
164

Martin Bauer's avatar
Martin Bauer committed
165
        shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
Frederik Hennig's avatar
Frederik Hennig committed
166
167
        assert len(set(lb_method.relaxation_rates)) == 1, (
            "In population space, guo only works for SRT, use Schiller instead")
Martin Bauer's avatar
Martin Bauer committed
168
169
        correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
        return [correction_factor * t for t in luo(lb_method)]
Martin Bauer's avatar
Martin Bauer committed
170

Frederik Hennig's avatar
Frederik Hennig committed
171
172
173
174
175
176
177
    def moment_space_forcing(self, lb_method):
        luo = Luo(self._force)
        q = len(lb_method.stencil)
        correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
        moments = correction_factor * (lb_method.moment_matrix * sp.Matrix(luo(lb_method)))
        return moments.expand()

Martin Bauer's avatar
Martin Bauer committed
178
    def macroscopic_velocity_shift(self, density):
Martin Bauer's avatar
Martin Bauer committed
179
        return default_velocity_shift(density, self._force)
Martin Bauer's avatar
Martin Bauer committed
180

181
182
183
    def macroscopic_momentum_density_shift(self, density):
        return default_momentum_density_shift(self._force)

Martin Bauer's avatar
Martin Bauer committed
184
185
    def equilibrium_velocity_shift(self, density):
        return default_velocity_shift(density, self._force)
Martin Bauer's avatar
Martin Bauer committed
186
187


Michael Kuron's avatar
Michael Kuron committed
188
189
class Schiller:
    r"""
190
    Force model by Schiller  :cite:`schiller2008thermal`, equation 4.67
Michael Kuron's avatar
Michael Kuron committed
191
192
    Equivalent to Guo but not restricted to SRT.
    """
Frederik Hennig's avatar
Frederik Hennig committed
193

Michael Kuron's avatar
Michael Kuron committed
194
195
196
197
198
199
    def __init__(self, force):
        self._force = force

    def __call__(self, lb_method):
        u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)
        force = sp.Matrix(self._force)
Frederik Hennig's avatar
Frederik Hennig committed
200

Michael Kuron's avatar
Michael Kuron committed
201
202
203
        uf = u.dot(force) * sp.eye(len(force))
        omega = get_shear_relaxation_rate(lb_method)
        omega_bulk = get_bulk_relaxation_rate(lb_method)
Michael Kuron's avatar
Michael Kuron committed
204
        G = (u * force.transpose() + force * u.transpose() - uf * sp.Rational(2, lb_method.dim)) * sp.Rational(1, 2) * \
205
            (2 - omega) + uf * sp.Rational(1, lb_method.dim) * (2 - omega_bulk)
Michael Kuron's avatar
Michael Kuron committed
206
207
208
209
210
211
212

        result = []
        for direction, w_i in zip(lb_method.stencil, lb_method.weights):
            direction = sp.Matrix(direction)
            tr = sp.trace(G * (direction * direction.transpose() - sp.Rational(1, 3) * sp.eye(len(force))))
            result.append(3 * w_i * (force.dot(direction) + sp.Rational(3, 2) * tr))
        return result
Frederik Hennig's avatar
Frederik Hennig committed
213

Michael Kuron's avatar
Michael Kuron committed
214
215
216
    def macroscopic_velocity_shift(self, density):
        return default_velocity_shift(density, self._force)

217
218
219
    def macroscopic_momentum_density_shift(self, density):
        return default_momentum_density_shift(self._force)

Michael Kuron's avatar
Michael Kuron committed
220

221
class Buick:
222
223
224
225
226
227
228
229
    r"""
    This force model :cite:`buick2000gravity` has a force term with zero second moment. 
    It is suited for incompressible lattice models.
    """

    def __init__(self, force):
        self._force = force

Martin Bauer's avatar
Martin Bauer committed
230
    def __call__(self, lb_method, **kwargs):
231
232
        simple = Simple(self._force)

Martin Bauer's avatar
Martin Bauer committed
233
        shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
234
        assert len(set(lb_method.relaxation_rates)) == 1, "Buick only works for SRT"
Martin Bauer's avatar
Martin Bauer committed
235
236
        correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
        return [correction_factor * t for t in simple(lb_method)]
237

Martin Bauer's avatar
Martin Bauer committed
238
    def macroscopic_velocity_shift(self, density):
Martin Bauer's avatar
Martin Bauer committed
239
        return default_velocity_shift(density, self._force)
240

241
242
243
    def macroscopic_momentum_density_shift(self, density):
        return default_momentum_density_shift(self._force)

Martin Bauer's avatar
Martin Bauer committed
244
245
    def equilibrium_velocity_shift(self, density):
        return default_velocity_shift(density, self._force)
246
247


248
class EDM:
Martin Bauer's avatar
Martin Bauer committed
249
250
251
252
253
    r"""Exact differencing force model"""

    def __init__(self, force):
        self._force = force

Martin Bauer's avatar
Martin Bauer committed
254
255
256
257
    def __call__(self, lb_method, **kwargs):
        cqc = lb_method.conserved_quantity_computation
        rho = cqc.zeroth_order_moment_symbol if cqc.compressible else 1
        u = cqc.first_order_moment_symbols
Martin Bauer's avatar
Martin Bauer committed
258

Martin Bauer's avatar
Martin Bauer committed
259
260
261
262
        shifted_u = (u_i + f_i / rho for u_i, f_i in zip(u, self._force))
        eq_terms = lb_method.get_equilibrium_terms()
        shifted_eq = eq_terms.subs({u_i: su_i for u_i, su_i in zip(u, shifted_u)})
        return shifted_eq - eq_terms
Martin Bauer's avatar
Martin Bauer committed
263

Martin Bauer's avatar
Martin Bauer committed
264
    def macroscopic_velocity_shift(self, density):
Martin Bauer's avatar
Martin Bauer committed
265
        return default_velocity_shift(density, self._force)
Martin Bauer's avatar
Martin Bauer committed
266

267
268
269
    def macroscopic_momentum_density_shift(self, density):
        return default_momentum_density_shift(self._force)

Martin Bauer's avatar
Martin Bauer committed
270

Martin Bauer's avatar
Martin Bauer committed
271
272
273
# --------------------------------  Helper functions  ------------------------------------------------------------------


Martin Bauer's avatar
Martin Bauer committed
274
def default_velocity_shift(density, force):
275
    return [f_i / (2 * density) for f_i in force]
276
277
278
279


def default_momentum_density_shift(force):
    return [f_i / 2 for f_i in force]