Coeff 11.f90 - calculate the coefficients
From CFD-Wiki
(Difference between revisions)
(Added QUICK implementation) |
|||
Line 83: | Line 83: | ||
Sp(i,j)= 0. | Sp(i,j)= 0. | ||
- | !--------------------------------QUICK SCHEME------------------------- | + | |
+ | !-------------------------------- 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).and.(j.GT. | + | 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 | end if | ||
- | + | 800 continue | |
+ | |||
+ | !-------------------------------- QUICK SCHEME---------------------------- | ||
!-------------------------------- HLPA SCHEME---------------------------- | !-------------------------------- HLPA SCHEME---------------------------- |
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