# Coeff 11.f90 - calculate the coefficients

(Difference between revisions)
 Revision as of 01:57, 19 September 2005 (view source)Michail (Talk | contribs)← Older edit Revision as of 20:44, 20 September 2005 (view source)Michail (Talk | contribs) Newer edit → Line 1: Line 1: +
+
+                                                                                             Sample program for solving Smith-Hutton Test using different schemes of covective terms approximation - Coefficients computing modul
+                                                                                             Copyright (C) 2005  Michail Kirichkov
+
+                                                                                             This program is free software; you can redistribute it and/or
+                                                                                             modify it under the terms of the GNU General Public License
+
+                                                                                             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.)

Line 102:                                                                                                                                                               Line 78:

-                                                                                            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.
Line 132:                                                                                                                                                               Line 108:

700 continue                                                                                                                                                             700 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 145:                                                                                                                                                               Line 118:
!------------------ 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 159:                                                                                                                                                               Line 129:

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 172:                                                                                                                                                               Line 139:
!------------------ 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 187:                                                                                                                                                               Line 150:

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 211:                                                                                                                                                               Line 170:

!----------------------------------------------------------------                                                                                                        !----------------------------------------------------------------
-
NImax = NXmaxp                                                                                                                                                           NImax = NXmaxp
-
NJmax = NYmaxp                                                                                                                                                           NJmax = NYmaxp

Line 226:                                                                                                                                                               Line 183:

Return                                                                                                                                                                   Return
+                                                                                             End

-                                                                                            End                                                                        +

## Revision as of 20:44, 20 September 2005

```
Sample program for solving Smith-Hutton Test using different schemes of covective terms approximation - Coefficients computing modul

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License

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 700

if( (i.GT.2).AND.(i.LT.NXmax).and.(j.GT.4).AND.(j.LT.NYmax) ) 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) )

if(Con_w.GT.0.)  Sp(i,j) = Sp(i,j) + Con_w * (-1.) * (-0.125 * F(i-2,j  ,nf) - 0.25 * F(i-1,j  ,nf) + 0.375 * F(i-2,j  ,nf) )

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) )

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) )

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) )

if(Con_w.LT.0.)  Sp(i,j) = Sp(i,j) + Con_w * (-1.) * ( 0.375 * F(i-1,j  ,nf) - 0.25 * F(i  ,j  ,nf) - 0.125 * F(i+1,j  ,nf) )

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) )

! 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) )

end if

700 continue

!-------------------------------- 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

```