1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032013 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/osl3.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 componentbycomponent 
27 
vectorvector 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_0c*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^{(n1)}u^{(n2)} + 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^{(n1)}+h$ and $a^{(n)}$ is the solution of 
53 
\begin{eqnarray} \label{LUMPING WAVE VALET 2} 
54 
a^{(n)}=c^2 u^{(n1)}_{,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 timestep. 
58 
In the notation of equation~\ref{LINEARPDE.SINGLE.1} we thus set $D=1$ and 
59 
$X=c^2 u^{(n1)}_{,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 wavefront 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^{(n1)}u^{(n2)} + (dt\cdot c)^2 u^{(n1)}_{,ii} \; . 
110 
\end{eqnarray} 
111 
This can also be interpreted as a PDE that must be solved at each timestep, 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^{(n1)}u^{(n2)}$; and $X=(h\cdot c)^2 u^{(n1)}_{,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 wavefront 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 twostep TaylorGalerkin scheme 
141 
\index{TaylorGalerkin 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^{(n1)})_i \\ 
146 
du^{(n)} = dt (v_i (u^{(n1)}+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^{(n1)} + \frac{dt}{2} (v_i u^{(n1)})_i \\ 
152 
u^{(n)} = u^{(n1)} + 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 wavefront 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 
TaylorGalerkin 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 
TaylorGalerkin 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 
wavefront 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 wavefront 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 
TaylorGalerkin 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 
TaylorGalerkin 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 
TaylorGalerkin 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 
