CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > Coeff 11.f90 - calculate the coefficients

Coeff 11.f90 - calculate the coefficients

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
(Added QUICK implementation)
 
(5 intermediate revisions not shown)
Line 1: Line 1:
 +
<pre>
 +
 +
!Sample program for solving Smith-Hutton Test using different schemes
 +
!of covective terms approximation - Coefficients computing modul
 +
!Copyright (C) 2005  Michail Kiričkov
 +
 +
!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 Coef_1(nf)
Subroutine Coef_1(nf)
-
 
include 'icomm_1.f90'   
include 'icomm_1.f90'   
Dimension F_out(nx,ny),CheckFlux(nx,ny)
Dimension F_out(nx,ny),CheckFlux(nx,ny)
-
 
Character  Filename*10
Character  Filename*10
-
 
!  calculation of fluxes
!  calculation of fluxes
-
 
!  all geometry has rectangular 2D notation
!  all geometry has rectangular 2D notation
-
 
Do 100 I= 2,NXmax
Do 100 I= 2,NXmax
-
 
Do 100 J= 2,NYmax
Do 100 J= 2,NYmax
-
 
-
 
Gam_e = ( Gam(i+1,j  ) + Gam(i  ,j  ) ) * 0.5
Gam_e = ( Gam(i+1,j  ) + Gam(i  ,j  ) ) * 0.5
-
 
Gam_w = ( Gam(i-1,j  ) + Gam(i  ,j  ) ) * 0.5
Gam_w = ( Gam(i-1,j  ) + Gam(i  ,j  ) ) * 0.5
-
 
Gam_s = ( Gam(i  ,j-1) + Gam(i  ,j  ) ) * 0.5
Gam_s = ( Gam(i  ,j-1) + Gam(i  ,j  ) ) * 0.5
-
 
Gam_n = ( Gam(i  ,j+1) + Gam(i  ,j  ) ) * 0.5
Gam_n = ( Gam(i  ,j+1) + Gam(i  ,j  ) ) * 0.5
-
 
-
 
Ro_e = ( Ro(i+1,j  ) + Ro(i  ,j  ) ) * 0.5
Ro_e = ( Ro(i+1,j  ) + Ro(i  ,j  ) ) * 0.5
-
 
Ro_w = ( Ro(i-1,j  ) + Ro(i  ,j  ) ) * 0.5
Ro_w = ( Ro(i-1,j  ) + Ro(i  ,j  ) ) * 0.5
-
 
Ro_s = ( Ro(i  ,j-1) + Ro(i  ,j  ) ) * 0.5
Ro_s = ( Ro(i  ,j-1) + Ro(i  ,j  ) ) * 0.5
-
 
Ro_n = ( Ro(i  ,j+1) + Ro(i  ,j  ) ) * 0.5
Ro_n = ( Ro(i  ,j+1) + Ro(i  ,j  ) ) * 0.5
Area_w = Y_xi(i-1,j-1)
Area_w = Y_xi(i-1,j-1)
-
 
Area_e = Y_xi(i  ,j-1)
Area_e = Y_xi(i  ,j-1)
-
 
Area_s = X_et(i-1,j-1)
Area_s = X_et(i-1,j-1)
-
 
Area_n = X_et(i-1,j  )
Area_n = X_et(i-1,j  )
Del_e  = Del_X_et(i  ,j  )  
Del_e  = Del_X_et(i  ,j  )  
-
 
Del_w  = Del_X_et(i-1,j  )
Del_w  = Del_X_et(i-1,j  )
-
 
Del_n  = Del_Y_xi(i  ,j  )
Del_n  = Del_Y_xi(i  ,j  )
-
 
Del_s  = Del_Y_xi(i  ,j-1)
Del_s  = Del_Y_xi(i  ,j-1)
! upwind differencing (all other will be included into the source term)  
! upwind differencing (all other will be included into the source term)  
-
Con_e = Area_e *  ( F(i,j,1) + F(i+1,j  ,1) ) * 0.5
+
Con_e = Area_e *  ( F(i,j,1) + F(i+1,j  ,1) ) * 0.5
Con_w = Area_w *  ( F(i,j,1) + F(i-1,j  ,1) ) * 0.5
Con_w = Area_w *  ( F(i,j,1) + F(i-1,j  ,1) ) * 0.5
-
 
