Changes between Initial Version and Version 1 of cypress/Programming/HeatMassTransfer/ex3.f90

Timestamp:

05/14/15 15:12:10 (10 years ago)

Author:

cmaggio

Comment:

migrated content from ccs wiki

cypress/Programming/HeatMassTransfer/ex3.f90

@@ +1 @@
+= Ex3.f90 =
+{{{#!fortran
+!
+!  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
+}}}