[aspect-devel] Internal heating in aspect (Ludovic Jeanniot)

Tue Aug 28 08:59:19 PDT 2018

All,
Thank you for your help with these issues. To summarize the discussion so
far, it seems that four ideas have been suggested, in order of increasing
complexity

0) Try using the DG elements for the energy equation, as no stabilization
is needed
     I tried this and ran into failures to converge within the first 100
timesteps. I suspect that I need to adjust solver parameters but have not
done so yet.
1) Increase global resolution
     Currently running, but will take some time.
2) Disable the stabilization of the temperature solution by setting cR
and/or beta=0
    This is an intriguing idea. I have a calculation in progress without
stabilization. Note that there is also some difference in the temperature
evolution because I also set a maximum timestep. By default, the
shell_simple_3d cookbook contained no initial perturbation so the initial
time step was very large in the calculation shown in black at GR3 with
default stabilization parameters.
[image: image.png]

3) Adjust the mesh to achieve more uniform resolution in the radial
direction
4) Compute the boundary heat flux using CBF

I am in the process of trying each of these things. Some will take quite a
bit of time. The global models that we were running when we identified the
discrepancy between thermal evolution predicted by ASPECT and CitcomS were
run at either global refinement 4 without any adaptive refinement or with
global refinement 3 and 2 levels of adaptive refinement. I can provide full
.prm files but I don't think that it's relevant to the immediate issue,
which is the apparent lack of energy conservation in much simpler, purely
thermal convection calculations.

Before moving forward with higher resolution thermal convection
calculations, I thought that it would be informative to look at a simpler
problem. I set up an ASPECT calculation to model the cooling of a spherical
shell without convection. To achieve this, I disabled gravity in the input
file for shell_simple_3d.prm. I also enabled the 'conduction timestep' to
prevent very large timesteps from compromising the accuracy of the
solution. The models were run at global refinement levels 3, 4, and 5. I
verified that the RMS velocities are zero everywhere. The spherical shell
is subject to constant temperature boundary conditions at the surface and
CMB. The initial condition is uniformly equal to the CMB temperature. There
is an analytic solution to the steady conduction problem. I ran the
calculations for 200 Gyr, long enough for the temperature field and heat
fluxes to reach steady state. The heat flow at the surface and CMB should
be equal at steady state to conserve energy. In this problem, the
temperature gradient is very small, so I would expect the agreement between
the analytic steady state solution and the numerical solution to be
exceedingly good, even at relatively coarse resolutions. What I found is
actually very slow convergence towards reasonable heat flow values and
significant deviation from the analytic solution. There is also a very
significantly different rate of evolution towards the analytic solution,
seen clearly in the time evolution of the surface and CMB heat flow. There
is also a very large mismatch between the surface and CMB heat flows. I
would be very very surprised if such a large disagreement could be
explained by not using the CBF method for computing the heat flux. The
thermal gradients in this problem are on the order of mK per meter.

[image: image.png]
The second figure below shows the temperature profile (computed using the
depth average postprocessor), along with the analytic solution.

[image: image.png]

Clearly, I must be doing something wrong, because this is one of the
simplest problems to solve numerically. Does anyone have an idea what the
problem may be? I have attached the .prm file (set up for the highest
resolution run shown here). Based on what I am seeing in the calculation
with no stabilization, I am wondering whether the default stabilization
parameters result in a far too diffusive solution. I will repeat the
cooling sphere calculations next without stabilization and see how it
compares with the analytic solution.

Best
Max

On Mon, Aug 27, 2018 at 9:36 AM Scott King <sking07 at vt.edu> wrote:

> Thanks.  We will check this out.
>
> We did not change the cells along the circumference but changed the radial
> distribution be uniform.  I think it required a code modification not a
> parameter change but I could be wrong.
>
> Scott
>
> Sent from my iPhone
>
> > On Aug 26, 2018, at 11:03 AM, Jeanniot, L. (Ludovic) <l.jeanniot at uu.nl>
> wrote:
> >
> > Rene, Scott, Ian and Max,
> >
> > I am currently benchmarking thermal convection in a 3D hollow sphere
> following Arrial et al. 2014 setup to compare the obtained results (average
> temperature, V_rms, Nusselt number and thermal shape at steady-state) in
> ASPECT with CitcomS and a RBF-PS model (Arrial et al. 2014). The whole
> setup is made to be nondimensional.
> > I did encountered similar issues of lowered average temperature than
> expected at low resolution (GMR = 3 (level 4)) and found a certain solution
> although I still need to understand why now it works.
> >
> >     The first issue is about the default mesh:
> > 1.- The default mesh is (I think) depth-dependent as the mesh is finer
> beneath a certain depth and coarser above. What is happening is that at low
> resolution, the temperature solution at the boundary between the two
> finer/coarser meshes is misbehaving (under-solved?). In CitcomS, the mesh
> is composed of about 12 'caps' of identical number of cells and I therefore
> changed the parameter "Cells along circumference" to 12 in the geometry
> model subsection to reproduce a similar and more uniform grid. But I must
> admit that what Cells along circumference does remain obscure to me.
> >
> >     The second issue may be what you are looking for:
> > 2.- I attached some experiments with changing the sensitivity of the
> beta factor in the entropy viscosity stabilization and the cR factor in the
> artificial viscosity stabilization (in the manual: A.41 Parameters in
> section Discretization/Stabilization parameters). By default in ASPECT, the
> values of the two factors are 0.078 and 0.33 respectively. In the figure,
> the 3 blue lines, which were obtained when turning the two factor values to
> 0, or one of the two to 0, are identical and they correspond at what I
> should get (note the level 4 resolution). Keeping at low resolution, by
> default with two factor values beta = 0.078 and cR = 0.33 results in
> dropping the average temperature from 0.218 to 0.19. Running the models
> with higher resolution then results in the average temperature I was
> expecting (similar to CitcomS and RBF-PS model)
> >
> > Maybe this can help too?
> >
> > Ludovic
> >
> >
> > -----Original Message-----
> > From: Aspect-devel [mailto:aspect-devel-bounces at geodynamics.org] On
> Behalf Of aspect-devel-request at geodynamics.org
> > Sent: 25 August 2018 22:17
> > To: aspect-devel at geodynamics.org
> > Subject: Aspect-devel Digest, Vol 81, Issue 7
> >
> > Send Aspect-devel mailing list submissions to
> >    aspect-devel at geodynamics.org
> >
> > To subscribe or unsubscribe via the World Wide Web, visit
> >    http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
> > or, via email, send a message with subject or body 'help' to
> >    aspect-devel-request at geodynamics.org
> >
> > You can reach the person managing the list at
> >    aspect-devel-owner at geodynamics.org
> >
> > When replying, please edit your Subject line so it is more specific than
> "Re: Contents of Aspect-devel digest..."
> >
> >
> > Today's Topics:
> >
> >   1. Re: Internal heating in aspect (Scott King)
> >   2. Re: Internal heating in aspect (Ian Rose)
> >
> >
> > ----------------------------------------------------------------------
> >
> > Message: 1
> > Date: Sat, 25 Aug 2018 15:50:25 -0400
> > From: Scott King <sdk at vt.edu>
> > To: aspect-devel <aspect-devel at geodynamics.org>
> > Subject: Re: [aspect-devel] Internal heating in aspect
> > Message-ID: <1E881A37-48D1-46FA-A2F1-9FEF1FAC8308 at vt.edu>
> > Content-Type: text/plain; charset="utf-8"
> >
> > Hi Rene, Max, Wolfgang, and Diogo;
> >
> > Grant found when looking at the Zhong 2008 “Benchmark” calculations that
> he gets better agreement with mean temperature and mean RMS velocity than
> he does with surf/cmb Nusselt number.  We find generally better results
> with a nearly uniform radial mesh than the aspect default mesh pretty much
> across the board.  As we go to higher and higher resolution (5 or 6
> refinements), we get convergence of the surf and cmb heatflow and good
> agreement with other codes.  We are puzzling about it because the mean
> temperature and mean velocity profiles (fn of r) are spot on with the ones
> from CitcomS, so the fluxes should really be correct.  The numerical
> diffusion still being included sounds like a possibility. On the other
> hand, that explanation puzzles me because the 2D Cartesian cases are spot
> on.  It made me wonder if there is a scaling issue that only shows up in
> the spherical case?  On my list of things to do is to check into that but
> life interrupted my to do list.
> >
> > Interestingly we get a significantly higher mean temperatures
> (essentially the profile shifted by 10%) for the C cases (Ra=10^5) and
> we’ve never quite figured out why.  That sort of concerns me because all
> the recent spherical papers I’ve seen only address the Ra=7000 cases and
> assume this tell us enough about how these codes behave at higher Ra, which
> seems to me a bit optimistic.   We’ve also never seen any codes try to
> reproduce the 10^5 results or to use more than a single grid either :)   Of
> course that makes me wonder which code is correct. (LOL)  Convergence
> studies seem to be out of favor.
> >
> > Scott
> >
> >> On Aug 25, 2018, at 1:58 PM, Rene Gassmoeller <
> rene.gassmoeller at mailbox.org> wrote:
> >>
> >> Hi Max,
> >>
> >> let me cc this to the mailing list, as it might be of interest to
> others.
> >> I found this old email thread on the mailing list that seems to discuss
> the same issue:
> http://lists.geodynamics.org/pipermail/aspect-devel/2012-November/000138.html
> <
> http://lists.geodynamics.org/pipermail/aspect-devel/2012-November/000138.html>
> (there are some additional messages when you scroll through the thread). In
> short: At the time Ian found that this is a problem when the boundary layer
> is underresolved. Since the element at the surface are larger than the
> elements at the CMB they do resolve a different temperature gradient for
> the boundary layer if the boundary layer is not sufficiently resolved, and
> we compute the heat flux based on the temperature gradient.
> >> However, thinking about it again, I do not think this can be the only
> reason. If the difference were caused by numerical resolution alone, the
> average temperature should change, until the two heat fluxes are in balance
> as you mentioned. But another idea I had was that the computed heat flux is
> a postprocessor that only considers physical diffusion (i.e. diffusion due
> to thermal conductivity), while in the assembly we also include the
> numerical diffusion to stabilize the equation. So maybe the difference
> between the heat fluxes is caused by numerical diffusion (and it is so
> large because the cells are relatively coarse)? If that were the case the
> difference should increase/decrease with higher/lower resolution, as the
> stabilization scales with the element size. Unfortunately that is not a
> criterion to distinguish between the two possible causes (obviously
> resolution of the temperature gradient across the boundary layer also
> scales with element size). You could try to re-run the model with
> discontinuous temperature elements, that should disable the numerical
> diffusion. Let me know if that makes sense.
> >> Best,
> >>
> >> Rene
> >>
> >>> On 08/24/2018 06:46 PM, Max Rudolph wrote:
> >>> Rene and Wolfgang,
> >>> We're running some models in 3D spherical geometry in aspect with
> internal heating and are unable to reproduce thermal histories consistent
> with what we saw in identically configured calculations in CitcomS. It is
> possible that there is some aspect of our model setup that is inconsistent,
> but we are using the same depth profile of viscosity, same boundary
> conditions, and same initial conditions. The ASPECT models cool much faster
> (which also drives up viscosity). We are trying to get understand these
> results, and I was looking at the heat transport in the 3D spherical model
> discussed in section 5.3.2 of the aspect manual 'simple convection in the
> 3D spherical shell'. In figure 40, it appears that this calculation has
> been run to nearly statistically steady state. Is this correct? I am
> wondering, basically, how to interpret the imbalance between reported
> surface and CMB heat flow. To follow up on this, I ran the cookbook .prm
> file from 5.3.2 (shell_simple_3d) for 5 Gyr, which appears to be very close
> to or at statistically steady state at the end of the calculation. Due to
> resource constraints, I only ran at global refinement level 3 (definitely
> under-resolved, mea culpa). However, at the end of the calculation, there
> is still a significant imbalance between the surface and core mantle
> boundary heat flows. This is a model with no internal heating, so at
> statistically steady state, the surface and cmb heat flows should balance
> to conserve energy. Do you have any insight into why we might be seeing
> this behavior? I've attached the input .prm file, the statistics file, and
> a plot showing the average temperature and surface and cmb heat flows.
> Thank you for thinking about this!
> >>>
> >>> Best
> >>> Max
> >>>
> >>> p.s. I sent this to both of you because it looked to me like Wolfgang
> contributed the cookbook. Please let me know if there is someone else that
> might be better positioned to help with this.
> >>>
> >>> <image.png>
> >>>
> >>>
> >>>
> >>
> >> --
> >> Rene Gassmoeller
> >> https://gassmoeller.github.io/
> >> <https://gassmoeller.github.io/>______________________________________
> >> _________
> >> Aspect-devel mailing list
> >> Aspect-devel at geodynamics.org
> >> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
> >
> > -------------- next part --------------
> > An HTML attachment was scrubbed...
> > URL: <
> http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20180825/b72f5484/attachment-0001.html
> >
> >
> > ------------------------------
> >
> > Message: 2
> > Date: Sat, 25 Aug 2018 16:17:13 -0400
> > From: Ian Rose <ian.rose at berkeley.edu>
> > To: aspect-devel <aspect-devel at geodynamics.org>
> > Subject: Re: [aspect-devel] Internal heating in aspect
> > Message-ID:
> >    <CAEm=8TR0FX7wiMDYXGh7zWtVouQ0cgH1GYnJgNQzz7JPnjQfmg at mail.gmail.com>
> > Content-Type: text/plain; charset="utf-8"
> >
> > It is also possible that the accuracy could be improved by using the
> consistent-boundary-flux method (CBF) for calculating the fluxes in
> postprocessing. The CBF can be seen as a form of projection for the
> solution into the flux space. It essentially constructs a new linear
> system, the solution to which is the fluxes at the boundaries. That system
> can be made to be diagonal, so it is cheap to solve.
> >
> > I found that using CBF in calculating the dynamic topography
> significantly improved the convergence of the solution (quadratic rather
> than linear, if memory serves). The mathematics for computing the heat flux
> is essentially identical, so a lot of the dynamic topography code could be
> copied whole-cloth into the heat flux postprocessor, in case anyone wants
> to give it a shot.
> >
> > Best,
> > Ian
> >
> >> On Sat, Aug 25, 2018 at 4:01 PM Scott King <sdk at vt.edu> wrote:
> >>
> >> Hi Rene, Max, Wolfgang, and Diogo;
> >>
> >> Grant found when looking at the Zhong 2008 “Benchmark” calculations
> >> that he gets better agreement with mean temperature and mean RMS
> >> velocity than he does with surf/cmb Nusselt number.  We find generally
> >> better results with a nearly uniform radial mesh than the aspect
> >> default mesh pretty much across the board.  As we go to higher and
> >> higher resolution (5 or 6 refinements), we get convergence of the surf
> >> and cmb heatflow and good agreement with other codes.  We are puzzling
> >> about it because the mean temperature and mean velocity profiles (fn
> >> of r) are spot on with the ones from CitcomS, so the fluxes should
> >> really be correct.  The numerical diffusion still being included
> >> sounds like a possibility. On the other hand, that explanation puzzles
> >> me because the 2D Cartesian cases are spot on.  It made me wonder if
> >> there is a scaling issue that only shows up in the spherical case?  On
> >> my list of things to do is to check into that but life interrupted my
> to do list.
> >>
> >> Interestingly we get a significantly higher mean temperatures
> >> (essentially the profile shifted by 10%) for the C cases (Ra=10^5) and
> >> we’ve never quite figured out why.  That sort of concerns me because
> >> all the recent spherical papers I’ve seen only address the Ra=7000
> >> cases and assume this tell us enough about how these codes behave at
> higher Ra, which seems to me a bit
> >> optimistic.   We’ve also never seen any codes try to reproduce the 10^5
> >> results or to use more than a single grid either :)   Of course that
> makes
> >> me wonder which code is correct. (LOL)  Convergence studies seem to be
> >> out of favor.
> >>
> >> Scott
> >>
> >> On Aug 25, 2018, at 1:58 PM, Rene Gassmoeller <
> >> rene.gassmoeller at mailbox.org> wrote:
> >>
> >> Hi Max,
> >>
> >> let me cc this to the mailing list, as it might be of interest to
> others.
> >>
> >> I found this old email thread on the mailing list that seems to
> >> discuss the same issue:
> >> http://lists.geodynamics.org/pipermail/aspect-devel/2012-November/0001
> >> 38.html (there are some additional messages when you scroll through
> >> the thread). In
> >> short: At the time Ian found that this is a problem when the boundary
> >> layer is underresolved. Since the element at the surface are larger
> >> than the elements at the CMB they do resolve a different temperature
> >> gradient for the boundary layer if the boundary layer is not
> >> sufficiently resolved, and we compute the heat flux based on the
> temperature gradient.
> >>
> >> However, thinking about it again, I do not think this can be the only
> >> reason. If the difference were caused by numerical resolution alone,
> >> the average temperature should change, until the two heat fluxes are
> >> in balance as you mentioned. But another idea I had was that the
> >> computed heat flux is a postprocessor that only considers physical
> >> diffusion (i.e. diffusion due to thermal conductivity), while in the
> >> assembly we also include the numerical diffusion to stabilize the
> >> equation. So maybe the difference between the heat fluxes is caused by
> >> numerical diffusion (and it is so large because the cells are
> >> relatively coarse)? If that were the case the difference should
> >> increase/decrease with higher/lower resolution, as the stabilization
> >> scales with the element size. Unfortunately that is not a criterion to
> >> distinguish between the two possible causes (obviously resolution of
> >> the temperature gradient across the boundary layer also scales with
> >> element size). You could try to re-run the model with discontinuous
> >> temperature elements, that should disable the numerical diffusion. Let
> me know if that makes sense.
> >>
> >> Best,
> >>
> >> Rene
> >>
> >> On 08/24/2018 06:46 PM, Max Rudolph wrote:
> >>
> >> Rene and Wolfgang,
> >> We're running some models in 3D spherical geometry in aspect with
> >> internal heating and are unable to reproduce thermal histories
> >> consistent with what we saw in identically configured calculations in
> >> CitcomS. It is possible that there is some aspect of our model setup
> >> that is inconsistent, but we are using the same depth profile of
> >> viscosity, same boundary conditions, and same initial conditions. The
> >> ASPECT models cool much faster (which also drives up viscosity). We
> >> are trying to get understand these results, and I was looking at the
> >> heat transport in the 3D spherical model discussed in section 5.3.2 of
> >> the aspect manual 'simple convection in the 3D spherical shell'. In
> >> figure 40, it appears that this calculation has been run to nearly
> >> statistically steady state. Is this correct? I am wondering,
> >> basically, how to interpret the imbalance between reported surface and
> >> CMB heat flow. To follow up on this, I ran the cookbook .prm file from
> >> 5.3.2
> >> (shell_simple_3d) for 5 Gyr, which appears to be very close to or at
> >> statistically steady state at the end of the calculation. Due to
> >> resource constraints, I only ran at global refinement level 3
> >> (definitely under-resolved, mea culpa). However, at the end of the
> >> calculation, there is still a significant imbalance between the
> >> surface and core mantle boundary heat flows. This is a model with no
> >> internal heating, so at statistically steady state, the surface and
> >> cmb heat flows should balance to conserve energy. Do you have any
> >> insight into why we might be seeing this behavior? I've attached the
> >> input .prm file, the statistics file, and a plot showing the average
> temperature and surface and cmb heat flows.
> >> Thank you for thinking about this!
> >>
> >> Best
> >> Max
> >>
> >> p.s. I sent this to both of you because it looked to me like Wolfgang
> >> contributed the cookbook. Please let me know if there is someone else
> >> that might be better positioned to help with this.
> >>
> >> <image.png>
> >>
> >>
> >>
> >>
> >> --
> >> Rene Gassmoellerhttps://gassmoeller.github.io/
> >>
> >> _______________________________________________
> >> Aspect-devel mailing list
> >> Aspect-devel at geodynamics.org
> >> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
> >>
> >>
> >> _______________________________________________
> >> Aspect-devel mailing list
> >> Aspect-devel at geodynamics.org
> >> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
> > -------------- next part --------------
> > An HTML attachment was scrubbed...
> > URL: <
> http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20180825/41f8abc5/attachment.html
> >
> >
> > ------------------------------
> >
> > Subject: Digest Footer
> >
> > _______________________________________________
> > Aspect-devel mailing list
> > Aspect-devel at geodynamics.org
> > http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
> >
> > ------------------------------
> >
> > End of Aspect-devel Digest, Vol 81, Issue 7
> > *******************************************
> > <test stabilization.png>
> > _______________________________________________
> > Aspect-devel mailing list
> > Aspect-devel at geodynamics.org
> > http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20180828/14f75630/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 24255 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20180828/14f75630/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 28810 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20180828/14f75630/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 28271 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20180828/14f75630/attachment-0005.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cooling_spherical_shell.prm
Type: application/octet-stream
Size: 2584 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20180828/14f75630/attachment-0001.obj>