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


Ignore:
Timestamp:
05/14/15 15:19:48 (9 years ago)
Author:
cmaggio
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • cypress/Programming/HeatMassTransfer/ex17.f90

    v1 v1  
     1= Ex17.f90 =
     2{{{#!fortran
     3!
     4!  Large Scale Computing
     5!  2D Heat/Mass Transfer
     6!  ex17.f90 : Numerical Analysis
     7!  Solve for
     8!  d2c/dx2 + d2c/dy2 + 1 = 0
     9!  With the boundary conditions of c = 0
     10!  along lines of x=1,-1, and y=1, and -1
     11!
     12program diffusion2d
     13  implicit none
     14  integer, parameter ::  NUM=20  ! number of points for x and y
     15  double precision, parameter :: TOL=1.0e-12 ! convergence tolerance
     16  integer ::  i,j,k,n,ie,iw,in,is
     17  double precision :: x,y,dx,w,err
     18  double precision, dimension(NUM * NUM) :: cnc ! define 1D array of length NUM * NUM
     19 
     20  !assuming dx = dy : domain is 2x2
     21  dx = 2.d0 / dble(NUM - 1)
     22 
     23  ! Initialize
     24  do n = 1, NUM * NUM
     25     cnc(n) = 0.0;
     26  end do
     27 
     28  !computing for C with Gauss-Seidel
     29  k = 0
     30  do
     31    err = 0.d0
     32    do n = 1, NUM * NUM
     33       i = mod(n - 1,NUM) + 1
     34       j = (n - 1) / NUM + 1
     35       if (i == 1) cycle !skip left end
     36       if (i == NUM) cycle !skip right end
     37       if (j == 1) cycle !skip bottom end
     38       if (j == NUM) cycle !skip top end
     39 
     40       ie = n + 1;
     41       iw = n - 1;
     42       in = n + NUM;
     43       is = n - NUM;
     44       w = 0.25 * (dx * dx + cnc(ie) + cnc(iw) + cnc(in) + cnc(is))
     45       err = err + (w - cnc(n)) * (w - cnc(n))
     46       cnc(n) = w
     47    end do
     48    k = k + 1
     49    if (err < TOL) exit
     50 end do
     51 print *,"# of Iteration=",k
     52 
     53 ! Output Result
     54 open(10,FILE='res.dat')
     55 do i = 1,NUM
     56    do j = 1,NUM
     57       x = -1.0 + 2.0 * dble(i-1) / dble(NUM - 1)
     58       y = -1.0 + 2.0 * dble(j-1) / dble(NUM - 1)
     59       write(10,'(3e16.8)') x,y,cnc(i + NUM * (j - 1))
     60    end do
     61 end do
     62 close(10)
     63 
     64 print *,"Done"
     65end program diffusion2d
     66}}}