aboutsummaryrefslogtreecommitdiff
path: root/content.tex
blob: a1286b6295ebaf387884e0e0c92eca7d1f533d98 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
\section{Einführung}

\subsection{Wieso Strömungen simulieren?}

Das hoch technisierte Lebensumfeld des modernen Menschen ist ohne ein detailliertes Verständnis des Verhaltens der, durch ihn in wachsendem Maße kontrollierten, Natur undenkbar. Eine wichtige Komponente dieses Naturverständnisses ist das Wissen um das Verhalten von Flüssigkeiten und Gasen, die sich als Strömungen bewegen. Ohne dieses Verständnis führe kein Automobil, flöge kein Flugzeug und drehten sich weder ein Windrad um seine Achse noch ein Satellit um unseren Planeten.

Experimente zur Bestimmung des Verhaltens von Strömungen -- z. B. in Wind- und Wasserkanälen -- sind möglich, liefern naturnaheste Ergebnisse und stellen eine zentrale Komponente der Entwicklung eben genannter Errungenschaften dar. Leider sind reale Experimente im Allgemeinen nicht nur sehr aufwendig in Aufbau und Durchführung, sondern stoßen auch insbesondere bei der Betrachtung mikroskopischer Probleme -- etwa im Bereich der Medizin, deren menschliches Subjekt zu einem großen Teil ebenfalls eine \emph{Strömungsmaschine} bildet -- an Grenzen von Messmethoden und Ethik.

\begin{figure}[h]
\centering
\includegraphics[scale=0.35]{img/static/nose.png}
\caption{Strömung in der komplexen Geometrie der menschlichen Nase \cite{olbNose}}
\end{figure}

So trifft es sich, dass der Wunsch nach theoretischer Lösung komplexer und nur schwer analytisch zugänglicher Strömungsprobleme mit der Entwicklung von immer leistungsfähigeren Rechenmaschinen nicht nur einherging, sondern auch eine der Triebfedern in deren initialen Entwicklung war. Moderne numerische Verfahren zur Simulation von Strömungen versprechen eine zunehmende Reduzierung benötigter realer Experimente und sind heute gängiges Werkzeug in Forschung und Maschinenbau.

\subsection{Weshalb mit Lattice Boltzmann Methoden?}

Während Finite Elemente Methoden die wohl verbreitetsten Zugänge zur Lösung der strömungsbeschreibenden partiellen Differentialgleichungen und somit der numerischen Strömungsmechanik bilden, erfreut sich auch die Herangehensweise der Lattice Boltzmann Methoden in den letzten Jahrzehnten wachsender Nutzbarkeit und Verbreitung. Im Gegensatz zu anderen Lösungsmethoden werden hier die Navier-Stokes Gleichungen nicht direkt numerisch gelöst. Lösungen ergeben sich vielmehr aus der Simulation des Fluidverhaltens auf \emph{mesoskopischer} Ebene -- d. h. aus der Betrachtung nicht aus Sicht der Kollision einzelner Fluidpartikel und nicht aus Sicht der analytischen Strömungsbeschreibung, sondern aus Sicht der Wahrscheinlichkeit, dass sich das Fluid zu bestimmter Zeit an einem bestimmten Ort mit bestimmter Geschwindigkeit bewegt.

\bigskip
Ein Vorteil dieses, auf den Arbeiten von Ludwig Eduard Boltzmann im Bereich der statistischen Physik aufbauenden, Ansatzes, ist seine Eignung für komplexe Geometrien mit verschiedensten Randbedingungen \cite{Aidun10}. Weiterhin gewinnt in den letzten Jahren auch die sehr gute Parallelisierbarkeit von Lattice Boltzmann Methoden in Hinblick auf einen technischen Fortschritt an Anziehungskraft, nach welchem die Leistungsfähigkeit von Großrechnern eher aus deren Parallelität als aus individueller Prozessorleistung erwächst.

\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 -- beispielsweise in komplexen Geometrien und an Rändern -- 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, diese 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. Tabelle~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~\cite[Kap.~6]{AmannEscherII}:
\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 Simulationsverlauf ohne aufwendige Restrukturierung verschoben werden können, um beispielsweise komplexere Strömungsstrukturen zu \emph{verfolgen}. Nachteil ist, dass das Einsparpotentiale in Speicher- und Rechenbedarf wegen mehrfacher Abdeckung von Teilen der Simulationsdomäne durch überlagerte 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.

\subsubsection{Koinzidierende und versetzte Gitter}

Neben der Unterscheidung zwischen Multi-Grid und Multi-Domain Ansätzen kann auch die Positionierung des groben im Bezug auf das feine Gitter varieren. Gitterkonstellationen, in denen grobe Knoten soweit möglich mit feinen Knoten übereinstimmen -- wie auch in den beiden vorherigen Abbildung zu sehen -- werden als koinzidierend bezeichnet:

\begin{figure}[h]
\centering
\input{img/coinciding_grid.tikz}
\caption{Koinzidierende Gitter}
\label{fig:CoincidingGrid}
\end{figure}

\noindent
Lösen wir uns von dieser Einschränkung und positionieren das feine Gitter abseits der groben Gitterknoten, sprechen wir über zueinander versetzte -- im angelsächsischen Sprachraum als \emph{staggered} beschriebene -- Gitter:

\begin{figure}[h]
\centering
\input{img/staggered_grid.tikz}
\caption{Zueinander versetzte Gitter ohne überlappende Knoten}
\label{fig:StaggeredGrid}
\end{figure}

\newpage
\section{Verfeinerungsmethode nach Lagrava et al.}\label{kap:Lagrava}

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 dem Multi-Domain Ansatz \cite[Kap.~3.1]{Lagrava12} mit koinzidierenden Gittern. Es werden also die feiner aufgelösten Teilbereiche der Simulationsdomäne so aus dem gröber aufgelösten Gitter ausgeschlossen, dass sie 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}

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

\bigskip
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.

\newpage
\subsection{Gitterdiskretisierung}

\begin{Definition}[Grobe und feine Simulationsdomänen]
\label{def:SimDomain}
Sei \(D \subseteq \R^2\) die physikalische Simulationsdomäne unabhängig der verwendeten Gitterauflösung. Diese Teilmenge von \(\R^2\) sei dabei bereits in Hinblick auf die angestrebte kartesische Diskretisierung gewählt, d. h. als Vereinigung zusammenhängender Rechtecke.

Den durch ein grobes Gitter abzubildenden Teilbereich bezeichnen wir mit \(D_g \subset D\). Weiter ist \(D_f \subset D\) mit \(D_f \cap D_g \neq \emptyset\) der zu verfeinernde Teilbereich. O.B.d.A. nehmen wir an, dass \(D_g \cup D_f = D\) gilt, die Simulationsdomäne also in genau zwei Auflösungsstufen aufgeteilt wird.
\end{Definition}

\begin{figure}[h]
\centering
\input{img/simulation_domain.tikz}
\caption{Schematische Darstellung der Simulationsdomänen der Gitter}
\label{fig:SimDomain}
\end{figure}

Die Annahme von zwei Verfeinerungsstufen ist berechtigt, da mehrfach verfeinerte Gitter sich durch erneute -- rekursive -- Anwendung von Definition~\ref{def:SimDomain} mit \(\tilde{D} := D_f\) modellieren lassen. Auch Verfeinerung eines groben Gitters in mehreren disjunkten Bereichen kann unter dieser Annahme betrachtet werden, da für die Koppelung zwischen Gittern nur Knoten in \(D_f\) relevant sind.

\begin{Definition}[Diskretisierung der Simulationsdomänen]
\label{def:DiskretRefinedGitter}
Wir betrachten kartesische Gitter \(\G\) und \(\F\) als Diskretisierungen der Simulationsdomänen \(D_g\) bzw. \(D_f\). Diese seien so gewählt, dass sie gerade die konvexen Hüllen ihrer koinzidierenden Diskretisierungsgitter beschreiben.
\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}

Zur Entwicklung der Gitterkopplung fordern wir, dass sich \(\G\) und \(\F\) um mindestens 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\). Formal sei dabei das Innere der groben Domäne \(D_g\) ohne den Übergangsbereich der Breite \(\delta x_g\) gegeben als:
\[ D_g^\circ := \{ x \in D_g | \forall y \in \R^2 \setminus D_g : \|x-y\| > \delta x_g \} \]
Vergleiche hierzu den unschraffierten Bereich der Darstellung von \(D_g\) in Abbildung~\ref{fig:SimDomain}.
Unter der Annahme, dass \(D_g\) die feine Simulationsdomäne \(D_f\) komplett umschließt, also der komplette Rand des feinen Gitters mit dem groben gekoppelt werden soll, fordern wir dann für den Rand des feinen Gitters:
\[ \partial D_f \stackrel{!}{\subset} \partial D_g^\circ \numberthis\label{eq:FineBorderIntersectCoarseGrid} \]
Diese Einschränkung garantiert, dass die äußersten Gitterknoten des feinen Gitters sich soweit möglich mit groben Gitterknoten schneiden und nicht etwa zwischen zwei groben \emph{Gitterreihen} liegen.

\begin{figure}[h]
\centering
\input{img/invalid_overlap_area.tikz}
\caption{Mit Forderung (\ref{eq:FineBorderIntersectCoarseGrid}) ausgeschlossener Übergangsbereich \cite[vgl. Abb.~9]{Lagrava12}}
\label{fig:InvalidOverlapArea}
\end{figure}

\noindent
Wir können die Gitterknoten der Übergangsbereiche nun detailliert 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.

\bigskip
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}

\bigskip
Zusammenfassend wird die Aufgabe der im kommenden Kapitel zu erarbeitenden Skalierungs-, Restriktions- und Interpolationsschritte also \emph{nur} darin bestehen, die 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 dabei 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 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 \overline{\tau_f} \delta t_f \mathrm{Q}_i : \mathrm{S} = \alpha \overline{\tau_g} \delta t_g \mathrm{Q}_i : \mathrm{S} \\
&\iff \alpha = \frac{\delta t_f}{\delta t_g} \frac{\overline{\tau_f}}{\overline{\tau_g}}\\
\end{align*}
Einsetzen der Relationen (\ref{eq:gridTime}) und (\ref{eq:gridTauShift}) reduziert den Skalierungsfaktor auf einen nur von der groben Relaxionszeit abhängigen Ausdruck:
\begin{align*}
\alpha &= \frac{\delta t_f}{\delta t_g} \frac{\overline{\tau_f}}{\overline{\tau_g}} \\
&= \frac{1}{2} \frac{2\overline{\tau_g} - \frac{1}{2}}{\overline{\tau_g}} \\
&= 1 - \frac{1}{4\overline{\tau_g}} \numberthis\label{eq:scaleFactor}
\end{align*}
Schließlich erhalten wir so die folgende Relation der Nicht-Equilibriumsverteilungen:
\[f_{f,i}^\text{neq} = \left( 1-\frac{1}{4\overline{\tau_g}} \right) f_{g,i}^\text{neq} \numberthis\label{eq:scaleFneq}\]

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(1-\frac{1}{4\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(1-\frac{1}{4\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,i}^\text{neq}(x_{f \to g} + \delta x_f \xi_j) \numberthis\label{eq:neqAvgRestrictionF2G}\]

\newpage
\subsubsection{Interpolation}\label{kap: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)]{AmannEscherI} 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}\label{kap: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,i}^\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: Randfälle der Restriktion ausarbeiten, analog zu Interpolation (fehlt im Paper)
% ToDo: Experimentelle Begründung, warum Kopplungsformel immer auf alle Richtungen angewandt wird
% ToDo: Bemerkungen zu MPI Unterstützung

