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

Max Rudolph maxrudolph at ucdavis.edu
Wed Aug 29 16:11:37 PDT 2018 

Gerry,
Could you please send me the talk slides?

We have had some problems with over/undershoots using DG, even when using
the bound preserving limiter in aspect. These are not small overshoots
either. A recent case had min/max values of -0.36 and 1.37 even though the
bound preserving limiter was set to [0,1].

Max

On Wed, Aug 29, 2018 at 2:17 PM Elbridge Gerry Puckett <egp at math.ucdavis.edu>
wrote:

> Hi Everyone,
>
> I'm sorry to come to this discussion so late.  I have been meaning to read
> the entire sequence of posts before posting myself.
>
> Ying He implemented DG with a Bound Preserving (BP) Limiter for the
> temperature equation (with incomprehensible flow).  We (Ying, Magali,
> Louise, and myself) submitted a poster to the 2016 AGU fall meeting and it
> was chosen for a talk, which Ying gave.
>
> I have the talk slides and they appear to be relevant to to this
> discussion, even though we did not study internal heating.  Again, I would
> have liked to read the entire discussion before posting, but Max just
> mentioned DG and hence it seems appropriate for me to offer Ying's talk
> slides to the [aspect-devel] mailing list as well as the rest of the
> community.
> A few words about what Ying / we did.  Ying was interested in studying the
> effectiveness of this new methodology as a function of the  Péclet number
> (Pe) number.  For many of the problems we were computing as a group, such
> as our recent LLSVP paper in the special issue of PEPI to associated with
> the 15th SEDI conference, we (Ying) determined that at that  Péclet number
> the changes to our results were negligible when one used the original
> algorithm in ASPECT with EV and the new DGBP algorithm.
>
> However, for other problems, Magali had an argument that for some problems
> one might want to compute in ASPECT there would be an *effective local
> Péclet number, *for which there are decided differences between the two
> algorithms.  This difference is shown in the talk slides.  (There is a
> movie of this, which I can probably find.)
>
> I don't see this material in the manual, which is unfortunate. I can put
> this on my to-do list if others think it would be helpful.
>
> Two things to note:
>
>
>    1. Part of Ying's to-do list was to add the ability to use AMR to
>    algorithm with this DGBP temperature advection-diffusion algorithm, which I
>    think remains to be done.
>
>
>    1. Another issue to note is that that Ying found it necessary to use a
>    second-order accurate Strong Stability Preserving time stepping algorithm
>    (SSP2) instead of the usual BDF2 for the DGBP algorithm. This is from
>    Ying's 2017 PEPI paper with Magali and me.
>
> "For the case when the FEM is used for the spatial discretization, we use
> the same semi-implicit BDF2 scheme as shown in Kronbichler et al. (2012)
> except we fix the time-step size dt for simplicity. However, if the DG
> method is used, then a Strong Stability-Preserving (SSP) high order time
> discretization method (Gottlieb, 2005; Gottlieb et al., 2001) is applied in
> order to maintain the bound preserving property from the limiter."
>
> This does not seem to be mentioned in the talk slides, probably due to
> lack of time.  A search on SSP2 in the ASPECT manual generated June 22,
> 2018 yields zero hits and and even a search on DGBP yields only three hits
> in thre parameter section.  There is probably room for more documentation
> there. : -(
>
> Anyway, the point is that one might keep (1) and (2) in mind while using
> DGBP for the temperature equation.
>
> In summary, I didn't intend to write such a long post, I hope some of you
> find it is worthwhile.
>
> Shall I post the talk slides and if so, where?  I can put them on
> math.ucdavis.edu/~egp/... or post to this discussion with the slides
> appended, or ..
>
> Sincerely,
>
> - Gerry
>
> On 08/29/2018 12:49 PM, Max Rudolph wrote:
>
> Rene,
> Thank you for summarizing the current state of affairs. I would suggest
> that when you make the changes in (1)-(2), you also make ASPECT print at
> least a warning when the entropy viscosity exceeds, say, 1% of the thermal
> conductivity. It might be even better for the code to exit in with an error
> unless a parameter is explicitly set to ignore such a condition.
>
> It also seems that one way around (1) is to use DG elements for
> composition. Is this what most users are doing anyways?
>
> Thank you for addressing these issues so quickly. I am happy to help with
> testing as needed.
>
> Best,
> Max
>
> On Wed, Aug 29, 2018 at 12:16 PM Rene Gassmoeller <
> rene.gassmoeller at mailbox.org> wrote:
>
>> Hi all,
>>
>> I think by now we have a pretty good understanding of the problem,
>> however there is no clear path forward yet. I have done some further
>> testing with the shell_simple_3d cookbook:
>>
>>  - As Max pointed out the artificial viscosity at resolution <= 2 global
>> refinement in spherical models is bigger than the natural diffusion,
>> although it reduces drastically with resolution (compare to a natural
>> conductivity of ~4 W/m*K).
>>
>> Global refinement / Maximum artificial viscosity timestep 0 / timestep 2:
>>
>> 1 (4 radial elements): 94.1 W/m*K  /  59.9 W/m*K
>>
>> 2 (8 radial elements): 48.9 W/m*K / 6.12 W/m*K
>>
>> 3 (16 radial elements): 24.8 W/m*K / 0.92 W/m*K
>>
>> There are some pitfalls here, in particular the artificial viscosity in
>> the output of timestep 0 is in general much bigger than in later timesteps
>> (and it only scales linearly with cell size), because we do not have a
>> velocity solution yet. Timestep 1 is somewhat unreliable as well, because
>> the BDF2 scheme is not fully initialized. Nevertheless, for later timesteps
>> the artificial viscosity scales at least with cell size squared (h^2) as
>> expected (actually slightly better, because it is not only h^2, but also
>> based on the residual of the equation, which reduces as well, ideally with
>> h^3). This means higher resolutions should significantly reduce the total
>> amount of artificial diffusion, although it might still be significant (>1%
>> of total diffusion) up to resolution 4 or so, which seems unacceptably
>> expensive. Moreover, as Scott and Cedric mentioned there will always be an
>> imbalance between the heat fluxes at the top and bottom boundary with this
>> method, unless we use radially uniformly spaced elements, or disable the
>> stabilization, although the magnitude of the imbalance should reduce with
>> increased resolution. Also the thermal conductivity within the domain will
>> be unequally distributed (artificial conductivity in the upper mantle will
>> be higher than in the lower mantle).
>>
>> While the above point suggests that the artificial viscosity scheme is
>> correctly implemented, there still exist three points for possible
>> improvements:
>>
>>  - As Wolfgang mentioned the parameters that were originally chosen for
>> the artificial viscosity scheme were different from what they are today,
>> because we noticed that the original smoothing was too weak to stabilize
>> oscillations in compositional fields. It is completely possible that the
>> original values were sufficient for the temperature (which already has
>> natural diffusion), and in hindsight I should have thought of just using
>> different parameters for composition and temperature (though that was 6
>> years ago, and I was just a 2nd year PhD student at the time I worked on
>> modifying the parameters for stabilizing the composition equation). I can
>> easily make a change that allows for different parameters for temperature
>> and composition, which would allow everyone to test their favorite values,
>> without risking oscillations in compositional fields.
>>
>> -  Even if we can reduce the parameters it is of course possible that as
>> Max pointed out the anisotropic nature of the diffusion in SUPG is more
>> appropriate / "less wrong" than the isotropic entropy viscosity for our
>> problem. As Juliane and Timo pointed out we could reimplement the content
>> of https://github.com/geodynamics/aspect/pull/412 and see if SUPG
>> performs better for low resolutions. Does anyone know if the "artificial
>> viscosity" of SUPG also scales with h^2? Because if it is linear, we might
>> implement something that is better for low resolutions, but worse for high
>> resolutions.
>>
>> - As Max mentioned: We should not need a stabilization for pure
>> conduction problems (where velocity is 0), and should modify the algorithm
>> accordingly.
>>
>> So the next steps could be:
>>
>> 1. Allow for different stabilization parameters for temperature and
>> composition, and check which values are still stable.
>> 2. Do not stabilize advection/diffusion solutions where the velocity is
>> zero (because it is only a diffusion equation).
>> 3. Reimplement the SUPG based on
>> https://github.com/geodynamics/aspect/pull/412 and see how it performs
>> (at low and high resolutions).
>>
>> Does that summarize our discussion appropriately? I can easily make the
>> code adjustments 1 and 2 (they are easy and useful in any case), and could
>> also look into creating an initial version of 3 (although it would take a
>> bit of time), but I currently do not have the time for much testing of the
>> methods, so I would be greatful if someone else could do the testing and
>> benchmarking of the methods.
>>
>> Best,
>> Rene
>>
>>
>> On 08/29/2018 09:19 AM, Max Rudolph wrote:
>>
>> On Tue, Aug 28, 2018 at 8:01 PM Wolfgang Bangerth <bangerth at colostate.edu>
>> wrote:
>>
>>> On 08/28/2018 05:33 PM, Max Rudolph wrote:
>>> >  From this, it is very obvious why the solution to the convection
>>> problem at
>>> > low resolution is very diffusive and also why the interior temperature
>>> is much
>>> > closer to the surface temperature than to the CMB temperature because
>>> the
>>> > artificial viscosity is on the order of 20 times larger than the
>>> thermal
>>> > conductivity near the surface.
>>>
>>> Would it be easy to verify whether the artificial viscosity ("artificial
>>> conductivity") decreases at the expected rate with mesh refinement?
>>>
>>
>> What is the most helpful way for me to show this? Visualization of a
>> couple of slices from the 3D conduction model? I tried to get the depth
>> average of artificial viscosity but the postprocessor is not implemented.
>>
>>
>>> > For the conduction problem, the default values of the artificial
>>> viscosity are
>>> > also much larger than the thermal conductivity.
>>>
>>> I think that's the point worth investigating. Since in this case the
>>> velocity
>>> is zero, one would expect the artificial viscosity to also be at least
>>> quite
>>> small. Why is it not?
>>>
>>
>> *Maybe the spherical and 2D annulus geometry models are returning an
>> unhelpful length scale, like planetary radius instead of layer depth?*
>>
>> aspect/source/simulator/entropy_viscosity.cc (starting line 191):
>> // If the velocity is 0 we have to assume a sensible velocity to calculate
>> // an artificial diffusion. We choose similar to nondimensional
>> // formulations: v ~ thermal_diffusivity / length_scale, which cancels
>>     // the density and specific heat from the entropy formulation. It
>> seems
>>     // surprising at first that only the conductivity remains, but
>> remember
>>     // that this actually *is* an additional artificial diffusion.
>>     if (std::abs(global_u_infty) < 1e-50)
>>       return parameters.stabilization_beta *
>>              max_conductivity / geometry_model->length_scale() *
>>              cell_diameter;
>>
>>
>>
>>>
>>> Best
>>>   W.
>>>
>>> --
>>> ------------------------------------------------------------------------
>>> Wolfgang Bangerth          email:                 bangerth at colostate.edu
>>>                             www:
>>> http://www.math.colostate.edu/~bangerth/
>>>
>>>
>>
>> _______________________________________________
>> Aspect-devel mailing listAspect-devel at geodynamics.orghttp://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>>
>>
>> --
>> 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 listAspect-devel at geodynamics.orghttp://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/20180829/d9567d20/attachment-0001.html>

More information about the Aspect-devel mailing list