A new variant of solvePDE (renaming the original variant)

[mhvwerts]

A new way of calling solvePDE is implemented, preserving also the existing solvePDE implementation. I am not a big fan of programming multiple dispatch mechanisms in Python, but I used it here to allow for two ways of calling solvePDE.

The new variant hides the matrix algebra, and works directly on the solution variable and the various terms. Look at the update example code in README.md to get an idea. I think it looks quite logical. Supplying the boundary condition term as a separate parameter enforces its essential presence, and makes that the list of equation terms only contains, well, equation terms. The solution variable is updated automatically (if you want to keep the "old" value, you should store it separately before calling this variant of solvePDE).

The implementation is simple, and did not require any other changes to the module code. There are some potential problems with this simple implementation, but will probably be easy to fix if we come across them. I left some comments in the code how this may be achieved.

As a test, I worked with my favorite 1D convection-diffusion equation (test_pdesolver_mason_weaver.py), which should some day be converted into a nice example Notebook. Mass is indeed well conserved by PyFVTool.

P.S. PyFVTool is now version 0.2.1!

[mhvwerts]

P.S. The potential problem will occur when there is a minus sign before a term which is a tuple of an M matrix and a RHS vector. This could be simply solved by introducing a 'signed tuple' subclass of tuple. But I am not sure if we will ever run into this problem.

[simulkade]

Having boundary conditions as an input to the solver is an excellent idea. The solver is also an excellent step in the right direction. The reason I initially wrote it to deal with matrices directly was a personal need: I had a system of coupled PDEs that I wanted to solve simultaneously. I ended up with two linearized equations that could be combined, solved, and separated as follows (PDE 1 and 2 discretized for variables p and sw):

M = vstack(
                    [
                        hstack([Mdiffp1 + Mbcp, Mconvsw1]),
                        hstack([Mdiffp2, Mconvsw2 + Mtranssw2 + Mbcsw]),
                    ]
                )
RHS = np.hstack([RHS1 + RHSbcp, RHS2 + RHStrans2 + RHSbcsw])
x = spsolve(M, RHS)
p_new = np.reshape(x[0 : (Nxyz + 2)], (Nxyz + 2))
sw_new = np.reshape(x[(Nxyz + 2) :], (Nxyz + 2))

As you can see, assembling the matrix of coefficients of the coupled system follows a simple rule and can be done automatically. However, I never thought about having a solver that can accept these equations and return the solution. I think we can have a different solver in mind for the coupled system (e.g. solvePDESystem and also potentially think about a new PDE or DiscretizedPDE class that contains the matrices of coefficients and RHS vectors for each equation for each primary variable. In the meantime, I am in favor of keeping the simple matrix-based solver perhaps with a different name to reflect its rather low-level implementation (solveLinearizedPDE?). The new solver then can be used under the name solvePDE, because its usage feels much less awkward than the matrix-based solver. What do you think, @mhvwerts ?

[mhvwerts]

With your remark about coupled equations, you actually started answering an other question that is emerging in the back of my mind about coupled systems and to which I hope to come back soon. Thanks!

With this PR, it is not my intention to completely hide the matrix algebra, on the contrary. Direct access to the matrix equation allows to create specialized numerical schemes, e.g. for coupled equations. Also, it is instructive to work directly with the matrix equations, even calling the matrix solver (spsolve) like you do without passing through a convenience function like solvePDE.

For the sake of readability for new users and ease-of-use for the simplest PDEs, I thought of this new variant of solvePDE. I completely agree that there are many interesting possibilities for creating specialized routines and classes (as a long term goal, with also the trade-off between coding effort and ease of use). In particular, solvePDESystem is a nice short-term perspective.

If I understand you correctly, you are proposing a way out of the multiple dispatch maze for solvePDE:

• use solvePDE exclusively for the new-style 'easy' solvePDE
• rename the previous matrix version to solveMatrixPDE (I propose to use this name instead of solveLinearizedPDE. @simulkade, do you agree?)

I can easily change this with some well-targeted find & replace in this PR.

From there, we can imagine a entire collection of new and improved solvers, but also still work directly with the matrix equations.

[simulkade]

You read my mind. I was thinking about a way out of multiple dispatch. It is OK to have several solvePDE functions/classes with different names. It feels a bit matlabesque to constantly check for input type (or maybe it is my feeling) :-). Agree with everything else including the new name solveMatrixPDE name for the old solvePDE.

[mhvwerts]

Multiple dispatch can have its utility in some use cases, justifying the use of clunky constructs, but is best avoided in Python because it leads to code that is not always very readable. You are right. Especially with these solver routines it is much clearer and straightforward to have differently named routines. More routines can now easily be added in the future.

Anyway, I think this is now ready to be merged. Thanks for the fruitful discussion.

[simulkade]

Thanks, @mhvwerts. It was indeed a fruitful discussion.