\newpage
\section{Implementierung in OpenLB}

OpenLB~\cite{olb12} ist ein umfangreiches frei verfügbares C++ Framework zur Implementierung von LBM basierenden Simulationen. Schwerpunkte sind dabei eine große Flexibilität in Hinblick auf die unterstützten Lattice Boltzmann Modelle sowie weitreichende Modularisierung zur einfachen Umsetzung neuer Anwendungen bei gleichzeitiger Eignung für hochperformante Berechnungen durch Skalierbarkeit auf parallele Großrechner.

\bigskip
Eine LBM Bibliothek mit idealer Unterstützung für Gitterverfeinerung würde es erlauben, dass zu lösende physikalische Problem zunächst ohne Rücksicht auf die konkrete Zusammensetzung des Gitters zu modellieren. Entstünde dann im Konflikt zwischen Rechenressourcen und angestrebter Ergebnisqualität ein Bedarf für Gitterverfeinerung, sollte diese a posteriori ohne weitere Anpassungen hinzugeschalten werden können. In einer solchen Umgebung wäre es dann sogar denkbar, die Entscheidung über Position und Ausmaß der Verfeinerung an ein automatisches Kriterium \cite{Lagrava15} zu übertragen.

Diese Vision ist selbstredend von großem Anspruch und weit außerhalb der Grenzen dieser Arbeit. Ein erster Schritt muss jedoch die Ergänzung von OpenLB um grundlegende Unterstützung für manuelle Gitterverfeinerung sein -- d. h. Position, Größe und Auflösung des zu verfeinernden Bereichs müssen durch den Nutzer vorgegeben werden können.

\subsection{Architekturübersicht}

OpenLB verwaltet die Diskretisierung der zu modellierenden Simulationsdomäne in einer \class{CuboidGeometry2D} Instanz. Diese Klasse teilt eine gegebene physikalische Domäne in eine beliebige Anzahl von, durch \class{Cuboid2D} Instanzen dargestellten, achsenparallelen Rechtecken mit definierter Auflösung ein. Diese Aufteilung dient der Parallelisierung, da die Gitter dieser einzelnen Quader außerhalb von kommunizierenden Randbereichen vollständig nebenläufig verarbeitet werden können. Entsprechend werden diese Teilbereiche von einer Implementierung der \class{LoadBalancer} Schnittstelle auf die zu Verfügung stehenden Prozessoren aufgeteilt.

Auf Grundlage dieses räumlichen Umrisses wird dann das eigentliche Gitter in einer \class{SuperLattice2D} Instanz abhängig des gewählten Lattice Boltzmann Modells aufgebaut. Die einzelnen Gitterzellen werden dabei von \class{Cell} Klassen abgebildet, welche anhand des spezifischen, durch einen sogenannten Deskriptor \emph{beschriebenen} und von \class{Dynamics} Instanzen durchgeführten, Kollisionsschrittes der lokalen Umsetzung des LBM Modells Sorge tragen.

Diese \class{Dynamics} beschreiben dabei zusammen mit zellübergreifenden Postprozessoren nicht nur das einfache Strömungsverhalten, sondern modellieren auch die zahlreichen Randkonditionen, welche die Abbildung komplexer Geometrien in LBM erst ermöglichen.
Zur einfacheren Zuordnung dieser zellspezifischen Klassen verwendet OpenLB sogenannte Materialzahlen, welche von einer \class{SuperGeometry2D} Instanz verwaltet werden.

Parallel zu diesen Strukturen kapselt und berechnet die \class{UnitConverter} Klasse etwaig benötigte Konstanten zur Konvertierung zwischen physikalischen Einheiten und den, diese modellierenden, Lattice-Einheiten sowie Relaxionszeiten und Fluidkonstanten wie die Reynolds-Nummer.

\bigskip
Im Allgemeinen ergibt sich aus diesen Komponenten folgende übliche Struktur von OpenLB-basierenden Anwendungen \cite[Kap.~2.1]{olb12userguide}:
\begin{enumerate}
	\item Erstellung des \class{UnitConverter} mit den beabsichtigten Gitterkonstanten
	\item Beschreibung der Simulationsdomäne durch Konstruktion einer \class{CuboidGeometry2D}
	\item Bereitstellung eines \class{LoadBalancer} zur Instantierung einer \class{SuperGeometry2D}
	\item Definition der Materialzahlen in einer \method{prepareGeometry} Methode
	\item Konstruktion der \class{SuperLattice2D} Instanz aus der \class{SuperGeometry2D}
	\item Instantierung der benötigten \class{Dynamics} und etwaigen Randkonditionen
	\item Bindung von \class{Dynamics} und Randkonditionen an die, von \class{SuperLattice2D} verwalteten, \class{Cell} Objekte anhand der Materialzahlen in einer \method{prepareLattice} Methode
	\item Starten der Simulationsschleife zum Aufruf von \method{SuperLattice2D::collideAndStream}
\end{enumerate}

In letzterem, die eigentliche Simulation durchführendem, Schritt, werden weiter durch kanonisch benannte Funktionen wie \method{getResults} und \method{error} die Ergebnisse zur Analyse in Dateien~\cite{vtkGuide10} geschrieben, Fehlernormen berechnet und Konvergenzkriterien bestimmt.

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

Ein erster Gedanke zur Integration von Gitterverfeinerung in OpenLB ist die Nutzung der bestehenden Dekomposition der Simulationsdomäne in achsenparallele Rechtecke. Insbesondere aus Sicht des Einfügens von Gitterverfeinerung in die existierende Architektur, sowie der unveränderten Weiterverwendung der \class{LoadBalancer} Algorithmen zur Steuerung der Parallelisierung, scheint ein solcher Ansatz sinnvoll.

Bei Variation der Auflösung einzelner Quader im Rahmen der \class{CuboidGeometry2D} Struktur handelte es sich zwangsweise um einen Multi-Domain Ansatz. Gingen wir diesen Weg, benötigten wir zunächst \class{Cuboid2D} spezifische \class{UnitConverter} Instanzen zur Verwaltung der auflösungsabhängigen Konstanten. Dies müsste dann im Rahmen von \method{prepareLattice} zur Setzung der dann ebenfalls quaderspezifischen \class{Dynamics} und Randkonditionen beachtet werden.

Zur Ermöglichung von Parallelisierung berücksichtigt die, der Gitterverwaltung in \class{SuperLattice2D} zugrundeliegende, Aufteilung der Domäne durch \class{CuboidGeometry2D} bereits Übergangsbereiche, deren Funktion mit zusätzlicher Auflösungskopplung in Einklang zu bringen wäre.

Weiter würde das Problemfeld eben dieser Dekomposition um die Restriktion auf Auflösungsübergänge im Verhältnis \(1:2\) erweitert. So müsste ein guter Algorithmus zur Aufteilung der Simulationsdomäne die Anforderungen an Parallelisierung, Verteilung der Rechenressourcen, Geometrie und Verfeinerung sinnvoll auflösen und zugleich manuelle Eingriffe erlauben. Diese zusätzliche starke Einschränkung sowie der dann bei Anpassung der Verfeinerungsstruktur unumgängliche komplette Neuaufbau der Simulation bilden ein starkes Gegenargument zu diesem ersten Gedanken.

\bigskip
Der tatsächlich umgesetzte Ansatz ergibt sich aus dem Verständnis von Gitterverfeinerung als Kopplung von ansonsten komplett allein stehenden Simulationen. Die Übergangsbereiche wären in diesem Modell mit Randkonditionen vergleichbar, wie sie für Ein- und Ausflüsse verwendet werden. Gitterverfeinerung könnte so weitestgehend von der bestehenden Architektur getrennt ergänzt werden, was insbesondere auch die veränderungsfreie Unterstützung existierender Anwendungen begünstigen würde.

Eine solche nebenläufige Überlagerung von Simulationen mit jeweils komplett eigenständig verwalteten Gittern gebietet sich bei erster Betrachtung als klarer Umriss eines Multi-Grid Verfahrens. Beachten wir jedoch, dass es einfach möglich ist, die überlagerten Gitterflächen durch \emph{Nullen} der entsprechenden Materialzahlen effizient aus der Verarbeitung auszuschließen und trotzdem bei Bedarf -- z.B. in Hinblick auf Verschiebung von verfeinerten Bereiche während des Simulationsverlaufs -- zu reaktivieren, stellen sich auch Multi-Domain Ansätze in diesem Modell als sinnvoll implementierbar heraus. Vorteil ist hier also gerade auch, dass prinzipiell beide Ansätze zur Gitterverfeinerung umgesetzt werden können und wir nicht durch Festhalten an der existierenden Struktur auf Multi-Domain Verfahren beschränkt sind. Da die Positionierung der Gitter in diesem Ansatz komplett frei wäre, ließen sich aus Architektursicht auch nicht-koinzidierende oder sogar zueinander rotierte Verfeinerungsgitter abbilden.

Ein Vorbild für dieses Konzept zur Umsetzung von Gitterverfeinerung existiert in Form der Optimierungskomponente von OpenLB, welche ebenfalls komplette Simulationen in einem sogenannten \class{Solver} kapselt \cite[vgl.~Abb.~4.1]{Krause10}. Langfristig könnten mit diesem Ansatz also beide gitterübergreifenden Module in einem gemeinsamen Konzept abgebildet werden.

