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