Solve Pressure Correction.f90 - Solution of pressure-correction equation and correction of U,V and P

<pre>

!Sample program for solving Lid-driven cavity flow test using SIMPLE-algorithm
! solution of pressure-correction quation and correction U,V, and P
!Copyright (C) 2010  Michail Kiričkov
!Copyright (C) 2016  Michail Kiričkov, Kaunas University for Technology

!This program is free software; you can redistribute it and/or
!modify it under the terms of the GNU General Public License
!as published by the Free Software Foundation; either version 2
!of the License, or (at your option) any later version.

!This program is distributed in the hope that it will be useful,
!but WITHOUT ANY WARRANTY; without even the implied warranty of
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!GNU General Public License for more details.

!You should have received a copy of the GNU General Public License
!along with this program; if not, write to the Free Software
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

!**********************************************************************
Subroutine Solve_Pressure_Correction
include 'icomm_1.f90'
integer test
    DpU(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,1)
    DpV(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,2) 
!This is 1./AP coefficient involved in the derivation of the pressure equation 
!and correction velocities. Check for simple algorithm in dedicated books.
!AP is the central coefficent there is one for U momentum equation 
!    --> APU and one for V-momentum APV.
! in fact for this code APU(i,j) = AP(i,j,1)
! APV(i,j) = AP(i,j,2)
! because he works with vector componnents style
!*****************************************************************************   
! vertical  faces East - e
   Do 101 I=1,NXmax
   Do 101 J=1,NYmax-1
 	     DXc_e = 0.
		 S_e   = 0.
		 VOL_e = 0.
		 APU_e = 0.
		 Ul_e  = 0.
		 DPl_e = 0.   
		 DPx_e = 0.
		 U_e   = 0.
	If(i.ne.NXmax)then
 	     !this is the distance between node P and node E
         DXc_e = Xc(i+1,j+1) - Xc(i,j+1) 
         !This is the east surface S_e
         S_e   =  Y(i  ,j+1) -  Y(i,j  ) 
		 VOL_e = DXc_e * S_e 
	!-------------------------------------------------------------
        !this is the iterpolation of 1/APU at east location         
         APU_e =       0.5 * (   DpU(i,j+1)   +   DpU(i+1,j+1)   ) 
        !this is the interpolation of velocity U on east location
         Ul_e  =       0.5 * (     F(i,j+1,1) +     F(i+1,j+1,1) ) 
        ! the same for interpolation of DPx
         DPl_e =       0.5 * ( DPx_c(i,j+1)   + DPx_c(i+1,j+1)   ) 
        ! this is the gradient of pressure at east location usint central difference
         DPx_e =             ( F(i+1,j+1,4) - F(i,j+1,4) ) / DXc_e 
        ! This is the interpolated velocity using at east location 
        !         using the interpolation of Rhie and Chow 	
         U_e   = Ul_e  + APU_e * VOL_e * ( DPl_e - DPx_e) 
    !--------------------------------------------------------------    
		 Con_we(i,j) =   U_e  * S_e ! this is the convective flux.
	End If
 	     Con_we(1,j) = 0.
  	     Con_we(NXmax,j) = 0.