\bigskip
Nachdem nun das grobe Umfeld eines Gitterverfeinerungsframeworks feststeht, gilt es, ein geeignetes Verfahren zur Umsetzung in und Nutzung mit eben diesem Framework zu wählen. Das von Lagrava et al. in \citetitle{Lagrava12}~\cite{Lagrava12} beschriebene Verfahren, welches insbesondere auf \cite{DupuisChopard03} und \cite{Filippova98} einen anpassungsfähigen Multi-Domain Gitterverfeinerungsalgorithmus für BGK LBM auf koinzidierenden D2Q9 Gittern aufbaut, erscheint hier als guter Kandidat. Die anfängliche Beschränkung auf zwei Dimensionen passend zur Einschränkung dieser Arbeit sowie die Flexibilität in Hinblick auf die verwendeten Restriktions- und Interpolationsoperatoren bilden hier eine gute Grundlage für eine erste und doch ausbaufähige Umsetzung von Gitterverfeinerung in OpenLB.

\newpage
\subsection{Struktur des Gitterverfeinerungsframework}

Wie im vorangehenden Kapitel erläutert, soll das Framework zur Gitterverfeinerung nicht tief in vorhandene Strukturen integriert, sondern viel mehr über diesen stehend angesiedelt werden. Eine erste Hürde zu diesem Ziel ist die, größtenteils aus zwangfreien Konventionen bestehende, Struktur von OpenLB Anwendungen. So sind zwar die einzelnen Komponenten der Simulation wie \class{CuboidGeometry2D} und \class{SuperLattice2D} vorgegeben, deren Konstruktion und Verknüpfung erfolgt jedoch größtenteils manuell.

Für sich ist diese Herangehensweise des flexiblen Zusammensetzens von Modulen durchaus erhaltenswert und bildet eine der Stärken von OpenLB. Zur übergreifenden und für den Nutzer möglichst bequemen Einbindung von Gitterverfeinerung -- wir erinnern uns: das Ziel ist es, Gitter erst im Nachhinein mit einem einzigen Funktionsaufruf zu verfeinern -- muss jedoch zumindest die Konstruktion des auflösungseigenen \class{SuperLattice2D} mit dem dazugehörigen Umfeld aus \class{UnitConverter}, \class{LoadBalancer}, \class{CuboidGeometry2D} und \class{SuperGeometry2D} soweit wie möglich gekapselt werden.

\bigskip
Entsprechend besteht das Framework aus zwei Komponenten: Einer \class{Grid2D} Klasse, die in einem Konstruktoraufruf ein \class{SuperLattice2D} zusammen mit dem benötigten Umfeld instanziiert und einer \class{Coupler2D} Klasse zur Abbildung der Übergänge zwischen mehreren \class{Grid2D} Instanzen. Die Gitterklasse stellt dabei eine Methode \method{Grid2D::refine} bereit, die anhand einer Parametrisierung der zu verfeinernden Domäne ein neues \method{Grid2D} konstruiert und über entsprechende \class{Coupler2D} Objekte mit sich selbst verknüpft. Funktionen wie \method{prepareGeometry} und \method{prepareLattice} können in diesem Umfeld dann durch entsprechende \emph{Getter} mit \class{Grid2D} zusammenarbeiten. Arbeiten diese Funktionen bereits auf Grundlage von analytischen Indikatoren, d. h. unabhängig der Auflösung, können sie sogar ohne Anpassung für alle Gitterauflösungen verwendet werden.

\begin{listing}[H]
\begin{minted}{cpp}
template <typename T, template<typename> class DESCRIPTOR>
RefiningGrid2D<T,DESCRIPTOR>::RefiningGrid2D(
	Grid2D<T,DESCRIPTOR>& parentGrid, Vector<T,2> origin, Vector<T,2> extend):
	Grid2D<T,DESCRIPTOR>(
		std::unique_ptr<IndicatorF2D<T>>(new IndicatorCuboid2D<T>(extend, origin)),
		2*parentGrid.getConverter().getResolution(),     // Auflösungsübergang $2:1$
		2*parentGrid.getConverter().getLatticeRelaxationTime() - 0.5, // Siehe $(\ref{eq:gridTauShift})$
		parentGrid.getConverter().getReynoldsNumber()),               // $\text{Re}_g = \text{Re}_f$
	_origin(origin),
	_extend(extend),
	_parentGrid(parentGrid)
{ }
\end{minted}
\caption{Konstruktor der verfeinernden Gitter}
\label{lst:RefiningGrid}
\end{listing}

Die Konstruktion von \class{Grid2D} erfolgt anhand einer indikatorgegebenen Beschreibung der Simulationsdomäne sowie der gewünschten Auflösung zusammen mit der Relaxionszeit und der modellierenden Reynolds-Nummer:
\begin{minted}{cpp}
Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF, int resolution, T tau, int re);
\end{minted}
Während sich die Realisierung dieser Signatur als einfache Konstruktion der erläuterten OpenLB Struktur erweist, gestaltet sich die Konstruktion der vererbten \class{RefiningGrid2D} Klasse in Listing~\ref{lst:RefiningGrid} interessanter, da hier dann Kraft der Ergebnisse von Kapitel~\ref{kap:Skalierung} die Vorgabe des groben Gitters zusammen mit dem verfeinerungsbedürftigen Teilbereich zur Erstellung eines neuen Gitters ausreicht.

\begin{listing}[H]
\inputminted{cpp}{code/grid2d_collide_and_stream.cpp}
\caption{Rekursiver Kollisions- und Strömungsschritt mit Gitterkopplung}
\label{lst:GridCollideAndStream}
\end{listing}

Wie in Kapitel~\ref{kap:Algorithmus} dargelegt, müssen zur Gitterkopplung nach jedem Kollisions- und Strömungsschritt verschiedene Arbeiten durchführt werden. So ist die Ausführung von Kollisions- und Strömungsschritten auf dem feinen Gitter zusammen mit der jeweiligen Vor- und Nacharbeit strikt an die Nacharbeit von Kollisions- und Strömungsschritten auf dem groben Gitter gebunden.

Wurde das grobe Gitter um einen Zeitschritt weiterentwickelt, muss der Zustand des feinen Gitters ebenfalls um entsprechend zwei feine Zeitschritte weiterentwickelt werden, damit die groben Verteilungsfunktionen vervollständigt werden können. Es liegt somit nahe, die Aufrufe von \method{SuperLattice2D::collideAndStream} in einer Methode \method{Grid2D::collideAndStream} zu kapseln, um auf diese Weise die Aufrufe von \class{Coupler2D} an den korrekten Stellen durchzuführen.

Konkret erhalten wir in Listing~\ref{lst:GridCollideAndStream} bei gesammelter Verwaltung der von \method{Grid2D::refine} erstellten feinen Gitter und den zugehörigen Kopplern eine, der Algorithmenübersicht in Abbildung~\ref{fig:AlgorithmBirdsEye} nicht unähnliche, Implementierung von \method{Grid2D::collideAndStream}.
Zu bemerken ist, dass die Konstellation aus dieser Methode zusammen mit \method{Grid2D::refine} durch Selbstaufruf bereits die freie Schachtelung von Verfeinerungsbereichen erlaubt.

\begin{listing}[H]
\begin{minted}{cpp}
// Initialisiere gröbstes Gitter mit gewünschten Fluidkonstanten
Grid2D<T,DESCRIPTOR> coarseGrid(coarseDomain, resolution, tau, Re);
prepareGeometry(coarseGrid);

// Verfeinere ein Einheitsquadrat beginnend bei $(0.5,0.5) \in \R^2$
auto fineGrid = coarseGrid.refine({0.5, 0.5}, {1.0, 1.0});
prepareGeometry(fineGrid);

// Schließe den feinen Bereich auf dem groben Gitter aus (optional)
auto refinedOverlap = fineGrid.getRefinedOverlap();
coarseGrid->getSuperGeometry().rename(1,0,*refinedOverlap);

// Binde Dynamics und Randkonditionen an die beiden Gitter
prepareLattice(coarseGrid);
prepareLattice(fineGrid);

// Simulationsschleife mit Ausgabe
for (int iT = 0; iT < coarseGrid->getConverter().getLatticeTime(100); ++iT) {
	coarseGrid->collideAndStream();

	getResults(coarseGrid, iT);
	getResults(fineGrid,   iT);
}
\end{minted}
\caption{Beispielhafte Nutzung von \class{Grid2D}}
\label{lst:RefinementUsageExample}
\end{listing}

Wie in diesem Listing zu sehen, kann Gitterverfeinerung innerhalb des beschriebenen Frameworks schon erfreulich kompakt formuliert werden. Tatsächlich fehlt zum Etappenziel der einzeilig aktivierbaren manuell positionierten Gitterverfeinerung nur eine weitere Abstraktion von von \method{prepareGeometry} und \method{prepareLattice}.

\bigskip
Zur umfassenden Beschreibung des Gitterverfeinerungsframework fehlt uns jetzt nur noch die Definition der \class{Coupler2D} Objekte. Da diese die Einzelheiten des Verfeinerungsverfahrens umsetzen, werden sie im Rahmen des folgenden Kapitels zur Implementierung des Verfahrens von Lagrava et al. näher beleuchtet werden.

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

Grundsätzlich implementiert jede Instanz von \class{Coupler2D} die Kopplung zweier Gitter in einer Richtung entlang einer, durch Ursprung und Ausdehnung charakterisierten, Linie innerhalb der physikalischen Simulationsdomäne. Für die Kopplung einer rechteckigen \class{RefiningGrid2D} Instanz werden von \method{Grid2D::refine} in diesem Fall acht Kopplungsobjekte erzeugt -- vier Seiten mit jeweils zwei Kopplern.

\begin{listing}[H]
\inputminted{cpp}{code/coupler2d.cpp}
\caption{Gemeinsame Struktur beider Kopplungsklassen}
\label{lst:Coupler2D}
\end{listing}

Die im Zuge dieser Arbeit entwickelte Version von \class{Coupler2D} beschränkt sich hierbei auf zu einem Einheitsvektor parallele Gitterübergänge. Sowohl für die Kopplung der mit \class{CuboidGeometry2D} modellierbaren Aufteilungen als auch für das umzusetzende Verfahren ist diese Einschränkung kein Hindernis, da Lagrava et al. ebenfalls nur von horizontalen bzw. vertikalen Gitterübergängen ausgehen.

Da \class{Grid2D} Methoden zur Diskretisierung physikalischer Koordinaten auf das Gitter bereitstellt, besteht das Fundament der beiden benötigten Kopplungsklassen größtenteils nur aus der Bestimmung aller zu setzenden Kopplungsknoten entlang der Übergangslinie. Die abgeleiteten Klassen \class{FineCoupler2D} und \class{CoarseCoupler2D} können auf die Liste dieser Knoten dann mittels \method{getFineLatticeR} und \method{getCoarseLatticeR} zugreifen und so ihre eigene Implementierung auf das Wesentliche beschränken.