Con_s = Area_s *  ( F(i,j,2) + F(i  ,j-1,2) ) * 0.5
Con_s = Area_s *  ( F(i,j,2) + F(i  ,j-1,2) ) * 0.5
-
 
Con_n = Area_n *  ( F(i,j,2) + F(i  ,j+1,2) ) * 0.5
Con_n = Area_n *  ( F(i,j,2) + F(i  ,j+1,2) ) * 0.5
Diff_e = Area_e * Gam_e / Del_e
Diff_e = Area_e * Gam_e / Del_e
-
 
  Diff_w = Area_w * Gam_w / Del_w  
  Diff_w = Area_w * Gam_w / Del_w  
-
 
Diff_s = Area_s * Gam_s / Del_s
Diff_s = Area_s * Gam_s / Del_s
-
 
Diff_n = Area_n * Gam_n / Del_n  
Diff_n = Area_n * Gam_n / Del_n  
                 Flux_e = Area_e * Con_e * Ro_e
                 Flux_e = Area_e * Con_e * Ro_e
-
 
Flux_w = Area_w * Con_w * Ro_w
Flux_w = Area_w * Con_w * Ro_w
-
 
Flux_s = Area_s * Con_s * Ro_s
Flux_s = Area_s * Con_s * Ro_s
 +
Flux_n = Area_n * Con_n * Ro_n
-
Flux_n = Area_n * Con_n * Ro_n
 
Aw(i,j) = Diff_w + max(    Con_w,0.)
Aw(i,j) = Diff_w + max(    Con_w,0.)
-
 
Ae(i,j) = Diff_e + max(-1.* Con_e,0.)
Ae(i,j) = Diff_e + max(-1.* Con_e,0.)
-
 
As(i,j) = Diff_s + max(    Con_s,0.)
As(i,j) = Diff_s + max(    Con_s,0.)
-
 
An(i,j) = Diff_n + max(-1.* Con_n,0.)
An(i,j) = Diff_n + max(-1.* Con_n,0.)
 +
         CheckFlux(i,j) =  Flux_e - Flux_w + Flux_s - Flux_n
         CheckFlux(i,j) =  Flux_e - Flux_w + Flux_s - Flux_n
-
Ap(i,j)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j) + CheckFlux(i,j)
+
 
 +
Ap(i,j)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j) + CheckFlux(i,j)
Sp(i,j)= 0.
Sp(i,j)= 0.
-
!--------------------------------QUICK SCHEME-------------------------
+
 +
!-------------------------------- QUICK SCHEME----------------------------
 +
    go to 800 !(now quick is "off")
-
  go to 700
+
! Subroutine QUICK(Con_f,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)
-
  if( (i.GT.2).AND.(i.LT.NXmax).and.(j.GT.4).AND.(j.LT.NYmax) ) then  
+
  if( (i.GT.2).AND.(i.LT.NXmax-0).and.(j.GT.2).AND.(j.LT.NYmax-0) ) then  
-
  if(Con_e.GT.0.) Sp(i,j) = Sp(i,j) + Con_e * (-1.) * (-0.125 * F(i-1,j ,nf) - 0.25 * F(i  ,j ,nf) + 0.375 * F(i-1,j ,nf) )
+
  !------------------ w face -------------------
 +
Fi_U = F(i-2,j,nf)
 +
Fi_C  = F(i-1,j,nf)
 +
Fi_D  = F(i  ,j,nf)
 +
Fi_DD = F(i+1,j,nf)
-
  if(Con_w.GT.0.) Sp(i,j) = Sp(i,j) + Con_w * (-1.) * (-0.125 * F(i-2,,nf) - 0.25 * F(i-1,,nf) + 0.375 * F(i-2,j  ,nf) )
+
call QUICK(Con_w,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)
-
  if(Con_n.GT.0.)  Sp(i,j) = Sp(i,j) + Con_n * (-1.) * (-0.125 * F(i  ,j-1,nf) - 0.25 * F(i  ,j  ,nf) + 0.375 * F(i  ,j-1,nf) )
+
Sp(i,j) = Sp(i,j) + Con_w * Delta_f
-
   if(Con_s.GT.0.)  Sp(i,j) = Sp(i,j) + Con_s * (-1.) * (-0.125 * F(i  ,j-2,nf) - 0.25 * F(i  ,j-1,nf) + 0.375 * F(i  ,j-2,nf) )
+
   !------------------ e face--------------------
