1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
|
PROGRAM NS
IMPLICIT NONE
INTEGER::i,n,t,j
!dynamically declare arrays a,b,c,u,f
REAL,DIMENSION(:),ALLOCATABLE::a,b,c,u,f
REAL::Lx,dx,dpdx,h=1,upsilon,deltaT,temp !Lx - domain size, dx - step length, dpdx pressure gradient
!Assigned values,
!number of internal points:
n=23
!domain size
Lx = 2*h
!pressure gradient
dpdx = -2.0
!viscosity, mhu/rho
upsilon=0.2
!calculate dx
dx=Lx/(n+1)
!allocate arrays a,b,c,f of size n
ALLOCATE(a(n),b(n),c(n),f(n))
!allocate u - mind that is has entries indexed from 0 to n+1
ALLOCATE(u(0:n+1))
!initial values !??? profil ???
u(0)=0
DO i=1,n-1
u(i) = -1+ (-1+i*dx)**2
END DO
u(N+1)=0
OPEN(unit=111,file='output.dat',status='unknown')
!Outer*DO*loop*for*time*index*
DO t=0,100
!Inner*DO*for*space*index*
DO i=1,n
!calculate*a(i),b(i),c(i),*and*f(i)*using Crank-Nicolson
!calculate deltaT or assigned ??
deltaT=0.1
!put proper values to a,b,c
temp=(upsilon*deltaT)/(2*(dx**2))
a(i) = -temp
b(i) = (1+2*temp)
c(i) = -temp
!calculate right-hand-size f
f(1) =temp*u(2)+(1-2*temp)*u(1)+temp*u(0)-dpdx+temp*u(0)
DO j=1,n-1
f(j) =temp*u(j+1)+(1-2*temp)*u(j)+temp*u(j-1)-dpdx
END DO
f(n) =temp*u(2)+(1-2*temp)*u(1)+temp*u(0)-dpdx+temp*u(n+1)
END DO
!call subroutine TDMA
call TDMA(a,b,c,u(1:n),f,n)
!print results to a file and on the screen
DO i=0,n+1
WRITE(111,*) (dx*i-1),u(i)
!PRINT*,dx*i,u(i)
END DO
END DO
CLOSE(111)
DEALLOCATE(a,b,c,u,f)
CONTAINS
SUBROUTINE TDMA(a,b,c,u,f,n)
!*****************************************
! subroutine solves a system of linear equations
! with tridiagonal matrix
! a - size n matrix holding lower diagonal values
! b - size n matrix holding diagonal values
! c - size n matrix holding upper diagonal values
! f - size n matrix holding right-hand-sides of equation
! u - size n matrix holding solution for internal points of domain
! (without boundary values)
!
!mind that b and f will be overwritten during the execution of subroutine
!*****************************************
INTEGER,INTENT(IN)::n
REAL,DIMENSION(n),INTENT(INOUT)::a,b,c,f
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)
f(i) = f(i) - a(i)/b(i-1)*f(i-1)
END DO
!backward substitution
!mind the inverse direction of the do loop
u(n) = f(n)/b(n)
DO i=n-1,1,-1
u(i) = (f(i) - c(i)*u(i+1))/b(i)
END DO
END SUBROUTINE TDMA
END PROGRAM NS |
Partager