Modify the plastic dilation terms to reduce noises in effective viscosity

[YiminJin]

This PR reduces noises in the effective viscosity by modifying the plastic dilation terms. It greatly improves the results of Duretz's (2018) test:

I am preparing to board the plane. I will provide a detailed explaination when I arrive at Davis.

[YiminJin]

First, let us see the comparison between the results of the strip footing test and Duretz's (2018) test before and after the modifications in this PR (the computations have also applied the modifications in #6471):

Strip footing:

Duretz (2018, test3):

It can be observed that before the modifications, some areas near the shear bands have a slightly lower viscosity than the unyielding area. These "noises" in the viscosity field affects the velocity field and eventually lead to incorrect strain rate patterns (see the results of Duretz's test).

I have suspected that the "noises" were caused by the modifications in elastic rheology, but after two days' debugging it turns out that the mistake is mine (apologies to Anne). In #6373, the plastic dilation is calculated by the following steps:

1. For each composition, check if the material is yielding;
2. If the material is yielding, then compute the dilation terms; otherwise, set the dilation terms to 0;
3. At each quadrature point, calculate the arithmetic average of all compositional fields.

This implementation will cause a problem: assume that there are two compositional fields with the same mechanical parameters at one quadrature point. Further, assume that the volume fractions of the two compositional fields are 1/2 and 1/2, and one is yielding, while the other is not, then the left-hand side dilation term will be calculated by $\frac{1}{2}(\bar\alpha\alpha/\eta^{ve} + 0) = \frac{1}{2}\bar\alpha\alpha/\eta^{ve}$, which could be regarded as a "new" material with smaller dilation angle and friction angle. This is incorrect. In this PR I changed the implementation for dilation terms, and now we can get the correct result for Duretz's test. In the strip footing test the "noises" in the viscosity field still exist, but it is the best I can get for now.

[YiminJin]

The broken test is the one I added in #6373 for plastic dilation. On a refined grid, the results before and after the modification are:

The result is actually improved after modification. So I applied the changes directly.

[gassmoeller]

Hi Yimin this is pretty cool, thank you!

I have one comment on the code (just a request for documentation), and one question: I completely agree the new results look less noisy and more localized. However, it looks like the Duretz picture you posted has even more pronounced and spatially limited shear bands than the picture in Duretz et al 2018 (I am looking at figure 5e for example). So my question is: Do we even expect such a localized shear band? I.e. is our result now "better" than the one in Duretz or just "different" and we do not know which one is the correct result? I am ok merging this PR if you are convinced this result is better, but I wanted to ask what you think about this.

[gassmoeller]

Could you leave a comment here explaining why you need this if-condition and what the individual if-else branches compute?

[YiminJin]

I have added comments here. Please let me know if I need to expand them in more details.

[YiminJin]

However, it looks like the Duretz picture you posted has even more pronounced and spatially limited shear bands than the picture in Duretz et al 2018 (I am looking at figure 5e for example). So my question is: Do we even expect such a localized shear band? I.e. is our result now "better" than the one in Duretz or just "different" and we do not know which one is the correct result?

I am not sure about this question. I remember that I have obtained sharper shear bands with elASPECT too, but I do not have that figure at hand. The following is a similar test obtained by pseudo-transient method (https://arxiv.org/abs/2305.01701):

In subfigure (c) and (g) the patterns are not even symmetric. I think this reflects the ill-posedness of elasto-(visco)-plastic problems without regularization, and a small difference in model settings (like the viscosity averaging, or mesh resolution) may lead to large differences in the result.

The sharp shear bands do have a problem: the convergence rates of both the linear solver and the linear solver are low. When applying cellwise harmonic averaging to viscosity, the strain rate field looks like

the computing time is halved, but the widths of shear bands are almost doubled. I am trying to implement "project to Q1" for elastic rheology to see if it helps.

[YiminJin]

I forgot to say: the reason that the strain field in each "block" cut by the shear bands is nearly constant in my result is perhaps because the material is incompressible (the tests in Duretz et al. (2018) are compressible).

[YiminJin]

I give up implementing project_to_Q1 for elastic rheology: it is not suitable for plastic models. Moreover, I find that cellwise averaging affects the strain rate pattern greatly in my tests, while the computational cost is not reduced sufficiently. I will try to accelerate the computation through other means.

@gassmoeller The result of test1 (which is very similar to test3) produced by elASPECT looks like:

The shear bands are sharper than not only the results in Duretz's paper (Figure 5), but also the results of ASPECT. I think it is normal to get different results from different codes in this sort of problems.

[gassmoeller]

The shear bands are sharper than not only the results in Duretz's paper (Figure 5), but also the results of ASPECT. I think it is normal to get different results from different codes in this sort of problems.

Thanks, I agree that it is normal to see differences, I was mostly wondering if we have any way to say if the differences are reasonable.

I forgot to say: the reason that the strain field in each "block" cut by the shear bands is nearly constant in my result is perhaps because the material is incompressible (the tests in Duretz et al. (2018) are compressible).

Ah, that is good to know, this explains it well.

I will restart the testers and if they are ok we can go forward with this PR.

[gassmoeller]

/rebuild

[gassmoeller]

looks like the testers are fine. lets merge.