Fast Poisson Solver- Boundary values

[Sav__Joe]

Hello all, 

I am using the Intel MKL Fast Poisson Solver to solve a standard Poisson problem (q=0) with Dirichlet boundary conditions everywhere (BC = "DDDDDD"). This has thus far worked brilliantly. For the purposes of optimizing my solver however I break the calculation up into two steps. The first is for the "compact support" problem, where I define the forcing function "f" everywhere in the domain EXCEPT the boundary, where it is set zero, this is working fine.

The second step is carried out later and I solve the boundary with a given set of BC's, and forcing "f" defined ONLY on the boundary. The BC appears to influence the solution, however it appears as though the solver completely ignores the value of f on the boundary, or at least the solution it produces is exactly the same as the solution where the boundary has zero forcing f=0. 

Perhaps I have misunderstood the implementation of the FFT for such a solver and the second step is not permissible? This however to me seems very counter intuitive, and the solver has work perfectly otherwise.

Thanks in advance for any help.

[mecej4]

The Helmholtz/Poisson/Laplace solver requires a consistent sequence of calls to specify and solve the problem. The array f does double duty: it contains the source term values at the grid points on input, and delivers the solution for temperature on output. You cannot obtain correct results if you modify some of the arrays in unapproved ways midway through the sequence.

The manual page for Helmholtz_2D/Helmholtz_3D says, "Note that an attempt to substitute the original right-hand side vector, which was passed to the ?_commit_Helmholtz_2D/?_commit_Helmholtz_3D routine, at this point results in an incorrect solution. "

[Sav__Joe]

Hi @mecej4  and thanks for the quick reply,

I am of course aware that the f array acts both as the 'f' input and also the solution output, otherwise my other uses of the function would have been unsuccessful/confusing.

The "true" input array which is copied into the 'f' array prior to calling _init, _commit, _helmoltz etc is stored elsewhere, and is copied into the arrays prior to calling the solver, this is being carried out  correctly. Perhaps my first message was misleading. The second call to the poisson solver takes place completely independently of the first, different input, different BCs. I was only describing this to justify the calling the solver with only 'f' defined on the boundary.

Prior to calling the solver the 'f' array (heat flux q in your analogy) is set to zero everywhere except for the values directly on the boundary. Then comes the issue. The solver executes without any problems, and the results (for Temperature in your analogy) are of course as expected dependent upon the value of the boundary conditions (defined in bd_+ arrays), however appear to be independent of the value of 'f' set at input.

If it is useful I can attach an example data set + visualization of the output, but I think the description above probably suffices.

Thanks again for the quick response and any help you can provide.

[mecej4]

There is not enough detail in your description for me to understand what you are trying to do -- are you solving (i) one Poisson problem in several steps, or (ii) are you solving a sequence of different Poisson problems but using the solution of one problem as the trial solution to the next problem?

It would be useful to see the mathematical problem being defined before we can talk about obtaining approximate solutions.

I do not understand your phrase "compact support" in this context.

Please note that, depending on how the discretization is performed, the values of the source term f on the boundary nodes many never be needed or used. For example, consider a standard Poisson problem on a rectangle with Dirichlet B.C. The values of the dependent variable on the boundary nodes being known, the discrete approximation to the PDE is written down only at the interior nodes; thus the values of the source term on the boundary are neither needed nor is there any means to use any values that you specify on the boundary.

[Sav__Joe]

Hi again,

For all intents and purposes I solving 1) one Poisson Problem. Apologies for any unclarity in my explanation.

Yes I agree that the dependent variable being known at the boundary needn't be solved for, however the discrete approximation within the boundary depends on f on the boundary. This dependency I do not see in my output. When I change the value of f on my boundary (remember the problem I am solving has f zero everywhere except ON the boundary), the solution within the domain remains unchanged.

[mecej4]

Sav__Joe, you wrote: "...however the discrete approximation within the boundary depends on f on the boundary". 

Sorry, that is not true at all. Consider the singular case where f = 1 on the boundary and f = 0 in the interior. Remember that f is a source term, representing some conserved quantity per unit volume. Since the volume of any boundary is zero, in terms of an overall balance, the contribution of f to mass, energy, charge, etc., is also zero.

To convince yourself about this issue, take a square domain, say x = -1..1, y = -1..1, and place one interior node, 'o', at (0, 0) and four boundary nodes at the midpoints of the four edges, and name these boundary nodes 'e', 'w', 'n' and 's'. The five-point approximation to the Poisson equation for the dependent variable u is then 

    u_e + u_w + u_n + u_s - 4 u_o = - f_o

Do you see any presence of f values on the boundaries there?

[Sav__Joe]

Ok thankyou for the enlightening illustration. So I can conclude that the boundary values of f in the input have no influence on the solution?

[mecej4 · Accepted solution]

There do exist high accuracy finite-difference approximations, some of which may use scores of nodes rather than just five. Some of these approximations may use f values at the boundary.

For the MKL solver, however, as your experimentation has shown, that does not seem to apply.

[Sav__Joe]

Hello again,

The problem herein lied in the fact that I had assumed that there was some sort of additional check on the boundary to ensure continuity with the value of the solution within the domain, aka within the one-sided differences used in the discrete Laplacian approximation to the values of f on the boundary were somehow being used.

This would imply the solver also ensured that the Poisson equation was being satisfied on the boundary of the domain in some sense. Upon closer inspection of the solver reference page however I see that the approximation is only correct for x_a < x < x_b, and hence not x_a <= x <= x_b. 

Thanks again @mecej4 for the helpful explanation and of course also for the very quick replies.

[Kirill_V_Intel]

Hi!

If I got it right, you wanted to impose the Poisson equation at the boundary of the domain. This is just not how a boundary value problem for a Laplace operator is formulated.

What the physical representation of your problem? Do you have "heat" (or whatever it is that you want to find) sources at the boundary of the domain? Are they local point-wise sources? If so, I think you need to first think whether the problem is well-posed. Typically, when solved numerically this at least causes convergence issues (your solution might have irregularities around such points).

Best,
Kirill