\begin{listing}[H]
\inputminted{cpp}{code/fineCoupler2d.cpp}
\caption{Struktur des Kopplers von grob nach fein}
\label{lst:FineCoupler2D}
\end{listing}

Wie in Listing~\ref{lst:FineCoupler2D} zu sehen, benötigt das Setzen der feinen Verteilungsfunktionen in \method{FineCoupler2D::couple} einen Zwischenspeicher der groben Verteilungsmomente und der groben Nicht-Equilibriumsverteilung entlang der Kopplungslinie. Dieser wird in der Methode \method{FineCoupler2D::store} gesetzt und dient in \method{FineCoupler2D::interpolate} der linearen Zeitinterpolation der groben Verteilungsfunktionen zu Zeitpunkt \(t+\delta t_f\) entsprechend der Ausführungen in Kapitel~\ref{kap:Algorithmus}.

\begin{listing}[H]
\inputminted{cpp}{code/fineCoupler2d_couple_extract.cpp}
\caption{Ausschnitt der Methode \method{FineCoupler2D::couple}}
\label{lst:FineCoupler2D_couple_extract}
\end{listing}

Zur Beleuchtung des Herzstücks der feinen Kopplung sehen wir in Listing~\ref{lst:FineCoupler2D_couple_extract} einen Ausschnitt der Kopplungsfunktion für Knoten \(x_{g \to f}^f \in \F\) mit räumlicher Interpolation der benötigten groben Werte. Zusammen mit der zugehörigen Interpolationsfunktion in Listing~\ref{lst:ipol4ord} lässt sich dabei sehen, wie sich das Verfeinerungsvefahren durch Verwenden von OpenLB Modulen wie den \class{lbHelpers} und der \class{Vector} Datenstruktur sehr nah an seiner mathematischen Formulierung umsetzen lässt.

\begin{listing}[H]
\begin{minted}{cpp}
template <typename T, unsigned N>
Vector<T,N> order4interpolation(const std::vector<Vector<T,N>>& data, int y)
{
	return 9./16. * (data[y] + data[y+1]) - 1./16. * (data[y-1] + data[y+2]);
}
\end{minted}
\caption{Templatefunktion der Interpolationsformel (\ref{eq:ipol4ord})}
\label{lst:ipol4ord}
\end{listing}

Der Koppler \class{CoarseCoupler2D} zum Setzen der groben Verteilungsfunktionen gestaltet sich derweil einfacher, da kein Zwischenspeicher benötigt wird und nur eine Restriktionsoperation anzuwenden ist.

\begin{listing}[H]
\inputminted{cpp}{code/computeRestrictedFneq.cpp}
\caption{Umsetzung der Restriktionsoperation (\ref{eq:neqAvgRestrictionF2G})}
\label{lst:CoarseCoupler2D_restriction}
\end{listing}

Hiermit sind die zentralen Bestandteile der Umsetzung des Verfahrens von Lagrava et al. im Kontext des in dieser Arbeit entwickelten Gitterverfeinerungsframework für OpenLB beschrieben. Zum Abschluss verbleibt nun noch die Evaluierung der Qualität eben dieses Verfahrens anhand verschiedener Beispiele.

\newpage
\section{Evaluierung}

Die Auswahl von Beispielen und Kriterien zur Bewertung der Qualität eines Gitterverfeinerungsverfahrens gestaltet sich zunächst unklarer, als man annehmen könnte. So stellt OpenLB zwar eine gute Sammlung verschiedener Simulationsbeispiele bereit, aber nur wenige von ihnen besitzen analytische Lösungen oder einfach zu verwendende gesicherte Vergleichsdaten. Solche Referenzwerte sind -- abseits offensichtlicher Gütekriterien wie der Vermeidung divergierender Simulationen -- Voraussetzung für die Bewertung der Verfahrensqualität. Ohne diese ist beim Vergleich mit aus uniformen Gittern gewonnenen Lösungen nicht klar, welches Ergebnis besser ist.

Entsprechend beschränken wir uns je nach Beispiel auf die Betrachtung einer Auswahl der folgenden Gütekriterien:
\begin{enumerate}
	\item Subjektive Qualität des Strömungsbildes
	\item Stetigkeit der Momente im Gitterübergang
	\item Bestimmung des Fehlers zu einer analytischen Lösung
	\item Anwendung eines Gitterverfeinerungskriteriums
	\item Vergleich von lokal verfeinerten und uniformen Gittern mit gleicher Knotenanzahl
\end{enumerate}

Das hochwertigste Kriterium ist dabei der Fehlervergleich zur analytischen Lösung bei Beschränkung der Gesamtknotenanzahl. Existiert die dazu notwenige formalisierte Lösung nicht, sind auch Vergleiche mit Referenzwerten aus realen Experimenten oder gesicherten Simulationen denkbar. Die subjektive Qualitätsbewertung hingegen erlaubt nur die Bewertung von Extremfällen wie im Strömungsbild erkennbaren Übergängen zwischen Gittern oder möglicher Divergenzvermeidung durch gezielte Verfeinerung.

\newpage
\subsection{Rohrströmung}

Als Einstiegspunkt wollen wir zunächst die grundsätzliche Funktion des Verfeinerungsverfahrens an einem möglichst einfachen Beispiel gegen möglichst korrekte Vergleichsdaten testen. Ein solches Beispiel ist gegeben durch die laminare Strömung in einem Rohr mit kreisförmigem Querschnitt und ohne Hindernisse abseits der Wände.

\begin{figure}[h]
\centering
\input{img/poiseuille2d_overview.tikz}
\caption{Schematische Darstellung der Rohrströmung}
\label{fig:PoiseuilleOverview}
\end{figure}

Diese Rohrströmung stellt nicht nur eine der einfachsten Strömungssituationen dar, sondern besitzt als Poiseuille-Fluss auch eine analytische Lösung, so dass wir ideale Vergleichsbedingungen vorfinden. Lieferte unser Verfahren in diesem Beispiel keine guten Ergebnisse, wäre nicht davon auszugehen, dass dies sich in komplexeren Situationen verbessern würde.

\begin{Definition}[Analytische Lösung des Poiseuille-Flusses]
\label{def:analyticPoiseuille}
Seien \(L_x, L_y \in \R_+\) die räumlichen Rohrdimensionen, \(\nu\) die kinematische Viskosität und \(\Delta p\) der Druckgradient zwischen Ein- und Ausfluss. Dann ist die analytische Geschwindigkeit in \(x\)-Richtung gegeben als \cite[vgl.~Kap.~4]{Bao11}:
\[u_x(y) = \frac{1}{2\nu} \frac{\Delta p}{L_x} y (y-L_y)\]
Dies ergibt mit \(u_x(L_y/2):=u_\text{max}\) und \(\Delta p := -p_0\) die analytische Lösung des Drucks:
\[p_0 = \frac{8 L_x \nu u_\text{max} }{L_y^2}\]
Mit \(u_\text{max}:=1\) vereinfacht sich damit die analytische Lösung der \(x\)-Geschwindigkeit zu:
\[u_x(y) = -\frac{4}{L_y^2} y (y-L_y)\]
Das Geschwindigkeitsprofil des Poiseuille-Flusses ist also parabelförmig.
\end{Definition}

Wir wollen in einem \(1 \times 4\) Meter bemessenden Rohr einen Poiseuille-Fluss simulieren. Als Auflösung einer Längeneinheit sei dabei \(N=10\) gewählt, was in der Diskretisierung durch \(11 \times 21\) grobe und \(21 \times 43\) feine Knoten abgebildet wird.
In Abbildung~\ref{fig:PoiseuilleGridSetup} sehen wir das resultierende Gitter zusammen mit den zugewiesenen Materialzahlen. Wand- und Einflusszellen werden nach dieser Vorlage mit lokalen Geschwindigkeitsrandbedingungen umgesetzt. Während für den Einfluss dabei das Geschwindigkeitsprofil der analytischen Poiseuille-Lösung vorausgesetzt wird, erhält der Ausfluss eine Druckrandbedingung.

\begin{figure}[H]
\begin{adjustbox}{center}
\input{img/poiseuille2d_grid.tikz}
\end{adjustbox}
\caption{Gitterstruktur einer halbseitig verfeinerten Rohrströmung}
\label{fig:PoiseuilleGridSetup}
\end{figure}

Neben diesen knotenspezifischen Eigenschaften sei \(u=\num{0.01}\) die charateristische Geschwindigkeit in Lattice-Einheiten und \(\text{Re}=10\) die modellierte Reynolds-Zahl. Erstellen wir unsere grobe \class{Grid2D} Instanz mit diesen, die Relaxionszeit \(\tau\) fixierenden, Werten und führen die Simulation bis zur Konvergenz aus, erblicken wir nach geeigneter Aufbereitung der \class{VTK}-Ausgabe \cite[Kap.~19.3]{vtkGuide10} schließlich das in Abbildung~\ref{fig:PoiseuilleVelocityGrid} dargestellte Strömungsbild. Konvergenz bedeutet in diesem Kontext, dass die durchschnittliche Energie des feinen Gitters unter einen Residuumswert, hier \num{1e-5}, gefallen ist.

\begin{figure}[h]
\begin{adjustbox}{center}
\input{img/poiseuille2d_velocity_grid.tikz}
\end{adjustbox}
\caption{Geschwindigkeiten in \(x\)-Richtung im simulierten Poiseuille-Fluss}
\label{fig:PoiseuilleVelocityGrid}
\end{figure}

Bei erster Betrachtung lässt sich erkennen, dass die Strömung den Gitterübergang subjektiv ideal bestritten hat. Es treten also keine ungewöhnlichen Artefakte im Geschwindigkeitsbild auf und dieses setzt sich nach dem Übergang in, bis auf die neuen Zwischenwerte, unveränderter Weise fort. Tatsächlich ist bei Bildung einer geschlossenen Fläche durch Interpolation der Zwischenbereiche kein Gitterübergang erkennbar.

\newpage
\subsubsection{Vergleich mit der analytischen Lösung}

Zur formalen Qualitätsbewertung ziehen wir im Folgenden die analytische Lösung von Geschwindigkeit und Druck des Poiseuille-Flusses in Definition~\ref{def:analyticPoiseuille} heran. Wir können diese in \mbox{OpenLB} einfach mit Hilfe des \class{SuperRelativeErrorL2Norm2D} Funktors auf beiden Gittern mit der simulierten Lösung vergleichen:

\begin{figure}[h]
\centering
\input{img/poiseuille2d_error_comparison.tikz}
\caption{Geschwindigkeits- und Druckfehler in verschiedenen Gittern}
\label{fig:PoiseuilleErrorNorm}
\end{figure}

