\section{Einführung}

\subsection{Warum Gitterverfeinerung?}

Die einfachsten und zugleich am weitesten verbreiteten Umsetzungen von Simulationen mit Lattice Boltzmann Methoden basieren auf uniformen Gittern, in denen Zellen immer den gleichen Abstand zu ihren Nachbarzellen haben.

Die Genauigkeit von Lattice Boltzmann basierenden Simulationen hängt maßgeblich von der Auflösung des verwendeten Gitters ab. Bei Außerachtlassung weiterer wichtiger Faktoren wie dem verwendeten Kollisionsterm und Randkonditionen kann im Allgemeinen davon ausgegangen werden, dass eine feinere Auflösung des Gitters zu besseren Ergebnissen führt.

In praktischen Beispielen können innerhalb eines Modells große Unterschiede in der Strömungskomplexität existieren. So kann es große Gebiete eines Modells geben, die mit einem vergleichsweise groben Gitter gut simuliert werden können, während in anderen Gebieten \textendash{} beispielsweise in komplexen Geometrien und an Rändern \textendash{} ein vielfach feineres Gitter zur adäquaten Behandlung benötigt wird. In uniformen Gittern muss jedoch das gesamte Modell unabhängig der lokalen Situation mit der maximal benötigten Auflösung abgebildet werden.

Da die Anzahl der benötigten Gitterpunkte sich maßgeblich auf den Speicherbedarf und Rechenaufwand auswirkt, ist es wünschenswert die Anzahl der Gitterpunkte zu minimieren. Ein Ansatz, dies zu erreichen, ist die lokale Variation der Gitterauflösung.

\newpage
\section{Grundlagen}

In diesem Kapitel werden wir die, dem weiteren Verlauf dieser Arbeit zugrunde liegende, Lattice Boltzmann Methode in 2D nachvollziehen.

\subsection{Lattice Boltzmann Methode}
\label{kap:LBM}

Grundlage und Namensgeber von Simulationen mit Lattice Boltzmann Methoden ist die Boltzmann Gleichung. Sie beschreibt das Verhalten von Gasen auf mesoskopischer Ebene als Verteilungsfunktion der Masse von Partikeln in einer Raumregion mit gegebener Geschwindigkeit.

\begin{Definition}[Die Boltzmann-Gleichung]
Sei \(f(x,\xi,t)\) die Verteilungsfunktion der Partikelmasse zu Zeit \(t\) in Ort \(x \in \R^2\) mit Geschwindigkeit \(\xi \in \R^2\), \(\rho\) die Dichte und \(F \in \R^2\) eine etwaige äußere Kraft. Die Boltzmann-Gleichung beschreibt die zeitliche Veränderung der Verteilungsfunktion anhand des totalen Differential \(\Omega(f)\):
\[ \left( \partial_t + \xi \cdot \partial_x + \frac{F}{\rho} \cdot \partial_\xi \right) f = \Omega(f) \left( = \partial_x f \cdot \frac{dx}{dt} + \partial_\xi f \cdot \frac{d\xi}{dt} + \partial_t f \right) \]

Hierbei handelt es sich um eine Advektionsgleichung wobei der Term \(\partial_t f + \xi \cdot \partial_x f\) die Strömung der Partikelverteilung mit Geschwindigkeit \(\xi\) und \(\frac{F}{\rho} \cdot \partial_\xi f\) einwirkende Kräfte darstellt. Der Term \(\Omega(f)\) beschreibt, entsprechend als Kollisionsoperator bezeichnet, die kollisionsbedingte lokale Neuverteilung von \(f\).
\end{Definition}

Zentrale Anforderung an den Kollisionsoperator ist die Impuls- und Masseerhaltung. Die im Folgenden betrachtete Lattice Boltzmann Methode verwendet die übliche BGK Approximation der Boltzmann-Gleichung ohne äußere Kraft von Bhatnagar, Gross und Krook (siehe \citetitle{krueger17}~\cite[Kap.~3.5.3]{krueger17}).
Grundlegendes Element dieser Approximation ist der BGK Operator
\[\Omega(f) := -\frac{f-f^\text{eq}}{\tau} \Delta t\]
welcher die Partikelverteilung mit Rate \(\tau\) gegen eine Equilibriumsverteilung \(f^\text{eq}\) relaxiert. Ohne Beschränkung der Allgemeinheit setzen wir dabei  im Folgenden \(\Delta t = 1\).

Wenden wir den BGK Operator auf die Boltzmann-Gleichung ohne äußere Kräfte an, erhalten wir die BGK Approximation:

\begin{Definition}[BGK Approximation]
Sei \(\tau \in \R_{\geq 0}\) eine Relaxionszeit und \(f^\text{eq}\) die von der Maxwell-Boltzmann-Verteilung gegebene Equilibriumsverteilung.
\[ (\partial_t + \xi \cdot \nabla_x) f = -\frac{1}{\tau} (f(x,\xi,t) - f^\text{eq}(x,\xi,t)) \]
\end{Definition}

Analog zur Boltzmann-Gleichung ist auch bei deren BGK Approximation der beschriebene Ort \(x \in \R^2\) im Allgemeinen frei gewählt. Da unser Ziel jedoch gerade die Diskretisierung der Simulationsdomäne in einem Gitter ist, wollen wir \(x\) einschränken:

\begin{Definition}[Ortsdiskretisierung]
\label{def:SpatialDiscretizationLBM}
Sei die örtliche Simulationsdomäne \(D := \R^2\) diskretisiert als kartesisches Gitter mit Zellabstand \(\Delta x \in \R_+\). Dann ist die Einbettung der kartesischen Gitterdomäne \(L := \Z^2\) in \(D\) gegeben als:
\[d : L \to D, \ l \mapsto \Delta x \cdot l\]
Analog zur o.B.d.A. erfolgten Wahl von \(\Delta t = 1\) setzen wir auch hier \(\Delta x = 1\), sodass wir die Simulations- und Gitterdomäne -- bei Inklusion stetiger Übergänge zwischen Elementen aus \(L\) -- identisch identifizieren können. Der stetige Übergang zwischen zwei orthodonal benachbarten Gitterknoten erfolgt somit auf dem Einheitsintervall. Praktisch können wir also im Folgenden bei Annahme von \(x \in L \subset D\) die explizite Ausführung der Gittereinbettung vernachlässigen.
\end{Definition}

Zu erwähnen ist an dieser Stelle die Wichtigkeit einer klaren Unterscheidung zwischen Simulationsdomäne und dem durch diese zu modellierenden physikalischen System für die konkrete Interpretation des Simulationsergebnisses \cite[Kap.~7]{krueger17}.

Für die verbleibende Herleitung der LBM können wir diese Interpretation, d.h. die Unterscheidung zwischen physikalischen- und Lattice-Einheiten, jedoch außer Acht lassen, da sich die modellierten physikalischen Einheiten aus der Wahl der Relaxionszeit und der Skalierung der Lattice-Momente ergeben. Diese Wahl wird hier nicht weiter eingeschränkt.

\bigskip

Wir bemerken nun, dass die BGK Approximation nicht nur für beliebige Orte, sondern auch für beliebige Geschwindigkeiten \(\xi \in \R^2\) definiert ist. Da wir die LBM auf einem endlichen Rechner umsetzen wollen, müssen wir auch die Menge der betrachteten Geschwindigkeiten auf eine endliche Menge diskretisieren.

Eine übliche Menge diskreter Geschwindigkeiten in 2D ist \emph{D2Q9} wobei \emph{D2} die Anzahl der Dimensionen und \emph{Q9} die Anzahl der Geschwindigkeiten verschlüsselt.

\begin{Definition}[D2Q9 Modell]
\[ \{\xi_i\}_{i=0}^8 = \left\{ \V{0}{0}, \V{-1}{\phantom{-}1}, \V{-1}{\phantom{-}0}, \V{-1}{-1}, \V{\phantom{-}0}{-1}, \V{\phantom{-}1}{-1}, \V{1}{0}, \V{1}{1}, \V{0}{1} \right\} \]
\end{Definition}

