Ex17.f90
!
! Large Scale Computing
! 2D Heat/Mass Transfer
! ex17.f90 : Numerical Analysis
! Solve for
! d2c/dx2 + d2c/dy2 + 1 = 0
! With the boundary conditions of c = 0
! along lines of x=1,-1, and y=1, and -1
!
program diffusion2d
implicit none
integer, parameter :: NUM=20 ! number of points for x and y
double precision, parameter :: TOL=1.0e-12 ! convergence tolerance
integer :: i,j,k,n,ie,iw,in,is
double precision :: x,y,dx,w,err
double precision, dimension(NUM * NUM) :: cnc ! define 1D array of length NUM * NUM
!assuming dx = dy : domain is 2x2
dx = 2.d0 / dble(NUM - 1)
! Initialize
do n = 1, NUM * NUM
cnc(n) = 0.0;
end do
!computing for C with Gauss-Seidel
k = 0
do
err = 0.d0
do n = 1, NUM * NUM
i = mod(n - 1,NUM) + 1
j = (n - 1) / NUM + 1
if (i == 1) cycle !skip left end
if (i == NUM) cycle !skip right end
if (j == 1) cycle !skip bottom end
if (j == NUM) cycle !skip top end
ie = n + 1;
iw = n - 1;
in = n + NUM;
is = n - NUM;
w = 0.25 * (dx * dx + cnc(ie) + cnc(iw) + cnc(in) + cnc(is))
err = err + (w - cnc(n)) * (w - cnc(n))
cnc(n) = w
end do
k = k + 1
if (err < TOL) exit
end do
print *,"# of Iteration=",k
! Output Result
open(10,FILE='res.dat')
do i = 1,NUM
do j = 1,NUM
x = -1.0 + 2.0 * dble(i-1) / dble(NUM - 1)
y = -1.0 + 2.0 * dble(j-1) / dble(NUM - 1)
write(10,'(3e16.8)') x,y,cnc(i + NUM * (j - 1))
end do
end do
close(10)
print *,"Done"
end program diffusion2d