Fix some small Newton related problems

[MFraters]

There are some problems with the derivatives. This pull request intents to fix them, but it is currently a work in progress, as discussed with @bangerth.

[bangerth]

Did you want to include all of these other things into this pull request?

[MFraters]

not all of them, only the add derivatives pull reqest. I can see if I can clean it up a bit

[bangerth]

Hm, then I don't quite understand. Does this patch rely on functionality that is still pending in other PRs? Or why is it so much?

[MFraters]

yes, or part of my testing relies on the added utilities function in the #2030, but the order it is shown in here look quite random. There are the three important commits:

commit 5f3e6db406f80c8c2e406b065fc98899cd4bc081
Author: MFraters <menno.fraters@hotmail.com>
Date:   Mon Feb 19 21:50:46 2018 +0100

    commenting some some stuff out and changing a bit...

commit 5c491b87324e9150bc86a67597df58dade3442fd
Author: MFraters <menno.fraters@hotmail.com>
Date:   Mon Feb 19 19:47:32 2018 +0100

    add testing input file.

commit 474f401b1b7fd8062c8d0a287839be81cb5ba564
Author: MFraters <menno.fraters@hotmail.com>
Date:   Mon Feb 19 19:45:47 2018 +0100

    Add drucker-prager compositions material model.

commit 4528385018398401bde010a99401f8e259c37f9c
Author: MFraters <menno.fraters@hotmail.com>
Date:   Mon Feb 19 19:44:36 2018 +0100

    testing changes.

I can also just push them without the commits from #2030

[MFraters]

@bangerth, I tested turning the n_basis_for_tensor function into just 1's, like you suggested, but it is still different. The values look something like this:

(nonlinear channel flow with drucker_prager_compositions)

composition_viscosities[c] = 1.71197e+21, strain_rate_norm_square = 1.17187e-27, strain_rate = -1.91722e-30 2.42061e-14 2.42061e-14 -3.15544e-30, deviator(strain_rate) = 6.19111e-31 2.42061e-14 2.42061e-14 -6.19111e-31, deviator_strain_rate = -1.91722e-30 2.42061e-14 2.42061e-14 -3.15544e-30, composition_viscosities_derivatives[c] = 2.80083e+18 -3.53623e+34 -3.53623e+34 4.60972e+18
analytical = 2.80083e+18, fd = -1.99101e+27
analytical = 4.60972e+18, fd = -1.99101e+27
analytical = -3.53623e+34, fd = -5.86167e+34

(spiegelman with drucker_prager_compositions)

composition_viscosities[c] = 6.85975e+22, strain_rate_norm_square = 2.68196e-30, strain_rate = -1.12488e-15 -9.96315e-19 -9.96315e-19 1.19021e-15, deviator(strain_rate) = -1.15755e-15 -9.96315e-19 -9.96315e-19 1.15755e-15, deviator_strain_rate = -1.12488e-15 -9.96315e-19 -9.96315e-19 1.19021e-15, composition_viscosities_derivatives[c] = 2.87715e+37 2.54832e+34 2.54832e+34 -3.04426e+37
analytical = 2.87715e+37, fd = 9.03494e+36
analytical = -3.04426e+37, fd = -9.55971e+36
analytical = 2.54832e+34, fd = 1.60034e+34

[MFraters]

Found the problem. Forgot that this term was also present in the Spiegelman paper, and should also present in the derivative. This explains why the analytical derivative didn't converge well, while the finite difference did.

[bangerth]

OK, great. So everything now agrees?

[MFraters]

It only agrees when I use

 0    1/2
1/2    0

And only when some parts of the strain-rate are not zero, because when they are the analytical derivatives in that part of the tensor becomes zero, and the finite differences still has a value other then zero.

[bangerth]

I now think that the correct form for these units tensors is actually

  0            1/sqrt(2)
  1/sqrt(2) 0

That's because this is a tensor that has norm 1.

[bangerth]

Please rebase and let's see how this works then.

[MFraters]

So I managed to fix two problems. The first one has to do with the fact that the derivative of the regularization we computed actually needed the unregularized viscosity to compute the regularized derivative. The second problem was that I never added the analytical pressure derivative. Everthing seems to work fine, but I cheated a bit on the analytical pressure derivative, because I removed a 1/2 to get it working, but I do not know why yet. I will look into this tomorrow.

[MFraters]

ow, I forgot to mention, the analytical and the finite difference derivative are now exactly the same, except for the off diagonal terms, which are very slightly off, but it doesn't seem to hinder convergence:

strain-rate derivative tensor:
analytical = 9.11386e+37, fd = 9.11386e+37
analytical = -9.11377e+37, fd = -9.11377e+37
analytical = -1.08329e+32, fd = -1.11412e+32
pressure derivative
analytical = 6.73883e+13, fd = 6.73882e+13

