mkl poisson solver, boundary condtion problem

[cavaiola__mattia]

Hi, I'm trying to solve the Poisson equation with the Fast Poisson solver.  I'm started by imposing a fully periodic 3D domain, then my BCtype='PPPPPP', but I get strange results, with a warning message: 

Intel MKL POISSON SOLVER WARNING: 
The problem is degenerate due to rounding errors. The approximate solution
that provides the minimal Euclidean norm of the solution will be computed.

I looked at the values of ipar, by referring to this page https://software.intel.com/content/www/us/en/develop/documentation/mkl-developer-reference-fortran/top/partial-differential-equations-support/fast-poisson-solver-routines/common-parameters-for-the-poisson-solver/ipar.html

From this page, one can see that the integer values of the array ipar from 4 to 9,  takes this value according to BCtype. In my case, they could be all equal to 2.

Then I print the values of ipar and they have different values form 2. I get 0 0 0 2 2 2.

Here there is my code written in fortran

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine solvermkl(phi)
implicit none
real(rp) :: ax,bx,ay,by,az,bz
integer :: i,j,k,l,stat
integer, parameter :: nx=(imax-1), ny=(jmax-1),nz=(kmax-1)
real(rp), parameter :: q = 0. ! constant for Helmotz equation --> 0 for Poisson
integer, dimension(128) :: ipar
real(rp), dimension(13.*(nx+ny)/2.+9.) :: dpar
type(DFTI_DESCRIPTOR), pointer :: xhandle, yhandle
character(6) BCtype
real(rp), dimension(imax,jmax,kmax), intent(out) :: phi
real(rp) :: bd_ax(jmax,kmax), bd_bx(jmax,kmax), bd_ay(imax,kmax), bd_by(imax,kmax), bd_az(imax,jmax), bd_bz(imax,jmax)

ax= 0.
bx= lx
ay= 0.
by= ly
az= 0.
bz= lz

do l=1,128
   ipar(l)=0
enddo

BCtype = 'PPPPPP'

call d_init_Helmholtz_3D(ax, bx, ay, by, az, bz, nx, ny, nz, BCtype, q, ipar, dpar, stat)

do l=4,9
 print *, 'ipar=', ipar(l)
enddo

call d_commit_Helmholtz_3D(RHS, bd_ax, bd_bx, bd_ay, bd_by, bd_az, bd_bz, xhandle,yhandle, ipar, dpar, stat)

call d_Helmholtz_3D(RHS, bd_ax, bd_bx, bd_ay, bd_by, bd_az, bd_bz, xhandle, yhandle, ipar, dpar, stat)

call free_Helmholtz_3D(xhandle, yhandle, ipar, stat)

phi(:,:,:) = RHS(:,:,:)

end subroutine solvermkl

end module mod_solver_mkl

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

There is someone can help me??

[Kirill_V_Intel]

Hi,

Without running the code, I notice that you have a 3D Poisson equation (q = 0) with periodic boundary conditions everywhere. For this case, the solution to the problem is not unique. For example, adding a constant to a solution is also a solution. Hence the problem is degenerate and that's what the warning is saying.

So I think you should re-consider the statement of the problem you're trying to solve. Typically, one can avoid degeneracy by imposing additional condition or fixing a solution somewhere. 

We will check why the values of ipar are changed the way you see. Immediate guess is that's because we find the solution with the minimal norm (essentially, that's an additional condition MKL uses to make the solution unique).

Thanks,
Kirill

[cavaiola__mattia]

Thank you for your answer... But I already solved Poisson with periodic boundary conditions everywhere, with other direct solvers, without encountering this kind of problem

As regards the ipar, by running the tutorial case, such as d_Poisson_3D_f.f90, I noticed the same shift of index. Seems that the information about boundary conditions are stored in ipar from index 7 to 12. I think there is an error in https://software.intel.com/content/www/us/en/develop/documentation/mkl-developer-reference-c/top/partial-differential-equations-support/fast-poisson-solver-routines/common-parameters-for-the-poisson-solver/ipar.html

[mecej4]

In C, array indices start at 0. In Fortran, by default, array indices start at 1. If you wish to use 0-base indices in Fortran, you have to dimension ipar(0:127), rather than ipar(128). You could avoid some of this confusion by using the MKL documentation for Fortran, instead of MKL-C.

Before calling ?_commit_Helmholtz_3D, you have to set values into the array RHS, which you have left undefined in the code above.

The document for the 2020 version of MKL, at https://software.intel.com/content/www/us/en/develop/documentation/mkl-developer-reference-fortran/top/partial-differential-equations-support/fast-poisson-solver-routines/routines-for-the-cartesian-solver/init-helmholtz-2d-init-helmholtz-3d.html , supposedly for Fortran users, unfortunately displays only C syntax! You may need to look at the include file mkl_poisson.f90 to clarify the interface. Or is there a new expectation that Fortran users of MKL should start learning C?

[cavaiola__mattia]

RHS is not undefined... It is defined outside of this routine...

I know the difference in indexing between C and Fortran.

The problem with ipar is that the shift is about 2 positions when I declare ipar(0:127), and 3 positions when ipar(1:128)

[mecej4]

I am trying to troubleshoot this issue, and it would be helpful if you would state the mathematical problem being solved. Poisson equation in 3D box, with periodic boundary conditions, source term (RHS) specified -- this part I see. Can I take it that RHS = uniform, say 1.0, or does it have to have some specified spatial distribution?

[cavaiola__mattia]

This is a part of code to solve the incompressible Navier-Stokes equations. I need to solve the Pressure Poisson equation, where the right-hand side is:

rhs= 1/dt*(div(u*)). Then in general it is not uniform.

[mecej4]

Thanks for answering my query.

[Kirill_V_Intel]

Hello again,

(somehow my previous reply has not been published so I repeat it)

I've looked at the issue. Most likely it is a documentation bug, description of the ipar array is incorrect. We're going to  check it all carefully and fix it.

At the doc page https://software.intel.com/content/www/us/en/develop/documentation/mkl-developer-reference-fortran/top/partial-differential-equations-support/fast-poisson-solver-routines/common-parameters-for-the-poisson-solver/ipar.html (notice, it is for C) it is said that ipar[15]-[16] are unused and ipar[4]-ipar[9] descibed the boundary conditions. However, it should say instead that ipar[4]-ipar[5] are unused, and shift the definitions of ipar[4]-ipar[14] two items below, so the BC are described in ipar[6]-ipar[11]. In Fortran, hence, you should see BC in ipar(7)-ipar(12) and I believe in your case you will see all 2's there.

As for the warning, as I said earlier, it comes from the fact that the poblem is degenerate and I believe that this behavior is correct.

Hope this helps.

Best,
Kirill