aspect
aspect copied to clipboard
geodynamics
tolerance for the "exact" implementation of the Peierls creep
Hi Bob, @bobmyhill , thanks for your reply. I closed the previous issue since I found a more practical one to ask. So I am testing the "exact" implementation of the Peierls with the set of parameters taken from Idrissi. I know Daniel Douglas, @danieldouglas92 has tested it in his model, so it's feasible. But I am stuck with adjusting the "Peierls strain rate residual tolerance" parameter, where the "exact" implementation doesn't converge.
I am starting from the 0th step where I set the reference strain rate to 1e-15. These are the parameters in the ViscoPlastic material model. I started from the default tolerance of 1e-22 and increases to 1e-17, but the inconverge persists:
set Reference temperature = 273
set Maximum viscosity = 1e25
set Densities = 3300.0
set Thermal diffusivities = 1.0e-6
set Thermal expansivities = 0.0
set Heat capacities = 1250.0
set Reference strain rate = 1e-15
set Grain size = 1e-02
# rheology
set Viscosity averaging scheme = harmonic
set Include Peierls creep = true # Change to true for restart
set Viscous flow law = composite # change to composite for restart
set Reference viscosity = 1e20
set Minimum viscosity = 1e18 # change to 1e18 on restart
# diffusion creep
set Prefactors for diffusion creep = 0.5e-31
set Grain size exponents for diffusion creep = 0.0
set Activation energies for diffusion creep = 0.0
set Activation volumes for diffusion creep = 0.0
# dislocation creep
set Prefactors for dislocation creep = 0.5e-31
set Stress exponents for dislocation creep = 1.0
set Activation energies for dislocation creep = 0.0
set Activation volumes for dislocation creep = 0.0
# Parameters for Mei Low Temperature Plasticity
set Peierls creep flow law = exact
set Peierls strain rate residual tolerance = 1e-17
set Peierls glide parameters p = 0.5
set Peierls glide parameters q = 2.0
set Stress exponents for Peierls creep = 0.0
set Peierls stresses = 3.8e9
set Prefactors for Peierls creep = 1e6
set Activation energies for Peierls creep = 566e3
set Activation volumes for Peierls creep = 0.0
Here is the error
Exception 'ExcMessage("No convergence has been reached in the loop that determines " "the Peierls creep viscosity. Aborting! " "Residual is " + Utilities::to_string(strain_rate_residual) + " after " + Utilities::to_string(stress_iteration) + " iterations. " "You can increase the number of iterations by adapting the " "parameter 'Maximum Peierls strain rate iterations'.")' on rank 0 on processing:
An error occurred in line <230> of file </home/lochy/Softwares/aspect/source/material_model/rheology/peierls_creep.cc> in function
double aspect::MaterialModel::Rheology::PeierlsCreep<dim>::compute_exact_viscosity(double, double, double, unsigned int, const std::vector<double, std::allocator<double> >&, const std::vector<unsigned int>&) const [with int dim = 2]
The violated condition was:
!abort_newton_iteration
Additional information:
No convergence has been reached in the loop that determines the
Peierls creep viscosity. Aborting! Residual is -nan after 9
iterations. You can increase the number of iterations by adapting the
parameter 'Maximum Peierls strain rate iterations'.
--------------------------------------------------------
Aborting!
So, any instruction on how to set the tolerance?
@lhy11009, the Idrissi flow law (https://www.science.org/doi/full/10.1126/sciadv.1501671) is unphysical. You can see this by plugging zero stress into the flow law, which results in a non-zero strain rate (i.e. the viscosity approaches zero at small strain rates).
To deal with this, @danieldouglas92 implemented a stress cutoff in the Peierls flow law (https://github.com/geodynamics/aspect/pull/4010), which can be added to your prm file using the "Cutoff stresses for Peierls creep" parameter. You can see an example of the use of this parameter in the test file tests/visco_plastic_peierls_cutoff_test.prm.
Let us know if you encounter any other problems.
Hi @bobmyhill. This solves the issue but also brings another question: what is a good choice for this "Cutoff stress"? I am using 1Mpa in my prm. Let me know what you think.
@danieldouglas92 - Can you comment on what values you used for the Cutoff stress in the Hawaii flexure models?
Hey @lhy11009, in my models I used a value of 1000 Pa for the cutoff stress. I set it at this value without doing too much testing, but definitely the cutoff stress should be set lower rather than higher since the strain rate is approximated by a quadratic for stresses smaller than the cutoff stress.
My thought would be that as long as your cutoff stress is chosen to be sufficiently smaller than the maximum stresses in your models that whatever deviation exists between the approximate and exact solutions at these small stresses wouldn't impact results too much.
Hi @lhy11009, Thank you @naliboff and @danieldouglas92 for answering. @lhy11009, one thing you may wish to do is plot viscosity vs strain-rate for different cut-off stresses. That will give you an idea of sensible choices to use - you want to find a cutoff stress that is small, but is large enough to ensure that (d viscosity / d strain rate) is always negative (and that the viscosity is much much higher than the dominant creep mechanism at low strain rates, usually diffusion creep).
Thanks for your answers, @danieldouglas92 @bobmyhill @naliboff . I'll do more tests on that before using it. It makes sense to me that there is a trade-off between having a negative value of the derivative (if I get this right, that's because we are solving with a newton solver) and reproducing the flow law.
The negative d viscosity / d strain rate at low strain rate is a property of any well-behaved power law flow law (nothing to do with Newton solvers). See the figure below. If a flow law (like the unmodified Idrissi law) has a positive d viscosity / d strain rate at low strain rates, then at very low strain rates that flow law would be "more active" than diffusion creep (which has d viscosity / d strain rate = 0). That is (a) physically unreasonable and (b) leads to poor convergence behaviour in numerical models.
The cut-off stress implemented in ASPECT is just one way to ensure the Idrissi flow law does not become active at low stresses. We could also have written a few lines to turn off the flow law at low stress. The benefit of our "cutoff stress" implementation is that eta(epsilon) remains continuous, differentiable and cheap to evaluate.
(taken from https://www.science.org/doi/10.1126/sciadv.abm1821)
Hi Guys
I have put up this test for the Iddrissi flow law using the parameters in the test tests/visco_plastic_peierls_cutoff_test.prm
I expand this to a diagram by setting a linear variation of temperature in a box geometry: (x= 0, T = 273 K; x = 800 km, T= 2000 K) and a reference viscosity of 1e-15.
Other creep viscosities are set to a very big value so that only the Peierls viscosity is output.
Here is the test prm file:
original.prm.txt
Here is the result of the viscosity (in log value and with two contours of T at 1200 and 1400 K):
The issue I found with this test is the low viscosity in the hotter region. In this case, as the strain rate is set to 1e-15, it's clear that the stress value in this region is very small (the smallest values around are 1e-4 Pa).
- Could you assist me further with this issue by setting a better value for the cutoff stress?
- By reading the implementation in the function PeierlsCreep
::compute_exact_strain_rate_and_derivative, the cutoff stress is not really cutoff in terms of the value of the stress, but is modified with another expression. So the idea is, the negative derivatives are cut off and the stress could drop to a smaller value, is this correct?
While debugging issue myself, I played with the value of Cutoff stresses for Peierls creep and changed it from 1e3 (left figure) to 1e4 (middle), and 1e5 (right). The changes lie in the very hot region, bigger values of viscosity could be computed with higher values of cutoff stresses. Meanwhile, the values around T = 1400 K are still very small. Note that the lower limits of these plots are 1e11 (red color).
So the takeaway is, with this set of parameters, the values from the Peierls creep are weaker than 1e18 when T > 1400, which would lead to a very weak mantle material. @danieldouglas92 , have you observed something like this in your study? Or could it be that you are mainly focused on lithosphere deformation so the warmer temperature is not included in your model? Cause I remember you did a model with the bending curvature of the lithosphere.
@lhy11009 please see the PR description here: https://github.com/geodynamics/aspect/pull/4010/files Your tests agree with this description - if you turn the cutoff off, you can get negative viscosities, and higher values of the cutoff increases the viscosities when stresses are lower than the cutoff.
Could you assist me further with this issue by setting a better value for the cutoff stress?
No, because the "ideal" cutoff should reproduce real behaviour, but that "ideal" value will depend on the poor behaviour of the chosen Peierls flow law. Essentially you are recalibrating the low stress part of the flow law.
By reading the implementation in the function PeierlsCreep::compute_exact_strain_rate_and_derivative, the cutoff stress is not really cutoff in terms of the value of the stress, but is modified with another expression. So the idea is, the negative derivatives are cut off and the stress could drop to a smaller value, is this correct?
Yes, that is correct. The cutoff stress is the stress at which the implemented Peierls law gets turned off and replaced (at lower stresses) with a very simple flow law.
Another note: we did think a little bit about whether in composite theologies the low stress cutoff Peierls law might become dominant at low stress (this would be unwanted behaviour). I don't remember what our conclusion was, but if that is ever a problem for you, either (a) let us know and we'll think of a workaround or (b, preferred) refit the dodgy flow law with a more reasonable stress dependence.
Hi @bobmyhill
Thanks for your answers. I read the highlighted part of the code again and I am more sure that I understand what it does.
We did think a little bit about whether in composite theologies the low-stress cutoff Peierls law might become dominant at low stress (this would be unwanted behavior)
This is basically the issue I am pointing to. In my case, this is triggered from a medium-high T ( >1400 k). My current takeaway is , this is what I want to work around, as the values at these conditions (T >1400k, strain rate = 1e-15) are far smaller than the diffusion creep.
(a) let us know and we'll think of a workaround
I have already considered a workaround. My current workaround is to make the cutoff stress a real cutoff on the stress computed. So basically, I am just removing those very small values in the diagram I showed by limiting the minimum stress the Peierls creep could get to. I'll continue the conversation once I have more results to post.
The flat cutoff is probably a decent, simple workaround. I was considering a power law decay, but that's probably not necessary (at least for composite flow laws).
Hi @bobmyhill , first, thanks for your previous advice. After another look at the code, I have a question about implementation of the function
PeierlsCreep<dim>::compute_exact_strain_rate_and_derivative (const double stress, const double pressure, const double temperature, const PeierlsCreepParameters creep_parameters) const
In this line here, is there a reason that the b term ((E + PV)/ RT) is not in the "arrhenius_cutoff" term?
const double arrhenius_cutoff = std::exp(-d_cutoff);
(https://github.com/geodynamics/aspect/blob/00ee3d0567afbd84f90e5dd20a7f89c7798ca1ba/source/material_model/rheology/peierls_creep.cc#LL339C16-L339C16)
While a few lines after this, the "arrhenius" term actually contains the b term:
std::exp(-b*d_cutoff)
@lhy11009 The "Arrhenius cutoff" is calculated at a reference temperature (defined to be equal to E/R) and pressure (=0). Thus, b_cutoff = 1.
The strain rate and derivative are calculated using this reference "Arrhenius cutoff". We derived the strain rate and derivative with the intention that they should both be continuous with respect to stress at the cutoff, regardless of the pressure or temperature. If you think there is a mistake in this derivation, let us know.
The "Arrhenius cutoff" is calculated at a reference temperature (defined to be equal to E/R) and pressure (=0). Got it. That would answer my question for now. But the value is E / T = 5.66e5 / 8.314 = 68077.9408227 K? Okay, this is more detail than I need, and it's just an early value into the implementation, but I am just double-checking that I understand how that work.
We derived the strain rate and derivative with the intention that they should both be continuous with respect to stress at the cutoff Yes, I can see this point in the implementation.
@lhy11009 yes, E/R is a silly temperature. It is a meaningless normalisation factor, but maybe E/(1000*R) would have been more reasonable.
Hi all Another curious situation is if we have phase transitions (e.g. upper and lower mantle) and we want to turn down the Peierls flow law for some of the phases (e.g. for the lower mantle). Do you have any suggestions? I currently only have a walk-away solution.
@lhy11009 do you mean that the ability to assign different values for each phase has not been implemented for the Peierls flow law? If so, the solution is to implement that functionality (using another flow law as a basis). It may be one of those features that no one has yet needed.
Do you mean that the ability to assign different values for each phase has not been implemented for the Peierls flow law
Yes, I don't remember if we have included that ability for the Peierls flow in the master branch. But with my own branch, I have been using different values with the approximate form of the Peierls creep.
so this is a little bit in detail. What I am used to is assigning a very small prefactor to the phases where I wanted to turn off the rheology:
set Prefactors for dislocation creep = background:5.9076e-16|5.9076e-16|5.9076e-16|5.9076e-16|5.0000e-32|5.0000e-32|5.0000e-32|5.0000e-32
Here, the prefactors are set to 5e-32 in the lower mantle and the purpose is to turn off dislocation creep for the lower mantle (as the viscosity calculated will be very big, so it's turned off effectively).
Initially, I applied this trick for the exact form of the Idrissi flow law (setting Prefactors for Peierls creep instead). This doesn't work.
Again this is in detail, when looking into the code, I understood that the prefactor is in the computation of the gradient, and it cannot be such a small value (the real reason being the range of the float?). I observed the issue by making the "compute_exact_viscosity" function outputs more debug information (P, T, strain_rate, gradient...), and I saw that "stress is 1.915469625816906e+30, Derivative is -0" in the following stderr outputs.
My current walkaround in my model is to add another variable to the Material model to turn off the Peierls creep law.
It may be one of those features that no one has yet needed
I agree, but I expect this is needed when one includes the Peierls creep and also includes a lower mantle (high viscosity).
Also, this might work better if we can talk through it at a user meeting.
`--------------------------------------------------------
An error occurred in line <304> of file </home/lochy/Softwares/aspect/source/material_model/rheology/peierls_creep.cc> in function
double aspect::MaterialModel::Rheology::PeierlsCreep
@lhy11009 Thanks for the info; this makes sense. And you've nicely shown that (a) Peierls is indeed implemented to work with phases and (b) the problem is almost certainly due to dividing through by a very large stress in the depsilon/dsigma derivative function.
One solution is not to be so extreme with your choice of Peierls prefactor and/or vary other parameters so that the Peierls component of the strain rate is negligible over the stress-T range in the lower mantle, but doesn't produce zeros in the derivative function.
Happy to talk through your options at a user meeting. I think other solutions would involve making changes to the code.
@bobmyhill I have another thought into this issue: using the gradient of $\frac{dln(\dot{\epsilon})}{dln(\sigma)}$ instead of $\frac{d\dot{\epsilon}}{d\sigma}$. The idea is the logarithm values are in a safer range than the values themself (e.g. -15 vs 1e-15). Also, there is no A in the expression of $\frac{d\dot{\epsilon}}{d\sigma}$, so there won't be an issue if A is too small (at least not from the gradient and iterations).
I have done 2 tests with Python where I used both implementations of gradient and set the same tolerance (relative residue, 1e-8). In both tests, the iteration using logarithm values converge faster (iterations 8 vs 4, and iterations 11 vs 4):
(Among the two, function peierls_visc_from_edot_newton is implemented with the logarithm value, n is the number of iterations)
Let me know what you think
@lhy11009 this is a great idea. Something similar was used for iterating on grain size (i.e. iterating on ln(d) rather than d). I'm in favour of doing this.
@bobmyhill , Yeah, it makes to iterate the grain size that way. I'll add this change to my Hackathon todo list. My Idrissi flow law has been working well. I feel comfortable concluding this issue. I learned many things from all you guys (@bobmyhill, @naliboff @danieldouglas92 @Rombur), thanks for your help:
- The cutoff stress is important for the Idrissi flow lay, and it might work better with a strict cutoff scheme (use the cutoff stress as the minimum stress)
- The current implementation raises issues with very small prefactors (e.g. 1e-31). The reason is the limit of float numbers during iteration triggered from A in the gradient term.
- To account for the mantle PZs (turn off then peierls creep for perovskite), the current implementation could be improved.
And my todo list includes:
-
Adding an option to use the cutoff stress directly as the minimum stress value. This will solve issue number 1.
-
Use the log value in iteration, this will address 2 directly (A will be removed from the gradient) and will also address 3 with a small trick (Assign a vary small prefactor to phases one want to turn off the Peierls creep on).
I'll leave this issue open till the hackathon if people want to add their thoughts. And we'll discuss more at the Hackathon.
This issue is addressed in PR #5227 and #5229. THe iterative scheme we put forward in conversation (iterate log values) have better performance on the convergence. Moreover, The strict limiting stress would ensure a reasonable value from the Idrissi 16 flow law.