issue in morency_doin.cc when calculate the second invariant of strain rate

[QwZhang]

@bobmyhill I recently notice something wrong in aspect/source/material_model/morency_doin.cc when it deals with the second invariant of strain rate: the final value of strain rate should be taken the square root if one have used deal.ii's function second_invariant(), in which this operation is not included internally. In short, it should be somthing like this:

SymmetricTensor<2,dim> strain_rate = in.strain_rate[i];
double second_invariant_of_strain_rate = std::sqrt(second_invariant(strain_rate));

not:

SymmetricTensor<2,dim> strain_rate = in.strain_rate[i];
double second_invariant_of_strain_rate = second_invariant(strain_rate);

Or do I miss something that I don't knew?

[bobmyhill]

Dear Qingwen Zhang,

I have no experience with the Morency Doin model, so I'll leave that discussion to someone more qualified. There is a pull request waiting that addresses your issue for the diffusion dislocation model that should be ready and merged very soon.

Best wishes, Bob On 9 Jun 2015 07:16, "Qingwen Zhang" notifications@github.com wrote:

@bobmyhill https://github.com/bobmyhill I recently notice something wrong in aspect/source/material_model/morency_doin.cc when it deals with the second invariant of strain rate: the final value of strain rate should be taken the square root if one have used deal.ii's function second_invariant(), in which this operation is not included internally. In short, it should be somthing like this: ... SymmetricTensor<2,dim> strain_rate = in.strain_rate[i]; double second_invariant_of_strain_rate = std::sqrt(second_invariant(strain_rate)); ... not: ... SymmetricTensor<2,dim> strain_rate = in.strain_rate[i]; double second_invariant_of_strain_rate = second_invariant(strain_rate); ... Or do I omit something that I don't knew?

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/aspect/issues/527.

[jperryhouts]

