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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26