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