-
  if(Con_e.LT.0.) Sp(i,j) = Sp(i,j) + Con_e * (-1.) * ( 0.375 * F(i  ,j ,nf) - 0.25 * F(i+1,j ,nf) - 0.125 * F(i+2,j ,nf) )
+
Fi_U = F(i-1,j,nf)
 +
Fi_C  = F(i  ,j,nf)
 +
Fi_D  = F(i+1,j,nf)
 +
Fi_DD = F(i+2,j,nf)
-
  if(Con_w.LT.0.) Sp(i,j) = Sp(i,j) + Con_w * (-1.) * ( 0.375 * F(i-1,,nf) - 0.25 * F(i  ,j  ,nf) - 0.125 * F(i+1,j  ,nf) )
+
call QUICK(Con_e,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)
-
  if(Con_n.LT.0.)  Sp(i,j) = Sp(i,j) + Con_n * (-1.) * ( 0.375 * F(i  ,j  ,nf) - 0.25 * F(i  ,j+1,nf) - 0.125 * F(i  ,j+2,nf) )
+
    Sp(i,j) = Sp(i,j) - Con_e * Delta_f               
-
! if(Con_s.LT.0.) Sp(i,j) = Sp(i,j) + Con_s * (-1.) * ( 0.375 * F(i  ,j-1,nf) - 0.25 * F(i  ,j+0,nf) - 0.125 * F(i  ,j+1,nf) )
+
  !------------------ s face--------------------
-
+
Fi_U = F(i ,j-2,nf)
-
+
Fi_C  = F(i  ,j-1,nf)
-
end if
+
Fi_D  = F(i  ,j ,nf)
 +
Fi_DD = F(i  ,j+1,nf)
-
  700 continue
+
call  QUICK(Con_s,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)
-
!--------------------------------QUICK SCHEME-------------------------
+
  Sp(i,j) = Sp(i,j) + Con_s * Delta_f
-
!------------------------------------------------------------------------
+
  !------------------ n face--------------------
-
!------------------------------------------------------------------------
+
Fi_U  = F(i  ,j-1,nf)
 +
Fi_C  = F(i  ,j  ,nf)
 +
Fi_D  = F(i  ,j+1,nf)
 +
Fi_DD = F(i  ,j+2,nf)
-
!------------------------------------------------------------------------
+
call  QUICK(Con_n,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)
 +
 
 +
  Sp(i,j) = Sp(i,j) - Con_n * Delta_f           
 +
 
 +
 
 +
end if
 +
 
 +
  800 continue
 +
 
 +
!-------------------------------- QUICK SCHEME----------------------------
!-------------------------------- HLPA SCHEME----------------------------
!-------------------------------- HLPA SCHEME----------------------------
-
!    go to 600
+
    go to 600 (now HLPA is "off")
! Subroutine HLPA(Uw,Fww,Fw,Fp,Fe,Delta_f)
! Subroutine HLPA(Uw,Fww,Fw,Fp,Fe,Delta_f)
Line 136: Line 149:
   !------------------ w face -------------------
   !------------------ w face -------------------
Fww = F(i-2,j,nf)
Fww = F(i-2,j,nf)
-
 
Fw  = F(i-1,j,nf)
Fw  = F(i-1,j,nf)
-
 
Fp  = F(i  ,j,nf)
Fp  = F(i  ,j,nf)
-
 
Fe  = F(i+1,j,nf)
Fe  = F(i+1,j,nf)
Line 150: Line 160:
Fww = F(i-1,j,nf)
Fww = F(i-1,j,nf)
-
 
Fw  = F(i  ,j,nf)
Fw  = F(i  ,j,nf)
-
 
Fp  = F(i+1,j,nf)
Fp  = F(i+1,j,nf)
-
 
Fe  = F(i+2,j,nf)
Fe  = F(i+2,j,nf)
Line 163: Line 170:
   !------------------ s face--------------------
   !------------------ s face--------------------
Fww = F(i  ,j-2,nf)
Fww = F(i  ,j-2,nf)
-
 
Fw  = F(i  ,j-1,nf)
Fw  = F(i  ,j-1,nf)
-
 
Fp  = F(i  ,j  ,nf)
Fp  = F(i  ,j  ,nf)
-
 
Fe  = F(i  ,j+1,nf)
Fe  = F(i  ,j+1,nf)
call  HLPA(Con_s,Fww,Fw,Fp,Fe,Delta_f)
call  HLPA(Con_s,Fww,Fw,Fp,Fe,Delta_f)
-
 
   Sp(i,j) = Sp(i,j) + Con_s * Delta_f
   Sp(i,j) = Sp(i,j) + Con_s * Delta_f
