Calling external sparse direct solvers

[mhvwerts]

It is an important feature of PyFVTool that it works with "pure scientific Python" (i.e. just python, numpy, scipy, and matplotlib for visualization), especially for sharing work with students and colleagues.

However, it may be anticipated that for certain future users, it would be interesting to provide an easy method to call an external sparse solver instead of scipy.sparse.linalg.spsolve, without sacrificing the default "pure scientific Python" character of PyFVTool.

Indeed, I noticed when doing some work on more realistic models, that my PyFVTool only uses a single thread on my processor when solving . Not a serious issue now, the calculations still only take less than 5 min, and it teaches me to be thoughtful on planning calculations. BTW, the 2D cylindrical results obtained with PyFVTool are really nice, so I am happy.

Just opening this to collect ideas on how to provide for a transparent external solver calling mechanism, and which external solvers might be easily accessible. A quick search yielded, for instance, PyPardiso

[mhvwerts]

By default scipy.sparse.linalg.spsolve appears to use the SuperLU solver, which is always available on scipy. If UMFPACK is available in the Python environment (via SciKits), then it is expected to use that.

However, UMFPACK on Windows is not an easily accessible option:

• https://github.com/scikit-umfpack/scikit-umfpack/issues/85
• https://github.com/scikit-umfpack/scikit-umfpack/issues/64

I will stick with the vanilla scipy.sparse.linalg.spsolve for now.

[mhvwerts]

Easiest mechanism for accessing an external solver may be to add an optional keyword argument to solvePDE(), e.g. externalsolver with default None or otherwise a function having the same interface as scipy.sparse.linalg.spsolve.

I could look into that in some future.

[simulkade]

It is a very important point indeed. A solver I really like to add is pyamg. I will have to some tests on its installation on different platforms and how to properly call it.

[mhvwerts]

I quickly tried something while taking a short break. I got the external solver calling mechanism working, but ran into a problem with pypardiso. Probably not a big problem (something to do with scipy.sparse switching to arrays instead of matrices; some conversion will need to be done before calling pypardiso).

The external solver mechanism does work when providing scipy.sparse.linalg.spsolve as an "external solver".

See the corresponding example notebook in my fix_details branch: https://github.com/mhvwerts/PyFVTool/tree/fix_details

[simulkade]

I think the new solvePDE looks fine, and as you correctly said, does not add new dependencies. We can later make a list of direct linear solvers that can be installed by the user. I have had positive experiences with MUMPS too, if there is an up-to-date Python interface for it.