A possible tiny bug in the temperature field initialization of the 3D spherical shell?

[Shangxin-Liu]

Skipping the process how I detect this (accidentally), this possible tiny bug is described here:

When the 3D spherical shell geometry is used (I haven't tested this in 2D spherical shell), the temperature field at all the quadrature points along the north half polar axis is always 0. In the verification test, I fix the temperature field everywhere constantly to 0.5, including top (radius 1)/ bottom (radius 0.55) boundaries in the prm file, and run only one 0 time step.

I add the printout command to simple.cc file to print four columns data of every quadrature point into the log file and I found that all the quadrature points along the north half polar axis have 0 temperature values, which is not correct. The data output is like this:

x y z in.temperature[i]

0.000000 0.000000 1.000000 0.000000 0.000000 0.000000 0.999550 0.000000 0.000000 0.000000 0.999099 0.000000 0.000000 0.000000 0.998649 0.000000 0.000000 0.000000 0.998198 0.000000 0.000000 0.000000 0.997748 0.000000 ............ 0.000000 0.000000 0.801351 0.000000 0.000000 0.000000 0.800901 0.000000 0.000000 0.000000 0.800450 0.000000 0.000000 0.000000 0.800000 0.000000 0.000000 0.000000 0.799550 0.000000 ............ 0.000000 0.000000 0.552703 0.000000 0.000000 0.000000 0.552252 0.000000 0.000000 0.000000 0.551802 0.000000 0.000000 0.000000 0.551351 0.000000 0.000000 0.000000 0.550901 0.000000 0.000000 0.000000 0.550450 0.000000 0.000000 0.000000 0.550000 0.000000 ............

And all the other quadrature points have the correct constant 0.5 temperature values.

The bug seems due to an omission while assigning the input temperature field to each quadrature point within the 3D spherical shell domain. But for now I don't know where in the code assigns the temperature field to the geometry domain.

Because only the north half polar axis has the ghost zero temperature and hence wrong density values, I'm kind of doubting this can significantly change the solution of 3D spherical shell but it's still worth fixing this possible tiny bug:-)

Also, in my test of the time-dependent cases, this tiny issue only happens in the first time step when assigning the temperature field to the 3D spherical shell domain and doesn't appear in all the following time steps.

So how do I fix this issue? @bangerth @gassmoeller @tjhei @jdannberg

Best, Shangxin

[gassmoeller]

Hi Shangxin, the assignment of initial temperature to the solution happens in simulator/initial_conditions.cc:98 if you use a deal.II newer than 8.4 and in simulator/initial_conditions.cc:166 if you use deal.II 8.4 (this will be removed in the future). However, all of this is geometry independent, the mesh itself does not know if it is at a polar axis or not (it does not even know if it is a sphere or not, it is just a collection of coordinates). I would therefore assume the bug is inside the initial condition plugin, maybe that one does not calculate temperature correctly. Did you use any of the existing plugins, and which one? Do you have a simple parameter file that reproduces the problem?

[Shangxin-Liu]

@gassmoeller Hi Rene,

The attachment is the prm file I use to produce this issue. My observation is that the initial condition subsection does not matter. I've tried function, harmonic_perturbation, S40RTS, SAVANI, etc.... And as long as I use 3D spherical shell, this issue always happens. The prm file attached here uses function initial condition to set the temperature field to constant 0.5 everywhere and the boundary condition also set 0.5 to both top and bottom. So I assume the output temperature field would be 0.5 everywhere but the quadrature points along the north half polar axis are not, they are all 0. Due to only a single axis, I suppose this issue cannot be seen in paraview so hasn't been detected previously.

I know you have your own way to print out the temperature values of quadrature points. But to confirm, I use these two lines of the code in the beginning of simple.cc subroutine to print out the ghost temperature values along the polar axis:

