Finite relative residuals for 0th Newton iteration in time-dependent problem

dolfinx
1

Hello everyone,

solving a time-dependent, nonlinear problem with fenicsx v0.7.3, at some time step I’m facing that the solver at converges with 0 iterations, since the relative residual was below tolerance already at Newton iteration 0. I’m wondering how this can be, since at Newton iteration 0, I would expect r (rel) = inf.

Looking at the source code at dolfinx/NewtonSolver.cpp at main · FEniCS/dolfinx · GitHub as proposed by @dokken in this post, I find the absolute residual is initialized to 0, thus the 0th relative residual is inf. But apparently, this is only the case for the first time step and not any further.

This behavior can be reproduced with @dokken 's hyperelasticity tutorial. Setting solver.convergence_criterion = "residual" one gets the following output for the first two time steps:

2024-07-10 11:42:07.983 (   3.574s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 0: r (abs) = 0.163333 (tol = 1e-08) r (rel) = inf(tol = 1e-08)
2024-07-10 11:42:08.123 (   3.714s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:08.453 (   4.044s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 1: r (abs) = 1330.75 (tol = 1e-08) r (rel) = 8.03268(tol = 1e-08)
2024-07-10 11:42:08.617 (   4.208s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:08.865 (   4.456s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 2: r (abs) = 73.9037 (tol = 1e-08) r (rel) = 0.446096(tol = 1e-08)
2024-07-10 11:42:09.016 (   4.606s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:09.259 (   4.850s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 3: r (abs) = 0.471433 (tol = 1e-08) r (rel) = 0.00284566(tol = 1e-08)
2024-07-10 11:42:09.411 (   5.002s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:09.665 (   5.256s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 4: r (abs) = 0.904561 (tol = 1e-08) r (rel) = 0.0054601(tol = 1e-08)
2024-07-10 11:42:09.824 (   5.414s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:10.071 (   5.662s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 5: r (abs) = 0.0010116 (tol = 1e-08) r (rel) = 6.10621e-06(tol = 1e-08)
2024-07-10 11:42:10.223 (   5.814s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:10.469 (   6.060s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 6: r (abs) = 2.32809e-05 (tol = 1e-08) r (rel) = 1.40528e-07(tol = 1e-08)
2024-07-10 11:42:10.625 (   6.216s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:10.876 (   6.467s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 7: r (abs) = 1.95203e-10 (tol = 1e-08) r (rel) = 1.17828e-12(tol = 1e-08)
2024-07-10 11:42:10.876 (   6.467s) [main            ]       NewtonSolver.cpp:255   INFO| Newton solver finished in 7 iterations and 7 linear solver iterations.
Time step 1, Number of iterations 7, Load [ 0.   0.  -1.5]
2024-07-10 11:42:11.103 (   6.694s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 0: r (abs) = 0.163333 (tol = 1e-08) r (rel) = 0.000985911(tol = 1e-08)
2024-07-10 11:42:11.254 (   6.845s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:11.508 (   7.099s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 1: r (abs) = 1053.16 (tol = 1e-08) r (rel) = 7.16324(tol = 1e-08)
2024-07-10 11:42:11.661 (   7.252s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:11.913 (   7.504s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 2: r (abs) = 44.9774 (tol = 1e-08) r (rel) = 0.305922(tol = 1e-08)
2024-07-10 11:42:12.068 (   7.659s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:12.326 (   7.917s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 3: r (abs) = 1.59719 (tol = 1e-08) r (rel) = 0.0108636(tol = 1e-08)
2024-07-10 11:42:12.485 (   8.076s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:12.760 (   8.351s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 4: r (abs) = 10.0883 (tol = 1e-08) r (rel) = 0.0686175(tol = 1e-08)
2024-07-10 11:42:12.915 (   8.506s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:13.169 (   8.760s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 5: r (abs) = 0.0336323 (tol = 1e-08) r (rel) = 0.000228757(tol = 1e-08)
2024-07-10 11:42:13.324 (   8.915s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:13.575 (   9.166s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 6: r (abs) = 0.116157 (tol = 1e-08) r (rel) = 0.000790063(tol = 1e-08)
2024-07-10 11:42:13.730 (   9.321s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:13.985 (   9.576s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 7: r (abs) = 4.33836e-06 (tol = 1e-08) r (rel) = 2.95082e-08(tol = 1e-08)
2024-07-10 11:42:14.146 (   9.737s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-07-10 11:42:14.402 (   9.993s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 8: r (abs) = 2.48803e-09 (tol = 1e-08) r (rel) = 1.69228e-11(tol = 1e-08)
2024-07-10 11:42:14.402 (   9.993s) [main            ]       NewtonSolver.cpp:255   INFO| Newton solver finished in 8 iterations and 8 linear solver iterations.
Time step 2, Number of iterations 8, Load [ 0.  0. -3.]

Apparently, Newton iteration 0 of the second time steps does not yield an infinite relative residual.

  1. Is this behavior expected?
  2. What formula is used to compute the relative residuals?

Any help is much appreciated.
Best regards.

2

That does not seem true to me.

The code around the line you linked is.

    // Initialize _residual0
    if (_iteration == 1)
    {
      PetscReal _r = 0.0;
      VecNorm(_dx, NORM_2, &_r);
      _residual0 = _r;
    }

It is true that _r is initialized to zero, but the next call to VecNorm will replace the float contained in _r with the norm of _dx.

3

Honestly, I’m very confused about this C++ Code. I guess you’re right, but I’m still questioning if it is reasonable, that the relative residual is non-inf on the 0th iteration of any further time step. Can you comment on that too please?

4

This is because we are re-using the solver, and the residual is not reset at the beginning of a new solve step. Thus the initial residual (_residual0) used to compute r (rel) is not 0.
You can see that the residual is evaluated at:

where the convergence check is called: dolfinx/cpp/dolfinx/nls/NewtonSolver.cpp at main · FEniCS/dolfinx · GitHub
This in turn calls solver.residual0() which is initialized as 0: dolfinx/cpp/dolfinx/nls/NewtonSolver.cpp at main · FEniCS/dolfinx · GitHub
but it is updated after the first iterationL: dolfinx/cpp/dolfinx/nls/NewtonSolver.cpp at main · FEniCS/dolfinx · GitHub

Note that if you find the C++ Newton solver confusing, you can easily write out your own Newton solver where you can choose your own convergence criteria, as shown in:
https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html

5

This seems very counterintuitive. Does a new solve/time step not indicate solving a new problem?

Taking the other answer into account:

it seems that the relative residual is computed as r_\mathrm{rel} = \frac{r_\mathrm{abs}}{\Delta x}, with \Delta x being the increment of the first Newton iteration. Thus, e.g. for the hyperelastic problem, the relative residual would have the unit \frac{N}{m}. Is this physical? What’s the interpretation of the relative residual?

6

Feel free to file a pull request to dolfinx rectifying this. It would simply be adding _residual0 = 0.0; at:

which would then consistently give you inf for iteration 0 using “residual” method.

You would have to ask the original author (@garth?). As i have already stated it is easy to implement your own Newton solver if you are not satisfied with the convergence criteria or logic.

7

@dokken is there a way to “reset” the residual manually?

I’m facing a similar issue where, in a pulsed scenario, after a long period of “nothing happens” (ie. solved in zero iteration), a new pulse starts but the relative residual is already below the tolerance

8

Hi,
just to mention that the pull request New NewtonSolver._residual0 initialisation (residual) by SalzmanA · Pull Request #3611 · FEniCS/dolfinx · GitHub may address some of your concerns.
If by chance you can test it (need to apply the change to the c++ sources and compile dolfinx thought …) with your problem it might be interesting to see if it goes in the right direction.

9

I opened this simple pull request to try and implement a quick fix as suggested by Jorgen