Issue with strain values when using plastic strain

[erinheilman]

In some models you can get unrealistic values of strain (in this case plastic strain) when using strain-weakening and strain-healing. There are very negative values and really high positive values as well. It's occurred in some large 3D models as well as simpler cylindrical models. I've included the parameter file for the simpler cylindrical case.

strain_case.txt

[bangerth]

I don't know what might be wrong, but the pictures are cool!

[gassmoeller]

You are right negative strain values are certainly not realistic, but we will need some more details to find the reason for this. @naliboff may have some more ideas, but here is what I would suggest:

• I dont see where the negative values occur in your pictures, can you try to zoom in on the region, and pinpoint the time they first occur? Ideally we want to know the exact conditions under which the values become negative. I.e. does it happen in high-strain rate regions? Does it happen in regions where there is a large gradient in accumulated strain (this could indicate overshooting in the compositional field)? Is it related to negative pressures, or any other property of the solution that can influence the accumulation of plastic strain. If you do not see that in the current output, try rerunning with visual output at every timestep
• Concerning the large values: What do you mean by "unrealistic"? If I remember right at least in principle there is no upper bound for the accumulation of plastic strain. I agree that your box model looks weird because there is a single spike of strain of about 2000, so we should look into why strain can reach very large values in a very small region, but it will be hard to define a perfect criterion for what constitutes too large values.
• For easier debugging I would suggest the following: Try to find the simplest 2d box model in the coarsest resolution that can create the negative values, then it will be simpler to pinpoint the conditions that lead to the values. Try to remove all other complexities one by one from the model (e.g. strain healing, strain weakening) to find out which of them causes the problems.
• In your example parameter file you use very coarse solver tolerances for the nonlinear Stokes solver and for the temperature and advection solver. You probably did this for better performance, but it could also be the reason for the spikes, if Stokes or the strain advection equations are solved too inaccurately. Do the problems also appear with our default tolerances?

[bangerth]

Regarding a limit for the strain rate: From a physical perspective, the stress has to be a continuous function and because the stress is viscosity * strain rate, this means that the strain rate should be a function of the sort continuous function / viscosity. So if the viscosity becomes very small in some region, then the strain rate can become very large there.

Of course, this is exactly what happens in a fault: the viscosity becomes quite small and the strain rate becomes quite large.

Whether that is what happens here is a question to which I have nothing to contribute, though :-)

[bangerth]

That said, about the bottom left picture with the point of high strain rate: Given that strain rate = continuous function / viscosity, and given that you plot the viscosity with a log scale at the bottom right, it would make sense to show the strain rate on a log scale as well. If you do, I bet you will see an inverted pattern from the viscosity -- or at least that's what I would expect if the simulation is correct.

[gassmoeller]

@bangerth: The plots and problem exist in the integrated plastic strain (i.e. the time integral of the part of the strain rate that is accomodated by plastic deformation), not necessarily in the strain rate itself, although that is of course one possibility. All of your comments still apply, but we additionally advect and integrate this property over time, which could also be part of the problem (i.e. overshooting in the advection solution could cause the negative strain). Plotting the strain in log scale is certainly useful to check if the positive spike is really anomolous. Paraview will complain that there are negative values, but I think it will still display the values > 0.

[erinheilman]

@gassmoeller I'm looking into this now, @naliboff agreed that it may be an issue with the long timescales the models are run for and the advection causing over/undershoots as well

[naliboff]

@gassmoeller @bangerth - as mentioned above, @erinheilman, @bobmyhill, @anne-glerum, and I discussed this yesterday.

We all agree a likely source for the large negative strains is overshooting during the advection solution. Small negative strains are almost observed in the models @anne-glerum and I run, but the finite deformations are simply much larger in the models here.

@erinheilman - I agree with @gassmoeller that a good next step is to try a coarser 2D model with much stricter solver tolerances. I typically never adjust the advection solver tolerances, but the Stokes systems I would use a non-linear solver tolerance of 1e-4 or 1-e-5. The single Advection, iterated defect correction Stokes solver scheme should allow fairly faster convergence. See the following for an example of the parameters once could set for this scheme:

set Nonlinear solver scheme                = single Advection, iterated defect correction Stokes
set Nonlinear solver tolerance             = 1e-4
set Max nonlinear iterations               = 100

subsection Solver parameters
  subsection Stokes solver parameters
    set Number of cheap Stokes solver steps = 0
  end
  subsection Newton solver parameters
    set Max pre-Newton nonlinear iterations      = 100
    set Maximum linear Stokes solver tolerance   = 1e-2
    set Nonlinear Newton solver switch tolerance = 1e-4
    set Use Eisenstat Walker method for Picard iterations = true
  end
end

@erinheilman - In addition to the various testing suggested above, I think it would be great to add the strain healing law to the visco_plastic particle properties (larger rework of that code to happen later). We can also try doing some tests with simpler models to see how the strain values compare with compositional fields verses particles.