# Annotation of /branches/doubleplusgood/doc/user/lumping.tex

Revision 4345 - (hide annotations)
Fri Mar 29 07:09:41 2013 UTC (6 years ago) by jfenwick
File MIME type: application/x-tex
File size: 11840 byte(s)
Spelling fixes

 1 gross 3379 2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 jfenwick 4154 % Copyright (c) 2003-2013 by University of Queensland 4 jfenwick 3989 5 gross 3379 % 6 % Primary Business: Queensland, Australia 7 % Licensed under the Open Software License version 3.0 8 9 % 10 jfenwick 3989 % Development until 2012 by Earth Systems Science Computational Center (ESSCC) 11 % Development since 2012 by School of Earth Sciences 12 % 13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 14 gross 3379 15 16 \section{Some Remarks on Lumping} 17 gross 4052 \label{LUMPING} 18 gross 3379 19 ahallam 3381 Explicit time integration schemes (two examples are discussed later in this 20 section), require very small time steps in order to maintain numerical stability. 21 Unfortunately, these small time increments often result in a prohibitive 22 computational cost. 23 In order to minimise these costs, a technique termed lumping can be utilised. 24 Lumping is applied to the coefficient matrix, reducing it to a simple diagonal 25 matrix. This can significantly improve the computational speed, because the 26 solution updates are simplified to a simple component-by-component 27 vector-vector product. However, some care is required when making radical 28 approximations such as these. In this section, two commonly applied lumping 29 techniques are discussed, namely row sum lumping 30 gross 3379 \index{linear solver!row sum lumping}\index{row sum lumping} and HRZ 31 lumping \index{linear solver!HRZ lumping}\index{HRZ lumping}. 32 33 ahallam 3381 \subsection{Scalar wave equation} 34 One example where lumping can be applied to a hyperpolic problem, is 35 the scalar wave equation 36 gross 3379 \begin{eqnarray} \label{LUMPING WAVE} 37 u_{,tt}=c^2 u_{,ii} \; . 38 \end{eqnarray} 39 ahallam 3381 In this example, both of the lumping schemes are tested against the reference solution 40 gross 3379 \begin{eqnarray} \label{LUMPING WAVE TEST} 41 u=sin(5 \pi (x_0-c*t) ) 42 \end{eqnarray} 43 over the 2D unit square. Note that $u_{,i}n_i=0$ on faces $x_1=0$ and $x_1=1$. 44 ahallam 3381 Thus, on the faces $x_0=0$ and $x_0=1$ the solution is constrained. 45 gross 3379 46 ahallam 3381 To solve this problem the explicit Verlet scheme\index{Verlet scheme} was used 47 with a constant time step size $dt$ given by 48 gross 3379 \begin{eqnarray} \label{LUMPING WAVE VALET} 49 ahallam 3381 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + dt^2 a^{(n)} 50 gross 3379 \end{eqnarray} 51 ahallam 3381 for all $n=2,3,\ldots$ where the upper index ${(n)}$ refers to values at 52 time $t^{(n)}=t^{(n-1)}+h$ and $a^{(n)}$ is the solution of 53 gross 3379 \begin{eqnarray} \label{LUMPING WAVE VALET 2} 54 a^{(n)}=c^2 u^{(n-1)}_{,ii} \; . 55 \end{eqnarray} 56 ahallam 3381 This equation can be interpreted as a PDE for the unknown value $a^{(n)}$, 57 which must be solved at each time-step. 58 In the notation of equation~\ref{LINEARPDE.SINGLE.1} we thus set $D=1$ and 59 $X=-c^2 u^{(n-1)}_{,i}$. Furthermore, in order to maintain stability, 60 jfenwick 4345 the time step size needs to fulfill the Courantâ€“Friedrichsâ€“Lewy condition 61 ahallam 3381 (CFL condition). 62 \index{Courant condition} 63 \index{explicit scheme!Courant condition} 64 For this example, the CFL condition takes the form 65 gross 3379 \begin{eqnarray} \label{LUMPING WAVE CFL} 66 dt = f \cdot \frac{dx}{c} . 67 \end{eqnarray} 68 ahallam 3381 where $dx$ is the mesh size and $f$ is a safety factor. In this example, 69 we use $f=\frac{1}{6}$. 70 gross 3379 71 ahallam 3381 Figure~\ref{FIG LUMPING VALET A} depicts a temporal comparison between four 72 alternative solution algorithms: the exact solution; using a full mass matrix; 73 jfenwick 4345 using HRZ lumping; and row sum lumping. The domain utilised rectangular order 1 74 ahallam 3381 elements (element size is $0.01$) with observations taken at the point 75 $(\frac{1}{2},\frac{1}{2})$. 76 All four solutions appear to be identical for this example. This is not the case 77 for order $2$ elements, as illustrated in Figure~\ref{FIG LUMPING VALET B}. 78 For the order $2$ elements, the row sum lumping has become unstable. Row sum 79 lumping is unstable in this case because for order $2$ elements, a row sum can 80 result in a value of zero. HRZ lumping does not display the same problems, but 81 rather exhibits behaviour similar to the full mass matrix solution. When using 82 both the HRZ lumping and full mass matrix, the wave-front is slightly delayed 83 when compared with the analytical solution. 84 85 \begin{figure}[ht] 86 gross 3379 \centerline{\includegraphics[width=7cm]{lumping_valet_a_1}} 87 jfenwick 4345 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the acceleration formulation~\ref{LUMPING WAVE VALET 2} of the 88 gross 3379 Velet scheme with order $1$ elements, element size $dx=0.01$, and $c=1$.} 89 \label{FIG LUMPING VALET A} 90 \end{figure} 91 92 ahallam 3381 \begin{figure}[ht] 93 gross 3379 \centerline{\includegraphics[width=7cm]{lumping_valet_a_2}} 94 jfenwick 4345 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the acceleration formulation~\ref{LUMPING WAVE VALET 2} of the 95 gross 3379 Velet scheme with order $2$ elements, element size $0.01$, and $c=1$.} 96 \label{FIG LUMPING VALET B} 97 \end{figure} 98 99 ahallam 3381 \begin{figure}[ht] 100 gross 3379 \centerline{\includegraphics[width=7cm]{lumping_valet_u_1}} 101 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the displacement formulation~\ref{LUMPING WAVE VALET 3} of the 102 Velet scheme with order $1$ elements, element size $0.01$ and $c=1$.} 103 \label{FIG LUMPING VALET C} 104 \end{figure} 105 106 ahallam 3381 Alternatively, one can directly solve for $u^{(n)}$ by inserting 107 equation~\ref{LUMPING WAVE VALET} into equation~\ref{LUMPING WAVE VALET 2}: 108 gross 3379 \begin{eqnarray} \label{LUMPING WAVE VALET 3} 109 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + (dt\cdot c)^2 u^{(n-1)}_{,ii} \; . 110 \end{eqnarray} 111 ahallam 3381 This can also be interpreted as a PDE that must be solved at each time-step, but 112 for the unknown $u^{(n)}$. 113 As per equation~\ref{LINEARPDE.SINGLE.1} we set the general form coefficients to: 114 $D=1$; $Y=2u^{(n-1)}-u^{(n-2)}$; and $X=-(h\cdot c)^2 u^{(n-1)}_{,i}$. 115 For the full mass matrix, the acceleration ~\ref{LUMPING WAVE VALET 2} and displacement formulations ~\ref{LUMPING WAVE VALET 3} 116 are identical. 117 gross 3379 118 ahallam 3381 The displacement solution is depicted in Figure~\ref{FIG LUMPING VALET C}. The 119 ahallam 3384 domain utilised order $1$ elements (for order $2$, both 120 ahallam 3381 lumping methods are unstable). The solutions for the exact and the full mass 121 matrix approximation are almost identical while the lumping solutions, whilst 122 jfenwick 4345 identical to each other, exhibit a considerably faster wave-front propagation 123 ahallam 3381 and a decaying amplitude. 124 125 \subsection{Advection equation} 126 Consider now, a second example that demonstrates the advection equation 127 gross 3379 \begin{eqnarray} \label{LUMPING ADVECTIVE} 128 u_{,t}=(v_i u)_i \; . 129 \end{eqnarray} 130 ahallam 3381 where $v_i$ is a given velocity field. To simplify this example, set $v_i=(1,0)$ and 131 gross 3379 \begin{equation} \label{LUMPING ADVECTIVE TEST} 132 u(x,t)= 133 \left\{ 134 \begin{array}{cl} 135 1 & x_0 < t \\ 136 0 & x_0 \ge t 137 \end{array} 138 jfenwick 3382 \right\}. 139 gross 3379 \end{equation} 140 ahallam 3381 The solution scheme implemented, is the two-step Taylor-Galerkin scheme 141 \index{Taylor-Galerkin scheme} 142 gross 3379 (which is in this case equivalent to SUPG\index{SUPG}): 143 ahallam 3381 the incremental formulation is given as 144 gross 3379 \begin{eqnarray} \label{LUMPING SUPG 1} 145 du^{(n-\frac{1}{2})} = \frac{dt}{2} (v_i u^{(n-1)})_i \\ 146 du^{(n)} = dt (v_i (u^{(n-1)}+du^{(n-\frac{1}{2})}) )_i \\ 147 u^{(n)} = u^{(n)} + du^{(n)} 148 \end{eqnarray} 149 ahallam 3381 This can be reformulated to calculate $u^{(n)}$ directly: 150 gross 3379 \begin{eqnarray} \label{LUMPING SUPG 2} 151 u^{(n-\frac{1}{2})} = u^{(n-1)} + \frac{dt}{2} (v_i u^{(n-1)})_i \\ 152 u^{(n)} = u^{(n-1)} + dt (v_i u^{(n-\frac{1}{2})} )_i 153 \end{eqnarray} 154 ahallam 3381 In some cases it may be possible to combine the two equations to calculate 155 $u^{(n)}$ without the intermediate step. This approach is not discussed, because 156 it is inflexible when a greater number of terms (e.g. a diffusion term) are 157 added to the right hand side. 158 159 The advection problem is thus similar to the wave propagation problem, because 160 the time step also needs to satisfy the CFL condition 161 \index{Courant condition}\index{explicit scheme!Courant condition}. For the 162 advection problem, this takes the form 163 gross 3379 \begin{eqnarray} \label{LUMPING ADVECTION CFL} 164 dt = f \cdot \frac{dx}{\|v\|} . 165 \end{eqnarray} 166 ahallam 3381 where $dx$ is the mesh size and $f$ is a safty factor. 167 For this example, we again use $f=\frac{1}{6}$. 168 gross 3379 169 ahallam 3381 Figures~\ref{FIG LUMPING SUPG INC A} and~\ref{FIG LUMPING SUPG INC B} illustrate 170 the four incremental formulation solutions: the true solution; the exact mass matrix; 171 the HRZ lumping; and the row sum lumping. Observe, that for the order $1$ elements 172 case, there is little deviation from the exact solution before the wave front, 173 whilst there is a significant degree of osciallation after the wave-front has 174 passed. For the order $2$ elements example, all of the numerical techniques fail. 175 176 \begin{figure}[ht] 177 gross 3379 \centerline{\includegraphics[width=7cm]{lumping_SUPG_du_1}} 178 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the incremental formulation~\ref{LUMPING SUPG 1} of the 179 Taylor-Galerkin scheme with order $1$ elements, element size $dx=0.01$, $v=(1,0)$.} 180 \label{FIG LUMPING SUPG INC A} 181 \end{figure} 182 183 ahallam 3381 \begin{figure}[ht] 184 gross 3379 \centerline{\includegraphics[width=7cm]{lumping_SUPG_du_2}} 185 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the incremental formulation~\ref{LUMPING SUPG 1} of the 186 Taylor-Galerkin scheme with order $2$ elements, element size $0.01$, $v=(1,0)$.} 187 \label{FIG LUMPING SUPG INC B} 188 \end{figure} 189 190 ahallam 3381 Figure~\ref{FIG LUMPING SUPG A} depicts the results from the direct formulation 191 of the advection problem for an order $1$ mesh. Generally, the results have 192 improved when compared with the incremental formulation. The full mass matrix 193 still introduces some osciallation both before and after the arrival of the 194 wave-front at the observation point. The two lumping solutions are identical, and 195 have introduced additional smoothing to the solution. There are no oscillatory 196 effects when using lumping for this example. In Figure~\ref{FIG LUMPING SUPG Ab} 197 the mesh or element size has been reduced from 0.01 to 0.002 units. As predicted 198 by the CFL condition, this significantly improves the results when lumping is 199 applied. However, when utilising the full mass matrix, a smaller mesh size will 200 result in post wave-front oscilations which are higher frequency and slower to 201 decay. 202 gross 3379 203 ahallam 3381 Figure~\ref{FIG LUMPING SUPG B} illustrates the results when utilising elements 204 of order $2$. The full mass matrix and HRZ lumping formulations are unable to 205 correctly model the exact solution. Only the row sum lumping was capable of 206 producing a smooth and sensical result. 207 208 \begin{figure}[ht] 209 gross 3379 \centerline{\includegraphics[width=7cm]{lumping_SUPG_u_1}} 210 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the direct formulation~\ref{LUMPING SUPG 2} of the 211 Taylor-Galerkin scheme using order $1$ elements, element size $dx=0.01$, $v=(1,0)$.} 212 \label{FIG LUMPING SUPG A} 213 \end{figure} 214 ahallam 3381 215 \begin{figure}[ht] 216 gross 3379 \centerline{\includegraphics[width=7cm]{lumping_SUPG_u_1b}} 217 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the direct formulation~\ref{LUMPING SUPG 2} of the 218 Taylor-Galerkin scheme using order $1$ elements, element size $dx=0.002$, $v=(1,0)$.} 219 \label{FIG LUMPING SUPG Ab} 220 \end{figure} 221 222 ahallam 3381 \begin{figure}[ht] 223 gross 3379 \centerline{\includegraphics[width=7cm]{lumping_SUPG_u_2}} 224 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the direct formulation~\ref{LUMPING SUPG 2} of the 225 Taylor-Galerkin scheme using order $2$ elements, element size $0.01$, $v=(1,0)$.} 226 \label{FIG LUMPING SUPG B} 227 \end{figure} 228 229 gross 4052 \subsection{Summary} 230 ahallam 3381 The examples in this section have demonstrated the capabilities and limitations 231 of both HRZ and row sum lumping with comparisons to the exact and full mass 232 matrix solutions. Wave propagation type problems that utilise lumping, produce 233 results simular the full mass matrix at significantly 234 jfenwick 4345 lower computation cost. An acceleration based formulation, with HRZ lumping 235 should be implemented for such problems, and can be applied to both order $1$ and 236 ahallam 3381 order $2$ elements. 237 gross 3379 238 ahallam 3381 In transport type problems, it is essential that row sum lumping is used to 239 achieve a smooth solution. Additionally, it is not recommended that second order 240 elements be used in advection type problems. 241 gross 3379 242 243

 ViewVC Help Powered by ViewVC 1.1.26