* A. Salih, Dept. of Mechanical Engg., NIT - Trichy, India. *
*************************************************************************
* THIS PROGRAM SOLVES TRANSIENT HEAT EQUATION IN A 2-D DOMAIN *
* USING FINITE VOLUME METHOD (UNIFORM CV) *
* GENERAL IMPLICIT METHOD (BETA METHOD) *
*************************************************************************
*
program IMPLICIT_2D
implicit doubleprecision (a-h,o-z)
include 'scalars/integers.inc'
include 'scalars/parameters.inc'
open(unit=11,file='input.dat',status='unknown')
open(unit=21,file='output.dat',status='unknown')
open(unit=31,file='outputs/x.dat',status='unknown')
open(unit=32,file='outputs/y.dat',status='unknown')
open(unit=33,file='outputs/tmp.dat',status='unknown')
open(unit=34,file='outputs/time.dat',status='unknown')
open(unit=35,file='outputs/time2.dat',status='unknown')
open(unit=36,file='outputs/t_trans.dat',status='unknown')
open(unit=37,file='outputs/iterations.dat',status='unknown')
c
pi = 4*datan(1d0)
c
CALL READ_IN ! read-in data
close (unit=11)
CALL GRID ! setting up the grid points in the domain
CALL INITIAL_COND ! setting up initial condition
CALL BOUNDARY_COND ! setting up boundary values
CALL SOLVE ! setting up system of equations
CALL PRINTOUT ! printing out the computational results
CALL WRITE_OUT ! write-out the basic data used for computation
c
stop
end
*
***SUBROUTINE: READING IN THE BASIC DATA *******************************
*
subroutine READ_IN
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/solveparams.inc'
c
read(11,*) alx ! actual length of domain
read(11,*) aly ! actual length of domain
read(11,*) boundary_type ! flag for type of boundary conditions
read(11,*) alpha ! thermal diffusivity
read(11,*) Fo ! value of mean grid Fourier number for
! calculating the time step
read(11,*) solver ! flag for type of solver to be used
read(11,*) tolsor ! tolerance value for terminating the
! iteration (for SOR method)
read(11,*) wtsor ! relaxation factor for SOR
read(11,*) wtlsor ! relaxation factor for Line SOR
read(11,*) itermax ! maximum number of iterations (for solver)
read(11,*) tolstdy ! convergence criterion for steady-state
read(11,*) tfinal ! time at which the solution is desired
read(11,*) print_freq !
read(11,*) beta ! weight factor for implicit method
! 0-FTCS, 1-BTCS, 0.5-Crank-Nicolson
c
if (beta < 5d-1) then
Fourier = 1d0 / (4*(1-2*beta))
write(*,*)''
print*, 'The mean Fourier number supplied for the given value of'
print*, 'beta =',beta, 'must be less than', Fourier
write(*,*)''
print*, 'Press Enter to continue'
read(*,*)
endif
c
c calculation of time-step dt based on the stability condition
dx = alx /(m-2)
dy = aly /(n-2)
dt = 2*Fo *(dx**2 * dy**2)/(dx**2 + dy**2) /alpha
maxntimestp = tfinal /dt
dt = tfinal /maxntimestp ! revised time-step
Fox = alpha*dt/dx**2
Foy = alpha*dt/dy**2
Fo = (Fox + Foy)/2d0 ! revised mean Fourier number
if (tfinal < dt) then
print*, 'Warning: final time is less than time-step!'
stop
endif
c
return
end
*
***SUBROUTINE: WRITING OUT THE BASIC DATA USED FOR COMPUTATION**********
*
subroutine WRITE_OUT
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/boundaryvalues.inc'
c
if (boundary_type == 1) then
write(21,*) 'Dirichlet bc is specified on all boundaries'
elseif (boundary_type == 2) then
write(21,*) 'Dirichlet on W,E,N Neumann on S'
endif
write(21,*) 'Thermal diffusivity', alpha
write(21,*) 'Mean grid Fourier number, Fo =', Fo
write(21,*) 'Length of the domain in x-direction', alx
write(21,*) 'Length of the domain in y-direction', aly
write(21,*) 'Number of control volumes:', nx, ' X', ny
write(21,*) 'dx =', dx
write(21,*) 'dy =', dy
write(21,*) 'dt =', dt
write(21,*) 'No. of Time step = ', ntimestp
write(21,*) 'Final time = ', time
write(21,*) 'convergence criterion for steady-state solution:',
$tolstdy
write(21,*) 'Weight factor for the implicit scheme = ', beta
c
if (beta == 0d0) then
! nothing
elseif (solver == 1) then
write(21,*) 'The solver used is PSOR'
elseif (solver == 2) then
write(21,*) 'The solver used is LSOR'
endif
c
return
end
*
****SUBROUTINE: GRID****************************************************
*
SUBROUTINE GRID
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'arrays/grid.inc'
c
c m-2 = number of control volumes the in x-direction
c n-2 = number of control volumes the in y-direction
dx = alx /(m-2)
x(1) = 0
x(2) = x(1) + dx/2d0
do i = 3, m-1
x(i) = x(i-1) + dx
enddo
x(m) = x(m-1) + dx/2d0
c
dy = aly /(n-2)
y(1) = 0
y(2) = y(1) + dy/2d0
do i = 3, n-1
y(i) = y(i-1) + dy
enddo
y(n) = y(n-1) + dy/2d0
c
return
end
*
***SUBROUTINE: BOUNDARY CONDITIONS FOR TEMPERATURE**********************
*
SUBROUTINE BOUNDARY_COND
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/boundaryvalues.inc'
include 'scalars/parameters.inc'
include 'arrays/temperature.inc'
include 'arrays/grid.inc'
c
c defining the boundary conditions
include 'bc/bc.inc'
c
c boundary_type: 1 - dirichlet: all
c boundary_type: 2 - dirichlet: W,E,N; neumann: S;
c
if (boundary_type == 1) then
c west and east
do j = 2, n-1
TMP(1,j) = BTW(x(1),y(j))
TMP(m,j) = BTE(x(m),y(j))
enddo
c south and north
do i = 2, m-1
TMP(i,1) = BTS(x(i),y(1))
TMP(i,n) = BTN(x(i),y(n))
enddo
c corner boundary nodes
TMP(1,1) = BTW(x(1),y(1))
TMP(m,1) = BTE(x(m),y(1))
TMP(1,n) = BTW(x(1),y(n))
TMP(m,n) = BTE(x(m),y(n))
elseif (boundary_type == 2) then
c west and east
do j = 2, n-1
TMP(1,j) = BTW(x(1),y(j))
TMP(m,j) = BTE(x(m),y(j))
enddo
c south and north
do i = 2, m-1
qs = BQS(x(i),y(1))
TMP(i,n) = BTN(x(i),y(n))
enddo
c corner boundary nodes
TMP(1,1) = BTW(x(1),y(1))
TMP(m,1) = BTE(x(m),y(1))
TMP(1,n) = BTW(x(1),y(n))
TMP(m,n) = BTE(x(m),y(n))
endif
c
return
end
*
***SUBROUTINE: INITIAL CONDITIONS FOR TEMPERATURE**********************
*
SUBROUTINE INITIAL_COND
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/parameters.inc'
include 'arrays/temperature.inc'
include 'arrays/grid.inc'
c
c defining the initial conditions
TINIT(x,y,t) = 0 ! dsin(pi*x/alx)
c
do j = 2, n-1
do i = 2, m-1
TMP(i,j) = TINIT(x(i),y(j),0)
enddo
enddo
c
return
end
*
***SUBROUTINE: CONSTRUCTION OF SYSTEM OF EQUATION AND ITS SOLUTION******
*
subroutine SOLVE
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/boundaryvalues.inc'
include 'scalars/parameters.inc'
include 'arrays/grid.inc'
include 'arrays/temperature.inc'
include 'arrays/coeff.inc'
doubleprecision TN(0:m1,0:n1), maxdif
c
include 'arithmeticfunctions/source.inc' ! source.inc file
! contains the definition of source term
c
c storing the values of temperature at nth time level
TN = TMP
c
ntimestp = 0
maxdif = 1d9
time = 0d0
write(34,101) time ! continuous time
write(35,101) time ! periodic time
monitorx = nx/4
monitory = ny/4
write(36,*) TMP(monitorx,monitory)
c
write(37,*) 'Number of iterations by solver at each time step'
c
10 if ((maxdif > tolstdy).and.(ntimestp < maxntimestp)) then
c construction of coefficient matrix
do j = 3, n-2
do i = 3, m-2
AS(i,j) = -beta *Foy
AW(i,j) = -beta *Fox
AE(i,j) = -beta *Fox
AN(i,j) = -beta *Foy
AP(i,j) = 1d0 - (AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
enddo
enddo
c construction of RHS vector
do j = 3, n-2
do i = 3, m-2
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(TN(i-1,j) - 2*TN(i,j) + TN(i+1,j)) +
$ Foy*(TN(i,j-1) - 2*TN(i,j) + TN(i,j+1)) ) + S(x(i),y(j))*dt
enddo
enddo
c
c incorporation of boundary conditions:
c
c boundary_type: 1 - dirichlet: all
c boundary_type: 2 - dirichlet: W,E,N; neumann: S;
c
if (boundary_type == 1) then
c west
i = 2
do j = 3, n-2
AS(i,j) = -beta *Foy
AW(i,j) = -8d0/3d0 *beta *Fox
AE(i,j) = -4d0/3d0 *beta *Fox
AN(i,j) = -beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(8*TN(i-1,j) - 12*TN(i,j) + 4*TN(i+1,j))/3d0 +
$ Foy*(TN(i,j-1) - 2*TN(i,j) + TN(i,j+1)) ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AW(i,j)*TMP(i-1,j)
AW(i,j) = 0
enddo
c east
i = m-1
do j = 3, n-2
AS(i,j) = -beta *Foy
AW(i,j) = -4d0/3d0 *beta *Fox
AE(i,j) = -8d0/3d0 *beta *Fox
AN(i,j) = -beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(4*TN(i-1,j) - 12*TN(i,j) + 8*TN(i+1,j))/3d0 +
$ Foy*(TN(i,j-1) - 2*TN(i,j) + TN(i,j+1)) ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AE(i,j)*TMP(i+1,j)
AE(i,j) = 0
enddo
c south
j = 2
do i = 3, m-2
AS(i,j) = -8d0/3d0 *beta *Foy
AW(i,j) = -beta *Fox
AE(i,j) = -beta *Fox
AN(i,j) = -4d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(TN(i-1,j) - 2*TN(i,j) + TN(i+1,j)) +
$ Foy*(8*TN(i,j-1) - 12*TN(i,j) + 4*TN(i,j+1))/3d0 ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AS(i,j)*TMP(i,j-1)
AS(i,j) = 0
enddo
c north
j = n-1
do i = 3, m-2
AS(i,j) = -4d0/3d0 *beta *Foy
AW(i,j) = -beta *Fox
AE(i,j) = -beta *Fox
AN(i,j) = -8d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(TN(i-1,j) - 2*TN(i,j) + TN(i+1,j)) +
$ Foy*(4*TN(i,j-1) - 12*TN(i,j) + 8*TN(i,j+1))/3d0 ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AN(i,j)*TMP(i,j+1)
AN(i,j) = 0
enddo
c south-west
i = 2
j = 2
AS(i,j) = -8d0/3d0 *beta *Foy
AW(i,j) = -8d0/3d0 *beta *Fox
AE(i,j) = -4d0/3d0 *beta *Fox
AN(i,j) = -4d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(8*TN(i-1,j) - 12*TN(i,j) + 4*TN(i+1,j))/3d0 +
$ Foy*(8*TN(i,j-1) - 12*TN(i,j) + 4*TN(i,j+1))/3d0 ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AS(i,j)*TMP(i,j-1) - AW(i,j)*TMP(i-1,j)
AS(i,j) = 0
AW(i,j) = 0
c south-east
i = m-1
j = 2
AS(i,j) = -8d0/3d0 *beta *Foy
AW(i,j) = -4d0/3d0 *beta *Fox
AE(i,j) = -8d0/3d0 *beta *Fox
AN(i,j) = -4d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(4*TN(i-1,j) - 12*TN(i,j) + 8*TN(i+1,j))/3d0 +
$ Foy*(8*TN(i,j-1) - 12*TN(i,j) + 4*TN(i,j+1))/3d0 ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AS(i,j)*TMP(i,j-1) - AE(i,j)*TMP(i+1,j)
AS(i,j) = 0
AE(i,j) = 0
c north-west
i = 2
j = n-1
AS(i,j) = -4d0/3d0 *beta *Foy
AW(i,j) = -8d0/3d0 *beta *Fox
AE(i,j) = -4d0/3d0 *beta *Fox
AN(i,j) = -8d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(8*TN(i-1,j) - 12*TN(i,j) + 4*TN(i+1,j))/3d0 +
$ Foy*(4*TN(i,j-1) - 12*TN(i,j) + 8*TN(i,j+1))/3d0 ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AW(i,j)*TMP(i-1,j) - AN(i,j)*TMP(i,j+1)
AW(i,j) = 0
AN(i,j) = 0
c north-east
i = m-1
j = n-1
AS(i,j) = -4d0/3d0 *beta *Foy
AW(i,j) = -4d0/3d0 *beta *Fox
AE(i,j) = -8d0/3d0 *beta *Fox
AN(i,j) = -8d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(4*TN(i-1,j) - 12*TN(i,j) + 8*TN(i+1,j))/3d0 +
$ Foy*(4*TN(i,j-1) - 12*TN(i,j) + 8*TN(i,j+1))/3d0 ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AE(i,j)*TMP(i+1,j) - AN(i,j)*TMP(i,j+1)
AE(i,j) = 0
AN(i,j) = 0
elseif (boundary_type == 2) then
c west
i = 2
do j = 3, n-2
AS(i,j) = -beta *Foy
AW(i,j) = -8d0/3d0 *beta *Fox
AE(i,j) = -4d0/3d0 *beta *Fox
AN(i,j) = -beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(8*TN(i-1,j) - 12*TN(i,j) + 4*TN(i+1,j))/3d0 +
$ Foy*(TN(i,j-1) - 2*TN(i,j) + TN(i,j+1)) ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AW(i,j)*TMP(i-1,j)
AW(i,j) = 0
enddo
c east
i = m-1
do j = 3, n-2
AS(i,j) = -beta *Foy
AW(i,j) = -4d0/3d0 *beta *Fox
AE(i,j) = -8d0/3d0 *beta *Fox
AN(i,j) = -beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(4*TN(i-1,j) - 12*TN(i,j) + 8*TN(i+1,j))/3d0 +
$ Foy*(TN(i,j-1) - 2*TN(i,j) + TN(i,j+1)) ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AE(i,j)*TMP(i+1,j)
AE(i,j) = 0
enddo
c south
j = 2
do i = 3, m-2
AS(i,j) = 0
AW(i,j) = -beta *Fox
AE(i,j) = -beta *Fox
AN(i,j) = -beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(TN(i-1,j) - 2*TN(i,j) + TN(i+1,j)) +
$ Foy*(TN(i,j+1) - TN(i,j)) ) - qs*Foy*dy +
$ S(x(i),y(j))*dt
enddo
c north
j = n-1
do i = 3, m-2
AS(i,j) = -4d0/3d0 *beta *Foy
AW(i,j) = -beta *Fox
AE(i,j) = -beta *Fox
AN(i,j) = -8d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(TN(i-1,j) - 2*TN(i,j) + TN(i+1,j)) +
$ Foy*(4*TN(i,j-1) - 12*TN(i,j) + 8*TN(i,j+1))/3d0 ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AN(i,j)*TMP(i,j+1)
AN(i,j) = 0
enddo
c south-west
i = 2
j = 2
AS(i,j) = 0
AW(i,j) = -8d0/3d0 *beta *Fox
AE(i,j) = -4d0/3d0 *beta *Fox
AN(i,j) = -4d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(8*TN(i-1,j) - 12*TN(i,j) + 4*TN(i+1,j))/3d0 +
$ Foy*(TN(i,j+1) - TN(i,j)) ) - qs*Foy*dy +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AW(i,j)*TMP(i-1,j)
AW(i,j) = 0
c south-east
i = m-1
j = 2
AS(i,j) = 0
AW(i,j) = -4d0/3d0 *beta *Fox
AE(i,j) = -8d0/3d0 *beta *Fox
AN(i,j) = -4d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(4*TN(i-1,j) - 12*TN(i,j) + 8*TN(i+1,j))/3d0 +
$ Foy*(TN(i,j+1) - TN(i,j)) ) - qs*Foy*dy +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AE(i,j)*TMP(i+1,j)
AE(i,j) = 0
c north-west
i = 2
j = n-1
AS(i,j) = -4d0/3d0 *beta *Foy
AW(i,j) = -8d0/3d0 *beta *Fox
AE(i,j) = -4d0/3d0 *beta *Fox
AN(i,j) = -8d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(8*TN(i-1,j) - 12*TN(i,j) + 4*TN(i+1,j))/3d0 +
$ Foy*(4*TN(i,j-1) - 12*TN(i,j) + 8*TN(i,j+1))/3d0 ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AW(i,j)*TMP(i-1,j) - AN(i,j)*TMP(i,j+1)
AW(i,j) = 0
AN(i,j) = 0
c north-east
i = m-1
j = n-1
AS(i,j) = -4d0/3d0 *beta *Foy
AW(i,j) = -4d0/3d0 *beta *Fox
AE(i,j) = -8d0/3d0 *beta *Fox
AN(i,j) = -8d0/3d0 *beta *Foy
AP(i,j) = 1d0 -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = TN(i,j) + (1d0 - beta)* (
$ Fox*(4*TN(i-1,j) - 12*TN(i,j) + 8*TN(i+1,j))/3d0 +
$ Foy*(4*TN(i,j-1) - 12*TN(i,j) + 8*TN(i,j+1))/3d0 ) +
$ S(x(i),y(j))*dt
C(i,j) = C(i,j) - AE(i,j)*TMP(i+1,j) - AN(i,j)*TMP(i,j+1)
AE(i,j) = 0
AN(i,j) = 0
endif
c
c solution of system of algebraic equation
if (beta == 0d0) then ! no system of eqns to be solved
do j = 2, n-1
do i = 2, m-1
TMP(i,j) = C(i,j)
enddo
enddo
elseif (solver == 1) then
CALL SOR(TMP) ! point SOR
elseif (solver == 2) then
CALL LSOR(TMP) ! line SOR
endif
c
time = time + dt
ntimestp = ntimestp + 1
c
maxdif = 0
do j = 2, n-1
do i = 2, m-1
dif = TMP(i,j) - TN(i,j)
if (dabs(dif) > maxdif) then
maxdif = dabs(dif)
endif
enddo
enddo
c
c updating the temperature
TN = TMP
c
write(34,101) time ! continuous time
write(6,102) ntimestp, maxdif
c
c printing out transient temperature at selected point of domain
j = jcounter + 1
jcounter = mod(j,print_freq)
if (jcounter == 0) then
write(35,101) (time + dt) ! periodic time
write(36,*) TMP(monitorx,monitory)
endif
go to 10
endif
c
101 format(e14.7)
102 format(i7,5x,e14.7)
c
c calculation of boundary temperatures where Neumann BC is specified
if (boundary_type == 1) then
!
elseif (boundary_type == 2) then
c south
j = 1
do i = 1, m
TMP(i,j) = (9*TMP(i,j+1) - TMP(i,j+2) - 3*qs*dy) /8d0
enddo
endif
c
return
end
*
**SUBROUTINE: POINT SOR FOR TRIDIAGONAL MATRIX**************************
*
SUBROUTINE SOR(X)
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/solveparams.inc'
include 'arrays/coeff.inc'
doubleprecision X(0:m1,0:n1)
c
mm = m-1
nn = n-1
iter = 1
wt = wtsor
relresidual = 1d9 ! maximum norm of relative residual
c
c calculation of maximum norm of rhs vector
rhsnorm = 0
do j = 2, nn
do i = 2, mm
rhsnorm = max(rhsnorm, dabs(C(i,j)))
enddo
enddo
rhsnorm = rhsnorm + 1d-30 ! avoid overflow if rhsnorm is zero
c
10 if ((relresidual > tolsor).and.(iter < itermax)) then
do j = 2, nn
do i = 2, mm
xold = X(i,j)
xnew = (C(i,j)- AS(i,j)*X(i,j-1) - AW(i,j)*X(i-1,j)
$ - AE(i,j)*X(i+1,j) - AN(i,j)*X(i,j+1))/AP(i,j)
X(i,j) = wt*xnew + (1d0 - wt)*xold
enddo
enddo
c
c compute new max-norm of residual vector
residnorm = 0
do j = 2, nn
do i = 2, mm
res = C(i,j) -(AS(i,j)*X(i,j-1) + AW(i,j)*X(i-1,j) + AP(i,j)
$ *X(i,j) + AE(i,j)*X(i+1,j) + AN(i,j)*X(i,j+1))
residnorm = max(residnorm, dabs(res))
enddo
enddo
relresidual = residnorm /rhsnorm
* write(333,101) iter, relresidual
iter = iter + 1
go to 10
endif
c
if (relresidual > tolsor) then
write(21,*) 'Tolerance condition for SOR is not met'
endif
write(37,*) iter
*101 format('',3x,i4,4x,e10.3)
c
return
end
*
***SUBROUTINE: LINE SOR FOR NONADJACENT PENTADIAGONAL MATRIX************
*
subroutine LSOR(X)
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/solveparams.inc'
include 'arrays/coeff.inc'
doubleprecision X(0:m1,0:n1)
parameter( mn = max(m,n)+1 )
doubleprecision LO(0:mn),DI(0:mn),UP(0:mn),CO(0:mn),XX(0:mn)
c
mm = m-1
nn = n-1
iter = 1
wt = wtlsor
relresidual = 1d9 ! maximum norm of relative residual
c
c calculation of maximum norm of rhs vector
rhsnorm = 0
do j = 2, nn
do i = 2, mm
rhsnorm = max(rhsnorm, dabs(C(i,j)))
enddo
enddo
rhsnorm = rhsnorm + 1d-30 ! avoid overflow if rhsnorm is zero
c
10 if ((relresidual > tolsor).and.(iter < itermax)) then
do j = 2, nn
do i = 2, mm
LO(i) = wt*AW(i,j)
DI(i) = AP(i,j)
UP(i) = wt*AE(i,j)
CO(i) = wt*(C(i,j) - AS(i,j)*X(i,j-1) - AN(i,j)*X(i,j+1))
$ + (1 - wt)*AP(i,j)*X(i,j)
* write(333,*) AS(i,j), AN(i,j)
enddo
CALL TDMA(LO,DI,UP,CO,XX)
do i = 2, mm
X(i,j) = XX(i)
enddo
enddo
c
c compute new max-norm of residual vector
residnorm = 0
do j = 2, nn
do i = 2, mm
res = C(i,j) -(AS(i,j)*X(i,j-1) + AW(i,j)*X(i-1,j) + AP(i,j)
$ *X(i,j) + AE(i,j)*X(i+1,j) + AN(i,j)*X(i,j+1))
residnorm = max(residnorm, dabs(res))
enddo
enddo
relresidual = residnorm /rhsnorm
* write(333,101) iter, relresidual
iter = iter + 1
go to 10
endif
c
if (relresidual > tolsor) then
write(21,*) 'Tolerance condition for LSOR is not met'
endif
write(37,*) iter
*101 format('',3x,i4,4x,e10.3)
c
return
end
*
***SUBROUTINE: TDMA*****************************************************
*
subroutine TDMA(L,D,U,C,X)
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
parameter( mn = max(m,n)+1 )
doubleprecision L(0:mn),D(0:mn),U(0:mn),C(0:mn),X(0:mn)
doubleprecision P(-1:mn),Q(-1:mn)
c
P(1) = 0d0
Q(1) = 0d0
c
c forward elimination
do i = 2, m-1
denom = D(i) + L(i)*P(i-1)
P(i) = -U(i) /denom
Q(i) = (C(i) - L(i)*Q(i-1)) /denom
enddo
c
c back substitution
do i = m-1, 2, -1
X(i) = P(i)*X(i+1) + Q(i)
enddo
*
return
end
*
***SUBROUTINE: Printing-out the output data*****************************
*
SUBROUTINE PRINTOUT
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'arrays/temperature.inc'
include 'arrays/grid.inc'
c
write(31,101) (x(i), i=1,m)
write(32,101) (y(i), i=1,n)
write(33,102) ((TMP(i,j), i=1,m), j=1,n)
101 format(e13.6)
102 format(34(e13.6,1x))! format number = no. of CV in x-direction + 2
c
return
end
|
