# Request algorithm recommomendation for 2D generalised Poisson solution

I need to examine the (static) electric field distribution around various electrode configurations in the presence of a dielectric, and have been using a finite difference approach to the 2D generalised Poisson equation do this. James Nagel's nice paper at http://www.drjamesnagel.com/notes/Nagel%20-%20Numerical%20Poisson.pdf has been my guide and suggests SOR on a discretized integral form of Poisson's equation. But Nagel's paper also notes that for large dielectric constants - which I have - SOR converges slowly or fails: then he uses the Matlab/Octave "\" as a blackbox operation on the sparse system matrix. This works fine for me on problems up to 300x300 grid or thereabouts - but it is slow and I really want to use C or something similar to integrate with other routines. So I am seeking advice on a suitable algorithm to apply to this problem. I think ADI may be what I need - but before learning how to code and use it, I'd value guidance and perhaps alternative recommendations. Perhaps I just need to call the LAPACK version of "\"? Of course, I want unconditional stability, high accuracy and instant convergence! Thanks for reading!

There are many excellent codes out there that solve the Poisson equation. I may point you at step-6 of the deal.II tutorial at http://dealii.org/developer/doxygen/deal.II/step_6.html as just one example. (Disclaimer: That's my own code.) Most other finite element libraries will have similar example codes.

• Many thanks for your advice, Wolfgang. I'd rather avoid FEM as I am not familiar with the underlying math and the learning threshold is real - though I understand it is a "better" method in many ways and I've bookmarked your tutorial to work through in due course. – user2093190 Jul 7 '13 at 16:28

A good way to approach this is by using existing powerful numerical packages such as PETSc. Look for PETSc examples of solving the Poisson equation, and that will be a good starting point.

• Thanks Maxim, I didn't know about PETSc - a quick look shows it is well worth looking further. I'll get to grips with a tutorial. – user2093190 Jul 7 '13 at 17:16

The FDM paper (which I personally wrote, BTW) is intended more for basic education than it is for large-scale implementations. As you may have noticed, the major bottleneck for large domains is setting up the system matrix and inverting it. In my experience, Matlab simply runs out of memory when trying to run the inversion outright, so you have to be clever how this step is implemented. I did my best to try and streamline this step in the code, but getting a solid matrix inversion for this kind of computation is a numerical challenge unto itself. If anyone has some clever updates for allowing my codes to process larger matrices, I'd love to hear about it.

... or of course you can always just use FEM instead, and simply sample the grid more efficiently. :)

EDIT: I'm not sure about the differences between Matlab and Octave when dealing with sparse matrices. In Matlab, setting up the matrix has always been the easy part for me. It has always been the inversion step that kills it by running out of memory.

• James, thanks for your paper. As reported, direct use of the umfpack library is as good as it gets, AFAIK. Those routines seem to be extremely well optimised and memory efficient. – user2093190 Oct 6 '14 at 16:31

I suggest using the "Run and Time" option in MATLAB to see exactly where the time is being spent in your 300x300 grid example. I am almost certain it is not in the \ operator. Generally, MATLAB has very fast, state-of-the-art implementations for sparse matrix factorizations.

I took at look at James Nagel's MATLAB code at the site referred to in your link. I didn't try to run it since I couldn't find an example case. But I did notice that there is a loop over all points in the grid (e.g. 300x300) to fill the sparse matrix prior to calling . I strongly suspect this is where most of the time is being spent. It is generally more efficient in MATLAB to write the code so that operations are performed on long vectors rather than scalars.

You can find many references on this topic if you search for "vectorized MATLAB code".

I should add that if you are interested in a fast, basic finite element code for the Poisson equation written in MATLAB, you can find one at the following link

http://www.isima.fr/~jkoko/codes.html

along with a discussion of how it was designed.

• Bill, thanks for your constructive comments. I am sure you are right - much of the time is spent in setting up rather than solving the system. I am using Octave rather than Matlab (being poor!) and believe the "\" functionality is based on a classification of the matrix type followed by a call to the correct LAPACK routine. So understanding which LAPACK method to call, could be an easy route ahead. – user2093190 Jul 7 '13 at 16:41
• I agree vectorising the Octave code would be a good route to speed up the calculation. It's not a style that comes immediately and I am sure it can improve on my present code which uses traditional loops. But I'd rather migrate to C which I am more familiar with. – user2093190 Jul 7 '13 at 16:53

So, finding a good solution has been instructive. First point is that as Bill suggested, setting up the system in Octave (rather than solving it) takes most of the time - but it is not easy code to vectorise. Also, for my problem formulation with mixed boundary conditions, the system matrix is not symmetric, which rules out some efficient approaches. I looked at different libraries and now use the sparse LU solver in umfpack, included in my Linux system as part of the suitesparse library. I gather both Octave and Matlab use umfpack which I take as a recommendation - it is a traditional library, easy to call from C as well as C++, and has excellent, clear documentation. A straightforward C routine solves a 320x320 FDE grid >100 times faster than in Octave and easily copes with 1000x1000 in 30s on a cheap laptop.

FYI, I've recently published a more comprehensive tutorial with the IEEE APS magazine. Here is the direct link. This one includes more material and has been formally peer-reviewed, so publication references finally have something solid to point to.

It's a great starter paper to get the basic idea of FDM working. It also seems to work very well with domain sizes on the order of 300x300 or so, depending on your memory. As we've all seen, the only real bottleneck for larger domains is setting up that pesky system matrix and inverting it.

• Thanks for that, James - a clear paper at a valuable level of detail. Another good reference for other approaches to solution of the Poisson eqn is Lawrence Dworsky's book Introduction to Numerical Electrostatic Using MATLAB. Use of a library like umfpack has its own educational value. You learn little about the algorithm, but it's worth learning that sometimes a black box is the best route to an efficient solution. One aspect missing from most literature is the handling of anisotropic dielectrics - important for liquid crystals as well as some circuit substrates, I understand. Any hints? – user2093190 Feb 22 '15 at 13:10
• Anisotropic materials add a whole new level of complexity. Basically, you have to go all the way back to square-1 and modify the dielectric function in Gauss' law. Instead of a simple dielectric function in space, it has its own directional components that have to be tracked. Then you have to propagate that little change all the way back through the finite-difference algorithm. For simple uniaxial materials, it's probably not so bad, but it can easily get complicated for varying types of anisotropy. – James Nagel Mar 25 '15 at 19:52