| * { dg-do run } |
| |
| program main |
| ************************************************************ |
| * program to solve a finite difference |
| * discretization of Helmholtz equation : |
| * (d2/dx2)u + (d2/dy2)u - alpha u = f |
| * using Jacobi iterative method. |
| * |
| * Modified: Sanjiv Shah, Kuck and Associates, Inc. (KAI), 1998 |
| * Author: Joseph Robicheaux, Kuck and Associates, Inc. (KAI), 1998 |
| * |
| * Directives are used in this code to achieve paralleism. |
| * All do loops are parallelized with default 'static' scheduling. |
| * |
| * Input : n - grid dimension in x direction |
| * m - grid dimension in y direction |
| * alpha - Helmholtz constant (always greater than 0.0) |
| * tol - error tolerance for iterative solver |
| * relax - Successice over relaxation parameter |
| * mits - Maximum iterations for iterative solver |
| * |
| * On output |
| * : u(n,m) - Dependent variable (solutions) |
| * : f(n,m) - Right hand side function |
| ************************************************************* |
| implicit none |
| |
| integer n,m,mits,mtemp |
| include "omp_lib.h" |
| double precision tol,relax,alpha |
| |
| common /idat/ n,m,mits,mtemp |
| common /fdat/tol,alpha,relax |
| * |
| * Read info |
| * |
| write(*,*) "Input n,m - grid dimension in x,y direction " |
| n = 64 |
| m = 64 |
| * read(5,*) n,m |
| write(*,*) n, m |
| write(*,*) "Input alpha - Helmholts constant " |
| alpha = 0.5 |
| * read(5,*) alpha |
| write(*,*) alpha |
| write(*,*) "Input relax - Successive over-relaxation parameter" |
| relax = 0.9 |
| * read(5,*) relax |
| write(*,*) relax |
| write(*,*) "Input tol - error tolerance for iterative solver" |
| tol = 1.0E-12 |
| * read(5,*) tol |
| write(*,*) tol |
| write(*,*) "Input mits - Maximum iterations for solver" |
| mits = 100 |
| * read(5,*) mits |
| write(*,*) mits |
| |
| call omp_set_num_threads (2) |
| |
| * |
| * Calls a driver routine |
| * |
| call driver () |
| |
| stop |
| end |
| |
| subroutine driver ( ) |
| ************************************************************* |
| * Subroutine driver () |
| * This is where the arrays are allocated and initialized. |
| * |
| * Working varaibles/arrays |
| * dx - grid spacing in x direction |
| * dy - grid spacing in y direction |
| ************************************************************* |
| implicit none |
| |
| integer n,m,mits,mtemp |
| double precision tol,relax,alpha |
| |
| common /idat/ n,m,mits,mtemp |
| common /fdat/tol,alpha,relax |
| |
| double precision u(n,m),f(n,m),dx,dy |
| |
| * Initialize data |
| |
| call initialize (n,m,alpha,dx,dy,u,f) |
| |
| * Solve Helmholtz equation |
| |
| call jacobi (n,m,dx,dy,alpha,relax,u,f,tol,mits) |
| |
| * Check error between exact solution |
| |
| call error_check (n,m,alpha,dx,dy,u,f) |
| |
| return |
| end |
| |
| subroutine initialize (n,m,alpha,dx,dy,u,f) |
| ****************************************************** |
| * Initializes data |
| * Assumes exact solution is u(x,y) = (1-x^2)*(1-y^2) |
| * |
| ****************************************************** |
| implicit none |
| |
| integer n,m |
| double precision u(n,m),f(n,m),dx,dy,alpha |
| |
| integer i,j, xx,yy |
| double precision PI |
| parameter (PI=3.1415926) |
| |
| dx = 2.0 / (n-1) |
| dy = 2.0 / (m-1) |
| |
| * Initialize initial condition and RHS |
| |
| !$omp parallel do private(xx,yy) |
| do j = 1,m |
| do i = 1,n |
| xx = -1.0 + dx * dble(i-1) ! -1 < x < 1 |
| yy = -1.0 + dy * dble(j-1) ! -1 < y < 1 |
| u(i,j) = 0.0 |
| f(i,j) = -alpha *(1.0-xx*xx)*(1.0-yy*yy) |
| & - 2.0*(1.0-xx*xx)-2.0*(1.0-yy*yy) |
| enddo |
| enddo |
| !$omp end parallel do |
| |
| return |
| end |
| |
| subroutine jacobi (n,m,dx,dy,alpha,omega,u,f,tol,maxit) |
| ****************************************************************** |
| * Subroutine HelmholtzJ |
| * Solves poisson equation on rectangular grid assuming : |
| * (1) Uniform discretization in each direction, and |
| * (2) Dirichlect boundary conditions |
| * |
| * Jacobi method is used in this routine |
| * |
| * Input : n,m Number of grid points in the X/Y directions |
| * dx,dy Grid spacing in the X/Y directions |
| * alpha Helmholtz eqn. coefficient |
| * omega Relaxation factor |
| * f(n,m) Right hand side function |
| * u(n,m) Dependent variable/Solution |
| * tol Tolerance for iterative solver |
| * maxit Maximum number of iterations |
| * |
| * Output : u(n,m) - Solution |
| ***************************************************************** |
| implicit none |
| integer n,m,maxit |
| double precision dx,dy,f(n,m),u(n,m),alpha, tol,omega |
| * |
| * Local variables |
| * |
| integer i,j,k,k_local |
| double precision error,resid,rsum,ax,ay,b |
| double precision error_local, uold(n,m) |
| |
| real ta,tb,tc,td,te,ta1,ta2,tb1,tb2,tc1,tc2,td1,td2 |
| real te1,te2 |
| real second |
| external second |
| * |
| * Initialize coefficients |
| ax = 1.0/(dx*dx) ! X-direction coef |
| ay = 1.0/(dy*dy) ! Y-direction coef |
| b = -2.0/(dx*dx)-2.0/(dy*dy) - alpha ! Central coeff |
| |
| error = 10.0 * tol |
| k = 1 |
| |
| do while (k.le.maxit .and. error.gt. tol) |
| |
| error = 0.0 |
| |
| * Copy new solution into old |
| !$omp parallel |
| |
| !$omp do |
| do j=1,m |
| do i=1,n |
| uold(i,j) = u(i,j) |
| enddo |
| enddo |
| |
| * Compute stencil, residual, & update |
| |
| !$omp do private(resid) reduction(+:error) |
| do j = 2,m-1 |
| do i = 2,n-1 |
| * Evaluate residual |
| resid = (ax*(uold(i-1,j) + uold(i+1,j)) |
| & + ay*(uold(i,j-1) + uold(i,j+1)) |
| & + b * uold(i,j) - f(i,j))/b |
| * Update solution |
| u(i,j) = uold(i,j) - omega * resid |
| * Accumulate residual error |
| error = error + resid*resid |
| end do |
| enddo |
| !$omp enddo nowait |
| |
| !$omp end parallel |
| |
| * Error check |
| |
| k = k + 1 |
| |
| error = sqrt(error)/dble(n*m) |
| * |
| enddo ! End iteration loop |
| * |
| print *, 'Total Number of Iterations ', k |
| print *, 'Residual ', error |
| |
| return |
| end |
| |
| subroutine error_check (n,m,alpha,dx,dy,u,f) |
| implicit none |
| ************************************************************ |
| * Checks error between numerical and exact solution |
| * |
| ************************************************************ |
| |
| integer n,m |
| double precision u(n,m),f(n,m),dx,dy,alpha |
| |
| integer i,j |
| double precision xx,yy,temp,error |
| |
| dx = 2.0 / (n-1) |
| dy = 2.0 / (m-1) |
| error = 0.0 |
| |
| !$omp parallel do private(xx,yy,temp) reduction(+:error) |
| do j = 1,m |
| do i = 1,n |
| xx = -1.0d0 + dx * dble(i-1) |
| yy = -1.0d0 + dy * dble(j-1) |
| temp = u(i,j) - (1.0-xx*xx)*(1.0-yy*yy) |
| error = error + temp*temp |
| enddo |
| enddo |
| |
| error = sqrt(error)/dble(n*m) |
| |
| print *, 'Solution Error : ',error |
| |
| return |
| end |
