SNES solver interface for standard, block and nest systems

[jorgensd]

A "functional" approach to the snes problem interface.
This reduced the code to:

• One factory function: create_snes_solver, that returns the the user a SNESObject and a vector to solve into x
• A SNESSolver, that handles destruction of PETSc objects. It also store convenient wrappers to transfer data from your dolfinx.fem.Function|list[dolfinx.fem.Function (u) to the PETSc vector used for solving (x)

The tests now show examples of:

1. How to use SNESSolver, the highest level interface
2. How to use create_snes_solver, the mid-level interface
3. How to use the functions supplied in dolfinx.nls.petsc and dolfinx.fem.petsc to interact with your own snes solver instance.

[jorgensd]

Why not just a dictionary that maps AssemblyType to a vector and matrix factories and let user call these? Maybe as a user I want to only create the matrix, because I already have a vector in the correct layout. And maybe, why the need for a helper at all

If you want to use the lowest level API the users can call the appropriate create_vector_*, create_matrix_* and shouldn’t really need the create_snes_data_structures.

[jorgensd]

There should be one-- and preferably only one --obvious way to do it."

The others wanted a low level API. We can of course remove the mid-layer, ie not expose create_snes_solver, and only expose the wrappers for F,J and moving data, and then let anyone that wants to do something else than SNESSolver do everything manually.

[jorgensd]

like the layered approach and I like the high-level SNESSolver class. But I am missing the point in this specific design of the lower- or mid- level API. What is the typical use case there? There should be an example. In the current code you help users with creating matrices create_snes_matrices_and_vectors and you provide some dolfinx.nls.petsc.F_nest that user controls what data is passed into it. But when would this very API be needed?

The low level API can for instance be used to set your own F_nest,J_nest that takes into account the external operator. However, you might not want to manually control all the other things by hand, such as matrices/vec creation.

[jhale]

I think I struggle with create_snes_matrices_and_vectors in particular. Why not just a dictionary that maps AssemblyType to a vector and matrix factories and let user call these? Maybe as a user I want to only create the matrix, because I already have a vector in the correct layout. And maybe, why the need for a helper at all? The more layers in-between you expose to public API, the more ways to do things appear, and "There should be one-- and preferably only one --obvious way to do it."

Seems like a good suggestion to break this out further.

I like the layered approach and I like the high-level SNESSolver class. But I am missing the point in this specific design of the lower- or mid- level API. What is the typical use case there? There should be an example. In the current code you help users with creating matrices create_snes_matrices_and_vectors and you provide some dolfinx.nls.petsc.F_nest that user controls what data is passed into it. But when would this very API be needed?

My view is that SNESSolver will end up being re-born as NonlinearProblem (in the same way LinearProblem wraps KSP), but that's another discussion to be had.

So then in terms of using petsc4py SNES, there will be the create_snes_solver which will be the high-level API which returns a fully setup petsc4py SNES, and the low-level components that can be put together as the user wishes by interacting with petsc4py SNES directly. The idea here is that we reduce boilerplate - for example, I have applications where we simply need to redefine the function J.

Edit: @michalhabera There is an "example" of using SNES directly at https://github.com/FEniCS/dolfinx/pull/3648/files#diff-f9b6597717f57bdb47379a64d265490d51aa09a0320d4d47511b09fba2371b1bR304

[jhale]

If we can agree that SNESSolver should be NonlinearProblem then this call should go inside solve and also solve should do any necessary ghost updates, in the same way that LinearProblem works.

[jhale]

Let's do that in this PR so that we can review this important set of changes as a whole.

[michalhabera]

I think I struggle with create_snes_matrices_and_vectors in particular. Why not just a dictionary that maps AssemblyType to a vector and matrix factories and let user call these? Maybe as a user I want to only create the matrix, because I already have a vector in the correct layout. And maybe, why the need for a helper at all? The more layers in-between you expose to public API, the more ways to do things appear, and "There should be one-- and preferably only one --obvious way to do it."

Seems like a good suggestion to break this out further.

I like the layered approach and I like the high-level SNESSolver class. But I am missing the point in this specific design of the lower- or mid- level API. What is the typical use case there? There should be an example. In the current code you help users with creating matrices create_snes_matrices_and_vectors and you provide some dolfinx.nls.petsc.F_nest that user controls what data is passed into it. But when would this very API be needed?

My view is that SNESSolver will end up being re-born as NonlinearProblem (in the same way LinearProblem wraps KSP), but that's another discussion to be had.

So then in terms of using petsc4py SNES, there will be the create_snes_solver which will be the high-level API which returns a fully setup petsc4py SNES, and the low-level components that can be put together as the user wishes by interacting with petsc4py SNES directly. The idea here is that we reduce boilerplate - for example, I have applications where we simply need to redefine the function J.

Edit: @michalhabera There is an "example" of using SNES directly at https://github.com/FEniCS/dolfinx/pull/3648/files#diff-f9b6597717f57bdb47379a64d265490d51aa09a0320d4d47511b09fba2371b1bR304

But I am missing the point in this specific design of the lower- or mid- level API. What is the typical use case there? There should be an example

Please see the current tests. They show all different usages of the API.

I've seen the tests and "how" to use the API. But my question was more about "why" to use it, in the tests all could be neatly done with the high level SNESSolver.

[jorgensd]

I think I struggle with create_snes_matrices_and_vectors in particular. Why not just a dictionary that maps AssemblyType to a vector and matrix factories and let user call these? Maybe as a user I want to only create the matrix, because I already have a vector in the correct layout. And maybe, why the need for a helper at all? The more layers in-between you expose to public API, the more ways to do things appear, and "There should be one-- and preferably only one --obvious way to do it."

Seems like a good suggestion to break this out further.

I like the layered approach and I like the high-level SNESSolver class. But I am missing the point in this specific design of the lower- or mid- level API. What is the typical use case there? There should be an example. In the current code you help users with creating matrices create_snes_matrices_and_vectors and you provide some dolfinx.nls.petsc.F_nest that user controls what data is passed into it. But when would this very API be needed?

My view is that SNESSolver will end up being re-born as NonlinearProblem (in the same way LinearProblem wraps KSP), but that's another discussion to be had.
So then in terms of using petsc4py SNES, there will be the create_snes_solver which will be the high-level API which returns a fully setup petsc4py SNES, and the low-level components that can be put together as the user wishes by interacting with petsc4py SNES directly. The idea here is that we reduce boilerplate - for example, I have applications where we simply need to redefine the function J.
Edit: @michalhabera There is an "example" of using SNES directly at https://github.com/FEniCS/dolfinx/pull/3648/files#diff-f9b6597717f57bdb47379a64d265490d51aa09a0320d4d47511b09fba2371b1bR304

But I am missing the point in this specific design of the lower- or mid- level API. What is the typical use case there? There should be an example

Please see the current tests. They show all different usages of the API.

I've seen the tests and "how" to use the API. But my question was more about "why" to use it, in the tests all could be neatly done with the high level SNESSolver.

This was a request of Jack. As he stated, he wants to work directly with the SNES-object, and handle memory desctruction, but maybe only replace the F or J.

[jhale]

I've seen the tests and "how" to use the API. But my question was more about "why" to use it, in the tests all could be neatly done with the high level SNESSolver.

Three (edit) points here:

1. The basic principle in DOLFINx is that we return first-class petsc4py objects, except at a couple of high-level entry points (LinearProblem and the new proposed NonlinearProblem). I think this has been very successful. Personally, as I've discussed with Jorgen, I will probably not use LinearProblem and NonlinearProblem as I prefer to work with petsc4py directly, but I appreciate that this is not very welcoming.
2. Your point is correct, these tests could be implemented using SNESSolver/NonlinearProblem. They are not really meant to be examples of use, but experiments + tests of how the lower-level API looks + works. I'm inclined to agree with you that create_snes_matrices_and_vectors could be removed from the public API.
3. As an example of use, I will be able to strip down what we are providing in ExternalOperator to three function-based implementations of the residual computation function F (standard, nest, block). I no longer have to maintain or write classes/subclasses, I can just use the first-class functionality to make 95% of a SNES, and instruct users they need to call setFunction with a special F.

[jhale]

The thinking behind returning x is that users must not under any circumstances use the data behind u to call SNES solve. This may be signposting that might not be necessary for users of the lower-level interface. On the other hand, easy mistake to make.

[garth-wells]

A first pass of comments.

I would prefer to see this built up in stages, perhaps by starting with bottom-up functionality that simplifies the test examples. An example to start could be a single function copies data from DOLFINx Functions in/out of blocked PETSc data structures, and get that as simple as possible. This could be merged aheads of a SNES interface.

[garth-wells]

It's not clear from the docs what is 'computed'. It should say something like it computes the derivative of a UFL form to generate a UFL form for the Jacobian.

Need to document the arguments.

[jhale]

I'm going to review and edit all of the documentation in due course. Many issues, also with consistency with respect to LinearProblem.

[garth-wells]

I think assembly_type should be a PETSc Mat/Vec type. This would allow the caller to not only specify 'nest`, but also other PETSc matrix types.

[jhale]

Already discussed this, it was also my preference, but @jorgensd told me it's not possible to discern between nest and block types. I do see MATNEST and MATBLOCK though: https://petsc.org/release/manualpages/Mat/MatType/ - would these work @jorgensd ?

[jorgensd]

We don't use createBlockMat (https://petsc.org/release/manualpages/Mat/MatCreateBlockMat/).
We use createMat and attach an index set to it based on the concatenated index maps.

[garth-wells]

We don't need to discern between block and nest types. If the forms are a list, then it's block or nest. Nest is used if the user asks for nest, otherwise MATAIJ or whatever else they want.

The interface should probably have options for the vector and the matrix (matrices for nest), since they can differ. Even for a single matrix or a matrix in a nest there are options, e.g. block CSR. See

dolfinx/python/dolfinx/fem/petsc.py

Line 174 in 192222e

def create_matrix(a: Form, mat_type=None) -> PETSc.Mat:

where the users passes the MatType. The other "create_matrix functions should have that option too.

[jorgensd]

So you want:

1. Send in two arguments, "MatType" and "VecType" to the SNESSolver-constructor.
2. If the user sends in any mat-type but nest, check if the residual is a list or not.
   2a. If it is a list use assemble_block and corresponding methods
   2b. OTherwise use assemble etc.
   What should happen if:

• A user sends in snes as argument, but a single form (non-blocked).
• The user sends in a non-matching argument to MatType and VecType.

[garth-wells]

Focusing on the matrix, the user should do one of the following:

1. Send in None for the matrix type. In this case the default PETSc matrix type will be used. If the forms are blocked, the assembly will be blocked; or
2. Send in MATNEST as the matrix type. In this case each block will have the default matrix type; or
3. Send in a list-of-lists with the matrix type for each block in a MATNEST matrix.

The C++ interfaces for creating the matrices support these options. The code should handle a PETSc prefix to support setting of matrix types from the command line.

[garth-wells]

The F_foo function names looks weird.

[jorgensd]

What should we call them then?
Just calling them weird isn't very helpful.

[jorgensd]

A first pass of comments.

I would prefer to see this built up in stages, perhaps by starting with bottom-up functionality that simplifies the test examples. An example to start could be a single function copies data from DOLFINx Functions in/out of blocked PETSc data structures, and get that as simple as possible. This could be merged aheads of a SNES interface.

I find this an unreasonable request.
This has been now reviewed several times by several core developers. This does not interfere with the core of DOLFINx. The API can be refined prior to release and it simplifies and enhances the user workflow. There are way bigger internal changes that get merged on a monthly basis with way less reviewing.

[garth-wells]

A first pass of comments.
I would prefer to see this built up in stages, perhaps by starting with bottom-up functionality that simplifies the test examples. An example to start could be a single function copies data from DOLFINx Functions in/out of blocked PETSc data structures, and get that as simple as possible. This could be merged aheads of a SNES interface.

I find this an unreasonable request. This has been now reviewed several times by several core developers. This does not interfere with the core of DOLFINx. The API can be refined prior to release and it simplifies and enhances the user workflow. There are way bigger internal changes that get merged on a monthly basis with way less reviewing.

There are concerns on the design. It will be easier, and faster, to get the design right by building it up in steps, bottom up.

[jhale]

My impression is that we should first tackle the comment made by @garth-wells https://github.com/FEniCS/dolfinx/pull/3648/files#r1973997836 in a separate PR.

This can then be refactored to use #3659 which should remove many of the if statements regarding assigning functions to vectors and vice-versa. We can then think about if there are any smart ways to deal with Mat and Vec types in the SNES solver setup.