Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Poisson Solver: issue with corner values

ivanp
New Contributor I
358 Views

I'm attempting to use the Fast Poisson Solver routines from the oneMKL library. The driver, in Fortran, is provided in the attachment.

The solver works, however the corner values at locations (1,1) and (1,0) seem to be wrong, when compared with the output of a different solver. 

poisson_solution.png


What could be the problem?

Labels (1)
0 Kudos
1 Solution
ivanp
New Contributor I
268 Views

I realized it was an error on my side.  The boundary arrays are of length nx+1 and ny+1, but I was wrongly only setting the first nx, ny values. Hence, there was by default a 0.0 in the last value.

 

The boundary routine is supposed to read

   subroutine fbd(nx,ny,x,y,bd_ax,bd_bx,bd_ay,bd_by)
      integer, intent(in) :: nx, ny
      real(dp), intent(in) :: x(nx+1), y(ny+1)
      real(dp), intent(out) :: bd_ax(ny+1), bd_bx(ny+1)
      real(dp), intent(out) :: bd_ay(nx+1), bd_by(nx+1)

      ! boundaries are given by
      ! u(x,y) = cos( pi x ) - sin( 2 pi y )

      bd_ax = cos(pi*0.0_dp) - sin(2*pi*y) ! x = 0
      bd_bx = cos(pi*1.0_dp) - sin(2*pi*y) ! x = 1

      bd_ay = cos(pi*x) - sin(2*pi*0.0_dp) ! y = 0
      bd_by = cos(pi*x) - sin(2*pi*1.0_dp) ! y = 1
   end subroutine


Now the solution matches the multigrid solver in the link in my original post.


Poisson solution (correct)Poisson solution (correct)

 

View solution in original post

0 Kudos
3 Replies
Gennady_F_Intel
Moderator
271 Views

Ivan, I am not quite sure which results you are expecting to see here.

running with the current version of oneMKL, I see 

1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 

.....
1.00000000000000 1.00000000000000 0.000000000000000E+000 

 

Which results do you see with another solvers? 
 

 

0 Kudos
ivanp
New Contributor I
269 Views

I realized it was an error on my side.  The boundary arrays are of length nx+1 and ny+1, but I was wrongly only setting the first nx, ny values. Hence, there was by default a 0.0 in the last value.

 

The boundary routine is supposed to read

   subroutine fbd(nx,ny,x,y,bd_ax,bd_bx,bd_ay,bd_by)
      integer, intent(in) :: nx, ny
      real(dp), intent(in) :: x(nx+1), y(ny+1)
      real(dp), intent(out) :: bd_ax(ny+1), bd_bx(ny+1)
      real(dp), intent(out) :: bd_ay(nx+1), bd_by(nx+1)

      ! boundaries are given by
      ! u(x,y) = cos( pi x ) - sin( 2 pi y )

      bd_ax = cos(pi*0.0_dp) - sin(2*pi*y) ! x = 0
      bd_bx = cos(pi*1.0_dp) - sin(2*pi*y) ! x = 1

      bd_ay = cos(pi*x) - sin(2*pi*0.0_dp) ! y = 0
      bd_by = cos(pi*x) - sin(2*pi*1.0_dp) ! y = 1
   end subroutine


Now the solution matches the multigrid solver in the link in my original post.


Poisson solution (correct)Poisson solution (correct)

 

0 Kudos
Gennady_F_Intel
Moderator
253 Views

ok, then the thread is closing.


0 Kudos
Reply