\begin{figure}
\centering
\input{img/d2q9.tikz}
\caption{Umgebung einer Zelle in D2Q9}
\end{figure}

Mithilfe einer solchen endlichen Menge diskreter Geschwindigkeiten lässt sich die BGK Approximation bezüglich der Geschwindigkeit diskretisieren:

\begin{Definition}[BGK Geschwindigkeitsdiskretisierung]
\label{def:disVelBGK}
Seien \(\xi_i\) Vektoren einer Menge mikroskopischer Geschwindigkeiten wie z.B. D2Q9 und \(f_i(x,t) \equiv f(x,\xi_i,t)\). Dann ist
\[ (\partial_t + \xi_i \cdot \nabla_x) f_i(x,t) = -\frac{1}{\tau} (f_i(x,t) - f_i^\text{eq}(x,t)) \]
die Diskretisierung der BGK Approximation entlang der Geschwindigkeiten in diskreten Gitterknoten \(x \in L\). Die Geschwindigkeiten müssen hier dank der Wahl von \(\Delta x = 1\) nicht weiter skaliert werden.
\end{Definition}

Hierbei ist die diskrete Equilibriumsverteilung \(f_i^\text{eq}\) wie folgt definiert:

\begin{Definition}[Diskrete Equilibriumsverteilung]
\label{def:fieq}
Seien \(\rho \in \R_{\geq 0}\) die Dichte, \(u \in \R^2\) die Gesamtgeschwindigkeit, \(\xi_i\) die \(i\)-te diskrete Geschwindigkeitskomponente, \(w_i\) das Gewicht jener Komponente bzgl. des Lattice und \(c_s\) die Lattice-Schallgeschwindigkeit.
\[f_i^\text{eq} = w_i \rho \left( 1 + \frac{u \cdot \xi_i}{c_s^2} + \frac{(u \cdot \xi_i)^2}{2c_s^4} - \frac{u \cdot u}{2c_s^2} \right)\]
\end{Definition}

Die Werte von \(u = u(x,t)\) und \(\rho = \rho(x,t)\) in Ort \(x\) zu Zeit \(t\) ergeben sich dabei aus den \emph{Momenten} der Verteilungsfunktion \(f_i\):

\begin{Definition}[Momente der Verteilungsfunktion]
\label{def:Momente}
\begin{align*}
\rho(x,t) &= \sum_{i=0}^{q-1} f_i(x,t) \\
\rho u(x,t) &= \sum_{i=0}^{q-1} \xi_i f_i(x,t)
\end{align*}
\end{Definition}

Für D2Q9 ergeben sich nach \cite[Gl.~3.60 bzw. Tab.~3.3]{krueger17} die Gewichte:
\[w_0 = \frac{4}{9}, \ w_{2,4,6,8} = \frac{1}{9}, \ w_{1,3,5,7} = \frac{1}{36}\]

Weiter folgt zusammen mit der Bedingung \(\sum_{i=1}^{q-1} w_i (\xi_i)_a (\xi_i)_b = c_s^2 \delta_{a,b}\) aus \cite[Gl.~3.60]{krueger17} die Schallgeschwindigkeit \(c_s = \sqrt{1/3}\) des Lattice. Konditionen zur Bestimmung dieser gitterspezifischen Konstanten sind hierbei die Erhaltung von Impuls und Masse sowie die Forderung von \emph{Rotationsisotropie}.

Zur Entwicklung einer \emph{implementierbaren} expliziten BGK Gleichung können wir nun die Geschwindigkeitsdiskretisierung~\ref{def:disVelBGK} integrieren:
\[ f_i(x+\xi_i, t+1) - f_i(x,t) = \int_0^1 \Omega_i(x+\xi_i s,t+s) ds \]

Wobei \(\Omega_i(x,t)\) hier die diskrete Formulierung des BGK Kollisionsoperators darstellt:
\[\Omega_i(x,t) := -\frac{1}{\tau} ( f_i(x,t) - f_i^\text{eq}(x, t) )\]

Da sich die exakte Lösung des Integrals auf der rechten Seite schwierig gestaltet, wird es in der Praxis nur approximiert. Während es dazu vielfältige Ansätze gibt, beschränken wir uns an dieser Stelle auf Anwendung der Trapezregel:
\begin{align*}
f_i(x+\xi_i,t+1) - f_i(x,t) &= \frac{1}{2} \left( \Omega_i(x,t) + \Omega_i(x+\xi_i,t+1) \right) \\
&= -\frac{1}{2\tau} \left( f_i(x+\xi_i,t+1) + f_i(x,t) - f_i^\text{eq}(x+\xi_i,t+1) - f_i^\text{eq}(x,t) \right)
\end{align*}

Zur expliziten Lösung dieser impliziten Gleichung benötigen wir nun nur noch eine geeignete Verschiebung von \(f_i\) und \(\tau\):

\begin{Definition}[Diskrete LBM BGK Gleichung]
\label{def:LBGKeq}
Seien \(\overline{f_i}\) und \(\overline\tau\) definiert:
\begin{align*}
\overline{f_i} &= f_i + \frac{1}{2\tau}(f_i - f_i^\text{eq}) \\
\overline\tau &= \tau + \frac{1}{2}
\end{align*}

Setzen wir diese verschobenen Variablen in das Ergebnis der Trapezregel ein, erhalten \cite[Kap.~A.5 mit \(\Delta t=1\)]{krueger17} wir die die vollständig diskretisierte LBM BGK Gleichung:

\[\overline{f_i}(x+\xi_i,t+1) = \overline{f_i}(x,t) - \frac{1}{\overline\tau} (\overline{f_i}(x,t) - f_i^\text{eq}(x,t))\]
\end{Definition}

Bemerkenswert ist an dieser Stelle, dass die Momente der Verteilungen mit \(\overline{f_i}\) analog zu \ref{def:Momente} berechnet werden können:
\begin{align*}
\sum_{i=0}^{q-1} \overline{f_i} = \sum_{i=0}^{q-1} f_i + \frac{1}{2\tau} \sum_{i=0}^{q-1} (f_i - f_i^\text{eq}) &= \rho \\
\sum_{i=0}^{q-1} \xi_i \overline{f_i} = \sum_{i=0}^{q-1} \xi_i f_i + \frac{1}{2\tau} \sum_{i=0}^{q-1} (f_i - f_i^\text{eq}) &= \rho u
\end{align*}

\newpage
\subsubsection{Algorithmus}
\label{kap:LBMimpl}

Bei der Implementierung der {diskreten LBM BGK Gleichung}~\ref{def:LBGKeq} auf einem Computer ist die Aufteilung in Kollisions- und Strömungsschritt üblich.

\begin{Definition}[Kollisionsschritt]
Annäherung der Verteilungsfunktion an die lokal berechnete Equilibriumsverteilung entsprechend dem BGK Kollisionsoperator.
\[f_i^\text{out}(x,t) = f_i(x,t) - \frac{1}{\tau}(f_i(x,t) - f_i^\text{eq}(x,t))\]
\end{Definition}

\begin{Definition}[Strömungsschritt]
Strömen der neuen Verteilungen auf die benachbarten Zellen entsprechend der jeweiligen diskreten Geschwindigkeit.
\[f_i(x+\xi_i,t+1) = f_i^\text{out}(x,t)\]
\end{Definition}

Bemerkenswert ist hierbei, dass der Kollisionsschritt nur lokale Informationen der jeweiligen Zelle benötigt und sich somit sehr gut zur parallelen Verarbeitung eignet.

\subsubsection{Chapman-Enskog Analyse}

Ziel der beschriebenen Lattice Boltzmann Methode ist die möglichst gute Approximation der schwach-kompressiblen Navier-Stokes Gleichungen auf der Simulationsdomäne.

