0
c---|
c..A POINT/LINE ITERATIVE SOLVER FOR ELLIPTIC PDEs |
c Course: AE305 | |
c---|
program ELLIPTIC
parameter (imx=101, jmx=201)
common /var/ beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)
common /stress/ tauzx(imx,jmx),tauzy(imx,jmx)
real error, error_tol, error0
integer i,j,k
c..Read the input data, generate the grid data and initialize the solution
call INIT
write(*,*) 'enter error_tol'
read(*,*) error_tol
c..Start the iterative solution loop
k = 0
m=alog10(error/error0)
DO WHILE (m.gt.error_tol)
k = k + 1
c..Point iterative solutions
call ITERATIVE
c..Line iterative solution
c call ADI(... )
c..Update phik array and evaluate residual
rl2=0.
do j = 1,jmax
do i = 1,imax
rl2 = rl2 + (phikp(i,j) - phik(i,j))**2
phik(i,j) = phikp(i,j)
enddo
enddo
rl2 = sqrt(rl2/(imax*jmax))
print*, 'Residual@k =',k,rl2
c... error definition
error=0
do j=2,jmax
do i=2,imax
error=error+(phik(i,j)-phik(i,j))**2
enddo
enddo
error=error/float(imax*jmax)
error=sqrt(error)
c...
call sSTRESS
c..Output intermediate solutions
if( mod(k,iostep).eq.0 .and. k.ne.kmx) call IO(k)
ENDDO
call IO(k)
stop
end
c---
subroutine INIT
parameter (imx=101, jmx=201)
common /var/ beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)
common /stress/ tauzx(imx,jmx),tauzy(imx,jmx)
xl=0.05
yl=0.10
Torque = 100
c..set the second moment of area, J and the rhs
sma = xl*yl*(xl**2+yl**2)/3.
rhs = -2*Torque/sma
c..Read input parameters like imax, jmax
print*, 'Enter Imax and Jmax :'
read(*,*) imax, jmax
dx = xl/float(imax-1)
dy = yl/float(jmax-1)
beta2 = (dx/dy)**2
c..Grid generation
do j=1,jmax
do i=1,imax
x(i,j)= dx*float(i-1)
y(i,j)= dy*float(j-1)
enddo
enddo
c..Initialize the phi arrays
do j = 1,jmax
do i = 1,imax
phik(i,j) = 0.
tauzx(i,j)=0.
tauzy(i,j)=0.
enddo
enddo
error0=1.
error=error0
return
end
c---
subroutine ITERATIVE
parameter (imx=101, jmx=201)
common /var/ beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)
c..Solve for phi^k+1
do j = 2,jmax-1
do i = 2,imax-1
c..Implement, Jacobi, Gauss-Seidel and SOR methods
phikp(i,j) =(1./2.*(1+beta2))*(-rhs+phik(i-1,j)+phik(i+1,j)+beta2*
&(phik(i,j-1)+phik(i,j+1)))
enddo
enddo
return
end
c---
subroutine sSTRESS
parameter (imx=101, jmx=201)
common /var/beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)
common /stress/ tauzx(imx,jmx),tauzy(imx,jmx)
do j=2,jmax-1
do i=2,imax-1
phik(i,j)=phikp(i,j)
tauzx(i,j)=-(phik(i+1,j)-phik(i-1,j))/2.*dx
tauzy(i,j)=(phik(i,j+1)-phik(i,j-1))/2.*dy
enddo
enddo
return
end
c---
subroutine IO(k)
c..Output solution in tecplot format
parameter (imx=101, jmx=201)
common /var/ beta2,rhs, phik(imx,jmx), phikp(imx,jmx)
common /grid/ imax, jmax,dx,dy,x(imx,jmx),y(imx,jmx)
common /stress/ tauzx(imx,jmx),tauzy(imx,jmx)
character fname*32,string*8,ext*5
write(string,'(f8.5)') float(k)/100000
read(string,'(3x,a5)') ext
fname = 'phi-'//ext//'.tec'
open(1,file=fname, form='formatted')
write(1,*) ' variables="x","y","phi","tauzx","tauzy" '
write(1,*) ' zone i=',imax, ', j=',jmax
do j = 1,jmax
do i = 1,imax
write(1,*) x(i,j),y(i,j),phik(i,j),tauzx(i,j),tauzy(i,j)
enddo
enddo
close(1)
return
end
c---
c subroutine THOMAS(AA,BB,CC,FF,n)
---
! Solution of a tridiagonal system of n equations of the form
! aa(i)*x(i-1) + bb(i)*x(i) + cc(i)*x(i+1) = ff(k) for i=1,n
! the solution x(k) is stored in ff(k)
! aa(1) and cc(n) is not used.
! aa,bb,cc,ff are arrays to be filled by the caller program
!______________________________________________________________
c real aa(1), bb(1), cc(1), ff(1)
c
c n2 = 2
c n1n = n+1
c..forward sweep
c
c bb(1) = 1./bb(1)
c aa(1) = ff(1)*bb(1)
c do k = n2,n
c k1 = k-1
c cc(k1) = cc(k1)*bb(k1)
c bb(k) = 1./(bb(k)-aa(k)*cc(k1))
c aa(k) = (ff(k)-aa(k)*aa(k1))*bb(k)
c enddo
c
c..back substitution
c
c ff(n) = aa(n)
c do k1 = n2,n
c k = n1n-k1
c ff(k) = aa(k)-cc(k)*ff(k+1)
c enddo
c
c return
c end
Tümünü Göster