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 
When using explicit time integration schemes (we will discuss two examples later in this section) 
19 
very small time steps need to be used in order to maintain stability. In orer to reduce 
20 
computational costs at each time step lumping of the coefficient matrix to a diagonal matrix 
21 
is appled reducing costs for updating the solution to a simple componentbycomponent 
22 
vectorvector product. Although lumping can speedup simulations dramatically 
23 
it needs to be applied with some care. In this section we will briefly discuss 
24 
two commonly used techiques, namely row sum lumping 
25 
\index{linear solver!row sum lumping}\index{row sum lumping} and HRZ 
26 
lumping \index{linear solver!HRZ lumping}\index{HRZ lumping}. 
27 

28 
As an example for an hyperpolic problem we take a look 
29 
at the scalar wave equation 
30 
\begin{eqnarray} \label{LUMPING WAVE} 
31 
u_{,tt}=c^2 u_{,ii} \; . 
32 
\end{eqnarray} 
33 
The schemes are tested against the reference solution 
34 
\begin{eqnarray} \label{LUMPING WAVE TEST} 
35 
u=sin(5 \pi (x_0c*t) ) 
36 
\end{eqnarray} 
37 
over the 2D unit square. Note that $u_{,i}n_i=0$ on faces $x_1=0$ and $x_1=1$. 
38 
On faces $x_0=0$ and $x_0=1$ the solution is constraints. 
39 

40 
The solution scheme is explicit Verlet scheme\index{Verlet scheme} with constant time step size $dt$ given by 
41 
\begin{eqnarray} \label{LUMPING WAVE VALET} 
42 
u^{(n)}=2u^{(n1)}u^{(n2)} + dt^2 a^{(n)} \\ 
43 
\end{eqnarray} 
44 
for all $n=2,3,\ldots$ where the upper index ${(n)}$ refers to values at time $t^{(n)}=t^{(n1)}+h$ and 
45 
$a^{(n)}$ is the solution of 
46 
\begin{eqnarray} \label{LUMPING WAVE VALET 2} 
47 
a^{(n)}=c^2 u^{(n1)}_{,ii} \; . 
48 
\end{eqnarray} 
49 
The equation is read as a PDE for the unknown $a^{(n)}$ which needs to be solved in each timestep. 
50 
In the notation of equation~\ref{LINEARPDE.SINGLE.1} we need to set $D=1$ and 
51 
$X=c^2 u^{(n1)}_{,i}$. In order to maintain stability the time step size needs to fullfill 
52 
Courant–Friedrichs–Lewy condition (CFL condition)\index{Courant condition}\index{explicit scheme!Courant condition} which 
53 
in the case discussed here takes the form 
54 
\begin{eqnarray} \label{LUMPING WAVE CFL} 
55 
dt = f \cdot \frac{dx}{c} . 
56 
\end{eqnarray} 
57 
where $dx$ is the mesh size and $f$ is a safty factor. Here we use $f=\frac{1}{6}$. 
58 

59 
\begin{figure} 
60 
\centerline{\includegraphics[width=7cm]{lumping_valet_a_1}} 
61 
\caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the accelaraton formulation~\ref{LUMPING WAVE VALET 2} of the 
62 
Velet scheme with order $1$ elements, element size $dx=0.01$, and $c=1$.} 
63 
\label{FIG LUMPING VALET A} 
64 
\end{figure} 
65 

66 
\begin{figure} 
67 
\centerline{\includegraphics[width=7cm]{lumping_valet_a_2}} 
68 
\caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the accelaraton formulation~\ref{LUMPING WAVE VALET 2} of the 
69 
Velet scheme with order $2$ elements, element size $0.01$, and $c=1$.} 
70 
\label{FIG LUMPING VALET B} 
71 
\end{figure} 
72 

73 
\begin{figure} 
74 
\centerline{\includegraphics[width=7cm]{lumping_valet_u_1}} 
75 
\caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the displacement formulation~\ref{LUMPING WAVE VALET 3} of the 
76 
Velet scheme with order $1$ elements, element size $0.01$ and $c=1$.} 
77 
\label{FIG LUMPING VALET C} 
78 
\end{figure} 
79 

