|
0
|
1 * { dg-do run }
|
|
|
2
|
|
|
3 program main
|
|
|
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.
|
|
|
9 *
|
|
|
10 * Modified: Sanjiv Shah, Kuck and Associates, Inc. (KAI), 1998
|
|
|
11 * Author: Joseph Robicheaux, Kuck and Associates, Inc. (KAI), 1998
|
|
|
12 *
|
|
|
13 * Directives are used in this code to achieve paralleism.
|
|
|
14 * All do loops are parallized with default 'static' scheduling.
|
|
|
15 *
|
|
|
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
|
|
|
22 *
|
|
|
23 * On output
|
|
|
24 * : u(n,m) - Dependent variable (solutions)
|
|
|
25 * : f(n,m) - Right hand side function
|
|
|
26 *************************************************************
|
|
|
27 implicit none
|
|
|
28
|
|
|
29 integer n,m,mits,mtemp
|
|
|
30 include "omp_lib.h"
|
|
|
31 double precision tol,relax,alpha
|
|
|
32
|
|
|
33 common /idat/ n,m,mits,mtemp
|
|
|
34 common /fdat/tol,alpha,relax
|
|
|
35 *
|
|
|
36 * Read info
|
|
|
37 *
|
|
|
38 write(*,*) "Input n,m - grid dimension in x,y direction "
|
|
|
39 n = 64
|
|
|
40 m = 64
|
|
|
41 * read(5,*) n,m
|
|
|
42 write(*,*) n, m
|
|
|
43 write(*,*) "Input alpha - Helmholts constant "
|
|
|
44 alpha = 0.5
|
|
|
45 * read(5,*) alpha
|
|
|
46 write(*,*) alpha
|
|
|
47 write(*,*) "Input relax - Successive over-relaxation parameter"
|
|
|
48 relax = 0.9
|
|
|
49 * read(5,*) relax
|
|
|
50 write(*,*) relax
|
|
|
51 write(*,*) "Input tol - error tolerance for iterative solver"
|
|
|
52 tol = 1.0E-12
|
|
|
53 * read(5,*) tol
|
|
|
54 write(*,*) tol
|
|
|
55 write(*,*) "Input mits - Maximum iterations for solver"
|
|
|
56 mits = 100
|
|
|
57 * read(5,*) mits
|
|
|
58 write(*,*) mits
|
|
|
59
|
|
|
60 call omp_set_num_threads (2)
|
|
|
61
|
|
|
62 *
|
|
|
63 * Calls a driver routine
|
|
|
64 *
|
|
|
65 call driver ()
|
|
|
66
|
|
|
67 stop
|
|
|
68 end
|
|
|
69
|
|
|
70 subroutine driver ( )
|
|
|
71 *************************************************************
|
|
|
72 * Subroutine driver ()
|
|
|
73 * This is where the arrays are allocated and initialzed.
|
|
|
74 *
|
|
|
75 * Working varaibles/arrays
|
|
|
76 * dx - grid spacing in x direction
|
|
|
77 * dy - grid spacing in y direction
|
|
|
78 *************************************************************
|
|
|
79 implicit none
|
|
|
80
|
|
|
81 integer n,m,mits,mtemp
|
|
|
82 double precision tol,relax,alpha
|
|
|
83
|
|
|
84 common /idat/ n,m,mits,mtemp
|
|
|
85 common /fdat/tol,alpha,relax
|
|
|
86
|
|
|
87 double precision u(n,m),f(n,m),dx,dy
|
|
|
88
|
|
|
89 * Initialize data
|
|
|
90
|
|
|
91 call initialize (n,m,alpha,dx,dy,u,f)
|
|
|
92
|
|
|
93 * Solve Helmholtz equation
|
|
|
94
|
|
|
95 call jacobi (n,m,dx,dy,alpha,relax,u,f,tol,mits)
|
|
|
96
|
|
|
97 * Check error between exact solution
|
|
|
98
|
|
|
99 call error_check (n,m,alpha,dx,dy,u,f)
|
|
|
100
|
|
|
101 return
|
|
|
102 end
|
|
|
103
|
|
|
104 subroutine initialize (n,m,alpha,dx,dy,u,f)
|
|
|
105 ******************************************************
|
|
|
106 * Initializes data
|
|
|
107 * Assumes exact solution is u(x,y) = (1-x^2)*(1-y^2)
|
|
|
108 *
|
|
|
109 ******************************************************
|
|
|
110 implicit none
|
|
|
111
|
|
|
112 integer n,m
|
|
|
113 double precision u(n,m),f(n,m),dx,dy,alpha
|
|
|
114
|
|
|
115 integer i,j, xx,yy
|
|
|
116 double precision PI
|
|
|
117 parameter (PI=3.1415926)
|
|
|
118
|
|
|
119 dx = 2.0 / (n-1)
|
|
|
120 dy = 2.0 / (m-1)
|
|
|
121
|
|
|
122 * Initilize initial condition and RHS
|
|
|
123
|
|
|
124 !$omp parallel do private(xx,yy)
|
|
|
125 do j = 1,m
|
|
|
126 do i = 1,n
|
|
|
127 xx = -1.0 + dx * dble(i-1) ! -1 < x < 1
|
|
|
128 yy = -1.0 + dy * dble(j-1) ! -1 < y < 1
|
|
|
129 u(i,j) = 0.0
|
|
|
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)
|
|
|
132 enddo
|
|
|
133 enddo
|
|
|
134 !$omp end parallel do
|
|
|
135
|
|
|
136 return
|
|
|
137 end
|
|
|
138
|
|
|
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
|
|
|
145 *
|
|
|
146 * Jacobi method is used in this routine
|
|
|
147 *
|
|
|
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
|
|
|
156 *
|
|
|
157 * Output : u(n,m) - Solution
|
|
|
158 *****************************************************************
|
|
|
159 implicit none
|
|
|
160 integer n,m,maxit
|
|
|
161 double precision dx,dy,f(n,m),u(n,m),alpha, tol,omega
|
|
|
162 *
|
|
|
163 * Local variables
|
|
|
164 *
|
|
|
165 integer i,j,k,k_local
|
|
|
166 double precision error,resid,rsum,ax,ay,b
|
|
|
167 double precision error_local, uold(n,m)
|
|
|
168
|
|
|
169 real ta,tb,tc,td,te,ta1,ta2,tb1,tb2,tc1,tc2,td1,td2
|
|
|
170 real te1,te2
|
|
|
171 real second
|
|
|
172 external second
|
|
|
173 *
|
|
|
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
|
|
|
178
|
|
|
179 error = 10.0 * tol
|
|
|
180 k = 1
|
|
|
181
|
|
|
182 do while (k.le.maxit .and. error.gt. tol)
|
|
|
183
|
|
|
184 error = 0.0
|
|
|
185
|
|
|
186 * Copy new solution into old
|
|
|
187 !$omp parallel
|
|
|
188
|
|
|
189 !$omp do
|
|
|
190 do j=1,m
|
|
|
191 do i=1,n
|
|
|
192 uold(i,j) = u(i,j)
|
|
|
193 enddo
|
|
|
194 enddo
|
|
|
195
|
|
|
196 * Compute stencil, residual, & update
|
|
|
197
|
|
|
198 !$omp do private(resid) reduction(+:error)
|
|
|
199 do j = 2,m-1
|
|
|
200 do i = 2,n-1
|
|
|
201 * Evaluate residual
|
|
|
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
|
|
|
205 * Update solution
|
|
|
206 u(i,j) = uold(i,j) - omega * resid
|
|
|
207 * Accumulate residual error
|
|
|
208 error = error + resid*resid
|
|
|
209 end do
|
|
|
210 enddo
|
|
|
211 !$omp enddo nowait
|
|
|
212
|
|
|
213 !$omp end parallel
|
|
|
214
|
|
|
215 * Error check
|
|
|
216
|
|
|
217 k = k + 1
|
|
|
218
|
|
|
219 error = sqrt(error)/dble(n*m)
|
|
|
220 *
|
|
|
221 enddo ! End iteration loop
|
|
|
222 *
|
|
|
223 print *, 'Total Number of Iterations ', k
|
|
|
224 print *, 'Residual ', error
|
|
|
225
|
|
|
226 return
|
|
|
227 end
|
|
|
228
|
|
|
229 subroutine error_check (n,m,alpha,dx,dy,u,f)
|
|
|
230 implicit none
|
|
|
231 ************************************************************
|
|
|
232 * Checks error between numerical and exact solution
|
|
|
233 *
|
|
|
234 ************************************************************
|
|
|
235
|
|
|
236 integer n,m
|
|
|
237 double precision u(n,m),f(n,m),dx,dy,alpha
|
|
|
238
|
|
|
239 integer i,j
|
|
|
240 double precision xx,yy,temp,error
|
|
|
241
|
|
|
242 dx = 2.0 / (n-1)
|
|
|
243 dy = 2.0 / (m-1)
|
|
|
244 error = 0.0
|
|
|
245
|
|
|
246 !$omp parallel do private(xx,yy,temp) reduction(+:error)
|
|
|
247 do j = 1,m
|
|
|
248 do i = 1,n
|
|
|
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
|
|
|
253 enddo
|
|
|
254 enddo
|
|
|
255
|
|
|
256 error = sqrt(error)/dble(n*m)
|
|
|
257
|
|
|
258 print *, 'Solution Error : ',error
|
|
|
259
|
|
|
260 return
|
|
|
261 end
|
