geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
223 stars 236 forks source link

Add ln(strain rate) and d ln(strain rate) / d ln(stress) functions in potentially composite rheologies #5322

Open bobmyhill opened 1 year ago

bobmyhill commented 1 year ago

After some testing, @lhy11009 has found that iterating ln(stress norm) to minimize the residual of ln(strain rate norm) is more efficient than iterating stress norm to minimize the residual of strain rate norm. This makes sense, as ln(strain rate) is close to linear in ln(stress). This mirrors the iterative scheme in the grain size model (ln(grain size)) and possible future temperature iterations in the multicomponent entropy model (entropy at high temperature is linear in ln(T)).

So that we don't have to duplicate / replace all the functions in the individual rheology modules, we should simply use the following identity in the iterative loops: d ln(strain rate) / d ln(stress) = (stress / strain rate) * (d (strain rate) / d (stress))

Here's a list of places we should implement the change. @lhy11009 has kindly offered to do the diffusion-dislocation model.

bobmyhill commented 1 year ago

Follow-up: To avoid writing a lot of new functions, we could as a stopgap use the identity: dlny / dlnx = (x / y) * dy / dx

lhy11009 commented 1 year ago

@bobmyhill I just looked at the composite visco plastic material model and It seems to me the issue it has is different. It computes the composite viscosity directly by harmonic average among components so there is no internal iteration from strain rate to stress. However, the compute_viscosity_derivatives function computes the derivatives of viscosity to strain rate and pressure, this seems to be what the Newton solver needs to assemble. I guess here are two questions:

  1. Do we still need to modify the composite viscoplastic as well under the current workflow mentioned in this issue?

  2. Could we modify the compute_viscosity_derivatives and let it computes d_log_viscosity / d_log_strain_rate and d_log_viscosity / d P and change the way the newton solver assembles? This is what I intend to do. Yes, we are switching, but if the Newton solver we have is easy to tweak, it could be a fantastic testing ground for this. Even if I can make converge a little faster, it would directly help my research. Let me know what you think.

FYI @mibillen

mibillen commented 1 year ago

@Haoyuan, thanks for for looping me on this. From our discussion on this, I do think it is worth trying this out to see the effect on the Newton Solver. The performance of the newton Solver with our subduction models has always been very disappointing, and it was not clear why it should be. @bobmyhill - your input and guidance on this is for Haoyuan is very much appreciated.

bobmyhill commented 1 year ago

@lhy11009 thanks for looking at this.

1) The composite VP scheme is this one: source/material_model/rheology/composite_visco_plastic.cc rather than the one currently implemented as visco_plastic. 2) This is an interesting idea. Let's see if the log scheme is generally useful first, then consider implementing other changes. I think using SUNDIALS might be a good next thing to implement - its internal line search algorithm is likely to stabilise these internal solves.

lhy11009 commented 1 year ago

@MFraters you may also be interested in what we are doing here.