* A. Salih, Dept. of Mechanical Engg., NIT - Trichy, India. *
*************************************************************************
* THIS PROGRAM SOLVES POISSON EQUATION IN A ONE-DIMENSIONAL DOMAIN *
* USING FINITE DIFFERENCE METHOD (UNIFORM GRID) AND THEN DETERMINE *
* GLOBAL ERROR BY COMPARING THE NUMERICAL AND ANALYTICAL SOLUTION *
*************************************************************************
*
program FD1
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/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
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,*) itermax ! maximum number of iterations (for SOR)
c
write(21,*) 'Length of the domain', alx
if (boundary_type == 1) then
write(21,*) 'Dirichlet bc is specified on all boundaries'
elseif (boundary_type == 2) then
write(21,*) 'West-side Dirichlet, east-side Neumann conditions'
elseif (boundary_type == 3) then
write(21,*) 'Neumann bc is specified on all boundaries'
endif
if (solver == 1) then
write(21,*) 'The solver used is TDMA'
elseif (solver == 2) then
write(21,*) 'The solver used is SOR'
endif
if ((boundary_type == 3).and.(solver == 1)) then
write(21,*)
write(21,*) 'WARNING: TDMA solver will not work with a "pure Ne
$umann problem".'
write(21,*) 'Try SOR solver.'
elseif (boundary_type == 3) then
write(21,*)
write(21,*) '1) No unique solution exists for a pure Neumann pr
$oblem.'
write(21,*) '2) Solutions can only be determined within an arbi
$trary constant.'
write(21,*) '3) For a pure Neumann problem, the source and the
$net boundary'
write(21,*) ' flux must be compatible.'
write(21,*)
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 n = no. of grid points (including boundary nodes) in x-direction
dx = alx /(n-1)
x(1) = 0
do i = 2, n
x(i) = x(i-1) + dx
enddo
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 - left; Neumann - right
c boundary_type: 3 - Neumann - all
c
if (boundary_type == 1) then
TMP(1) = tw
TMP(n) = te
elseif (boundary_type == 2) then
TMP(1) = tw
elseif (boundary_type == 3) then
! nothing
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
c construction of coefficient matrix
do i = 2, n-1
AW(i) = 1
AE(i) = 1
AP(i) = -(AW(i) + AE(i))
enddo
c
c construction of RHS vector
do i = 2, n-1
C(i) = -S(x(i))*dx**2
enddo
c
c incorporation of boundary conditions
if (boundary_type == 1) then
c west
C(2) = C(2) - AW(2)*TMP(1)
AW(2) = 0
c east
C(n-1) = C(n-1) - AE(n-1)*TMP(n)
AE(n-1) = 0
elseif (boundary_type == 2) then
c west
C(2) = C(2) - AW(2)*TMP(1)
AW(2) = 0
c east
AW(n-1) = AW(n-1) - (1d0/3d0)*AE(n-1)
AP(n-1) = AP(n-1) + (4d0/3d0)*AE(n-1)
C(n-1) = C(n-1) - (2d0/3d0)*AE(n-1) *qe *dx
AE(n-1) = 0
elseif (boundary_type == 3) then
c west
AE(2) = AE(2) - (1d0/3d0)*AW(2)
AP(2) = AP(2) + (4d0/3d0)*AW(2)
C(2) = C(2) + (2d0/3d0)*AW(2) *qw *dx
AW(2) = 0
c east
AW(n-1) = AW(n-1) - (1d0/3d0)*AE(n-1)
AP(n-1) = AP(n-1) + (4d0/3d0)*AE(n-1)
C(n-1) = C(n-1) - (2d0/3d0)*AE(n-1) *qe *dx
AE(n-1) = 0
endif
if (solver == 1) then
CALL TDMA(AW,AP,AE,C,TMP) ! TDMA
elseif (solver == 2) then
CALL SOR(TMP) ! SOR
endif
c
c calculation of heat fluxes at boundaries (flux is proportional to
c the negative of slope of the temperature curve)
if (boundary_type == 1) then
qw = -(-TMP(3) + 4*TMP(2) - 3*TMP(1)) /(2*dx)
qe = -(3*TMP(n) - 4*TMP(n-1) + TMP(n-2)) /(2*dx)
elseif (boundary_type == 2) then
qw = -(-TMP(3) + 4*TMP(2) - 3*TMP(1)) /(2*dx)
TMP(n) = (4*TMP(n-1) - TMP(n-2) + 2*qe*dx) /3d0
elseif (boundary_type == 3) then
TMP(1) = (4*TMP(2) - TMP(3) - 2*qw*dx) /3d0
TMP(n) = (4*TMP(n-1) - TMP(n-2) + 2*qe*dx) /3d0
endif
write(21,*) ''
write(21,*) 'The temperature at West boundary is:', TMP(1)
write(21,*) 'The temperature at East boundary is:', TMP(n)
write(21,*) 'The heat flux at West boundary is:', qw
write(21,*) 'The heat flux at East boundary is:', qe
*
return
end
*
***SUBROUTINE: TDMA*****************************************************
*
subroutine TDMA(L,D,U,C,X)
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
doubleprecision L(0:n1),U(0:n1),D(0:n1),C(0:n1),X(0:n1)
doubleprecision P(-1:n1),Q(-1:n1)
c
P(1) = 0d0
Q(1) = 0d0
c
c forward elimination
do i = 2, n-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 = n-1, 2, -1
X(i) = P(i)*X(i+1) + Q(i)
enddo
*
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:n1)
c
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 i = 2, nn
rhsnorm = max(rhsnorm, dabs(C(i)))
enddo
rhsnorm = rhsnorm + 1d-30 ! avoid overflow if rhsnorm is zero
c
c initial guess of the solution
do i = 2, nn
X(i) = 0
enddo
c
10 if ((relresidual > tolsor).and.(iter < itermax)) then
do i = 2, nn
xold = X(i)
xnew = ( C(i)- AW(i)*X(i-1) - AE(i)*X(i+1) ) /AP(i)
X(i) = wt*xnew + (1d0 - wt)*xold
enddo
c
c compute new max-norm of residual vector
residnorm = 0
do i = 2, nn
res = C(i) -( AW(i)*X(i-1) + AP(i)*X(i) + AE(i)*X(i+1) )
residnorm = max(residnorm, dabs(res))
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: 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'
c
include 'arithmeticfunctions/analyt_sol.inc' ! analyt_sol.inc file
! contains the arithmetic fn for analytical solution
c
open(unit=33,file='outputs/tanal.dat',status='unknown')
open(unit=34,file='outputs/errors.dat',status='unknown')
c
do i = 1, n
if (boundary_type == 1) then
tana = TA(x(i),tw,te)
elseif (boundary_type == 2) then
tana = TA(x(i),tw,qe)
elseif (boundary_type == 3) then
tana = TA(x(i),qw,qe)
endif
error = 0d0
if (dabs(tana) > 1d-12) then
error = (tana - TMP(i))/tana *100
endif
write(33,101) tana
write(34,101) error
enddo
101 format(e13.6)
*
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,n)
write(32,102) (TMP(i), i=1,n)
101 format(e13.6)
102 format(e13.6)
*
return
end
|
