/[escript]/trunk/doc/user/lumping.tex
ViewVC logotype

Annotation of /trunk/doc/user/lumping.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3382 - (hide annotations)
Thu Nov 25 00:43:18 2010 UTC (8 years, 9 months ago) by jfenwick
File MIME type: application/x-tex
File size: 11732 byte(s)
I've added better credits for dev team.
There is a new page just before the contents page with current dev team 
and a more detailed list in the user guide appendix.

Please check this and let me know if you think there are errors or 
omissions.

Also fixed a typo.
Removed reference to python style files we no longer use from the 
license file.


1 gross 3379
2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3     %
4     % Copyright (c) 2003-2010 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/osl-3.0.php
11     %
12     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13    
14    
15     \section{Some Remarks on Lumping}
16     \label{REMARKS ON LUMPING}
17    
18 ahallam 3381 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 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 gross 3379 \index{linear solver!row sum lumping}\index{row sum lumping} and HRZ
30     lumping \index{linear solver!HRZ lumping}\index{HRZ lumping}.
31    
32 ahallam 3381 \subsection{Scalar wave equation}
33     One example where lumping can be applied to a hyperpolic problem, is
34     the scalar wave equation
35 gross 3379 \begin{eqnarray} \label{LUMPING WAVE}
36     u_{,tt}=c^2 u_{,ii} \; .
37     \end{eqnarray}
38 ahallam 3381 In this example, both of the lumping schemes are tested against the reference solution
39 gross 3379 \begin{eqnarray} \label{LUMPING WAVE TEST}
40     u=sin(5 \pi (x_0-c*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 ahallam 3381 Thus, on the faces $x_0=0$ and $x_0=1$ the solution is constrained.
44 gross 3379
45 ahallam 3381 To solve this problem the explicit Verlet scheme\index{Verlet scheme} was used
46     with a constant time step size $dt$ given by
47 gross 3379 \begin{eqnarray} \label{LUMPING WAVE VALET}
48 ahallam 3381 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + dt^2 a^{(n)}
49 gross 3379 \end{eqnarray}
50 ahallam 3381 for all $n=2,3,\ldots$ where the upper index ${(n)}$ refers to values at
51     time $t^{(n)}=t^{(n-1)}+h$ and $a^{(n)}$ is the solution of
52 gross 3379 \begin{eqnarray} \label{LUMPING WAVE VALET 2}
53     a^{(n)}=c^2 u^{(n-1)}_{,ii} \; .
54     \end{eqnarray}
55 ahallam 3381 This equation can be interpreted as a PDE for the unknown value $a^{(n)}$,
56     which must be solved at each time-step.
57     In the notation of equation~\ref{LINEARPDE.SINGLE.1} we thus set $D=1$ and
58     $X=-c^2 u^{(n-1)}_{,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 gross 3379 \begin{eqnarray} \label{LUMPING WAVE CFL}
65     dt = f \cdot \frac{dx}{c} .
66     \end{eqnarray}
67 ahallam 3381 where $dx$ is the mesh size and $f$ is a safety factor. In this example,
68     we use $f=\frac{1}{6}$.
69 gross 3379
70 ahallam 3381 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 gross 3379 \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 ahallam 3381 \begin{figure}[ht]
92 gross 3379 \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 ahallam 3381 \begin{figure}[ht]
99 gross 3379 \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 ahallam 3381 Alternatively, one can directly solve for $u^{(n)}$ by inserting
106     equation~\ref{LUMPING WAVE VALET} into equation~\ref{LUMPING WAVE VALET 2}:
107 gross 3379 \begin{eqnarray} \label{LUMPING WAVE VALET 3}
108     u^{(n)}=2u^{(n-1)}-u^{(n-2)} + (dt\cdot c)^2 u^{(n-1)}_{,ii} \; .
109     \end{eqnarray}
110 ahallam 3381 This can also be interpreted as a PDE that must be solved at each time-step, 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^{(n-1)}-u^{(n-2)}$; and $X=-(h\cdot c)^2 u^{(n-1)}_{,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 gross 3379
117 ahallam 3381 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     \subsection{Advection equation}
125     Consider now, a second example that demonstrates the advection equation
126 gross 3379 \begin{eqnarray} \label{LUMPING ADVECTIVE}
127     u_{,t}=(v_i u)_i \; .
128     \end{eqnarray}
129 ahallam 3381 where $v_i$ is a given velocity field. To simplify this example, set $v_i=(1,0)$ and
130 gross 3379 \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 jfenwick 3382 \right\}.
138 gross 3379 \end{equation}
139 ahallam 3381 The solution scheme implemented, is the two-step Taylor-Galerkin scheme
140     \index{Taylor-Galerkin scheme}
141 gross 3379 (which is in this case equivalent to SUPG\index{SUPG}):
142 ahallam 3381 the incremental formulation is given as
143 gross 3379 \begin{eqnarray} \label{LUMPING SUPG 1}
144     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 \\
146     u^{(n)} = u^{(n)} + du^{(n)}
147     \end{eqnarray}
148 ahallam 3381 This can be reformulated to calculate $u^{(n)}$ directly:
149 gross 3379 \begin{eqnarray} \label{LUMPING SUPG 2}
150     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
152     \end{eqnarray}
153 ahallam 3381 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 gross 3379 \begin{eqnarray} \label{LUMPING ADVECTION CFL}
163     dt = f \cdot \frac{dx}{\|v\|} .
164     \end{eqnarray}
165 ahallam 3381 where $dx$ is the mesh size and $f$ is a safty factor.
166     For this example, we again use $f=\frac{1}{6}$.
167 gross 3379
168 ahallam 3381 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 gross 3379 \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     Taylor-Galerkin 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 ahallam 3381 \begin{figure}[ht]
183 gross 3379 \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     Taylor-Galerkin scheme with order $2$ elements, element size $0.01$, $v=(1,0)$.}
186     \label{FIG LUMPING SUPG INC B}
187     \end{figure}
188    
189 ahallam 3381 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     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 gross 3379
202 ahallam 3381 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 gross 3379 \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     Taylor-Galerkin scheme using order $1$ elements, element size $dx=0.01$, $v=(1,0)$.}
211     \label{FIG LUMPING SUPG A}
212     \end{figure}
213 ahallam 3381
214     \begin{figure}[ht]
215 gross 3379 \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     Taylor-Galerkin scheme using order $1$ elements, element size $dx=0.002$, $v=(1,0)$.}
218     \label{FIG LUMPING SUPG Ab}
219     \end{figure}
220    
221 ahallam 3381 \begin{figure}[ht]
222 gross 3379 \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     Taylor-Galerkin scheme using order $2$ elements, element size $0.01$, $v=(1,0)$.}
225     \label{FIG LUMPING SUPG B}
226     \end{figure}
227    
228 ahallam 3381 \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 gross 3379
237 ahallam 3381 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 gross 3379
241    
242    

  ViewVC Help
Powered by ViewVC 1.1.26