multicolumn right-hand sides

[gdmcbain]

This is a nonissue in a sense; multicolumn right-hand sides just work with no modification required to the usual idiom of e.g. ex14 36d7ace

u[dofs] = solve(*condense(A, 0.*b, u, dofs))

where dofs is still a one-dimensional array of indices but now u and b are two-dimensional with u.shape[0] == b.shape[0] == A.shape[0] and u.shape[1] == b.shape[1] >= 1.

However, as that might not be obvious and as it is so useful in applications, this is being noted here as worth documenting and worth demonstrating in an example. (I hope to provide one shortly; I have it working in complicated cases but would like to strip it down to essentials before publishing.)

[kinnala]

I'm curious that what kind of applications you have in mind?

[gdmcbain]

Yes, sorry, a month has passed without my having come up with a nice potted example here.

So the ex14 referred to above was renamed ex13 in 366392b. It solves the Laplace equation with mixed boundary conditions: (i) unit potential, (ii) zero potential, (iii) zero current, and integrates the Joule heating over the domain to get the resistance R between the positive and ground edges. The resistance of something with two ports is usually thought of as a single quantity, giving the ratio of the signed voltage drop to the signed current. Here the role of the finite element method is condensation: to reduce the continuous multidimensional model to a lumped equivalent, a resistor. But what happens if there are more than two ports? Another way to think about the resistance is that there is a conductance matrix that multiplies the vector of port potentials to give the corresponding vector of port currents. In the case of a single (two-port) resistor, this is just [[+1, -1], [-1, +1]] / R. In general, the [i, j] coefficient of the conductance matrix gives the current into the i - th port given unit potential at the j - th port with the rest grounded. Thus one can calculate the matrix column-wise by solving the n boundary value problems which charge each port in turn. Instead of a single resistor, the reduced order model corresponds to a circuit that looks like a 'complete graph', with a resistor between each pair of distinct nodes with resistance minus one on the off-diagonal coefficient. It can be shown that the matrix has zero row and column sums, so actually only n – 1 of these problems needs to be solved. In the case n = 2, that's one, as in the given example, so there's not much point to this matrix formulation, but for larger n, being able to specify the n – 1 right-hand sides all at once is at least very tidy but also possibly more efficient if one is using a Krylov solver since the Krylov subspace can be shared by the n – 1 right-hand sides and solutions. (I haven't actually demonstrated the latter claim, but the tidiness is a plus enough for me so far.)

[kinnala]

I see now. This is an interesting use case. I'm glad if you are able to contribute.

[gdmcbain]

Perhaps this should be documented with the example from ‘Three ways to compute multiport inertance’ (hal-02159726), now in press.

[gdmcbain]

Meanwhile, the next question is whether this will still ‘just work with no modification required to the usual idiom’ if skfem.utils.solve is passed a solver from skfem.utils.solver_iter_krylov.

[gdmcbain]

No, it doesn't work with iterative solvers from scipy.sparse.linalg like cg or lgmres; they raise
ValueError: A and b have incompatible dimensions:

https://github.com/scipy/scipy/blob/e0cfa29e2fbe86f994924c0c7514ff5bbfffd514/scipy/sparse/linalg/isolve/utils.py#L95-L101

[gdmcbain]

Alternative implementations of GMRES like

• pyamg.krylov.gmres
• krypy.linsys.Gmres

only accept single-column right-hand sides. (And similarly for other Krylov solvers.) Perhaps there's a good reason for this…

[gdmcbain]

Revisting some code that did work back in 2018-06-28, I find that it doesn't any more; multicolumn right-hand sides no longer seem to be supported by skfem.utils.condense (or enforce). A few things have changed since then. And since it was subsequently learned that the idea of multicolumn right-hand sides isn't generally supported by Krylov solvers, it probably isn't worth pursuing further.

(Were this still an issue, I'd close it; not sure what the equivalent for a Discussion is. Deletion seems excessive.)

[gdmcbain]

No, no need to lock it, I shouldn't think. I guess it'll just sink as other discussions are voted up. Thanks.

[kinnala]

And let me add that going forwards I'm trying to keep issues at least somewhat "actionable". So after there is something that can be done to the codebase and explained in few sentences, then we can "Create issue from a discussion".

[kinnala]

Maybe moving this issue was a failure on my behalf, I did some mass moving of issues and caught something that are more "issue-like" than the others.

[gdmcbain]

Yes definitely. The new Discussions list was a really good innovation.

[gdmcbain]

Oops, I mean keeping only actionable issues is good. This issue would have been a Discussion, had there been Discussions back then.

[gdmcbain]

Revisting some code that did work back in 2018-06-28, I find that it doesn't any more; multicolumn right-hand sides no longer seem to be supported by skfem.utils.condense (or enforce).

Strike that, sorry, the code with the usual minimal modifications required by the intervening changes to condense (#228, &c.) still works fine.