for (unsigned int i=0; i < in.position.size(); ++i) { const Point position = in.position[i]; if (in.temperature[i] == 0.) printf("%f %f %f %f\n",position[0],position[1],position[2],in.temperature[i]);

a1.prm.txt

[gassmoeller]

Hi Shangxin, I checked your input file and could reproduce the issue. However, it is slightly different than you think. The initial condition is applied correctly to all quadrature points of the geometry, but the material model's evaluate function is also called for computing the reference temperature/pressure/density profile that is necessary for the boussinesq approximation (that happens from adiabatic_conditions/initial_profile.cc), and since the Adiabatic surface temperature of your model is not specified it is computed for a default of 0. For the reference profile we need to assume some positions for provide in.position for the material model. This is done by calling the Geometry models representative_point function, which in the case of the spherical shell simply provides positions along the positive z-axis with the corresponding depth. So that is why you see calls to the material model for positions along the positive z-axis with an input temperature of 0, these material output are only used for the reference profile, not for setting the initial conditions.

[Shangxin-Liu]

Hi Rene,

Thanks for the detailed explanation. So in the simple.cc file, I noticed that after my printout lines of code, in.temperature[i] will be used to subtract the reference temperature to get the temperature anomaly. And this temperature anomaly will further go into the temperature-dependent viscosity equation and density buoyancy equation. So do you mean that although this 0 positive z-axis reference profile will still go into the later temperature-dependent viscosity and density equations, the viscosity and density from this 0 reference temperature profile will be identified and will not be put into the Stokes and energy equations' solver matrix?

And if I set the parameter Adiabatic surface temperature also to 0.5, the values of this positive z-axis reference temperature profile will also be 0.5, yes?

I have one more question, what's the difference of this reference temperature profile along positive z-axis and the input parameter Reference temperature of simple model? In the case of incompressibility Boussinesq approximation, we will use a constant reference temperature throughout mantle. So what's the point of using one more reference temperature profile in the code? Is this common for every material model in ASPECT? @gassmoeller

[gassmoeller]

Thanks for the detailed explanation. So in the simple.cc file, I noticed that after my printout lines of code, in.temperature[i] will be used to subtract the reference temperature to get the temperature anomaly. And this temperature anomaly will further go into the temperature-dependent viscosity equation and density buoyancy equation. So do you mean that although this 0 positive z-axis reference profile will still go into the later temperature-dependent viscosity and density equations, the viscosity and density from this 0 reference temperature profile will be identified and will not be put into the Stokes and energy equations' solver matrix?

I believe you think about the process the wrong way around. Let me try to draw a little flowchart to visualize the relationship between the objects:

1. Initialization of model (once in the beginning) -> compute adiabatic reference profile -> call material model with reference position/temperature/pressure (ignore the viscosity output, because that one is not needed in the reference profile)

2. Assembly of Energy matrix (once per timestep) -> call material model with actual position/temperature/pressure (but for boussinesq approximation replace all occurences of the density in the assembly with the reference density stored in the adiabatic reference profile computed in 1.)

3. Assembly of Stokes matrix (once per timestep) -> call material model with actual position/temperature/pressure (use all of the material outputs for the matrix assembly)

The calls to the material model for computing the reference profile are not used for assembling the matrix.

And if I set the parameter Adiabatic surface temperature also to 0.5, the values of this positive z-axis reference temperature profile will also be 0.5, yes?

Yes.

I have one more question, what's the difference of this reference temperature profile along positive z-axis and the input parameter Reference temperature of simple model? In the case of incompressibility Boussinesq approximation, we will use a constant reference temperature throughout mantle. So what's the point of using one more reference temperature profile in the code? Is this common for every material model in ASPECT?

Well you are right this situation is unfortunate. There are two different temperatures for historic reasons (simple was created before we had put much work in the reference profile), but you are right, it can lead to confusion. I can see situations in which having two temperatures can be useful, e.g. if you have material properties measured in a laboratory at certain pressure-temperature conditions and want to use them in the material model, but want to use a different temperature for your reference profile. However this is definitely not the usual case, and maybe we should simply remove the reference temperature from the simple material model (if somebody wants to have the special case he can allows write a new material model for it). Do you want to give this a try at the hackathon?

[Shangxin-Liu]

Hi Rene, @gassmoeller

Thanks for a detailed workflow explanation. This wrongly reported issue makes me learn a lot:) Yes. I can surely give that a try at the hackathon. But there is a good news that I have successfully benchmarked my geoid code recently (still traditional spherical harmonic calculations but working with AMR). The numerical solution of my code is the same as the analytic solution calculated by hand in the benchmark cases. The benchmark process is somewhat complicated (containing two parts, density contribution and dynamic topography contribution) and I will briefly document that and bring it to hackathon for our discussion. I think this time is the point that we should put geoid postprocessor into ASPECT.

