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

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??

Link Copied

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

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

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

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

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

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?

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

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)

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

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?

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

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.

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

Thanks for answering my query.

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

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/t... (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

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