Line 178: Line 181:
Fww = F(i  ,j-1,nf)
Fww = F(i  ,j-1,nf)
-
 
Fw  = F(i  ,j  ,nf)
Fw  = F(i  ,j  ,nf)
-
 
Fp  = F(i  ,j+1,nf)
Fp  = F(i  ,j+1,nf)
-
 
Fe  = F(i  ,j+2,nf)
Fe  = F(i  ,j+2,nf)
call  HLPA(Con_n,Fww,Fw,Fp,Fe,Delta_f)
call  HLPA(Con_n,Fww,Fw,Fp,Fe,Delta_f)
-
 
   Sp(i,j) = Sp(i,j) + Con_n * Delta_f *(-1.)
   Sp(i,j) = Sp(i,j) + Con_n * Delta_f *(-1.)
Line 202: Line 201:
   
   
!----------------------------------------------------------------
!----------------------------------------------------------------
-
 
NImax = NXmaxp
NImax = NXmaxp
-
 
NJmax = NYmaxp
NJmax = NYmaxp
Line 217: Line 214:
Return
Return
 +
End
-
End
+
</pre>

Latest revision as of 17:58, 20 April 2012


!Sample program for solving Smith-Hutton Test using different schemes 
!of covective terms approximation - Coefficients computing modul
!Copyright (C) 2005  Michail Kiričkov

!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 Coef_1(nf)
include 'icomm_1.f90'  

Dimension F_out(nx,ny),CheckFlux(nx,ny)
Character  Filename*10

!  calculation of fluxes
!  all geometry has rectangular 2D notation

	Do 100 I= 2,NXmax
	Do 100 J= 2,NYmax

		Gam_e = ( Gam(i+1,j  ) + Gam(i  ,j  ) ) * 0.5
		Gam_w = ( Gam(i-1,j  ) + Gam(i  ,j  ) ) * 0.5
		Gam_s = ( Gam(i  ,j-1) + Gam(i  ,j  ) ) * 0.5
		Gam_n = ( Gam(i  ,j+1) + Gam(i  ,j  ) ) * 0.5
		
		Ro_e = ( Ro(i+1,j  ) + Ro(i  ,j  ) ) * 0.5
		Ro_w = ( Ro(i-1,j  ) + Ro(i  ,j  ) ) * 0.5
		Ro_s = ( Ro(i  ,j-1) + Ro(i  ,j  ) ) * 0.5
		Ro_n = ( Ro(i  ,j+1) + Ro(i  ,j  ) ) * 0.5

		Area_w = Y_xi(i-1,j-1)
		Area_e = Y_xi(i  ,j-1)
		Area_s = X_et(i-1,j-1)
		Area_n = X_et(i-1,j  )

		Del_e  = Del_X_et(i  ,j  ) 
		Del_w  = Del_X_et(i-1,j  )
		Del_n  = Del_Y_xi(i  ,j  )
		Del_s  = Del_Y_xi(i  ,j-1)

! upwind differencing (all other will be included into the source term) 

		Con_e = Area_e *  ( F(i,j,1) + F(i+1,j  ,1) ) * 0.5	
		Con_w = Area_w *  ( F(i,j,1) + F(i-1,j  ,1) ) * 0.5
		Con_s = Area_s *  ( F(i,j,2) + F(i  ,j-1,2) ) * 0.5
		Con_n = Area_n *  ( F(i,j,2) + F(i  ,j+1,2) ) * 0.5

		Diff_e = Area_e * Gam_e / Del_e
 		Diff_w = Area_w * Gam_w / Del_w 
		Diff_s = Area_s * Gam_s / Del_s
		Diff_n = Area_n * Gam_n / Del_n 

                Flux_e = Area_e * Con_e * Ro_e
		Flux_w = Area_w * Con_w * Ro_w
		Flux_s = Area_s * Con_s * Ro_s
		Flux_n = Area_n * Con_n * Ro_n

		
		Aw(i,j) = Diff_w + max(     Con_w,0.)
		Ae(i,j) = Diff_e + max(-1.* Con_e,0.)
		As(i,j) = Diff_s + max(     Con_s,0.)
		An(i,j) = Diff_n + max(-1.* Con_n,0.)


        CheckFlux(i,j) =  Flux_e - Flux_w + Flux_s - Flux_n


	Ap(i,j)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j) + CheckFlux(i,j)
	
		Sp(i,j)= 0.

 
