* A. Salih, Dept. of Mechanical Engg., NIT - Trichy, India. *
*************************************************************************
* THIS PROGRAM SOLVES POISSON EQUATION IN A TWO-DIMENSIONAL DOMAIN *
* USING FINITE VOLUME METHOD (UNIFORM CV) AND THEN DETERMINE GLOBAL *
* ERROR BY COMPARING THE NUMERICAL AND ANALYTICAL SOLUTION *
*************************************************************************
*
program FV2
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/t.dat',status='unknown')
c
pi = 4*datan(1d0)
c
CALL READ_WRITE ! read-in and write-out data
close (unit=11)
CALL GRID ! setting up the grid points in the domain
CALL BOUNDARY_COND ! setting up boundary values
CALL SOLVE ! setting up system of equations
CALL ANALYT ! analytical solution
CALL PRINTOUT ! printing out the results
c
stop
end
*
***SUBROUTINE: READING IN AND WRITING OUT THE BASIC DATA****************
*
subroutine READ_WRITE
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, x-direction
read(11,*) aly ! actual length of domain, y-direction
read(11,*) boundary_type ! flag for type of boundary conditions
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,*) wtadisor ! relaxation factor for ADI SOR
read(11,*) itermax ! maximum number of iterations (for solver)
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'
elseif (boundary_type == 3) then
write(21,*) 'Neumann on W,E,N Dirichlet on S'
endif
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
if (solver == 1) then
write(21,*) 'The solver used is PSOR'
elseif (solver == 2) then
write(21,*) 'The solver used is LSOR'
elseif (solver == 3) then
write(21,*) 'The solver used is ADISOR'
endif
*
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 boundary_type: 3 - neumann: W,E,N; dirichlet: 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))
elseif (boundary_type == 3) then
c west and east
do j = 2, n-1
qw = BQW(x(1),y(j))
qe = BQE(x(m),y(j))
enddo
c south and north
do i = 2, m-1
TMP(i,1) = BTS(x(i),y(1))
qn = BQN(x(i),y(n))
enddo
c corner boundary nodes
TMP(1,1) = BTS(x(1),y(1))
TMP(m,1) = BTS(x(m),y(1))
endif
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'
c
include 'arithmeticfunctions/source.inc' ! source.inc file
! contains the definition of source term
c
hx = dy/dx
hy = dx/dy
c
c construction of coefficient matrix for interior nodes
do j = 2, n-1
do i = 2, m-1
AS(i,j) = hy
AW(i,j) = hx
AE(i,j) = hx
AN(i,j) = hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
enddo
enddo
c
c construction of RHS vector
do j = 2, n-1
do i = 2, m-1
C(i,j) = -S(x(i),y(j))*dx*dy
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 boundary_type: 3 - neumann: W,E,N; dirichlet: S;
c
if (boundary_type == 1) then
c west
i = 2
do j = 3, n-2
AW(i,j) = 8d0/3d0 *hx
AE(i,j) = 4d0/3d0 *hx
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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
AW(i,j) = 4d0/3d0 *hx
AE(i,j) = 8d0/3d0 *hx
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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 *hy
AN(i,j) = 4d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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 *hy
AN(i,j) = 8d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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 *hy
AW(i,j) = 8d0/3d0 *hx
AE(i,j) = 4d0/3d0 *hx
AN(i,j) = 4d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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 *hy
AW(i,j) = 4d0/3d0 *hx
AE(i,j) = 8d0/3d0 *hx
AN(i,j) = 4d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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 *hy
AW(i,j) = 8d0/3d0 *hx
AE(i,j) = 4d0/3d0 *hx
AN(i,j) = 8d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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 *hy
AW(i,j) = 4d0/3d0 *hx
AE(i,j) = 8d0/3d0 *hx
AN(i,j) = 8d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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
AW(i,j) = 8d0/3d0 *hx
AE(i,j) = 4d0/3d0 *hx
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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
AW(i,j) = 4d0/3d0 *hx
AE(i,j) = 8d0/3d0 *hx
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) + qs*dx
enddo
c north
j = n-1
do i = 3, m-2
AS(i,j) = 4d0/3d0 *hy
AN(i,j) = 8d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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 *hx
AE(i,j) = 4d0/3d0 *hx
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) + qs*dx - 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 *hx
AE(i,j) = 8d0/3d0 *hx
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) + qs*dx - AE(i,j)*TMP(i+1,j)
AE(i,j) = 0
c north-west
i = 2
j = n-1
AS(i,j) = 4d0/3d0 *hy
AW(i,j) = 8d0/3d0 *hx
AE(i,j) = 4d0/3d0 *hx
AN(i,j) = 8d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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 *hy
AW(i,j) = 4d0/3d0 *hx
AE(i,j) = 8d0/3d0 *hx
AN(i,j) = 8d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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 == 3) then
c west
i = 2
do j = 3, n-2
AW(i,j) = 0
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) + qw*dy
enddo
c east
i = m-1
do j = 3, n-2
AE(i,j) = 0
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) - qe*dy
enddo
c south
j = 2
do i = 3, m-2
AS(i,j) = 8d0/3d0 *hy
AN(i,j) = 4d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
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, n-2
AN(i,j) = 0
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) - qn*dx
enddo
c south-west
i = 2
j = 2
AS(i,j) = 8d0/3d0 *hy
AW(i,j) = 0
AN(i,j) = 4d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) - AS(i,j)*TMP(i,j-1) + qw*dy
AS(i,j) = 0
c south-east
i = m-1
j = 2
AS(i,j) = 8d0/3d0 *hy
AE(i,j) = 0
AN(i,j) = 4d0/3d0 *hy
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) - AS(i,j)*TMP(i,j-1) - qe*dy
AS(i,j) = 0
c north-west
i = 2
j = n-1
AW(i,j) = 0
AN(i,j) = 0
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) + qw*dy - qn*dx
c north-east
i = m-1
j = n-1
AE(i,j) = 0
AN(i,j) = 0
AP(i,j) = -(AS(i,j) + AW(i,j) + AE(i,j) + AN(i,j))
C(i,j) = C(i,j) - qe*dy - qn*dx
endif
c
c solution of system of algebraic equation
if (solver == 1) then
CALL SOR(TMP) ! point SOR
elseif (solver == 2) then
CALL LSOR(TMP) ! line SOR
elseif (solver == 3) then
CALL ADISOR(TMP) ! line SOR with alternating direction
endif
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
elseif (boundary_type == 3) then
c west
i = 1
do j = 1, n
TMP(i,j) = (9*TMP(i+1,j) - TMP(i+2,j) - 3*qw*dx) /8d0
enddo
c east
i = m
do j = 1, n
TMP(i,j) = (9*TMP(i-1,j) - TMP(i-2,j) + 3*qe*dx) /8d0
enddo
c north
j = n
do i = 1, m
TMP(i,j) = (9*TMP(i,j-1) - TMP(i,j-2) + 3*qn*dy) /8d0
enddo
endif
* write(21,*) ''
* write(21,*) 'The heat flux at West boundary is:', qw
* write(21,*) 'The heat flux at East boundary is:', qe
* write(21,*) 'The heat flux at South boundary is:', qs
* write(21,*) 'The heat flux at North boundary is:', qn
*
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
c initial guess of the solution
do j = 2, nn
do i = 2, mm
X(i,j) = 0
enddo
enddo
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(21,*) 'Number of iterations performed by solver', iter
*101 format('',3x,i4,4x,e10.3)
*
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
c initial guess of the solution
do j = 2, nn
do i = 2, mm
X(i,j) = 0
enddo
enddo
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(21,*) 'Number of iterations performed by solver', iter
*101 format('',3x,i4,4x,e10.3)
*
return
end
*
***SUBROUTINE: ADISOR FOR NONADJACENT PENTADIAGONAL MATRIX****************
*
SUBROUTINE ADISOR(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 = wtadisor
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
c initial guess of the solution
do j = 2, nn
do i = 2, mm
X(i,j) = 0
enddo
enddo
c
10 if ((relresidual > tolsor).and.(iter < itermax)) then
c sweep by rows
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)
enddo
CALL TDMA(LO,DI,UP,CO,XX)
do i = 2, mm
X(i,j) = XX(i)
enddo
enddo
c
c sweep by columns
do i = 2, mm
do j = 2, nn
LO(j) = wt*AS(i,j)
DI(j) = AP(i,j)
UP(j) = wt*AN(i,j)
CO(j) = wt*(C(i,j) - AW(i,j)*X(i-1,j) - AE(i,j)*X(i+1,j))
$ + (1 - wt)*AP(i,j)*X(i,j)
enddo
CALL TDMA(LO,DI,UP,CO,XX)
do j = 2, nn
X(i,j) = XX(j)
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 ADISOR is not met'
endif
write(21,*) 'Number of iterations performed by solver', iter
*101 format('',3x,i4,4x,e10.3)
*
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: ANALYTICAL SOLUTION**************************************
*
subroutine ANALYT
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'
doubleprecision TANA(0:m1,0:n1), ER(0:m1,0:n1)
c
include 'arithmeticfunctions/analyt_sol.inc' ! analyt_sol.inc file
! contains the arithmetic fn for analytical solution
c
open(unit=34,file='outputs/tanal.dat',status='unknown')
open(unit=35,file='outputs/errors.dat',status='unknown')
c
if (boundary_type == 1) then
do j = 1, n
do i = 1, m
sum = 0
do kk = 1, 100
sum = sum + TA(x(i),y(j),kk)
enddo
TANA(i,j) = TMP(m/2,1) + (TMP(m/2,n) - TMP(m/2,1)) *2/pi*sum
ER(i,j) = 0d0
if (dabs(TANA(i,j)) > 1d-12) then
ER(i,j) = (TANA(i,j) - TMP(i,j))/TANA(i,j) *100
endif
enddo
enddo
endif
write(34,101) ((TANA(i,j), i=1,m), j=1,n)
write(35,101) ((ER(i,j), i=1,m), j=1,n)
101 format(66(e13.6,1x))
*
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(66(e13.6,1x))! format number = no. of CV in x-direction + 2
*
return
end
|
