1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 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/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

14 

15 
\section{Some Remarks on Lumping} 
16 
\label{REMARKS ON LUMPING} 
17 

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 componentbycomponent 
26 
vectorvector 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}. 
31 

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_0c*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. 
44 

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^{(n1)}u^{(n2)} + 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^{(n1)}+h$ and $a^{(n)}$ is the solution of 
52 
\begin{eqnarray} \label{LUMPING WAVE VALET 2} 
53 
a^{(n)}=c^2 u^{(n1)}_{,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 timestep. 
57 
In the notation of equation~\ref{LINEARPDE.SINGLE.1} we thus set $D=1$ and 
58 
$X=c^2 u^{(n1)}_{,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}$. 
69 

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 wavefront is slightly delayed 
82 
when compared with the analytical solution. 
83 

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$.} 
88 
\label{FIG LUMPING VALET A} 
89 
\end{figure} 
90 

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$.} 
95 
\label{FIG LUMPING VALET B} 
96 
\end{figure} 
97 

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$.} 
102 
\label{FIG LUMPING VALET C} 
103 
\end{figure} 
104 

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^{(n1)}u^{(n2)} + (dt\cdot c)^2 u^{(n1)}_{,ii} \; . 
109 
\end{eqnarray} 
110 
This can also be interpreted as a PDE that must be solved at each timestep, 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^{(n1)}u^{(n2)}$; and $X=(h\cdot c)^2 u^{(n1)}_{,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. 
116 

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 wavefront propagration 
122 
and a decaying amplitude. 
123 

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 twostep TaylorGalerkin scheme 
140 
\index{TaylorGalerkin 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^{(n1)})_i \\ 
145 
du^{(n)} = dt (v_i (u^{(n1)}+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^{(n1)} + \frac{dt}{2} (v_i u^{(n1)})_i \\ 
151 
u^{(n)} = u^{(n1)} + 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. 
157 

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}$. 
167 

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 wavefront has 
173 
passed. For the order $2$ elements example, all of the numerical techniques fail. 
174 

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 
TaylorGalerkin scheme with order $1$ elements, element size $dx=0.01$, $v=(1,0)$.} 
179 
\label{FIG LUMPING SUPG INC A} 
180 
\end{figure} 
181 

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 
TaylorGalerkin scheme with order $2$ elements, element size $0.01$, $v=(1,0)$.} 
186 
\label{FIG LUMPING SUPG INC B} 
187 
\end{figure} 
188 

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 
wavefront 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 wavefront oscilations which are higher frequency and slower to 
200 
decay. 
201 

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. 
206 

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 
TaylorGalerkin scheme using order $1$ elements, element size $dx=0.01$, $v=(1,0)$.} 
211 
\label{FIG LUMPING SUPG A} 
212 
\end{figure} 
213 

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 
TaylorGalerkin scheme using order $1$ elements, element size $dx=0.002$, $v=(1,0)$.} 
218 
\label{FIG LUMPING SUPG Ab} 
219 
\end{figure} 
220 

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 
TaylorGalerkin scheme using order $2$ elements, element size $0.01$, $v=(1,0)$.} 
225 
\label{FIG LUMPING SUPG B} 
226 
\end{figure} 
227 

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. 
236 

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. 
240 

241 

242 
