Ex3.f90

!
!  Large Scale Computing
!  Heat/Mass Transfer
!  ex3.f90
!
program diffusion1d
  implicit none
  double precision :: a,b,d
  double precision,dimension(:),allocatable :: phi,phi0
  double precision,dimension(:),allocatable :: rhs
  double precision :: dx
  double precision :: err,merr
  double precision,parameter :: tol = 1.0e-8
  integer :: i,num,itc
 
  character(len=32) :: arg
 
  ! Newer version of Fortran can get commadline arguments
  do i=1,4
     call get_command_argument(i,arg)
     if (len_trim(arg) == 0) then
        call get_command_argument(0,arg)
        print *,arg," [num] [A] [B] [D]"
        call exit(-1)
     end if
 
     ! get inputed depth
     if (i == 1) read(arg,*) num
     if (i == 2) read(arg,*) a
     if (i == 3) read(arg,*) b
     if (i == 4) read(arg,*) d
  end do
 
  print *,"num=",num,"A=",a,"B=",b,"D=",d
 
  ! Memory Allocation
  allocate(phi(num))
  allocate(phi0(num))
  allocate(rhs(num))
 
  ! Setup 
  dx = 1.d0 / dble(num - 1)
  do i = 1,num
     rhs(i) = -dx * dx / d * (dx * dble(i))
     phi(i) = 0.d0
     phi0(i) = 0.d0
  end do
 
  ! Boundary Condition
  phi(1) = a
  phi0(1) = a
  phi(num) = b
  phi0(num) = b
 
  ! Solve with Jacobi Method
  itc = 0;
  do 
     ! update phi 
     do i = 2,num-1
        phi(i) = -0.5 * (rhs(i) - phi0(i + 1) - phi0(i - 1))
     end do
 
    ! Check Convergence 
     merr = 0.d0
     do i = 1, num
        err = abs(phi(i) - phi0(i))
        if (merr < err) then
           merr = err
        end if
        phi0(i) = phi(i)
     end do
 
     itc = itc + 1
     if (merr < tol) exit
     if (mod(itc,10000) == 0) print *,"itc=",itc,"err=",merr
  end do
 
  print *,"Number of Iteration=",itc
 
  ! Output Result 
  open(10,FILE='res.dat')
  do i = 1,num
     write(10,'(2e16.8)') dx * dble(i-1),phi(i)
  end do
  close(10)
 
  ! END : Deallocate Memory
  deallocate(phi)
  deallocate(phi0)
  deallocate(rhs)
  print *,"Done"
end program diffusion1d