Installing FORTRAN 90 Compilier
Fortran Coding Tutorials:
1-A code example that reads data in from a text file and prints out 3 values to check that data was read in.
1-A code example that reads data in from a text file and prints out 3 values to check that data was read in.
Link on how to instal the fortran90 compiler on a windows machine.
Youtube link on how to instal FORTRAN90
Programming Using FORTRAN 90
The following coe reads in data from a text file.
PROGRAM DF
INTEGER::N
REAL:: B,C,num,I
REAL,DIMENSION(100)::A,BB,CC
OPEN(unit=111,file='my_data.txt',status='unknown')
DO I=1,100
READ(111,*)
CC(I)
END DO
close(111)
print*, CC(25), CC(50), CC(99)
END PROGRAM DF
INTEGER::N
REAL:: B,C,num,I
REAL,DIMENSION(100)::A,BB,CC
OPEN(unit=111,file='my_data.txt',status='unknown')
DO I=1,100
READ(111,*)
CC(I)
END DO
close(111)
print*, CC(25), CC(50), CC(99)
END PROGRAM DF
The following tutorial guides you to write data file:
The following program converts from degress to radians, in addition it uses a function:
PROGRAM GH
IMPLICIT NONE
REAL::CC
REAL::ANGLEDEG
ANGLEDEG=30.0
CC = CONVERSION(ANGLEDEG)
PRINT*,CC
CONTAINS REAL FUNCTION CONVERSION(ANGLEDEG)
REAL::C,PI,ANGLEDEG
INTEGER::I
PI=4.0*ATAN(1.0)
C=PI/180
CONVERSION = C*ANGLEDEG
END
FUNCTION CONVERSION
END PROGRAM GH
IMPLICIT NONE
REAL::CC
REAL::ANGLEDEG
ANGLEDEG=30.0
CC = CONVERSION(ANGLEDEG)
PRINT*,CC
CONTAINS REAL FUNCTION CONVERSION(ANGLEDEG)
REAL::C,PI,ANGLEDEG
INTEGER::I
PI=4.0*ATAN(1.0)
C=PI/180
CONVERSION = C*ANGLEDEG
END
FUNCTION CONVERSION
END PROGRAM GH
It would be advantagous if there is some material on Solving ODEs, Solving PDEs, Making initilization profiles, Domain creation, working on cross section planes in domains, takeing samples pf domain cells, finite diffrence schemes, finite volume schemes,using trignometric functions, cubic spline interpolation, Lagrange polynomial interpolation, creating 2D shapes, algebric operations on matrcies, using IF statements, using subroutines, having control over your output, the application of counters.
Online Material on CFD Fortran Coding
To write your CFD code it is recommended by the scientfic community to use FORTRAN 90 as your programming languadge:
1-Boston University Information Servicies & Technology
2-University of Oklahoma Fortran Coding by Dr.Ming Xue:
3-Uppsala Universitet Numerical Hydrodynamics and Radiative Transfer.
4-Universite De Bordeaux, Index
5-Code Paralization Material.
6-Code Paralization Material from the University of Texas Austin.
7-Fortran Source Codes.
8-Cornell Theory Center Parallel Processing Concepts Saleh Elmohamd.
9-UPMC SORBONNE UNIVERSITIES Stephane Zaleski
10-Finite Volume Method Tutorials by G Recktenwald
11-Computational Hydraulics by Dr David Apsley
1-Boston University Information Servicies & Technology
2-University of Oklahoma Fortran Coding by Dr.Ming Xue:
3-Uppsala Universitet Numerical Hydrodynamics and Radiative Transfer.
4-Universite De Bordeaux, Index
5-Code Paralization Material.
6-Code Paralization Material from the University of Texas Austin.
7-Fortran Source Codes.
8-Cornell Theory Center Parallel Processing Concepts Saleh Elmohamd.
9-UPMC SORBONNE UNIVERSITIES Stephane Zaleski
10-Finite Volume Method Tutorials by G Recktenwald
11-Computational Hydraulics by Dr David Apsley
Excellent Website for Working in a Supercomputing Center
The link provides lots of tips using supercomputting facilities : http://www.mcs.anl.gov/hs/software/index.php
With Notes on the use of Fortran 90 : http://www.mcs.anl.gov/~itf/dbpp/text/node82.html
With Notes on the use of Fortran 90 : http://www.mcs.anl.gov/~itf/dbpp/text/node82.html
Michal Kopera CFD Lab Notes
The lessons have been put up with permission from Michal Kopera.
Exmple One-Dimensional Case
PROGRAM TDMA
IMPLICIT NONE
REAL,DIMENSION(:),ALLOCATABLE::A,B,C,D,U0,f,U
REAL::AA,BB,CC,DELTA_Y,DELTA_T,L,H,ID,Y,T,DP
INTEGER::I,N,P,TIME
N=32
ALLOCATE(U0(N+2),U(N+2),D(N),C(N),A(N),B(N))
H=1
L=2*H
DELTA_T=0.1
DELTA_Y=L/N
AA=-1/(2*(DELTA_Y)**2)
BB=(1/DELTA_T-(2*AA))
CC=AA
!INITILAIZATION PROFILE
P=1
U0(1)=0
DO I=-15,15
P=P+1
DP=(-(I*DELTA_Y)**4)
U0(P)=-(2.13)**(-1)+(2.13)**DP
END DO
U0(N)=0
DO I=1,N
C(I)=AA
A(I)=AA
B(I)=BB
END DO
D(1)=1+BB*U0(1)+AA*U0(2)
DO I=2,N-1
D(I)=1+AA*U0(I-1)+BB*U0(I)+AA*U0(I+1)
END DO
D(N)=1+AA*U0(N-1)+BB*U0(N)
CALL TD(A,B,C,U,D,N)
open (unit=111, file='VISH.dat' ,status='unknown')
DO I=1,N+1
WRITE(111,*) U(I)
END DO
CLOSE(111)
!-----------------------------------------------
!TIME STEPPING
DO TIME=1,100
DO I=1,N
C(I)=AA
A(I)=AA
B(I)=BB
END DO
D(1)=1+BB*U(1)+AA*U(2)
DO I=2,N-1
D(I)=1+AA*U(I-1)+BB*U(I)+AA*U(I+1)
END DO
D(N)=1+AA*U(N-1)+BB*U(N)
CALL TD(A,B,C,U,D,N)
END DO
open (unit=222, file='VISH.dat' ,status='unknown')
DO I=1,N
WRITE(222,*) U(I)
END DO
CLOSE(222)
DEALLOCATE(A,B,C,U,D,U0)
!END TIME STEPPING
!-----------------------------------------------
CONTAINS
SUBROUTINE TD(A,B,C,U,D,N)
INTEGER,INTENT(IN)::N
REAL,DIMENSION(N),INTENT(INOUT)::A,B,C,D
REAL,DIMENSION(N),INTENT(OUT)::U
INTEGER::I
!forward elimination
DO I=2,N
B(I) = B(I) - A(I)/B(I-1)*C(I-1)
D(I) = D(I) - A(I)/B(I-1)*D(I-1)
END DO
!backward substitution
!mind the inverse direction of the do loop
U(N) = D(N)/B(N)
DO I=N-1,1,-1
U(I) = (D(I) - c(I)*U(I+1))/B(I)
END DO
END SUBROUTINE TD
END PROGRAM TDMA
IMPLICIT NONE
REAL,DIMENSION(:),ALLOCATABLE::A,B,C,D,U0,f,U
REAL::AA,BB,CC,DELTA_Y,DELTA_T,L,H,ID,Y,T,DP
INTEGER::I,N,P,TIME
N=32
ALLOCATE(U0(N+2),U(N+2),D(N),C(N),A(N),B(N))
H=1
L=2*H
DELTA_T=0.1
DELTA_Y=L/N
AA=-1/(2*(DELTA_Y)**2)
BB=(1/DELTA_T-(2*AA))
CC=AA
!INITILAIZATION PROFILE
P=1
U0(1)=0
DO I=-15,15
P=P+1
DP=(-(I*DELTA_Y)**4)
U0(P)=-(2.13)**(-1)+(2.13)**DP
END DO
U0(N)=0
DO I=1,N
C(I)=AA
A(I)=AA
B(I)=BB
END DO
D(1)=1+BB*U0(1)+AA*U0(2)
DO I=2,N-1
D(I)=1+AA*U0(I-1)+BB*U0(I)+AA*U0(I+1)
END DO
D(N)=1+AA*U0(N-1)+BB*U0(N)
CALL TD(A,B,C,U,D,N)
open (unit=111, file='VISH.dat' ,status='unknown')
DO I=1,N+1
WRITE(111,*) U(I)
END DO
CLOSE(111)
!-----------------------------------------------
!TIME STEPPING
DO TIME=1,100
DO I=1,N
C(I)=AA
A(I)=AA
B(I)=BB
END DO
D(1)=1+BB*U(1)+AA*U(2)
DO I=2,N-1
D(I)=1+AA*U(I-1)+BB*U(I)+AA*U(I+1)
END DO
D(N)=1+AA*U(N-1)+BB*U(N)
CALL TD(A,B,C,U,D,N)
END DO
open (unit=222, file='VISH.dat' ,status='unknown')
DO I=1,N
WRITE(222,*) U(I)
END DO
CLOSE(222)
DEALLOCATE(A,B,C,U,D,U0)
!END TIME STEPPING
!-----------------------------------------------
CONTAINS
SUBROUTINE TD(A,B,C,U,D,N)
INTEGER,INTENT(IN)::N
REAL,DIMENSION(N),INTENT(INOUT)::A,B,C,D
REAL,DIMENSION(N),INTENT(OUT)::U
INTEGER::I
!forward elimination
DO I=2,N
B(I) = B(I) - A(I)/B(I-1)*C(I-1)
D(I) = D(I) - A(I)/B(I-1)*D(I-1)
END DO
!backward substitution
!mind the inverse direction of the do loop
U(N) = D(N)/B(N)
DO I=N-1,1,-1
U(I) = (D(I) - c(I)*U(I+1))/B(I)
END DO
END SUBROUTINE TD
END PROGRAM TDMA
Example One Dimensional Advection (need to check it)
PROGRAM ADVECTION
IMPLICIT NONE
INTEGER :: iNx,iNt,i,j
DOUBLE PRECISION :: cfl=1,u,rDx,rDt,F(500),f1(500),f2(500),w1(500),w2(500),l1(500),l2(500),X(100)
! Input the variables required
iNx=60
iNt=100
u=1
rDx=80./(iNx-1)
rDt=cfl*rDx/u
! initialise
x(1)=-0.5
x(iNx)=0.5
DO i=2,iNx-1
F(i)=0.5*EXP(-(X(i)/0.1)**4)
END DO
OPEN(unit=111,file="cfl1_exp_s olution_init.dat",status="unknown")
!WRITE(1,*) 'TITLE="Numeric Solution for the Advection Equation"'
!WRITE(1,*) 'VARIABLES="x" "Amplitude"'
!WRITE(1,*) "ZONE I=",iNx," F=POINT"
DO i=2,iNx-1
WRITE(111,*)F(i)
END DO
close(111)
! solution with lax-Wendroff formulation
l1=F
l2=l1
! start time loop
DO j=1,iNt
DO i=1,iNx
l2(i)=l1(i)-cfl/2.*(l1(i+1)-l1(i-1))+cfl**2./2.*(l1(i+1)-2.*l1(i)+l1(i-1))
END DO
l1=l2
END DO
OPEN (unit=222,file="cfl1_exp_lax_wendroff.dat",status="unknown")
!WRITE(1,*) 'TITLE="Numeric Solution for the Advection Equation"'
!WRITE(1,*) 'VARIABLES="x" "Amplitude"'
!WRITE(1,*) "ZONE I=",iNx," F=POINT"
DO i=1,iNx
WRITE(222,*)l1(i)!-40.+(i-1)*rDx, l1(i)
END DO
CLOSE(222)
END PROGRAM ADVECTION
IMPLICIT NONE
INTEGER :: iNx,iNt,i,j
DOUBLE PRECISION :: cfl=1,u,rDx,rDt,F(500),f1(500),f2(500),w1(500),w2(500),l1(500),l2(500),X(100)
! Input the variables required
iNx=60
iNt=100
u=1
rDx=80./(iNx-1)
rDt=cfl*rDx/u
! initialise
x(1)=-0.5
x(iNx)=0.5
DO i=2,iNx-1
F(i)=0.5*EXP(-(X(i)/0.1)**4)
END DO
OPEN(unit=111,file="cfl1_exp_s olution_init.dat",status="unknown")
!WRITE(1,*) 'TITLE="Numeric Solution for the Advection Equation"'
!WRITE(1,*) 'VARIABLES="x" "Amplitude"'
!WRITE(1,*) "ZONE I=",iNx," F=POINT"
DO i=2,iNx-1
WRITE(111,*)F(i)
END DO
close(111)
! solution with lax-Wendroff formulation
l1=F
l2=l1
! start time loop
DO j=1,iNt
DO i=1,iNx
l2(i)=l1(i)-cfl/2.*(l1(i+1)-l1(i-1))+cfl**2./2.*(l1(i+1)-2.*l1(i)+l1(i-1))
END DO
l1=l2
END DO
OPEN (unit=222,file="cfl1_exp_lax_wendroff.dat",status="unknown")
!WRITE(1,*) 'TITLE="Numeric Solution for the Advection Equation"'
!WRITE(1,*) 'VARIABLES="x" "Amplitude"'
!WRITE(1,*) "ZONE I=",iNx," F=POINT"
DO i=1,iNx
WRITE(222,*)l1(i)!-40.+(i-1)*rDx, l1(i)
END DO
CLOSE(222)
END PROGRAM ADVECTION
Numerical Integration using the Simpson Rule
PROGRAM DF
REAL::L,N,DELTA_X,AREA,A,I,ODD_AREA,EVEN_AREA
REAL,DIMENSION(100)::F,B,C
L=1.0
N=100
A=1.0/3.0
!PRINT*, A
DELTA_X=L/N
F(1)=(DELTA_X)**2+2*DELTA_X+2
!PRINT*, F(1)
F(N)=(N*DELTA_X)**2+2*N*DELTA_X+2
!PRINT*, F(N)
DO I=2,N-1
IF (MOD(I,2.0)==0) THEN B(I)=(I*DELTA_X)**2+2*I*DELTA_X+2
!PRINT*, B(I)
ELSE IF (MOD(I,2.0)>0) THEN C(I)=(I*DELTA_X)**2+2*I*DELTA_X+2
!PRINT*, C(I)
END IF
END DO
ODD_AREA=SUM(B)
!PRINT*, ODD_AREA
EVEN_AREA=SUM(C)
!PRINT*, EVEN_AREA
AREA=A*DELTA_X*(F(1)+F(N)+ODD_AREA+EVEN_AREA)
PRINT*,'TOTAL AREA', AREA
END PROGRAM DF
REAL::L,N,DELTA_X,AREA,A,I,ODD_AREA,EVEN_AREA
REAL,DIMENSION(100)::F,B,C
L=1.0
N=100
A=1.0/3.0
!PRINT*, A
DELTA_X=L/N
F(1)=(DELTA_X)**2+2*DELTA_X+2
!PRINT*, F(1)
F(N)=(N*DELTA_X)**2+2*N*DELTA_X+2
!PRINT*, F(N)
DO I=2,N-1
IF (MOD(I,2.0)==0) THEN B(I)=(I*DELTA_X)**2+2*I*DELTA_X+2
!PRINT*, B(I)
ELSE IF (MOD(I,2.0)>0) THEN C(I)=(I*DELTA_X)**2+2*I*DELTA_X+2
!PRINT*, C(I)
END IF
END DO
ODD_AREA=SUM(B)
!PRINT*, ODD_AREA
EVEN_AREA=SUM(C)
!PRINT*, EVEN_AREA
AREA=A*DELTA_X*(F(1)+F(N)+ODD_AREA+EVEN_AREA)
PRINT*,'TOTAL AREA', AREA
END PROGRAM DF
Numerical Integration using the Trapzodial Rule
PROGRAM DF
REAL::L,N,DELTA_X,AREA
REAL,DIMENSION(100)::F,B
L=1
N=10
DELTA_X=L/N
F(I)=(I*DELTA_X)**2+2*I*DELTA_X+2
DO I=2,N-1
B(I)=(I*DELTA_X)**2+2*I*DELTA_X+2
END DO
AREA=DELTA_X*(0.5*F(1)+0.5*F(N)+SUM(B))
PRINT*,AREA
END PROGRAM DF
REAL::L,N,DELTA_X,AREA
REAL,DIMENSION(100)::F,B
L=1
N=10
DELTA_X=L/N
F(I)=(I*DELTA_X)**2+2*I*DELTA_X+2
DO I=2,N-1
B(I)=(I*DELTA_X)**2+2*I*DELTA_X+2
END DO
AREA=DELTA_X*(0.5*F(1)+0.5*F(N)+SUM(B))
PRINT*,AREA
END PROGRAM DF
Creating Tables
PROGRAM GH
IMPLICIT NONE
REAL::ANGLEDEG
INTEGER::I,N,T
REAL,DIMENSION(1000)::CC,VV
T=0
N=100
ANGLEDEG=30.0
PRINT* ,'----------NUMBER--ANGLE_DEG-----ANGLE_RAD-----'
DO I=1,N
T=T+1
ANGLEDEG=ANGLEDEG+0.01
VV(I)=ANGLEDEG
CC(I) = CONVERSION(ANGLEDEG)
PRINT*,T,'---',VV(I),CC(I)
END DO
OPEN(UNIT=111,FILE='DATA.DAT',STATUS='UNKNOWN')
WRITE(111,*) '----------NUMBER--ANGLE_DEG-----ANGLE_RAD-----'
T=0
DO I=1,N
T=T+1
WRITE(111,*) T,'---',VV(I),CC(I) END DO
CLOSE(111)
CONTAINS REAL
FUNCTION CONVERSION(ANGLEDEG)
REAL::C,PI,ANGLEDEG
PI=4.0*ATAN(1.0)
C=PI/180
CONVERSION = C*ANGLEDEG END FUNCTION CONVERSION
END PROGRAM GH
IMPLICIT NONE
REAL::ANGLEDEG
INTEGER::I,N,T
REAL,DIMENSION(1000)::CC,VV
T=0
N=100
ANGLEDEG=30.0
PRINT* ,'----------NUMBER--ANGLE_DEG-----ANGLE_RAD-----'
DO I=1,N
T=T+1
ANGLEDEG=ANGLEDEG+0.01
VV(I)=ANGLEDEG
CC(I) = CONVERSION(ANGLEDEG)
PRINT*,T,'---',VV(I),CC(I)
END DO
OPEN(UNIT=111,FILE='DATA.DAT',STATUS='UNKNOWN')
WRITE(111,*) '----------NUMBER--ANGLE_DEG-----ANGLE_RAD-----'
T=0
DO I=1,N
T=T+1
WRITE(111,*) T,'---',VV(I),CC(I) END DO
CLOSE(111)
CONTAINS REAL
FUNCTION CONVERSION(ANGLEDEG)
REAL::C,PI,ANGLEDEG
PI=4.0*ATAN(1.0)
C=PI/180
CONVERSION = C*ANGLEDEG END FUNCTION CONVERSION
END PROGRAM GH
Using the sum Command
PROGRAME FG
IMPLICIT NONE
INTEGER::N,I
REAL,DIMENSION(100)::A
N=10
DO I=1,N
A(I)=1
END DO
print *, sum(A)
END PROGRAME FG
IMPLICIT NONE
INTEGER::N,I
REAL,DIMENSION(100)::A
N=10
DO I=1,N
A(I)=1
END DO
print *, sum(A)
END PROGRAME FG
Using the mod Command
PROGRAM df
IMPLICIT NONE
REAL:: I
DO I=1,10
if (mod(I,2.0)>0) then
print*, 'yes'
else if (mod(I,2.0)==0) then
print*, 'no'
end if
END DO
END PROGRAM df
IMPLICIT NONE
REAL:: I
DO I=1,10
if (mod(I,2.0)>0) then
print*, 'yes'
else if (mod(I,2.0)==0) then
print*, 'no'
end if
END DO
END PROGRAM df
Using Functions
PROGRAM GH
IMPLICIT NONE
REAL::CC
REAL::ANGLEDEG
ANGLEDEG=30.0
CC = CONVERSION(ANGLEDEG)
PRINT*,CC
CONTAINS
REAL FUNCTION CONVERSION(ANGLEDEG)
REAL::C,PI,ANGLEDEG
INTEGER::I
PI=4.0*ATAN(1.0)
C=PI/180
CONVERSION = C*ANGLEDEG
END FUNCTION CONVERSION
END PROGRAM GH
IMPLICIT NONE
REAL::CC
REAL::ANGLEDEG
ANGLEDEG=30.0
CC = CONVERSION(ANGLEDEG)
PRINT*,CC
CONTAINS
REAL FUNCTION CONVERSION(ANGLEDEG)
REAL::C,PI,ANGLEDEG
INTEGER::I
PI=4.0*ATAN(1.0)
C=PI/180
CONVERSION = C*ANGLEDEG
END FUNCTION CONVERSION
END PROGRAM GH
Ideal Gas Table Creation
PROGRAM GH
IMPLICIT NONE
REAL,DIMENSION(10000)::PRESSURE,TEMP,R,DENSITY,T_C,T_K,MASS
INTEGER::I,N,NN
N=100 T_C(1)=0
!PRINT*,'1234567891234567891234567891234567891234567891234567891234567'
PRINT*,'------------N-----PRESSURE---------DENSITY----------R---------'
NN=0
DO I=1,N
NN=NN+1
T_C(I+1)=T_C(I)+0.1
T_K(I)=273+T_C(I)
DENSITY(I)=1.29
R(I)=8.314
MASS(I)=1.0
PRESSURE(I)=MASS(I)*DENSITY(I)*R(I)*T_K(I)
PRINT*,'-',NN,'-',PRESSURE(I),'-',DENSITY(I),'-',R(I)
END DO
!PRINT*,'1234567891234567891234567891234567891234567891234567891234567'
PRINT*,'------------N-----T_C--------------T_K--------------MASS------'
NN=0
DO I=1,N
NN=NN+1
PRINT*,'-',NN,'-',T_C(I),'-',T_K(I),'-',MASS(I)
END DO
OPEN(UNIT=111,FILE='TEMP_VS_PRESSUR.DAT',STATUS='UNKNOWN')
DO I=1,N
WRITE(111,*) PRESSURE(I),T_K(I)
END DO
CLOSE(111)
END PROGRAM GH
IMPLICIT NONE
REAL,DIMENSION(10000)::PRESSURE,TEMP,R,DENSITY,T_C,T_K,MASS
INTEGER::I,N,NN
N=100 T_C(1)=0
!PRINT*,'1234567891234567891234567891234567891234567891234567891234567'
PRINT*,'------------N-----PRESSURE---------DENSITY----------R---------'
NN=0
DO I=1,N
NN=NN+1
T_C(I+1)=T_C(I)+0.1
T_K(I)=273+T_C(I)
DENSITY(I)=1.29
R(I)=8.314
MASS(I)=1.0
PRESSURE(I)=MASS(I)*DENSITY(I)*R(I)*T_K(I)
PRINT*,'-',NN,'-',PRESSURE(I),'-',DENSITY(I),'-',R(I)
END DO
!PRINT*,'1234567891234567891234567891234567891234567891234567891234567'
PRINT*,'------------N-----T_C--------------T_K--------------MASS------'
NN=0
DO I=1,N
NN=NN+1
PRINT*,'-',NN,'-',T_C(I),'-',T_K(I),'-',MASS(I)
END DO
OPEN(UNIT=111,FILE='TEMP_VS_PRESSUR.DAT',STATUS='UNKNOWN')
DO I=1,N
WRITE(111,*) PRESSURE(I),T_K(I)
END DO
CLOSE(111)
END PROGRAM GH
Unless otherwise noted, all content on this site is @Copyright by Ahmed Al Makky 2012-2015- http://cfd2012.com