Problem with Poisson solver Neumann boundary conditions

[barnaclejunior]

I'm having problems solving a Poisson equation using MKL's s_Helmholtz_2D, Win32 binaries, 10.0.4.023...

When I use Dirichlet boundary conditions with zeros, the function executes and the results are ok. When I use Neumann boundary conditions, I invariably get this warning printed to the console when I call s_Helmholtz_2D:

MKL POISSON LIBRARY WARNING:
The problem is degenerate up to rounding errors! The approximate solution
that provides the minimal Eucledean (sic) norm of the residual will be computed...

The function call also fills the entire input/output array with -1.#IND000, except for the left boundary, which gets filled with 1.#QNAN00. I've attached my solver code below. I pass in the Laplacian (positive), and the x and y gradients of the solution. Because MKL wants the negative lap, and I have the positive lap, I negate the Bx and By boundary conditions rather than the Ax and Ay conditions, so I should get the negative solution, which I negate in the loop at the end.

Any ideas are much appreciated.

.sean

void PoissonSolve(const FloatMatrix& lap, const FloatMatrix& Gx, const FloatMatrix& Gy,
FloatMatrixPtr* ppSolution) {

IntSize2 size = lap.size;
float ax = 0, bx = size.width, ay = 0, by = size.height;
int nx = size.width - 1, ny = size.height - 1;
float q = 0;
int ipar[128];
std::vector spar(5 * nx / 2 + 7);
int stat;

s_init_Helmholtz_2D(&ax, &bx, &ay, &by, &nx, &ny, "NNNN", &q, ipar, &spar[0], &stat);

std::vector boundaryAx(ny + 1), boundaryBx(ny + 1),
boundaryAy(nx + 1), boundaryBy(nx + 1);
for(int row(0); row < size.height; ++row) {
boundaryAx[row] = Gx.Row(row)[0];
boundaryBx[row] = -Gx.Row(row)[size.width - 1];
}
for(int i(0); i < size.width; ++i) {
boundaryAy = Gy.Row(0);
boundaryBy = -Gy.Row(size.height - 1);
}

assert(0 == stat);

std::vector > data(size.Area());
for(int row(0); row < size.height; ++row)
memcpy(&data[row * size.width], lap.Row(row), sizeof(float) * size.width);

DFTI_DESCRIPTOR_HANDLE xHandle;
s_commit_Helmholtz_2D(&data[0], &boundaryAx[0], &boundaryAy[0], &boundaryBx[0],
&boundaryBy[0], &xHandle, ipar, &spar[0], &stat);
s_Helmholtz_2D(&data[0], &boundaryAx[0], &boundaryAy[0], &boundaryBx[0],
&boundaryBy[0], &xHandle, ipar, &spar[0], &stat);
free_Helmholtz_2D(&xHandle, ipar, &stat);

// copy the results into the solution matrix an d negate for the positive solution
FloatMatrixPtr solution(new FloatMatrix(size));
for(int row(0); row < size.height; ++row) {
memcpy(solution->Row(row), &data[row * size.width], sizeof(float) * size.width);
for(int col(0); col < size.width; ++col)
solution->Row(row)[col] *= -1;
}

*ppSolution = solution;
}

[barnaclejunior]

One other thing... It seems like the Dirichlet boundary conditions require arrays bigger than ny + 1, nx + 1 to return good results. I size them to ny + 1000, nx + 1000 and the results are fine... Is the solver actually reading beyond the specified sizes of the boundary arrays?

[barnaclejunior]

In fact, if I feed the solver Neumann boundary conditions with large boundary arrays (like ny + 10000), the solver executes without producing NAN's or anything like that, but just fills the output buffer with a smooth gradient... The solution is still rubbish, and I still get the warning with the mispelling of "Euclidean," but at least it executes.

[Sergey_G_Intel]

Hi Sean!

Thank you for using Intel MKL!

It would be great if you can provide us with the test cases to understand your issue. It is hard to say something about the root cause of your problems looking at the code of the routine. It might be that the size of the problem is computed incorrectly, say, 1 less than it should be. It also could be that some troubles come from the loops where arrays like boundaryAx[row] are filled. As you write ++row, you are filling the arrays starting from position 1, while you provide MKL routines with arrays starting from position 0.

As for the warning, you can remove it using proper setting of ipar array: ipar[2]=0.

Thank you for pointing out the misprinting in the message!

-Sergey G

[barnaclejunior]

The ++row doesn't matter; it's the same as row++. The boundary condition arrays are being filled correctly. Can you define the meaning of the ax, bx, ay, by terms passed to s_init_Helmholtz_2D? I'm setting them to 0, size.width, 0, size.height, respectively, and nx = size.width - 1, ny = size.height - 1... Is this correct for processing a size.width x size.height image? Also, am I specifying the correct signs on the boundary derivatives? The specs say that the ax and ay boundaries take the partials of the field wrt -x and -y, respectively... So that's just the .x and .y components of the gradient times -1?

thanks

[Sergey_G_Intel]

Hi Sean!

It seems youve encountered a problem in Helmholtz solver. Could you please check if your program will work fine with ipar[256] setting instead of ipar[128]?

The meaning of ax, ay, bx, by is to describe rectangular domain x=ax (left boundary), x=bx (right boundary), y=ay (lower boundary), y=by (upper boundary). With these settings the computations will be done for the domain ax

It looks to me that your settings of the boundary derivatives are correct.

[barnaclejunior]

Quoting - SERGEY G (Intel)

I gave ipar 4k elements and the results are the same.. Basically, if the boundary arrays are ny + 1, nx + 1 in size, respectively, i get INF and QNAN for results. If I make the boundary arrays very large (just allocate a lot of heap space and still only fill out the first ny + 1, nx + 1 elements), the solver returns stable results, but the results are competely wrong, just giving me a smooth gradient between corners.

thanks,

Hi Sean!

It seems youve encountered a problem in Helmholtz solver. Could you please check if your program will work fine with ipar[256] setting instead of ipar[128]?

The meaning of ax, ay, bx, by is to describe rectangular domain x=ax (left boundary), x=bx (right boundary), y=ay (lower boundary), y=by (upper boundary). With these settings the computations will be done for the domain ax

It looks to me that your settings of the boundary derivatives are correct.

[Sergey_G_Intel]

Hi Sean,

Sorry for the delay with the reply. It seems the system deleted my alerts for this forum thread. Could you provide us with the test case to uderstandthe root cause of theproblem?It seems that on our side everything works fine.

Kind regards,

Sergey