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

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

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

1 Like

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.

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

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

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

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.

Hi Rene,

Thank you!

– G

Hi Rene,

I’m trying to write in our paper your explanation of how the coefficients \alpha_{max} and \alpha_E, etc. were chosen. In particular, I am trying to understand what the coefficient or parameter that you call cR is.

But I can’t find a reference to cR in either the 2012 or 2017 papers.

I understand that what is called \alpha_{max} in the 2012 paper you are referring to as \beta. What are you now calling the parameter that is called \alpha_E in the 2012 paper?

Is that you refer to as cR? In other words is referred to as \alpha_{max} in the 2012 paper is what you are calling cR?

Thanks!

  • G

Gerry, have you checked the manual / parameter documentation in ASPECT?
Let me quote:

  • set cR = 0.33 [List list of <[Double 0…MAX_DOUBLE (inclusive)]> of length 0…4294967295 (inclusive)]

    • The c_R factor in the entropy viscosity stabilization. This parameter controls the part of the entropy viscosity that depends on the solution field itself and its residual in addition to the cell diameter and the maximum velocity in the cell. This parameter can be given as a single value or as a list with as many entries as one plus the number of compositional fields. In the former case all advection fields use the same stabilization parameters, in the latter case each field (temperature first, then all compositions) use individual parameters. This can be useful to reduce the stabilization for the temperature, which already has some physical diffusion. (For historical reasons, the name used here is different from the one used in the 2012 paper by Kronbichler, Heister and Bangerth that describes ASPECT, see \cite{KHB12}. This parameter corresponds to the factor \alpha_E in the formulas following equation (15) of the paper. After further experiments, we have also chosen to use a different value than described there.) Units: None.
  • set beta = 0.078 [List list of <[Double 0…MAX_DOUBLE (inclusive)]> of length 0…4294967295 (inclusive)]

    • The \beta factor in the artificial viscosity stabilization. This parameter controls the maximum dissipation of the entropy viscosity, which is the part that only scales with the cell diameter and the maximum velocity in the cell, but does not depend on the solution field itself or its residual. An appropriate value for 2d is 0.078 and 0.117 for 3d. (For historical reasons, the name used here is different from the one used in the 2012 paper by Kronbichler, Heister and Bangerth that describes ASPECT, see \cite{KHB12}. This parameter can be given as a single value or as a list with as many entries as one plus the number of compositional fields. In the former case all advection fields use the same stabilization parameters, in the latter case each field (temperature first, then all compositions) use individual parameters. This can be useful to reduce the stabilization for the temperature, which already has some physical diffusion. This parameter corresponds to the factor \alpha_{\text{max}} in the formulas following equation (15) of the paper. After further experiments, we have also chosen to use a different value than described there: It can be chosen as stated there for uniformly refined meshes, but it needs to be chosen larger if the mesh has cells that are not squares or cubes.) Units: None.

Hi Timo,

That was going to be my next step, which I left for today. :slight_smile:

I didn’t expect the manual entry to be so detailed that it would state:

For historical reasons, the name used here is different from the one used in the 2012
paper by Kronbichler, Heister and Bangerth that describes ASPECT, see \cite{KHB12}.
This parameter corresponds to the factor αE in the formulas following equation (15) of the
paper.

That is why I wanted to double check with Rene. These manual entries are great!

May I now ask you an unrelated question? How did you get into the math / LaTeX mode so
that cR and αE printed correctly, with R and E as subscripts?

Thank you!

  • Gerry

By writing latex between dollar signs: \sqrt{\frac{1}{\pi}}

Rene has explained how the default value of c_R was chosen; namely,

The entropy viscosity parameters have indeed changed since the first ASPECT paper, 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 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 c_R that led to a constant or decreasing maximum and minimum composition value over time.

Could anyone tell me how the default value of \beta was chosen?

I think the description in the results section of step-31 is still valid:
https://dealii.org/developer/doxygen/deal.II/step_31.html#Numericalexperimentstodetermineoptimalparameters
The exact value of beta may later have been modified, but the methodology to
choose them was similar to the one discussed in the link above.

Best
W.

1 Like

Thanks for the quick response!

I plan to reference that in my paper. I am trying a BibTeX entry of the following form

@MISC{BETA-CR:2019,
TITLE = {Numerical experiments to determine optimal parameters},
YEAR = {2019},
howpublished = {{R}eference documentation for deal.II},
note = {{T}he step-31 tutorial program},
}

However, the entry in the bibliography doesn’t look good:

References
[1] , 2019. Numerical experiments to determine optimal parameters. Reference
documentation for deal.II, The step-31 tutorial program.

Note the comma after [1]. I think if I add an author or authors it would be better. If this is OK who should I add as author(s)?

I have previously used the following entries for step-31 and step-32:

@Misc{dealiistep31,
author = {M. Kronbichler and W. Bangerth},
title = {deal.{II} tutorial program step-31,
\url{http://www.dealii.org/developer/doxygen/deal.II/step_31.html}},
year = 2011
}

@Misc{dealiistep32,
author = {M. Kronbichler and T. Heister and W. Bangerth},
title = {deal.{II} tutorial program step-32,
\url{http://www.dealii.org/developer/doxygen/deal.II/step_32.html}},
year = 2011
}

1 Like

Thank you! :heart:

Post must be at least 20 characters
Have you tried the :heart: button?

I’m not sure why the heart button is the preferred button, but what the heck! :smiley:

I forwarded the post-length issue to Tyler and he fixed it. Any length should be fine from now on.

Thanks! (Just testing …) :slight_smile: