[aspect-devel] Help needed with setting up an EBA benchmark

[Wolfgang Bangerth]

All,
I need a bit of help from people who understand the 
nondimensionalization procedure in geodynamics better than I do. I want 
to set up the Davies et al. benchmark 3.1 outlined in

https://github.com/bangerth/aspect/blob/more-davies-benchmarks/benchmark/davies_et_al/Cylindrical_Benchmark.pdf?raw=true

It uses the extended Boussinesq approximation (EBA) that includes 
adiabatic heating and friction heating terms. The temperature equation 
stated in this document (equation (17), p. 3) looks like this:

[1]
   cp DT/Dt  -  Di alpha g . u (T+T0) - div(k grad T) = Di/Ra Phi

where
   Phi = 2 eta eps(u):eps(u)
is the friction heating term. The problem, of course and as always, is 
that the benchmark is stated in terms of Dissipation and Rayleigh 
numbers Di and Ra, instead of the physical quantities we use in ASPECT. 
In particular, the equation we use has the form

[2]
   rho cp DT/Dt - alpha T u . grad p - div(k grad T) = Phi

So I need to finagle [1] into form [2]. I believe that the authors of 
[1] have in mind that in their form cp=alpha=|g|=k=1. They prescribe 
Ra=10^4 and Di=0.25.

Then, if I use these values for the physical quantities and multiply 
their equation by Ra/Di, I get

[1']
   Ra/Di DT/Dt  -  Ra g . u (T+T0) - div(Ra/Di grad T) = Phi

So, to put this into form [2], I do
   g=10^10
   alpha=10^-6
      --> these two guarantee that rho=1 to good accuracy because T=0..1
   cp=Ra/Di=4.10^4
   k=Ra/Di=4.10^4
With this, I can match the first and third terms of [1'] and [2]. 
However, I'm stuck at the second term. I do have that
   grad p  =  -rho g
(I assuming that Davies et al. define g to point upward, whereas we 
define it as pointing downward) so the terms look like this:

[adiabatic heating in 1' -- Davies et al.]
   - Ra g . u (T+T0)               (g points upward)
[adiabatic heating in 2  -- ASPECT]
   + alpha rho g . u T             (g points downward)

With the values I have chosen above, alpha*rho*g=Ra, so this part is 
good. But how do I get the factor T0 into my equations? The benchmark 
defines it as T0=0.091 (page 4). This term simply doesn't appear 
anywhere in ASPECT. Anyone's got an idea how to finagle it into ASPECT 
without having to change the code?

Thanks in advance
  Wolfgang


PS: Model setup file is here:
https://github.com/bangerth/aspect/blob/more-davies-benchmarks/benchmark/davies_et_al/case-3.1.prm

-- 
------------------------------------------------------------------------
Wolfgang Bangerth               email:            bangerth at math.tamu.edu
                                 www: http://www.math.tamu.edu/~bangerth/