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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3379 - (show annotations)
Wed Nov 24 04:48:49 2010 UTC (7 years, 10 months ago) by gross
File MIME type: application/x-tex
File size: 10262 byte(s)
some clarification on lumping
1
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 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 component-by-component
22 vector-vector product. Although lumping can speed-up 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_0-c*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^{(n-1)}-u^{(n-2)} + 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^{(n-1)}+h$ and
45 $a^{(n)}$ is the solution of
46 \begin{eqnarray} \label{LUMPING WAVE VALET 2}
47 a^{(n)}=c^2 u^{(n-1)}_{,ii} \; .
48 \end{eqnarray}
49 The equation is read as a PDE for the unknown $a^{(n)}$ which needs to be solved in each time-step.
50 In the notation of equation~\ref{LINEARPDE.SINGLE.1} we need to set $D=1$ and
51 $X=-c^2 u^{(n-1)}_{,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^{(n-1)}-u^{(n-2)} + (dt\cdot c)^2 u^{(n-1)}_{,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 time-step.
96 In the notation of equation~\ref{LINEARPDE.SINGLE.1} we need to set $D=1$, $Y=2u^{(n-1)}-u^{(n-2)}$ and
97 $X=-(h\cdot c)^2 u^{(n-1)}_{,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 two-step Taylor-Galerkin scheme \index{Taylor-Galerkin 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^{(n-1)})_i \\
124 du^{(n)} = dt (v_i (u^{(n-1)}+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^{(n-1)} + \frac{dt}{2} (v_i u^{(n-1)})_i \\
130 u^{(n)} = u^{(n-1)} + 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 Taylor-Galerkin 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 Taylor-Galerkin 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 Taylor-Galerkin 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 Taylor-Galerkin 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 Taylor-Galerkin 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

  ViewVC Help
Powered by ViewVC 1.1.26