Question about Rayleigh number calculation in basic statistics postprocessor

ASPECT
1

Hi team,

I have a question about the Rayleigh number that the basic statistics postprocessor returns. To demonstrate the question I’ve run the shell_simple_3d.prm cookbook with slight modification (I’m pasting the content of the parameter file below since as a ‘new user’ of the form I don’t seem to be allowed to add attachments). Following the equation in the manual, the material model should result in a Rayleigh number of 0.910e+07.
However, when running the parameter file it returns a Rayleigh number of 1.00e+07 (see log file). The reason for this is that the Reference density is returned as 3458.4 instead of 3300, which is the value of the material model. Looking at the simple material model, the reference density is returned as an average:

out.densities[i] = MaterialUtilities::average_value(volume_fractions, eos_outputs.densities, MaterialUtilities::arithmetic);

Given that the reference temperature is 1600, a density of 3458.4 is obtained at T = 0. I’m not really sure how this can be the average over any part of the model domain and is not reflective of the total average density. I don’t know the averaging function too well so wanted to check with you all if you see what’s going wrong. Why is the reference density not closer to 3300?

The difference in reference density also propagates into a reference thermal diffusivity that differs between what the log prints vs. what would be calculated based on the parameter file.

Thank you for your help!! The relevant part of the log file and the parameter file are pasted below.

Best,
Jacky

relevant part of log file

– This is ASPECT, the Advanced Solver for Problems in Earth’s ConvecTion.
– . version 2.2.0-pre (master, 729c918e1)
– . using deal.II 9.0.0
– . with 32 bit indices and vectorization level 1 (128 bits)
– . using Trilinos 12.12.1
– . using p4est 2.0.0
– . running in DEBUG mode
– . running with 1 MPI process

Number of active cells: 768 (on 2 levels)
Number of degrees of freedom: 28,690 (20,790+970+6,930)

*** Timestep 0: t=0 years
Solving temperature system… 0 iterations.
Rebuilding Stokes preconditioner…
Solving Stokes system… 118+0 iterations.

Postprocessing:

 Model domain depth (m):                        2.89e+06
 Temperature contrast across model domain (K):  3000
 Reference depth (m):                           0
 Reference temperature (K):                     0
 Reference pressure (Pa):                       0
 Reference gravity (m/s^2):                     9.81
 Reference density (kg/m^3):                    3458.4
 Reference thermal expansion coefficient (1/K): 3e-05
 Reference specific heat capacity (J/(K*kg)):   1250
 Reference thermal conductivity (W/(m*K)):      4.125
 Reference viscosity (Pa*s):                    7.72399e+21
 Reference thermal diffusivity (m^2/s):         9.54198e-07
 Rayleigh number:                               1e+07

[…]

Parameter file

A simple setup for convection in a 3d shell. See the

manual for more information.

set Dimension = 3
set Use years in output instead of seconds = true
set End time = 0
set Output directory = output-shell_simple_3d

subsection Material model
set Model name = simple

subsection Simple model
set Reference density = 3300
set Viscosity = 7.72399e21
set Thermal expansion coefficient = 3e-5
set Reference temperature = 1600
set Thermal conductivity = 4.125
set Reference specific heat = 1250
end
end

subsection Geometry model
set Model name = spherical shell

subsection Spherical shell
set Inner radius = 3481000
set Outer radius = 6371000
end
end

subsection Boundary velocity model
set Zero velocity boundary indicators = bottom
set Tangential velocity boundary indicators = top
end

subsection Boundary temperature model
set Fixed temperature boundary indicators = top, bottom
set List of model names = spherical constant

subsection Spherical constant
set Inner temperature = 3300
set Outer temperature = 300
end
end

subsection Initial temperature model
set Model name = function
subsection Function
set Function expression = 1.473e3
end
end

subsection Gravity model
set Model name = radial constant
subsection Radial constant
set Magnitude = 9.81
end
end

subsection Mesh refinement
set Initial global refinement = 1
set Initial adaptive refinement = 0
set Strategy = temperature
set Time steps between mesh refinement = 15
end

subsection Postprocess
set List of postprocessors = basic statistics
end

2

Hi Jacky,

Yes attaching files is not allowed on the first (or first 3?) posts, but I manually allowed it for your next posts. It is a spam filter thing.

I think you ran into an inconsistency in the input parameter file (that exists in a slightly different form in shell_simple_3d.prm as well). The basic statistics postprocessor that computes the Rayleigh number and all the other output does not actually compute the average temperature in the model. Instead it uses the adiabatic reference profile as a representative value for the model. However, in shell_simple the adiabatic surface temperature is not set, and the default is set to 0. Thus, when computing the density it computes it for a T of 0 K, which leads to the bigger density. Try setting the following outside of any subsection:

set Adiabatic surface temperature = 1600

Personally, I find it confusing that there are two different reference temperatures (the material model should just use the adiabatic reference profile), but removing the parameter from the material model might break a lot of .prm files.

Hope that helps,
Rene

3

Hi Rene,

That makes sense, thanks! One follow-up question: I assume this issue only affects the postprocessor and not the actual calculation, is that right? Or is the adiabatic surface temperature used somewhere else in an incompressible calculation?
Would it make sense to set the default Adiabatic surface temperature to be the reference temperature in case an adiabatic surface temperature is not specified in order to correctly calculate the Rayleigh number here?

Thanks!

Best,
Jacky

4

Hi Jacky,

I think under incompressibility the adiabatic surface temperature is also used to set up the reference density profile (though constant throughout the mantle from the simple material model) and will be used as the density term at the LHS of the energy equation. I was once also confused on this set-up a little bit. The density term in momentum equation and energy equation are set up from the different places under the incompressibility. I will let Rene to confirm this further:)

In my model, I set the adiabatic surface temperature same as my reference temperature in incompressible model.

Best,

Shangxin

5

Yes, the ‘Simple’ material model was the first one we ever implemented,
and it predates the availability of the adiabatic temperature profile.

There is of course also the issue that the “basic statistics”
postprocessor really only works for the ‘Simple’ material model because
it makes assumptions about how the material model is parameterized. I’ve
previously thought that we should just remove it because it doesn’t (and
can’t) produce useful information for more realistic situation.

@Jacky: Out of curiosity, what kind of information were you trying to
get out of this postprocessor?

Best
W.

6

Thanks Shangxin for your response! I understand.
Wolfgang, I’m actually not sure exactly why I had the basic statistics postprocessor included (maybe because I copied things over from a cookbook file in the first place). I wasn’t really using it to gain information but when I looked at the log output I realized the difference in Ra number and just wanted to make sure I understood where this comes.
Thanks everyone!

Best,
Jacky