# Diff of /trunk/doc/user/lumping.tex

revision 3379 by gross, Wed Nov 24 04:48:49 2010 UTC revision 3381 by ahallam, Wed Nov 24 23:37:55 2010 UTC
# Line 15  Line 15
15  \section{Some Remarks on Lumping}  \section{Some Remarks on Lumping}
16  \label{REMARKS ON LUMPING}  \label{REMARKS ON LUMPING}
17
18  When using explicit time integration schemes (we will discuss two examples later in this section)  Explicit time integration schemes (two examples are discussed later in this
19  very small time steps need to be used in order to maintain stability. In orer to reduce  section), require very small time steps in order to maintain numerical stability.
20  computational costs at each time step lumping of the coefficient matrix to a diagonal matrix  Unfortunately, these small time increments often result in a prohibitive
21  is appled reducing costs for updating the solution to a simple component-by-component  computational cost.
22  vector-vector product. Although lumping can speed-up simulations dramatically  In order to minimise these costs, a technique termed lumping can be utilised.
23  it needs to be applied with some care. In this section we will briefly discuss  Lumping is applied to the coefficient matrix, reducing it to a simple diagonal
24  two commonly used techiques, namely row sum lumping  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  \index{linear solver!row sum lumping}\index{row sum lumping} and HRZ
30  lumping \index{linear solver!HRZ lumping}\index{HRZ lumping}.  lumping \index{linear solver!HRZ lumping}\index{HRZ lumping}.
31
32  As an example for an hyperpolic problem we take a look  \subsection{Scalar wave equation}
33  at the scalar wave equation  One example where lumping can be applied to a hyperpolic problem, is
34    the scalar wave equation
35  \begin{eqnarray} \label{LUMPING WAVE}  \begin{eqnarray} \label{LUMPING WAVE}
36  u_{,tt}=c^2 u_{,ii} \; .  u_{,tt}=c^2 u_{,ii} \; .
37  \end{eqnarray}  \end{eqnarray}
38  The schemes are tested against the reference solution  In this example, both of the lumping schemes are tested against the reference solution
39  \begin{eqnarray} \label{LUMPING WAVE TEST}  \begin{eqnarray} \label{LUMPING WAVE TEST}
40  u=sin(5 \pi (x_0-c*t) )  u=sin(5 \pi (x_0-c*t) )
41  \end{eqnarray}  \end{eqnarray}
42  over the 2D unit square. Note that $u_{,i}n_i=0$ on faces $x_1=0$ and $x_1=1$.  over the 2D unit square. Note that $u_{,i}n_i=0$ on faces $x_1=0$ and $x_1=1$.
43  On faces $x_0=0$ and $x_0=1$ the solution is constraints.  Thus, on the faces $x_0=0$ and $x_0=1$ the solution is constrained.
44
45  The solution scheme is explicit Verlet scheme\index{Verlet scheme} with constant time step size $dt$ given by  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}  \begin{eqnarray} \label{LUMPING WAVE VALET}
48  u^{(n)}=2u^{(n-1)}-u^{(n-2)} + dt^2 a^{(n)} \\  u^{(n)}=2u^{(n-1)}-u^{(n-2)} + dt^2 a^{(n)}
49  \end{eqnarray}  \end{eqnarray}
50  for all $n=2,3,\ldots$ where the upper index ${(n)}$ refers to values at time $t^{(n)}=t^{(n-1)}+h$ and  for all $n=2,3,\ldots$ where the upper index ${(n)}$ refers to values at
51  $a^{(n)}$ is the solution of  time $t^{(n)}=t^{(n-1)}+h$ and $a^{(n)}$ is the solution of
52  \begin{eqnarray} \label{LUMPING WAVE VALET 2}  \begin{eqnarray} \label{LUMPING WAVE VALET 2}
53  a^{(n)}=c^2 u^{(n-1)}_{,ii} \; .  a^{(n)}=c^2 u^{(n-1)}_{,ii} \; .
54  \end{eqnarray}  \end{eqnarray}
55  The equation is read as a PDE for the unknown $a^{(n)}$ which needs to be solved in each time-step.  This equation can be interpreted as a PDE for the unknown value $a^{(n)}$,
56  In the notation of equation~\ref{LINEARPDE.SINGLE.1} we need to set $D=1$ and  which must be solved at each time-step.
57  $X=-c^2 u^{(n-1)}_{,i}$. In order to maintain stability the time step size needs to fullfill  In the notation of equation~\ref{LINEARPDE.SINGLE.1} we thus set $D=1$ and
58  Courant–Friedrichs–Lewy condition (CFL condition)\index{Courant condition}\index{explicit scheme!Courant condition} which  $X=-c^2 u^{(n-1)}_{,i}$. Furthermore, in order to maintain stability,
59  in the case discussed here takes the form  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}  \begin{eqnarray} \label{LUMPING WAVE CFL}
65  dt = f \cdot \frac{dx}{c} .  dt = f \cdot \frac{dx}{c} .
66  \end{eqnarray}  \end{eqnarray}
67  where $dx$ is the mesh size and $f$ is a safty factor. Here we use $f=\frac{1}{6}$.  where $dx$ is the mesh size and $f$ is a safety factor. In this example,
68    we use $f=\frac{1}{6}$.
69
70  \begin{figure}  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.
83
84    \begin{figure}[ht]
85  \centerline{\includegraphics[width=7cm]{lumping_valet_a_1}}  \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  \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$.}  Velet scheme with order $1$ elements, element size $dx=0.01$, and $c=1$.}
88  \label{FIG LUMPING VALET A}  \label{FIG LUMPING VALET A}
89  \end{figure}  \end{figure}
90
91  \begin{figure}  \begin{figure}[ht]
92  \centerline{\includegraphics[width=7cm]{lumping_valet_a_2}}  \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  \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$.}  Velet scheme with order $2$ elements, element size $0.01$, and $c=1$.}
95  \label{FIG LUMPING VALET B}  \label{FIG LUMPING VALET B}
96  \end{figure}  \end{figure}
97
98  \begin{figure}  \begin{figure}[ht]
99  \centerline{\includegraphics[width=7cm]{lumping_valet_u_1}}  \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  \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$.}  Velet scheme with order $1$ elements, element size $0.01$ and $c=1$.}
102  \label{FIG LUMPING VALET C}  \label{FIG LUMPING VALET C}
103  \end{figure}  \end{figure}
104
105  Figure~\ref{FIG LUMPING VALET A} shows a comparison between  Alternatively, one can directly solve for $u^{(n)}$ by inserting
106  the exact solution, using a full mass matrix, using HRZ lumping and row sum lumping  equation~\ref{LUMPING WAVE VALET} into equation~\ref{LUMPING WAVE VALET 2}:
for rectangular, order 1 elements (element size is $0.01$) at the point $(\frac{1}{2},\frac{1}{2})$
over time. All four graphs are allmost identical. As shown in
Figure~\ref{FIG LUMPING VALET B} the situation is changing when using order $2$ elements.
In fact one can observe that the row sum lumping is unstable (Notice that for
row sum lumping for order $2$ elements zero row sums can be produced) while using HRZ lumping
shows similar behaviour like using using the full mass matrix with the a
slightly slower wave propagation speed.

