Governing equation approximations

[ian-r-rose]

Hi all,

So we have discussed this in one form or another a couple of times, but I thought that I would bring it up again. As we know, there are several approximations to the equations for mass, momentum, and energy that are commonly used in mantle convection (Boussinesq, anelastic, etc). Some are more thermodynamically consistent than others. Discussing this with Juliane and Rene this week, I think we agreed that this could be clarified in Aspect. This recently reared its head when trying to reproduce the Zhong et al 2008 benchmark, as the "simple" material model is not Boussinesq in the strictest sense (does not use the reference density for thermal diffusivity, as I recall? I have not verified, but this is what Rene indicated to me). I think it is important to be able to use these approximations with confidence.

I see two different ways of addressing this:

1. The easier solution. Try to introduce a Boussinesq material model which does exactly that and nothing more. The reference density is used everywhere except for the terms which multiply gravity. Additionally, it would be nice for Aspect to query the material model about its compressibility, look at whether viscous heating and adiabatic heating are included, and then infer which approximation is being used. This would be useful output information at the start of a run.

2. Fully support the relevant approximations that are used (Boussinesq, extended Boussinesq, anelastic) so that it can be specified in a parameter file. This would require some extensive changes to the assembly of the matrix, as well as some checks for compressibility in the material model, but it may be worth it.

I am willing to work on this, but I'd like to hear thoughts from people as to what would be most appropriate before starting.

[gassmoeller]

I agree that we should provide some more controls / reasonable defaults on this. Not only can one set up inconsistent models at the moment, but we also might need to change the formulation slightly to reproduce benchmarks, they often ignore dynamic pressure for adiabatic heating. However I do not know whether the material model would be the most appropriate place to put it in. They are already quite long. Also it is not intrinsically a material property, but rather a different way to formulate the equations so the assembly seems the most reasonable place for it.
In general I see 3 different ways to formulate the heating in a model:

1. No heating and no shear dissipation (for incompressible models using Boussinesq approximation)
2. Shear heating and adiabatic heating with adiabatic heating being simplified by neglecting dynamic pressure (replacing \Delta P by \rho g)
3. Shear heating and adiabatic heating with full pressure (as it is implemented today)
   The shear heating hereby needs to be calculated differently in the compressible and incompressible case as it is already implemented.
   Additionally a material model might either provide latent heating information or not.

Maybe we can simply introduce a function in the material model interface that returns the kind of heating assumed by / implemented in the material model (as an enum, as for example for the NonlinearSolver) and then accordingly compute the heating term in the assembly? This would remove the need to have them as model input parameters and prevent users from using settings that are not handled consistently in the material model. Material models that support different approximations might still allow for the heating formulation to be a user input.

Setting up a very simple material model fully implementing the Boussinesq approximation (and others for EBA, TALA) still would be a great thing and would be possible independently.

[tjhei]

Any news on this? We really do need an easy way for the user to specify BA, EBA, TALA and also see which method is used for a given computation.

[bangerth]

The changes really have to happen in the simulator proper because there one has to decide whether one wants to use the actual density or the reference density for a particular term.

The thing that makes me unhappy is that we've pretty much gotten rid of using any of the reference_* functions of the material models and I'd be much happier if we got rid of them altogether. On the other hand, if people want to run simulations that employ the various simplified models, then I don't see any other options than reintroducing them or using the workarounds we use right now (e.g., g=10^10, alpha=10^-6, so that the density used in the temperature equation is essentially constant).

[ian-r-rose]

Just so I'm clear, the only change that would need to be made to the matrix assembly is using a reference density rather than the total density in a few places? That doesn't sound too bad, overall. I am still in favor of directly supporting these approximations (after all, the Stokes equation itself is an approximation).

The workarounds that you describe are confusing; I had seen them before and not understood why those parameters were so weird until just now. Also, those extreme values for certain parameters caused issues with the solver during the hackathon. In theory, we've solved that, but perhaps there are other corner cases that we've missed.

As far as I can see, there are two major changes that need to be made in order to support BA, EBA, TALA, ALA:

1. Define a reference density profile. Trivial for Boussinesq and friends, as it is constant. Trickier for anelastic, as we don't just want a constant reference density, but the depth-dependent reference density. Perhaps we can reuse some of the machinery from adiabatic_conditions? Or require that a material model implement a depth-dependent reference_density function if they want to use ALA?

2. Provide an easy interface for the user. Personally, I'd be interested in just specifying the governing equations at the top of the parameter file, and then ditching parameters like include_shear_heating and include_adiabatic_heating.

[gassmoeller]

To Ians point 1): Would it be fair to expect the density function to return the 'reference density' whenever it is called with a 'reference temperature' and 'reference pressure'? Because then we would not need separate reference_density functions, we could simply use the existing function to also calculate a value 'reference_density', which could then be used in the assembly instead of the density in case we ordered a BA/EBA/TALA/ALA model. We would need to think about what kind of reference temperature and pressure would be appropriate, but in the compressible or EBA case this could be the adiabatic pressure and temperature (from adiabatic conditions), and in the BA case, pressure should not matter (density should be independent) and temperature should be constant (Maybe 'Adiabatic surface temperature'? Or the 'Reference temperature' of the material model, but this is not necessarily defined in all material models). We could then think about either define a reference_density_profile inside adiabatic_conditions (as it is done for pressure and temperature, and we have the density at Tadi and Padi there anyway), or let the material model return another variable in materialmodeloutputs that is called reference_density. Actually I think the first one might be more efficient, since we really only need to compute the reference values once.

To 2): There should definitely be an option to select the different approximations, but I am not sure ditching the include_shear_heating and include_adiabatic_heating is a good idea. There might be cases where one could want to use them separately, be it for testing purposes (comparing temperature increase with heating amount), be it because a growing number of people uses Aspect also for lithosphere scale modelling and I am not aware of the kind of approximations that are common in that community. Think for example of a 2D simplified plate, driven by velocity boundary conditions and open side walls to investigate the shear heating in the asthenosphere. There would be not vertical flow (so no need for compressibility or adiabatic heating) but the shear heating is an essential part of the model. I would therefore argue that we at least keep the possibility to set the approximation to 'none' or 'individual' and set the shear heating and adiabatic heating individually.

[jdannberg]

Okay, we have thought about this a bit more and we have come to the conclusion that (1) does not work in the way we thought, e.g. in the EBA case, the reference temperature should be adiabatic, but the reference density should be constant, so we can not get the reference density by just calling the density in the material model with the reference temperature and pressure.

This means, if we want to support these approximations, we need to define a reference density somewhere, which depends on what approximation we want to use. As the adiabatic temperature is the reference temperature for EBA, the adiabatic conditions might be a good place to do this. The only additional change would be in the assembly, to ask whether to use the 'normal' density everywhere or just in the gravity term, and the reference density everywhere else.

[gassmoeller]

I have uploaded a first implementation of how the reference density inside of the adiabatic conditions could work (#212). Feel free to criticize. This does not mean we need to go this way, I just wanted to illustrate, what I imagined. There is a complication at the moment about the meaning of the adiabatic pressure profile that is not directly related to this issue but came up, when I implemented this. We do not clarify whether the adiabatic pressure profile is meant as reference value for lookup tables (using compressible density even in incompressible models) or as initial guess for the stokes solver (using the same density as the model approximation). Sometimes this conflicts.

[gassmoeller]

Closed by #1275. Just took 2 and a half years since opening this issue 😄.

[bangerth]

👍