Grid structure
* A cartesian or \(\lambda\)-\(\phi\)-\(z\) grid is used.
* The orography and other obstacles (buildings) are cut out.
* This technique can be loosely described as ”cut cell” or ”shaved cell” or ”partial step” approaches.
* The cut grid information is stored by the partial face and volume areas of each grid cell.
* This information is used for the spatial discretization.
* For the \(\lambda\)-\(\phi\)-\(z\) grid also the shape deformation is within these area information.
* To parallelize the code the grid is decomposed in blocks.
* Blocks can have different resolution.
$$\begin{matrix}
F_W & =& F_E= r_C \Delta \varphi \Delta z \\
F_S & =& r_C \Delta \lambda \cos{\varphi_S} \Delta z\\
F_N & =& r_C \Delta \lambda \cos{\varphi_N} \Delta z\\
F_B & =& r_B^2 \Delta \lambda \Delta \varphi \cos{\varphi_C} \\
F_T & =& r_T^2 \Delta \lambda \Delta \varphi \cos{\varphi_C}\\
V_C & =& r_C^2 \Delta \lambda \Delta \varphi \cos{\varphi_C} \Delta z
\end{matrix}$$
Spacial discretization
* Staggered/collocated arrangement of momentum components. (Staggered components are the main ones. Collocated components are used only as intermediate values.)
* The collocated (cell centered) momentum components are obtained by shifting the face values in both directions, \(U\,C_L=U\,F_L\) and \(U\,C_R=U\,F_R\).
* Equations are discretized in the cell centers.
* Most of the discretization is performed with the face and cell area information.
Discretization of the continuity equation:
$$V_C\dot{\rho}+F_E U_E – F_W U_W + F_T W_T – F_B W_B=0$$
and the advective part of any other equation:
$$V_C\dot{c}+F_E U_E(c/\rho)_E – F_W U_W(c/\rho)_W F_T \\[0,5cm]
W_T(c/\rho)_T – F_B W_B(c/\rho)_B=0$$
The face values \((c/\rho)_{Face}\) are determined by higher order interpolation.
Temporal discretization
The spatial discretized system is abbreviated as \(y’=F(y)\) and integrated in time by a Rosenbrock-Method, which is a higher order linear implicit integration method. The following method is proposed by [[Main_Page#Einzelnachweise|Lanser et al]] {{Anker|Lanser}}. and is of order 3.
\(\begin{matrix}
w^{n+1} &=& w^n+\frac{5}{4}k_1+\frac{3}{4}k_2 \\[0,2cm]
Sk_1 &=& \tau F(w^n) \\[0,2cm]
Sk_2 &=& \tau F(w^n+\frac{2}{3}k_1)-\frac{4}{3}k_1\\[0,2cm]
S &=&I -\gamma\,\tau\,J,\quad J=F’(w^n)
\end{matrix} \\[0,2cm] \)
with \(\gamma=\frac{1}{2}+\frac{1}{6}\sqrt{3}\).
Another representation of the method:
\(\begin{matrix}
Sk_1 &=& \tau F(w^n) \\[0,2cm]
w^n_1 &=& w^n+\frac{2}{3}k_1\\[0,2cm]
Sk_2 &=& \tau F(w^n+\frac{2}{3}k_1)-\frac{4}{3}k_1\\[0,2cm]
w^{n+1} &=& w^n_1+\frac{7}{12}k_1+\frac{3}{4}k_2
\end{matrix}\)
Error estimation by difference of two approximations of different
order.
\(w_1=w^n+k_1\, \quad w_2 = w^{n+1}=w^n+\frac{5}{4}k_1+\frac{3}{4}k_2\, \\
e=w_2-w_1=\frac{1}{4}k_1+\frac{3}{4}k_2\)
The described Rosenbrock method allows a simplified solution of the linear systems without loosing the order. If the Jacobian matrix is decomposed in two parts \(J=J_A+J_B\) the matrix \(S\) can be replaced by \(S=(I-\gamma\tau J_A)(I-\gamma\tau J_B)\). A further simplification can be reached by omitting some parts of the Jacobian or replacement of the derivatives by the same derivatives of a simplified operator: \(\tilde{F}(w^n)\). A typical example is the replacement of higher order interpolation formula by simple constant upwind value. The Jacobian: \(F’\) has the following block structure:
\(J=
\begin{pmatrix}
\frac{\partial F{_\rho}}{\partial \rho} & \frac{\partial F{_\rho}}{\partial \mathbf{V}} & 0 \\[0,2cm]
\frac{\partial F_{\mathbf{}{V}}}{\partial \rho} & \frac{\partial F_{\mathbf{V}}}{\partial \mathbf{V}} & \frac{\partial F_{\mathbf{V}}}{\partial \Theta} \\[0,2cm]
0 & \frac{\partial F_\Theta}{\partial \mathbf{V}} & \frac{\partial
F_\Theta}{\partial \Theta}
\end{pmatrix}\)
where \(0\) indicate, that this block is not included in the Jacobian. The derivative with respect to \(\rho\) is only taken for the Buoyancy term in the vertical momentum equation.
The matrix is decomposed as:
\(J=J_T+J_P=
\begin{pmatrix}
\frac{\partial F{_\rho}}{\partial \rho} & 0 & 0 \\[0,2cm]
\frac{\partial F_{\mathbf{}{V}}}{\partial \rho} & \frac{\partial F_{\mathbf{V}}}{\partial \mathbf{V}} & 0 \\[0,2cm]
0 & 0 & \frac{\partial F_\Theta}{\partial \Theta}
\end{pmatrix}
+
\begin{pmatrix}
0 & \frac{\partial F{_\rho}}{\partial \mathbf{V}} & 0 \\[0,2cm]
0 & 0 & \frac{\partial F_{\mathbf{V}}}{\partial \Theta} \\[0,2cm]
0 & \frac{\partial F_\Theta}{\partial \mathbf{V}} & 0
\end{pmatrix}\)
The first part \(J_T\) is called the transport/source part because it contains the advection, diffusion and other source terms like Coriolis, curvature, Buoyancy, latent heat, and so on. The second matrix is called the pressure part. Both systems are solved by preconditioned CG–like methods. The first linear system:
$$(I-\gamma \tau J_{AD}-\gamma \tau J_S)\Delta w =R$$
is preconditioned from the right with the matrix:
$$P_r=(I-\gamma \tau J_{AD})^{-1}$$
and from the left with the matrix
$$P_l=(I-\gamma \tau J_S)^{-1}$$
The matrix \(J_{AD}\) stands for advection/diffusion, its elements are coupled between grid cells and is the “same” for each component. The matrix \(J_S\) assembles the source terms, here the coupling is between different components in each grid cell. The new matrix
$$P_l(I-\gamma \tau J_{AD}-\gamma \tau J_S)P_r$$
can be written in the form:
$$(I-\gamma \tau P_lJ_{AD})P_r=(I+P_l((I-\gamma \tau J_{AD})+I))P_r$$
Therefore only the LU-decomposition of the matrix:
$$(I-\gamma \tau J_S)$$
has to be stored and not the matrix \(J_S\) itself. The matrix [/latex](I- \gamma \tau J_{AD})[/latex] is inverted approximately by a fixed number of Gauss-Seidel iterations.
The coefficient matrix of the second equation has the form:
$$(I-\gamma \tau J_P)=
\begin{pmatrix}
V_F & \gamma \tau Grad D_\Theta \\[0,2cm]
\gamma \tau Div D_\mathbf{V}& V_C
\end{pmatrix}$$
where \(V_F\), \(V_C\), \(D_\mathbf{V}\) and \(D_\Theta\) are diagonal matrices. Elimination of the momentum part gives a Helmholtz equation for the increment of the potential temperature. This equation is solved by a CG-method with multigrid as a preconditioner.