forcemodels.py 9.37 KB
 Martin Bauer committed Oct 17, 2017 1 r"""  Martin Bauer committed Feb 09, 2017 2 3 4 .. module:: forcemodels :synopsis: Collection of forcing terms for hydrodynamic LBM simulations  Martin Bauer committed Oct 17, 2017 5 6  Get started:  Martin Bauer committed Apr 30, 2018 7 ------------  Martin Bauer committed Oct 17, 2017 8   Martin Bauer committed Apr 13, 2018 9 This module offers different models to introduce a body force in the lattice Boltzmann scheme.  Michael Kuron committed Sep 04, 2020 10 If you don't know which model to choose, use :class:lbmpy.forcemodels.Schiller.  Martin Bauer committed Oct 17, 2017 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 committed Apr 13, 2018 29  \sum_q \pmb{c}_q F_q = \pmb{a}  Martin Bauer committed Oct 17, 2017 30 31   Martin Bauer committed Apr 13, 2018 32 The second order moment is different for the forcing models - if it is zero the model is suited for  Martin Bauer committed Oct 17, 2017 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 committed Apr 13, 2018 55 For all force models the computation of the macroscopic velocity has to be adapted (shifted) by adding a term  Martin Bauer committed Oct 17, 2017 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 committed Apr 13, 2018 65 Some models also shift the velocity entering the equilibrium distribution.  Martin Bauer committed Oct 17, 2017 66 67 68 69 70 71 72  Comparison ---------- Force models can be distinguished by 2 options: **Option 1**:  Martin Bauer committed Apr 13, 2018 73 74  :math:C_F = 1 and equilibrium velocity is not shifted, or :math:C_F=(1 - \frac{\omega}{2}) and equilibrum is shifted.  Martin Bauer committed Oct 17, 2017 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 committed Feb 09, 2017 87 """  Martin Bauer committed Oct 17, 2017 88   Martin Bauer committed Feb 09, 2017 89 import sympy as sp  Martin Bauer committed Jul 11, 2019 90   Michael Kuron committed Jun 25, 2020 91 from lbmpy.relaxationrates import get_bulk_relaxation_rate, get_shear_relaxation_rate  Martin Bauer committed Feb 09, 2017 92 93   Martin Bauer committed Jul 26, 2019 94 class Simple:  Martin Bauer committed Feb 09, 2017 95  r"""  Martin Bauer committed Feb 09, 2017 96  A simple force model which introduces the following additional force term in the  Martin Bauer committed Oct 17, 2017 97  collision process :math:\frac{w_q}{c_s^2} c_{qi} \; a_i (often: force = rho * acceleration)  Martin Bauer committed Feb 09, 2017 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 committed Aug 05, 2021 101   Martin Bauer committed Feb 09, 2017 102 103 104  def __init__(self, force): self._force = force  Martin Bauer committed Apr 10, 2018 105 106  def __call__(self, lb_method, **kwargs): assert len(self._force) == lb_method.dim  Martin Bauer committed Feb 13, 2017 107   Martin Bauer committed Apr 10, 2018 108  def scalar_product(a, b):  Martin Bauer committed Feb 13, 2017 109 110  return sum(a_i * b_i for a_i, b_i in zip(a, b))  Martin Bauer committed Apr 10, 2018 111  return [3 * w_i * scalar_product(self._force, direction)  Martin Bauer committed Apr 10, 2018 112  for direction, w_i in zip(lb_method.stencil, lb_method.weights)]  Martin Bauer committed Feb 09, 2017 113   Frederik Hennig committed Aug 05, 2021 114 115 116  def moment_space_forcing(self, lb_method): return (lb_method.moment_matrix * sp.Matrix(self(lb_method))).expand()  Martin Bauer committed Apr 10, 2018 117  def macroscopic_velocity_shift(self, density):  Martin Bauer committed Apr 10, 2018 118  return default_velocity_shift(density, self._force)  Martin Bauer committed Feb 21, 2017 119   Helen Schottenhamml committed Mar 29, 2021 120 121 122  def macroscopic_momentum_density_shift(self, density): return default_momentum_density_shift(self._force)  Martin Bauer committed Feb 09, 2017 123   Martin Bauer committed Jul 26, 2019 124 class Luo:  Martin Bauer committed Apr 10, 2018 125  r"""Force model by Luo :cite:luo1993lattice.  Martin Bauer committed Feb 09, 2017 126 127 128  Shifts the macroscopic velocity by F/2, but does not change the equilibrium velocity. """  Frederik Hennig committed Aug 05, 2021 129   Martin Bauer committed Feb 09, 2017 130 131 132  def __init__(self, force): self._force = force  Martin Bauer committed Apr 10, 2018 133 134  def __call__(self, lb_method): u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)  Martin Bauer committed Feb 09, 2017 135 136 137  force = sp.Matrix(self._force) result = []  Martin Bauer committed Apr 10, 2018 138  for direction, w_i in zip(lb_method.stencil, lb_method.weights):  Martin Bauer committed Feb 09, 2017 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 committed Aug 05, 2021 143 144 145  def moment_space_forcing(self, lb_method): return (lb_method.moment_matrix * sp.Matrix(self(lb_method))).expand()  Martin Bauer committed Apr 10, 2018 146  def macroscopic_velocity_shift(self, density):  Martin Bauer committed Apr 10, 2018 147  return default_velocity_shift(density, self._force)  Martin Bauer committed Feb 09, 2017 148   Helen Schottenhamml committed Mar 29, 2021 149 150 151  def macroscopic_momentum_density_shift(self, density): return default_momentum_density_shift(self._force)  Martin Bauer committed Feb 09, 2017 152   Martin Bauer committed Jul 26, 2019 153 class Guo:  Martin Bauer committed Feb 09, 2017 154  r"""  Martin Bauer committed Oct 17, 2017 155  Force model by Guo :cite:guo2002discrete  Martin Bauer committed Feb 09, 2017 156 157  Adapts the calculation of the macroscopic velocity as well as the equilibrium velocity (both shifted by F/2)! """  Frederik Hennig committed Aug 05, 2021 158   Martin Bauer committed Feb 09, 2017 159  def __init__(self, force):  Martin Bauer committed Feb 09, 2017 160 161  self._force = force  Martin Bauer committed Apr 10, 2018 162  def __call__(self, lb_method):  Martin Bauer committed Feb 09, 2017 163  luo = Luo(self._force)  Martin Bauer committed Mar 05, 2017 164   Martin Bauer committed Apr 10, 2018 165  shear_relaxation_rate = get_shear_relaxation_rate(lb_method)  Frederik Hennig committed Aug 05, 2021 166 167  assert len(set(lb_method.relaxation_rates)) == 1, ( "In population space, guo only works for SRT, use Schiller instead")  Martin Bauer committed Apr 10, 2018 168 169  correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate) return [correction_factor * t for t in luo(lb_method)]  Martin Bauer committed Feb 09, 2017 170   Frederik Hennig committed Aug 05, 2021 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 committed Apr 10, 2018 178  def macroscopic_velocity_shift(self, density):  Martin Bauer committed Apr 10, 2018 179  return default_velocity_shift(density, self._force)  Martin Bauer committed Feb 09, 2017 180   Helen Schottenhamml committed Mar 29, 2021 181 182 183  def macroscopic_momentum_density_shift(self, density): return default_momentum_density_shift(self._force)  Martin Bauer committed Apr 10, 2018 184 185  def equilibrium_velocity_shift(self, density): return default_velocity_shift(density, self._force)  Martin Bauer committed Feb 09, 2017 186 187   Michael Kuron committed Jun 25, 2020 188 189 class Schiller: r"""  Michael Kuron committed Sep 04, 2020 190  Force model by Schiller :cite:schiller2008thermal, equation 4.67  Michael Kuron committed Jun 25, 2020 191 192  Equivalent to Guo but not restricted to SRT. """  Frederik Hennig committed Aug 05, 2021 193   Michael Kuron committed Jun 25, 2020 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 committed Aug 05, 2021 200   Michael Kuron committed Jun 25, 2020 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 committed Sep 04, 2020 204  G = (u * force.transpose() + force * u.transpose() - uf * sp.Rational(2, lb_method.dim)) * sp.Rational(1, 2) * \  Michael Kuron committed Sep 07, 2020 205  (2 - omega) + uf * sp.Rational(1, lb_method.dim) * (2 - omega_bulk)  Michael Kuron committed Jun 25, 2020 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 committed Aug 05, 2021 213   Michael Kuron committed Jun 25, 2020 214 215 216  def macroscopic_velocity_shift(self, density): return default_velocity_shift(density, self._force)  Helen Schottenhamml committed Mar 29, 2021 217 218 219  def macroscopic_momentum_density_shift(self, density): return default_momentum_density_shift(self._force)  Michael Kuron committed Jun 25, 2020 220   Martin Bauer committed Jul 26, 2019 221 class Buick:  Martin Bauer committed Oct 17, 2017 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 committed Apr 10, 2018 230  def __call__(self, lb_method, **kwargs):  Martin Bauer committed Oct 17, 2017 231 232  simple = Simple(self._force)  Martin Bauer committed Apr 10, 2018 233  shear_relaxation_rate = get_shear_relaxation_rate(lb_method)  Michael Kuron committed Jun 26, 2020 234  assert len(set(lb_method.relaxation_rates)) == 1, "Buick only works for SRT"  Martin Bauer committed Apr 10, 2018 235 236  correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate) return [correction_factor * t for t in simple(lb_method)]  Martin Bauer committed Oct 17, 2017 237   Martin Bauer committed Apr 10, 2018 238  def macroscopic_velocity_shift(self, density):  Martin Bauer committed Apr 10, 2018 239  return default_velocity_shift(density, self._force)  Martin Bauer committed Oct 17, 2017 240   Helen Schottenhamml committed Mar 29, 2021 241 242 243  def macroscopic_momentum_density_shift(self, density): return default_momentum_density_shift(self._force)  Martin Bauer committed Apr 10, 2018 244 245  def equilibrium_velocity_shift(self, density): return default_velocity_shift(density, self._force)  Martin Bauer committed Oct 17, 2017 246 247   Martin Bauer committed Jul 26, 2019 248 class EDM:  Martin Bauer committed Feb 08, 2018 249 250 251 252 253  r"""Exact differencing force model""" def __init__(self, force): self._force = force  Martin Bauer committed Apr 10, 2018 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 committed Feb 08, 2018 258   Martin Bauer committed Apr 10, 2018 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 committed Feb 08, 2018 263   Martin Bauer committed Apr 10, 2018 264  def macroscopic_velocity_shift(self, density):  Martin Bauer committed Apr 10, 2018 265  return default_velocity_shift(density, self._force)  Martin Bauer committed Feb 08, 2018 266   Helen Schottenhamml committed Mar 29, 2021 267 268 269  def macroscopic_momentum_density_shift(self, density): return default_momentum_density_shift(self._force)  Martin Bauer committed Feb 08, 2018 270   Martin Bauer committed Feb 09, 2017 271 272 273 # -------------------------------- Helper functions ------------------------------------------------------------------  Martin Bauer committed Apr 10, 2018 274 def default_velocity_shift(density, force):  Martin Bauer committed Feb 09, 2017 275  return [f_i / (2 * density) for f_i in force]  Helen Schottenhamml committed Mar 29, 2021 276 277 278 279  def default_momentum_density_shift(force): return [f_i / 2 for f_i in force]