!-------------------------------- QUICK SCHEME----------------------------
     go to 800 !(now quick is "off")

! Subroutine QUICK(Con_f,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)

 if( (i.GT.2).AND.(i.LT.NXmax-0).and.(j.GT.2).AND.(j.LT.NYmax-0) ) then 

   !------------------ w face -------------------
	Fi_U  = F(i-2,j,nf)
	Fi_C  = F(i-1,j,nf)
	Fi_D  = F(i  ,j,nf)
	Fi_DD = F(i+1,j,nf)

	call  QUICK(Con_w,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)

	Sp(i,j) = Sp(i,j) + Con_w * Delta_f

  !------------------ e face--------------------

	Fi_U  = F(i-1,j,nf)
	Fi_C  = F(i  ,j,nf)
	Fi_D  = F(i+1,j,nf)
	Fi_DD = F(i+2,j,nf)

	call  QUICK(Con_e,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)

    Sp(i,j) = Sp(i,j) - Con_e * Delta_f                

  !------------------ s face--------------------
	Fi_U  = F(i  ,j-2,nf)
	Fi_C  = F(i  ,j-1,nf)
	Fi_D  = F(i  ,j  ,nf)
	Fi_DD = F(i  ,j+1,nf)

	call  QUICK(Con_s,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)

   Sp(i,j) = Sp(i,j) + Con_s * Delta_f

  !------------------ n face--------------------

	Fi_U  = F(i  ,j-1,nf)
	Fi_C  = F(i  ,j  ,nf)
	Fi_D  = F(i  ,j+1,nf)
	Fi_DD = F(i  ,j+2,nf)

	call  QUICK(Con_n,Fi_U,Fi_C,Fi_D,Fi_DD,Delta_f)

   Sp(i,j) = Sp(i,j) - Con_n * Delta_f             


 end if	

   800 continue	

!-------------------------------- QUICK SCHEME----------------------------

!-------------------------------- HLPA SCHEME----------------------------
    go to 600 (now HLPA is "off")

! Subroutine HLPA(Uw,Fww,Fw,Fp,Fe,Delta_f)

 if( (i.GT.2).AND.(i.LT.NXmax-0).and.(j.GT.2).AND.(j.LT.NYmax-0) ) then 

   !------------------ w face -------------------
	Fww = F(i-2,j,nf)
	Fw  = F(i-1,j,nf)
	Fp  = F(i  ,j,nf)
	Fe  = F(i+1,j,nf)

	call  HLPA(Con_w,Fww,Fw,Fp,Fe,Delta_f)

	Sp(i,j) = Sp(i,j) + Con_w * Delta_f

  !------------------ e face--------------------

	Fww = F(i-1,j,nf)
	Fw  = F(i  ,j,nf)
	Fp  = F(i+1,j,nf)
	Fe  = F(i+2,j,nf)

	call  HLPA(Con_e,Fww,Fw,Fp,Fe,Delta_f)

    Sp(i,j) = Sp(i,j) + Con_e * Delta_f * (-1.)

  !------------------ s face--------------------
	Fww = F(i  ,j-2,nf)
	Fw  = F(i  ,j-1,nf)
	Fp  = F(i  ,j  ,nf)
	Fe  = F(i  ,j+1,nf)

	call  HLPA(Con_s,Fww,Fw,Fp,Fe,Delta_f)

   Sp(i,j) = Sp(i,j) + Con_s * Delta_f

  !------------------ n face--------------------

	Fww = F(i  ,j-1,nf)
	Fw  = F(i  ,j  ,nf)
	Fp  = F(i  ,j+1,nf)
	Fe  = F(i  ,j+2,nf)

	call  HLPA(Con_n,Fww,Fw,Fp,Fe,Delta_f)

   Sp(i,j) = Sp(i,j) + Con_n * Delta_f *(-1.)


 end if	

   600 continue	

!-------------------------------- HLPA SCHEME----------------------------

!------------------------------------------------------------------------
		
100 continue
 
!----------------------------------------------------------------
	NImax = NXmaxp
	NJmax = NYmaxp

	F_out     = CheckFlux

 	Filename  ='0_Flx.txt' 

    Call  Out_array(F_out,NImax,NJmax,Filename)

!-------------------------------------------------------------------


Return
End

My wiki