101 continue
!********************************************************************
! horisontal  faces
   Do 102 I=1,NXmax-1
   Do 102 J=1,NYmax
 	     DYc_n = 0.
	     S_n   = 0.
	     VOL_n = 0.
         APV_n = 0.
         Vl_n  = 0.
         DPl_n = 0.
         DPy_n = 0.
   		 V_n   = 0.
	If((j.ne.NYmax))then
 	     DYc_n = Yc(i+1,j+1) - Yc(i+1,j)
	     S_n   =  X(i+1,j  ) -  X(i  ,j)
	     VOL_n = DYc_n * S_n 
	!-----------------------------------------------------------
         APV_n =       0.5 * (   DpV(i+1,j)   +   DpV(i+1,j+1)   )
         Vl_n  =       0.5 * (     F(i+1,j,2) +     F(i+1,j+1,2) )
         DPl_n =       0.5 * ( DPy_c(i+1,j)   + DPy_c(i+1,j+1)   )
         DPy_n = ( F(i+1,j+1,4) - F(i+1,j,4) ) / DYc_n
  		 V_n   = Vl_n  + APV_n * VOL_n * ( DPl_n - DPy_n ) 	
    !-----------------------------------------------------------
		 Con_ns(i,j) =   V_n  * S_n 
    End If
         Con_ns(i,1) = 0.
 	     Con_ns(i,NYmax) = 0.
 102 continue
  Summ = 0.
	Do 200 I=2,NXmax
	Do 200 J=2,NYmax
		 S_e   =  Y(i  ,j) -  Y(i,j-1)
		 S_n   =  X(i  ,j) -    X(i-1,j)
    An(i,j) = 0.5 * (   DpV(i,j)   +   DpV(i,j+1)   )  * S_n * S_n
    As(i,j) = 0.5 * (   DpV(i,j)   +   DpV(i,j-1)   )  * S_n * S_n
    Ae(i,j) = 0.5 * (   DpU(i,j)   +   DpU(i+1,j)   )  * S_e * S_e
    Aw(i,j) = 0.5 * (   DpU(i,j)   +   DpU(i-1,j)   )  * S_e * S_e
    Ap(i,j,3) = (An(i,j) + As(i,j) + Ae(i,j) +	Aw(i,j))
    Sp(i,j,3) = -1. * ( Con_ns(i-1,j) - Con_ns(i-1,j-1) &
                      + Con_we(i,j-1) - Con_we(i-1,j-1) )
    Summ = Summ + Sp(i,j,3) ! 
 200 continue
    write(*,*)'Sum',summ
    write(*,*)'solve PP'
 niter = 0
20 continue
        call Convergence_Criteria(3,Res_sum_before,niter)
        niter= niter + 1
        Call TDMA_1(3)
		call Convergence_Criteria(3,Res_sum_After,niter)

    if((abs(Res_sum_before-Res_sum_After).Ge.1.0E-10) &
    .and.(niter.le.1500).and.(abs(Res_sum_After).ge.1.0E-06))go to 20		 
	write(*,*)'Pressure correction Res_sum_before-Res_sum_After' &
    ,Res_sum_before-Res_sum_After,niter
!  pressure correction	
	Do 301 I=2,NXmax
	Do 301 J=2,NYmax
        F(i,j,4) = F(i,j,4) + 0.2 *   F(i,j,3) 
 301 continue
!  velocities and pressure correction
	Do 302 I=2,NXmax
	Do 302 J=2,NYmax
		DY = Y(i,j)-Y(i,j-1)
		DX = X(i,j)-X(i-1,j)
		PPe = 0.5 * ( F(i,j,3) + F(i+1,j,3) )
		PPw = 0.5 * ( F(i,j,3) + F(i-1,j,3) )
		PPn = 0.5 * ( F(i,j,3) + F(i,j+1,3) )
		PPs = 0.5 * ( F(i,j,3) + F(i,j-1,3) )
		F(i,j,1) = F(i,j,1) + (PPw - PPe)  * DpU(i,j)  * Dy 
		F(i,j,2) = F(i,j,2) + (PPs - PPn)  * DpV(i,j)  * Dx 
 302 continue
!**********************************************************************
	Do 501 I=1,NXmaxC
       F(i,1,4) = F(i,2,4)            !   
       F(i,NYmaxC,4) = F(i,NYmaxC-1,4)!   
 501 continue

 	Do 502 J=1,NYmaxC
       F(1,j,4) = F(2,j,4)             !  
       F(NXmaxC,j,4) = F(NXmaxC-1,j,4) !  
 502 continue
 F(:,:,4) = F(:,:,4) - F(3,4,4) 
999 continue

Return
End

</pre>