\begin{Definition}[Schwach-kompressible Navier-Stokes Gleichungen]
Sei \(\rho\) die Dichte, \(u\) die Geschwindigkeit und \(p\) der Druck zu Zeit \(t\) sowie \(\nu\) die kinematische Viskosität und \(\mathrm{S}\) der Verzerrungstensor.
\begin{align*}
\partial_t  \rho + \nabla \cdot (\rho u) &= 0 \\
\partial_t u + (u \cdot \nabla) u &= -\frac{1}{\rho} \nabla p + 2\nu\nabla \cdot (\mathrm{S})
\end{align*}

Dabei sind Druck \(p\), kinetische Viskosität \(\nu\) und Verzerrungstensor \(\mathrm{S}\) definiert als:
\begin{align*}
p &= c_s^2 \rho \\
\nu &= c_s^2 \tau \\
\mathrm{S} &= \frac{1}{2} (\nabla u + (\nabla u)^\top) \numberthis\label{eq:Verzerrungstensor}
\end{align*}
\end{Definition}

Nach \cite[Kap.~4.1]{krueger17} kann die asymptotische Äquivalenz von LBM BGK Gleichung und schwach-kompressiblen Navier-Stokes Gleichungen mit der Entwicklung von Chapman-Enskog gezeigt werden.

\begin{Definition}[Chapman-Enskog Ansatz]
\label{def:ChapmanEnskog}
Der Chapman-Enskog Ansatz besteht in der Annahme, dass die Verteilungsfunktion \(f_i\) als leicht gestörte Equilibriumsverteilung dargestellt werden kann: \[f_i = f_i^\text{eq} + \epsilon f_i^{(1)} + \mathcal{O}(\epsilon^2)\]
Hierbei ist \(\epsilon f_i^{(1)}\) mit \(\epsilon \ll 1\) der Störterm 1. Ordnung. Dieser ist gegeben als:
\[\epsilon f_i^{(1)} = \frac{w_i}{2 c_s^4} \mathrm{Q}_i : \mathrm{\Pi}^{(1)} \numberthis\label{eq:firstOrderPertubation}\]
Wobei der Geschwindigkeitstensor \(\mathrm{Q}_i\) und das Störmoment \(\mathrm\Pi^{(1)}\) nach \cite[Kap.~4.1.3]{krueger17} definiert sind als:
\begin{align*}
\mathrm{Q}_i &= \xi_i \xi_i - c_s^2 I \numberthis\label{eq:velocityTensor} \\
\mathrm\Pi^{(1)} &= \sum_{i=0}^{q-1} \xi_i \xi_i \epsilon f_i^{(1)} = -2 c_s^2 \rho \overline\tau \mathrm{S} \numberthis\label{eq:pertubationMoment}
\end{align*}
\end{Definition}

\begin{Definition}[Nicht-Equilibriumsverteilung]
Die Relaxion der Verteilungsfunktion \(f_i\) gegen die Equilibriumsverteilung \(f_i^\text{eq}\) impliziert eine Dekomposition in Equilibriums- und Nicht-Equilibriumsverteilung: \[f_i := f_i^\text{eq} + f_i^\text{neq}\]
\end{Definition}

Unter Vernachlässigung von Störtermen der Ordnung \(\mathcal{O}(\epsilon^2)\) ergibt sich eine Näherung der Nicht-Equilibriumsverteilung:
\[f_i^\text{neq} \cong \epsilon f_i^{(1)}\]

Diese Darstellung können wir unter Verwendung von Definition~\ref{def:ChapmanEnskog} ausführen als:
\[f_i^\text{neq} \cong -\frac{w_i \rho \overline\tau}{c_s^2} \mathrm{Q}_i : \mathrm{S} \numberthis\label{eq:approxFneq}\]

\newpage
\subsection{Herangehensweisen an Gitterverfeinerung}

Grundsätzlich existieren mit der Multi-Grid bzw. Multi-Domain Herangehenweise zwei verschiedene Ansätze für Gitterverfeinerung in Lattice Boltzmann Methoden \cite[Kap.~3.1]{lagrava12}. Im Westenlichen unterscheiden die Ansätze sich in der Struktur der variabel aufgelösten Teilgitter der Simulationsdomäne. Weitere Unterschiede folgen dann aus dieser grundlegenden Struktur.

\subsubsection{Multi-Grid Ansatz}

Bei Verfahren des Multi-Grid Ansatzes existiert das gröbste Gitter über der gesamten Domäne. Feinere Teilgitter werden über gröberen Gittern platziert und nicht aus deren Verarbeitung ausgeschlossen. Somit existieren über der gesamten \emph{Fläche} feinerer Gitter auch die Knoten gröberer Gitter.

\begin{figure}[h]
\centering
\input{img/multi_grid.tikz}
\caption{Teilgitter in der Multi-Grid Herangehensweise}
\end{figure}

Vorteil dieser Herangehensweise ist es, dass feinere Teilgitter im Verlauf der Simulation ohne aufwendige Restrukturierung verschoben werden können, um beispielsweise komplexere Strömungsstrukturen zu \emph{verfolgen}. Nachteil ist, dass das Einsparpotentiale in Speicher- und Rechenbedarf durch Mehrfachabdeckung von Teilen der Simulationsdomäne mittels mehrerer Gitter nicht voll ausgenutzt werden können.

\subsubsection{Multi-Domain Ansatz}

Kern des Multi-Domain Ansatzes ist es, außerhalb von etwaigen verfahrensbedingten Übergangsbereichen, jede Position in der Simulationsdomäne durch genau ein Teilgitter abzubilden. Konkret werden also bereits durch feinere Gitter abgedeckte Bereiche aus gröberen Teilgittern ausgespart.

\begin{figure}[h]
\centering
\input{img/multi_domain_plain.tikz}
\caption{Teiligitter in der Multi-Domain Herangehensweise}
\end{figure}

\noindent
Vorteil gegenüber des Multi-Grid Ansatzes ist hier der weiter reduzierte Speicherbedarf sowie erwartete Einsparungen in der benötigten Rechenzeit. Erkauft werden diese Vorteile durch aufwendigere Kopplung \cite[Kap.~3.1]{lagrava12} der verschiedenen Teilgitter in den Übergangsbereichen.

\newpage
\section{Verfeinerungsmethode nach Lagrava et al.}

Wie in Kapitel~\ref{sec:olbRefinementChoice} noch näher begründet werden wird, bieten sich der Multi-Domain Ansatz als Grundlage für Gitterverfeinerung in OpenLB an. Passend zu dieser Wahl sowie der, im Rahmen dieser Arbeit getroffenen, Einschränkung auf zweidimensionale LBM mit BGK-Kollisionsoperator haben Lagrava et al. 2012 in \citetitle{lagrava12}~\cite{lagrava12} ein solches Verfeinerungsverfahren entwickelt. Die Stuktur dieses Verfahrens, mit potenziell austauschbaren Restriktions- und Interpolationsverfahren im zentralen Kopplungsschritt, erscheint dabei sogleich auch als Fundament eines Multi-Domain Gitterverfeinerungsframeworks in OpenLB.

\subsection{Übersicht}

Das Verfahren basiert auf der Multi-Domain Herangehensweise \cite[Kap.~3.1]{lagrava12} , bei welcher die feiner aufgelösten Teilbereiche der Simulationsdomäne aus dem gröber aufgelösten Gitter ausgeschlossen werden und sich nur in Übergangsbereichen überlappen.

\begin{figure}[h]
\centering
\input{img/multi_domain.tikz}
\caption{Multi-Domain Herangehensweise mit Übergangsbereich \cite[vgl.~Abb.~3]{lagrava12}}
\label{fig:MultiDomainOverlap}
\end{figure}

In diesen Übergangsbereichen, welche eine Breite von mindestens einer Einheit des gröberen Zellabstands haben, liegt die Hauptarbeit des Verfeinerungsverfahrens.

Während der Übergang vom feinen zum groben Gitter sich im Wesentlichen auf eine skalierte und gefilterte Restriktion der Verteilungen beschränkt, gestaltet sich der Übergang vom groben zum feinen Gitter aufwendiger, da feine Knoten, für deren Position kein grober Knoten existiert, aus den übrigen Daten interpoliert werden müssen.

