[WIP] Revert one change of PR #5580 to improve Newton solver convergence

[gassmoeller]

This is a result of #6159. I would like to discuss reverting this change that happened originally in #5580. The difference between the strain_rate and effective_strain_rate is for most models (that dont use elasticity) that effective_strain_rate is simply the deviatoric part of the strain rate. I do not understand why we use the deviatoric strain rate here (it will be used to assemble the Newton matrix, i.e. the Jacobian). Maybe I just didnt follow the discussion in #5580 closely enough. From my understanding this should be the full strain rate, and at least the benchmark in #6159 massively benefits from my change here:

These are the relevant results from the benchmark in #6159 pre #5580:

These are the relevant results post #5580. You can see you the stabilized version / dots converge faster, but the unstabilized version / lines are much worse:

Now these are the relevant results post #5580 but with the change I propose in this PR. The unstabilized models are reverted to their original (better) convergence, and even the stabilized models gain another small boost:

I think this change may be responsible for some of the changes in #6159 that I couldnt explain. I will have to rerun more benchmarks, but in the meantime @YiminJin and @MFraters could you remind me why we did this change originally?

[YiminJin]

@gassmoeller Sorry for the late reply. The reason that I use deviatoric strain rate in the Jacobian matrix in PR #5580 is that there were only few rheologies (ViscoPlastic, CompositeViscoPlastic, ...) providing viscosity derivatives at that time, and I remember that all of them use the deviatoric strain rate when calculating the viscosity. As a result, the strain rate in the Jacobian matrix must be deviatoric in order to maintain the consistency. I have carried out experiments with material model ViscoPlastic and the convergence rate with deviatoric strain rate was much faster that with full strain rate.

I noticed that in the Spiegelman fail test, the parameter Use deviator of strain-rate is set to false, so the material model use the full strain rate in calculating the viscosity. In this case, the assembler should use the full strain rate to assemble the Jacobian matrix.

The modifications I made in PR #5580 was some sort of hard coding. We may have to find some way to tackle this problem.

[gassmoeller]

@YiminJin Thanks for the feedback! I agree with your statement that using the deviatoric strain rate provides much better convergence rate for models that use the stabilization, this is also what I observed in my tests in #6159 and which is the reason I kept using the deviatoric strain rate if the PD stabilization is applied. However, I am still in favor of including my modification in this PR. I have the the following two reasons (tell me if I am wrong here):

• This modification will only use the full strain-rate if no PD stabilization is applied. As I showed in Update spiegelman newton benchmark #6159 and in the figures above, using the deviatoric strain rate destroys the quadratic convergence of the Newton solver for unstabilized models. In practice unstabilized models are not common, because they lead to crashes with most of our rheologies, but if we cannot show the quadratic convergence of the Newton solver at least in simple benchmark models, we cannot really claim to have a correctly implemented Newton solver. As far as I can see your test models in Fixup newton solver and elastic rheology #5580 always activated the stabilization, this is why you didnt notice that the unstabilized models actually became worse (and we have very few tests for not stabilized models, this is why we also didnt notice it in the test suite).

• The first reason is independent of whether the material model uses the deviatoric strain or the full strain rate to compute the viscosity. I have tested both, and in both cases the full strain rate allows for better convergence of the unstabilized models than the deviatoric strain rate. This is a bit unintuitive, because as you mention, most of our material models use the deviatoric strain rate to compute the derivative, and it would seem natural to use the deviatoric strain rate in the Jacobian as well. However, consider the following case (with $\eta$ being viscosity and $\epsilon$ being strain-rate:

  We want to compute the term $\frac{\partial}{\partial\epsilon} (\eta(\epsilon) \epsilon)$, then the product rule is

  $\frac{\partial\eta(\epsilon)}{\partial \epsilon} \epsilon + \eta(\epsilon) \frac{\partial \epsilon}{\epsilon}$. I would expect the $\epsilon$ in the first term to be the full strain-rate, no matter how $\frac{\partial\eta(\epsilon)}{\partial \epsilon}$ is computed inside the material model.

Would you agree with this? Let me know in case I got something wrong, I only really started looking into this in the past few weeks.

[YiminJin]

Thank you for the detailed reply. With the help of your explanation, I realized that I made two mistakes in my previous comment:

1. The strain rate in the system Jacobian should be the full strain rate;
2. The experiments I carried out in my memory was not for the deviatoric strain rate in the system Jacobian, but for the full strain rate in the system residual.

I have tested your modifications with the modified continental extension cookbook in Issue #6046, and found that the convergence rate was much better than the main branch when no stabilization is applied (the main branch takes 101 iterations with a lot of line search steps to converge to 10-7 in the first timestep, while the present PR takes only 52 iterations and no line search step).
Therefore, I agree with your statements. Thank you for pointing out my mistake. I have also taken a look at #6159 and I will think about the reason for 1. why the PD stabilization needs deviatoric strain rate; 2. why the convergence rate with GMG solver is worse that AMG.

[gassmoeller]

Thanks for checking! Glad to know you come to the same conclusion. I will close this PR then, as the change I proposed here was merged as part of #6159 already.

I will think about the reason for 1. why the PD stabilization needs deviatoric strain rate;

Yes that would be good to understand better. When I tried using the full strain rate with the PD stabilization I got warnings that the matrix was no longer positive definite (these checks/warnings are only enabled in debug mode).

2. why the convergence rate with GMG solver is worse that AMG.

Also a good question, but I suspect it is related to material model averaging. One could for example run the Spiegelman benchmark set of #6159 with AMG and harmonic material model averaging and compare that against GMG with harmonic material model averaging. If they more or less agree, then the difference is likely just that GMG requires the averaging, while AMG does not.