Zunächst ist die halbseitig verfeinerte Lösung aus Abbildung~\ref{fig:PoiseuilleGridSetup} also um etwa eine Größenordnung schlechter als eine gleichmäßig mit \(n=20\) aufgelöste Berechnung. Auch im Vergleich mit dem uniform \(n=10\) aufgelösten Gitter erweist sich der Fehler der verfeinerten Lösung als deutlich größer -- zumindest in diesem speziellen Beispiel ist Gitterverfeinerung also der Simulationsgenauigkeit messbar abträglich.

Nicht vergessen werden sollte jedoch, dass die untersuchte halbseitig verfeinerte Rohrströmung als Beispiel sehr konstruiert und nicht realitätsnah ist. Auch die noch folgenden Beobachtungen in Abbildung~\ref{fig:PoiseuilleOutflowProfile}, nach welchen die lineare Interpolation zu kleineren Geschwindigkeitsfehlern führt, deutet auf eine beschränkte Aussagekraft dieses Beispiels hin.
Eine Besonderheit ist in diesem Kontext auch die Existenz von Randbedingungen im Übergangsbereich. Eine etwaige Behandlung dieser wurde weder von Lagrava et al. angesprochen, noch in der dem Leser vorliegenden Arbeit näher untersucht. Tatsächlich verschwindet der \emph{Verfeinerungsfehler} fast vollständig, wenn die Wände aus dem verfeinerten Bereich ausgespart werden.

\bigskip

Abschließend erscheint es beispielübergreifend intuitiv erwartbar, dass eine nicht aus dem konkreten Strömungsproblem informierte Anwendung von Gitterverfeinerung -- und damit eine unbegründete Erhöhung der Simulationskomplexität gegenüber einem uniformen Gitter -- einer Verbesserung der Präzision nicht zuträglich ist.

\subsubsection{Vergleich der Interpolationsverfahren}

Den halbseitig verfeinerten Poiseuille Simulationsaufbau können wir an dieser Stelle auch zur Nachvollziehung des, von Lagrava et al. für die Verwendung eines Verfahrens vierter Ordnung in der räumlichen Interpolation feiner Übergangsknoten ohne koinzidierende grobe Stützpunkte hervorgebrachten, Arguments verwenden.

\begin{figure}[h]
\begin{adjustbox}{center}
\input{img/massloss_interpolation_plot.tikz}
\end{adjustbox}
\caption{Druckverlauf bei linearer und kubischer Interpolation \cite[vgl.~Abb.~11]{Lagrava12}}
\label{fig:PoiseuilleMassloss}
\end{figure}

Wir setzen dazu \(N=50\) als Auflösung der Längeneinheit, \(\text{Re}=100\) als Reynolds-Zahl und eine Geschwindigkeitsrandbedingung mit Poiseuilleprofil im Ausfluss. Führen wir dann die Simulation mit dem linearen Interpolationsverfahren (\ref{eq:ipol2ord}) aus und vergleichen den Verlauf des physikalischen Drucks auf einer vertikal zentrierten horizontalen Linie mit den, aus einem Durchlauf mit dem Verfahren vierter Ordnung (\ref{eq:ipol4ord}) gewonnen, Daten, erhalten wir den in Abbildung~\ref{fig:PoiseuilleMassloss} zu sehenden Plot.

Entsprechend der Beobachtungen in \cite[Kap.~3.7]{Lagrava12} sehen auch wir bei linearer Interpolation einen prominenten Abfall des physikalischen Drucks im Übergangsbereich der Gitter. Bei kubischer Interpolation tritt dieser Fehler nicht auf, der Druckverlauf ist in diesem Fall so glatt, dass der Übergang nicht mehr zu erkennen ist.

\begin{figure}[h]
\begin{adjustbox}{center}
\input{img/poiseuille2d_velocity_outflow.tikz}
\end{adjustbox}
\caption{Vergleich des vertikalen Geschwindigkeitsprofil bei \(x=3\)}
\label{fig:PoiseuilleOutflowProfile}
\end{figure}

Zum Abschluss dieses Beispiels vergleichen wir in Abbildung~\ref{fig:PoiseuilleOutflowProfile} das Geschwindigkeitsprofil entlang einer die gesamte Rohrbreite durchgeschreitenden vertikalen Linie bei \(x=3\) im halbseitig verfeinerten Fall. Wir sehen hier zunächst die subjektiv beurteilt gute Qualität des Strömungsbildes erneut bestätig, erkennen aber auch den größeren Fehler im Vergleich mit der uniform fein aufgelösten Simulation. Interessant ist hier, dass die den Druckabfall ausgleichende kubische Interpolation zugleich den Geschwindigkeitsfehler gegenüber der linearen Interpolation leicht zu erhöhen scheint.

Legen wir jedoch erneut ein randlos halbseitig verfeinertes Gitter wie beim Vergleich der L2-Normen in Abbildung~\ref{fig:PoiseuilleErrorNorm} an, verschwindet dieses Problem und wir erhalten bei dieser Messung sogar einen leicht geringeren Fehler als im uniform fein aufgelösten Gitter.

\bigskip

Fassen wir die Ergebnisse der zurückliegenden ersten Analyse des Gitterverfeinerungsverfahrens unter dem Eindruck der beschränkten Aussagekraft bezogen auf die Wahl der Rohrströmung zusammen, bestätigt sich das Verfahren als grundsätzlich anwendbar. Es treten also sowohl bei subjektiver Betrachtung als auch bei Vergleich der analytischen Lösung keine extremen Fehler im Strömungsbild auf. Der beobachtete negative Einfluss von Randbedingungen im Übergangsbereich der Gitter sollte hier aber nicht unerwiedert bleiben -- bis die korrekte Behandlung solcher Randbedingungen nicht geklärt ist, sollten Gitterübergänge nur im offenen Fluid platziert werden, denn nur dessen Skalierung wird in Kapitel~\ref{kap:Lagrava} und der zugrundeliegendem Arbeit \cite{Lagrava12} behandelt.

\newpage
\subsection{Umströmter Zylinder}

Bei dem \emph{umströmten Zylinder} handelt es sich um ein verbreitetes Strömungsbeispiel welches entsprechend in der Menge der OpenLB Beispielanwendungen enthalten ist. Grundsätzlich ähnelt es dem Aufbau der Rohrströmung -- simuliert wird die von zwei Wänden begrenzte Strömung zwischen Ein- und Ausfluss ergänzt um ein zylindrisches Hindernis im Eingangsbereich.

\begin{figure}[h]
\centering
\input{img/cylinder2d_overview.tikz}
\caption{Schematische Darstellung des umströmten Zylinder \cite[vgl.~Abb.~1]{SchaeferTurek96}}
\label{fig:cylinder2d_overview}
\end{figure}

Während für diese Strömungssituation noch keine analytische Lösung gefunden wurde, stehen in \citetitle{SchaeferTurek96}~\cite{SchaeferTurek96} detaillierte, hochwertige und mit verschiedenen Verfahren berechnete Vergleichsdaten zur Verfügung. Bevor wir diese jedoch zur Evaluation heranziehen, wollen wir uns vorerst der subjektiven Qualität der Ergebnisse versichern und uns nähere Gedanken zur Lokalisierung der Verfeinerunsbereiche machen.

\begin{figure}[H]
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_unrefined_n20_re100_16s.pdf}
\end{adjustbox}
\caption{Uniform mit \(N=20\) aufgelöstes Strömungsbild zu \(t=16s\)}
\label{fig:UniformCylinderVelocity16s}
\end{figure}

Für die Umsetzung in OpenLB parametrisieren wir die Geometrie bezogen auf den Zylinderdurchmesser \(D\) und dimensionalisieren diesen wiederum als \(D := 0.1m\), was zugleich der charakteristischen Länge entspreche. Auflösungsangaben beziehen sich im Folgenden also auf den Durchmesser des Zylinders in groben Gitterweiten. Hinblickend auf die Vorgaben zum instationären Testfall \cite[Kapitel~2.2b]{SchaeferTurek96} sei \(\text{Re}:=100\) die Reynolds-Zahl und für den Einfluss sei ein Poiseuille-Geschwindigkeitsprofil angelegt. Wände und Ausflüsse werden analog zur hindernisfreien Rohrströmung durch lokale Geschwindigkeits- bzw. Druckrandbedingungen konstruiert, während der Zylinder den Fluss durch Bounce-Back hindere. Eine Relaxionszeit \(\overline\tau_g := 0.51\) des gröbsten Gitters vervollständigt unser Modell.

\begin{figure}[h]
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_single_refinement_n20_re100_16s.pdf}
\end{adjustbox}
\caption{Einfach verfeinertes Strömungsbild zu \(t=16s\)}
\label{fig:SingleLevelRefinementCylinderVelocity16s}
\end{figure}

Als Grundlage für den subjektiven Vergleich des Strömungsbildes simulieren wir zunächst in Abbildung~\ref{fig:UniformCylinderVelocity16s} auf einem unverfeinert mit \(N=20\) aufgelösten Gitter. Charakteristisch ist hier direkt die Bildung einer Kármánschen Wirbelstraße zu beobachten.

\noindent
Vergleichen wir diese Grundsituation mit der in Abbildung~\ref{fig:SingleLevelRefinementCylinderVelocity16s} zu sehenden, um den Zylinder herum verfeinerten, Simulation, wirkt das Strömungsbild subjektiv gut aber unterschiedlich: Während positiv auffällt, dass der der Gitterübergang im Geschwindigkeitsbild trotz komplexerer Strömungsstruktur nicht zu erkennen ist, unterscheidet sich die Position der Wirbel trotz gleichem Zeitpunkt \(t=16s\) erkennbar.

\begin{figure}[H]
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_unrefined_n40_re100_16s.pdf}
\end{adjustbox}
\caption{Uniform mit \(N=40\) aufgelöstes Strömungsbild zu \(t=16s\)}
\label{fig:UniformCylinderVelocityN4016s}
\end{figure}

Ergänzen wir unsere Auswahl von Geschwindigkeitsbildern jedoch um Abbildung~\ref{fig:UniformCylinderVelocityN4016s}, welche aus einem unverfeinerten \(N=40\) Gitter hervorgeht, ist oberflächlich gegenüber dieser kein Unterschied zur einfach verfeinerten Variante erkennbar. Beachtenswert ist dabei, dass das uniform mit \(N=40\) aufgelöste Gitter \(\sim 145000\) Knoten enthält, während das subjektiv identische lokal verfeinerte \(N=20\) Gitter mit \(\sim 66000\) nicht einmal halb so viele Knoten benötigt.

\bigskip

