= Ex17.f90 = {{{#!fortran ! ! 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 }}}