1 | % $Id: numerics_bc.tex 915 2012-05-30 15:11:11Z witha $ |
---|
2 | \input{header_tmp.tex} |
---|
3 | %\input{header_lectures.tex} |
---|
4 | |
---|
5 | \usepackage[utf8]{inputenc} |
---|
6 | \usepackage{ngerman} |
---|
7 | \usepackage{pgf} |
---|
8 | \usetheme{Dresden} |
---|
9 | \usepackage{subfigure} |
---|
10 | \usepackage{units} |
---|
11 | \usepackage{multimedia} |
---|
12 | \usepackage{hyperref} |
---|
13 | \newcommand{\event}[1]{\newcommand{\eventname}{#1}} |
---|
14 | \usepackage{xmpmulti} |
---|
15 | \usepackage{tikz} |
---|
16 | \usetikzlibrary{shapes,arrows,positioning} |
---|
17 | \def\Tiny{\fontsize{4pt}{4pt}\selectfont} |
---|
18 | \usepackage{amsmath} |
---|
19 | \usepackage{amssymb} |
---|
20 | \usepackage{multicol} |
---|
21 | \usepackage{pdfcomment} |
---|
22 | |
---|
23 | \institute{Institut fÌr Meteorologie und Klimatologie, Leibniz UniversitÀt Hannover} |
---|
24 | \date{last update: \today} |
---|
25 | \event{PALM Seminar} |
---|
26 | \setbeamertemplate{navigation symbols}{} |
---|
27 | |
---|
28 | \setbeamertemplate{footline} |
---|
29 | { |
---|
30 | \begin{beamercolorbox}[rightskip=-0.1cm]& |
---|
31 | {\includegraphics[height=0.65cm]{imuk_logo.pdf}\hfill \includegraphics[height=0.65cm]{luh_logo.pdf}} |
---|
32 | \end{beamercolorbox} |
---|
33 | \begin{beamercolorbox}[ht=2.5ex,dp=1.125ex, |
---|
34 | leftskip=.3cm,rightskip=0.3cm plus1fil]{title in head/foot} |
---|
35 | {\leavevmode{\usebeamerfont{author in head/foot}\insertshortauthor} \hfill \eventname \hfill \insertframenumber \; / \inserttotalframenumber} |
---|
36 | \end{beamercolorbox} |
---|
37 | \begin{beamercolorbox}[colsep=1.5pt]{lower separation line foot} |
---|
38 | \end{beamercolorbox} |
---|
39 | } |
---|
40 | %\logo{\includegraphics[width=0.3\textwidth]{luhimuk_logo.pdf}} |
---|
41 | |
---|
42 | \title[Numerics and Boundary Conditions]{Numerics and Boundary Conditions\\ |
---|
43 | (used in PALM) |
---|
44 | } |
---|
45 | \author{Siegfried Raasch} |
---|
46 | |
---|
47 | \begin{document} |
---|
48 | |
---|
49 | % Folie 1 |
---|
50 | \begin{frame} |
---|
51 | \titlepage |
---|
52 | \end{frame} |
---|
53 | |
---|
54 | % Folie 2 |
---|
55 | \begin{frame} |
---|
56 | \frametitle{Overview} |
---|
57 | \scriptsize PALM is (almost) using simple, standard and fast numerical schemes: |
---|
58 | \begin{itemize} |
---|
59 | \scriptsize |
---|
60 | \item<2-> \textbf{Spatial and temporal discretization by finite differences}\\ |
---|
61 | \item<3-> \textbf{Explicit timestep methods:}\\ |
---|
62 | - Euler\\ |
---|
63 | - Leapfrog\\ |
---|
64 | - \underline{Runge-Kutta}, second or \underline{third order} |
---|
65 | \item<4-> \textbf{Advection method}\\ |
---|
66 | - Upstream\\ |
---|
67 | - Piacsek-Williams (second order central finite differences)\\ |
---|
68 | - Upstream-spline\\ |
---|
69 | - Bott-Chlond-scheme (monotone, positiv definit, for scalars only)\\ |
---|
70 | - \underline{5th-order scheme of Wicker and Skamarock}, (as used in WRF model) |
---|
71 | \item<5-> \textbf{Poisson-equation for pressure}\\ |
---|
72 | - \underline{Direct FFT-method}\\ |
---|
73 | - Multigrid-method |
---|
74 | \item<6-> \textbf{Lagrangian particle model included} |
---|
75 | \item<7-> \textbf{Boundary conditions:}\\ |
---|
76 | - \underline{Cyclic} and non-cyclic horizontal boundary conditions\\ |
---|
77 | - Surface layer with Monin-Obukhov similarity\\ |
---|
78 | - Topography\\ |
---|
79 | - Turbulent inflow (for non-cyclic horizontal boundary conditions) |
---|
80 | |
---|
81 | \end{itemize} |
---|
82 | \end{frame} |
---|
83 | |
---|
84 | |
---|
85 | \section{Numerics} |
---|
86 | \subsection{Numerics} |
---|
87 | |
---|
88 | % Folie 3 |
---|
89 | \begin{frame} |
---|
90 | \frametitle{Numerical Grid} |
---|
91 | \footnotesize |
---|
92 | \vspace{2mm} |
---|
93 | \includegraphics[width=\textwidth]{numerics_bc_figures/numerical_grid.png} |
---|
94 | \begin{itemize} |
---|
95 | \item<1->Equations are spatially discretized on an Arakawa-C grid. |
---|
96 | \item<2->All scalar variables s (e.g. , $p^*$, $e$, $K_{\mathrm{m}}$, |
---|
97 | $K_{\mathrm{h}}$) are defined at the cell centers. |
---|
98 | \item<3->Velocity components ($u$, $v$, $w$) are shifted by half of the grid spacing. |
---|
99 | \item<4->Spacings are equidistant, stretching along $z$ is possible. |
---|
100 | \end{itemize} |
---|
101 | |
---|
102 | \tikzstyle{plain} = [rectangle, draw, fill=white!20, text width=0.33\textwidth, font=\small] |
---|
103 | \begin{tikzpicture}[remember picture, overlay] |
---|
104 | \node at (current page.north west){% |
---|
105 | \begin{tikzpicture}[overlay] |
---|
106 | \node[plain, draw,anchor=west] at (88mm,-55mm) { |
---|
107 | \noindent \scriptsize general definition (cylic):\\ |
---|
108 | $\Psi$(0:nz+1,-1:ny+1,-1:nx+1)\\ |
---|
109 | $\Psi$(:,-1,:) $=\Psi$(:,ny,:)\\ |
---|
110 | $\Psi$(:,ny+1,:) $=\Psi$(:,0,:) |
---|
111 | |
---|
112 | }; |
---|
113 | \end{tikzpicture} |
---|
114 | }; |
---|
115 | \end{tikzpicture} |
---|
116 | \end{frame} |
---|
117 | |
---|
118 | % Folie 4 |
---|
119 | \begin{frame} |
---|
120 | \frametitle{Timestep Methods (I)} |
---|
121 | \footnotesize |
---|
122 | \begin{itemize} |
---|
123 | \item<1->\textbf{Euler}\\ |
---|
124 | \vspace{1mm} |
---|
125 | $\dfrac{\partial \psi(t)}{\partial t} = F (\psi(t)) \rightarrow |
---|
126 | \dfrac{\psi(t + \Delta t) - \psi(t)}{\Delta t} |
---|
127 | \approx F (\psi(t))$ \hspace{8mm} |
---|
128 | \onslide<2-> $u\dfrac{\Delta t}{\Delta x}=C<1$\\ |
---|
129 | \begin{flushright} |
---|
130 | for stability |
---|
131 | \end{flushright} |
---|
132 | |
---|
133 | \onslide<1->$\psi (t+\Delta t) = \psi(t) |
---|
134 | + \Delta t \cdot F(\psi(t)) \hspace{28mm} \mathcal{O}(\Delta t)$\\ |
---|
135 | (used for SGS-TKE in special cases) |
---|
136 | |
---|
137 | \vspace{3mm} |
---|
138 | \item<3-> \textbf{Leapfrog}\\ |
---|
139 | $\psi(t + \Delta t) = \psi(t - \Delta t)+2 \Delta t \cdot F(\psi(t))$ \hspace{16mm} $ |
---|
140 | \mathcal{O}(\Delta t^2)$ \hspace{3mm} $C \le 0.1$\\ |
---|
141 | Time-splitting requires a weak time filter (Asselin filter) |
---|
142 | |
---|
143 | \vspace{3mm} |
---|
144 | \item<4-> \textbf{Runge-Kutta, third-order}\\ |
---|
145 | $k_1=F(\psi(t))$\\ |
---|
146 | \vspace{1mm} |
---|
147 | $k_2=F \left( \psi(t) + \frac{1}{3} \Delta t \cdot k_1 \right)$\\ |
---|
148 | \vspace{1mm} |
---|
149 | $k_3=F \left( \psi(t) - \frac{3}{16} \Delta t \cdot k_1 |
---|
150 | + \frac{15}{16} \Delta t \cdot k_2 \right)$\\ |
---|
151 | \vspace{1mm} |
---|
152 | $\psi(t + \Delta t) = \psi(t) + \frac{1}{30}\Delta t (5 k_1 + 9 k_2 + 16 k_3)$ |
---|
153 | \hspace{12mm} $\mathcal{O}(\Delta t^2)$ \hspace{3mm} $C \le 0.9$\\ |
---|
154 | \end{itemize} |
---|
155 | \end{frame} |
---|
156 | |
---|
157 | % Folie 5 |
---|
158 | \begin{frame} |
---|
159 | \frametitle{Timestep Methods (II)} |
---|
160 | \footnotesize |
---|
161 | \onslide<1->In the PALM code, the different timestep schemes are treated by one |
---|
162 | equation using switches: |
---|
163 | $\psi (t + \Delta t ) = (1 - c_1) \cdot \psi (t - \Delta t ) + c_1 \cdot \psi (t) |
---|
164 | + \Delta t \cdot \left[ c_2 \cdot F (\psi (t) ) + c_3 \cdot F (\psi (t - |
---|
165 | \Delta t ) ) \right]$ |
---|
166 | \vspace{1mm} |
---|
167 | |
---|
168 | \onslide<2-> |
---|
169 | \begin{centering} |
---|
170 | \begin{table} |
---|
171 | \begin{tabular}{cccc} |
---|
172 | \bf{Scheme} & \bf{c$_1$} & \bf{c$_2$} & \bf{c$_3$}\\ |
---|
173 | Euler & 1 & 1 & 0\\ |
---|
174 | Leapfrog & 0 & 2 & 0\\ |
---|
175 | RK (1st step) & 1 & 1/3 & 0\\ |
---|
176 | RK (2nd step) & 1 & 15/16 & -25/48\\ |
---|
177 | RK (3rd step) & 1 & 8/15 & 1/15\\ |
---|
178 | \end{tabular} |
---|
179 | \end{table} |
---|
180 | \end{centering} |
---|
181 | |
---|
182 | \onslide<3-> |
---|
183 | \begin{align*} |
---|
184 | \psi (t - \Delta t) &= \psi (t) \hspace{15mm} \textbf{after each RK substep}\\ |
---|
185 | \psi (t) &= \psi (t + \Delta t) |
---|
186 | \end{align*} |
---|
187 | \end{frame} |
---|
188 | |
---|
189 | % Folie 6 |
---|
190 | \begin{frame} |
---|
191 | \frametitle{Advection Methods (I)} |
---|
192 | \small |
---|
193 | \begin{itemize} |
---|
194 | \item<1-> Piacsek Williams C3 (1970, J. Comput. Phy., 6, 392). |
---|
195 | \item<2-> Scheme of 2nd order accuracy. |
---|
196 | \item<3-> Conserves integrals of linear and quadratic quantities. |
---|
197 | \item<4-> Requires incompressibility $\rightarrow$ flux form of advection term. |
---|
198 | \onslide<4-> \includegraphics[width=0.8\textwidth]{numerics_bc_figures/advection_methods.png} |
---|
199 | \end{itemize} |
---|
200 | $$\left.\frac{\partial (u \psi)}{\partial x}\right\vert_i = \frac{1}{2 \Delta x} |
---|
201 | \left( u_{i+\frac{1}{2}} \psi_{i+1} - u_{i-\frac{1}{2}} \psi_{i-1} \right)$$ |
---|
202 | \begin{itemize} |
---|
203 | \item<5-> In case of momentum advection (e.g. $\psi=u$), $u_{i-1}$ and |
---|
204 | $u_{i+1}$ have to be obtained by linear interpolation. |
---|
205 | \item<5-> May cause $2 \Delta x$ wiggles in case of sharp gradients. |
---|
206 | \end{itemize} |
---|
207 | \end{frame} |
---|
208 | |
---|
209 | % Folie 7 |
---|
210 | \begin{frame} |
---|
211 | \frametitle{Advection Methods (II)} |
---|
212 | \begin{itemize} |
---|
213 | \item<1-> \small Upstream-spline\\ \scriptsize |
---|
214 | \onslide<1-> - Requires Euler timestep\\ |
---|
215 | \onslide<2-> - Nonlocal scheme: produces heavy load on communication network\\ |
---|
216 | \onslide<3-> - Numerically unstable under stable stratification |
---|
217 | \item<4-> \small Bott-Chlond\\ \scriptsize |
---|
218 | \onslide<4-> - Chlond (1994)\\ |
---|
219 | \onslide<5-> - Monotone, positive definit. Can only be used for scalars\\ |
---|
220 | \onslide<6-> - Conserves sharp gradients\\ |
---|
221 | \onslide<7-> - Numerically expensive\\ |
---|
222 | \onslide<8-> - Not optimized for use on cache-based machines. |
---|
223 | \item<9-> \small Default: Wicker and Skamarock scheme (5th order)\\ \scriptsize |
---|
224 | \onslide<9-> - Much better accuracy than Piacsek Williams\\ |
---|
225 | \onslide<10-> - Much simpler algorithm than Bott-Chlond\\ |
---|
226 | \onslide<11-> - Requires additional ghost layers\\ |
---|
227 | \onslide<12-> - Adds additional numerical dissipation |
---|
228 | |
---|
229 | \end{itemize} |
---|
230 | \end{frame} |
---|
231 | |
---|
232 | % Folie 8 |
---|
233 | \begin{frame} |
---|
234 | \frametitle{Advection Methods â Wicker/Skamarock (I)} |
---|
235 | \footnotesize |
---|
236 | \begin{itemize} |
---|
237 | \item Wicker and Skamarock (2002, Mon. Wea. Rev. 130, 2088 â 2097). |
---|
238 | \item Based on flux form of advection term |
---|
239 | \item Difference of fluxes at the edge of the grid cell is used to |
---|
240 | discretise advection term |
---|
241 | \end{itemize} |
---|
242 | |
---|
243 | \begin{columns}[T] |
---|
244 | \begin{column}{0.55\textwidth} |
---|
245 | \hspace{8mm}\includegraphics[width=0.8\textwidth]{numerics_bc_figures/numerical_grid_small.png} |
---|
246 | \end{column} |
---|
247 | \begin{column}{0.45\textwidth} |
---|
248 | $\frac{ \partial \psi}{\partial t} = - \nabla (u_i \psi) \approx |
---|
249 | - \frac{F_{i+\frac{1}{2}} - F_{i-\frac{1}{2}}}{\Delta x}$ |
---|
250 | \end{column} |
---|
251 | \end{columns} |
---|
252 | |
---|
253 | \tikzstyle{plain} = [rectangle, text width=0.1\textwidth, font=\small] |
---|
254 | \begin{tikzpicture}[remember picture, overlay] |
---|
255 | \node at (current page.north west){% |
---|
256 | |
---|
257 | \begin{tikzpicture}[overlay] |
---|
258 | \node[plain, anchor=west] at (2mm,-68mm) { |
---|
259 | \tikz |
---|
260 | { |
---|
261 | \draw[blue, -latex', line width=5pt] (1,0) -- (2,0); |
---|
262 | } |
---|
263 | |
---|
264 | $F_{i-\frac{1}{2}}$ |
---|
265 | }; |
---|
266 | \end{tikzpicture} |
---|
267 | |
---|
268 | \begin{tikzpicture}[overlay] |
---|
269 | \node[plain, anchor=west] at (62mm,-68mm) { |
---|
270 | \tikz |
---|
271 | { |
---|
272 | \draw[blue, -latex', line width=5pt] (,0) -- (2,0); |
---|
273 | } |
---|
274 | $F_{i+\frac{1}{2}}$ |
---|
275 | }; |
---|
276 | \end{tikzpicture} |
---|
277 | |
---|
278 | }; |
---|
279 | \end{tikzpicture} |
---|
280 | \end{frame} |
---|
281 | |
---|
282 | % Folie 9 |
---|
283 | \begin{frame} |
---|
284 | \frametitle{Advection Methods â Wicker/Skamarock (II)} |
---|
285 | |
---|
286 | |
---|
287 | |
---|
288 | \textbf{Finite difference approxiamtion of 6$^{\text{th}}$ order} |
---|
289 | \begin{tikzpicture}[scale=2] |
---|
290 | \tikzstyle{ann} = [draw=none,fill=none,right] |
---|
291 | \matrix[nodes={draw, thick, fill=blue!20}, row sep=0.3cm,column sep=0.5cm]{ |
---|
292 | \node[rectangle, rounded corners]{ |
---|
293 | $F^{\text{6th}}_{i-\frac{1}{2}} = \frac{1}{60} u_{i-\frac{1}{2}} \left( 37 |
---|
294 | (\Psi_i + \Psi_{i-1}) - 8 (\Psi_{i+1} + \Psi_{i-2}) + (\Psi_{i+2} |
---|
295 | + \Psi_{i-3}) \right)$ |
---|
296 | };\\ |
---|
297 | }; |
---|
298 | \end{tikzpicture} |
---|
299 | |
---|
300 | \vspace{5mm} |
---|
301 | |
---|
302 | \textbf{Artificially added numerical dissipation term} |
---|
303 | \begin{tikzpicture}[scale=2] |
---|
304 | \tikzstyle{ann} = [draw=none,fill=none,right] |
---|
305 | \matrix[nodes={draw, thick, fill=blue!40}, row sep=0.3cm,column sep=0.5cm]{ |
---|
306 | \node[rectangle, rounded corners]{ |
---|
307 | $\frac{-1}{60} \left| u_{i-\frac{1}{2}} \right| \left( 10 (\Psi_i - |
---|
308 | \Psi_{i-1}) - 5 (\Psi_{i+1} - \Psi_{i-2}) + (\Psi_{i+2} - \Psi_{i-3}) \right)$ |
---|
309 | };\\ |
---|
310 | }; |
---|
311 | \end{tikzpicture} |
---|
312 | \end{frame} |
---|
313 | |
---|
314 | % Folie 10 |
---|
315 | \begin{frame} |
---|
316 | \frametitle{Advection Methods â Wicker/Skamarock (III)} |
---|
317 | |
---|
318 | \begin{tikzpicture}[scale=2] |
---|
319 | \tikzstyle{ann} = [draw=none,fill=none,right] |
---|
320 | \matrix[nodes={draw, thick, fill=blue!20}, row sep=0.3cm,column sep=0.5cm]{ |
---|
321 | \node[rectangle, rounded corners]{ |
---|
322 | $F^{\text{6th}}_{i-\frac{1}{2}} = \frac{1}{60} u_{i-\frac{1}{2}} |
---|
323 | \left( 37 (\Psi_i + \Psi_{i-1}) - 8 (\Psi_{i+1} + \Psi_{i-2}) + (\Psi_{i+2} + |
---|
324 | \Psi_{i-3}) \right)$ |
---|
325 | };\\ |
---|
326 | }; |
---|
327 | \end{tikzpicture} |
---|
328 | |
---|
329 | \begin{columns}[T] |
---|
330 | \begin{column}{0.7\textwidth} |
---|
331 | \includegraphics[width=1\textwidth]{numerics_bc_figures/numerical_oscillations.png} |
---|
332 | \end{column} |
---|
333 | \begin{column}{0.3\textwidth} |
---|
334 | Centered Finite Differences produces numerical oscillations (''wiggles'') |
---|
335 | near sharp gradients. |
---|
336 | \end{column} |
---|
337 | \end{columns} |
---|
338 | \end{frame} |
---|
339 | |
---|
340 | % Folie 11 |
---|
341 | \begin{frame} |
---|
342 | \frametitle{Advection Methods â Wicker/Skamarock (IV)} |
---|
343 | \footnotesize |
---|
344 | |
---|
345 | \begin{columns}[T] |
---|
346 | \begin{column}{0.1\textwidth} |
---|
347 | \begin{tikzpicture}[scale=2] |
---|
348 | \tikzstyle{ann} = [draw=none,fill=none,right] |
---|
349 | \matrix[nodes={draw, thick, fill=blue!40}, row sep=0.3cm,column sep=0.5cm]{ |
---|
350 | \node[draw=none,fill=none]{ |
---|
351 | $F^{\text{5th}}_{i-\frac{1}{2}} = F^{\text{6th}}_{i-\frac{1}{2}}$ |
---|
352 | };\\ |
---|
353 | }; |
---|
354 | \end{tikzpicture} |
---|
355 | \end{column} |
---|
356 | \begin{column}{0.8\textwidth} |
---|
357 | \begin{tikzpicture}[scale=2] |
---|
358 | \tikzstyle{ann} = [draw=none,fill=none,right] |
---|
359 | \matrix[nodes={draw, thick, fill=blue!40}, row sep=0.3cm,column sep=0.5cm]{ |
---|
360 | \node[rectangle, rounded corners]{ |
---|
361 | $- \frac{1}{60} \left| u_{i-\frac{1}{2}} \right| \left( 10 (\Psi_i - \Psi_{i-1}) - |
---|
362 | 5 (\Psi_{i+1} - \Psi_{i-2}) + (\Psi_{i+2} - \Psi_{i-3}) \right)$ |
---|
363 | };\\ |
---|
364 | }; |
---|
365 | \end{tikzpicture} |
---|
366 | \end{column} |
---|
367 | \end{columns} |
---|
368 | |
---|
369 | \begin{columns}[T] |
---|
370 | \begin{column}{0.7\textwidth} |
---|
371 | \includegraphics[width=1\textwidth]{numerics_bc_figures/numerical_oscillations_2.png} |
---|
372 | \end{column} |
---|
373 | \begin{column}{0.3\textwidth} |
---|
374 | \vspace{3mm} |
---|
375 | \underline{Advantage}\\ |
---|
376 | Numerical Dissipation damps small scale oscillations.\\ |
---|
377 | \vspace{3mm} |
---|
378 | \underline{Disadvantage}\\ |
---|
379 | In a turbulent flow numerical dissipation removes energy from small scales. |
---|
380 | |
---|
381 | \end{column} |
---|
382 | \end{columns} |
---|
383 | \end{frame} |
---|
384 | |
---|
385 | % Folie 12 |
---|
386 | \begin{frame} |
---|
387 | \frametitle{Advection Methods â Wicker/Skamarock (V)} |
---|
388 | |
---|
389 | \begin{columns}[T] |
---|
390 | \begin{column}{0.6\textwidth} |
---|
391 | \includegraphics[width=1\textwidth]{numerics_bc_figures/numerical_properties.png} |
---|
392 | \end{column} |
---|
393 | \begin{column}{0.4\textwidth} |
---|
394 | \includegraphics[width=1\textwidth]{numerics_bc_figures/pw_ws.png} |
---|
395 | \end{column} |
---|
396 | \end{columns} |
---|
397 | |
---|
398 | \begin{itemize} |
---|
399 | \item Better resolution of larger scales $(> 8\,\Delta x)$ and hence less |
---|
400 | numerical energy transfer from larger to smaller scales compared to lower |
---|
401 | order schemes. |
---|
402 | \item Less spectral energy at smaller scales. |
---|
403 | \end{itemize} |
---|
404 | |
---|
405 | \end{frame} |
---|
406 | |
---|
407 | % Folie 13 |
---|
408 | \begin{frame} |
---|
409 | \frametitle{Pressure Solver (I)} |
---|
410 | \footnotesize |
---|
411 | \begin{itemize} |
---|
412 | \item<1-> Governing equations of PALM require incompressibility |
---|
413 | \item<2-> Incompressibility is reached by a predictor-corrector method\\ |
---|
414 | \scriptsize |
---|
415 | 1. Momentum equations are solved without the pressure term giving a |
---|
416 | provisional velocity field which is not free of divergence.\\ |
---|
417 | $\overline{u}^{t+\Delta t}_{i_{\mathrm{prov}}} = \overline{u}^t_i + |
---|
418 | \Delta t \left( - \frac{\partial}{\partial x_k} \overline{u}^t_k |
---|
419 | \overline{u}^t_i - (\varepsilon_{ijk} f_j \overline{u}^t_k |
---|
420 | - \varepsilon_{i3k} f_3 u_{\mathrm{g}_k}) |
---|
421 | + g \frac{\overline{\theta^*}^t}{\theta_0} \delta_{i3} |
---|
422 | - \frac{\partial}{\partial x_k} \overline{u'_k u'_i}^t \right)$\\ |
---|
423 | \vspace{3mm} |
---|
424 | \onslide<3-> 2. Assign all remaining divergences to the (perturbation) |
---|
425 | pressure $p^*$ so that the new corrected velocity field is the sum of the |
---|
426 | provisional, divergent field and the perturbation pressure term.\\ |
---|
427 | $\overline{u}^{t+\Delta t}_{i} = |
---|
428 | \overline{u}^{t+\Delta t}_{i_{\mathrm{prov}}} - |
---|
429 | \frac{\Delta t}{\rho_0} \frac{\partial \overline{p^*}^t}{\partial x_i}$\\ |
---|
430 | \vspace{3mm} |
---|
431 | \onslide<4-> 3. The divergence operator is applied to this equation. |
---|
432 | Demanding a corrected velocity field free of divergence, this leads to a |
---|
433 | Poisson equation for the perturbation pressure.\\ |
---|
434 | $\frac{\partial^2 \overline{p^*}^t}{\partial x_i^2} = \frac{\rho_0} |
---|
435 | {\Delta t} \frac{\partial \overline{u}_{i_{\mathrm{prov}}}^{t + \Delta t}} |
---|
436 | {\partial x_i}$\\ |
---|
437 | \vspace{3mm} |
---|
438 | \onslide<5-> 4. After solving the Poisson equation, the final velocity field |
---|
439 | is \\ |
---|
440 | calculated as given in step 2.\\ |
---|
441 | |
---|
442 | \end{itemize} |
---|
443 | \end{frame} |
---|
444 | |
---|
445 | % Folie 14 |
---|
446 | \begin{frame} |
---|
447 | \frametitle{Pressure Solver (II)} |
---|
448 | \small |
---|
449 | |
---|
450 | \begin{itemize} |
---|
451 | \item FFT-solver\\ |
---|
452 | \onslide<1-> 1. Discretization of the Poisson-equation by central differences\\ |
---|
453 | \onslide<2-> 2. 2D discrete FFT in both horizontal directions\\ |
---|
454 | \onslide<3-> 3. Solving the resulting tridiagonal set of linear equations\\ |
---|
455 | \onslide<4-> 4. Inverse 2D discrete FFT in both horizontal directions leading |
---|
456 | to the perturbation pressure |
---|
457 | |
---|
458 | \begin{itemize} |
---|
459 | \item<5-> Very fast and accurate, $\mathcal{O}(n \log n)$, $n$: |
---|
460 | number of gridpoints |
---|
461 | \item<6-> CPU requirement $<$ 50\% of total CPU time |
---|
462 | \item<7-> Due to non-locality of the FFT, transpositions are required |
---|
463 | on parallel computers |
---|
464 | \item<8-> Usage requires periodic boundary conditions and uniform grids |
---|
465 | along $x$ and $y$ |
---|
466 | \end{itemize} |
---|
467 | \end{itemize} |
---|
468 | \end{frame} |
---|
469 | |
---|
470 | % Folie 15 |
---|
471 | \begin{frame} |
---|
472 | \frametitle{Pressure Solver (III)} |
---|
473 | \scriptsize |
---|
474 | \begin{columns} |
---|
475 | \begin{column}{1.03\textwidth} |
---|
476 | \begin{itemize} |
---|
477 | \item<1-> Multigrid-method\\ |
---|
478 | \begin{itemize}\scriptsize |
---|
479 | \item Iterative solver\\ |
---|
480 | \scriptsize |
---|
481 | basic idea: Poisson equation is transformed to a fixed point problem:\\ |
---|
482 | $\vec{p}^{k+1} = T \cdot \vec{p}^k + \vec{c}^k$\\ |
---|
483 | \vspace{1mm} |
---|
484 | \onslide<2-> starting from a first guess, the solution will be |
---|
485 | improved by repeated execution of the fixed point problem:\\ |
---|
486 | $\begin{array}{rcl} |
---|
487 | \vec{p}^{1} &=&T \cdot \vec{p}^0 + \vec{c}^0\\ |
---|
488 | \vec{p}^{2} &=&T \cdot \vec{p}^1 + \vec{c}^1 |
---|
489 | \vspace{-2mm}\\ |
---|
490 | |
---|
491 | &\vdots& |
---|
492 | \vspace{-1.5mm}\\ |
---|
493 | |
---|
494 | \vec{p}^{k} &=&T \cdot \vec{p}^{k-1} + \vec{c}^{k-1}\\ |
---|
495 | \vec{p}^{k+1} &=&T \cdot \vec{p}^k + \vec{c}^k\\ |
---|
496 | \end{array}$\\ |
---|
497 | \vspace{1mm} |
---|
498 | \onslide<3-> Depending on the structure of the matrix $T$ and vector |
---|
499 | $c$ different iterative solvers can be defined, e.g.: Jacobi-scheme |
---|
500 | (here on 2D-uniform grid):\\ |
---|
501 | $p^{k+1}_{i,j} = \frac{1}{4} \cdot \left( p^k_{i-1,j} + p^k_{i+1,j} |
---|
502 | + p^k_{i,j-1} + p^k_{i,j+1} - \Delta x^2 f(i,j,k) \right)$\\ |
---|
503 | \scriptsize \vspace{2mm} |
---|
504 | \item<4-> With each iteration step $k$ the improved solution converges towards |
---|
505 | the exact solution. |
---|
506 | \item<5-> Iterative schemes are 'local schemes' $\rightarrow$ information |
---|
507 | is needed \\ only from neighboring grid-points. |
---|
508 | \item<6-> Very low convergence: $\mathcal{O}(n^2)$. |
---|
509 | \end{itemize} |
---|
510 | \end{itemize} |
---|
511 | \end{column} |
---|
512 | \end{columns} |
---|
513 | \end{frame} |
---|
514 | |
---|
515 | % Folie 16 |
---|
516 | \begin{frame} |
---|
517 | \frametitle{Pressure Solver (IV)} |
---|
518 | \begin{itemize} |
---|
519 | \item<1-> Multigrid-method\\ |
---|
520 | \begin{itemize} |
---|
521 | \item Due to their locality, iterative solvers show a frequency-dependent |
---|
522 | reduction of the residual: low frequencies are reduced slower than high |
---|
523 | frequencies. |
---|
524 | \item<2-> The main idea of the multigrid method is to reduce errors of different |
---|
525 | frequencies on grids with different grid spacing: |
---|
526 | \begin{itemize} |
---|
527 | \item errors of high frequency are reduced on fine grids |
---|
528 | \item errors of low frequency are reduced on coarse grids. |
---|
529 | \end{itemize} |
---|
530 | \end{itemize} |
---|
531 | \end{itemize} |
---|
532 | \onslide<2-> |
---|
533 | \begin{figure}[htp] |
---|
534 | \centering |
---|
535 | \includegraphics[scale=0.35]{numerics_bc_figures/errors.png} |
---|
536 | \end{figure} |
---|
537 | \end{frame} |
---|
538 | |
---|
539 | % Folie 17 |
---|
540 | \begin{frame} |
---|
541 | \frametitle{Pressure Solver (V)} |
---|
542 | |
---|
543 | \begin{columns}[T] |
---|
544 | \begin{column}{0.65\textwidth} |
---|
545 | \begin{itemize} |
---|
546 | \item<1-> Multigrid-method\\ |
---|
547 | \begin{itemize} |
---|
548 | \footnotesize |
---|
549 | \item On each grid-level an approximate solution of the fixed point |
---|
550 | equation is obtained performing a few iterations. |
---|
551 | \item<2-> The solution is transmitted to the next coarser grid-level |
---|
552 | where it is used as the first guess to solve the fixed point problem. |
---|
553 | \item<3-> This procedure is performed up to the coarsest grid-level |
---|
554 | containing two grid-points in each direction. |
---|
555 | \item<4-> From the coarsest grid-level the procedure is passed in |
---|
556 | backward order to get the final solution. |
---|
557 | \item<5-> For large grids faster than FFT method. |
---|
558 | \item<6-> V- and W-cycles are implemented. |
---|
559 | \end{itemize} |
---|
560 | \end{itemize} |
---|
561 | \end{column} |
---|
562 | \begin{column}{0.5\textwidth} |
---|
563 | \onslide<2-> |
---|
564 | \includegraphics[width=1\textwidth]{numerics_bc_figures/multigrid.png} |
---|
565 | \end{column} |
---|
566 | \end{columns} |
---|
567 | \end{frame} |
---|
568 | |
---|
569 | |
---|
570 | |
---|
571 | \section{Boundary Conditions} |
---|
572 | \subsection{Boundary Conditions} |
---|
573 | |
---|
574 | % Folie 18 |
---|
575 | \begin{frame} |
---|
576 | \frametitle{Boundary Conditions (I)} |
---|
577 | |
---|
578 | \begin{itemize} |
---|
579 | \item<1-> Lateral $(xy)$ boundary conditions:\\ |
---|
580 | \begin{itemize} |
---|
581 | \item Cyclic by default, allowing undisturbed evolution / advection of turbulence. |
---|
582 | \begin{columns}[T] |
---|
583 | \begin{column}{0.2\textwidth} |
---|
584 | %leer |
---|
585 | \end{column} |
---|
586 | \begin{column}{0.5\textwidth} |
---|
587 | \includegraphics[width=1\textwidth]{numerics_bc_figures/lateral_bc.png}\\ |
---|
588 | \vspace{2mm} |
---|
589 | \end{column} |
---|
590 | \begin{column}{0.4\textwidth} |
---|
591 | $\begin{array}{rcl} |
---|
592 | \Psi(-1) &=& \Psi(n)\\ |
---|
593 | \Psi(n+1) &=& \Psi(0)\\ |
---|
594 | \end{array}$ |
---|
595 | \end{column} |
---|
596 | \end{columns} |
---|
597 | |
---|
598 | \item<2-> Dirichlet (inflow) and radiation (outflow) conditions are allowed along |
---|
599 | either $x$- or $y$-direction. |
---|
600 | |
---|
601 | \item<3-> In case of a Dirichlet condition, the inflow is laminar (by default) and |
---|
602 | the domain has to be extended to allow for the development of a turbulent |
---|
603 | state, if neccessary. |
---|
604 | \item<4-> Non-cyclic lateral conditions require the use of the multigrid-method |
---|
605 | for solving the Poisson-equation. |
---|
606 | |
---|
607 | \end{itemize} |
---|
608 | \end{itemize} |
---|
609 | \end{frame} |
---|
610 | |
---|
611 | % Folie 19 |
---|
612 | \begin{frame} |
---|
613 | \frametitle{Boundary Conditions (II)} |
---|
614 | |
---|
615 | \begin{columns}[T] |
---|
616 | \begin{column}{0.7\textwidth} |
---|
617 | \scriptsize |
---|
618 | \begin{itemize} |
---|
619 | \item<1-> Surface boundary condition: |
---|
620 | \begin{itemize} |
---|
621 | \scriptsize |
---|
622 | \item<1-> Monin-Obukhov-similarity is used by default, i.e. a Prandtl-layer is |
---|
623 | assumed between the surface and the first grid layer.\\ |
---|
624 | $\frac{\partial \overline{u}}{\partial z} = \frac{u_{*}}{\kappa z} |
---|
625 | \Phi_{\mathrm{m}}$; \hspace{3mm} $u_{*} = \sqrt{- \overline{w' u'_0}} |
---|
626 | = \sqrt{\frac{\tau_0}{\overline{\rho}}}$\\ |
---|
627 | $\frac{\partial \overline{\theta}}{\partial z} = \frac{\vartheta_{*}} |
---|
628 | {\kappa z} \Phi_{\mathrm{h}}$; \hspace{3mm} $\vartheta_{*} = |
---|
629 | \frac{\overline{w' \theta'_0}}{u_{*}}$\\ |
---|
630 | \vspace{2mm} |
---|
631 | \item<2-> Integration between $z=z_0$ (roughness height) and $z=z_{\mathrm{p}}$ |
---|
632 | (top of Prandtl-layer, $k=1$) gives the only unknowns $u_{*}$ and |
---|
633 | $\theta_{*}$ which then define the surface momentum and heat flux, |
---|
634 | used as the real boundary conditions.\\ |
---|
635 | \vspace{2mm} |
---|
636 | \item<3-> $\Phi_{\mathrm{m}}$, $\Phi_{\mathrm{h}}$: Dyer-Businger functions\\ |
---|
637 | \onslide<4->$\Phi_{\mathrm{m}} = \left\{ \begin{array}{cc} |
---|
638 | 1+5\,\mathrm{Rif} & \text{stable}\\ |
---|
639 | 1 & \text{neutral}\\ |
---|
640 | (1-16\,\mathrm{Rif})^{-1/4} & \text{unstable}\\ |
---|
641 | \end{array} \right. $ |
---|
642 | |
---|
643 | \end{itemize} |
---|
644 | \end{itemize} |
---|
645 | \end{column} |
---|
646 | \begin{column}{0.3\textwidth} |
---|
647 | \onslide<1-> |
---|
648 | \includegraphics[width=1\textwidth]{numerics_bc_figures/surface_bc.png}\\ |
---|
649 | \vspace{-2mm} |
---|
650 | Prandtl-layer\\ |
---|
651 | \vspace{8mm} |
---|
652 | \onslide<5-> |
---|
653 | $\mathrm{Rif} = \frac{\frac{g}{\tilde{\theta}} |
---|
654 | \overline{w' \theta'_0}}{\overline{w' u'} |
---|
655 | \frac{\partial \overline{u}}{\partial z}}$\\ |
---|
656 | \scriptsize Richardson flux number |
---|
657 | \end{column} |
---|
658 | \end{columns} |
---|
659 | \end{frame} |
---|
660 | |
---|
661 | % Folie 20 |
---|
662 | \begin{frame} |
---|
663 | \frametitle{Boundary Conditions (III)} |
---|
664 | |
---|
665 | \begin{itemize} |
---|
666 | \item<1-> Surface boundary condition: |
---|
667 | \begin{itemize} |
---|
668 | \footnotesize |
---|
669 | \item<1-> Monin-Obukhov-similarity is only valid for a |
---|
670 | horizontal surface with homogeneous conditions. |
---|
671 | \item<2-> The surface temperature has to be prescribed. |
---|
672 | Alternatively, the surface heat flux can be prescribed. |
---|
673 | \item<3-> Instead of MO-similarity, no-slip conditions or |
---|
674 | free-slip conditions can be used |
---|
675 | \begin{equation*} |
---|
676 | u(z=0) = 0, \quad v(z=0)=0 \hspace{15mm} |
---|
677 | \frac{\partial u}{\partial z} = 0, \quad |
---|
678 | \frac{\partial v}{\partial z} = 0 \qquad |
---|
679 | \end{equation*} |
---|
680 | realized by |
---|
681 | \begin{flalign*} |
---|
682 | u(k=0)=-u(k=1) \hspace{18mm} u(k=0)=u(k=1)\\ |
---|
683 | v(k=0)=-v(k=1) \hspace{18mm} v(k=0)=v(k=1) |
---|
684 | \end{flalign*} |
---|
685 | \item<4-> Pressure boundary condition: |
---|
686 | $\dfrac{\partial p}{\partial z} = 0$ |
---|
687 | in order to guarantee $w(z=0)=0$ |
---|
688 | \item<5-> SGS-TKE condition $\dfrac{\partial e}{\partial z}=0$ |
---|
689 | \end{itemize} |
---|
690 | \end{itemize} |
---|
691 | \end{frame} |
---|
692 | |
---|
693 | % Folie 21 |
---|
694 | \begin{frame} |
---|
695 | \frametitle{Boundary Conditions (IV)} |
---|
696 | |
---|
697 | \begin{itemize} |
---|
698 | \item<1-> Boundary conditions at the top (default) |
---|
699 | \begin{itemize} |
---|
700 | \item<1-> Dirichlet conditions for velocities: $u=u_{\mathrm{g}}, |
---|
701 | \quad v=v_{\mathrm{g}}, \quad w=0$ |
---|
702 | \item<2-> Neumann conditions (temporal constant gradients) for scalars: |
---|
703 | $$\frac{\partial \theta}{\partial z} = |
---|
704 | \left. \frac{\partial \theta}{\partial z} \right\vert_{t=0}$$ |
---|
705 | \item<3-> Pressure: Dirichlet $p=0$ |
---|
706 | or Neumann $\dfrac{\partial p}{\partial z} = 0$ |
---|
707 | \item<4-> SGS-TKE: Neumann $\dfrac{\partial e}{\partial z} = 0$ |
---|
708 | \item<5-> A damping layer can be switched on in order to absorb |
---|
709 | gravity waves. |
---|
710 | \end{itemize} |
---|
711 | \end{itemize} |
---|
712 | \end{frame} |
---|
713 | |
---|
714 | % Folie 22 |
---|
715 | \begin{frame} |
---|
716 | \frametitle{Initial Conditions} |
---|
717 | |
---|
718 | All 3D-arrays are initialized with vertical profiles (horizontally homogeneous).\\ |
---|
719 | \quad \\ |
---|
720 | Two different profiles can be chosen: |
---|
721 | \begin{itemize} |
---|
722 | \item<2-> \textbf{constant (piecewise linear) profiles} |
---|
723 | \begin{itemize} |
---|
724 | \footnotesize |
---|
725 | \item \textbf{e.g.} $u=0, v=0, \dfrac{\partial \theta}{\partial z}=0$ |
---|
726 | \textbf{up to} $z=\unit[1000]{m}$, |
---|
727 | $\dfrac{\partial \theta}{\partial z}=+1.0$ \textbf{up to top} |
---|
728 | \end{itemize} |
---|
729 | \item<3-> \textbf{velocity profiles calculated by a 1D-model (which is a part of PALM)} |
---|
730 | \begin{itemize} |
---|
731 | \footnotesize |
---|
732 | \item \textbf{constant (piecewise linear) temperature profile is used |
---|
733 | for the 1D-model} |
---|
734 | \end{itemize} |
---|
735 | \end{itemize} |
---|
736 | \onslide<4-> |
---|
737 | \underline{Under horizontally homogeneous initial conditions, random}\\ |
---|
738 | \underline{fluctuations have to be added in order to generate turbulence!} |
---|
739 | \end{frame} |
---|
740 | |
---|
741 | \end{document} |
---|