dolfinx.fem.petsc.LinearProblem for nest and block systems

[jorgensd]

With the unification of create_*, assemble_* and assign for various PETSc operators, I suggest we let LinearProblem support non-blocked and blocked problems, i.e.

problem = LinearProblem(a, L, [bc0, bc1], petsc_options={"ksp_type": "preonly",
                                                         "pc_type": "lu",
                                                         "pc_factor_mat_solver_type":
                                                           "mumps"})
problem = LinearProblem([[a00,a01],[None, a11]], [L0, L1], bcs=[bc0, bc1],
                        u=[uh0, uh1])

For now, I've chosen to force the user to specify u for blocked problems, as it is hard extract the correct function spaces (dolfinx.fem.FunctionSpace, python class) from the nested list of forms with the current functions (extract_function_spaces work on the C++ function-spaces).

This also introduces an immediate solution vector _x, in the same fashion as within the proposed SNES solver.

It also adds the option of sending in a preconditioner to LinearProblem, as a ufl.Form or a nested sequence of Forms, which is then set to the krylov-solver.

[francesco-ballarin]

Looks like a reasonable interface to me. The only issue I see: what if the user wants a MatNest assembly?

[jorgensd]

Looks like a reasonable interface to me. The only issue I see: what if the user wants a MatNest assembly?

MatNest can't be used for direct solves (as far as I am aware).

Edit:
Looking at: https://petsc.org/release/overview/linear_solve_table/#direct-solvers
it seems like it could be possible. @jhale , @garth-wells what do you think?

[jorgensd]

Looks like a reasonable interface to me. The only issue I see: what if the user wants a MatNest assembly?

Context: firedrakeproject/firedrake#431 (comment)
even if the tables says it is supported, it is unclear if that means that it is supported only on each sub-matrix.

[francesco-ballarin]

MatNest can't be used for direct solves (as far as I am aware).

I do understand that, but in principle the user could choose something different from

petsc_options={"ksp_type": "preonly",
                                                         "pc_type": "lu",
                                                         "pc_factor_mat_solver_type":
                                                           "mumps"}

that is compatible with MatNest. I do understand that it will be a kind of advanced user that probably will have no trouble to set up the KSP object on their own, but still it would feel odd to have Nest-compatible SNES and Nest-incompatible linear solver.

[jhale]

On MatNest + LU/MUMPS, try it and see what happens? We could simply return an error message if the user wants a nest solve and passes no KSP options - clearly this requires problem specific help that we cannot automate.

Passing u isn't ideal - could we write something that would make it work without?

[garth-wells]

The logic around kind looks complicated. I'd suggest passing kind to initialiser. If the user asks for a matrix type and a solver that are not compatible, let PETSc raise the exception.

[jorgensd]

I've made "kind" input to LinearProblem.
It now accepts "nest" as well, and it creates the functions that stores the solution, regardless of the type of system (if not supplied by the user).

I've added a test that check that all modes work in serial/parallel (here we can use mumps+nest as it can invert each of the systems individually).

[garth-wells]

Looking at his more closely, it's not clear to me what need this change addresses. Block/nest systems are usually constructed to support block-type preconditioners. But, it's not obvious that the changes support this. E.g., how should the standard block-preconditioned solvers implemented in the Stokes demo be implemented via the LinearProblem interface?

[jorgensd]

Looking at his more closely, it's not clear to me what need this change addresses.

For mixed dimensional/submesh problems, one has to use blocked problems (either manually or through ufl.MixedFunctionSpace), where one doesn't necessarily use blocked-preconditioning.

s. But, it's not obvious that the changes support this. E.g., how should the standard block-preconditioned solvers implemented in the Stokes demo be implemented via the LinearProblem interface

We can of course add a set_preconditioner function to the LinearProblem, which creates the P matrix, and sets it to the krylov solver. I guess it could even be an input argument to the constructor.

[jorgensd]

Looking at his more closely, it's not clear to me what need this change addresses. Block/nest systems are usually constructed to support block-type preconditioners. But, it's not obvious that the changes support this. E.g., how should the standard block-preconditioned solvers implemented in the Stokes demo be implemented via the LinearProblem interface?

I've now added the option to pass in a preconditioner to the LinearProblem at construction.