80 
Figure~\ref{FIG LUMPING VALET A} shows a comparison between 
81 
the exact solution, using a full mass matrix, using HRZ lumping and row sum lumping 
82 
for rectangular, order 1 elements (element size is $0.01$) at the point $(\frac{1}{2},\frac{1}{2})$ 
83 
over time. All four graphs are allmost identical. As shown in 
84 
Figure~\ref{FIG LUMPING VALET B} the situation is changing when using order $2$ elements. 
85 
In fact one can observe that the row sum lumping is unstable (Notice that for 
86 
row sum lumping for order $2$ elements zero row sums can be produced) while using HRZ lumping 
87 
shows similar behaviour like using using the full mass matrix with the a 
88 
slightly slower wave propagation speed. 
89 

90 
Alternatively, one can directly solve for $u^{(n)}$ by inserting equation~(\ref{LUMPING WAVE VALET}) 
91 
into equations~\ref{LUMPING WAVE VALET 2}: 
92 
\begin{eqnarray} \label{LUMPING WAVE VALET 3} 
93 
u^{(n)}=2u^{(n1)}u^{(n2)} + (dt\cdot c)^2 u^{(n1)}_{,ii} \; . 
94 
\end{eqnarray} 
95 
Again we can read this update as a PDE for the unknown $u^{(n)}$ which needs to be solved in each timestep. 
96 
In the notation of equation~\ref{LINEARPDE.SINGLE.1} we need to set $D=1$, $Y=2u^{(n1)}u^{(n2)}$ and 
97 
$X=(h\cdot c)^2 u^{(n1)}_{,i}$. Notice that for the case of using the full mass matrix 
98 
the accelaraton formulation~\ref{LUMPING WAVE VALET 2} and displacement formulation~\ref{LUMPING WAVE VALET 3} 
99 
are identical. Figure~\ref{FIG LUMPING VALET C} shows the result for order $1$ elements (for order $2$ both 
100 
lumping methods are not stable). The curves for the exact and the full mass matrix approximation 
101 
are almost identical while the  in this case identical  curves for row sum and 
102 
HRZ lumping are showing faster period with a decaying amplitude. 
103 

104 
As a second example 
105 
we look at the advection equation 
106 
\begin{eqnarray} \label{LUMPING ADVECTIVE} 
107 
u_{,t}=(v_i u)_i \; . 
108 
\end{eqnarray} 
109 
where $v_i$ is a given velocity field. For our simple test we set $v_i=(1,0)$ and 
110 
\begin{equation} \label{LUMPING ADVECTIVE TEST} 
111 
u(x,t)= 
112 
\left\{ 
113 
\begin{array}{cl} 
114 
1 & x_0 < t \\ 
115 
0 & x_0 \ge t 
116 
\end{array} 
117 
\right. 
118 
\end{equation} 
119 
As a solution scheme we use the twostep TaylorGalerkin scheme \index{TaylorGalerkin scheme} 
120 
(which is in this case equivalent to SUPG\index{SUPG}): 
121 
which in the incremental formulation is given as 
122 
\begin{eqnarray} \label{LUMPING SUPG 1} 
123 
du^{(n\frac{1}{2})} = \frac{dt}{2} (v_i u^{(n1)})_i \\ 
124 
du^{(n)} = dt (v_i (u^{(n1)}+du^{(n\frac{1}{2})}) )_i \\ 
125 
u^{(n)} = u^{(n)} + du^{(n)} 
126 
\end{eqnarray} 
127 
This can be reformulated calculating $u^{(n)}$ directly: 
128 
\begin{eqnarray} \label{LUMPING SUPG 2} 
129 
u^{(n\frac{1}{2})} = u^{(n1)} + \frac{dt}{2} (v_i u^{(n1)})_i \\ 
130 
u^{(n)} = u^{(n1)} + dt (v_i u^{(n\frac{1}{2})} )_i 
131 
\end{eqnarray} 
132 
Notice that one can insert the first equation into the second to calculate $u^{(n)}$ with out the intermediate step. We 
133 
will not discuss this approach has it does not very flexible when more terms (e.f. a diffusion) into the right hand side. 
134 
Similar to the wave propagation problem the time step size needs to fullfill a 
135 
Courant–Friedrichs–Lewy condition (CFL condition)\index{Courant condition}\index{explicit scheme!Courant condition} which 
136 
in the case of the advection problem takes the form 
137 
\begin{eqnarray} \label{LUMPING ADVECTION CFL} 
138 
dt = f \cdot \frac{dx}{\v\} . 
139 
\end{eqnarray} 
140 
where $dx$ is the mesh size and $f$ is a safty factor. Again we use $f=\frac{1}{6}$. 
141 

142 
\begin{figure} 
143 
\centerline{\includegraphics[width=7cm]{lumping_SUPG_du_1}} 
144 
\caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the incremental formulation~\ref{LUMPING SUPG 1} of the 
145 
TaylorGalerkin scheme with order $1$ elements, element size $dx=0.01$, $v=(1,0)$.} 
146 
\label{FIG LUMPING SUPG INC A} 
147 
\end{figure} 
148 