Based on your explanation here and also Wolfgang/Timo's on the mailing list, I carefully went though the code in the _adiabatic_conditions/initialprofile.cc, and I think now I am close to the correct concepts. Maybe this is also related to our time-dependent 3D spherical shell benchmark struggles. To make sure now I understand the whole process correctly, I'd like to confirm with you again about my concepts. Feel free to point out if I still think anything wrong:-)

From the file _adiabatic_conditions/initialprofile.cc, I noticed that the code integrates the reference temperature downward from the input adiabatic surface temperature and if it's incompressible (like our BA cases), all the temperature values along the reference profile will be adiabatic surface temperature constantly; if it's compressible, an adiabatic temperature profile with gradient will be computed. That's why in my BA case, the default 0 adiabatic surface temperature will results in all 0 reference temperature profile, yes? And in the initial temperature module, let's say S40RTS/SAVANI, the background temperature will be either this initial temperature profile if the model is compressible, or the input Reference temperature of the initial temperature subsection if the model is incompressible (although a constant initial temperature profile in BA is still calculated but is not used), yes?

For the initial pressure profile, the hydrostatic pressure profile will be used to evaluate the material model in Initialization of model (ignoring the viscosity) and also be used to evaluate the material model in the first time step's matrix assembly before the actual full pressure is firstly solved, yes?

As for the calculation of the initial reference density profile, I'm still a bit confused. In the file _adiabatic_conditions/initialprofile.cc, I noticed that every density value along the initial density profile is calculated from a loop evaluating the density by the material model (simple.cc in my case) with every representative_point, temperature, and pressure along the initial profile. I have two questions about the usage of the first quadrature point in.position[0] as the bridge. Firstly, when the loop index i is 1, since the material property hasn't been evaluated, where does the out.densities[0] in the line 89 come from? Secondly, while the first quadrature point in.position[0] is used as the bridge to evaluate reference density profile from material model, how about other quadrature points in this process with 1,2,3, etc index? Does the code just ignore them because here we only care about the the first quadrature point in.position[0] as the bridge? Maybe I still think things wrong here.

The last thing I want to confirm is the density term in the assembly of Energy matrix. From your explanation, commonly, this density term is from the material model. However, for the BA cases (Timo also mentioned this on mailing list), you mention that ASEPCT will replace all occurrences of the density in the assembly with the reference density stored in the adiabatic reference profile computed in the Initialization of model. However, based on this process, I have a concern about the usage of dynamic pressure instead of full pressure. When I use a material model in which the density has no constant term and only contains perturbation term to make the full pressure equals the dynamic pressure, in the compressible cases, it seems that the density term in the assembly of Energy matrix will also be only the density perturbation, which is not correct; And in my BA cases, although the density term comes from _adiabatic_conditions/initialprofile.cc, because this initial reference density profile still comes from the evaluate function in the material model (simple.cc in our cases), it will still be only the density perturbation form, which makes me a little worried, although we can set a value of adiabatic surface temperature to make this reference density profile to a constant value we want. Do I still understand things wrong here?

Sorry for the long message but I think I'm closer than before and only several more points to confirm with you:)

[bangerth]

our discussion. I think this time is the point that we should put geoid postprocessor into ASPECT.

Great, we look forward to that!

