PCA - An introduction

Introduction 

The origin of PCA is just analogous to TLS (Total Least Squares ) ,used for regression but it is extended to data compression (reducing dimensionality of data), data imputation (imputation of missing data),Iterative PCA (IPCA) (Used for finding the errors in the data )  , Non-Negative Matrix Factorisation  (NMF) (Used for extracting clean data from noisy data - The famous Cock tail party problem), Network Component Analysis (NCA) (Used generally in Bio-Informatics) , finally  Kernel PCA (KPCA) (to model any physical data). PCA is such a powerful tool , which can be used to do all this provided no consideration is taken on the computational power , that is assuming we have enough computational power we can device an algorithm for all the above mentioned .

Introduction to SVD 

PCA uses an mathematical tool in Matlab called a Singular Valued Decomposition (SVD) . This is very much analogous to any other factorisation technique such as the QR , Eigen Value decomposition ,etc . 

Let us define a matrix A , its SVD is given by :

$A_{m*n}  = U_{m*n}S_{n*m}V_{n*m}^{'}$ 

The Properties of the matrix :  The matrix U , V are the left singular ,right singular vectors which are orthonormal to each other , S is a diagonal matrix consisting of singular values (Generic : first few singular values are non zero arranged in descending order, followed by zero singular values). 

(ie $s_{1} \geq s_{2} \geq s_{3} \cdots  s_{r} \geq s_{r+1} \approx 0$ .(Here singular values after the rth value are approximately equal to 0).

It can now be easily seen that :

$AA^{'} = (USV^{'} )*(VS^{'}U^{'} ) $

$\implies AA^{'} = U(SS^{'})U^{'}$

Hence it can be seen that $ U $, $SS^{'}$ is set of eigen vectors ,eigen values respectively of the matrix $AA^{'}$  . It can also be noted that if  $A$ is a data matrix that contains each data points as a column vector (ie the matrix sizing is Variable x No of datapoints ) , then the matrix  $AA^{'}$ kind of acts like a covariance matrix across the variables . 

Similarly ,

$A^{'}A = (VS^{'}U^{'} )*(USV^{'} ) $

$\implies A^{'}A = V(S^{'}S)V^{'}$

Hence it can be seen that $ V $, $S^{'}S$ is set of eigen vectors ,eigen values respectively of the matrix $A^{'}A$  . It can also be noted that if  $A$ is a data matrix that contains each data points as a column vector (ie the matrix sizing is Variable x No of datapoints ) , then the matrix  $AA^{'}$ kind of acts like a covariance matrix across the datapoints. 

Now it can also be seen that 

$A = \sum _{i = 1} ^{k} U_{i}s_{i}V_{i}^{'}$

Where $ U_{i},V_{i}$ is the ith column of the U matrix , V matrix respectively ,$s_{i}$ represents the ith non zero singular value  ,where $ k = min(m,n) $ 

Now that adequate details have been given about the SVD , the properties of the decomposed matrices  , I hope you can yourself formulate the method that is required to perform the decomposition (ie SVD) .

Algorithm 

A setup of the problem :

For all purposes let us define some basic conventions that will be used through out the series of blogs under MVDA , 

Let A  - $A_{m*n}$be the data matrix (Variable x No of Datapoints ) .

Then matrix $Z = A - \bar{A}$ , here $\bar{A}$ is the mean of the data matrix .

for all purposes we shall now be dealing with Z matrix. 

Now some inferences we shall make on Z matrix :

First let us assume the Z matrix to be a matrix with no errors , the system following an exact linear relation then the number of constraints  = m - row rank of matrix (Z).  

Algorithm Approach  :

Then , let us perform an SVD on Z ,

then if we take the eigen vector corresponding ($U_{m}$) to the last singular value gives we get  :
Note : when $m \leq n$ (a well posed data matrix) then $k = m$ ,hence : 

$ U_{k}^{'}  Z = U_{k}^{'}* (\sum _{i = 1} ^{k} U_{i}s_{i}V_{i}^{'})$

$\implies U_{k}^{'}  Z =  s_{k}V_{k}^{'} \approx  0$. (pre-multiply with the vector)

Now you may see how the TLS algorithm discussed in my previous blog even originated .

If the data matrix is ill posed $m \geq n$ then $m \geq k$ hence

$ U_{m}^{'}  Z = U_{m}^{'}* (\sum _{i = 1} ^{k} U_{i}s_{i}V_{i}^{'})$

$\implies U_{m}^{'}  Z =  0$. (pre-multiply with the vector)

Now if there are multiple constraints then we obviously take the last few vectors of the U matrix that can now be called a regression set . 

Now as mentioned that the Z matrix consists of each and every data entry in the form of a column , if the data matrix that you are using is the other way around the you can post multiply with the last few vectors of $V$ corresponding to least singular values .

Matlab Function

Since we have seen a brief intro to SVD  , I shall now introduce to you the matlab function that would do the decomposition for you :

Matlab Function : [U,S,V] = svd(A).

The above function can always be used to do the SVD  , to obtain a regression set.
Since it was already mentioned that $k = min(m,n) $, lets consider a case where $m \leq n$ , then the extra vectors of the $V$ matrix (ie columns of the V matrix from $k+1$ to $n$ columns) are unnecessarily computed , thereby increasing computational power.

Hence we can use a specialized version of the same function under certain conditions:

[U,S,V] = svd(A,'econ')

This version basically ignores the above mentioned columns of the $V$ matrix  , also returns $S$ as a square matrix.
Note : Assuming blog convention of data matrix :
           The econ version is only used when $m \leq n$ .
            If $m \geq n$ then econ shouldn't be used.

Lid Driven Cavity - Post 3

Code Files$^1$

1. Main.m (The file to be run)
2. streamfunc.m (Solves the stream vorticity function)
3. velocity.m (Solves the velocity function)
4. omega.m (Iterates the vorticity transport equation)
5. BC.m (Implements the Boundary conditions)
6. Upwind.m (Implements the upwind scheme)
7. rmse psi.m (Computes Errors)

Terminologies used in code$^2$

1. w = $\omega$(Vorticity)
2. psi = $\psi$(StreamFunction)
3. u = (U −V elocity)
4. v = (V −V elocity)
5. u0 = $u_0$(Lid Velocity)
6. psi_1 = $c_1$
7. gamma = $\gamma$ (Kinematic viscosity)
8. dx = ∆x
9. dy = ∆y
10. alpha = Relaxation parameter for stream vorticity function
11. alpha1 = Relaxation parameter for vorticity transport equation
12. x = (X −length) of cavity
13. y = (Y −length) of cavity
14. Nx = $N_x$(No of grid points along x-direction
15. Ny = $N_y$(No of grid points along y-direction
16. err - Percentage error in $\psi$
17. err1 - Percentage error in $\omega$
18. ERR1 - Error in stream vorticity equation
19. ERR2 - Error in vorticity transport equation

-------------------------------------------------------------------------------------------------------------------------
$^1$ Link to code given under section Matlab Code URL
$^2$ Variable Names used in the Matlab code implementing the Stream Vorticity approach for the Modified Lid Driven Cavity Problem

Lid Driven Cavity - Post 2

(You can find the PDF version of the below blogpost here. )
(Citation : R, AAKASH. "Lid Driven Cavity." Classroom fundaes. Aakash R, 2 May 2017. Web. http://instafundae.blogspot.in/2017/05/lid-driven-cavity-post-2.html.)

Abstract

The lid driven cavity is a standard benchmark problem for any CFD software/Code, Lid driven cavity flows are important in many industrial processing applications such as short-dwell and flexible blade coaters.This problem consists of a 2-D cavity , whose top side moves with  a uniform velocity $u_0$.There are numerous results form the literature for the above problem , which have been used to benchmark the code,its respective results presented in this paper.The benchmarking references used in this paper are A. AbdelMigid [1] ,O. BOTELLA [3],Ghia $et$ $al.$[4]
This paper presents steady state solutions to Modified Lid Driven Cavity(Same as lid driven cavity problem but with inflow and outflow at the left and right hand side boundary) introduced in the book by Hoffman [5]  for low to Middle  $Re$s which can be extended to produce results for the Lid Driven Cavity . It is to be duly noted that the results presented in this paper are just a small part of what the Code (The Matlab code written to solve the Modified LDC) has been tested for,the same code can be used for middle to high level $Re$s too with appropriate grid sizes.
The approach used in solving the above mentioned problem is the standard stream vorticity approach . The numerical schemes used for producing results in this include central ,forward ,backward differencing schemes,Guass seidel PSOR(Point Successive Over Relaxation) methods in solving equation,etc.

The schematic diagram of the Modified Lid Driven Cavity is as follows :

Governing Equations

The governing equations of the stream vorticity approach are Vorticity transport equation,stream function  equation , continuity equations which are coupled .The below equations can be easily derived .
\begin{eqnarray}
\frac{\delta \omega}{\delta t} + u \frac{\delta \omega}{\delta x} +v\frac{\delta \omega}{\delta y} &=& \gamma \left(\frac{\delta ^2 \omega}{\delta x ^2} + \frac{\delta ^2 \omega}{\delta y ^2}\right)\\
  \nabla ^2\psi&=& - \omega\\
u &=&  \frac{\delta \psi}{\delta y}\\
v &=& -\frac{\delta \psi}{\delta x}
\end{eqnarray}
Note the continuity equation is implicitly used in deriving the stream vorticity formulations.

Boundary Conditions
Modified Lid Driven Cavity

Boundary conditions for $u,v,\psi$

Bottom side :
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_1
\end{eqnarray}
Top side :
\begin{eqnarray}
u &=& 1\\
v &=& 0\\
\psi &=& c_2
\end{eqnarray}
Right side:
For $1:(j_3 -1)$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_1
\end{eqnarray}
 For $j_3:j_4$
\begin{eqnarray}
\frac{\delta \psi}{\delta x} &=& 0
\end{eqnarray}
 For $(j4+1):end$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_2
\end{eqnarray}

Left side :
 For $1:(j_1 -1)$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_1
\end{eqnarray}
 For $j_1:j_2$
\begin{eqnarray}
\frac{\delta \psi}{\delta x} &=& 0
\end{eqnarray}
 For $(j2+1):end$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_2
\end{eqnarray}

Boundary conditions for $\omega$

\begin{equation}
  \omega  = -\left(  \frac{\delta ^2 \psi}{\delta x^2} + \frac{\delta ^2 \psi}{\delta y^2}\right)
\end{equation}

Discretized Equations

Vorticity Transport Equation

Inner Grid Formulation
Central Differencing Scheme
\begin{equation}
u(i,j)\frac{\omega(i+1,j)-\omega(i-1,j) }{2\Delta x}+v(i,j)\frac{\omega(i,j+1)-\omega(i,j-1) }{2\Delta y} =  \gamma \left( \frac{\omega(i+1,j)-2\omega(i,j)+\omega(i-1,j)}{(\Delta x)^2} + \frac{\omega(i,j+1)-2\omega(i,j)+\omega(i,j-1)}{(\Delta y)^2}\right)
\end{equation}
Rearranging , equating $\Delta x = \Delta y$,we get:
\begin{equation}
\omega(i,j) = \frac{1}{4} \left( \omega(i+1,j)\left(1 - \frac{u(i,j)\Delta x}{2\gamma}\right) + \omega(i-1,j)\left(1 + \frac{u(i,j)\Delta x }{2\gamma}\right) \\ + \omega(i,j+1)\left(1 - \frac{v(i,j)\Delta x }{2\gamma}\right) +  \omega(i,j-1)\left(1 + \frac{v(i,j)\Delta x}{2\gamma}\right)\right)
\end{equation}
Finally applying the PSOR
\begin{equation}
\omega(i,j) =(1-\alpha)\omega(i,j) +  \frac{\alpha}{4} \left( \omega(i+1,j)\left(1 - \frac{u(i,j)\Delta x}{2\gamma}\right) + \omega(i-1,j)\left(1 + \frac{u(i,j)\Delta x }{2\gamma}\right) \\ + \omega(i,j+1)\left(1 - \frac{v(i,j)\Delta x }{2\gamma}\right) +  \omega(i,j-1)\left(1 + \frac{v(i,j)\Delta x}{2\gamma}\right)\right)
\end{equation}
Boundary Formulation

Left Boundary 
\begin{eqnarray}
\omega(1,j )&=& - \left. \frac{\delta ^2 \psi}{\delta x^2} \right \vert _{(1,j)} - \left.\frac{\delta u}{\delta y}\right \vert _{(1,j)}\\
\implies \omega(1,j) &=& 2\frac{\psi(1,j)-\psi(2,j)}{\Delta x^2} - 2\frac{v(1,j)}{\Delta x} -\frac{u(1,j+1)-u(1,j-1)}{(2\Delta y)}
\end{eqnarray}
Right Boundary 
\begin{eqnarray}
\omega(im,j )&=& - \left. \frac{\delta ^2 \psi}{\delta x^2} \right \vert _{(1,j)} - \left.\frac{\delta u}{\delta y}\right \vert _{(1,j)}\\
\implies \omega(1,j) &=& 2\frac{\psi(1,j)-\psi(2,j)}{\Delta x^2} - 2\frac{v(1,j)}{\Delta x} -\frac{u(1,j+1)-u(1,j-1)}{(2\Delta y)}
\end{eqnarray}

Stream Vorticity Equation

Inner Grid Formulation
Central Differencing Scheme
\begin{eqnarray}
\omega = - \left( \frac{\psi(i+1,j)-2\omega(i,j)+\psi(i-1,j)}{(\Delta x)^2} + \frac{\psi(i,j+1)-2\psi(i,j)+\psi(i,j-1)}{(\Delta y)^2}\right)\\
\implies \psi(i,j) = \frac{1}{4}\left(\psi(i,j-1)+\psi(i,j+1)+\psi(i-1,j)+\psi(i+1,j)+ \omega (i,j)\Delta x ^2 \right)
\end{eqnarray}
Finally applying PSOR we get
\begin{equation}
 \psi(i,j) = (1-\alpha)\psi(i,j)+\frac{\alpha}{4}(\psi(i,j-1)+\psi(i,j+1)+\psi(i-1,j)+\psi(i+1,j)+ \omega (i,j)\Delta x ^2 )
\end{equation}

Boundary Formulation

Left Boundary ($j_1 : j_2$):
\begin{eqnarray}
\left. \frac{\delta \psi}{\delta x} \right \vert _{(1,j)} = v
\end{eqnarray}
Using 2nd order accurate scheme we get
\begin{eqnarray}
\psi(1,j) = \frac{4\psi(2,j)-\psi(3,j) + 2\Delta x v(1,j)}{3}
\end{eqnarray}
 Right Boundary ($j_3 : j_4$):
Using 2nd order accurate scheme we get
\begin{eqnarray}
\psi(im,j) = \frac{4\psi(im-1,j)-\psi(im-2,j) - 2\Delta x v(1,j)}{3}
\end{eqnarray}

Velocity Equations

Inner Grid Formulation

Central Differencing Scheme
\begin{eqnarray}
u = \frac{\psi(i,j+1)-\psi(i,j-1)}{(2\Delta  y)}\\
v = \frac{\psi(i+1,j)-\psi(i-1,j)}{(2\Delta  x)}
\end{eqnarray}

Boundary Formulation

Forward Differencing Scheme
\begin{eqnarray}
u &=& \frac{(\psi(j+1,i) - \psi(j,i))}{\Delta y}\\
v &=& 0
\end{eqnarray}

Parameters used for convergence measure
Before we see the algorithm used in solving the problem ,I would introduce few parameters used in checking the convergence of the solution.
The parameters are as follows :

1.  $e_1$: RMSE $\left ( \frac{\psi  ^{k+1} _{i,j} -\psi  ^k _{i,j}}{\psi  ^k _{i,j}} \right )$ (where k is the iteration number)
   Variable name used in the program ("err").
2.  $e_2$ :RMSE$\left ( \frac{\omega  ^{k+1} _{i,j} -\omega  ^k _{i,j}}{\omega  ^k _{i,j}}\right )$ .
   Variable name used in the program ("err1")
3.  $e_3$ : Numerical solution accuracy of Vorticity equation .
   Variable name used in the program ("ERR1")
4. $e_4$ : Numerical solution accuracy of Navier stokes equation .
   Variable name used in the program ("ERR2")

It is to be noted that $e_1<e_2<e_3<e_4$ for the algorithm used in this paper , hence convergence is measured using $e_4$ , divergence measured using $e_1$.

Algorithm

1. Assuming a specific grid size $N_x,N_y $  .
2.  Construct matrices for  $\psi,\omega,u,v$  using the above mentioned size using random values .  
3. Assume appropriate boundary conditions for  $\psi,\omega,u,v$ .
4. Solve for  $\psi$  using the PSOR method using Stream Vorticity function .
5.  Compute  $u,v$  using central difference scheme .
6. Iterate the PSOR method for  $\omega$   using vorticity transport Equation once (using the pre solved  $\psi,u,v $ .
   (Note : Don't solve for  $\omega$ run the PSOR over the entire grid once.)
7. Compute necessary Error Functions, check for convergence.
8. If converged , STOP. Else go to Step 4.

Results

Presented below are few results which can be used for used for benchmarking purposes . 

129x129 Grid 

S.No	U Velocity  (m/s)					V velocity  (m/s)
Distance (m)	Re 100	Re 400	Re 800	Re 1000	Re 3200	Re 100	Re 400	Re 800	Re 1000	Re 3200
1.00	1.0000	1.0000	1.0000	1.0000	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000
0.9760	0.8431	0.7568	0.6799	0.6537	0.5189	-0.0465	-0.0885	-0.1349	-0.1655	-0.3331
0.9687	0.7911	0.6825	0.5955	0.569	0.4747	-0.0621	-0.1224	-0.1885	-0.2311	-0.4441
0.9609	0.7402	0.6155	0.5283	0.5068	0.4609	-0.0777	-0.1577	-0.2434	-0.2967	-0.5218
0.9531	0.6906	0.5567	0.4775	0.4632	0.4583	-0.0930	-0.1936	-0.2974	-0.3586	-0.5537
0.8515	0.2357	0.2883	0.3211	0.3350	0.3465	-0.2388	-0.4452	-0.4241	-0.4061	-0.3679
0.7344	0.0037	0.1603	0.1806	0.1892	0.2000	-0.2134	-0.2522	-0.2295	-0.2356	-0.2300
0.6172	-0.1388	0.0198	0.0483	0.0551	0.0775	-0.0738	-0.0900	-0.0981	-0.0999	-0.1057
0.5000	-0.2087	-0.1159	-0.0705	-0.0620	-0.0367	0.0574	0.0526	0.0295	0.0258	0.0143
0.4531	-0.2134	-0.1722	-0.1161	-0.1089	-0.0800	0.0970	0.1110	0.0798	0.0771	0.0605
0.2813	-0.1570	-0.3244	-0.2982	-0.2808	-0.2405	0.1752	0.2855	0.2758	0.2690	0.2361
0.1719	-0.1013	-0.2376	-0.3564	-0.3852	-0.3434	0.1696	0.2862	0.3515	0.3715	0.3556
0.1016	-0.0641	-0.1418	-0.2456	-0.2923	-0.4297	0.1323	0.2357	0.3070	0.3370	0.4291
0.0703	-0.0464	-0.1001	-0.1783	-0.2236	-0.4036	0.1031	0.1947	0.2625	0.2958	0.4071
0.0625	-0.0418	-0.0897	-0.1610	-0.2029	-0.3798	0.0944	0.1812	0.2475	0.2804	0.3914
0.0547	-0.0371	-0.0792	-0.1435	-0.1817	-0.3499	0.0850	0.1661	0.2302	0.2625	0.3719
0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000

Comparison of my code with Ghia et al (129x129 Grid) 

Grid Points	U-Velocity		Grid Points	V-Velocity
	Present	Ghia et al		Present	Ghia et al
	Re 100			Re 100
129	1	1	129	0	0
126	0.84311	0.84123	125	-0.0621	-0.05906
125	0.79114	0.78871	124	-0.07766	-0.07391
124	0.74015	0.73722	123	-0.09305	-0.08864
123	0.69055	0.68717	122	-0.10817	-0.10313
110	0.23571	0.23151	117	-0.17657	-0.16914
95	0.00375	0.00332	111	-0.23289	-0.22445
80	-0.1388	-0.13641	104	-0.2526	-0.24533
65	-0.2087	-0.20581	65	0.057388	0.05454
59	-0.2134	-0.2109	31	0.178879	0.17527
37	-0.157	-0.15662	30	0.178671	0.17507
23	-0.1013	-0.1015	21	0.164132	0.16077
14	-0.0641	-0.06434	13	0.125833	0.12317
10	-0.0464	-0.04775	11	0.11128	0.1089
9	-0.0418	-0.04192	10	0.103134	0.10091
8	-0.0371	-0.03717	9	0.094379	0.09233
1	0	0	1	0	0

251x251 Grid

S.No	U Velocity  (m/s)				V velocity  (m/s)
Distance (m)	Re 100	Re 400	Re 1000	Re 3200	Re 100	Re 400	Re 1000	Re 3200
1.00	1.0000	1.0000	1.0000	1.0000	0.0000	0.0000	0.0000	0.0000
0.9760	0.8431	0.7568	0.6537	0.5189	-0.0465	-0.0885	-0.1655	-0.3331
0.9687	0.7911	0.6825	0.569	0.4747	-0.0621	-0.1224	-0.2311	-0.4441
0.9609	0.7402	0.6155	0.5068	0.4609	-0.0777	-0.1577	-0.2967	-0.5218
0.9531	0.6906	0.5567	0.4632	0.4583	-0.0930	-0.1936	-0.3586	-0.5537
0.8515	0.2357	0.2883	0.3350	0.3465	-0.2388	-0.4452	-0.4061	-0.3679
0.7344	0.0037	0.1603	0.1892	0.2000	-0.2134	-0.2522	-0.2356	-0.2300
0.6172	-0.1388	0.0198	0.0551	0.0775	-0.0738	-0.0900	-0.0999	-0.1057
0.5000	-0.2087	-0.1159	-0.0620	-0.0367	0.0574	0.0526	0.0258	0.0143
0.4531	-0.2134	-0.1722	-0.1089	-0.0800	0.0970	0.1110	0.0771	0.0605
0.2813	-0.1570	-0.3244	-0.2808	-0.2405	0.1752	0.2855	0.2690	0.2361
0.1719	-0.1013	-0.2376	-0.3852	-0.3434	0.1696	0.2862	0.3715	0.3556
0.1016	-0.0641	-0.1418	-0.2923	-0.4297	0.1323	0.2357	0.3370	0.4291
0.0703	-0.0464	-0.1001	-0.2236	-0.4036	0.1031	0.1947	0.2958	0.4071
0.0625	-0.0418	-0.0897	-0.2029	-0.3798	0.0944	0.1812	0.2804	0.3914
0.0547	-0.0371	-0.0792	-0.1817	-0.3499	0.0850	0.1661	0.2625	0.3719
0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000

Comparison with Ghia etal ,and O .BOTELLA:

Grid Points	U-Velocity					Grid Points	V-Velocity
	Present	Ghia et al	Present	Ghia et al	O.BOTELLA		Present	Ghia et al	Present	Ghia et al	O.BOTELLA
	Re 400		Re 1000				Re 400		Re 1000
251	1	1	1	1	1	251	0	0	0	0	0
245	0.75389422	0.75837	0.653659036	0.65928	0.6644227	243	-0.127997501	-0.12146	-0.2311	-0.21388	-0.22792
243	0.67893746	0.68439	0.569419488	0.57492	0.5808359	241	-0.16480199	-0.15663	-0.2967	-0.27669	-0.29369
241	0.61153269	0.61756	0.506788795	0.51117	0.5169277	239	-0.202194664	-0.19254	-0.3586	-0.33714	-0.35532
239	0.5527114	0.55892	0.463245781	0.46604	0.4723329	237	-0.239373072	-0.22847	-0.4134	-0.39188	-0.41038
214	0.29147429	0.29093	0.334964485	0.33304	0.3372212	228	-0.382661909	-0.23827	-0.5226	-0.5155	-0.52644
185	0.16391904	0.16256	0.189150105	0.18719	0.1886747	216	-0.452317228	-0.44993	-0.426	-0.42665	-0.42645
155	0.01930732	0.02135	0.055104972	0.05702	0.0570178	202	-0.383530806	-0.38598	-0.3172	-0.31966	-0.32021
126	-0.1152734	-0.11477	-0.061988317	-0.0608	-0.062056	126	0.052196894	0.05188	0.02584	0.02526	0.0258
114	-0.1730285	-0.17119	-0.108910482	-0.10648	-0.1082	60	0.302110802	0.30174	0.32157	0.32235	0.32536
71	-0.3276167	-0.32726	-0.280799838	-0.27805	-0.28037	58	0.30257526	0.30203	0.33037	0.33075	0.33399
44	-0.2423208	-0.24299	-0.385150406	-0.38289	-0.388569	40	0.282024218	0.28124	0.37335	0.37095	0.37692
26	-0.1428726	-0.14612	-0.292346467	-0.2973	-0.300456	24	0.228366656	0.22965	0.32684	0.32627	0.33304
19	-0.1047669	-0.10338	-0.223618183	-0.2222	-0.222896	21	0.212590596	0.2092	0.30931	0.30353	0.30991
17	-0.0938794	-0.09266	-0.202882065	-0.20196	-0.20233	19	0.200473862	0.19713	0.29582	0.29012	0.29627
15	-0.0829239	-0.08186	-0.181688508	-0.18109	-0.181288	17	0.1868535	0.1836	0.28043	0.27485	0.28071
1	0	0	0	0	0	1	0	0	0	0	0

501x501 Grid:

Note : The Underlined values show some discrepancies with the results of Ghia etal.

Grid Points	U-Velocity			Grid Points	V-Velocity
	Present	Ghia et al	AbdelMigid		Present	Ghia et al	AbdelMigid
	Re 3200	Re 3200	Re 3200		Re 3200	Re 3200	Re 3200
501	1	1	1	501	0	0	0
489	0.518868	0.53236	0.536	485	-0.44412	-0.39017	−0.4106
485	0.474687	0.48296	0.4886	481	-0.52184	-0.47425	−0.5049
481	0.460908	0.46547	0.4661	478	-0.55368	-0.52357	−0.5560
478	0.458339	0.46101	0.461	474	-0.56399	-0.54053	−0.5671
427	0.346469	0.34682	0.3491	454	-0.44061	-0.44307	−0.4462
368	0.200026	0.19791	0.2035	431	-0.37699	-0.37401	−0.3799
310	0.077457	0.07156	0.078	403	-0.31086	-0.31184	−0.3148
251	-0.03674	-0.04272	−0.0369	251	0.014275	0.00999	0.0142
228	-0.07997	-0.86636	−0.0808	118	0.286802	0.28188	0.2864
142	-0.24047	-0.24427	−0.2414	114	0.295331	0.2903	0.2954
87	-0.34339	-0.34323	−0.3445	79	0.37531	0.37119	0.3753
52	-0.42972	-0.41933	−0.4314	48	0.4294	0.42768	0.4318
36	-0.40356	-0.37827	−0.4081	40	0.418656	0.41906	0.4224
32	-0.37977	-0.35344	−0.3839	36	0.407058	0.40917	0.4109
28	-0.34989	-0.32407	−0.3530	32	0.391374	0.3956	0.3949
1	0	0	0	1	0	0	0

Matlab Code URL

Link : https://github.com/RAAKASH/Intro-to-CFD-/tree/master/LDC
Check LDC - Post 3 for description of the Matlab code .

References

[1] Tamer A. AbdelMigid, Khalid M. Saqr, Mohamed A. Kotb, and Ahmed A. Aboelfarag. Revisiting the lid-driven cavity flow problem: Review and new steady state benchmarking results using {GPU} accelerated code. Alexandria Engineering Journal, 56(1):123 – 135, 2017.

[2] C. K. Aidun, N. G. Triantafillopoulos, and J. D. Benson. Global stability of a liddriven cavity with throughflow: Flow visualization studies. Physics of Fluids A: Fluid Dynamics, 3(9):2081–2091, 1991.

[3] O Botella and R Peyret. Benchmark spectral results on the lid-driven cavity flow. Computers & Fluids, 27(4):421–433, 1998.

[4] UKNG Ghia, Kirti N Ghia, and CT Shin. High-re solutions for incompressible flow using the navier-stokes equations and a multigrid method. Journal of computational physics, 48(3):387–411, 1982.

[5] K.A. Hoffman. Computational Fluid Dynamics for Engineers. Number v. 2. Engineering Education System, 1989.

Lid Driven Cavity - Mini Project (Post 1)

Lid Driven Cavity - Post 1

The following series of posts starting from this one gives you my approach in solving the Lid Driven Cavity , the numerical methods used and my collection of results up to Reynolds Number = 3200 . 

The series of posts also gives you my results in comparison with standard results that already exist for the problem such as the Ghia Ghia and Shin.  
This post gives you an outline of the approach used in solving the problem.
Before discussing the methods involved ,let me give you a brief introduction to what the problem actually is about. This problem is a standard benchmarking problem for testing any new software,code,etc. This consists of a rectangular or square cavity (2D or 3D) ,for now consider it to be a 2D square cavity whose top side is of infinite length driven with a constant velocity.

Image Courtesy : http://www.nacad.ufrj.br/~rnelias/gallery/cavity.html

Numerical Approach 

This problem has been solved using using the Stream Vorticity approach.

•  The Navier Stokes Equation has been solved using  Point Successive Over  Relaxation method.
• The Vorticity function has also been solved using the PSOR method.
• The Velocity functions have been solved by straight forward Central differencing scheme . 

  Algorithm used for Solving

1. Assuming a specific grid size  $N_x,N_y$  (best if they are equal).
2. Construct matrices for $\psi , \omega ,u,v$ using the above mentioned size using random values (best if values are zeros.)  
3. Assume appropriate boundary conditions for $\psi , \omega ,u,v$ . 
4. Solve for $\psi$ using the PSOR method using Vorticity function .
5.  Compute $u,v$ using central difference scheme .
6. Iterate the PSOR method for $\omega$  using Navier Stokes Equation (using the presolved $\psi,u,v$ .
   (Note : Don't solve for $\omega$ just run the PSOR over the entire grid just once.) 
7. Compute necessary Error Functions. Check for convergence.
8. If converged , STOP. Else go to Step 4.  

Note : In the above algorithm I haven't Explicitly mentioned what those Error functions are,there are many methods to check convergence . But I have mainly used four Parameters ,namely:

• Error 1 :RMSE $\left ( \frac{\psi  ^{k+1} _{i,j} -\psi  ^k _{i,j}}{\psi  ^k _{i,j}} \right )$ (where k is the iteration number)
  Variable name used in my code ("err").
• Error 2 :RMSE$\left ( \frac{\omega  ^{k+1} _{i,j} -\omega  ^k _{i,j}}{\omega  ^k _{i,j}}\right )$ .
  Variable name used in my code ("err1")
• Error 3 : Numerical solution accuracy of Vorticity equation .
  Variable name used in my code ("ERR1")
• Error 4 : Numerical solution accuracy of Navier stokes equation .
  Variable name used in my code ("ERR2")

It is necessary to note that out of these error parameters , Error 4 is the dominant parameter , this parameter would be used in monitoring the accuracy of the numerical scheme as well as convergence.

The following post would consist of the derivation of the stream vorticity functions , the math of the numerical schemes deployed. For now as a motivation for you to continue reading my blog I give you a beautiful graph of the vorticity contour solution for the $R_e = 600$.

Note: I still haven't introduced as of now the numerical schemes employed in Elliptic,hyperbolic PDEs (check timeline of blog),  which I soon hope to,but as of now I would explicitly mention the Numerical schemes used in the code.

1-D Parabolic Equation PDE schemes

1-D Heat conduction equation (with no internal heat generation)

$\frac{\delta T}{\delta t} = \alpha \left ( \frac{\delta ^2 T}{\delta x^2} \right )$

Before we dwell any further in the discretization of the above PDE we need to know 3 important terms.

1. Consistency
2. Stability 
3. Convergence

Here a finite difference scheme is said to be consistent when as $\Delta x,\Delta t \rightarrow 0$ , then the scheme  $\rightarrow$ to the original PDE.

A scheme is said to be stable when the ratio of the error propagation from one time step to another is less than $1$.

A scheme is said to converge when the above to conditions are met as well as the solution is close to the actual solution (analytic or numerical) of the PDE .

  

Note : Implicit schemes are mostly unconditionally stable ,while the explicit schemes are conditionally stable or unconditionally unstable. 

Explicit schemes

FTCS scheme

$\frac{T_{i} ^{n+1} -T_{i} ^{n}}{\Delta t} = \frac{T_{i+1} ^{n} -2T_{i} ^{n}+T_{i-1} ^{n}}{(\Delta x)^2}  + O(\Delta t,(\Delta x)^2)$

This scheme is conditionally stable , consistent ,convergent ,second order accurate in space,1st order accurate in time.
For more information on the stability criterion please check the tutorial 1 under the CFD tab.

CTCS scheme

$\frac{T_{i} ^{n+1} -T_{i} ^{n-1}}{2 \Delta t} = \frac{T_{i+1} ^{n} -2T_{i} ^{n}+T_{i-1} ^{n}}{(\Delta x)^2}  + O((\Delta t)^2,(\Delta x)^2)$

This scheme is unconditionally unstable. Hence this scheme is absolutely useless .Therefore the Dufort -Frankel scheme emerged.

Dufort Frankel scheme

$\frac{T_{i} ^{n+1} -T_{i} ^{n-1}}{2 \Delta t} = \frac{T_{i+1} ^{n} -2\left (\frac{T_{i} ^{n+1} +T_{i} ^{n-1}}{2}\right )+T_{i-1} ^{n}}{(\Delta x)^2}  + O((\Delta t)^2,(\Delta x)^2,(\frac{\Delta t}{\Delta x})^2)$

This scheme is unconditionally stable. Hence this scheme is brilliant . The error term can be easily derived , do this as an exercise.

Implicit schemes
BTCS scheme (laasonen scheme)

$\frac{T_{i} ^{n+1} -T_{i} ^{n}}{\Delta t} = \frac{T_{i+1} ^{n+1} -2T_{i} ^{n+1}+T_{i-1} ^{n+1}}{(\Delta x)^2}  + O(\Delta t,(\Delta x)^2)$

Unconditionally stable scheme . This scheme forms a tridiagonal matrix , which can be solved using matrix inversion ,or the Gauss elimination method or the TDMA (thomas) algorithm. For more information on the TDMA algorithm check tutorial 2 under the CFD tab .

Crank Nicolson scheme 

$\frac{T_{i} ^{n+1} -T_{i} ^{n}}{\Delta t} =  \frac{1}{2}\left (\frac{T_{i+1} ^{n+1} -2T_{i} ^{n+1}+T_{i-1} ^{n+1}}{(\Delta x)^2} +\frac{T_{i+1} ^{n} -2T_{i} ^{n}+T_{i-1} ^{n}}{(\Delta x)^2}    \right ) + O((\Delta t )^2,(\Delta x)^2)$

This is a really interesting scheme centered around the time step $n + 1/2$ . As can be seen the above scheme is second order accurate in time as well.
Unconditionally stable scheme . This scheme forms a tridiagonal matrix , which can be solved using matrix inversion ,or the Gauss elimination method or the TDMA (thomas) algorithm. 

Beta Formulation

Finally comes the $\beta $ formulation , this combines all the above mentioned schemes in a single shot .

$\frac{T_{i} ^{n+1} -T_{i} ^{n}}{\Delta t} =  \beta \left (\frac{T_{i+1} ^{n+1} -2T_{i} ^{n+1}+T_{i-1} ^{n+1}}{(\Delta x)^2}  \right )+(1-\beta) \left ( \frac{T_{i+1} ^{n} -2T_{i} ^{n}+T_{i-1} ^{n}}{(\Delta x)^2}   \right )  $

This scheme behaves like an FTCS scheme for $\beta < 1/2$
This scheme behaves like an BTCS scheme for $\beta > 1/2$
This scheme behaves like an Crank Nicolson scheme for $\beta = 1/2$ 

2D Parabolic PDEs

Here let us now consider a plate whose initial conditions are specified .
FTCS scheme
So Let us discretize the governing 2 -Dimensional Parabolic PDE using FTCS scheme .
$$\frac{T_{i,j}^{n+1}-T_{i,j}^{n}}{\Delta t} = \alpha \left (\frac{T_{i+1,j}^{n}-2T_{i,j}^{n}+T_{i-1,j}^{n}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n}-2T_{i,j}^{n}+T_{i,j-1}^{n}}{\Delta y ^2} \right ) $$
Since this is an explicit scheme it becomes easy to calculate the solutions of $T$ as a function of time.(ie)

$$T_{i,j}^{n+1} = T_{i,j}^{n}+ \Delta t \alpha \left (\frac{T_{i+1,j}^{n}-2T_{i,j}^{n}+T_{i-1,j}^{n}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n}-2T_{i,j}^{n}+T_{i,j-1}^{n}}{\Delta y ^2} \right )$$
For computational ease $\Delta x = \Delta y$ can be used. Now let us see the stability criterion of the scheme below.

Stability Condition

It can be derived to show that the stability for the 2 D heat conduction equation is similar to that of the 1-D heat equation ,hence 

$$ \alpha \Delta (t) \left (\frac{1}{(\Delta x )^2} + \frac{1}{(\Delta x )^2}  \right  )\leq \frac{1}{2}$$

Check Computational Fluid Dynamics by  Klaus A. Hoffmann (Author), Steve T. Chiang (Author) for reference.

This can be further reduced to 

$$ \alpha \Delta (t) \left  (\frac{1}{(\Delta x )^2} \right )\leq \frac{1}{4}$$ when $\Delta x = \Delta y$   

BTCS scheme

The same can be done in BTCS scheme to obtain 

$$\frac{T_{i,j}^{n+1}-T_{i,j}^{n}}{\Delta t} = \alpha \left (\frac{T_{i+1,j}^{n+1}-2T_{i,j}^{n+1}+T_{i-1,j}^{n+1}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n+1}-2T_{i,j}^{n+1}+T_{i,j-1}^{n+1}}{\Delta y ^2} \right ) $$

The above equation as can be seen that it forms a pentadiagonal matrix which can be solved by matrix inversion which is computationally expensive ,hence this scheme wont be used most of the time. Hence emerges the ADI scheme

ADI scheme
The ADI (alternate direction implicit method) is the execution of the implicit schemes at one direction at time steps  $n,n +1/2$
$$\frac{T_{i,j}^{n+1/2}-T_{i,j}^{n}}{\Delta t} = \alpha \left (\frac{T_{i+1,j}^{n+1/2}-2T_{i,j}^{n+1/2}+T_{i-1,j}^{n+1/2}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n}-2T_{i,j}^{n}+T_{i,j-1}^{n}}{\Delta y ^2} \right ) $$
$$\frac{T_{i,j}^{n+1}-T_{i,j}^{n+1/2}}{\Delta t} = \alpha \left (\frac{T_{i+1,j}^{n+1}-2T_{i,j}^{n+1}+T_{i-1,j}^{n+1}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n+1/2}-2T_{i,j}^{n+1/2}+T_{i,j-1}^{n+1/2}}{\Delta y ^2} \right ) $$

The above equations as can be seen that it forms a tridiagonal matrix which can easily be solved by TDMA algorithm . The first equation is solved for finding the solutions for intermediate time step $n +1/2$ then this is used to find the solution at time step $n+1$ .

CFD - Finite Difference schemes (Intro)

First Order Derivatives

Forward Difference

$\frac{\delta u}{\delta x} =  \frac{u_{i+1} - u_{i}}{\Delta x} + O(\Delta x)$

Derivation :
  $u_{i+1}=u_{i} + \left .\frac{\Delta x}{1 !} \frac{\delta u}{\delta x} \right \vert _{i} + \left .\frac{(\Delta x)^2}{2 !} \frac{\delta ^2u}{\delta x ^2} \right \vert _{i} + \left .\frac{(\Delta x)^3}{3 !} \frac{\delta ^3u}{\delta x ^3} \right \vert _{i} +O((\Delta x)^4) $
$\implies \frac{u_{i+1}-u_{i}}{\Delta x} = \left .\frac{1}{1 !} \frac{\delta u}{\delta x} \right \vert _{i} + \left .\frac{\Delta x}{2 !} \frac{\delta ^2u}{\delta x ^2} \right \vert _{i} + \left .\frac{(\Delta x)^2}{3 !} \frac{\delta ^3u}{\delta x ^3} \right \vert _{i} +O((\Delta x)^3) $
$\implies \frac{u_{i+1}-u_{i}}{\Delta x} = \left .\frac{1}{1 !} \frac{\delta u}{\delta x} \right \vert _{i} + O((\Delta x)) $

Backward Difference

$\frac{\delta u}{\delta x} =  \frac{u_{i} - u_{i-1}}{\Delta x} + O(\Delta x)$

Derivation :

  $u_{i-1}=u_{i} - \left .\frac{\Delta x}{1 !} \frac{\delta u}{\delta x} \right \vert _{i} + \left .\frac{(\Delta x)^2}{2 !} \frac{\delta ^2u}{\delta x ^2} \right \vert _{i} - \left .\frac{(\Delta x)^3}{3 !} \frac{\delta ^3u}{\delta x ^3} \right \vert _{i} +O((\Delta x)^4) $
$\implies \frac{u_{i}-u_{i-1}}{\Delta x} = \left .\frac{1}{1 !} \frac{\delta u}{\delta x} \right \vert _{i} + O((\Delta x)) $

Central Difference

$\frac{\delta u}{\delta x} =  \frac{u_{i+1} - u_{i-1}}{ 2\Delta x } + O((\Delta x)^2)$

Derivation :

$-u_{i+1}+u_{i} + \left .\frac{\Delta x}{1 !} \frac{\delta u}{\delta x} \right \vert _{i} + \left .\frac{(\Delta x)^2}{2 !} \frac{\delta ^2u}{\delta x ^2} \right \vert _{i} + \left .\frac{(\Delta x)^3}{3 !} \frac{\delta ^3u}{\delta x ^3} \right \vert _{i} +O((\Delta x)^4)  =0      -(1) $ 

  $-u_{i-1}+u_{i} - \left .\frac{\Delta x}{1 !} \frac{\delta u}{\delta x} \right \vert _{i} + \left .\frac{(\Delta x)^2}{2 !} \frac{\delta ^2u}{\delta x ^2} \right \vert _{i} - \left .\frac{(\Delta x)^3}{3 !} \frac{\delta ^3u}{\delta x ^3} \right \vert _{i} +O((\Delta x)^4) =0-(2)$ 

$(\alpha eqn_1 +\beta eqn_2) = ..... + 0*\left .\frac{(\Delta x)^2}{2 !} \frac{\delta ^2u}{\delta x ^2} \right \vert _{i}  +O((\Delta x)^3)  =0 -(3)$

Solving for $\alpha,\beta $ substituting in equation 3 , rearranging we get the above result. 

Similarly using the above method, we can obtain higher order accuracies for  forward, backward difference schemes as well.

Second Order Derivatives 

Forward Difference

$\frac{\delta^2 u}{\delta x^2} =  \frac{u_{i+2}-2u_{i+1} + u_{i}}{ \Delta x ^2} + O((\Delta x))$

Central Difference

$\frac{\delta^2 u}{\delta x^2} =  \frac{u_{i}-2u_{i-1} + u_{i-2}}{ \Delta x ^2} + O((\Delta x))$

Central Difference

$\frac{\delta^2 u}{\delta x^2} =  \frac{u_{i+1}-2u_{i} + u_{i-1}}{ \Delta x ^2} + O((\Delta x)^2)$

Note : Derive the above equations for your own practice.

Other methods of derivation

There are many methods using which the above schemes can be derived . One is that of the polynomial method ,  the other one is discussed for a simple case below :

 $$\frac{\delta ^n f}{\delta x ^n} = \frac{\delta ^{n-1} }{\delta x ^{n-1}} \left (\frac{\delta f}{\delta x }\right )    $$

Hence for a second order derivative (using forward difference) we get

$$\frac{\delta ^2 f}{\delta x ^2} = \frac{\delta  }{\delta x  } \left (\frac{f_{i+1} - f_{i}}{\Delta x} \right )   $$

$$\frac{\delta ^2 f}{\delta x ^2} = \left (\frac{\delta f_{i+1}  }{\delta x}-\frac{\delta f_{i}  }{\delta x}\right )\frac{1}{\Delta x}   $$

$$\frac{\delta ^2 f}{\delta x ^2} = \left (\frac{ f_{i+1}  -f_{i+1}  }{\Delta x}-\frac{ f_{i+1} - f_{i}  }{\Delta x}\right )\frac{1}{\Delta x}   $$

$$\frac{\delta ^2 f}{\delta x ^2} = \left (\frac{ f_{i+1}  - 2f_{i+1} +f_{i} }{\Delta x ^2} \right )$$

Polynomial method

Let there be a polynomial $f(x)   = ax^2+bx+c$ .

$\frac{\delta f}{\delta x} = 2ax+b    -(1)$

$\frac{\delta ^2 f}{\delta x ^2} = 2a  -(2)$

Let at $x =0 ,f(x) = f_{i},x = \Delta x ,f(x) = f_{i+1} ,x =2 \Delta x ,f(x) = f_{i+2} $.

Now solve for $a , b,c$ in terms of $f_{i},f_{i+1},f_{i+2},\Delta x$ 

Substituting $a , b$ back in equation 1, 2 will give you the appropriate derivatives.

Try the above calculations for practice.