[bangerth]

On 02/21/2018 12:20 AM, MFraters wrote: |analytical = 9.11386e+37, fd = 9.11386e+37 analytical = -9.11377e+37, fd = -9.11377e+37 analytical = -1.08329e+32, fd = -1.11412e+32|

The difference is about 0.03e32=3e30 out of a total of almost 1e38. This difference is zero.

[MFraters]

@bangerth, to summarize, there are two things we need to decide.

1. The simple nonlinear in here has an extra 2 in front of the strain-rate, but in the paper we don't have that. I think it is good to have them the same. Which one should we choose?
2. The drucker_prager.cc in the material model also contains a regularization (but slightly different from the spiegelman). I could add a pure drucker prager to the benchmark folder for the paper?

[bangerth]

About 1: Do you mean this line?

const double edot_ii = 2.0 * std::max(edot_ii_strict, min_strain_rate[c]);

This looks wrong. Both of the definitions we discuss in the paper on p. 6 use edot_ii = 1/2 eps:eps or the negative of this. This is also how edot_ii_strict is defined in the formula copied above. The factor of 2 seems wrong to me. Where did it come from? Are you trying to reproduce an existing benchmark that has this factor?

About 2: Is drucker_prager.cc used in ASPECT anywhere, or is it a file that we only use for the purposes of the benchmark of this paper? I have to admit that I don't understand the context of your question -- none of the material models we discuss in the paper uses pure DP. I think that for the benchmarks in the paper, we should implement exactly the material model we state in the paper. If ASPECT does not, out of the box, already provide such a material model, then we should implement it separately from ASPECT. Or do I completely misunderstand where you are going with the question?

[MFraters]

about 1: It comes I think from that people often in rheologies a term 2.0 * edot_ii appears somewhere, see for example the drucker prager (line 106) and the visco plastic material model (line 246 and 291). But I am totally fine with removing it here, it would make more sense to me.

about 2: the material model drucker_prager.cc was already there before the paper, and I think it is used in the crustal cookbooks. You are right that non of the benchmarks we show in the paper use pure drucker prager, but we discuss pure ducker prager in appendix B2. That is why I ask whether we should put it somewhere in aspect, so that people can for themself see that alpha is always a half if you implement it in a pure way.

[bangerth]

On 03/01/2018 01:07 PM, MFraters wrote: about 1: It comes I think from that people often in rheologies a term |2.0 * edot_ii| appears somewhere, see for example the drucker prager (line 106) and the visco plastic material model (line 246 and 291). But I am totally fine with removing it here, it would make more sense to me.

I think I'm still not quite clear about how the factor of 2 got there. If the formula reproduces a material model that has been used in a different paper, then we should have the same factors that they have as well. If we are discussing a material model without reference to other works, then the factor of two seems wrong. (I can't look up right now where the lines you reference are -- I assume you are referencing lines in the paper -- because I don't have an internet connection: I'm somewhere over the Aleutians in the middle of the North Pacific.)

about 2: the material model |drucker_prager.cc| was already there before the paper, and I think it is used in the crustal cookbooks. You are right that non of the benchmarks we show in the paper use pure drucker prager, but we discuss pure ducker prager in appendix B2. That is why I ask whether we should put it somewhere in aspect, so that people can for themself see that alpha is always a half if you implement it in a pure way.

I see. Yes, I think that would actually be quite nice. It reminds me that we should get better at providing material models *as building blocks* for other material models. For example, the regularized Drucker Prager model should not reimplement things, but just *build* on a pure DP model. That is something to fix for a later time, but if you were willing to add such a pure DP model, then that would be nice!

[MFraters]

about 1: I talked to @cedrict and we looked it up in the place where the benchmark originally comes from, which is the book of Taras Gerya (equation 16.4). In there the 2 is present, so I think we should keep it.

about 2: Should I add a pure Drucker Prager in this or an other pull request?

[bangerth]

On 03/07/2018 09:42 AM, MFraters wrote: about 1: I talked to @cedrict <https://github.com/cedrict> and we looked it up in the place where the benchmark originally comes from, which is the book of Taras Gerya (equation 16.4). In there the 2 is present, so I think we should keep it.

OK, let's keep it then.

about 2: Should I add a pure Drucker Prager in this or an other pull request?

Separate is always easier.

[MFraters]

I will add a reference to the comments

[MFraters]

To summarize, the wasn't really an spd problem in the end, we just didn't account for the regularization. In the mean time we found three other issues, which we solved here.

1. We now use the exact powerlaw description from the book of Taras Gerya for the simple nonlinear.
2. optimalize the drucker prager derivatives and add the deviator change for the derivative.
3. Align the nontation of the spd computation with the paper.