Entsprechend liegt der Fokus des von Lagrava et al. entwickelten Algorithmus auf der Auswahl des Interpolationsverfahrens sowie der Skalierung der physikalischen Werte zwischen den unterschiedlich aufgelösten Verteilungen. Um diese Kopplung der verschiedenen Gitterauflösungen theoretisch erfassen zu können, müssen wir zunächst die Gitter selbst konkreter definieren:

\begin{Definition}[Diskretisierung der Gitter]
\label{def:DiskretRefinedGitter}
Wir betrachten kartesische Gitter \(\G\) und \(\F\) als Diskretisierung der physikalischen Simulationsdomänen \(D_g\) bzw. \(D_f\). Diese seien so gewählt, dass sie gerade die konvexen Hüllen ihrer Diskretisierungsgitter darstellen.
\begin{align*}
\G &\subset D_g \cap \{ x \in \R^2 | \exists i \in \Z^2 : x = \delta x_g \cdot i \} && \text{Gröberes Gitter} \\
\F &\subset D_f \cap \{ x \in \R^2 | \exists i \in \Z^2 : x = \delta x_f \cdot i \} && \text{Feineres Gitter}
\end{align*}
\(\delta x_g = \delta x_g / 2 \in \R_+\) seien die Diskretisierungsauflösungen im Verhältnis \(1:2\).
\end{Definition}

Anders als noch in Definition~\ref{def:SpatialDiscretizationLBM} betrachten wir die Gitter jetzt also nicht mehr unabhängig des darzustellenden physikalischen Modells, sondern unterscheiden anhand der physikalischen Auflösung \(\delta x\). Während im Kontext der LBM an sich weiterhin für beide Gitter \(\Delta x = 1\) gesetzt wird, führt die Relation von \(\delta x_g\) und \(\delta x_f\) im kommenden Kapitel~\ref{kap:Skalierung} u.a. zu einer Relation zwischen grober und feiner Relaxionszeit.

Eine stringente Behandlung von Gitterkopplung in diesem Modell benötigt Abbildungen der physikalisch eingebetteten Knoten aus \(\G\) und \(\F\) in die zugehörigen \emph{implementierenden} Gitter mit uniformer Auflösung \(\Delta x = 1\). Diese stellen gerade die Definitionsmengen der groben bzw. feinen Verteilungsfunktionen dar.