Yep, this still needs more testing, but it looks like the paper actually uses even a different definition for their second invariant (maybe incorrect, I'm not sure). The deal.ii second_invariant function calcuates something like: const SymmetricTensor<2,dim> e = in.strain_rate[i]; const double e2inv = second_invariant(e); // = e[0][0]*e[1][1] - e[0][1]*e[0][1];

But what they call second invariant in the paper is defined as: const double e2inv = std::sqrt((e[0][0]*e[0][0] + e[0][1]*e[0][1])/2.0);

This will be fixed soon. Sorry for the inconvenience.

[QwZhang]

Hi all, Thank you for your reply. I have one thing for sure about which method is correct. I run a plastic material to demonstrate this. If I use:

SymmetricTensor<2,dim> strain_rate = in.strain_rate[i];
const SymmetricTensor<2,dim> strain_rate_dev = deviator(strain_rate);
double second_invariant_of_strain_rate_dev = std::sqrt(0.5) * strain_rate_dev.norm(); //CORRECT

I get result (the screenshot when timestep=5) like this: or:

double second_invariant_of_strain_rate_dev = std::sqrt(second_invariant(strain_rate_dev));//CORRECT

which yields: These two results are exactly the same and more reasonable for plastic deformation. However, if use:

double second_invariant_of_strain_rate_dev = second_invariant(strain_rate_dev);

then I get wired result like this:

Best, Qingwen

[bobmyhill]

Hi again,

Thanks for sending us your test results! Jonathan and I will fix these two material models. That being said, ASPECT is a community code, so if you find anything else you think is wrong, or if you implement something new, we would be very happy for you to submit a pull request yourself! If you need any help in doing this, we can certainly give you some tips.

Best wishes, Bob On 9 Jun 2015 23:13, "Qingwen Zhang" notifications@github.com wrote:

Hi all, Thank you for your reply. I have one thing for sure about which method is correct. I run a plastic material to demonstrate this. If I use SymmetricTensor<2,dim> strain_rate = in.strain_rate[i]; const SymmetricTensor<2,dim> strain_rate_dev = deviator(strain_rate); double second_invariant_of_strain_rate_dev = std::sqrt(0.5) * strain_rate_dev.norm();//CORRECT Here is the screenshot when timestep=5: [image: fig 1_sqrt 0 5 strain_rate_dev norm] https://cloud.githubusercontent.com/assets/12026068/8074155/1d52215e-0f61-11e5-8377-fbaa4993349b.png or: double second_invariant_of_strain_rate_dev = std::sqrt(second_invariant(strain_rate_dev));//CORRECT [image: fig 2_sqrt abs second_invariant strain_rate_dev] https://cloud.githubusercontent.com/assets/12026068/8074158/274b2340-0f61-11e5-8bae-ca8f5d4d366c.png

These two results are exactly the same and more reasonable. If use double second_invariant_of_strain_rate_dev = second_invariant(strain_rate_dev); I get wired result like this: [image: fig 3_second_invariant strain_rate_dev] https://cloud.githubusercontent.com/assets/12026068/8074166/44ce1526-0f61-11e5-871f-64b62b1624cb.png

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/aspect/issues/527#issuecomment-110570984.

[bangerth]

All, just as a point of reference: it surprises me that anyone would define the second invariant involving the square root. The first, second, and third invariants of a tensor are commonly defined so that they are linear, quadratic, and cubic functions of the tensor elements. See, for example, http://en.wikipedia.org/wiki/Cauchy_stress_tensor#Principal_stresses_and_stress_invariants and http://en.wikipedia.org/wiki/Invariants_of_tensors .

If Morency and Doin really use a definition that involves the square root, can you please leave a comment in the code making it clear that their definition deviates from the common one, to avoid confusion by readers of the code at a later time?

[QwZhang]

Yes, this definition and implementation per se is correct. I will recheck where the problem is as soon as possible.

[jperryhouts]

I had been sort of focused on this problem with e[0][0]*e[1][1] - e[0][1]^2 vs e[0][0]^2 + e[0][1]^2 but the square root thing is weirder I guess. @QwZhang do you know where this definition comes from?

And yes, once this is all sorted out I will document this awkward definition in the comments and in the manual cookbook (which I am also working on improving). After this week I will finally have time to follow through with those fixes.

[jperryhouts]

For the record though, I am still kind of in favor of removing this material model from Aspect because I don't use it, and can't really vouch for its usefulness in most models.

I think that the cookbook is instructive, but the material model should maybe only be distributed as part of the cookbook, and not one of the built-in material models because it is so specific to this one paper's model setup and I'm not convinced that it's appropriate for general purpose use.

[QwZhang]

I think I find the problems. @anne-glerum First, if one want to get the second invariant of a tensor like (deviatoric) stress or (deviatoric) stress rate, etc., the function second_invariant() is a better choice for its being concise, elegant and more importantly, intuitive; Secondly, for plasticity (yielding), we generally get the effective viscosity by

effective_viscosity = strength  / (2 * sqrt(second_invariant_of_deviatoric_stress_rate) )

So, for the following lines of implementation of plastic rheology:

const SymmetricTensor<2,dim> strain_rate_dev = deviator(strain_rate);
double strain_rate_dev_inv = std::sqrt(0.5) * strain_rate_dev.norm();
strength = (std::max(pressure,0.0) * std::sin(phi) + cohesion * std::cos(phi));
const double viscosity = strength / (2.0 * strain_rate_dev_inv);

the final result is correct but this strain_rate_dev_inv here actually equals to the square root of second invariant of deviatoric stress rate tensor, not the second invariant of deviatoric stress rate tensor as indicated by its name. Therefore, it's not very intuitive and can be replaced by:

const SymmetricTensor<2,dim> strain_rate_dev = deviator(strain_rate);
const double strain_rate_dev_inv = second_invariant(strain_rate_dev);
strength = (std::max(pressure,0.0) * std::sin(phi) + cohesion * std::cos(phi));
const double viscosity = strength / (2.0 * std::sqrt(strain_rate_dev_inv));

This actually yield same results as shown in the figures above.

For viscocity, the effective viscosity can be calculated in various way, in which both (1) the square root of second invariant of (deviatoric) stress rate tensor or (2) the second invariant of (deviatoric) stress rate tensor can be involved to get the effective viscosity. So, it depends on the practical situation.

@jperryhouts, perhaps I have been misleading you, sorry for the inconvenience!

Regards, Qingwen

[bangerth]

On 06/11/2015 11:24 PM, Jonathan Perry-Houts wrote:

I had been sort of focused on this problem with e[0][0]*e[1][1] - e[0][1]^2 vs e[0][0]^2 - e[0][1]^2, but maybe this is okay so long as div(u) = 0? Or would there be a sign issue??

If div(u)=0, then e[0][0]=-e[1][1], so indeed the sign appears wrong.

And yes, once this is all sorted out I will document this awkward definition in the comments and in the manual cookbook (which I am also working on improving). After this week I will finally have time to follow through with those fixes.

Thanks.

[bangerth]

On 06/11/2015 11:30 PM, Jonathan Perry-Houts wrote:

For the record though, I am still kind of in favor of removing this material model from Aspect because I don't use it, and can't really vouch for its usefulness in most models.

I think that the cookbook is instructive, but the material model should maybe only be distributed as part of the cookbook, and not one of the built-in material models because it is so specific to this one setup and I'm not convinced that it's appropriate for general purpose use.

I like the idea of having many material models. Nobody is forced to use them, and we can clearly state limitations and problems in the documentation of every model. I generally also favor being backward compatible, which would suggest keeping models that have been distributed in previous releases. But I'm also not an expert in this particular model, so will defer to your judgement.