- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

http://software.intel.com/sites/products/documentation/hpc/compilerpro/en-us/cpp/win/mkl/refman/appe...

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

Link Copied

5 Replies

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

it's seem that you need to set

do ix=1,nx+1

bd_ay

bd_by

end do

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

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Alexander.

Thank you, that did the trick indeed

Myron

Thank you, that did the trick indeed

Myron

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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).

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page