\begin{Definition}[Abbildung auf implementierende Gitter]
\label{def:BijImplGitter}
Sei \(\# \in \{f, g\}\) Symbol des feinen oder groben Gitters.
Dann können wir o.B.d.A. einen beliebigen physikalischen Knoten \(x_{0,\#}^\text{phys}\) mit dem Knoten \(x_0^\text{impl} = 0 \in L\) identifizieren. Eine Bijektion zwischen physikalischen und implementierenden Gittern ist damit schon eindeutig definiert:
\begin{align*}
x_\#^\text{impl}(x^\text{phys}) &:= \frac{1}{\delta x_\#} (x^\text{phys} - x_{0,\#}^\text{phys}) \\
x_\#^\text{phys}(x^\text{impl}) &:= x_{0,\#}^\text{phys} + \delta x_\# \cdot x^\text{impl}
\end{align*}
Diese Abbildung der physikalisch eingebetteten Gitterknoten in die Definitionsmenge der Verteilungsfunktionen nehmen wir dabei zur Vereinfachung der Notation implizit an, wann immer die Verteilung in Elementen aus \(\G\) oder \(\F\) betrachtet wird.
\end{Definition}

Da die Einbettungen von Knoten verschiedener Auflösungen in deren Übergangsbereichen nicht disjunkt sind, benötigen wir im Weiteren gitterspezifische Bezeichnungen für die Verteilungen und deren Momente:

\begin{Definition}[Gitterspezifische Verteilungsfunktionen und Momente]
Sei \(\# \in \{f, g\}\) wieder Symbol des feinen oder groben Gitters. Wir bezeichnen dann mit \(f_{\#,i}\) die \(i\)-te Verteilungsfunktion des entspechenden Gitters. Analog formulieren sich mit Definition~\ref{def:Momente} die Momente für \(x \in \G\) respektive \(x \in \F\):
\begin{align*}
\rho_\#(x) &:= \sum_{j=0}^{q-1} f_{\#,j}(x) \\
u_\#(x) &:= \frac{1}{\rho_\#(x)} \sum_{j=0}^{q-1} \xi_j f_{\#,j}(x)
\end{align*}
\end{Definition}

Zur Entwicklung der Gitterkopplung fordern wir, dass sich \(\G\) und \(\F\) um eine grobe Gitterweite \(\delta x_g\) überlappen -- vgl. Abbildungen \ref{fig:MultiDomainOverlap} und \ref{fig:OverlapZone}. Die Seitenlängen der konvexen Hüllen \(D_g\) und \(D_f\) sind ganzzahlige Vielfache von \(\delta x_g\) und \(\delta x_f\). Wir benötigen diese, in Form der Gitter diskretisierten, Mengen, um die Gitterknoten der Übergangsbereiche näher zu klassifizieren:

\begin{figure}[h]
\centering
\input{img/overlap_zone.tikz}
\caption{Skizze des Übergangsbereich \cite[vgl.~Abb.~4]{lagrava12}}
\label{fig:OverlapZone}
\end{figure}

\begin{Definition}[Gitterknoten der Übergangsbereiche]
\label{def:OverlapGridNodes}
\begin{align*}
\U_g &:= \G \cap \F && \text{Grobe Knoten im Übergangsbereich} \\
\U_f &:= \F \cap (D_f \cap D_g) && \text{Feine Knoten im Übergangsbereich} \\
\U_{g \to f} &:= \partial D_f \cap \U_f && \text{Knoten mit Übertragung von grob nach fein} \\
\U_{f \to g} &:= \partial D_g \cap \U_g && \text{Knoten mit Übertragung von fein nach grob}
\end{align*}
\end{Definition}

Die Übertragungsrichtungen in \(\U_{g \to f}\) und \(\U_{f \to g}\) ergeben sich aus den jeweils fehlenden Verteilungsfunktionen an den Rändern der Gitter. So fehlen beispielsweise zur Kollision der groben Gitterknoten in \(\U_{f \to g}\) Verteilungsfunktionen in Richtung des feinen Gitters, während die feinen Zellen in dieser Menge noch vollständig definiert sind, da sie im Inneren des feinen Gitters liegen.

Mit diesem Argument lässt sich auch die Notwendigkeit eines Übergangsbereiches \(\U_g \cup \U_f\) der Mindestbreite \(\delta x_g\) erklären: Gäbe es diesen nicht, so fehlten an der Grenze zwischen grobem und feinem Gitter Verteilungsfunktionen in beide Richtungen zugleich.

Die Aufgabe der Skalierungs-, Restriktions- und Interpolationsschritte des Verfahrens besteht also \emph{nur} darin, diese jeweils fehlenden Verteilungsfunktionen möglichst gut zu rekonstruieren.

\newpage
\subsection{Komponenten der Gitterkopplung}
\subsubsection{Skalierung}
\label{kap:Skalierung}

Während die Skalierung räumlicher Größen durch die Festlegung des Verfahrens auf Übergänge im Verhältnis \(1:2\) definiert ist, eröffnen sich für die zeitliche Skalierung zwei Möglichkeiten: Konvektive oder diffusive Skalierung. Unterschied der beiden Ansätze ist dabei das jeweilige Verhältnis zwischen räumlicher und zeitlicher Auflösung.

\begin{Definition}[Konvektive Skalierung]
Sei \(\delta t > 0\) die zeitliche und \(\delta x > 0\) die räumliche Diskretisierung. Dann gilt bei konvektiver Skalierung das Verhältnis:
\[ \delta t \sim \delta x \]
Es besteht also eine lineare Proportionalität.
\end{Definition}

\begin{Definition}[Diffusive Skalierung]
Sei \(\delta t > 0\) die zeitliche und \(\delta x > 0\) die räumliche Diskretisierung. Dann gilt bei diffusiver Skalierung das Verhältnis:
\[ \delta t \sim \delta x^2 \]
Es besteht hier also eine quadratische Proportionalität. Im Vergleich zu einer konvektiven Skalierung ist die zeitliche Auflösung somit um eine Ordnung feiner.
\end{Definition}

Es ist klar zu erkennen, dass diffusive Skalierung einen deutlich größeren numerischen Aufwand gegenüber der konvektiven Skalierung nach sich zieht. Vorteil der bei diffusiver Skalierung erhöhten Schrittanzahl pro Zeiteinheit sind kleinere Fehler.

Für die Autoren des hier erörterten Gitterverfeinerungsverfahrens überwog jedoch das Argument der numerischen Effizienz, weshalb auch wir hier nun die konvektive Skalierung betrachten wollen. Die Austauschbarkeit des Skalierungsverfahrens sollte jedoch bei der Implementierung in OpenLB beachtet werden, da dieser Aspekt eine weitere prinzipiell flexible Komponente des Verfahrens darstellt.

\bigskip

Aus der Wahl von konvektiver Skalierung ergibt sich zunächst:
\[\frac{\delta t_g}{\delta x_g} = \frac{\delta t_f}{\delta x_f} \land \delta x_f = \frac{\delta x_g}{2} \implies \delta t_f = \frac{\delta t_g}{2} \numberthis\label{eq:gridTime}\]
Auf dem feinen Gitter müssen also doppelt so viele Zeitschritte wie auf dem groben Gitter durchgeführt werden. Geschwindigkeit, Druck und Dichte sind stetig im Gitterübergang. Dies gilt nicht für die kinetische Viskosität \(\nu = c_s^2 \tau\), was wir bei der Herleitung der feinen Relaxionszeit \(\tau_f\) aus \(\tau_g\) beachten müssen.

\begin{Definition}[Physikalische Reynolds-Zahl]
\label{def:PhysicalReynoldsNumber}
Seien \(U\) die charakteristische Geschwindigkeit, \(L\) die charakteristische Länge und \(\nu\) die kinetische Viskosität in physikalischen Einheiten. Dann ist die Reynolds-Zahl definiert als: \[\text{Re} := \frac{U L}{\nu}\]
\end{Definition}

\begin{Definition}[Lattice Reynolds-Zahl]
\label{def:LatticeReynoldsNumber}
Sei \(\# \in \{f, g\}\) Symbol des feinen oder groben Gitters.
Seien \(U_\# := \delta t_\# / \delta x_\# \cdot U\) die charakteristische Geschwindigkeit, \(L_\# := 1 / \delta x_\# \cdot L\) die charakteristische Länge und \(\nu_\# := c_s^2 \tau_\#\) die kinetische Viskosität in Lattice-Einheiten. Dann ist die \emph{Lattice} Reynolds-Zahl des feinen bzw. groben Gitters definiert als: \[ \text{Re}_\# := \frac{U_\# L_\#}{\nu_\#} = \frac{\delta t_\# U L}{(\delta x_\#)^2 \nu_\#} \]
\end{Definition}

Wir erzwingen nun mit \(\text{Re}_g = \text{Re}_f\) die Unabhängigkeit von Reynolds-Zahl und Gitterauflösung. Diese Gleichsetzung ist sinnvoll, da die Reynolds-Zahl gerade die Vergleichbarkeit von Strömungen verschiedener Modellgrößen ermöglicht. Durch Einsetzen von Definition~\ref{def:LatticeReynoldsNumber} erhalten wir eine Verknüpfung der Relaxionszeiten \(\tau_f\) und \(\tau_g\):
\begin{align*}
\text{Re}_g = \text{Re}_f &\iff \frac{\delta t_g U L}{(\delta x_g)^2 \nu_g} = \frac{\delta t_f U L}{(\delta x_f)^2 \nu_f} \\
&\iff \frac{\delta t_g}{(\delta x_g)^2 \nu_g} = \frac{\delta t_f}{(\delta x_f)^2 \nu_f} \\
&\iff \frac{\delta t_g}{(\delta x_g)^2 c_s^2 \tau_g} = \frac{\delta t_f}{(\delta x_f)^2 c_s^2 \tau_f} \\
&\iff \tau_f = \frac{\delta t_f (\delta x_g)^2}{(\delta x_f)^2 \delta t_g} \tau_g \\
&\iff \tau_f = 2 \tau_g \numberthis\label{eq:gridTau}\\
\end{align*}

Für die zur expliziten Lösung der diskreten LBM BGK Gleichung in Definition~\ref{def:LBGKeq} verschobenen Relaxionszeiten ergibt sich somit:
\[\overline{\tau_f} = 2 \overline{\tau_g} - \frac{1}{2} \numberthis\label{eq:gridTauShift}\]

Die Equilibriumsverteilung \(f_i^\text{eq}\) ergibt sich nach Definition~\ref{def:fieq} aus Geschwindigkeit \(u\) und Dichte \(\rho\). Sie sind also explizit unabhängig der Gitterauflösung und, wie erwähnt, stetig im Gitterübergang. Diese Aussage gilt nicht für die Nicht-Equilibriumsverteilung \(f_i^\text{neq}\), da diese nach (\ref{eq:approxFneq}) von dem Geschwindigkeitsgradienten \(\nabla u\) abhängt.

Bezeichnen nun \(f_{f,i}^\text{neq}\) und \(f_{g,i}^\text{neq}\) die gitterspezifischen Nicht-Equilibriumanteile und \(\mathrm{S}_f\) sowie \(\mathrm{S}_g\) die entsprechenden Verzerrungstensoren. Zur Skalierung von \(f_{f,i}^\text{neq}\) suchen wir ein \(\alpha \in \R\) s.d. gilt: \[f_{f,i}^\text{neq} = \alpha f_{g,i}^\text{neq} \numberthis\label{eq:scaleFneqReq}\]

Mit Hilfe von \ref{eq:approxFneq} lässt sich diese Relation nun nach \(\alpha\) auflösen:
\begin{align*}
f_{f,i}^\text{neq} = \alpha f_{g,i}^\text{neq} &\iff -\frac{w_i \rho \overline{\tau_f}}{c_s^2} \mathrm{Q}_i : \mathrm{S}_f = -\alpha \left( \frac{w_i \rho \overline{\tau_g}}{c_s^2} \mathrm{Q}_i : \mathrm{S}_g \right) \\
&\iff \overline{\tau_f} \mathrm{Q}_i : \mathrm{S}_f = \alpha \overline{\tau_g} \mathrm{Q}_i : \mathrm{S}_g \\
&\iff \frac{\overline{\tau_f}}{\delta t_f} \mathrm{Q}_i : \mathrm{S} = \alpha \frac{\overline{\tau_g}}{\delta t_g} \mathrm{Q}_i : \mathrm{S} \\
&\iff \alpha = \frac{\delta t_g}{\delta t_f} \frac{\overline{\tau_f}}{\overline{\tau_g}}\\
\end{align*}

Auflösen dieses \(\alpha\) in (\ref{eq:scaleFneqReq}) und Einsetzen der Relationen (\ref{eq:gridTime}) sowie (\ref{eq:gridTauShift}) ergibt dann:
\begin{align*}
f_{f,i}^\text{neq} &= \frac{\delta t_g}{\delta t_f} \frac{\overline{\tau_f}}{\overline{\tau_g}} f_{g,i}^\text{neq} \\
&= 2 \frac{2\overline{\tau_g} - \frac{1}{2}}{\overline{\tau_g}} f_{g,i}^\text{neq} \\
&= \left( 4 - \frac{1}{\overline{\tau_g}} \right) f_{g,i}^\text{neq} \numberthis\label{eq:scaleFneq}
\end{align*}

Insgesamt haben wir hiermit die Skalierung der Diskretisierungen in Raum und Zeit, der Relaxionszeit sowie der Nicht-Equilibriumsverteilung zwischen den Gittern \(\F\) und \(\G\) hergeleitet.

\bigskip

Seien \(x_{f \to g} \in \U_{f \to g}\) und \(x_{g \to f} \in \U_{g \to f}\) die Knoten aus dem Übergangsbereich mit Übertragung von fein nach grob bzw. von grob nach fein. Dann gelten bei Erinnerung an die implizite Knotenabbildung~\ref{def:BijImplGitter}:
\begin{align}
f_{g,i}(x_{f \to g}) &= f_i^\text{eq}(\rho(x_{f \to g}), u(x_{f \to g})) + \left(4-\frac{1}{\overline{\tau_g}}\right)^{-1} f_{f,i}^\text{neq}(x_{f \to g}) \label{eq:basicF2G} \\
f_{f,i}(x_{g \to f}) &= f_i^\text{eq}(\rho(x_{g \to f}), u(x_{g \to f})) + \left(4-\frac{1}{\overline{\tau_g}}\right) f_{g,i}^\text{neq}(x_{g \to f}) \label{eq:basicG2F}
\end{align}

Die zusammengesetzten Verteilungsfunktionen von Übergangsknoten des einen Gitters lassen sich also durch Skalierung des Nicht-Equilibriumanteils der Verteilungsfunktionen des jeweils anderen Gitters rekonstruieren. Leider reicht dies noch nicht zur vollständigen Beschreibung eines Gitterverfeinerungsverfahrens, da nicht für alle feinen Gitterknoten im Übergangsbereich passende grobe Gitterknoten existieren -- vgl. dazu Abbildung~\ref{fig:OverlapZone}. Auch der Übergang von fein nach grob gestaltet sich trotz passenden feinen Knoten potenziell komplizierter, als eine bloße Skalierung, wie wir im nächsten Kapitel zeigen werden.

\newpage
\subsubsection{Restriktion}

Kraft seiner höheren Auflösung enthält das feine Gitter mehr Informationen als das umgebende grobe Gitter. Der Übergang von fein nach grob stellt also eine Restriktion der Verteilungsfunktionen dar.

Konkret suchen wir nach einer sinnvollen Definition der in \(x_{f \to g} \in \U_{f \to g}\) fehlenden Verteilungsfunktionen \(f_{g,i}\). Eine Solche ergibt sich aus der skalierten Dekomposition \ref{eq:basicF2G} durch Ersetzen der einfachen Nicht-Equilibriumsverteilung \(f_{f,i}^\text{neq}(x_{f \to g})\) mit einer restringierten Variante ebendieser.
\[f_{g,i}(x_{f \to g}) = f_i^\text{eq}(\rho(x_{f \to g}), u(x_{f \to g})) + \alpha^{-1} \ \resarg{i}{x_{f \to g}} \numberthis\label{eq:restrictedF2G}\]
Wir bemerken an dieser Stelle, dass nur die Nicht-Equilibriumsverteilung durch eine Restriktionsoperation eingeschränkt wird, während der Equilibriumanteil unangetastet bleibt. Dies ist damit zu begründen, dass Dichte und Geschwindigkeit bei der von uns verwendeten konvektiven Skalierung im Gitterübergang stetig bleiben.

Die skalierte Dekomposition \ref{eq:basicF2G} lässt sich in der Schreibweise von \ref{eq:restrictedF2G} formulieren, wenn die Identität als Restriktionsoperation eingesetzt wird: \[\resarg{i}{x_{f \to g}} := f_{f,i}^\text{neq}(x_{f \to g})\]

Die für unser Verfahren \cite[Kap.~3.3]{lagrava12} beschriebene Restriktion ist der Mittelwert aller umliegenden gerichteten Nicht-Equilibriumanteilen:
\[\resarg{i}{x_{f \to g}} := \frac{1}{q} \sum_{j=0}^{q-1} f_{f,j}^\text{neq}(x_{f \to g} + \delta x_f \xi_j) \numberthis\label{eq:neqAvgRestrictionF2G}\]

\newpage
\subsubsection{Interpolation}

Zunächst ergänzen wir die Gitterteilmengen aus Definition~\ref{def:OverlapGridNodes} um eine Unterscheidung zwischen alleinstehenden feinen Knoten und solchen für die ein, der Übertragung von grob nach fein dienlicher, grober Knoten existiert.
\begin{Definition}[Gitterknoten mit Übertragung von grob nach fein]
\begin{align*}
\U_{g \to f}^g &:= \U_{g \to f} \cap \U_g && \text{Doppelte Knoten mit Übertragung von grob nach fein} \\
\U_{g \to f}^f &:= \U_{g \to f} \setminus \U_g && \text{Alleinstehende feine Knoten mit Übertragung von grob nach fein}
\end{align*}
\end{Definition}

Für \(x_{g \to f}^g \in \U_{g \to f}^g\) gilt insbesondere \(x_{g \to f}^g \in \U_g \cap \, \U_f\). Es existieren in diesen Gitterpunkten also vollständig definierte grobe Verteilungsfunktionen, die wir zur Bestimmung der Momente \(\rho\) und \(u\) in (\ref{eq:basicG2F}) heranziehen können. Entsprechend besitzen wir zur Interpolation des gesuchten Wertes geschickterweise eine Stützstelle an eben dessen Position. Die \emph{Interpolation} von \(x_{g \to f}^g\) beschränkt sich folglich auf die Skalierung~(\ref{eq:basicG2F}):
\[f_{f,i}(x_{g \to f}^g) = f_i^\text{eq}(\rho_g(x_{g \to f}^g), u_g(x_{g \to f}^g)) + \alpha f_{g,i}^\text{neq}(x_{g \to f}^g) \numberthis\label{eq:expandedDirectG2F}\]

Für \(x_{g \to f}^f \in \U_{g \to f}^f\) gilt insbesondere \(x_{g \to f}^f \notin \U_g\). Es existieren in diesen Gitterpunkten also im Gegensatz zur Situation \ref{eq:expandedDirectG2F} keine groben Verteilungsfunktionen. Die fehlenden Werte zur Bestimmung der Momente sowie des Nicht-Equilibriumanteils in (\ref{eq:basicG2F}) müssen hier also aus den umliegenden groben Verteilungsfunktionen interpoliert werden:
\[f_{f,i}(x_{g \to f}^f) = f_i^\text{eq}(\ipolarg{\rho_g}{x_{g \to f}^f}, \ipolarg{u_g}{x_{g \to f}^f}) + \alpha \ipolarg{f_{g,i}^\text{neq}}{x_{g \to g}^f} \numberthis\label{eq:expandedInterpolG2F}\]
Die unbekannten Werte der Moment- und Nicht-Equilibriumfunktionen werden in diesem Ausdruck durch eine Interpolationsoperation \(\ipol\) genähert. Neben der Art der Restriktion \(\res\) stellt die Wahl des Interpolationsverfahrens einen weiteren zentralen und flexiblen Bestandteil des, auf diesen Seiten nachvollzogenen, Gitterverfeinerungsverfahren dar.

\bigskip

Stützstellen für die Interpolation seien hier die parallel zum Gitterübergang liegenden groben Nachbarknoten des gesuchten Punktes \(x_{g \to f}^f\). Wir adressieren diese, parallel zu einem Einheitsvektor \(v\) positionierten, Stützknoten mit:
\[\mathcal{N}_{x_{g \to f}^f} := \left\{x \in \G \middle| x = x_{g \to f}^f + j \, \delta x_f \, v,\ j \in \Z \right\} \subseteq \U_{g \to f}^g\]
Bekannte Stützwerte von \(\ipolarg{\star}{x}\) befinden sich also relativ zum gesuchten Knoten \(x_{g \to f}^f\) in, um ungerade Vielfache der feinen Schrittweite \(\delta x_f\) skalierten, Verschiebungen entlang der normierten Übergangsparallele \(v\) (vgl. Abbildungen \ref{fig:InterpolationBasis} und \ref{fig:InterpolationDetail}). Die Einschränkung der hinzugezogenen Stützen auf Knoten, welche parallel zum Übergang liegen, wurde von Lagrava et al. so gewählt~\cite[Kap.~3.6]{lagrava12}, um in 2D ein eindimensionales Interpolationsproblem zu erhalten. Prinzipiell spricht nichts gegen eine Einbeziehung umfangreicherer Teilmengen der groben Knotennachbarschaft.

\begin{figure}[h]
\centering
\input{img/interpolation_basis.tikz}
\caption{Stützstellen der Interpolation im Übergangsbereich}
\label{fig:InterpolationBasis}
\end{figure}

Um die kommenden Ausführungen auf das Wesentliche -- namentlich das Verfahren unabhängig der konkret zu interpolierenden Funktion -- zu konzentrieren, sei definiert:
\[\sipolarg{h} := \ipolarg{\star}{x_{g \to f}^f + h \, \delta x_f \, v} \text{ für Zielfunktion } \star \in \{\rho_g, u_g, f_{g,i}^\text{neq}\}\]
In dieser Formulierung suchen wir also eine möglichst gute Interpolation des Wertes in \(\sipolarg{0}\) anhand der Stützstellen \(\sipolarg{h}\) für kleine \(h \in \Z \setminus 2\Z\). Ein naheliegender Ansatz hierfür ist das arithmetische Mittel der beiden engsten Nachbarn:
\[\ipolarg{\star}{x_{g \to f}^f} = \sipolarg{0} = \frac{\sipolarg{-1} + \sipolarg{1}}{2} \numberthis\label{eq:ipol2ord}\]
Vorteil eines solch einfachen Verfahrens wäre, dass die benötigten groben Nachbarn auch an den Ecken des Übergangsbereiches existieren (vgl. Abbildung~\ref{fig:InterpolationEdgeCase}) und daher keine Sonderbehandlung erforderlich wird.
\begin{figure}[h]
\centering
\input{img/interpolation_detail.tikz}
\caption{Stützstellen der Interpolation im Detail}
\label{fig:InterpolationDetail}
\end{figure}

Bessere Näherungen können unter Einsatz weiterer Stützknoten erzielt werden. Wir berechnen dazu mit dem Schema der dividierten Differenzen die Faktoren der Newtonschen Interpolationsformel \cite[IV.3~(3.10)]{amann_escher} auf den in Abbildung~\ref{fig:InterpolationDetail} dargestellten Stützpunkten \((-3,\sipolarg{-3})\), \((-1,\sipolarg{-1})\), \((1,\sipolarg{1})\) und \((3,\sipolarg{3})\):
\begin{align*}
\sipolarg{x} :&= \sipolarg{-3} \\
&+ \frac{1}{2}(\sipolarg{-1}-\sipolarg{-3})(x+3)\\
&+ \frac{1}{8}(\sipolarg{1}-\sipolarg{-1}+\sipolarg{3})(x+3)(x+1) \\
&+ \frac{1}{48}(\sipolarg{3}-3\sipolarg{1}+3\sipolarg{-1}-\sipolarg{-3})(x+3)(x+1)(x-1)
\end{align*}
Ausgewertet in \(0\) erhalten wir folgenden Ausdruck als Näherung von \(\ipolarg{\star}{x_{g \to f}^f}\):
\[\sipolarg{0} = \frac{9}{16}(\sipolarg{-1} + \sipolarg{1}) - \frac{1}{16}(\sipolarg{-3} + \sipolarg{3}) \numberthis\label{eq:ipol4ord}\]
Hierbei handelt es sich um ein Verfahren vierter Ordnung, wie sich durch Einsetzen der Taylor-Entwicklung von \(\sipol\) um \(0\) in die Auswertung~(\ref{eq:ipol4ord}) zeigen lässt:
\begin{align*}
&\sipolarg{h} = \sipolarg{0} + \sipolderivarg{1}{0}h + \frac{1}{2}\sipolderivarg{2}{0}h^2 + \frac{1}{6}\sipolderivarg{3}{0}h^3 + \mathcal{O}(h^4) \numberthis\label{eq:sipolTaylorOrder4} \\
\implies &\frac{9}{16}(\sipolarg{-1} + \sipolarg{1}) - \frac{1}{16}(\sipolarg{-3} + \sipolarg{3}) \stackrel{(\ref{eq:sipolTaylorOrder4})}{=} \sipolarg{0} + \mathcal{O}(h^4)
\end{align*}

In Abbildung~\ref{fig:InterpolationEdgeCase} erkennen wir zwei mögliche Randfälle des Gitterübergangs, welche zu Problemen bei Nutzung des Interpolationsverfahren vierter Ordnung führen können.
\begin{figure}[h]
\centering
\input{img/interpolation_edge_case.tikz}
\caption{Interpolation in Ecken und Enden des feinen Gitters}
\label{fig:InterpolationEdgeCase}
\end{figure}

\noindent
So kann es an den Außengrenzen der Simulationsdomäne dazu kommen, dass nur drei der vier benötigten groben Nachbarknoten zur Verfügung stehen, wie in der grün markierten Situation dargestellt.
Je nach Implementierung der Kommunikation zwischen den Gittern, kann die gleiche Einschränkung auch in Ecken des feinen Gitters -- hier blau markiert -- dazu kommen, dass eine Interpolation auf Grundlage von nur drei Nachbarn benötigt wird. Entsprechend erhalten wir nach Berechnung der dividierten Differenzen ein Interpolationspolynom auf drei Stützstellen:
\begin{align*}
\sipolarg{x} :&= \sipolarg{-1} \\
&+ \frac{1}{2}(\sipolarg{1}-\sipolarg{-1})(x+1)\\
&+ \frac{1}{8}(\sipolarg{3}-2\sipolarg{1}+\sipolarg{-1})(x+1)(x-1)
\end{align*}
Auch in dieser Situation ist nur der Wert in \(0\) zur Näherung von \(\ipolarg{\star}{x_{g \to f}^f}\) von Interesse:
\[\sipolarg{0} = \frac{3}{8}\sipolarg{-1} + \frac{3}{4}\sipolarg{1} - \frac{1}{8}\sipolarg{3} \numberthis\label{eq:ipol3ord}\]
Passend zur Anzahl der Stützstellen präsentiert sich dieses Verfahren nach erneutem Einsetzen der Taylor-Entwicklung (\ref{eq:sipolTaylorOrder4}) als eine Interpolationsformel dritter Ordnung:
\[\frac{3}{8}\sipolarg{-1} + \frac{3}{4}\sipolarg{1} - \frac{1}{8}\sipolarg{3} \stackrel{(\ref{eq:sipolTaylorOrder4})}{=} \sipolarg{0} + \mathcal{O}(h^3)\]

Trotz Behandlung dieses Sondernfalls werden die höheren Approximationsordnungen gegenüber (\ref{eq:ipol2ord}) weiterhin und unumgänglich durch zusätzliche Stützstellen erkauft, welche bei parallelisierten LBM-Umsetzungen kommuniziert werden müssen. Es gilt also, zwischen Güte der Interpolation und Anzahl sowie Position der Stützstellen abzuwiegen.

\newpage
\subsection{Algorithmus}

In der zurückliegenden Sektion haben wir, aufbauend auf der Skalierung von Verteilungsfunktionen zwischen Gitterauflösungen, die Restriktion von fein nach grob sowie die Interpolation von grob zu fein nachvollzogen. Die somit erfassten wesentlichen Grundlagen des Verfeinerungsverfahrens gilt es nun in einem \emph{Kopplungsalgorithmus}~\cite[Kap.~3.5]{lagrava12} zusammenzuführen.

\bigskip

Entsprechend (\ref{eq:gridTime}) müssen für jeden groben Zeitschritt \(\delta t_g\) zwei feine Zeitschritte \(\delta t_f\) durchgeführt werden. Die alternierenden Kollisions- und Strömungsschritte der beiden Gitter sind also strikt gekoppelt und werden als eine Schleifeneinheit betrachtet. Als Schleifeninvariante definieren wir dabei die vollständige Bekanntheit aller Verteilungsfunktionen aller Knoten in beiden Gittern.

\begin{figure}[h]
\input{img/algorithm_birds_eye.tikz}
\caption{Verfeinerungsalgorithmus mit Invariante aus der Vogelperspektive}
\label{fig:AlgorithmBirdsEye}
\end{figure}
\noindent
Aufbauend auf dieser Invariante ergibt sich die, in Abbildung~\ref{fig:AlgorithmBirdsEye} dargelegte, Reihenfolge der erforderlichen Schritte direkt aus den, für die einzelnen Komponenten der Gitterkopplung benötigen, Informationen. So sind zu Beginn alle Verteilungsfunktionen vollständig bekannt, was die Ausführung eines üblichen Kollisions- und Strömungsschritts (vgl. Kapitel~\ref{kap:LBMimpl}) in beiden Gittern ohne weitere Zuarbeit erlaubt. Nach diesen beiden Schritten fehlen Verteilungsfunktionen \(f_{g,i}(x_{f \to g})\) zur Wiederherstellung der Invariante des groben Gitters. Auch der benötigte zweite Simulationsschritt, um \(\F\) auf Zeitpunkt \(t+\delta t_g=t+2\delta t_f\) zu bringen, scheitert zunächst an der Unbestimmtheit von Verteilungsfunktionen \(f_{f,i}(x_{g \to f})\).

\begin{figure}[h]
\centering
\input{img/algo_completion_overview.tikz}
\caption{Übersicht der zu vervollständigenden Knoten}
\end{figure}

\begin{description}[style=unboxed,leftmargin=0cm]
\item[Vervollständigung von \(\F\) zu Zeitpunkt \(t+\delta t_f\):] Zur Vervollständigung des feinen Gitters nach dem ersten Zeitschritt müssen die fehlenden Verteilungen aus dem groben Gitter rekonstruiert werden. Um die dazu erarbeiteten Kopplungen (\ref{eq:expandedDirectG2F}) und (\ref{eq:expandedInterpolG2F}) anzuwenden, fehlen jedoch Werte der groben Stützstellen \(\smash{f_{g,i}(x_{g \to f}^g})\) zu Zeitpunkt \(t+\delta t_f\). Diese sind zwar in den gesuchten Punkten, dank Trennung der Kopplungsrichtungen durch den Übergangsbereich, nach jedem Simulationsschritt direkt vollständig vorhanden -- jedoch nur zu Zeit \(t\) und \(t+\delta t_g\). Hier findet sich eine Anwendung des Interpolationsverfahrens zweiter Ordnung (\ref{eq:ipol2ord}) zur linearen Zeitinterpolation der benötigten Werte von \(\rho_g, u_g\) und \(f_{g,i}^\text{neq}\):
\[\star(x,t+\delta t_f) \approx \frac{\star(x,t+\delta t_g) + \star(x,t)}{2} \text{ für } \star \in \{\rho_g,u_g,f_{g,i}^\text{neq}\}, x \in \G\]
Aufbauend darauf steht der Anwendung der Kopplungsformeln (\ref{eq:expandedDirectG2F}) und (\ref{eq:expandedInterpolG2F}) zur Vervollständigung von \(\F\) nichts mehr im Wege:
\begin{align*}
f_{f,i}(x_{g \to f}^g,t+\delta t_f) &= f_i^\text{eq}(\rho_g(x_{g \to f}^g,t+\delta t_f), u_g(x_{g \to f}^g,t+\delta t_f)) + \alpha f_{g,i}^\text{neq}(x_{g \to f}^g,t+\delta t_f) \\
f_{f,i}(x_{g \to f}^f,t+\delta t_f) &= f_i^\text{eq}(\ipolarg{\rho_g}{x_{g \to f}^f}, \ipolarg{u_g}{x_{g \to f}^f}) + \alpha \ipolarg{f_{g,i}^\text{neq}}{x_{g \to g}^f}
\end{align*}
Der Interpolationsoperator vierter Ordnung (\ref{eq:ipol4ord}) löst sich dabei für die Zielfunktionen \(\star \in \{\rho_g,u_g,f_{g,i}^\text{neq}\}\) und die Übergangsparallele \(v\) wie folgt auf:
\begin{align*}
\ipolarg{\star}{x_{g \to f}^f} = &\frac{9}{16}(\star(x_{g \to f}^f-\delta x_f v, t+\delta t_f) + \star(x_{g \to f}^f+\delta x_f v, t+\delta t_f))\\
+ &\frac{1}{16}(\star(x_{g \to f}^f-3\delta x_f v, t+\delta t_f) + \star(x_{g \to f}^f+3\delta x_f v, t+\delta t_f))
\end{align*}

\item[Vervollständigung von \(\F\) zu Zeitpunkt \(t+\delta t_g\):] Dieser zweite Rekonstruktionsschritt auf dem feinen Gitter gestaltet sich einfacher, da die benötigten groben Verteilungen in \(\U_{g \to f}\) zur Zeitpunkt \(t+\delta t_g\) bereits durch den initialen Simulationschritt auf dem groben Gitter bekannt sind. Entsprechend können die Kopplungsformeln (\ref{eq:expandedDirectG2F}) und (\ref{eq:expandedInterpolG2F}) direkt zur Vervollständigung von \(\F\) angewandt werden.

\item[Vervollständigung von \(\G\) zu Zeitpunkt \(t+\delta t_g\):] Nach zweimaliger Vervollständigung des feinen Gitters verbleibt zur Wiederherstellung der Schleifeninvariante der Abschluss des eingehenden Kollisions- und Strömungsschritts auf dem groben Gitter durch Restriktion der aus Richtung des feinen Gitters eingehenden Verteilungsfunktionen. Hierzu erlaubt die, durch die zuvorkommenden Schritte garantierte, Vollständigkeit des feinen Gitters zu Zeitpunkt \(t+\delta t_g\), die direkte Anwendung der Kopplungsformel (\ref{eq:restrictedF2G}) mit Restriktionsoperator (\ref{eq:neqAvgRestrictionF2G}) auf die Knoten in \(\U_{f \to g}\).
\begin{align*}
f_{g,i}(x_{f \to g},t+\delta t_g) &= f_i^\text{eq}(\rho_f(x_{f \to g},t+\delta t_g), u_f(x_{f \to g},t+\delta t_g))\\
&+ \frac{1}{\alpha} \frac{1}{q} \sum_{j=0}^{q-1} f_{f,j}^\text{neq}(x_{f \to g} + \delta x_f \xi_j,t+\delta t_g)
\end{align*}
\end{description}
Zu erwähnen bleibt, dass wir aus Konsistenzgründen alle Kopplungsformeln immer auf alle -- und nicht nur die fehlenden -- diskreten Richtungen \(i \in [q-1]\) einer betrachteten Zelle \(x\) anwenden.

\bigskip
Nach Durchführung der drei Vervollständigungsschritte haben wir die Invariante für \(t+\delta t_g\) wieder hergestellt und können von vorne beginnen. Wir haben damit an dieser Stelle das Verfeinerungsverfahren von Lagrava et al. vollständig nachvollzogen und können mit der Implementierung in OpenLB fortfahren.

% ToDo: Revisit Multi-Grid: Kopplung über der ganzen Domäne
	% Sollte so sein - Kopplung nur am Rand macht Inneres irrelevant, d.h. Domain==Grid
% ToDo: Einschränkungen der Gitterpositionierung (keine hängenden feinen Knoten) ausarbeiten
% ToDo: Randfälle der Restriktion ausarbeiten, analog zu Interpolation (fehlt im Paper)
% ToDo: Experimentelle Begründung, warum Kopplungsformel immer auf alle Richtungen angewandt wird

\newpage
\section{Implementierung in OpenLB}

OpenLB~\cite{olb12}

\subsection{Auswahl der Verfeinerungsmethode}\label{sec:olbRefinementChoice}

\subsection{Struktur des Gitterverfeinerungsframework}

\subsection{Umsetzung des Verfahrens von Lagrava et al.}

\newpage
\section{Evaluierung}

\subsection{Wahl der Beispiele}

\subsection{Rohrströmung}