Alternatively, one can directly solve for $u^{(n)}$ by inserting equation~(\ref{LUMPING WAVE VALET})
into equations~\ref{LUMPING WAVE VALET 2}:
107  \begin{eqnarray} \label{LUMPING WAVE VALET 3}  \begin{eqnarray} \label{LUMPING WAVE VALET 3}
108  u^{(n)}=2u^{(n-1)}-u^{(n-2)} + (dt\cdot c)^2 u^{(n-1)}_{,ii} \; .  u^{(n)}=2u^{(n-1)}-u^{(n-2)} + (dt\cdot c)^2 u^{(n-1)}_{,ii} \; .
109  \end{eqnarray}  \end{eqnarray}
110  Again we can read this update as a PDE for the unknown $u^{(n)}$ which needs to be solved in each time-step.  This can also be interpreted as a PDE that must be solved at each time-step, but
111  In the notation of equation~\ref{LINEARPDE.SINGLE.1} we need to set $D=1$, $Y=2u^{(n-1)}-u^{(n-2)}$ and  for the unknown $u^{(n)}$.
112  $X=-(h\cdot c)^2 u^{(n-1)}_{,i}$. Notice that for the case of using the full mass matrix  As per equation~\ref{LINEARPDE.SINGLE.1} we set the general form coefficients to:
113  the accelaraton formulation~\ref{LUMPING WAVE VALET 2} and displacement formulation~\ref{LUMPING WAVE VALET 3}  $D=1$; $Y=2u^{(n-1)}-u^{(n-2)}$; and $X=-(h\cdot c)^2 u^{(n-1)}_{,i}$.
114  are identical. Figure~\ref{FIG LUMPING VALET C} shows the result for order $1$ elements (for order $2$ both  For the full mass matrix, the acceleration ~\ref{LUMPING WAVE VALET 2} and displacement formulations ~\ref{LUMPING WAVE VALET 3}
115  lumping methods are not stable). The curves for the exact and the full mass matrix approximation  are identical.
116  are almost identical while the - in this case identical - curves for row sum and
117  HRZ lumping are showing faster period with a decaying amplitude.  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.
123
124  As a second example  \subsection{Advection equation}
125  we look at the advection equation  Consider now, a second example that demonstrates the advection equation
127  u_{,t}=(v_i u)_i \; .  u_{,t}=(v_i u)_i \; .
128  \end{eqnarray}  \end{eqnarray}
129  where $v_i$ is a given velocity field. For our simple test we set $v_i=(1,0)$ and  where $v_i$ is a given velocity field. To simplify this example, set $v_i=(1,0)$ and
131  u(x,t)=  u(x,t)=
132  \left\{  \left\{
# Line 114  u(x,t)= Line 134  u(x,t)=
134     1 &  x_0 < t \\     1 &  x_0 < t \\
135     0 &  x_0 \ge t     0 &  x_0 \ge t
136     \end{array}     \end{array}
137  \right.  \right}.
138
139  As a solution scheme we use the two-step Taylor-Galerkin scheme \index{Taylor-Galerkin scheme}  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}):  (which is in this case equivalent to SUPG\index{SUPG}):
142  which in the incremental formulation is given as  the incremental formulation is given as
143  \begin{eqnarray} \label{LUMPING SUPG 1}  \begin{eqnarray} \label{LUMPING SUPG 1}
144  du^{(n-\frac{1}{2})} = \frac{dt}{2} (v_i u^{(n-1)})_i \\  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 \\  du^{(n)} = dt (v_i (u^{(n-1)}+du^{(n-\frac{1}{2})}) )_i \\
146  u^{(n)} = u^{(n)} + du^{(n)}  u^{(n)} = u^{(n)} + du^{(n)}
147  \end{eqnarray}  \end{eqnarray}
148  This can be reformulated calculating $u^{(n)}$ directly:  This can be reformulated to calculate $u^{(n)}$ directly:
149  \begin{eqnarray} \label{LUMPING SUPG 2}  \begin{eqnarray} \label{LUMPING SUPG 2}
150  u^{(n-\frac{1}{2})} = u^{(n-1)} + \frac{dt}{2} (v_i u^{(n-1)})_i \\  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  u^{(n)} =  u^{(n-1)} + dt (v_i u^{(n-\frac{1}{2})} )_i
152  \end{eqnarray}  \end{eqnarray}
153  Notice that one can insert the first equation into the second to calculate $u^{(n)}$ with out the intermediate step. We  In some cases it may be possible to combine the two equations to calculate
154  will not discuss this approach has it does not very flexible when more terms (e.f. a diffusion) into the right hand side.  $u^{(n)}$ without the intermediate step. This approach is not discussed, because
155  Similar to the wave propagation problem the time step size needs to fullfill a   it is inflexible when a greater number of terms (e.g. a diffusion term) are
156  Courant–Friedrichs–Lewy condition (CFL condition)\index{Courant condition}\index{explicit scheme!Courant condition} which   added to the right hand side.
157  in the case of the advection problem takes the form
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
163  dt = f \cdot \frac{dx}{\|v\|} .  dt = f \cdot \frac{dx}{\|v\|} .
164  \end{eqnarray}  \end{eqnarray}
165  where $dx$ is the mesh size and $f$ is a safty factor. Again we use $f=\frac{1}{6}$.  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  \begin{figure}  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.
174
175    \begin{figure}[ht]
176  \centerline{\includegraphics[width=7cm]{lumping_SUPG_du_1}}  \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  \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)$.}  Taylor-Galerkin scheme with order $1$ elements, element size $dx=0.01$, $v=(1,0)$.}
179  \label{FIG LUMPING SUPG INC A}  \label{FIG LUMPING SUPG INC A}
180  \end{figure}  \end{figure}
181
182  \begin{figure}  \begin{figure}[ht]
183  \centerline{\includegraphics[width=7cm]{lumping_SUPG_du_2}}  \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  \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)$.}  Taylor-Galerkin scheme  with order $2$ elements, element size $0.01$, $v=(1,0)$.}
186  \label{FIG LUMPING SUPG INC B}  \label{FIG LUMPING SUPG INC B}
187  \end{figure}  \end{figure}
188
189  Figures~\ref{FIG LUMPING SUPG INC A} and~\ref{FIG LUMPING SUPG INC B} show the value of the true solution,  Figure~\ref{FIG LUMPING SUPG A} depicts the results from the direct formulation
190  using the exact mass matrix, the HRZ lumping and the row sum lumping for the incermental formulation. In the case of order $1$ elements  of the advection problem for an order $1$ mesh. Generally, the results have
191  osciallation prior to the arrivial of the front and -in particular- after the front has passed.  improved when compared with the incremental formulation. The full mass matrix
192  For order $2$ elements all techniques fail.  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.
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}  \begin{figure}[ht]
208  \centerline{\includegraphics[width=7cm]{lumping_SUPG_u_1}}  \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  \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)$.}  Taylor-Galerkin scheme using order $1$ elements, element size $dx=0.01$, $v=(1,0)$.}
211  \label{FIG LUMPING SUPG A}  \label{FIG LUMPING SUPG A}
212  \end{figure}  \end{figure}
213  \begin{figure}
214    \begin{figure}[ht]
215  \centerline{\includegraphics[width=7cm]{lumping_SUPG_u_1b}}  \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  \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)$.}  Taylor-Galerkin scheme using order $1$ elements, element size $dx=0.002$, $v=(1,0)$.}
218  \label{FIG LUMPING SUPG Ab}  \label{FIG LUMPING SUPG Ab}
219  \end{figure}  \end{figure}
220
221  \begin{figure}  \begin{figure}[ht]
222  \centerline{\includegraphics[width=7cm]{lumping_SUPG_u_2}}  \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  \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)$.}  Taylor-Galerkin scheme  using order $2$ elements, element size $0.01$, $v=(1,0)$.}
225  \label{FIG LUMPING SUPG B}  \label{FIG LUMPING SUPG B}
226  \end{figure}  \end{figure}
227
228  For the order $1$ case results are better when using the direct formulation. As shown in Figure~\ref{FIG LUMPING SUPG A}  \subsection{Sumamry}
229  using the full mass matrix introduces still osciallation around the time when the front hits the  The examples in this section have demonstrated the capabilities and limitations
230  observation point. Lumping (HRZ and row sum lumping are identical in our case) introduces additional  of both HRZ and row sum lumping with comparisons to the exact and full mass
231  smoothing into the solution avoiding any rippels. As confirmed by Figure~\ref{FIG LUMPING SUPG Ab}  matrix solutions. Wave propagation type problems that utilise lumping, produce
232  the approximation around the time of impact improves when a smaller mesh size and lumping is used.  results simular the full mass matrix at significantly
233  However, using the full mass matrix a smaller mesh size will produces stronger riples.  lower computation cost. An accelleration based formulation, with HRZ lumping
234  Figure~\ref{FIG LUMPING SUPG B} shows that  should be implemented for such problems, and can be appied to both order $1$ and
235  for order $2$ elements nor HRZ lumping neither using the full mass matrix is able to produce a stable solution   order $2$ elements.
236  while row sum lumping can smooth some of the osciallations.
237    In transport type problems, it is essential that row sum lumping is used to
238  In summary: For wave propagation type problem lumping produces results simular to the usage of the full mass matrix as significantly  achieve a smooth solution. Additionally, it is not recommended that second order
239  lower costs. An accelleration based formulation with HRZ lumping should be used which can be appied to order $1$ and order $2$  elements be used in advection type problems.
element. For transport type problems using row sum lumping is essential to achieve a smooth solution. It is not recommended to
use second order elements.
240
241
242

Legend:
 Removed from v.3379 changed lines Added in v.3381