LBM section: first readable version

\input{content/meshes/appendix.tex}

\chapter{Velocity Sets for LBM}

\label{appendix:LBM:velocity sets}

A Python script solving the system of equations that the given velocity set must satisfy.

The reference mentioned in the code is \cite[Section~3.4.7.2]{kruger:2017lattice}.

\inputminted[

    linenos=true,

    fontsize=\scriptsize,

    % background color + frame

    bgcolor=codebg,

    frame=single,

    framesep=2\fboxsep,

    framerule=0.01pt,

]{python}{data/lbm/lattice_weights.py}

 No newline at end of file

#! /usr/bin/env python3

# Reference: [1] Timm Krüger - The Lattice Boltzmann Method, Section 3.4.7.2, Eq. (3.60)

import itertools

from sympy import *

init_printing()

# input parameters

D = 3    # space dimension

Q = 27   # number of discrete velocities

#cs = 1/Rational(2)   # speed of sound (D3Q7)

cs = 1/sqrt(3)        # speed of sound (D1Q3, D2Q5, D2Q9, D3Q7*, D3Q15, D3Q19, D3Q27)

# all discrete velocity vectors

c_1_3 = Matrix([[0, 1, -1]])

c_2_9 = Matrix([

    [0, 1, 0, -1, 0,   1, 1,-1,-1],

    [0, 0, 1, 0, -1,   1,-1, 1,-1],

])

c_3_27 = Matrix([

    [0, 1, 0, 0, -1, 0, 0,   0, 0, 0, 0, 1, 1,-1,-1, 1, 1,-1,-1,   1, 1, 1, 1,-1,-1,-1,-1],

    [0, 0, 1, 0, 0, -1, 0,   1, 1,-1,-1, 1,-1, 1,-1, 0, 0, 0, 0,   1, 1,-1,-1, 1, 1,-1,-1],

    [0, 0, 0, 1, 0, 0, -1,   1,-1, 1,-1, 0, 0, 0, 0, 1,-1, 1,-1,   1,-1, 1,-1, 1,-1, 1,-1],

])

# select the discrete velocity vectors for D and Q

if D == 1:

    c = c_1_3

elif D == 2:

    c = c_2_9.extract(range(D), range(Q))

elif D == 3 and Q == 15:

    c = c_3_27.extract(range(D), [i for i in range(27) if i < 7 or 19 <= i])

elif D == 3:

    c = c_3_27.extract(range(D), range(Q))

else:

    raise NotImplementedError

assert sum(c) == 0

print(f"Matrix of discrete velocity vectors for D{D}Q{Q}")

pprint(c)

def sym_vector(dim, name):

    syms = []

    for i in range(dim):

        syms.append(symbols(f"{name}_{i}"))

    return Matrix(syms)

def delta(a, b):

    return 1 if a == b else 0

# solve weights according to [1]

w = sym_vector(Q, "w")

eqs = []

# eq 1

eqs.append(Eq(sum(w), 1))

# eq 2

for d1 in range(D):

    eqs.append(Eq(w.dot(c.row(d1)), 0))

# eq 3

for d1, d2 in itertools.product(range(D), repeat=2):

    rhs = cs ** 2 if d1 == d2 else 0

    s = Add(0)

    for i in range(Q):

        s += w[i] * c[d1, i] * c[d2, i]

    # avoid adding identical equation "0=0" (it would add a boolean constant instead of an expression)

    if s != 0:

        eqs.append(Eq(s, rhs))

# eq 4

for d1, d2, d3 in itertools.product(range(D), repeat=3):

    s = Add(0)

    for i in range(Q):

        s += w[i] * c[d1, i] * c[d2, i] * c[d3, i]

    # avoid adding identical equation "0=0" (it would add a boolean constant instead of an expression)

    if s != 0:

        eqs.append(Eq(s, 0))

# avoid overdetermined system