Noch weiter kann diese Beobachtung hier jedoch nicht bewertet werden, da nicht klar ist, welche Wirbelkonfiguration in dieser konkreten Strömungssituation korrekt ist. Aussagekräftiger für die formale Qualitätsbewertung wird der Vergleich der von Schäfer und Turek in \cite{SchaeferTurek96} zusammengestellten Referenzwerte sein.

\subsubsection{Anwendung eines formalen Kriteriums zur Gitterverfeinerung}

An dieser Stelle wollen wir die Gelegenheit nutzen, uns näher mit der Wahl der Verfeinerungsbereiche auseinanderzusetzen. So scheint es auf der einen Seite intuitiv sinnvoll, einen Verfeinerungsbereich in besonders komplexen Regionen einer Simulationsdomäne zu platzieren -- in diesem Fall also um den Zylinder, dessen Geometrie dann u.a. feiner aufgelöst würde. Auf der anderen Seite ist eine solch intuitive Entscheidung aber leicht fehlbar und schwer qualitativ zu beurteilen oder zu optimieren. Entsprechende Gedanken bewogen auch die Autoren unseres Verfeinerungsverfahrens, in \citetitle{Lagrava15}~\cite{Lagrava15} ein auf der Knudsen-Zahl basierendes automatisches Kriterium zur Gitterverfeinerung herzuleiten.

\begin{Definition}[Knudsen-Zahl]
Sei \(\lambda\) die mittlere freie Weglänge, \(L\) die charakteristische Länge des Systems. Die Knudsen-Zahl \(\text{Kn}\) ist definiert als:
\[ \text{Kn} := \frac{\lambda}{L} \quad \left(\overset{\text{\cite[Kap.~1.4.2]{Haenel04}}}{=} \frac{\text{Ma}}{\text{Re}} \frac{\rho^2 c_s \lambda}{\nu} \quad \overset{\text{\cite[Gl.~(7.22)]{Krueger17}}}{\simeq} \frac{\text{Ma}}{\text{Re}}\right) \]
Diese Knudsen-Zahl ist eine dimensionlose Kennzahl der Strömungslehre. Ihr Wert beschreibt, ob ein Gas als strömungsmechanisches Kontinuum nach Navier-Stokes oder als Bewegung einzelner Teilchen betrachtet werden kann \cite[S.~14]{Krueger17}:
\begin{description}
\item[\(\text{Kn} \ll 1\)] Strömungsmechanisches Kontinuum nach Navier-Stokes
\item[\(\text{Kn} \gtrsim 1\)] Betrachtung einzelner Teilchen
\end{description}
\end{Definition}

\begin{Definition}[Auftreten der Knudsen-Zahl in LBM]
Die Knudsen-Zahl wird im LBM-Kontext nach \cite[vgl.~(21)]{Lagrava15} durch den Quotienten der Nicht-Equilibriums- und Equilibriumsverteilung in Richtung \(i \in \{0,\dots,q-1\}\) genähert:
\[\text{Kn} \sim \frac{f_i^\text{neq}}{f_i^\text{eq}} \overset{!}{\ll} 1 \]
Die Mittelung über alle Richtungen ergibt so die lokal in einem Knoten des Gitters genäherte Knudsen-Zahl:
\[C(x) := \frac{1}{q} \sum_{i=0}^{q-1} \left|\frac{f_i^\text{neq}}{f_i^\text{eq}}\right|\]
\end{Definition}

\begin{Definition}[Lokaler Verfeinerungsfaktor]
Nach \cite[vgl.~(29)]{Lagrava15} führt die folgende Transformation der lokal genäherten Knudsen-Zahl zur Berechnung eines diskreten Verfeinerungsfaktors:
\[ R(x) := \text{round}\left( \log_2 \left( \frac{C(x)}{\text{Kn}} \right) \right) \in \Z \]
Dieser Faktor beschreibt die Anzahl der empfohlenen Auflösungsverdoppelungen.
\end{Definition}

Dieses, die theoretische mit der tatsächlich simulierten Knudsen-Kennzahl des Fluids vergleichende, Kriterium liefert auf diese Weise bis auf die Ebene einzelner Zellen auflösbare Informationen zur lokalen Simulationsqualität. Beschränken wir dessen Wertebereich zur besseren Unterscheidung auf eine diskrete Menge, erhalten wir folgende Darstellung der unverfeinerten Simulation:

\begin{figure}[H]
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_unrefined_n20_re100_16s_knudsen.pdf}
\end{adjustbox}
\caption{Verfeinerungskriterium in einem uniform mit \(N=20\) aufgelösten Gitter}
\label{fig:UnrefinedCylinderKnudsen60s}
\end{figure}

Der lokale Vergleich der Knudsen-Zahlen eröffnet einen neuen, interessanten, Blick auf die Fluidstruktur -- klar zu erkennen sind zunächst die wechselseitigen Wirbel sowie die Strömungskomplexität um den Zylinder im Eingangsbereich. Abseits dieses Hindernisses finden sich größere Bereiche mit guter Auflösung, also einem Verfeinerungsfaktor von größer oder gleich Null. Darüber hinaus weisen große Teile des Gitters mit Verfeinerungsfaktoren bis eins eine noch akzeptable Auflösung auf. Den größten Mangel an Simulationsqualität attestiert das Verfeinerungskriterium direkt um den Zylinder und in den wandnahen Wirbelbereichen.

\begin{figure}[H]
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_unrefined_n40_re100_16s_knudsen.pdf}
\end{adjustbox}
\caption{Verfeinerungskriterium in einem uniform mit \(N=40\) aufgelösten Gitter}
\label{fig:UnrefinedCylinderKnudsenN4060s}
\end{figure}

Bei der Interpretation der Verfeinerungsfaktoren ist zu beachten, dass einzelne kleine Bereiche mit großen Faktoren keine global mangelnde Auflösung beschreiben. Das auf den Faktoren aufbauende automatische Gitterverfeinerungskriterium betrachtet dazu jeweils die durchschnittliche Qualität von a priori kartierten möglichen Verfeinerungsdomänen -- nur wenn die extremalen Bereiche ausreichend aufgelöste Gebiete dominieren, ist demnach eine Verfeinerung empfohlen.

\begin{figure}[h]
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_single_refinement_n20_re100_16s_knudsen.pdf}
\end{adjustbox}
\caption{Verfeinerungskriterium in einem einfach verfeinerten \(N=20\) Gitter}
\label{fig:SingleLevelRefinementCylinderKnudsen16s}
\end{figure}

Insgesamt war unter dieser formaleren Analyse unsere intuitive Wahl des in Abbildung~\ref{fig:SingleLevelRefinementCylinderVelocity16s} verfeinerten Bereichs akzeptabel. Dies bestätigt sich auch bei Berechnung des Knudsen-Kriterium für das einfach verfeinerte Gitter in Abbildung~\ref{fig:SingleLevelRefinementCylinderKnudsen16s} -- die angemahnten Bereiche im Umfeld des Zylinders sind hier deutlich reduziert. Gehen wir weiterhin davon aus, dass das Kriterium belastbare Aussagen zur Simulationsqualität liefert, so ist die Reduzierung des des Verfeinerungsfaktors in den verfeinerten Gebieten um genau eins ein Indiz für die Qualität des Verfeinerungsalgorithmus.

\begin{figure}[h]
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_improved_grid_n20_re100_16s_knudsen.pdf}
\end{adjustbox}
\caption{Verbesserte Verfeinerungsstruktur um den Zylinder}
\end{figure}

Typischerweise sind LBM-Simulationen des umströmten Zylinders für wirbelbildende Reynolds-Zahlen bei kleineren Auflösungen empfindlich gegenüber Divergenz im Ausflussbereich. Dies deutet sich auch bei noch ausreichender Auslösung in den Verfeinerungsfaktoren des Ausflussbereiches an. So divergiert ein uniform mit \(N=5\) aufgelöstes Gitter mit den verwendeten Parametern und Randbedingungen schon nach wenigen Schritten.
Es wäre vorteilhaft, wenn sich dieses Problem unter Einsatz möglichst weniger zusätzlicher Gitterknoten beheben ließe. Und tatsächlich genügt schon die Verfeinerung eines schmalen Stücks des Ausgangsbereiches zur Stabilisierung der Simulation bei dieser niedrigen Auflösung. In Abbildung~\ref{fig:CylinderOptimizedGridN5} sehen wir den Geschwindigkeitsplot eines solchen Gitters.

\bigskip

Dieses Beispiel gewinnt insbesondere an Eindruck, wenn wir die Gesamtanzahl der Gitterzellen mit den für ein stabiles uniformes Gitter benötigten Zellen vergleichen: Während das im Ausgang verfeinerte Gitter nur 3608 aktive Zellen beinhaltet, benötigt ein stabiles uniformes Gitter eine Auflösung von mindestens \(N=12\), was 13500 Zellen entspricht. Interessieren wir uns an dieser Stelle also nur für ein subjektiv korrektes und nicht-divergentes Strömungsbild, so ermöglicht uns Gitterverfeinerung ein Auskommen mit nur etwa einem Viertel der für ein unverfeinertes Gitter benötigten Zellen.

\begin{figure}[h]
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_low_resolution_outflow_refine_n5_re100_16s.pdf}
\end{adjustbox}
\caption{Stabile Zylinderumströmung mit \(N=5\) durch Ausgangsverfeinerung}
\label{fig:CylinderOptimizedGridN5}
\end{figure}

Alternativ ist es möglich über das Festhalten der Anzahl der Freiheitsgrade, d. h. der Knotenanzahl, beschränkte Rechenressourcen besser auf die zu simulierende Strömungssituation aufzuteilen. Beispielsweise fällt auf, dass der zentrale Zylinder in Abbildung~\ref{fig:CylinderOptimizedGridN5} nur sehr grob modelliert wird. Entsprechend wollen wir versuchen, bei Einschränkung auf die Anzahl der Gitterknoten eines uniformen Gitters mit \(N=12\), ein besseres Gitter zu erzeugen.

\begin{figure}[H]
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_unrefined_n12_re100_16s.pdf}
\end{adjustbox}
\begin{adjustbox}{center}
\includegraphics[width=1.2\textwidth]{img/static/cylinder2d_optimized_refinement_n5_re100_16s.pdf}
\end{adjustbox}
\caption{Vergleich von uniform und problembezogen verteilten Freiheitsgraden}
\label{fig:CylinderOptimizedGridComparison}
\end{figure}

\newpage
Wir sehen das Geschwindigkeitsbild dieser Bemühungen in Abbildung~\ref{fig:CylinderOptimizedGridComparison}. Die dort dargestellten Gitter beinhalten beide jeweils maximal 13500 Zellen. Der kleine Unterschied in der Knotenanzahl ist dabei der Einschränkung auf quaderförmige Gitter geschuldet, welche eine exakte Fixierung der Knotenanzahl erschwert.