From the file /adiabatic_conditions/initial_profile.cc/, I noticed that the code integrates the reference temperature downward from the input adiabatic surface temperature and if it's incompressible (like our BA cases), all the temperature values along the reference profile will be adiabatic surface temperature constantly; if it's compressible, an adiabatic temperature profile with gradient will be computed. That's why in my BA case, the default 0 adiabatic surface temperature will results in all 0 reference temperature profile, yes? And in the initial temperature module, let's say S40RTS/SAVANI, the background temperature will be either this initial temperature profile if the model is compressible, or the input Reference temperature of the initial temperature subsection if the model is incompressible (although a constant initial temperature profile in BA is still calculated but is not used), yes?

I think this is correct.

For the initial pressure profile, the hydrostatic pressure profile will be used to /evaluate/ the material model in Initialization of model (ignoring the viscosity) and also be used to /evaluate/ the material model in the first time step's matrix assembly before the actual full pressure is firstly solved, yes?

Yes, correct.

As for the calculation of the initial reference density profile, I'm still a bit confused. In the file /adiabatic_conditions/initial_profile.cc/, I noticed that every density value along the initial density profile is calculated from a loop evaluating the density by the material model (/simple.cc/ in my case) with every representative_point, temperature, and pressure along the initial profile. I have two questions about the usage of the first quadrature point /in.position[0]/ as the bridge. Firstly, when the loop index /i/ is 1, since the material property hasn't been evaluated, where does the out.densities[0] in the line /89/ come from?

That was computed in line 127 in the previous iteration (i=0).

Secondly, while the first quadrature point /in.position[0]/ is used as the bridge to /evaluate/ reference density profile from material model, how about other quadrature points in this process with 1,2,3, etc index? Does the code just ignore them because here we only care about the the first quadrature point /in.position[0]/ as the bridge? Maybe I still think things wrong here.

The loop here walks down the reference profile, and for each depth, we only set up a structure with a single quadrature point (namely a reference point at the current depth).

You are confused by the fact that in.position is an array. It is an array because it is usually used to evaluate material properties at the quadrature points of cells, but here, it is only an array of length one. (See the initialization in line 56.)

The last thing I want to confirm is the density term in the assembly of Energy matrix. From your explanation, commonly, this density term is from the material model. However, for the BA cases (Timo also mentioned this on mailing list), you mention that ASEPCT will replace all occurrences of the density in the assembly with the reference density stored in the adiabatic reference profile computed in the Initialization of model.

I believe that this is correct -- you should be able to find this in the file source/simulator/assemblers/stokes.cc for the various formulations.

However, based on this process, I have a concern about the usage of dynamic pressure instead of full pressure. When I use a material model in which the density has no constant term and only contains perturbation term to make the full pressure equals the dynamic pressure, in the compressible cases, it seems that the density term in the assembly of Energy matrix will also be only the density perturbation, which is not correct;

That depends on the formulation you use. If you use the fully compressible formulation, then that would make no sense because the density that appears in the rhs of the Stokes equation is the same one that is used in the energy equation. Clearly, in that case the density better not be negative.

On the other hand, if you use ALA/TALA/BA, then the densities that appear in these two places are unrelated and can be chosen independently.

And in my BA cases, although the density term comes from /adiabatic_conditions/initial_profile.cc/, because this initial reference density profile still comes from the /evaluate/ function in the material model (/simple.cc/ in our cases), it will still be only the density perturbation form, which makes me a little worried, although we can set a value of adiabatic surface temperature to make this reference density profile to a constant value we want. Do I still understand things wrong here?

Then you haven't set things up correctly. For the BA, you can (and need to) independently choose the reference profile and the density in the material model. They can not be the same.

[jdannberg]

From the file /adiabatic_conditions/initial_profile.cc/, I noticed that the code integrates the reference temperature downward from the input adiabatic surface temperature and if it's incompressible (like our BA cases), all the temperature values along the reference profile will be adiabatic surface temperature constantly; if it's compressible, an adiabatic temperature profile with gradient will be computed. That's why in my BA case, the default 0 adiabatic surface temperature will results in all 0 reference temperature profile, yes? And in the initial temperature module, let's say S40RTS/SAVANI, the background temperature will be either this initial temperature profile if the model is compressible, or the input Reference temperature of the initial temperature subsection if the model is incompressible (although a constant initial temperature profile in BA is still calculated but is not used), yes?

