Have the default values of the entropy viscosity parameters in ASPECT changed?


#1

On behalf of Gerry Puckett and received on the deprecated mailing list:

Hi,

I’m writing a paper with Jon Robey for Computers & Fluids that describes the implementation of his VOF interface tracking algorithm in ASPECT in great detail. One of the reviewers wrote

  • In Section 3.3, the authors describe the use of [the] entropy
    viscosity stabilization technique in the energy (or temperature)
    equation. … I believe the use of such stabilization method[s,]
    akin to artificial diffusion is justifiable only if it does not
    significantly change the engineering results of the problem. The
    authors need to demonstrate that.

  • On the related note, the authors also need to show the functional
    form of this entropy and its
    dependence to the local Péclet number. A comment on how they select
    the stabilization parameters
    values are also needed for completeness.

Following the discussion in

[40] Kronbichler, M., Heister, T., Bangerth, W., 2012. High accuracy
mantle convection simulation through modern numerical methods.
Geophysical Journal International 191 (1), 12–29.

we have gone into great detail to address the first sentence in 2. above. In particular, in reference [40] on page 18 it is stated that

1. ... the constant \alpha_max = 0.026d depends only on the spatial
   dimension d.

2. ... we choose the constant \alpha_max _E = 1, see also the
   discussion in Guermond et al. (2011).

However, it seems that Jon has found that in his computations

3. "For the temperature-dependent models in this paper, the values
   chosen for the stabilization parameters are the default values
   for ASPECT, \alpha_max = 0.078, and \alpha_E = 0.33."

I have asked Jon to double check that this is what is in the parameter file that is output with the values of the parameters that have been used in a particular computation.

Nevertheless, I thought it would be a good idea to ask on aspect-devel what current defaults are.

One other observation. I noticed that

if the dimension d = 3, then we have,

   \alpha_max = 0.078 = 0.026 * 3 = = 0.026 * d

and then the statement in 1. concerning how \alpha_max is chosen is
consistent with the statement in 3., /except that our d should be /d =2.

Also is there a simple way to address this request from the reviewer?

“A comment on how they select the stabilization parameters values
are also needed for completeness.”

Comments, thoughts, explanations … ?

Thanks for your attention to this issue. I’m trying to ensure that what we write in a refereed journal article about ASPECT is correct, at least at the time at which we wrote it.

Happy New Year!

Gerry


#2

Hi Gerry,

The entropy viscosity parameters have indeed changed since the first ASPECT paper, (that happened in 2012 already, see https://aspect.geodynamics.org/doc/doxygen/changes_between_0_81_and_0_82.html) mostly because the initial parameters were chosen for the temperature equation that contains some natural diffusion. When using these parameters for the composition equation we noticed that oscillations could grow. The new parameters were chosen (if I remember correctly) by running the van Keken Rayleigh-Taylor overturn problem (Keken, P.V., King, S.D., Schmeling, H., Christensen, U.R., Neumeister, D. and Doin, M.P., 1997. A comparison of methods for the modeling of thermochemical convection. Journal of Geophysical Research: Solid Earth , 102 (B10), pp.22477-22495.) and choosing the smallest values for cR that led to a constant or decreasing maximum and minimum composition value over time.

You are right that the value of beta is currently hardcoded and not changed in dependence of dimensionality as it could be (this would require a bit of technical adjustments). It is also set to 0.078, which was the old recommended value for a 3D model, but is now the recommended value for 2D models as mentioned in the description of this parameter in the manual (so you are fine using it for 2D models). Of course that means the default is actually not optimal for 3D models, but if I remember correctly, beta is more of a hard upper bound for the stabilization in regions with extremely sharp gradients. Nearly all cells do not reach this bound, and so the stabilization is instead governed by cR. It has been a while since I looked into this though, you might want to check this.

In short, while the basic principles of the current stabilization in ASPECT are still the same the parameters have been adjusted to not rely on natural diffusion to provide sufficient stabilization and the formula also has changed a bit (it now also includes a dependence on the strain rate, but that is switched off by default).

Hope that helps,
Rene


#3

Hi, I am trying to determine why (I think) Ying wrote in a previous paper

In our work we do not use a second-order extrapolation to treat the advection
term (u · ∇T, ψ) and the entropy viscosity term (ν h (T ) ∇T , ∇ψ) Ω explicitly
in time as described in [41]. We use the fully implicit adaptive Backward Differentiation Formula of order 2 (BDF2) [41, 78] to discretize the temperature equation in time.
Jon does not know where this came from, and I find it hard to believe that Ying recoded one small piece of the entropy viscosity.

Is this also a default change that was made after

[41] Kronbichler, M., Heister, T., Bangerth, W., 2012. High accuracy mantle convection simulation through modern numerical methods. Geophysical Journal International 191 (1), 12–29.

was published.


#4

Hi Gerry,
yes that is another change that was introduced after the publication of the 2012 paper. See Section 3.1 of the second ASPECT paper:

Heister, T., Dannberg, J., Gassmöller, R., & Bangerth, W. (2017). High accuracy mantle convection simulation through modern numerical methods–II: realistic models and problems. Geophysical Journal International , 210 (2), 833-851.

Best,
Rene


#5

You will need to ask Ying. I have no idea.
Magali


#6

Rene, May I use you response to my question concerning “How the new parameters were chosen” in my response to a reviewer’s question

    "A comment on how they select the stabilization parameters 
     values are also needed "for completeness."

Also, do you or anyone else know the answer to

    (if I remember correctly)

Are the statements following them correct, or should I put the question to the forum community?

Thanks!

  • G

#7

Hi Gerry,
I found an email thread between Wolfgang and me from 11/12/2012 that discussed these changes and yes you can use my answer for your reply, I remembered correctly.


#8

Hi Rene,

Thank you!

– G