Boundary conditions MKL Helholtz solver

[vdamme]

When creating a test program to solve the Poisson equation I obtained the wrong result when using the MKL Helmholtz solver. The basis for the test program was the example given in the link below:
http://software.intel.com/sites/products/documentation/hpc/compilerpro/en-us/cpp/win/mkl/refman/appendices/mkl_appC_PLR.html

I did make the following changes.
1: I saved the vector f as a 1-dimensional array instead of a 2D one as given in the example
2: In the example code, the following lines are used. xi=hx*(ix-1)/lx and yi=hy*(iy-1)/ly
Since ix is an integer, xi will always be rounded of to the nearest integer although it is defined as a double precision. Hence instead I used: xi=hx*real(ix-1)/lx and yi=hy*real(iy-1)/ly
3: In the example program q is set to 1. However to solve the Poisson equation q should be 0.0d0

The core program can be found below
When applying Dirichlet boundary conditions everywhere, I manage to obtain the correct answer, however when applying Von Neumann boundary conditions, even on one side, the error is very large (maximum value for the error is 0.89. The core program can be found below. Does anyone have an idea where I went wrong?


The core program is now
ax=0.0d0
bx=1.0d0
ay=0.0d0
by=1.0d0
nx=gridx !gridx=10
ny=gridz !gridz=10
q=0.0d0
ipar(:)=0
dpar(:)=0.0d0
spar(:)=0.0d0

dx=1.0d0/gridx
dz=1.0d0/gridz
lx=bx-ax
hx=lx/nx
ly=by-ay
hy=ly/ny


do iy=1,ny+1
do ix=1,nx+1
xi=hx*real(ix-1)/lx
yi=hy*real(iy-1)/ly
cx=sin(2.0d0*pi*xi)
cy=sin(2.0d0*pi*yi)
u(ix+(iy-1)*(nx+1))=cx*cy
f(ix+(iy-1)*(nx+1))=((8.0d0*pi*pi))*u(ix+(iy-1)*(nx+1))
u(ix+(iy-1)*(nx+1))=u(ix+(iy-1)*(nx+1))+1.0d0
end do
end do

bd_ax(:)=1.0d0
bd_bx(:)=1.0d0

!do iy=1,nx+1
!bd_ax=-2.0d0*pi*sin(2.0d0*pi*real(iy-1)/ny)
!bd_bx=2.0d0*pi*sin(2.0d0*pi*real(iy-1)/ny)
!end do
do ix=1,nx+1
bd_ay=-2.0d0*pi*sin(2.0d0*pi*real(ix-1)/nx)
bd_by=2.0d0*pi*sin(2.0d0*pi*real(ix-1)/nx)
end do

bd_ax(:)=1.0d0
bd_bx(:)=1.0d0
!bd_ay(:)=1.0d0
!bd_by(:)=1.0d0

BCtype='DDNN'
ipar(:)=0

! Initializing simple data structures of Poisson Library for 2D Helmholtz Solver
call d_init_Helmholtz_2D(ax, bx, ay, by, nx, ny, BCtype, q, ipar, dpar, stat)
if(stat.ne.0)then
write(*,*)'Double precision 2D Helmholtz example FAILED to compute the solution...'
end if

! Computing the approximate solution of 2D Helmholtz problem
! between the Commit step and the subsequent call to the Solver routine!
! Otherwise the results may be wrong.
call d_commit_Helmholtz_2D(f, bd_ax, bd_bx, bd_ay, bd_by, xhandle, ipar, dpar, stat)
if(stat.ne.0)then
write(*,*)'Double precision 2D Helmholtz example FAILED to compute the solution...'
end if

! NOTE: Boundary data stored in the arrays bd_ax, bd_bx, bd_ay, bd_by should not be changed
! between the Commit step and the subsequent call to the Solver routine!
call d_Helmholtz_2D(f, bd_ax, bd_bx, bd_ay, bd_by, xhandle, ipar, dpar, stat)
if(stat.ne.0)then
write(*,*)'Double precision 2D Helmholtz example FAILED to compute the solution...'
end if


ix=3

c1=0.0d0
do iy=1,ny+1
write(*,*) (ix-1)*hx,(iy-1)*hy,f(ix+(iy-1)*(ny+1))-u(ix+(iy-1)*(ny+1))
end do


do i=1,gridx+1
do j=1,gridz+1
write(*,*)i,j, 'ans',f((i)+(j-1)*(gridx+1)),u((i)+(j-1)*(gridx+1))
end do
end do

pause

[Alexander_K_Intel2]

Hi,

it's seem that you need to set

do ix=1,nx+1
bd_ay(i)=-2.0d0*pi*sin(2.0d0*pi*real(ix-1)/nx)
bd_by(i)=2.0d0*pi*sin(2.0d0*pi*real(ix-1)/nx)
end do

instead of

do ix=1,nx+1
bd_ay=-2.0d0*pi*sin(2.0d0*pi*real(ix-1)/nx)
bd_by=2.0d0*pi*sin(2.0d0*pi*real(ix-1)/nx)
end do

Nevertheless I will check it later today.

With best regards,

Alexander Kalinkin

[vdamme]

Alexander.

Thank you, that did the trick indeed

Myron

[vdamme]

Dear Alexander

I do now still have a problem when running the code with all Von Neumann boundary conditions. When I attempt this I get the following MKL Poisson Library warning:
The problem is degenerate up to rounding errors! the approximate solution that provides the mimimal Euclidian norm or the solution will be computed...

An error is expected since an infinite number of solutions is possible when using 4 Von Neumann boundary conditions, however the gradients of the numerical solution should match the gradients found with the analytical solution. This seems not to be the case at the bottom left and bottom right corner of the computational domain. Everywhere else, the solution is OK. Do you have any idea how this larger error in the lower corner points (Error value of 0.18 on a 10x10 grid) is initiated and is there a way for me to prevent it.

Thank you for the help
Myron

[mecej4]

When you obtain a finite-difference solution for a 2-D problem with discontinuities in the boundary conditions, you should go in with reasonable expectations. Typically, discontinuities are smoothed over a small region. As the mesh is made finer, this smoothed region becomes smaller. Just as an experiment cannot be set up with discontinuous boundary values, a finite-difference approximation will be inaccurate in regions of high (or infinite) gradients.

If you want accurate solutions near a discontinuity, you should consider alternative methods of solution.

It is often troublesome when contrived problems are solved numerically whose virtue is that a simple analytical solution is possible. The numerical solution should not be expected to have the same properties as a Fourier series solution (for example, convergence to the mean at a point of discontinuity).

[Alexander_K_Intel2]

Hi Myron,

In case when you use only zero Neumann boundary condition the kernel space is not empty. Indifferential equation it is mean that gradients of all solutions are equal but it is not clear how to calculate gradient of grid solution. You can calculate numeriacal aproximation of gradients in each internal point of solution but what about boundary point? Could you explain how you calculate gradient of solution in corner points?

With best regards,

Alexander Kalinkin