program laplace
interface
subroutine initial_u(u,nx,ny)
integer nx,ny
double precision, dimension (1:nx,1:ny) :: u
end subroutine
subroutine get_coordinates(xcoord,ycoord,nx,ny)
integer nx,ny
double precision, dimension(1:nx) :: xcoord
double precision, dimension(1:ny) :: ycoord
end subroutine
subroutine put_boundary_value(u,nx,ny,xcoord,ycoord)
double precision, dimension(1:nx,1:ny) :: u
integer nx,ny
double precision, dimension(1:nx) :: xcoord
double precision, dimension(1:ny) :: ycoord
end subroutine
subroutine jacobi(u,uold,nx,ny)
integer nx,ny
double precision, dimension(1:nx,1:ny) :: u,uold
end subroutine
subroutine anal_error(error,u,nx,ny,xcoord,ycoord)
double precision error
integer nx,ny
double precision, dimension(1:nx, 1:ny):: u
double precision, dimension(1:nx) :: xcoord
double precision, dimension(1:ny) :: ycoord
end subroutine
subroutine interior_avg(u,nx,ny)
integer nx,ny
double precision, dimension(1:nx,1:ny) :: u
end subroutine
end interface
parameter(nx = 21)
parameter(ny = 11)
integer i,count
double precision error
double precision, dimension(1:nx, 1:ny):: u,uold
double precision, dimension(1:nx) :: xcoord
double precision, dimension(1:ny) :: ycoord
write(*,*) 'MESH RESOLUTION INFO:'
write(*,*) 'dx is ',2.0d0/dble(nx-1)
write(*,*) 'dy is ',1.0d0/dble(ny-1)
! GET COORDINATES
call get_coordinates(xcoord,ycoord,nx,ny)
! SET U to 0, UOLD to 1.0
u = 0.0d0
uold = 1.0d0
count = 0
! PUT BOUNDARY VALUES ON THE BOUNDARIES
call put_boundary_value(u,nx,ny,xcoord,ycoord)
call interior_avg(u,nx,ny)
error = 1.0d0
do while( error > 3.0d-06 )
count = count + 1
call jacobi(u,uold,nx,ny)
error = maxval(abs(u-uold))
! write(*,*) 'infinity norm between u and uold after step ', &
! count,' is ',error
enddo
write(*,*) 'PROGRAM COMPLETE!'
write(*,*) 'total number of Jacobi iteration steps: ',count
write(*,*) ' '
write(*,*) 'infinity norm between u and uold for last ', &
'Jacobi step is',error
call anal_error(error,u,nx,ny,xcoord,ycoord)
write(*,*) ' '
write(*,*) 'infinity norm between u and the analytical ', &
'solution is ',error
end
!********************************************************************
!** INTERIOR_AVG
!********************************************************************
!*** This routine calculates the average of the boundary of the
!*** input grid u, and sets the interior points of the grid
!*** to that value
!********************************************************************
subroutine interior_avg(u,nx,ny)
integer nx,ny
double precision, dimension(1:nx,1:ny) :: u
! local var
double precision avg
avg = 0.0d0
do i=1,nx
avg = avg + u(i,1) + u(i,ny)
enddo
do i=1,ny
avg = avg + u(1,i) + u(nx,i)
enddo
avg = avg / dble(2 * (nx + ny))
do i=2,nx-1
do j=2,ny-1
u(i,j) = avg
enddo
enddo
end
!********************************************************************
!*** ANAL_ERROR
!********************************************************************
!*** This routine returns the infinity norm between the input
!*** grid u and the known analytic solution (sampled at the grid points)
!********************************************************************
subroutine anal_error(error,u,nx,ny,xcoord,ycoord)
double precision error
integer nx,ny
double precision, dimension(1:nx, 1:ny):: u
double precision, dimension(1:nx) :: xcoord
double precision, dimension(1:ny) :: ycoord
! local variables
double precision, dimension(1:nx,1:ny) :: utemp
do i=1,nx
do j=1,ny
call calc_u(utemp(i,j),xcoord(i),ycoord(j))
enddo
enddo
error = maxval(abs(u-utemp))
end
!********************************************************************
!*** JACOBI
!********************************************************************
!*** This routine performs one jacobi iteration on the
!*** input grid u. The new values are returned in u, and
!*** uold on output contains the input u.
!********************************************************************
subroutine jacobi(u,uold,nx,ny)
integer nx,ny
double precision, dimension(1:nx,1:ny) :: u,uold
uold = u
do i=2,nx-1
do j=2,ny-1
u(i,j) = (uold(i-1,j) + uold(i+1,j) + uold(i,j-1) + &
uold(i,j+1) )/4.0d0
enddo
enddo
end
!********************************************************************
!*** PUT_BOUNDARY_VALUE
!********************************************************************
!*** This routine puts the boundary values on the boundary of u
!********************************************************************
subroutine put_boundary_value(u,nx,ny,xcoord,ycoord)
interface
subroutine calc_u( uval,x,y)
double precision uval,x,y
end subroutine
end interface
double precision, dimension(1:nx,1:ny) :: u
integer nx,ny
double precision, dimension(1:nx) :: xcoord
double precision, dimension(1:ny) :: ycoord
do i=1,nx
call calc_u( u(i,1) ,xcoord(i),0.0d0)
call calc_u( u(i,ny), xcoord(i),1.0d0)
enddo
do i=1,ny
call calc_u( u(1,i), 0.0d0, ycoord(i))
call calc_u( u(nx,i), 2.0d0, ycoord(i))
enddo
end
!*******************************************************************
!*** CALC_U
!*******************************************************************
!*** This routine calculates an analytic function, and returns it
!*** in uval
!*******************************************************************
subroutine calc_u( uval,x,y)
double precision uval,x,y
uval = -4.7 + 0.2 * x - 0.3 * y + 0.09 * (x * x - y * y)
uval = uval + log(sqrt( ((x+0.2)**2.0) + ((y+0.2)**2.0) ) )
! write(*,*) uval,x,y
end subroutine
!********************************************************************
!*** GET_COORDINATES
!********************************************************************
!*** This routine calculates the coordinates for the X-Y mesh
!********************************************************************
subroutine get_coordinates(xcoord,ycoord,nx,ny)
integer nx,ny
double precision, dimension(1:nx) :: xcoord
double precision, dimension(1:ny) :: ycoord
do i = 1,nx
xcoord(i) = 2.0 * dble(i-1)/dble(nx-1)
enddo
do i = 1,ny
ycoord(i) = 1.0 * (dble(i-1)/dble(ny-1))
enddo
end
!********************************************************************
subroutine initial_u(u,nx,ny)
integer nx,ny
double precision, dimension(1:nx,1:ny) :: u
u = 0.0
end