ViewVC logotype

Contents of /trunk/doc/user/lumping.tex

Parent Directory Parent Directory | Revision Log Revision Log

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

  ViewVC Help
Powered by ViewVC 1.1.26