Update spiegelman newton benchmark

[gassmoeller]

So this started out as a follow-up to #6135 to update the Spiegelman 2016 benchmark for the Newton solver. It ended up being a lengthy dive into our version history and I think I need outside input on how to put this together (@MFraters can we speak some time?). Before the wall of data here are my conclusions from the work:

• There were significant changes in the solver behavior of the Newton solver, both between the original Fraters et al. 2019 paper and ASPECT 2.5 and between ASPECT 2.5 and 3.0. However, my changes in Move pressure scaling #6135 are not affecting the nonlinear solver very much it was mostly other changes.
• To summarize the changes is hard, because they are not very systematic, but to try:
  • The SPD stabilized version of the solver tends to converge faster now (and is more often identical to the non stabilized version), indicating better stabilization choices.
  • However, the solver seems to have lost the quadratic convergence behavior in some models (it still converges, and faster than the DC or simple Picard, but not quadratically accelerating over iterations). This happened even for non-stabilized models so maybe something changed in the Jacobian?.
• The GMG solver has significant effect on the Newton solver. While the linear solver converges fine, the nonlinear solver often converges worse (though sometimes better) than when using AMG. It is not the GMG solver itself that is the problem, but the material model averaging that is required by the GMG solver. This was already mentioned in Fixup newton solver and elastic rheology #5580 (comment). In consequence: should we disable the combination Newton/GMG? It still converges in many models, just slowly. The faster GMG linear solver may make up for the worse nonlinear convergence behavior.
• The Spiegelman benchmark in ASPECT needs to be rewritten. The existing plot file does not produce Fig. 4 from the paper, but something completely different. The split into metabash.sh and bash.sh is very confusing (I simplified it here). It should be easy to reproduce the figure from the paper to check this benchmark for correctness. Also we should include the standard input.prm benchmark case as a test.

So here is the output data I have, while trying to reproduce Fig. 4 of Fraters et al. 2019:

Original Figure 4 from paper:

My version run with ASPECT 2.5 (AMG). I can not guarantee this uses identical parameters to the paper, because the figure and scripts included in ASPECT where clearly different from the original paper (I tried to reproduce according to the description of the paper). This is already different from the paper, but still close:

This is the most relevant comparison to ASPECT 3.0 (AMG). You can see how the SPD stabilized models (dots) converge better than above, but the non-stabilized versions (lines) are much worse:

Here are the results of some other things I tried:

Results with ASPECT 3.0 (GMG), for comparison with AMG above. Some models converge better, some worse:

Results with ASPECT 2.6.pre (GMG). This version was ASPECT right before #6135 was merged to tease apart the influence of #6135. There are some changes, but in both directions (worse/better). On average I would say the behavior is similar:

As mentioned above, I am not sure what conclusion to draw from all of this. I tried to look for problems/changes in the code that caused the different convergence behavior, but apart from the obvious (different SPD stabilization factor) I have trouble understanding the following:

• Why is the DC Picard convergence in all models now worse than the Fixed Point Picard convergence? This was not the case in the original paper, but was already present in ASPECT 2.5. Maybe the residual is now computed differently than back then? this was a bug in my run script
• What was the change the destroyed the quadratic convergence for (some) of the unstabilized models? It must have happened between 2.5 and 3.0, independent of linear solver (/tolerance), independent of stabilization. Some change in the Jacobian matrix? see discussion below and [WIP] Revert one change of PR #5580 to improve Newton solver convergence #6160
• Some of the ASPECT 3.0 models show a sudden increase in nonlinear residual right when switching from defect-correction solver to Newton solver that was not present in 2.5 or earlier models. I can prevent this increase by allowing for more line-search iterations (=not moving in the direction of increasing residual), but this leads some models to never converge. this must have been introduced in Fixup newton solver and elastic rheology #5580 and is only present in the stabilized models. since the stabilized models are now better than before once they do converge, it is likely an acceptable change

I think it would be great if we could figure out some of these questions.

[gassmoeller]

Ok, @MFraters I think this is ready for a review and I could address most of the things I didnt understand originally.

This PR does 3 things:

1. It reverts a change to the strain rate used in the Jacobian that only affected the convergence rate of the Newton solver if the matrix was not stabilized to be SPD (this is the change originally introduced in Fixup newton solver and elastic rheology #5580 and discussed in [WIP] Revert one change of PR #5580 to improve Newton solver convergence #6160).
2. It reworks the Spiegelman benchmark of the Newton solver to be easier to execute and plot and extends its documentation. The results produced are now also much more similar to the original results of the paper than what I showed above. I will post figures below.
3. It adds two tests that test the benchmark prm file and an unstabilized version of the benchmark prm file. This should help prevent accidental changes to this solver or the benchmark in the future.

Just as a reminder, these are the original results from the paper:

And these are the results I could produce with this PR:

[MFraters]

Thanks @gassmoeller for working to fix the convergence and cleaning up the benchmark code. In hindsight we probably should have split up #5580 in three different pull request, so that all individual changes to the solver could be reviewed and tested separately.

It is nice to see that the both the unstabilized and stabilized versions are now generally faster and have a more stable convergence rate than in the paper, and in some cases can converge, where the paper couldn't (although for the hardest case, it is the other way around, but that was just one sub-case which may just have gotten lucky).

I will wait with merging for a day in case @bangerth or @YiminJin still want to take a look at it, but I think it is good to merge.