Make Newton iterations in the CompositeViscoPlastic rheology more efficient

[bobmyhill]

This PR does a few small things that facilitate efficient calculation of CompositeViscoPlastic viscosities:

• The approximate Peierls creep rheology has a new function compute_approximate_log_strain_rate_and_derivative, mirroring compute_exact_log_strain_rate_and_derivative already in the code.
• A new power law Drucker-Prager rheology is introduced, that replaces a strain-rate invariant yield stress with a yield stress that varies as a function of strain rate.
• CompositeViscoPlastic is modified to use the approximate Peierls creep rheology and the power law Drucker-Prager rheology.
• The internal Newton iteration scheme inside CompositeViscoPlastic is modified to iterate over the log stress to get the correct total log strain rate. The change from iterating on stress and strain to log stress and strain, and the use of pure power law rheologies (at fixed state) leads to much improved performance, typically requiring 1--6 relatively inexpensive iterations rather than 1--20+ more expensive iterations with the old code.

Further improvements in efficiency can be made in future PRs, particularly in the approximate Peierls creep rheology (see #5945).

I'm happy to split this PR up if necessary. All brought together here because I developed the features in tandem, and to demonstrate the eventual goal.

• I have followed the instructions for indenting my code.
• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[bobmyhill]

/rebuild

[gassmoeller]

Looks very promising and I am curious how it performs in practice. I left a bunch of comments, but mostly small renaming and documentation requests.

[gassmoeller]

So at low temperatures the changes lead to higher viscosities and stresses. Is that because of the power law of the drucker prager rheology? Do you approach the old values as your stress exponent goes to infinity? I just want to make sure the difference we see is caused by the functional change in the drucker prager rheology, and not by the change to the Newton iteration.

[bobmyhill]

Yes and yes.

• As the temperature goes down, the purely viscous stress at a given strain rate increases.
• The "yield stress" in the new Drucker Prager flow law is associated with a constant reference strain rate
• At a fixed strain rate, the ratio of plastic strain to viscous strain will increase as temperature decreases.

At even lower temperatures, plastic strain increasingly will dominate, and so the stress will tend towards a constant value. That value will depend on the "yield stress", reference strain rate and actual strain rate.

[bobmyhill]

Hi @gassmoeller , thanks for the detailed review!
I've now added the variable grain size we talked about at the user meeting. This PR is increasingly becoming a testing ground for all the features we need, so I'm happy to split it up when we're happy it's the way we want it to go.

[bobmyhill]

@gassmoeller, @anne-glerum

Today @anne-glerum pointed me at the stress limiter in the viscoplastic rheology:

aspect/source/material_model/rheology/visco_plastic.cc

Lines 334 to 342 in e7ff10c

	case stress_limiter:
	{
	//Step 5b-1: always rescale the viscosity back to the yield surface
	const double viscosity_limiter = yield_stress / (2.0 * ref_strain_rate)
	* std::pow((edot_ii/ref_strain_rate),
	1./exponents_stress_limiter[j] - 1.0);
	effective_viscosity = 1. / ( 1./viscosity_limiter + 1./non_yielding_viscosity);
	break;
	}

This is exactly the same as DruckerPragerPower:

aspect/source/material_model/rheology/drucker_prager_power.cc

Lines 101 to 112 in 61bc576

	DruckerPragerPower<dim>::compute_viscosity (const double cohesion,
	const double angle_internal_friction,
	const double pressure,
	const double effective_strain_rate,
	const double max_yield_stress) const
	{
	const double yield_stress = compute_yield_stress(cohesion, angle_internal_friction, pressure, max_yield_stress);
	const double stress = yield_stress * std::pow(effective_strain_rate/drucker_prager_edot_ref, 1./drucker_prager_stress_exponent);
	const double apparent_viscosity = stress / (2. * effective_strain_rate);
	return apparent_viscosity;
	}

What do you think about renaming DruckerPragerPower to StressLimiter and then using it in both ViscoPlastic and CompositeViscoPlastic?

The only thing that bothers me is that in visco_plastic, the parameter "Reference strain rate" is used both to initialise the strain rate in the first time step, and also to set the reference strain rate in the stress limiter. I don't know what the best solution is:

• to add a new parameter (changing existing behaviour),
• for ViscoPlastic to inherit the parameter value from StressLimiter (we agreed that this was a bad idea™ at the hack).
• To not use the new StressLimiter rheology in ViscoPlastic (seems like a wasted opportunity)

Any thoughts?

[gassmoeller]

What do you think about renaming DruckerPragerPower to StressLimiter and then using it in both ViscoPlastic and CompositeViscoPlastic?

Yes, reusing code is always good.

The only thing that bothers me is that in visco_plastic, the parameter "Reference strain rate" is used both to initialise the strain rate in the first time step, and also to set the reference strain rate in the stress limiter. I don't know what the best solution is:
Any thoughts?

That seems like an earlier oversight, in particular because the use in the stress limiter is not even documented in the input parameter comment (it only mentions that the strain rate is used during the first time step). I would say let's introduce the new parameter and use that one.

For a future PR: In order to use the DruckerPragerPower as stress limiter in ViscoPlastic do you need a different limiting stress for different compositional fields? (I just noticed that the stress limiter seems to allow different values for different fields, while DruckerPragerPower seems to use a single yield stress).

This PR is increasingly becoming a testing ground for all the features we need, so I'm happy to split it up when we're happy it's the way we want it to go.

Yes, let's focus on one thing at a time. As far as I can see there are three things mixed up in this PR:

• grain size for composite visco plastic
• log stress Newton iterations
• DruckerPragerPower

It would be ideal to separate this in 3 PRs, but you know better what can be cleanly separated and what not. I will leave that decision up to you, but please do not add anything else in this PR (like using DruckerPragerPower in visco_plastic).

[bobmyhill]

Thanks @gassmoeller , that all seems sensible.
I noticed the different number of stress exponents in VP, and raised an issue here: #5961. I don't think there's a particular need for that (maybe you disagree) but I do think I should implement different Max yield stresses. That isn't in the base DP rheology yet though. Is that what you meant?

[gassmoeller]

I don't think there's a particular need for that (maybe you disagree)

I do not disagree, but I would like to check with @anne-glerum and @naliboff if those different exponents were implemented on purpose and if they have use cases for them. Removing unnecessary functionality is good, but removing functionality that is needed and/or already in production use would be bad.

but I do think I should implement different Max yield stresses. That isn't in the base DP rheology yet though. Is that what you meant?

That wasnt what I meant, but it seems reasonable and more useful than the different exponents. Something for a follow-up PR though.

[bobmyhill]

Addressed by #5978.