Klar zu erkennen ist die in der verfeinerten Variante deutlich bessere Diskretisierung des Zylinders durch Konzentration der verfügbaren Gitterknoten in dessen Umfeld. Auch liegt dem Ausfluss des verfeinerten Gitters die Divergenz ferner als dem Ausfluss des uniformen Gitters, an welchem sich schon Artefakte abzeichnen. Eine formalere Analyse der Qualität dieses optimierten Gitters erwartet uns im nachfolgenden Unterkapitel.

\bigskip

Bevor wir dazu kommen, bemerken wir, dass sich mit dieser flexibleren Verteilung der Knotenfreiheitsgrade hier auch ein, bis jetzt lediglich in der Einführung erwähnter, Vorteil von Gitterverfeinerung illustriert: Selbst wenn ein Verfeinerungsverfahren bezogen auf den Fehler im Vergleich mit analytischen Lösungen keine Verbesserungen oder sogar leichte Einbußen produziert, kann es doch potenziell eingesetzt werden, um Probleme zu behandeln, welche anderweitig nicht oder nur mit deutlich höherem Speicher- und Rechenaufwand zugänglich wären.

\subsubsection{Vergleich von Widerstands- und Auftriebskoeffizienten}

Bis hier haben wir die den Einfluss von Gitterverfeinerung auf die Zylinderumströmung entweder subjektiv, durch Vergleich der Strömungsbilder, oder grob, durch einfachen Vergleich der für eine divergenzfreie Simulation benötigten Knotenanzahl, bewertet. Wie eingangs erwähnt, existiert für die vorliegende Strömungssituation keine analytische Lösung, weshalb wir uns für eine formal belastbarere Bewertung auf vertrauenswürdige aber ebenfalls simulierte Referenzwerte stützen wollen.

\citetitle{SchaeferTurek96}~\cite{SchaeferTurek96} liefert eine, dieser formalen Bewertung dienliche, Übersicht der, von 17 verschiedenen Forschungsgruppen beigetragenen und auf unterschiedliche Weisen gewonnenen, Lösungen einer klar definierten Zylinderumströmung. Diese Lösungen sind dabei anhand einer Auswahl von charakteristischen Messwerten wie dem Strömungswiderstands- und Auftriebskoeffizient sowie dem Druckunterschied zwischen Vorder- und Rückseite des Zylinders gegeben.

\begin{Definition}[Strömungswiderstands- und Auftriebskraft]
	Beschreite der Weg \(S\) den Rand des Zylinders und sei \(n \in \R^2\) dessen Normale, \(v_t\) die Geschwindigkeit entlang der Tangente \(t:=(n_1,-n_0)\), \(\rho_c\) die charakteristische Dichte der Flüssigkeit, \(P\) der Druck sowie \(\nu\) die kinematische Viskosität. Dann sind Widerstands- und Auftriebskraft des Zylinders gegeben als~\cite[Kap.~2.2]{SchaeferTurek96}:
\begin{align*}
F_w &= \int_S \left( \rho_c \nu \frac{\partial v_t}{\partial n} n_1 - P n_0 \right) dS && \text{Widerstandskraft}\\
F_a &= - \int_S \left( \rho_c \nu \frac{\partial v_t}{\partial n} n_0 + P n_1 \right) dS && \text{Auftriebskraft}
\end{align*}
\end{Definition}

\newpage
\begin{Definition}[Strömungswiderstands- und Auftriebskoeffizient]
Sei \(D\) der Zylinderdurchmesser, \(\overline{U}\) die durchschnittliche Fluidgeschwindigkeit und \(\rho_c\) die charakteristische Dichte. Dann sind die Widerstands- und Auftriebskoeffizienten gegeben als~\cite[Kap.~2.2]{SchaeferTurek96}:
\begin{align*}
c_w &= \frac{2F_w}{\rho_c \overline{U}^2 D} && \text{Widerstandskoeffizient} \\
c_a &= \frac{2F_a}{\rho_c \overline{U}^2 D} && \text{Auftriebskoeffizient}
\end{align*}
\end{Definition}

Vergleichen wollen wir diese Koeffizienten nun im Rahmen des unstetigen Testfalls \cite[2.2b]{SchaeferTurek96} mit Reynolds-Zahl \(\text{Re}=100\), maximaler Einflussgeschwindigkeit \(U = 1.5 \,\text{m}/\text{s}\) und charakteristischer Dichte \(\rho_c = 1.0 \,\text{kg}/\text{m}^3\). Dieses Ziel hatten wir dabei schon zu Beginn dieses Kapitels im Blick, so dass die OpenLB-basierende Modellierung der Zylinderumströmung inklusive des optimierten Gitters in Abbildung~\ref{fig:CylinderOptimizedGridComparison} direkt verwendet werden kann. Für die Berechnung der Koeffizienten stellt OpenLB dabei in Form des \class{SuperLatticePhysDrag2D} Funktors bereits ein passendes \emph{Messinstrument} bereit.

\bigskip

Die folgenden Abbildungen zeichnen demnach den zeitlichen Verlauf der durch das uniform mit \(N=12\) aufgelöste Gitter berechneten Charakteristiken zusammen mit den Resultaten des problembewusst verfeinerten \(N=5\) Gitters mit näherungsweise gleicher Anzahl Freiheitsgraden. Der anzustrebende formale Referenzwert ist dabei jeweils das Mittel der oberen und unteren Grenzwerte \cite[Tabelle~4]{SchaeferTurek96}.

\begin{figure}[H]
\centering
\input{img/cylinder2d_drag_comparison.tikz}
\caption{Zeitlicher Verlauf des Strömungswiderstandskoeffizients}
\end{figure}

\begin{figure}[H]
\centering
\input{img/cylinder2d_lift_comparison.tikz}
\caption{Zeitlicher Verlauf des Auftriebskoeffizients}
\end{figure}

\begin{figure}[H]
\centering
\input{img/cylinder2d_deltap_comparison.tikz}
\caption{Zeitlicher Verlauf der Druckdifferenz}
\end{figure}

\newpage
Klar zu erkennen ist, dass die Simulation auf dem verfeinerten Gitter über alle drei Kriterien sowohl die Referenzwerte besser trifft als auch unabhängig davon eine bessere Qualität bezogen auf Fluktuation und Konsistenz in der periodischen Strömungsphase aufweist. Da diese Ergebnisse bei näherungsweise beibehaltener Anzahl von Gitterknoten erzielt wurden -- tatsächlich verwendet das verfeinerte Gitter sogar 46 Knoten weniger als das uniform aufgelöste -- ist somit auch bei formaler Betrachtung der positive Nutzen von Gitterverfeinerung zur Simulation der Zylinderumströmung klar demonstriert. Wir erhalten bei leicht geringerer Anzahl von Knoten und damit näherungsweise gleichem Speicherbedarf ein deutlich besseres Simulationsergebnis.

\bigskip

Das diesen Resultaten zugrundeliegende Gitter orientiert sich in der Anzahl seiner Zellen an einem uniform aufgelösten Gitter, welches knapp an der Grenze zur Divergenz im Ausflussgebiet arbeitet. Die im Vergleich dazu sehr guten Ergebnisse des optimierten Gitters sind also nicht zu überraschend.
Als weiteren Vergleich betrachten wir deshalb die entscheidenden maximalen Widerstands- und Auftriebskoeffizienten \(\widehat{c_w}\) bzw. \(\widehat{c_a}\) eines uniform mit \(N=40\), d. h. der maximalen lokalen Auflösung unseres problembewusst verfeinerten Gitters entsprechend, aufgelösten Gitters.
Auch im Vergleich mit diesem 145314 Knoten umfassenden Gitter bewähren sich die, mit nur etwa einem Zehntel der Knoten gewonnenen, Ergebnisse der verfeinerten Simulation:

\begin{table}[h]
\centering
\sisetup{round-precision=4}
\begin{tabular}{l l l l l}
& Uniform & Verfeinert & Referenzintervall \cite{SchaeferTurek96} \\
\hline
\hline
\(\widehat{c_w}\) & \num{3.28433} & \num{3.21927} & \([\num[round-mode=off]{3.22},\ \num[round-mode=off]{3.24}]\) \\
\(|\widehat{c_w}-\num[round-mode=off]{3.23}|\) & \num{0.05433} & \num{0.01073} & [\num[round-mode=off]{0}, \num[round-mode=off]{0.01}] \\
\hline
\(\widehat{c_a}\) & \num{1.07046} & \num{1.10359} & \([\num[round-mode=off]{0.99},\ \num[round-mode=off]{1.01}]\) \\
\(|\widehat{c_a}-\num[round-mode=off]{1.0}|\) & \num{0.07046} & \num{0.10359} & [\num[round-mode=off]{0}, \num[round-mode=off]{0.01}]  \\
\hline
\(\Delta P\) & \num{2.5793} & \num{2.44285} & \([\num[round-mode=off]{2.46},\ \num[round-mode=off]{2.5}]\)  \\
\(|\Delta P-\num[round-mode=off]{2.48}|\) & \num{0.0993} & \num{0.03715} & [\num[round-mode=off]{0}, \num[round-mode=off]{0.02}] \\
\hline
\hline
\(N=\) & \num{40} & \num{5} (max: \num{40}) & \\
Knotenanzahl & \num{145314} & \num{13454} & \\
\end{tabular}
\caption{Vergleichswerte zweier maximal \(N=40\) aufgelöster Gitter}
\label{tab:cylinder2dComparison}
\end{table}

\noindent
Tatsächlich ist der Fehler des verfeinerten Gitters für Widerstandskoeffizient und Druckdifferenz sogar kleiner als der des uniformen Gitters mit ungleich mehr Knoten.

\bigskip

Wir haben an dieser Stelle also auch im formalen Vergleich bestätigt, dass sich Gitterverfeinerung zur besseren Verteilung beschränkter Rechenressourcen einsetzen lässt.
Die bestimmten Vergleichswerte bestehen bei geeigneter Variation der lokalen Gitterweiten auch in Konkurrenz mit uniformen Gittern, die auf der ganzen Simulationsdomäne der feinsten Gitterweite des heterogenen Gitters entsprechend aufgelöst sind. Es stellt sich daher die Frage, ob dieser klare Vorteil auch auf höher aufgelöste Gitter übertragen werden kann und sich die Ergebnisse in vergleichbarem Maße verbessern.
Dazu sehen wir in Abbildung~\ref{fig:cylinder2dHighResComparison} und zugehöriger Tabelle~\ref{tab:cylinder2dHighResComparison} die aerodynamischen Kennzahlen der uniformen \(N=48\) und \(N=80\) Gitter sowie eines problembezogen varierten \(N=20\) Gitters entsprechend der Struktur in Abbildung~\ref{fig:CylinderOptimizedGridComparison} ohne Ausflussverfeinerung.