I think this is correct.

This is almost correct. Whether the reference temperature is constant or increases with depth depends on the question if you use adiabatic heating in your model or not (see lines 101-105 in source/adiabatic_conditions/initial_profile.cc). As you would usually include adiabatic heating in a compressible model, in many cases this is the same thing. Also, if you start with a surface temperature of zero, your profile will always be zero all the way down to the bottom of your model (as the slope of the adiabat depends on the temperature, and is zero when the temperature is zero).

The background temperature in your initial temperature profile depends on the plugin you use. In for example the S40RTS model, it indeed depends on the model being compressible or not if a constant value or the adiabatic profile is used, but that might be defined differently by other plugins.

And in my BA cases, although the density term comes from /adiabatic_conditions/initial_profile.cc/, because this initial reference density profile still comes from the /evaluate/ function in the material model (/simple.cc/ in our cases), it will still be only the density perturbation form, which makes me a little worried, although we can set a value of adiabatic surface temperature to make this reference density profile to a constant value we want. Do I still understand things wrong here?

The way we intended this feature to work is that the initial_profile adiabatic conditions plugin automatically gives you the adiabatic profile you would want to use in most of the cases. However, using only a density perturbation for the density is not a case this would work for. Hence, there is the option to use other plugins to define your adiabatic profile, such as the function or the ascii data plugin. In your case, I think it would be the easiest way to use function, and then you can define your reference temperature, pressure and density profiles as anything you want. You would do this by going to the subsection Adiabatic conditions model in the input file, and then setting the Model name to function.

[Shangxin-Liu]

Hi Wolfang, Rene, and Juliane, @bangerth @gassmoeller @jdannberg

Thanks a lot for all the detailed descriptions. I've figured out the infrastructure of the code by your instruction:) Now only for the last dynamic pressure issue, I think in my last message I may not describe my concern/question clearly. Let me confirm myself again here. As Wolfgang pointed out, for ALA/TALA/BA, the densities that appear in Stokes matrix and Energy matrix are unrelated and can be chosen independently. But when I use initial profile in Adiabatic conditions model and set the density in material model (simple.cc in my case) to only density perturbation like this form (the reference_T in simple.cc is set to 0): out.densities[i] = reference_rho (- thermal_alpha in.temperature[i]) Although in ALA/TALA/BA cases, the density term in energy matrix is unrelated to the density term in Stokes matrix, because the density term in the energy matrix comes from _adiabatic_conditions/initialprofile.cc where the reference density profile is still evaluated from the material model by _this->get_materialmodel().evaluate(in, out) at line 127 in initial_profile.cc, i.e., the above equation by set in.temperature[i] to the reference temperature profile at each depth, if we use a previous zero reference temperature profile due to the default zero adiabatic surface temperature, the reference density profile will also be zero, which is still not what we want, since we'd like to set the density term in the energy equation to constant 1, same as non-dimensional equation from Zhong 2008. Do I understand things correctly here? And is this why Juliane pointed out that using only a density perturbation for the density in material model is not a case initial_profile would work for so I should use function to set my own constant 1 reference density profile? @jdannberg

Also, in our benchmark prm file we set reference_rho and thermal_alpha both to 1. To make the reference density profile to constantly 1, it seems that if I set the adiabatic surface temperature to -1 to get the -1 constant reference temperature profile (no adiabatic heating), we indeed can get a constant 1 reference density profile from the above equation, although I'm not sure whether ASPECT allow the users to set a negative -1 adiabatic surface temperature. This seems provide a possibility to still use the default _initialprofile.cc adiabatic conditions plugin in our dynamic pressure cases.

