/[escript]/branches/doubleplusgood/doc/user/lumping.tex
ViewVC logotype

Contents of /branches/doubleplusgood/doc/user/lumping.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4345 - (show annotations)
Fri Mar 29 07:09:41 2013 UTC (5 years, 11 months ago) by jfenwick
File MIME type: application/x-tex
File size: 11840 byte(s)
Spelling fixes
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2013 by University of Queensland
4 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Open Software License version 3.0
8 % http://www.opensource.org/licenses/osl-3.0.php
9 %
10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 % Development since 2012 by School of Earth Sciences
12 %
13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14
15
16 \section{Some Remarks on Lumping}
17 \label{LUMPING}
18
19 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 \index{linear solver!row sum lumping}\index{row sum lumping} and HRZ
31 lumping \index{linear solver!HRZ lumping}\index{HRZ lumping}.
32
33 \subsection{Scalar wave equation}
34 One example where lumping can be applied to a hyperpolic problem, is
35 the scalar wave equation
36 \begin{eqnarray} \label{LUMPING WAVE}
37 u_{,tt}=c^2 u_{,ii} \; .
38 \end{eqnarray}
39 In this example, both of the lumping schemes are tested against the reference solution
40 \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 Thus, on the faces $x_0=0$ and $x_0=1$ the solution is constrained.
45
46 To solve this problem the explicit Verlet scheme\index{Verlet scheme} was used
47 with a constant time step size $dt$ given by
48 \begin{eqnarray} \label{LUMPING WAVE VALET}
49 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + dt^2 a^{(n)}
50 \end{eqnarray}
51 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 \begin{eqnarray} \label{LUMPING WAVE VALET 2}
54 a^{(n)}=c^2 u^{(n-1)}_{,ii} \; .
55 \end{eqnarray}
56 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 the time step size needs to fulfill the Courant–Friedrichs–Lewy condition
61 (CFL condition).
62 \index{Courant condition}
63 \index{explicit scheme!Courant condition}
64 For this example, the CFL condition takes the form
65 \begin{eqnarray} \label{LUMPING WAVE CFL}
66 dt = f \cdot \frac{dx}{c} .
67 \end{eqnarray}
68 where $dx$ is the mesh size and $f$ is a safety factor. In this example,
69 we use $f=\frac{1}{6}$.
70
71 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 using HRZ lumping; and row sum lumping. The domain utilised rectangular order 1
74 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 \centerline{\includegraphics[width=7cm]{lumping_valet_a_1}}
87 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the acceleration formulation~\ref{LUMPING WAVE VALET 2} of the
88 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 \begin{figure}[ht]
93 \centerline{\includegraphics[width=7cm]{lumping_valet_a_2}}
94 \caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the acceleration formulation~\ref{LUMPING WAVE VALET 2} of the
95 Velet scheme with order $2$ elements, element size $0.01$, and $c=1$.}
96 \label{FIG LUMPING VALET B}
97 \end{figure}
98
99 \begin{figure}[ht]
100 \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 Alternatively, one can directly solve for $u^{(n)}$ by inserting
107 equation~\ref{LUMPING WAVE VALET} into equation~\ref{LUMPING WAVE VALET 2}:
108 \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 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
118 The displacement solution is depicted in Figure~\ref{FIG LUMPING VALET C}. The
119 domain utilised order $1$ elements (for order $2$, both
120 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 identical to each other, exhibit a considerably faster wave-front propagation
123 and a decaying amplitude.
124
125 \subsection{Advection equation}
126 Consider now, a second example that demonstrates the advection equation
127 \begin{eqnarray} \label{LUMPING ADVECTIVE}
128 u_{,t}=(v_i u)_i \; .
129 \end{eqnarray}
130 where $v_i$ is a given velocity field. To simplify this example, set $v_i=(1,0)$ and
131 \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 \right\}.
139 \end{equation}
140 The solution scheme implemented, is the two-step Taylor-Galerkin scheme
141 \index{Taylor-Galerkin scheme}
142 (which is in this case equivalent to SUPG\index{SUPG}):
143 the incremental formulation is given as
144 \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 This can be reformulated to calculate $u^{(n)}$ directly:
150 \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 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 \begin{eqnarray} \label{LUMPING ADVECTION CFL}
164 dt = f \cdot \frac{dx}{\|v\|} .
165 \end{eqnarray}
166 where $dx$ is the mesh size and $f$ is a safty factor.
167 For this example, we again use $f=\frac{1}{6}$.
168
169 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 \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 \begin{figure}[ht]
184 \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 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
203 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 \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
215 \begin{figure}[ht]
216 \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 \begin{figure}[ht]
223 \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 \subsection{Summary}
230 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 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 order $2$ elements.
237
238 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
242
243

  ViewVC Help
Powered by ViewVC 1.1.26