\begin{figure}[H]
\centering
\input{img/cylinder2d_high_res_comparsion.tikz}
\caption{Zeitlicher Verlauf der Messwerte höher aufgelöster Gitter}
\label{fig:cylinder2dHighResComparison}
\end{figure}

\begin{table}[h]
\centering
\sisetup{round-precision=4}

\begin{tabular}{l l l l l l l}
& Uniform & Uniform & Uniform & Verfeinert & Referenzintervall \cite{SchaeferTurek96} \\
\hline
\hline
\(\widehat{c_w}\) & \num{3.23526} & \num{3.24441} & \num{3.24897} & \num{3.19928} & \([\num[round-mode=off]{3.22},\ \num[round-mode=off]{3.24}]\) \\
\(|\widehat{c_w}-\num[round-mode=off]{3.23}|\) & \num{0.00526} & \num{0.01441} & \num{0.01897} & \num{0.03072} & [\num[round-mode=off]{0}, \num[round-mode=off]{0.01}] \\
\hline
\(\widehat{c_a}\) & \num{1.00194} & \num{1.01377} & \num{1.02121} & \num{0.995577} & \([\num[round-mode=off]{0.99},\ \num[round-mode=off]{1.01}]\) \\
\(|\widehat{c_a}-\num[round-mode=off]{1.0}|\) & \num{0.00194} & \num{0.01377} & \num{0.02121} & \num{0.004423} & [\num[round-mode=off]{0}, \num[round-mode=off]{0.01}]  \\
\hline
\(\Delta P\) & \num{2.49969} & \num{2.50671} & \num{2.44354} & \num{2.46411} & \([\num[round-mode=off]{2.46},\ \num[round-mode=off]{2.5}]\)  \\
\(|\Delta P-\num[round-mode=off]{2.48}|\) & \num{0.01969} & \num{0.02671} & \num{0.03646} & \num{0.01589} & [\num[round-mode=off]{0}, \num[round-mode=off]{0.02}] \\
\hline
\hline
\(N=\) & \num{160} & \num{80} & \num{48} & \num{20} (max: \num{160}) & \\
Knotenanzahl & \num{2298014} & \num{576758} & \num{207862} & \num{208031} & \\
\end{tabular}
\caption{Aerodynamische Kennzahlen höher aufgelöster Zylinderumströmungen}
\label{tab:cylinder2dHighResComparison}
\end{table}

Das geeignet verfeinerte Gitter liefert demnach in zwei von drei Messwerten einen kleineren Fehler als ein äquivalent aufgelöstes \(N=48\) Gitter und besteht im Vergleich mit dem uniformen \(N=80\) Gitter erneut gegenüber einer deutlich mehr Knoten nutzenden Simulation. Selbst gegenüber eines mehr als zehnmal so viele Knoten aufwendenden uniformen \(N=160\) Gitter erhalten wir mit der uniformen Variante noch einen leicht kleineren Fehler in der Druckdifferenz.

Da alle vier getesteten Gitter gute Übereinstimmung zu den Referenzdaten von Schäfer und Turek aufweisen, fällt der Vorteil der Gitterverfeinerung im Allgemeinen jedoch geringer aus, als noch im Vergleich der niedrig aufgelösten Gitter aus Abbildung~\ref{fig:CylinderOptimizedGridComparison} und Tabelle~\ref{tab:cylinder2dComparison}. Darüber hinaus liegt die Knotenanzahl des betrachteten \(N=160\) Gitters weit oberhalb der maximalen referenzstiftenden Knotenzahlen \cite[Tabelle~4]{SchaeferTurek96}, so dass wir hier an die Grenzen der Aussagekraft von Fehlern bezüglich dieser Werte stoßen.

\bigskip
Rückblickend hat die Evaluation unseres Gitterverfeinerungsverfahrens anhand der Zylinderumströmung dessen Vorteil spendende Verwendbarkeit eindrücklich demonstriert. Während bei der Betrachtung der einfachen Poiseuilleströmung Sorgen bezüglich der Ergebnisgenauigkeit verfeinerter Gitter geweckt wurden, konnten diese in der Betrachtung einer komplexeren und so einen Verfeinerungsbedarf besser begründenden Strömungssituation weitestgehend beigelegt werden. Während die aufgeworfene Frage nach der korrekten Behandlung von Randbedingungen im Übergangsbereich noch offen ist, konnten wir unter Vermeidung dieser unklaren Situation ein sinnvolles, durch ein formales Kriterium informiertes, lokal verfeinertes Gitter beschreiben.

Verfeinerte Simulationen der Zylinderumströmung konnten auf diese Weise im Vergleich ihrer aerodynamischen Kennzahlen mit belastungsfähigen Referenzwerten \cite{SchaeferTurek96} durchweg gegen ungleich größere Knotengrade aufweisende uniforme Gitter bestehen.

\newpage
\section{Fazit}

Zusammenfassend hinterlässt die Anwendung Gitterverfeinerung in Lattice Boltzmann Methoden einen vielversprechenden Eindruck mit klaren Fragestellungen für mögliche Wege des Aufbaus auf dieser Arbeit. So hat sich die theoretische Formulierung des Verfahrens von Lagrava et al. als praktisch gut umsetzbar erwiesen und die zugehörigen Ergebnisse in \citetitle{Lagrava12}~\cite{Lagrava12} konnten reproduziert werden.

Die existierenden Strukturen der LBM-Bibliothek OpenLB ließen sich gut mit einer lokalen Verfeinerung des Gitters vereinbaren und der gewählte Implementierungsansatz, der Betrachtung von verfeinernden Gittern als weitestgehend eigenständige Simulationen mit besonderen Randkonditionen in den Übergangsbereichen, bewährte sich sowohl in der initialen Umsetzung des Verfahrens als auch in der Verwendung der resultierenden Schnittstelle zur Entwicklung der Evaluationsbeispiele.

Die im Rahmen dieser Beispiele betrachtete Zylinderumströmung bestätigte das Verfahren im Vergleich von verfeinerten und uniformen Simulationen als gewinnbringend einsetzbar -- aerodynamische Kennzahlen mit guter Übereinstimmung zu hochwertigen Vergleichsdaten \cite{SchaeferTurek96} konnten durch gezielte Verfeinerung mit einer deutlich reduzierten Anzahl von Gitterknoten simuliert werden. Dies ist insbesondere auch im Hinblick darauf beachtlich, dass ein mit äquivalenter Anzahl von Gitterknoten uniform aufgelöstes Gitter signifikant schlechtere und zur Divergenz neigende Ergebnisse lieferte. In diesem Kontext offenbarte sich der Einsatz von Gitterverfeinerung zur gezielten Vermeidung von Divergenz im Ausflussbereich als zusätzliches Anwendungsgebiet.

Problematischer erwies sich das Verfahren in der Anwendung auf eine einfache und analytisch lösbare Poiseuilleströmung. Im Verlauf der Untersuchung dieser grundlegenden Strömungssituation kristallisierte sich das in \cite{Lagrava12} unbeachtete Zusammenspiel von Randkonditionen und Übergangsbereichen durch Vergrößerung des Fehlers um bis zu eine Größenordnung als kritischer und weiterer theoretischer Untersuchung bedürfender Aspekt heraus. Auch bei Vermeidung von Randbedingungen in Auflösungsübergängen konnte ein negativer Einfluss von Gitterverfeinerung auf die Reproduktion der analytischen Lösung festgestellt werden. Es wurde somit klar, dass Gitterverfeinerung in dem hier untersuchten Rahmen keinesfalls unvorsichtig und in der Abwesenheit konkreter Notwendigkeiten angewendet werden sollte. Die zusätzliche Komplexität und Fehlerquelle einer lokalen Verfeinerung sollte in sinnvollen Anwendungen eben dieser durch Vorteile wie bessere Geometriediskretisierung oder Zwänge wie beschränkte Rechenressourcen aufgewogen oder geringstenfalls begründet werden.

Im Kontext der sinnvollen Anwendung von Gitterverfeinerung sowie der konkreten Strukturierung des heterogenen Gitters erwies sich das in \citetitle{Lagrava15}~\cite{Lagrava15} entwickelte Gitterverfeinerungskriterium als sehr gutes Maß zur Bewertung der lokalen Simulationsqualität. Nicht nur konnten in ihrer Anzahl beschränkte Knotenfreiheitsgrade mittels dieses Kriteriums formal fundiert problembezogen umverteilt werden, sondern auch problematische und zur Divergenz neigende Bereiche der Simulation ließen sich frühzeitig erkennen und vermeiden.

\bigskip
\noindent
Abschließend ergeben sich somit die folgenden theoretischen Fragestellungen zur weiteren Verfolgung in absteigender Dringlichkeit:
\begin{enumerate}
	\item Wie müssen Randkonditionen in Übergangsbereichen behandelt werden?
	\item Wie kann die Rechenlast zur parallelen Simulation der verfeinernden Gitter besser verteilt werden?
	\item Wie verhält sich das betrachtete Verfahren bei Übertragung auf dreidimensionale Lattice Boltzmann Methoden?
	\item Was gilt es bei Variation der Verfeinerung im Verlauf der Simulation zu beachten?
\end{enumerate}

\noindent
Bezüglich der Weiterentwicklung des nun in OpenLB existierenden Gitterverfeinerungsframework bieten sich darüber hinaus Ansatzpunkte in Form der Fragen:
\begin{enumerate}
	\item Wie kann der Kommunikationsaufwand zur Gitterkopplung reduziert werden?
	\item Welche Möglichkeiten gibt es zur weiteren Abstraktion des Simulationsaufbaus?
	\item Wie lassen sich die existierenden Konzepte aus der Optimierung und Mehrphasenkopplung sinnvoll mit den neuen Möglichkeiten zur Gitterverfeinerung verbinden?
	\item Kann die Verfeinerungsschnittstelle unabhängig des konkreten Verfahrens direkt auf drei Dimensionen übertragen werden?
\end{enumerate}
Diese Auswahl von Möglichkeiten für anknüpfende Arbeiten beschließe an dieser Stelle unsere Betrachtung von gitterverfeinerten Lattice Boltzmann Methoden in OpenLB. Die angestrebten Ziele der theoretischen Ausarbeitung, flexiblen praktischen Umsetzung und experimentellen Evaluation wurden erreicht und Gitterverfeinerung verspricht den Funktionskatalog der freien Lattice-Boltzmann-Bibliothek OpenLB um ein nützliches und ausbaufähiges Werkzeug erweitert zu haben.