The code below implements Jacobi iteration for solving the linear system arising from the steady state heat equation with a simple application of parallel do loops using OpenMP.
Compare to:
The code:
! $CLASSHG/codes/openmp/jacobi1.f90
program jacobi1
use omp_lib
implicit none
integer :: n, nthreads
real(kind=8), dimension(:), allocatable :: x,u,uold,f
real(kind=8) :: alpha, beta, dx, tol, dumax
real(kind=8), intrinsic :: exp
integer :: i,iter,maxiter
! Specify number of threads to use:
nthreads = 2
!$ call omp_set_num_threads(nthreads)
!$ print "('Using OpenMP with ',i3,' threads')", nthreads
print *, "Input n ... "
read *, n
! allocate storage for boundary points too:
allocate(x(0:n+1), u(0:n+1), uold(0:n+1), f(0:n+1))
open(unit=20, file="heatsoln.txt", status="unknown")
! grid spacing:
dx = 1.d0 / (n+1.d0)
! boundary conditions:
alpha = 20.d0
beta = 60.d0
!$omp parallel do
do i=0,n+1
! grid points:
x(i) = i*dx
! source term:
f(i) = 100.*exp(x(i))
! initial guess:
u(i) = alpha + x(i)*(beta-alpha)
enddo
! tolerance and max number of iterations:
tol = 0.1 * dx**2
print *, "Convergence tolerance: tol = ",tol
maxiter = 100000
print *, "Maximum number of iterations: maxiter = ",maxiter
! Jacobi iteratation:
uold = u ! starting values before updating
do iter=1,maxiter
dumax = 0.d0
!$omp parallel do reduction(max : dumax)
do i=1,n
u(i) = 0.5d0*(uold(i-1) + uold(i+1) + dx**2*f(i))
dumax = max(dumax, abs(u(i)-uold(i)))
enddo
if (mod(iter,10000)==0) then
print *, iter, dumax
endif
! check for convergence:
if (dumax .lt. tol) exit
!$omp parallel do
do i=1,n
uold(i) = u(i) ! for next iteration
enddo
enddo
print *, "Total number of iterations: ",iter
write(20,*) " x u"
do i=0,n+1
write(20,222), x(i), u(i)
enddo
222 format(2e20.10)
print *, "Solution is in heatsoln.txt"
close(20)
end program jacobi1