Problems with shell_simple_2d/3d cookbook examples

[rclam]

Version: ASPECT ver. 2.3.0, deal.II 9.2.0, Trilinos 12.18.1, p4est 2.2.0, in OPTIMIZED mode, with 1 MPI process

We have noticed two separate issue with ASPECT: one of the built-in cookbooks (shell_simple_3d.prm) is always terminating near immediately; there seems to be some numerical instability in the solver for the 2D spherical shell.

Firstly, running an unedited version of the shell_simple_3d.prm does not yield any output vtu files. The terminal output is killed almost immediately at time step 0. It claims that there are zero iterations for solving temperature system for each instance. The Stokes preconditioner is rebuilt twice (with 200+15 iterations and 200+16 iterations respectively) before the run is “killed”. My computer has 64 GB memory and no other processes were running, so it would be surprising if the issue is a question of RAM. However if it is, it would be nice if ASPECT's error message clearly stated so.

Onto the second issue: Using an edited version of the shell_simple_2d.prm, we attempted to model a system where the initial temperature field is a with a linear function of depth (as opposed to the spherical hexagonal perturbation model in the original). The set-up is shown below.

subsection Initial temperature model
  set Model name = function

  subsection Function
    set Variable names = x,y
    set Function constants = k=1e-6, tau=3e16, r_in = 3481000, T_1=774, T_0=3974, r_outer=6336000

    set Function expression =  T_0 + (   ((T_0-T_1)/(r_outer - r_in)) * ( -(sqrt(x*x + y*y)) + r_in)) #- 273

  end
end

Looking at the Temperature output files below, we can see clearly (in time steps 3 and 4) multiple rings of sudden temperature fluctuations holding along the shell radially.

We decided to look at other variables and found fluctuations also appear in the velocity output (magnitude below, x- and y- later) files.

Here are the x- and y-components at time 0 for the linear case:

Upon discovering this strange behavior, we decided to try the original shell_simple_2d.prm file and see if this also occurred there. Looking at the temperature output closely, there seems to be similar fluctuations near the outer radius (see image below).

[tjhei]

Thank you for reporting this issue. This is indeed not how it should look like. We might have an underresolved mesh for the viscosity chosen here. You can try that first.

I am a bit surprised how this could have changed at some point in the past because the prm itself did not change. Could this be from

commit fedeb6a29559291016f80a3a982da5b7480e9f20
Author: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date:   Sat Jan 20 13:33:28 2018 +0100

    Remove radial earth-like

diff --git a/cookbooks/shell_simple_2d.prm b/cookbooks/shell_simple_2d.prm
index d9f4d1392..fe31cd0e3 100644
--- a/cookbooks/shell_simple_2d.prm
+++ b/cookbooks/shell_simple_2d.prm
@@ -57,7 +57,7 @@ end

 subsection Gravity model
-  set Model name = radial earth-like
+  set Model name = ascii data
 end

@gassmoeller , what do you think?

[jdannberg]

I would be surprised if the change in the gravity profile would have such a big effect. Both profiles have values of approximately 10 m/s^2 (the new one -- the ascii data -- is from PREM).

This looks like an issue of resolution to me. We designed the simple 2d shell model to be a model that can be run quickly on a laptop, so it has a limited resolution. If you now pick an initial condition that does not have any perturbations (using a linear temperature profile, rather than the spherical hexagonal perturbation), the velocities are approximately zero at the start of the model run (basically, there are no lateral density variations that would drive the flow). This means that instead of physical properties, it's now the discretization that decides where material starts to move. The structure of the mesh always influences the solution a little, and if the velocities in the model are only around 1e-4, this influence becomes really visible, and the influence of the mesh becomes really strong. This will improve with (a) choosing a higher resolution, and (b) starting with some sort of initial perturbation. You also see that the velocity looks fine in the later time steps, where the model has evolved to actually have lateral density variations.

The oscillations in the temperature also look like a problem of resolution to me. We apply a stabilization to the temperature equation (some amount of numerical diffusion), but if the resolution is really coarse, this is simply not enough and there are some over- and undershoots. In the original model, where the flow structures are on a larger length scale, the problem is much less pronounced.

One thing that has changed since we made the cookbook is that we now by default write the visualization output not only at the vertices of the cells, but we also evaluate the solution at additional points in between (since we use second-order finite element, this better reflects the actual number of points where we compute the solution on). However, that has the disadvantage that sharp temperature gradients within a cell can cause that interpolation of the solution within the cell to have over- and undershoots, which might be what you see in the output. So again, this should be fixed by a higher resolution, or you can see if it goes away if you set

set Interpolate output = false

in the Visualization subsection (which will then only output the solution at the vertices of the mesh).

The 3d cookbook not working is a separate issue we will need to look into. I haven't tried to run that, but I just looked through the section of the manual, and it mentions significant computational cost and that "Computations of this size typically run with 1000 MPI processes". I would only try to run that on a cluster. But you are right, we should be more clear about that. @bangerth I think you ran that model? Do you still remember on how many cores/for what time you ran that model? Could you add that as a note to the start of the cookbook? I know it was a long time ago, but just having approximate values would already be useful.

[elodie-kendall]

The interpolation and small change to the gravity profile didn't fix the shell_simple_2D.prm but adding an extra global refinement means this cookbook now works. Thanks @jdannberg ! :) Let me know if you want me to fix as a PR

[rclam]

Thanks for the advice! The cookbook now works well for me.

[bangerth]

For the record, I can reproduce the crash. It would be useful to bisect which change broke this cookbook.

[bangerth]

As for the question of how many cores I used -- I have no recollection at all.