A possible tiny bug in the temperature field initialization of the 3D spherical shell? #1443
Comments
|
Hi Shangxin, |
|
@gassmoeller 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) |
|
Hi Shangxin, |
|
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 |
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:
The calls to the material model for computing the reference profile are not used for assembling the matrix.
Yes.
Well you are right this situation is unfortunate. There are two different temperatures for historic reasons ( |
|
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/initial_profile.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/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? 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/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? 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/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? Sorry for the long message but I think I'm closer than before and only several more points to confirm with you:) |
|
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.
|
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 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.
The way we intended this feature to work is that the |
|
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): 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. 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:-) |
|
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.
|
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
The text was updated successfully, but these errors were encountered: