source: palm/trunk/DOC/tec/methods/upstream_spline/ups.tex @ 1485

Last change on this file since 1485 was 481, checked in by raasch, 15 years ago

file structure of technical documentation revised; documents concerning numerical methods added

File size: 21.1 KB
Line 
1\documentclass[11pt,a4paper,fleqn,abstracton]{scrartcl}
2\usepackage[T1]{fontenc}
3\usepackage{enumerate,sqrt,chimuk,graphicx}
4\usepackage[german]{babel}
5\usepackage[hang,bf]{caption}
6
7\bibliographystyle{ups}
8
9\def\version{1.0}
10
11\title{Das Upstream-Spline Advektionsverfahren}
12\author{Michael Schr"oter}
13\date{Version \version{}, \today}
14
15
16\begin{document}
17%
18  \maketitle
19%
20  \nocite{mahrer:78,price:73}
21%
22\begin{abstract}
23     Dokumentation des Upstream-Spline Advektionsverfahrens. Es behandelt die
24     allgemeine Berechnung des Splines und die Berechnung im speziellen bei
25     Verwendung bestimmter Randbedingungen. Zus"atzlich wird die L"osung
26     tridiagonaler linearer Gleichungssysteme vorgestellt. Abschlie"send
27     wird ein selektiver Filter vorgestellt, der es erm"oglicht sogenannte
28     $2\Delta$-Instabilit"aten, welche zum Beispiel bei der Upstream-Spline
29     Advektion k"unstlich angeregt werden, auszufiltern.
30\end{abstract}
31%
32Gesucht ist die L"osung der Advektionsterme in prognostischen Gleichungen, welche hier
33f"ur eine allgemeine Gr"o"se $\psi$ und die Raumrichtung $x$ formuliert:
34\begin{equation}\label{advek:eq}
35\frac{\6\psi}{\6 t}=-u\frac{\6\psi}{\6 x}\;
36\end{equation}
37lauten.
38Ist die advehierende Geschwindigkeit $u$ konstant, so entspricht die L"osung von
39(\ref{advek:eq}) der Verlagerung von $\psi$ um die Distanz $u\cdot \Delta t$ 
40innerhalb des Zeitintervalls $\Delta t$. Die numerische L"osung von
41(\ref{advek:eq}) setzt eine konstante Verlagerungsgeschwindigkeit  w"ahrend eines
42Zeitschritts $\Delta t$ voraus. Entspricht die Verlagerungsdistanz
43$u\cdot \Delta t$ nicht dem ganzzahligen Vielfachen der numerischen Gitterweite,
44so ist die Anwendung eines Interpolationsverfahrens notwendig, um Werte von $\psi$ 
45zwischen den Maschenpunkten zu erhalten.
46
47Der Unterschied zu anderen Advektionsverfahren, wie zum Beispiel dem
48\textit{Upstream}-Verfahren oder dem Verfahren von \textit{Piacsek und Williams},
49bei denen zur Beschreibung der Advektion eine Approximation von
50Differenzenquotienten vorgenommen wird, liegt darin, da"s bei dem
51\textit{Upstream-Spline Advektionsverfahren} die Str"omung durch eine
52interpolierende Funktion angen"ahert wird. Als interpolierende Funktion setzt das
53\textit{Upstream-Spline Advektionsverfahren} eine
54kubische Splinefunktion $s_{\Delta}$ an. Es handelt sich dabei um eine Funktion,
55die auf jedem Teilintervall $[x_{i},x_{i+1}]$, mit $i=0,\ldots,N-1$, zweimal stetig
56differenzierbar sind und dort mit einem Polynom dritten Grades "ubereinstimmen.
57Eine Splinefunktion ist somit st"uckweise aus $N$ kubischen Polynomen so zusammengesetzt,
58da"s die Funktion $s_{\Delta}$ selbst und ihre beiden ersten Ableitungen an den
59Stellen $x_i$, mit $i=1,\ldots,N-1$, keine Sprungstellen besitzen
60\comshortcite{stoer:72}. D.~h. an den Randpunkten der Teilintervalle werden die
61Splinefunktionen $s_{\Delta,i}$ durch folgende  Beziehungen miteinander verkn"upft:
62\begin{eqnarray}
63s_{\Delta,i}(x_i)  & = & s_{\Delta,i+1}(x_i) \nonumber \\
64s'_{\Delta,i}(x_{i}) & = & s'_{\Delta,i+1}(x_{i})
65                              \hspace{0.5cm}\mbox{f"ur}\hspace{0.25cm} i=0,\ldots,N\!-\!1\nonumber \\
66s''_{\Delta,i}(x_{i}) & = & s''_{\Delta,i+1}(x_{i})\;. \nonumber
67\end{eqnarray}
68$s_{\Delta,i}(x_i)$ ist hier die auf dem $i$-ten Teilintervall g"ultige
69Splinefunktion und $x_i$ sind die Knotenpunkte (St"utzstellen) des Splines.
70Sie entsprechen des Knotenpunkten des numerischen Gitters in die jeweilige
71betrachtete Raumrichtung.
72
73Unter Verwendung der Steigungen der Splinefunktionen $m_i$, kann die Splinefunktion
74als
75\begin{eqnarray}\label{spline_eq}
76s_\Delta(x) & = & m_{i-1}\frac{\left(x_i-x\right)^2 \left(x-x_{i-1}\right) }{h_i^2}
77                  -m_i\frac{\left(x-x_i\right)^2\left(x_i-x\right)}{h_i^2}\nonumber\\
78            &   & +\psi_{i-1}
79                      \frac{\left(x_i-x\right)^2\left[2\left( x-x_i\right) +h_i\right]}{h_i^3}
80                  +\psi_i
81                      \frac{\left(x-x_{i-1}\right)^2\left[2\left(x_i-x\right)+h_i\right]}{h_i^3}
82\end{eqnarray}
83formuliert werden \comshortcite{ahlberg:67}, mit $h_i=x_i-x_{i-1}=dx_i$.
84F"ur die advehierte Gr"o"se an einem Gitterpunkt $i$ zum Zeitpunkt $t+\Delta t$ gilt
85dann:
86\begin{equation}
87  \psi_i^{t+\Delta t} = s^{t}_\Delta(x+\alpha_i\Delta x)\;,
88\end{equation}
89mit $\alpha_i=u_i\frac{\Delta t}{\Delta x}$. F"ur $u_i\ge 0$ erh"alt man durch
90Einsetzen von $x-\alpha_i$ in (\ref{spline_eq}):
91\begin{eqnarray}
92  \psi_i^{t+\Delta t} & = &
93     \psi_i^{t} - m_i\alpha + \left[
94       3\left(\psi_{i-1}^{t}-\psi_i^{t} \right) + 2m_i + m_{i-1}\right]\alpha^2 \\
95                      &   & - \left[ m_i + m_{i-1} +
96           2\left( \psi_{i-1}^{t}-\psi_i^{t}\right)\right]\alpha^3\;.
97\end{eqnarray}
98F"ur $u_i<0$ folgt mit $x+\alpha_i$ (der Spline wird hier zwischen den Punkten $x_i$
99und $x_i+1$ ausgewertet):
100\begin{eqnarray}
101  \psi_i^{t+\Delta t} & = &
102     \psi_i^{t} + m_i\alpha + \left[
103       3\left(\psi_{i}^{t}-\psi_{i+1}^{t} \right) + 2m_i + m_{i+1}\right]\alpha^2 \\
104                      &   & - \left[ m_i + m_{i+1} +
105           2\left( \psi_{i}^{t}-\psi_{i+1}^{t}\right)\right]\alpha^3\;.
106\end{eqnarray}
107
108Unter der Voraussetzung, da"s die Splinefunktionen an ihren R"andern stetig
109ineinander "ubergehen, ergibt sich folgende Bestimmungsgleichung f"ur die Steigungen
110$m_i$:
111\begin{equation}
112\frac{1}{h_i}m_{i-1}+2\left(\frac{1}{h_i}+\frac{1}{h_{i+1}}\right)m_i+\frac{1}{h_{i+1}}m_{i+1}=
113 3\frac{\psi_i-\psi_{i-1}}{h_i^2}+3\frac{\psi_{i+1}-\psi_i}{h_{i+1}^2}\;.
114\end{equation}
115F"ur $i$ gilt hier $i=1,\ldots,N-1$. Mit
116\begin{eqnarray}
117\lambda_i & = & \frac{h_{i+1}}{h_i+h_{i+1}}\label{lambda:eq}, \\[8pt]
118\mu_i     & = & 1-\lambda_i \label{mu:eq}, \\[8pt]
119c_i       & = & 3\frac{\psi_i-\psi_{i-1}}{h_i^2}+3\frac{\psi_{i+1}-\psi_i}{h_{i+1}^2}
120\end{eqnarray}
121folgt
122\begin{equation}\label{spl:eq}
123\lambda_i m_{i-1} + 2 m_i + \mu_i m_{i+1} = c_i\;.
124\end{equation}
125
126Die R"ander des Modellgebietes m"ussen eine gesonderte Untersuchung
127erfahren. Hierzu ist die Formulierung von Randbedingungen notwendig, wobei zwischen
128lateralen und vertikalen Randbedingungen unterschieden werden mu"s. Lateral werden in
129der LES meist zyklische Randbedingungen angesetzt. Der Index $i$ l"auft hier nicht von 0 bis $N$,
130sondern es gilt $i=-1,\ldots ,nx+1$ (f"ur $x$-Richtung). Es gilt:
131\[
132\psi_{-1} = \psi_{nx}\hspace{1cm}\mbox{und}\hspace{1cm}\psi_{0} = \psi_{nx+1}\;.
133\]
134In vertikaler Richtung mu"s hingegen zwischen Neumannschen und Dirichletschen
135Randbedingungen unterschieden werden. Der Laufindex f"ur die vertikale
136Richtung bewegt sich in den Grenzen $0\le i\le nz$ und es gilt:
137\begin{enumerate}
138\item Neumann-R"ander:
139\[
140\frac{\6 \psi_{0,nz}}{\6 x_i} = 0
141\]
142\item Dirichlet-R"ander:
143\[
144\psi_{0,nz}=\psi_0=const
145\]
146\end{enumerate} 
147Jede Wahl der Randbedingung wirkt sich anders auf die Bestimmung der
148Splinekoeffizienten $m_i$ aus, wie im folgenden erl"autert wird.
149
150\subsection*{Nichtperiodische Randbedingungen, Advektion in vertikaler Richtung}
151F"ur nichtperiodische Splinefunktionen (Neumann oder Dirichlet Randbedingungen) ist
152die Angabe zweier zus"atzlicher Randbedingungen f"ur $c$ und $i=0$ bzw.
153$i=nz$ notwendig:
154\begin{eqnarray}
1552 m_{0}+\mu_{0} m_1 = c_{0} & \mbox{f"ur} & i = 0 \nonumber\\
156\lambda_{nz} m_{nz-1} + 2 m_{nz} = c_{nz} & \mbox{f"ur} & i = nz
157\end{eqnarray}
158mit
159\[
160\mu_{0}= \lambda_{nz}=1
161\]
162und
163\begin{equation}
164c_{0}=3\frac{\psi_1-\psi_{0}}{h_1}\hspace{1cm}\mbox{bzw.}\hspace{1cm}
165c_{nz}=3\frac{\psi_{nz}-\psi_{nz-1}}{h_{nz}}\;.
166\end{equation}
167Dies gilt jedoch nur unter der Voraussetzung, da"s die Kr"ummungen der
168Ausgangsfunktion $\psi''(x)$ bzw. der Splinefunktionen $s''_\Delta(x)$ an den
169R"andern verschwinden. F"ur $i=1,\ldots,nz-1$ gilt weiterhin:
170\begin{equation}
171c_i=3\lambda_i\frac{\psi_i-\psi_{i-1}}{h_i}+3\mu_i\frac{\psi_{i+1}-\psi_i}{h_{i+1}}\;.
172\end{equation}
173Die einzige Unbekannte sind hier die Steigungen der Splinefunktionen $m_i$. Zu ihrer
174Bestimmung ist die L"osung des folgenden linearen Gleichungssystems notwendig:
175\begin{equation}
176\begin{array}{ccccccc}
1772 m_{0}               & + & \mu_{0} m_{1} &   &              & = & c_{0} \\
178\lambda_1 m_{0}       & + & 2 m_1         & + & \mu_1 m_{0}  & = & c_\\
179                      &   & \vdots        &   &              & = & \vdots \\
180\lambda_{nz} m_{nz-1} & + & 2 m_{nz}      &   &              & = & c_{nz}
181\end{array}
182\end{equation}
183bzw.
184\begin{equation}
185\left(\begin{array}{ccccc}
1862            & \mu_{0}  &                &              &  \\
187\lambda_1    & 2        & \mu_1          &              &           \\
188             & \ddots   & \ddots         & \ddots       &           \\
189             &          & \lambda_{nz-1} & 2            & \mu_{nz-1}  \\
190             &          &                & \lambda_{nz} & 2   
191\end{array}\right)
192\cdot
193\left(\begin{array}{c}
194m_{0} \\ m_1 \\ \vdots \\ m_{nz-1} \\ m_{nz}
195\end{array}\right)
196=
197\left(\begin{array}{c}
198c_{0} \\ c_1 \\ \vdots \\ c_{nz-1} \\ c_{nz}
199\end{array}\right)
200\end{equation}
201(leere Pl"atze stehen f"ur Nullen).
202Bei diesem linearen Gleichungssystem ist die Koeffizientenmatrix eine Tridiagonalmatrix mit
203konstanten Koeffizienten, die ausschlie"slich von der Struktur des verwendeten
204Modellgitters abh"angen (siehe (\ref{lambda:eq}) und (\ref{mu:eq})).
205
206
207\subsection*{Periodische Randbedingungen, Advektion in lateraler Richtung}
208
209F"ur zyklische laterale Randbedingungen (periodischer Splineansatz) kann
210(\ref{spl:eq}) f"ur $i=0$ und $i=nx$ vollst"andig formuliert werden:
211\begin{equation}\label{spl-lgs:eq}
212  \begin{array}{ccccccc}
213    \lambda_0 m_{-1} & + & 2 m_& + & \mu_0 m_{1} & = & c_0 \\
214    \lambda_1 m_{0}  & + & 2 m_& + & \mu_1 m_{2} & = & c_\\
215                     &   & \vdots &   &             & = & \vdots \\
216    \lambda_{nx} m_{nx-1} & + & 2 m_{nx} & + & \mu_{nx} m_{nx+1} & = & c_{nx}
217  \end{array}
218\end{equation}
219Unter Verwendung der Bedingungen f"ur den zyklischen Rand
220\newlength{\savearracs}
221\setlength{\savearracs}{\arraycolsep}
222\setlength{\arraycolsep}{8pt}
223\[
224\begin{array}{cccc}
225s_\Delta(x_{-1}) = s_\Delta(x_{nx})\,, &
226m_{-1} = m_{nx}\,, &
227s_\Delta(x_{nx+1}) = s_\Delta(x_{0})\,, &
228m_{nx+1} = m_{0}
229\end{array}
230\]
231\setlength{\arraycolsep}{\savearracs}
232kann (\ref{spl-lgs:eq}) zu
233\begin{equation}
234\begin{array}{ccccccc}
235\lambda_0 m_{nx} & + & 2 m_& + & \mu_0 m_{1} & = & c_0 \\
236\lambda_1 m_{0}  & + & 2 m_& + & \mu_1 m_{2} & = & c_\\
237                 &   & \vdots &   &             & = & \vdots \\
238\lambda_{nx} m_{nx-1} & + & 2 m_{nx} & + & \mu_{nx} m_{0} & = & c_{nx}
239\end{array}
240\end{equation}
241umgeschrieben werden, was dem L"osen von
242\begin{equation}
243\left(\begin{array}{ccccc}
2442            & \mu_&              &              & \lambda_0 \\
245\lambda_1    & 2      & \mu_1        &              &           \\
246             & \ddots & \ddots       & \ddots       &           \\
247             &        & \lambda_{nx-1} & 2            & \mu_{nx-1}  \\
248\mu_{nx}     &        &              & \lambda_{nx} & 2   
249\end{array}\right)
250\cdot
251\left(\begin{array}{c}
252m_0 \\ m_1 \\ \vdots \\ m_{nx-1} \\ m_{nx}
253\end{array}\right)
254=
255\left(\begin{array}{c}
256c_0 \\ c_1 \\ \vdots \\ c_{nx-1} \\ c_{nx}
257\end{array}\right)
258\end{equation}
259"aquivalent ist. Wie auch f"ur den nichtperiodischen Fall sind hier die Koeffizienten
260der Matrix nur von der Struktur des numerischen Gitters abh"angig, jedoch ist
261die Koeffizientenmatrix keine Tridiagonalmatrix mehr, wodurch sich der Aufwand f"ur
262die L"osung des Gleichungssystems erh"oht.
263
264\subsection*{L"osen der linearen Gleichungssysteme}
265
266Sowohl im periodischen als auch im nichtperiodischen Fall kann die Koeffizientenmatrix
267vor jeder Simulation einmal zur Verf"ugung gestellt werden und f"ur den Rest
268der Simulation unver"andert bleiben. In \textsf{PALM-1} geschieht die
269Bereitstellung der Koeffizientenmatrix nach der "Uberpr"ufung der
270Eingabeparameter (\texttt{check\_parameters.f90}) in dem Unterprogramm
271\texttt{init\_advec.f90}
272
273
274\minisec{Die Koeffizientenmatrix im nichtperiodischen Fall}
275
276F"ur tridiagonale Gleichungssysteme verringert sich der Rechenaufwand f"ur das
277Standard-Gau"s-L"osungsverfahren (Thomasalgorithmus). Ausgangspunkt f"ur diesen Algorithmus
278ist das folgende lineare Gleichungssystem:
279\begin{equation}\label{lgs:eq}
280\mathsf{A}_{(nz+1)\times (nz+1)}\cdot \vc{m} = \vc{c}\;.
281\end{equation}
282Dem Standard-Gau"s-L"osungsverfahren folgend wird zun"achst die Matrix
283$\mathsf{A}_{(nz+1)\times (nz+1)}$ in eine linke untere Dreiecksmatrix
284$\mathsf{L}_{(nz+1)\times (nz+1)}$
285und eine rechte obere Dreiecksmatrix $\mathsf{R}_{(nz+1)\times (nz+1)}$ zerlegt.
286Da $\mathsf{A}_{(nz+1)\times (nz+1)}$ eine Tridiagonalmatrix ist und alle Elemente
287au"ser den Diagonal-, Subdiagonal- und Superdiagonalelementen gleich 0 sind, sind
288bei den Matrizen $\mathsf{L}$ und $\mathsf{R}$
289auch nur die Diagonal- und Subdiagonal bzw. die Diagonal- und Superdiagonalelemente
290von 0 verschieden. Die Matrizen $\mathsf{L}$ und $\mathsf{R}$ nehmen deshalb folgende
291einfache Gestalt an:
292\[
293\mathsf{L}_{(nz+1)\times (nz+1)} =\left(\begin{array}{lllll}
2941             & 0          &        &           &       \\
295l_{1}         & 1          & \ddots &           &  \\
296              & l_{2}      & 1      & \ddots    &  \\
297              &            & \ddots & \ddots    & 0 \\
298              &            &        & l_{nz}    & 1
299\end{array}\right)
300\]
301und
302\[
303\mathsf{R}_{(nz+1)\times (nz+1)} =\left(\begin{array}{lllll}
304d_{0}  & r_{0}     &         &             & \\
3050      & d_{1}     & r_{1}   &             &  \\
306       & \ddots    & d_{2}   & \ddots      &  \\
307       &           & \ddots  & \ddots      & r_{nz-1} \\
308       &           &         & 0           & d_{nz}
309\end{array}\right)\;.
310\]
311Die Zerlegung der Matrix kann, wie bereits angesprochen, einmalig vor
312Beginn der eigentlichen Simulationen vorgenommen werden (\texttt{init\_advec.f90}),
313wobei sich die Koeffizienten wie folgt berechnen:
314\begin{eqnarray}
315d_{0}       & = & a_{0,0}       \\[6pt]
316r_{i}       & = & a_{i,i+1}\hspace{1cm}\mbox{mit}\hspace{10pt}i=0,\ldots,nz-1 \\[6pt]
317l_{i}       & = & \frac{c_{i}}{d_{i-1}}
318                           \hspace{1cm}\mbox{mit}\hspace{10pt}i=1,\ldots,nz \\[6pt]
319d_{i}       & = & a_{i,i}-a_{i-1,i}\cdot l_{i}
320                           \hspace{1cm}\mbox{mit}\hspace{10pt}i=1,\ldots,nz
321\end{eqnarray}
322Gel"ost wird das Gleichungssystem (\ref{lgs:eq}) nun in zwei Schritten, denn es gilt:
323\[
324\mathsf{A}_{(nz+1)\times (nz+1)}\cdot\vc{m} =
325\mathsf{L}_{(nz+1)\times (nz+1)}\underbrace{\mathsf{R}_{(nz+1)\times (nz+1)}
326 \cdot\vc{m}}_{\dis\glr \vc{y}} = \vc{c}
327\]
328\begin{enumerate}[\textbf{Schritt} 1:]
329\item 
330\[
331\mathsf{L}_{(nz+1)\times (nz+1)}\cdot\vc{y}=\vc{c}
332\]
333\item
334\[
335\mathsf{R}_{(nz+1)\times (nz+1)}\cdot\vc{m}=\vc{y}\;.
336\]
337\end{enumerate}
338
339\minisec{Die Koeffizientenmatrix im periodischen Fall}
340
341Im periodischen Fall handelt es sich bei der Matrix $\mathsf{A}_{(nx+1)\times (nx+1)}$
342um eine sogenannte zyklische Tridiagonalmatrix der Form:
343\begin{equation}
344\mathsf{A}_{(nx+1)\times (nx+1)}=\left(\begin{array}{ccccccc}
345a_{0,0}   & a_{0,1} & 0       & 0      & \cdots      & \beta \\
346a_{1,0}   & a_{1,1} & a_{1,2} & 0      & \cdots      & 0           \\
3470         & \ddots  & \ddots  & \ddots &             & \vdots      \\
348\vdots    & \ddots  & \ddots  & \ddots & \ddots      & 0      \\
3490         &         & \ddots  & \ddots & a_{nx-1,nx-1} & a_{nx-1,nx}    \\
350\alpha    & 0       & \cdots  & 0      & a_{nx,nx-1}   & a_{nx,nx}
351\end{array}\right)\;,
352\end{equation}
353mit $\alpha=\mu_{nx}$ und $\beta=\lambda_0$.
354
355Gleichungssysteme dieser Form k"onnen mit Hilfe der \textsc{Sherman-Morrison} Formel
356\comshortcite{press:86} gel"ost werden, dabei gilt es zwei lineare Gleichungssysteme
357zu l"osen:
358\begin{eqnarray}
359\mathsf{A}_{(nx+1)\times (nx+1)} \cdot \vc{y} & = & \vc{b}\label{lgs-cyc1:eq}\;, \\
360\mathsf{A}_{(nx+1)\times (nx+1)} \cdot \vc{z} & = & \vc{u}\label{lgs-cyc2:eq}\;.
361\end{eqnarray}
362Die endg"ultige L"osung $\vc{m}$ ergibt sich anschlie"send zu:
363\begin{equation}
364\vc{m}=\vc{y}-\left( \frac{\vc{v}\cdot\vc{y}}{1+(\vc{v}\cdot\vc{z})}\right)\vc{z}\;.
365\end{equation}
366Die Vektoren $\vc{u}$ und $\vc{v}$ sind dabei wie folgt definiert:
367\[
368\vc{u} = \left(\begin{array}{c}
369          \gamma  \\
370           0      \\
371           \vdots \\
372           0      \\
373           \alpha\end{array}\right)\;,\hspace{1cm}
374\vc{v} = \left(\begin{array}{c}
375           1  \\
376           0      \\
377           \vdots \\
378           0      \\
379           \beta/\alpha\end{array}\right)\;.
380\]
381Der Parameter $\gamma$ kann dabei beliebig gew"ahlt werden. Er wird hier gem"a"s des
382Vorschlages von \shortciteN{press:86} zu $\gamma=-a_{0,0}$ gew"ahlt. Wird nun die
383Matrix $\mathsf{A}_{(nx+1)\times (nx+1)}$ derart modifiziert, da"s gilt:
384\[
385a'_{0,0}=a_{0,0}-\gamma\hspace{1cm}\mbox{und}\hspace{1cm}
386a'_{nx,nx}=a_{nx,nx}-\frac{\alpha\beta}{\gamma}
387\]
388und
389\[
390\mathsf{A}_{(nx+1)\times (nx+1)}=\left(\begin{array}{ccccccc}
391a'_{0,0}  & a_{0,1} &         &        &             &  \\
392a_{1,0}   & a_{1,1} & a_{1,2} &        &             &  \\
393          & \ddots  & \ddots  & \ddots &             &  \\
394          &         & \ddots  & \ddots & \ddots      &  \\
395          &         &         & \ddots & a_{nx-1,nx-1} & a_{nx-1,nx}    \\
396          &         &         &        & a_{nx,nx-1}   & a'_{nx,nx}
397\end{array}\right)\;,
398\]
399so k"onnen die Gleichungen (\ref{lgs-cyc1:eq}) und (\ref{lgs-cyc2:eq}) wie
400herk"ommliche tridiagonale Gleichungssysteme behandelt werden.
401
402\minisec{Der Thomasalgorithmus}
403
404Gel"ost werden die linearen Gleichungssysteme sowohl im nichtperiodischen als auch im periodischen
405Fall mit dem sogenannten Thomasalgorithmus.
406\begin{description}
407\item[1. Vorw"artssubstitution:]
408\begin{eqnarray}
409y_{0}   & = & c_{0}  \hspace{0.5cm}\mbox{\small nichtperiodisch}       \\[6pt]
410y_{0}   & = & c_{0}  \hspace{0.5cm}\mbox{\small periodisch}       \\[6pt]
411y_{i}   & = & c_{i} - l_{i,i-1}\cdot y_{i-1}
412                                \hspace{0.5cm}\mbox{mit}\,
413\left\{\begin{array}{lr}
414i=1,\ldots,nz & \mbox{\small nichtperiod.} \\
415i=0,\ldots,nx & \mbox{\small period.}
416\end{array}\right.
417\end{eqnarray}
418\item[2. R"uckw"artssubstitution:]
419\begin{eqnarray}
420m_{nz}   & = & \frac{y_{nz}}{d_{nz}}\hspace{0.5cm}\mbox{\small nichtperiodisch}  \\[6pt]
421m_{nx}   & = & \frac{y_{nx}}{d_{nx}}\hspace{0.5cm}\mbox{\small periodisch}  \\[6pt]
422m_{i}      & = & \frac{y_{i} - r_{i+1}\cdot m_{i+1}}{d_{i}}
423                                \hspace{0.5cm}\mbox{mit}
424\left\{\begin{array}{lr}
425i=nz-1,\ldots,0 & \mbox{\small nichtperiod.} \\
426i=nx-1,\ldots,0 & \mbox{\small period.}
427\end{array}\right.
428\end{eqnarray}
429\end{description}
430
431
432
433\section*{Filterung der Daten auf $\mathbf{2\Delta}$-Instabilit"aten}
434
435Die Verwendung des Upstream-Spline Advektionsverfahrens hat zur Folge, da"s sich
436unrealistische, numerisch bedingte $2\Delta$-Instabilit"aten ausbilden k"onnen.
437Um diese auszufiltern werden die advektiven "Anderungen der prognostischen
438Variablen im Anschlu"s an die Berechnung der Advektion in die jeweilige
439Raumrichtung einer Filterung auf diese $2\Delta$-Instabilit"aten unterzogen.
440Es wird ein sogenannter selektiver Filter von \shortciteN{mahrer:78} angewendet.
441Die ungefilterte Gr"o"se $\phi$ und die gefilterte Gr"o"se $\phi^\ast$ m"ussen
442dabei folgendes inhomogene, lineare Gleichungssystem mit symmetrischer,
443diagonaldominanter, tridiagonaler Koeffizientenmatrix erf"ullen:
444\begin{equation}\label{long-filter:eq}
445  (1-\delta)\,\phi_{i-1}^\ast + 2(1+\delta)\,\phi_i^\ast +
446  (1-\delta)\,\phi_{i+1}^\ast = \phi_{i-1} + 2\,\phi_{i} + \phi_{i+1}
447\end{equation}
448(der Laufindex $i$ m"oge f"ur die Betrachtungen von $0$ bis $n$ laufen).
449Der Wert des Gl"attungsfaktors bestimmt dabei, inwieweit von der Filterung
450auch St"orungen mit Wellenl"angen gr"o"ser als $2\Delta$ betroffen sind.
451Damit diese aber, da nicht numerisch bedingt, m"oglichst unbeeinflu"st
452bleiben sollen, sollte der Gl"attungsfaktor so klein wie m"oglich gew"ahlt werden
453($0 < \delta \le 0.1 $).
454
455Ein Problem bei der Filterung stellen die Randpunkte des Modellgebietes ($0$ und $n$)
456dar, da man f"ur sie das Gleichungssystem (\ref{long-filter:eq}) nicht formulieren
457kann. L"a"st man die Randpunkte jedoch ungefiltert, so gilt:
458\newlength{\templength}
459\setlength{\templength}{\mathindent}
460\setlength{\mathindent}{0pt}
461\begin{equation}
462  \begin{array}{rcl}
463    \multicolumn{3}{l}{\mbox{$i = 1$:}} \\
464    \hspace{0.5cm}2(1+\delta)\,\phi_1^\ast + (1-\delta)\phi_2^\ast     & = & 
465      \delta \phi_{0} + 2\,\phi_1 + \phi_2   \\[8pt]
466    \multicolumn{3}{l}{\mbox{$1\le i\le n$:}} \\[6pt]
467    \hspace{0.5cm}(1-\delta)\,\phi_{i-1}^\ast + 2(1+\delta)\,\phi_{i}^\ast +
468    (1-\delta)\,\phi_{i+1}^\ast                          & = & 
469      \phi_{i-1} + 2\,\phi_i + \phi_{i+1} \\[8pt]
470    \multicolumn{3}{l}{\mbox{$i = n-1$:}} \\   
471    \hspace{0.5cm}2(1+\delta)\,\phi_{n-1}^\ast + (1-\delta)\phi_n^\ast & = & 
472       \phi_{n-1} + 2\,\phi_{n-1} +  \delta\phi_{n}
473  \end{array}
474\end{equation}
475\setlength{\mathindent}{\templength}
476Gel"ost wird dieses Gleichungssystem mit dem bereits beschriebenen
477Thomasalgorithmus.
478
479Tests haben gezeigt, da"s die Wirkung des Filters abh"angig von der Anzahl der betrachteten
480Gitterpunkte, jedoch unabh"angig von der Formulierung f"ur die Randpunkte ist.
481
482\section*{Ablauf in PALM-1}
483
484Die Umsetzung des Upstream-Spline Advektionsverfahrens in Verbindung mit der
485Filterung auf $2\Delta$-Instabilit"aten wird durch das Flu"sdiagramm in
486Abbildung~\ref{ablauf:fig} verdeutlicht.
487\begin{figure}[p]
488\centering
489\includegraphics[width=\textwidth]{ablauf.eps}
490\caption{Ablaufplan f"ur das Upstream-Spline Advektionsverfahren}
491\label{ablauf:fig}
492\end{figure}
493
494\newpage
495\bibliography{ups}
496
497\end{document}
498
Note: See TracBrowser for help on using the repository browser.