if not (D == 3 and Q == 7 and cs == 1/Rational(2)):

    # eq 5

    for d1, d2, d3, d4 in itertools.product(range(D), repeat=4):

        rhs = cs ** 4 * (delta(d1,d2) * delta(d3,d4) + delta(d1,d3) * delta(d2,d4) \

                + delta(d1,d4) * delta(d2,d3))

        s = Add(0)

        for i in range(Q):

            s += w[i] * c[d1, i] * c[d2, i] * c[d3, i] * c[d4, i]

        # avoid adding identical equation "0=0" (it would add a boolean constant instead of an expression)

        if s != 0:

            eqs.append(Eq(s, rhs))

    # eq 6

    for d1, d2, d3, d4, d5 in itertools.product(range(D), repeat=5):

        s = Add(0)

        for i in range(Q):

            s += w[i] * c[d1, i] * c[d2, i] * c[d3, i] * c[d4, i] * c[d5, i]

        # avoid adding identical equation "0=0" (it would add a boolean constant instead of an expression)

        if s != 0:

            eqs.append(Eq(s, 0))

print("Total number of equations:", len(eqs))

for i, eq in enumerate(eqs):

    pprint(f"{i}-th equation: {eq}")

solution = linsolve(eqs, [wi for wi in w])

if not solution:

    raise Exception("linsolve did not return a solution")

solution = Matrix(list(solution)[0])

print("solution:")

pprint(Eq(w, solution))

# D3Q27 is underdetermined (see [1])

if D == 3 and Q == 27:

    w26 = 1/Rational(216)

    print(f"assuming {w[-1]}={w26}:")

    pprint(Eq(w, solution.subs({w[-1]: w26})))

    %\vspace{0.25\baselineskip}%

    \par\noindent%

}

% environment for algorithms

\usepackage{amsthm}

\usepackage{thmtools}

\declaretheoremstyle[

    spaceabove=\topsep,

    spacebelow=\topsep,

    headfont=\normalfont\bfseries,

    notefont=\mdseries,

    notebraces={(}{)},

    bodyfont=\normalfont,

    headpunct={},

    % newline has no effect on enumerate environment and conflicts with the \mbox{} workaround

%    postheadspace=\newline,

]{algstyle}

\declaretheorem[

    name=Algorithm,

    style=algstyle,

    refname=algorithm,

    Refname=Algorithm,

]{algorithm}

%\crefname{algorithm}{algorithm}{algorithms}

\crefname{algorithm}{Algorithm}{Algorithms}

\Crefname{algorithm}{Algorithm}{Algorithms}

\usepackage{enumitem}

% enum list for algorithm items (basic style, no special indentation)

\newlist{algitems}{enumerate}{3}

% enum list for algorithm steps (style to keep labels left-aligned and steps indented)

\newlist{algsteps}{enumerate}{3}

% style of all levels

\setlist*[algitems,algsteps]{

    itemsep=0ex,

    % force the first item on a new line

    before=\mbox{}\\[-\baselineskip],

    label*=\arabic*.,

}

% ref style of first level

\setlist*[algitems,algsteps,1]{ref=\arabic*}

% ref style of second level

\setlist*[algitems,2]{ref=\arabic{algitemsi}.\arabic*}

\setlist*[algsteps,2]{ref=\arabic{algstepsi}.\arabic*}

% ref style of third level

\setlist*[algitems,3]{ref=\arabic{algitemsi}.\arabic{algitemsii}.\arabic*}

\setlist*[algsteps,3]{ref=\arabic{algstepsi}.\arabic{algstepsii}.\arabic*}

% style of second and third levels

\setlist*[algitems,algsteps,2,3]{

    topsep=0ex,

}

% label alignment for the algsteps list

\setlist*[algsteps,1]{labelsep=1em}

\setlist*[algsteps,2]{labelsep=2.4em}

\setlist*[algsteps,3]{labelsep=3.6em, itemindent=-0.5em}

\makeatother

\begin{appendices}

    \input{content/meshes/appendix.tex}

    \input{content/appendix.tex}

\end{appendices}

% restore the ToC hierarchy