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

Revision 3379 - (show annotations)
Wed Nov 24 04:48:49 2010 UTC (8 years, 9 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 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 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