New NewtonSolver._residual0 initialisation (residual) #3611
Conversation
* Forces _residual0 to be set to the norm of the first residual vector when convergence_criterion=="residual". * Done after the first convergence test, so that this test can play a role in the case of multiple calls to the solve method. * Retain the same behaviour when convergence_criterion=="incremental".
|
To me, this makes sense. @garth-wells @jhale what do you think? |
jorgensd
added
enhancement
|
@jorgensd I'm afraid I'm not able to test this anytime soon. However in terms of functionality, I guess there was a reason to move away from how the residual was calculated in legacy-fenics? If that's the case, would it be better to expose this behaviour to the users and let them choose between the old behaviour (what is implemented in this PR) and the one currently available in dolfinx? |
| @@ -232,6 +228,9 @@ std::pair<int, bool> nls::petsc::NewtonSolver::solve(Vec x) | |||
| // set. | |||
| if (_iteration == 1) | |||
| { | |||
| PetscReal _r = 0.0; | |||
| @@ -232,6 +228,9 @@ std::pair<int, bool> nls::petsc::NewtonSolver::solve(Vec x) | |||
| // set. | |||
| if (_iteration == 1) | |||
| { | |||
| PetscReal _r = 0.0; | |||
| VecNorm(_dx, NORM_2, &_r); | |||
| _residual0 = _r; | |||
| _residual = 1.0; | |||
There was a problem hiding this comment.
Unrelated to the PR, but I think "residual" initialization _residual = 1.0 for the incremental criterion should happen above, outside of the loop just as it does for residual criterion. In the current code, even if the first ||dx|| < atol I would have to wait one additional iteration.
Besides the naming is not great and what we call "residual" for incremental criterion is actually increment norm...
There was a problem hiding this comment.
I think this piece of code could be completely refactored so that the convergence criterion can be switched out based on passing a different std::function?
There was a problem hiding this comment.
It can? You can do set_convergence_check?
https://github.com/search?q=repo%3AFEniCS%2Fdolfinx%20set_convergence_check&type=code
If you mean that we should have a function that includes the modification of the tolerances, I'm not sure that is worth-while.
Hello,
just in case ... It's almost my first PR, so I hope I'm not too far away from how these things should be done.
The origins:
In making comparisons with another library, I saw that FENICSX needed more iterations to solve the same non-linear test case in order to converge to the same atol/rtol.
And it wasn't possible (due to non-convergence) to reduce the rtol as much as in the other library.
It was the same with either the Python API or the C++ API.
After investigation, both libraries use the same test (with convergence_criterion==‘residual’), but in FENICSX _residual0 is set to the norm of the first variation of the solution, whereas in the other library it is set to the norm of the initial residual vector.
In FENICSX it is not adimentional, which explains why in my case, with my units and settings, I got |dx| far less than |r0| and the relative convergence criterion was really not the same as in the other code.
In https://fenicsproject.discourse.group/t/finite-relative-residuals-for-0th-newton-iteration-in-time-dependent-problem/15214 they seem to have made a similar observation about the _residual0 setting, plus some other questions about using the solve method in a time-dependent context.
How:
This PR now forces _residual0 to be set to |r0| when convergence_criterion=="residual". If convergence_criterion=="incremental", this PR does not change anything.
Now regarding the multiple call to solve (time-dependent context), the assignment of _residual0 is placed just after the test for:
Tests:
For my test case (a crude elasto-damaged asymmetric tension/compression constitutive law), I run the solve method twice, with a 'small' or 'large' variation in the volume load imposed between the two runs.
Only relative convergence is observed (i.e. atol is set to a very small value).
The results (see asym.txt) show that with the PR, the first resolution uses fewer iterations (in agreement with the other library) and the second iterates little or converges immediately.
With the original version of the code, for this specific test, |dx| is well below |r0|, so there are convergence problems.
This PR does not change the cpp/demo/hyperelastic behaviour.
It pass the python unit test.
With the test case demo_cahn-hilliard.py, the analysis of the system resolution at the first time step is shown in ch.pdf.
So with PR and convergence_criterion=="residual" there is one more iteration compared to the original code version (but less than with convergence_criterion=="incremental").
In this case |dx|>>|r0| and the PR look more conservative.
Counting all non-linear iterations of all 50 timesteps, we have 266 iterations for convergence_criterion=="incremental", 189 iterations for convergence_criterion=="residual" with original code and 205 iterations with PR.
With https://jsdokken.com/dolfinx-tutorial/chapter2/hyperelasticity.html, using solver.convergence_criterion = "residual", the PR gives a slightly modified number of iterations (61) over the 10 time steps compared to the original code (59).
In any case (PR or original) we have 69 iterations with solver.convergence_criterion = "incremental".
Some runs are given in hyperelastic_tuto.zip and in particular "ir" runs where an iregular load has been tuned to have a variation at a time step small enough that the "residual" relative convergence criterion is less than rtol with the original code and conservatively greater with PR, with an absolute convergence criterion greater than atol.
Future:
It is certainly possible to find cases where this initial test with the previous _residual0 value disturbs the simulation.
And it is still possible to assign _residual0 before the test, but as this requires an extra norm calculation and leads to at least one iteration, it may be at the user's discretion.
This means an API change that goes far beyond the original intent of this PR, which was to address the one solve case.