4 ************************************************************
5 * program to solve a finite difference
6 * discretization of Helmholtz equation :
7 * (d2/dx2)u + (d2/dy2)u - alpha u = f
8 * using Jacobi iterative method.
10 * Modified: Sanjiv Shah, Kuck and Associates, Inc. (KAI), 1998
11 * Author: Joseph Robicheaux, Kuck and Associates, Inc. (KAI), 1998
13 * Directives are used in this code to achieve paralleism.
14 * All do loops are parallized with default 'static' scheduling.
16 * Input : n - grid dimension in x direction
17 * m - grid dimension in y direction
18 * alpha - Helmholtz constant (always greater than 0.0)
19 * tol - error tolerance for iterative solver
20 * relax - Successice over relaxation parameter
21 * mits - Maximum iterations for iterative solver
24 * : u(n,m) - Dependent variable (solutions)
25 * : f(n,m) - Right hand side function
26 *************************************************************
29 integer n,m,mits,mtemp
31 double precision tol,relax,alpha
33 common /idat/ n,m,mits,mtemp
34 common /fdat/tol,alpha,relax
38 write(*,*) "Input n,m - grid dimension in x,y direction "
43 write(*,*) "Input alpha - Helmholts constant "
47 write(*,*) "Input relax - Successive over-relaxation parameter"
51 write(*,*) "Input tol - error tolerance for iterative solver"
55 write(*,*) "Input mits - Maximum iterations for solver"
60 call omp_set_num_threads (2)
63 * Calls a driver routine
71 *************************************************************
72 * Subroutine driver ()
73 * This is where the arrays are allocated and initialzed.
75 * Working varaibles/arrays
76 * dx - grid spacing in x direction
77 * dy - grid spacing in y direction
78 *************************************************************
81 integer n,m,mits,mtemp
82 double precision tol,relax,alpha
84 common /idat/ n,m,mits,mtemp
85 common /fdat/tol,alpha,relax
87 double precision u(n,m),f(n,m),dx,dy
91 call initialize (n,m,alpha,dx,dy,u,f)
93 * Solve Helmholtz equation
95 call jacobi (n,m,dx,dy,alpha,relax,u,f,tol,mits)
97 * Check error between exact solution
99 call error_check (n,m,alpha,dx,dy,u,f)
104 subroutine initialize (n,m,alpha,dx,dy,u,f)
105 ******************************************************
107 * Assumes exact solution is u(x,y) = (1-x^2)*(1-y^2)
109 ******************************************************
113 double precision u(n,m),f(n,m),dx,dy,alpha
117 parameter (PI=3.1415926)
122 * Initilize initial condition and RHS
124 !$omp parallel do private(xx,yy)
127 xx = -1.0 + dx * dble(i-1) ! -1 < x < 1
128 yy = -1.0 + dy * dble(j-1) ! -1 < y < 1
130 f(i,j) = -alpha *(1.0-xx*xx)*(1.0-yy*yy)
131 & - 2.0*(1.0-xx*xx)-2.0*(1.0-yy*yy)
134 !$omp end parallel do
139 subroutine jacobi (n,m,dx,dy,alpha,omega,u,f,tol,maxit)
140 ******************************************************************
141 * Subroutine HelmholtzJ
142 * Solves poisson equation on rectangular grid assuming :
143 * (1) Uniform discretization in each direction, and
144 * (2) Dirichlect boundary conditions
146 * Jacobi method is used in this routine
148 * Input : n,m Number of grid points in the X/Y directions
149 * dx,dy Grid spacing in the X/Y directions
150 * alpha Helmholtz eqn. coefficient
151 * omega Relaxation factor
152 * f(n,m) Right hand side function
153 * u(n,m) Dependent variable/Solution
154 * tol Tolerance for iterative solver
155 * maxit Maximum number of iterations
157 * Output : u(n,m) - Solution
158 *****************************************************************
161 double precision dx,dy,f(n,m),u(n,m),alpha, tol,omega
165 integer i,j,k,k_local
166 double precision error,resid,rsum,ax,ay,b
167 double precision error_local, uold(n,m)
169 real ta,tb,tc,td,te,ta1,ta2,tb1,tb2,tc1,tc2,td1,td2
174 * Initialize coefficients
175 ax = 1.0/(dx*dx) ! X-direction coef
176 ay = 1.0/(dy*dy) ! Y-direction coef
177 b = -2.0/(dx*dx)-2.0/(dy*dy) - alpha ! Central coeff
182 do while (k.le.maxit .and. error.gt. tol)
186 * Copy new solution into old
196 * Compute stencil, residual, & update
198 !$omp do private(resid) reduction(+:error)
202 resid = (ax*(uold(i-1,j) + uold(i+1,j))
203 & + ay*(uold(i,j-1) + uold(i,j+1))
204 & + b * uold(i,j) - f(i,j))/b
206 u(i,j) = uold(i,j) - omega * resid
207 * Accumulate residual error
208 error = error + resid*resid
219 error = sqrt(error)/dble(n*m)
221 enddo ! End iteration loop
223 print *, 'Total Number of Iterations ', k
224 print *, 'Residual ', error
229 subroutine error_check (n,m,alpha,dx,dy,u,f)
231 ************************************************************
232 * Checks error between numerical and exact solution
234 ************************************************************
237 double precision u(n,m),f(n,m),dx,dy,alpha
240 double precision xx,yy,temp,error
246 !$omp parallel do private(xx,yy,temp) reduction(+:error)
249 xx = -1.0d0 + dx * dble(i-1)
250 yy = -1.0d0 + dy * dble(j-1)
251 temp = u(i,j) - (1.0-xx*xx)*(1.0-yy*yy)
252 error = error + temp*temp
256 error = sqrt(error)/dble(n*m)
258 print *, 'Solution Error : ',error