fully compressible Stokes equation

[YiminJin]

This PR implements the fully compressible Stokes equation. The mass conservation equation is expressed as
$$-\nabla\cdot\mathbf{u} = \frac{1}{\rho}\nabla\rho\cdot\mathbf{u} = \frac{\beta(p - p^{\text{old}})}{\Delta t}. (1)$$
The second "=" has utilized the approximation
$$\beta = \frac{1}{\rho}\frac{\partial\rho}{\partial p} \approx \frac{1}{\rho}\frac{\text{d}\rho}{\text{d}p}. (2)$$
To maintain the efficiency of the linear solver, I modified the Schur complement preconditioner by (according to "Finite Elements and Fast Iterative Solvers", H. C. Elman et al., 2014, Eq. [8.11]):
$$\mathbf{S} = \mathbf{B}\mathbf{F}^{-1}\mathbf{B}^T + \mathbf{C}, (3)$$
where $\mathbf{C}$ is the bottom right block of the system matrix. I have tested the formulation by a simple mantle convection model (similar to the experiment in https://doi.org/10.1093/gji/ggaa078):
fully_compressible.prm.txt
and it works well. The result is a little different from the isentropic compression and hydrostatic compression approximations, but the computation time is almost the same.

The fully compressible mass conservation equation is necessary for transient problems and plastic dilation, and has been implemented in a number of geodynamic codes, such as ELVIS and SLIM3D. Comparing to the other codes, my implementation is simple and not strict (Eq. (2) may have some problems). Also, the fully compressible formulation has not been implemented in the GMG solver. Please check this implementation if you are interested in it @naliboff @anne-glerum @bobmyhill

[bobmyhill]

Dear @YiminJin, thank you for posting this for us to look at!

I wonder whether you have looked at the Projected Density Approximation (PDA) in @gassmoeller's paper from 2020 (https://academic.oup.com/gji/article/221/2/1264/5735442), described in Section 2.5.3, which is implemented in ASPECT. There are some important comments in that section about pressure and temperature variations in the density and oscillations when using the full density in convection simulations. There is also a brief description of the limitations of the PDA.

We have previously discussed trying some problems using the full density (and time-dependent variation of the density) in all governing equations. For example, we've discussed papers on lithospheric overpressure* (which is viewed by many as implausible and others as inevitable). If you are particularly interested in problems that need full compressibility, @gassmoeller will have some good thoughts.

*e.g https://www.researchgate.net/profile/Torgeir-Andersen/publication/260694635_An_alternative_model_for_ultra-high_pressure_in_the_Svartberget_Fe-Ti_garnet-peridotite_Western_Gneiss_Region_Norway/links/0c9605320366a6749f000000/An-alternative-model-for-ultra-high-pressure-in-the-Svartberget-Fe-Ti-garnet-peridotite-Western-Gneiss-Region-Norway.pdf

[gassmoeller]

I will need some time to review this (or #5929) since I am currently traveling. I just wanted to say we appreciate the contributions @YiminJin!

[bobmyhill]

P.S. I realised that you have obviously read the 2020 paper as you cite it in your first post. I wonder whether you might clarify a couple of things for me:

• Your expression doesn't include thermal expansivity (e.g. the thermal divergence in Duretz; https://insu.hal.science/insu-03283712/document, Equation C1). Is there a reason for this?
• What is the reason that plastic dilation requires a fully compressible formulation? In the Duretz paper I note that they try incompressible and fully compressible, but nothing in between, and I'd be surprised if the PDA + "plastic dilatancy" didn't do a good job without needing full compressibility.

[YiminJin]

@Rombur Thanks for your comments. About your questions:

• Your expression doesn't include thermal expansivity.

It is actually a problem that should be corrected afterwards. However, I've noticed in the ASPECT documentation (section "Coefficient self-consistency", one of your contributions) that the formulations are quite complicated, which may take me a few days to understand. Currently I am more interested in plastic dilation, and I plan to implement a version without thermal expansion as a first step. (I have met problems relevant to plasticity in my own work, so I want to do it fast.)

• What is the reason that plastic dilation requires a fully compressible formulation?

This is because when the trial stress is mapped onto the yield surface, the direction of plastic strain should be perpendicular to the "equipotential surface" of plastic potential. When the dilatancy angle is not 0, the return-mapping of stress will change the pressure directly. In this case, not only the top-left block, but also the top-right and bottom-right blocks of the system matrix should be supplemented with compressibility-relevant terms. If one uses approximated formulations, I doubt if the pressure would change in the correct way and the convergence rate would be as good as the fully-compressible model (like SLIM3D and Duretz's code).

[bobmyhill]

Hey @YiminJin, thanks for your answers! They are really useful.

I wouldn't worry about the "Coefficient self-consistency" (the isothermal and isentropic compressibilities are almost identical at lithospheric temperatures), but you're right, the thermodynamics of the problem does need some thought.

For compact materials (i.e. those with no void space) the change in volume can be decomposed (nonuniquely) into pressure, temperature (or entropy) and composition terms. Your equation (1) contains the pressure term, the temperature term involves the thermal volumetric expansivity $\alpha_P$, and the compositional term involves partial volumes of the components involved in the reaction. We sometimes fold the compositional terms into the P-T/S terms at constant bulk composition, see, for example the recent paper from @jdannberg, @gassmoeller and @RanpengLi with Carolina Lithgow-Bertelloni and Lars Stixrude; https://zenodo.org/records/6494835).

As I understand it, plastic dilation implies that materials are no longer compact - i.e., this is the modelling analogue of fault-generated porosity. If that is correct, (a) an extra "viscoplastic divergence" term is required in the mass conservation equation, (b) some evolution law will be required for the generation of porosity, and (c) an extra field will need to be tracked (either the porosity or some related property).

All of the above is to motivate another question to your comment:

"not only the top-left block, but also the top-right and bottom-right blocks of the system matrix should be supplemented with compressibility-relevant terms."

Do you think that the "compact-equivalent" beta and alpha terms are locally important (in mapping the trial stress onto the yield surface)? Or is their importance restricted to non local effects (loading and unloading across the shear zone)? Comparing Figure 3a and 3b of the Duretz et al. (2021) paper, finite bulk compressibility moves the ripple instabilities around, but it doesn't kill them, which leads me to think that maybe bulk (compact) compressibility is not so important on a local scale.

Whatever the answer, I think we do want to play with fully compressible formulations in ASPECT, so thank you for opening this PR, and taking the time to answer my perhaps naive questions :)

[YiminJin]

Thanks @bobmyhill . Your questions are not naive, but provide me a lot of useful information.

I made a few mistakes in the previous comment. First of all, I think the plastic dilation does not change the Stokes system in Picard iterations --- the Picard method does not provide a descent direction for the solution. The return-mapping procedure is only useful in Newton method. Thus, I think the plastic dilation does not affect the "compactness" of the material.
Secondly, after some mathematical work I found that all the 4 blocks of the Newton system need to be supplemented with dilation-relevant terms. The calculations are shown in the attached pdf:
dilation.pdf
The calculations are very primative --- they did not take cell averaging into account. I am trying to add the dilation-relevant terms to the Newton system and see what happens to the solution.
I have a small question: the comment for class PrescribedPlasticDilation mentioned a paper of ChoiPeterson (2015), but I failed to find the paper. Would you please provide me the doi of the paper @naliboff ?

[YiminJin]

@bobmyhill Sorry, another mistake: the plastic dilation DOES affect the Stokes RHS in Picard iteration (see the last term of $F_p$ in Eq. (16) of the pdf). And now I understand the PrescribedPlasticDilation.

[cedrict]

@YiminJin the paper is http://dx.doi.org/10.1016/j.tecto.2015.06.026

[YiminJin]

Thank you @cedrict for the paper!
@bobmyhill I shall apologize for the rush uploading of the pdf which is full of mistakes. I directly transformed the basic equations of elASPECT to the framework of ASPECT but did not check the physics carefully. I hope nobody has wasted time on reading it.
I carried out some more tests and found that the fully compressible equation causes severe constant drift when pressure normalization is off. I don't know how to solve the problem, so I decide to close this PR. The plastic dilation issue will be discussed in a new PR, which has passed the strip-footing test.

[gassmoeller]

Hi @YiminJin you are of course welcome to close this PR if you think it is not useful anymore. I wanted to say that your continuity approximation itself is probably correct for a subset of very special cases (i.e. no temperature or chemistry changes as Bob already pointed out), however whether you think this is relevant enough to be merged into the main ASPECT version is up to you (i.e. if you want to run models that need this approximation and it may be useful for others I am open to consider it). Let me know in case you want me to take a closer look at this PR at some point in the future.