149 
\begin{figure} 
150 
\centerline{\includegraphics[width=7cm]{lumping_SUPG_du_2}} 
151 
\caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the incremental formulation~\ref{LUMPING SUPG 1} of the 
152 
TaylorGalerkin scheme with order $2$ elements, element size $0.01$, $v=(1,0)$.} 
153 
\label{FIG LUMPING SUPG INC B} 
154 
\end{figure} 
155 

156 
Figures~\ref{FIG LUMPING SUPG INC A} and~\ref{FIG LUMPING SUPG INC B} show the value of the true solution, 
157 
using the exact mass matrix, the HRZ lumping and the row sum lumping for the incermental formulation. In the case of order $1$ elements 
158 
osciallation prior to the arrivial of the front and in particular after the front has passed. 
159 
For order $2$ elements all techniques fail. 
160 

161 
\begin{figure} 
162 
\centerline{\includegraphics[width=7cm]{lumping_SUPG_u_1}} 
163 
\caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the direct formulation~\ref{LUMPING SUPG 2} of the 
164 
TaylorGalerkin scheme using order $1$ elements, element size $dx=0.01$, $v=(1,0)$.} 
165 
\label{FIG LUMPING SUPG A} 
166 
\end{figure} 
167 
\begin{figure} 
168 
\centerline{\includegraphics[width=7cm]{lumping_SUPG_u_1b}} 
169 
\caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the direct formulation~\ref{LUMPING SUPG 2} of the 
170 
TaylorGalerkin scheme using order $1$ elements, element size $dx=0.002$, $v=(1,0)$.} 
171 
\label{FIG LUMPING SUPG Ab} 
172 
\end{figure} 
173 

174 
\begin{figure} 
175 
\centerline{\includegraphics[width=7cm]{lumping_SUPG_u_2}} 
176 
\caption{Amplitude at point $(\frac{1}{2},\frac{1}{2})$ using the direct formulation~\ref{LUMPING SUPG 2} of the 
177 
TaylorGalerkin scheme using order $2$ elements, element size $0.01$, $v=(1,0)$.} 
178 
\label{FIG LUMPING SUPG B} 
179 
\end{figure} 
180 

181 
For the order $1$ case results are better when using the direct formulation. As shown in Figure~\ref{FIG LUMPING SUPG A} 
182 
using the full mass matrix introduces still osciallation around the time when the front hits the 
183 
observation point. Lumping (HRZ and row sum lumping are identical in our case) introduces additional 
184 
smoothing into the solution avoiding any rippels. As confirmed by Figure~\ref{FIG LUMPING SUPG Ab} 
185 
the approximation around the time of impact improves when a smaller mesh size and lumping is used. 
186 
However, using the full mass matrix a smaller mesh size will produces stronger riples. 
187 
Figure~\ref{FIG LUMPING SUPG B} shows that 
188 
for order $2$ elements nor HRZ lumping neither using the full mass matrix is able to produce a stable solution 
189 
while row sum lumping can smooth some of the osciallations. 
190 

191 
In summary: For wave propagation type problem lumping produces results simular to the usage of the full mass matrix as significantly 
192 
lower costs. An accelleration based formulation with HRZ lumping should be used which can be appied to order $1$ and order $2$ 
193 
element. For transport type problems using row sum lumping is essential to achieve a smooth solution. It is not recommended to 
194 
use second order elements. 
195 

196 

197 
