Add matrix free demo using PETSc

[jorgensd]

Requires: #3937

A demo showcasing how to use PETSc's Python type matrix to create a matrix-free operator.
Would partially resolve: #3779

A separate question is where one should put the MatrixFreeOperator, as it doesn't change with the bi-linear form.

[michalhabera]

I don't like that you go for a mixed problem and in addition try to solve it with two approaches. It obscures the core of matrix-free demo.
Also, why bother with BCs for projection problem?

[jorgensd]

I don't like that you go for a mixed problem and in addition try to solve it with two approaches. It obscures the core of matrix-free demo. Also, why bother with BCs for projection problem?

Because it is non-trivial to set up, as seen in the demo. It already exposed a bug in DOLFINx (the bc one). It also serves as a general matrixfree solver, rather than one that is specialised for projections.

I can move this demo to my own web page if that is more suitable.

[michalhabera]

I don't like that you go for a mixed problem and in addition try to solve it with two approaches. It obscures the core of matrix-free demo. Also, why bother with BCs for projection problem?

Because it is non-trivial to set up, as seen in the demo. It already exposed a bug in DOLFINx (the bc one). It also serves as a general matrixfree solver, rather than one that is specialised for projections.

I can move this demo to my own web page if that is more suitable.

I think that PETSc's

I don't like that you go for a mixed problem and in addition try to solve it with two approaches. It obscures the core of matrix-free demo. Also, why bother with BCs for projection problem?

Because it is non-trivial to set up, as seen in the demo. It already exposed a bug in DOLFINx (the bc one). It also serves as a general matrixfree solver, rather than one that is specialised for projections.

I can move this demo to my own web page if that is more suitable.

I think that PETSc' matrix free solver and diagonal assembly should be demoed. But why not keep it as concise, readable and simple as possible now, and then maybe add some form of MatrixFreeOperator into our petsc namespace later?

[IgorBaratta]

I like the idea, much better than https://github.com/FEniCS/dolfinx/tree/main/cpp/demo/poisson_matrix_free.
But also think the demo is a bit complex, and would suggest to make the demo less general and more concise.
And perhaps add MatrixFreeOperator to the core later as @michalhabera suggested.

[jorgensd]

I think that PETSc' matrix free solver and diagonal assembly should be demoed. But why not keep it as concise, readable and simple as possible now, and then maybe add some form of MatrixFreeOperator into our petsc namespace later?

Adding it later means that I need to keep on maintaining it somewhere (as later might be a year from now).

I additionally thought that we wanted more advanced demos in DOLFINx, not just simple setups.

[schnellerhase]

Adding it later means that I need to keep on maintaining it somewhere (as later might be a year from now).

How about moving the general interface of the MatrixFreeOperator into the tests, as done for NonlinearPDEProblem and NonlinearPDE_SNESProblem?

[jorgensd]

Adding it later means that I need to keep on maintaining it somewhere (as later might be a year from now).

How about moving the general interface of the MatrixFreeOperator into the tests, as done for NonlinearPDEProblem and NonlinearPDE_SNESProblem?

Moving it to the tests just makes it tricky for people to find it (and keep track of changes). That’s how we ended up with many poor copies of the SNES interface before finally adding it to the «core» PETSC interface

[jhale]

The right approach for the demo depends on what you want to achieve:

1. An example showing how to interface DOLFINx with PETSc at a lower level - in this case, the example should probably be simplified.
2. An example showing how to use a fully featured and library based NonlinearOperator - in this case, the code should probably go in the library.

On the second point, some lower-level functions like https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/fem/petsc.py#L1032 might be a better way to expose some of the routines for assembling actions/diagonals - these could then be placed into a context within the user's scope.

[jorgensd]

The right approach for the demo depends on what you want to achieve:

1. An example showing how to interface DOLFINx with PETSc at a lower level - in this case, the example should probably be simplified.

2. An example showing how to use a fully featured and library based `NonlinearOperator` - in this case, the code should probably go in the library.

On the second point, some lower-level functions like https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/fem/petsc.py#L1032 might be a better way to expose some of the routines for assembling actions/diagonals - these could then be placed into a context within the user's scope.

The main complication is that setting bcs in the diagonalized mat-vec (mult) case for PETSc does not fit our current application of BCs in assemble_residual: https://github.com/FEniCS/dolfinx/pull/3938/files#diff-0b48cc54bef51011d79955f12bee4a1caae04c7ddf91564444e148a592049c5fR194-R216

I would argue that a 150 Line class (MatrixFreeOperator, https://github.com/FEniCS/dolfinx/pull/3938/files#diff-0b48cc54bef51011d79955f12bee4a1caae04c7ddf91564444e148a592049c5fR194-R216) is not very complex. It is an advanced demo, for advanced users that know what they are doing.

I struggle to understand where we want to draw the line on "more complicated"/"too complicated" for our demos.

[jhale]

Well, I would leave the complication of handling mixed or block spaces to the library version.