Absolute and relative tolerances in `NewtonSolver`

[RemDelaporteMathurin]

Hi @jorgensd @finsberg thanks for making this package! I hope to use it soon to replace some of the code in FESTIM.

I was wondering if we could add a relative (or absolute?) tolerance to the blocked Newton solver in order to keep the functionality compared to the default Newton solver. Actually I'm not sure what the current tolerance is.

https://github.com/scientificcomputing/scifem/blob/f64874c878865914c5dbac0bf3f789a75ab066d2/src/scifem/solvers.py#L205-L208

If it's easy enough I'm happy to make a PR but I don't really know where to start

[finsberg]

Hi @RemDelaporteMathurin. Thanks for the kind words.

I agree that having the option to both have an absolute and relative tolerance would be good. I guess with the current setup (@jorgensd please correct me if I am wrong), the tolerance would be equivalent to absolute tolerance (ref https://github.com/FEniCS/dolfinx/blob/main/cpp/dolfinx/nls/NewtonSolver.cpp).

I also think that passing max_iterations to the __init__ function and the tolerance to the solve function is a bit confusing. Perhaps we should go for a more traditional parameter object, e.g

class NewtonSolver:
    def __init__(
        self,
        F: list[dolfinx.fem.Form],
        J: list[list[dolfinx.fem.Form]],
        w: list[dolfinx.fem.Function],
        bcs: list[dolfinx.fem.DirichletBC] | None = None,
        parameters: dict[str, Any] | None = None
    ):
        self.parameters = type(self).default_parameters()
        if parameters is not None:
            self.parameters.update(parameters)
        ...

    @staticmethod
    def default_parameters() -> dict[str, Any]:
        return {
            max_iterations: 5,
            petsc_options: {},
            error_on_nonconvergence:  True,
            beta: 1.0,
            abs_tol: 1e-6,
            rel_tol: 1e-6,
        }

PRs are of course welcomed :)

[RemDelaporteMathurin]

@finsberg would this be the way to have a relative error then?

correction_norm = self.dx.norm(0)
logger.info(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < atol:
    return i
if (correction_norm - old_norm)/correction_norm < rtol:
    return i
old_norm = correction_norm

If so I'm happy to work on a PR

[finsberg]

This looks correct. Perhaps we could also output the relative error in the logs. Great if you could work on a PR 😊

[jorgensd]

There are many ways to measure the norm. Here we use the norm on the increment of the solution, while DOLFINx checks the norm of the residual (by default, it can be customized to your own needs with set_convergence_check: https://github.com/FEniCS/dolfinx/blob/8261439905011d5cc6b6b50060dbe0b43dec5997/cpp/dolfinx/nls/NewtonSolver.cpp#L26-L50 We should probably do the same

[RemDelaporteMathurin]

Ok so I should call self._solver.residual0 as in the C++ code you sent?

[jorgensd]

If we want to fully emulate the C++ solver, we need to refactor the code such that there is residual evaluator function.

This is because: "incremental" means: Measure the convergence in the du update parameter (needs one Newton step at least). "residual": Measure the convergence in the RHS vector. Can happen before the first Newton iteration as the initial guess can solve the problem.

To implement these we would have to follow the C++ code quite closely.

[RemDelaporteMathurin]

Oh sorry I thought that this NewtonSolver was built upon the C++ class... Any pointers on how to compute the residual from the python interface?

[jorgensd]

One would have to refactor the solver a bit if we want an exact copy of the C++ one in Python.

[RemDelaporteMathurin]

Ok sounds like a big lift for me. I'll try to have a look and fiddle around but can't promise anything 😢

[finsberg]

I belive residual0 is just the residual from the first iteration (see https://github.com/FEniCS/dolfinx/blob/8261439905011d5cc6b6b50060dbe0b43dec5997/cpp/dolfinx/nls/NewtonSolver.cpp#L226), which is initialised to 0.0 (https://github.com/FEniCS/dolfinx/blob/8261439905011d5cc6b6b50060dbe0b43dec5997/cpp/dolfinx/nls/NewtonSolver.cpp#L71).

[RemDelaporteMathurin]

I belive residual0 is just the residual from the first iteration

ok so then it's this line that we need to be able to get in this solver https://github.com/FEniCS/dolfinx/blob/8261439905011d5cc6b6b50060dbe0b43dec5997/cpp/dolfinx/nls/NewtonSolver.cpp#L30

Right?

[finsberg]

If we want to use the same norm (which I guess we want), then you could do the following

import petsc4py.typing

residual = self.dx.norm(petsc4py.typing.NormType.NORM_2)

This would also make it more explicit I think. By the way, you can find a reference here: https://web.cels.anl.gov/projects/petsc/vault/petsc-3.20/docs/petsc4py/reference/petsc4py.typing.NormTypeSpec.html

[RemDelaporteMathurin]

Fixed in #50