Considering the -1 adiabatic surface temperature looks weird and as Juliane suggested, we'd better use function to prescribe our reference temperature, pressure and density profiles in Adiabatic conditions model. Here I'd like to describe our settings a little bit to confirm further. Because we are using BA with no adiabatic heating, our initial temperature model, harmonic perturbation will directly use the input Reference temperature and our material model is simple.cc where the material properties don't depend on pressure. So does this mean that in our benchmark cases, the prescribed reference temperature and pressure profile will not be used hence it doesn't matter what kind of temperature and pressure profile we set in function subsection? (i.e, we can simply set them both to the default 0.?) And the only reference profile we need to set here is the density where in our cases it is constantly 1, yes?

My last question is about the pressure normalization. Previously when we use the common total pressure, we always set pressure normalization to surface with 0 surface average pressure, the default setting. Now when we use the dynamic pressure instead, is there any special care we need in pressure normalization? Or we can still do the same as before, setting the pressure normalization to surface with 0 surface average pressure?

I believe the above are the last bunch of stuff I want to confirm further:-)

[bangerth]

Thanks a lot for all the detailed descriptions. I've figured out the infrastructure of the code by your instruction:) Now only for the last dynamic pressure issue, I think in my last message I may not describe my concern/question clearly. Let me confirm myself again here. As Wolfgang pointed out, for ALA/TALA/BA, the densities that appear in Stokes matrix and Energy matrix are unrelated and can be chosen independently. But when I use /initial profile/ in /Adiabatic conditions model/ and set the density in material model (/simple.cc/ in my case) to only density perturbation like this form out.densities[i] = reference_rho (- thermal_alpha in.temperature[i]) Although in ALA/TALA/BA cases, the density term in energy matrix is unrelated to the density term in Stokes matrix, because the density term in the energy matrix comes from /adiabatic_conditions/initial_profile.cc/ where the reference density profile is still evaluated from the material model by /this->get_material_model().evaluate(in, out)/ at line 127 in initial_profile.cc, i.e., the above equation by set in.temperature[i] to the reference temperature profile at each depth, if we use a previous zero reference temperature profile due to the default zero adiabatic surface temperature, the reference density profile will also be zero, which is still not what we want, since we'd like to set the density term in the energy equation to constant 1, same as non-dimensional equation from Zhong

2008. Do I understand things correctly here?

Yes, that is correct. It makes no sense to choose the adiabatic conditions in the ALA/TALA/BA from the same material model. It really does need to be independent.

It is independent, for example, if you do this:

Also, in our benchmark prm file we set reference_rho and thermal_alpha both to 1. To make the reference density profile to constantly 1, it seems that if I set the adiabatic surface temperature to -1 to get the -1 constant reference temperature profile (no adiabatic heating), we indeed can get a constant 1 reference density profile from the above equation, although I'm not sure whether ASPECT allow the users to set a negative -1 adiabatic surface temperature. This seems provide a possibility to still use the default /initial_profile.cc/ adiabatic conditions plugin in our dynamic pressure cases.

Yes. If you make the surface temperature and the adiabatic surface temperature different, then you achieve your goal. Choosing -1 is not wrong -- it's just another value, and in at least some of the approximations, the 'temperature' is only a deviation from a reference temperature anyway, so -1 is not an impossible choice.

Now, as you already suggest, just because this approach works doesn't mean that it's the best way to do it.

Considering the -1 adiabatic surface temperature looks weird and as Juliane suggested, we'd better use function to prescribe our reference temperature, pressure and density profiles in Adiabatic conditions model. Here I'd like to describe our settings a little bit to confirm further. Because we are using BA with no adiabatic heating, our initial temperature model, harmonic perturbation will directly use the input Reference temperature and our material model is simple.cc where the material properties don't depend on pressure. So does this mean that in our benchmark cases, the prescribed reference temperature and pressure profile will not be used hence it doesn't matter what kind of temperature and pressure profile we set in function subsection? (i.e, we can simply set them both to the default 0.?) And the only reference profile we need to set here is the density where in our cases it is constantly 1, yes?

That I don't know and would rather let others answer.