Simple setups

5.2 Simple setups

5.2.1 Convection in a 2d box

In this first example, let us consider a simple situation: a 2d box of dimensions $[0, 1] \times [0, 1]$ that is heated from below, insulated at the left and right, and cooled from the top. We will also consider the simplest model, the incompressible Boussinesq approximation with constant coefficients $η, ρ_{0}, g, C_{p}, k$ , for this testcase. Furthermore, we assume that the medium expands linearly with temperature. This leads to the following set of equations:

\begin{align} - \nabla \cdot [2 η ε (u)] + \nabla p & = ρ_{0} (1 - α (T - T_{0})) g & in Ω, & (45) \\ \nabla \cdot u & = 0 & in Ω, & (46) \\ ρ_{0} C_{p} (\frac{\partial T}{\partial t} + u \cdot \nabla T) - \nabla \cdot k \nabla T & = 0 & in Ω . & (47) \end{align}

It is well known that we can non-dimensionalize this set of equations by introducing the Rayleigh number $R a = \frac{ρ_{0} g α Δ T h^{3}}{η κ}$ , where $h$ is the height of the box, $κ = \frac{k}{ρ C_{p}}$ is the thermal diffusivity and $Δ T$ is the temperature difference between top and bottom of the box. Formally, we can obtain the non-dimensionalized equations by using the above form and setting coefficients in the following way:

\begin{array}{l} ρ_{0} = C_{p} = κ = α = η = h = Δ T = 1, T_{0} = 0, g = R a, \end{array}

where $g = - g e_{z}$ is the gravity vector in negative $z$ -direction. We will see all of these values again in the input file discussed below. One point to note is that for the Boussinesq approximation, as described above, the density in the temperature equation is chosen as the reference density $ρ_{0}$ rather than the full density $ρ (1 - α (T - T_{0}))$ as we see it in the buoyancy term on the right hand side of the momentum equation. As ASPECT is able to handle different approximations of the equations (see Section 2.10), we also have to specify in the input file that we want to use the Boussinesq approximation. The problem is completed by stating the velocity boundary conditions: tangential flow along all four of the boundaries of the box.

This situation describes a well-known benchmark problem for which a lot is known and against which we can compare our results. For example, the following is well understood:

For values of the Rayleigh number less than a critical number $R a_{c} \approx 780$ , thermal diffusion dominates convective heat transport and any movement in the fluid is damped exponentially. If the Rayleigh number is moderately larger than this threshold then a stable convection pattern forms that transports heat from the bottom to the top boundaries. The simulations we will set up operates in this regime. Specifically, we will choose $R a = 1 0^{4}$ .
On the other hand, if the Rayleigh number becomes even larger, a series of period doublings starts that makes the system become more and more unstable. We will investigate some of this behavior at the end of this section.
For certain values of the Rayleigh number, very accurate values for the heat flux through the bottom and top boundaries are available in the literate. For example, Blankenbach et al. report a non-dimensional heat flux of $4.884409 \pm 0.00001$ , see [BBC $^{+}$ 89]. We will compare our results against this value below.

With this said, let us consider how to represent this situation in practice.

The input file. The verbal description of this problem can be translated into an ASPECT input file in the following way (see Section A for a description of all of the parameters that appear in the following input file, and the indices at the end of this manual if you want to find a particular parameter; you can find the input file to run this cookbook example in cookbooks/convection-_box.prm):

# At the top, we define the number of space dimensions we would like to
# work in:
set ??                              = 2

# There are several global variables that have to do with what
# time system we want to work in and what the end time is. We
# also designate an output directory.
set ?? = false
set ??                               = 0.5
set ??                       = output-convection-box

# Then there are variables that describe how the pressure should
# be normalized. Here, we choose a zero average pressure
# at the surface of the domain (for the current geometry, the
# surface is defined as the top boundary).
set ??                 = surface
set ??                       = 0

# Then come a number of sections that deal with the setup
# of the problem to solve. The first one deals with the
# geometry of the domain within which we want to solve.
# The sections that follow all have the same basic setup
# where we select the name of a particular model (here,
# the box geometry) and then, in a further subsection,
# set the parameters that are specific to this particular
# model.
subsection ??
  set ?? = box

  subsection ??
    set ?? = 1
    set ?? = 1
  end
end

# The next section deals with the initial conditions for the
# temperature (there are no initial conditions for the
# velocity variable since the velocity is assumed to always
# be in a static equilibrium with the temperature field).
# There are a number of models with the ’function’ model
# a generic one that allows us to enter the actual initial
# conditions in the form of a formula that can contain
# constants. We choose a linear temperature profile that
# matches the boundary conditions defined below plus
# a small perturbation:
subsection ??
  set ?? = function

  subsection ??
    set ??      = x,z
    set ??  = p=0.01, L=1, pi=3.1415926536, k=1
    set ?? = (1.0-z) - p*cos(k*pi*x/L)*sin(pi*z)
  end
end

# Then follows a section that describes the boundary conditions
# for the temperature. The model we choose is called ’box’ and
# allows to set a constant temperature on each of the four sides
# of the box geometry. In our case, we choose something that is
# heated from below and cooled from above.
# All other parts of the boundary are insulated (i.e., no heat flux through
# these boundaries; this is also often used to specify symmetry boundaries).
subsection ??
  set ??   = bottom, top
  set ?? = box

  subsection ??
    set ?? = 1
    set ??   = 0
    set ??  = 0
    set ??    = 0
  end
end

# We then also have to prescribe several other parts of the model such as
# which boundaries actually carry a prescribed boundary temperature, whereas
# The next parameters then describe on which parts of the
# boundary we prescribe a zero or nonzero velocity and
# on which parts the flow is allowed to be tangential.
# Here, all four sides of the box allow tangential
# unrestricted flow but with a zero normal component:
subsection ??
  set ?? = left, right, bottom, top
end

# The following two sections describe first the
# direction (vertical) and magnitude of gravity and the
# material model (i.e., density, viscosity, etc). We have
# discussed the settings used here in the introduction to
# this cookbook in the manual already.
subsection ??
  set ?? = vertical

  subsection ??
    set ?? = 1e4   # = Ra
  end
end

subsection ??
  set ?? = simple

  subsection ??
    set ??             = 1
    set ??       = 1
    set ??         = 0
    set ??          = 1
    set ?? = 1
    set ??                     = 1
  end
end

# We also have to specify that we want to use the Boussinesq
# approximation (assuming the density in the temperature
# equation to be constant, and incompressibility).
subsection ??
  set ?? = Boussinesq approximation
end

# The settings above all pertain to the description of the
# continuous partial differential equations we want to solve.
# The following section deals with the discretization of
# this problem, namely the kind of mesh we want to compute
# on. We here use a globally refined mesh without
# adaptive mesh refinement.
subsection ??
  set ??                = 4
  set ??              = 0
  set ??       = 0
end

# The final part is to specify what ASPECT should do with the
# solution once computed at the end of every time step. The
# process of evaluating the solution is called ‘postprocessing’
# and we choose to compute velocity and temperature statistics,
# statistics about the heat flux through the boundaries of the
# domain, and to generate graphical output files for later
# visualization. These output files are created every time
# a time step crosses time points separated by 0.01. Given
# our start time (zero) and final time (0.5) this means that
# we will obtain 50 output files.
subsection ??
  set ?? = velocity statistics, temperature statistics, ...
                                           ... heat flux statistics, visualization

  subsection ??
    set ?? = 0.01
  end
end

subsection ??
  set ?? = 1e-10
end

Running the program. When you run this program for the first time, you are probably still running ASPECT in debug mode (see Section 4.3) and you will get output like the following:

Number of active cells: 256 (on 5 levels)
Number of degrees of freedom: 3,556 (2,178+289+1,089)

*** Timestep 0:  t=0 seconds
   Solving temperature system... 0 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 31+0 iterations.

[... ...]

*** Timestep 1085:  t=0.5 seconds
   Solving temperature system... 0 iterations.
   Solving Stokes system... 5 iterations.

   Postprocessing:
     RMS, max velocity:                  43.5 m/s, 70.3 m/s
     Temperature min/avg/max:            0 K, 0.5 K, 1 K
     Heat fluxes through boundary parts: 0.01977 W, -0.01977 W, -4.787 W, 4.787 W

Termination requested by criterion: end time

+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |      66.5s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble Stokes system          |      1086 |      8.63s |        13% |
| Assemble temperature system     |      1086 |        32s |        48% |
| Build Stokes preconditioner     |         1 |    0.0225s |         0% |
| Build temperature preconditioner|      1086 |      1.52s |       2.3% |
| Solve Stokes system             |      1086 |       7.7s |        12% |
| Solve temperature system        |      1086 |     0.729s |       1.1% |
| Initialization                  |         1 |    0.0316s |         0% |
| Postprocessing                  |      1086 |      7.76s |        12% |
| Setup dof systems               |         1 |    0.0104s |         0% |
| Setup initial conditions        |         1 |   0.00621s |         0% |
+---------------------------------+-----------+------------+------------+

If you’ve read up on the difference between debug and optimized mode (and you should before you switch!) then consider disabling debug mode. If you run the program again, every number should look exactly the same (and it does, in fact, as I am writing this) except for the timing information printed every hundred time steps and at the end of the program:

+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |      25.8s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble Stokes system          |      1086 |      2.51s |       9.7% |
| Assemble temperature system     |      1086 |      9.88s |        38% |
| Build Stokes preconditioner     |         1 |    0.0271s |      0.11% |
| Build temperature preconditioner|      1086 |      1.58s |       6.1% |
| Solve Stokes system             |      1086 |      6.38s |        25% |
| Solve temperature system        |      1086 |     0.542s |       2.1% |
| Initialization                  |         1 |     0.219s |      0.85% |
| Postprocessing                  |      1086 |      2.79s |        11% |
| Setup dof systems               |         1 |      0.23s |      0.89% |
| Setup initial conditions        |         1 |     0.107s |      0.41% |
+---------------------------------+-----------+------------+------------+

In other words, the program ran more than 2 times faster than before. Not all operations became faster to the same degree: assembly, for example, is an area that traverses a lot of code both in ASPECT and in deal.II and so encounters a lot of verification code in debug mode. On the other hand, solving linear systems primarily requires lots of matrix vector operations. Overall, the fact that in this example, assembling linear systems and preconditioners takes so much time compared to actually solving them is primarily a reflection of how simple the problem is that we solve in this example. This can also be seen in the fact that the number of iterations necessary to solve the Stokes and temperature equations is so low. For more complex problems with non-constant coefficients such as the viscosity, as well as in 3d, we have to spend much more work solving linear systems whereas the effort to assemble linear systems remains the same.

Visualizing results. Having run the program, we now want to visualize the numerical results we got. ASPECT can generate graphical output in formats understood by pretty much any visualization program (see the parameters described in Section ??) but we will here follow the discussion in Section 4.4 and use the default VTU output format to visualize using the Visit program.

In the parameter file we have specified that graphical output should be generated every 0.01 time units. Looking through these output files (which can be found in the folder output-convection-box, as specified in the input file), we find that the flow and temperature fields quickly converge to a stationary state. Fig. 7 shows the initial and final states of this simulation.

Figure 7: Convection in a box: Initial temperature and velocity field (left) and final state (right).

There are many other things we can learn from the output files generated by ASPECT, specifically from the statistics file that contains information collected at every time step and that has been discussed in Section 4.4.2. In particular, in our input file, we have selected that we would like to compute velocity, temperature, and heat flux statistics. These statistics, among others, are listed in the statistics file whose head looks like this for the current input file:

# 1: Time step number
# 2: Time (seconds)
# 3: Time step size (seconds)
# 4: Number of mesh cells
# 5: Number of Stokes degrees of freedom
# 6: Number of temperature degrees of freedom
# 7: Iterations for temperature solver
# 8: Iterations for Stokes solver
# 9: Velocity iterations in Stokes preconditioner
# 10: Schur complement iterations in Stokes preconditioner
# 11: RMS velocity (m/s)
# 12: Max. velocity (m/s)
# 13: Minimal temperature (K)
# 14: Average temperature (K)
# 15: Maximal temperature (K)
# 16: Average nondimensional temperature (K)
# 17: Outward heat flux through boundary with indicator 0 ("left") (W)
# 18: Outward heat flux through boundary with indicator 1 ("right") (W)
# 19: Outward heat flux through boundary with indicator 2 ("bottom") (W)
# 20: Outward heat flux through boundary with indicator 3 ("top") (W)
# 21: Visualization file name
... lots of numbers arranged in columns ...

Fig. 8 shows the results of visualizing the data that can be found in columns 2 (the time) plotted against columns 11 and 12 (root mean square and maximal velocities). Plots of this kind can be generated with Gnuplot by typing (see Section 4.4.2 for a more thorough discussion):

plot "output-convection-box/statistics" using 2:11 with lines

Fig. 8 shows clearly that the simulation enters a steady state after about $t \approx 0.1$ and then changes very little. This can also be observed using the graphical output files from which we have generated Fig. 7. One can look further into this data to find that the flux through the top and bottom boundaries is not exactly the same (up to the obvious difference in sign, given that at the bottom boundary heat flows into the domain and at the top boundary out of it) at the beginning of the simulation until the fluid has attained its equilibrium. However, after $t \approx 0.2$ , the fluxes differ by only $5 \cdot 1 0^{- 5}$ , i.e., by less than 0.001% of their magnitude.¹⁹ The flux we get at the last time step, 4.787, is less than 2% away from the value reported in [BBC $^{+}$ 89] ( $\approx$ 4.88) although we compute on a $16 \times 16$ mesh and the values reported by Blankenbach are extrapolated from meshes of size up to $72 \times 72$ . This shows the accuracy that can be obtained using a higher order finite element. Secondly, the fluxes through the left and right boundary are not exactly zero but small. Of course, we have prescribed boundary conditions of the form $\frac{\partial T}{\partial n} = 0$ along these boundaries, but this is subject to discretization errors. It is easy to verify that the heat flux through these two boundaries disappears as we refine the mesh further.

Figure 8: Convection in a box: Root mean square and maximal velocity as a function of simulation time (left). Heat flux through the four boundaries of the box (right).

Furthermore, ASPECT automatically also collects statistics about many of its internal workings. Fig. 9 shows the number of iterations required to solve the Stokes and temperature linear systems in each time step. It is easy to see that these are more difficult to solve in the beginning when the solution still changes significant from time step to time step. However, after some time, the solution remains mostly the same and solvers then only need 9 or 10 iterations for the temperature equation and 4 or 5 iterations for the Stokes equations because the starting guess for the linear solver – the previous time step’s solution – is already pretty good. If you look at any of the more complex cookbooks, you will find that one needs many more iterations to solve these equations.

Figure 9: Convection in a box: Number of linear iterations required to solve the Stokes and temperature equations in each time step.

Play time 1: Different Rayleigh numbers. After showing you results for the input file as it can be found in cookbooks/convection-_box.prm, let us end this section with a few ideas on how to play with it and what to explore. The first direction one could take this example is certainly to consider different Rayleigh numbers. As mentioned above, for the value $R a = 1 0^{4}$ for which the results above have been produced, one gets a stable convection pattern. On the other hand, for values $R a < R a_{c} \approx 780$ , any movement of the fluid dies down exponentially and we end up with a situation where the fluid doesn’t move and heat is transported from the bottom to the top only through heat conduction. This can be explained by considering that the Rayleigh number in a box of unit extent is defined as $R a = \frac{ρ_{0} g α Δ T}{η k}$ . A small Rayleigh number means that the viscosity is too large (i.e., the buoyancy given by the product of the magnitude of gravity times the density anomalies caused by temperature – $ρ_{0} α Δ T$ – is not strong enough to overcome friction forces within the fluid).

On the other hand, if the Rayleigh number is large (i.e., the viscosity is small or the buoyancy large) then the fluid develops an unsteady convection period. As we consider fluids with larger and larger $R a$ , this pattern goes through a sequence of period-doubling events until flow finally becomes chaotic. The structures of the flow pattern also become smaller and smaller.

Figure 10: Convection in a box: Temperature fields at the end of a simulation for

R a = 1 0^{2}

where thermal diffusion dominates (left) and

R a = 1 0^{6}

where convective heat transport dominates (right). The mesh on the right is clearly too coarse to resolve the structure of the solution.

Figure 11: Convection in a box: Velocities (left) and heat flux across the top and bottom boundaries (right) as a function of time at

R a = 1 0^{6}

We illustrate these situations in Fig.s 10 and 11. The first shows the temperature field at the end of a simulation for $R a = 1 0^{2}$ (below $R a_{c}$ ) and at $R a = 1 0^{6}$ . Obviously, for the right picture, the mesh is not fine enough to accurately resolve the features of the flow field and we would have to refine it more. The second of the figures shows the velocity and heatflux statistics for the computation with $R a = 1 0^{6}$ ; it is obvious here that the flow no longer settles into a steady state but has a periodic behavior. This can also be seen by looking at movies of the solution.

To generate these results, remember that we have chosen $g = R a$ in our input file. In other words, changing the input file to contain the parameter setting

subsection ??
  subsection ??
    set ?? = 1e6   # = Ra
  end
end

will achieve the desired effect of computing with $R a = 1 0^{6}$ .

Play time 2: Thinking about finer meshes. In our computations for $R a = 1 0^{4}$ we used a $16 \times 16$ mesh and obtained a value for the heat flux that differed from the generally accepted value from Blankenbach et al. [BBC $^{+}$ 89] by less than 2%. However, it may be interesting to think about computing even more accurately. This is easily done by using a finer mesh, for example. In the parameter file above, we have chosen the mesh setting as follows:

subsection ??
  set ??                = 4
  set ??              = 0
  set ??       = 0
end

We start out with a box geometry consisting of a single cell that is refined four times. Each time we split each cell into its 4 children, obtaining the $16 \times 16$ mesh already mentioned. The other settings indicate that we do not want to refine the mesh adaptively at all in the first time step, and a setting of zero for the last parameter means that we also never want to adapt the mesh again at a later time. Let us stick with the never-changing, globally refined mesh for now (we will come back to adaptive mesh refinement again at a later time) and only vary the initial global refinement. In particular, we could choose the parameter Initial global refinement to be 5, 6, or even larger. This will get us closer to the exact solution albeit at the expense of a significantly increased computational time.

A better strategy is to realize that for $R a = 1 0^{4}$ , the flow enters a steady state after settling in during the first part of the simulation (see, for example, the graphs in Fig. 8). Since we are not particularly interested in this initial transient process, there is really no reason to spend CPU time using a fine mesh and correspondingly small time steps during this part of the simulation (remember that each refinement results in four times as many cells in 2d and a time step half as long, making reaching a particular time at least 8 times as expensive, assuming that all solvers in ASPECT scale perfectly with the number of cells). Rather, we can use a parameter in the ASPECT input file that let’s us increase the mesh resolution at later times. To this end, let us use the following snippet for the input file:

subsection ??
  set ??                = 3
  set ??              = 0
  set ??       = 0
  set ??              = 0.2, 0.3, 0.4
  set ??                      = 1.0
  set ??                      = 0.0
end

What this does is the following: We start with an $8 \times 8$ mesh (3 times globally refined) but then at times $t = 0.2, 0.3$ and $0.4$ we refine the mesh using the default refinement indicator (which one this is is not important because of the next statement). Each time, we refine, we refine a fraction 1.0 of the cells, i.e., all cells and we coarsen a fraction of 0.0 of the cells, i.e. no cells at all. In effect, at these additional refinement times, we do another global refinement, bringing us to refinement levels 4, 5 and finally 6.

Figure 12: Convection in a box: Refinement in stages. Total number of unknowns in each time step, including all velocity, pressure and temperature unknowns (left) and heat flux across the top boundary (right).

Fig. 12 shows the results. In the left panel, we see how the number of unknowns grows over time (note the logscale for the $y$ -axis). The right panel displays the heat flux. The jumps in the number of cells is clearly visible in this picture as well. This may be surprising at first but remember that the mesh is clearly too coarse in the beginning to really resolve the flow and so we should expect that the solution changes significantly if the mesh is refined. This effect becomes smaller with every additional refinement and is barely visible at the last time this happens, indicating that the mesh before this refinement step may already have been fine enough to resolve the majority of the dynamics.

In any case, we can compare the heat fluxes we obtain at the end of these computations: With a globally four times refined mesh, we get a value of 4.787 (an error of approximately 2% against the accepted value from Blankenbach, $4.884409 \pm 0.00001$ ). With a globally five times refined mesh we get 4.879, and with a globally six times refined mesh we get 4.89 (an error of almost 0.1%). With the mesh generated using the procedure above we also get 4.89 with the digits printed on the screen²⁰ (also corresponding to an error of almost 0.1%). In other words, our simple procedure of refining the mesh during the simulation run yields the same accuracy as using the mesh that is globally refined in the beginning of the simulation, while needing a much lower compute time.

Play time 3: Changing the finite element in use. Another way to increase the accuracy of a finite element computation is to use a higher polynomial degree for the finite element shape functions. By default, ASPECT uses quadratic shape functions for the velocity and the temperature and linear ones for the pressure. However, this can be changed with a single number in the input file.

Before doing so, let us consider some aspects of such a change. First, looking at the pictures of the solution in Fig. 7, one could surmise that the quadratic elements should be able to resolve the velocity field reasonably well given that it is rather smooth. On the other hand, the temperature field has a boundary layer at the top and bottom. One could conjecture that the temperature polynomial degree is therefore the limiting factor and not the polynomial degree for the flow variables. We will test this conjecture below. Secondly, given the nature of the equations, increasing the polynomial degree of the flow variables increases the cost to solve these equations by a factor of $\frac{22}{9}$ in 2d (you can get this factor by counting the number of degrees of freedom uniquely associated with each cell) but leaves the time step size and the cost of solving the temperature system unchanged. On the other hand, increasing the polynomial degree of the temperature variable from 2 to 3 requires $\frac{9}{4}$ times as many degrees of freedom for the temperature and also requires us to reduce the size of the time step by a factor of $\frac{2}{3}$ . Because solving the temperature system is not a dominant factor in each time step (see the timing results shown at the end of the screen output above), the reduction in time step is the only important factor. Overall, increasing the polynomial degree of the temperature variable turns out to be the cheaper of the two options.

Following these considerations, let us add the following section to the parameter file:

subsection ??
set ?? = 2
set ?? = 3
end

This leaves the velocity and pressure shape functions at quadratic and linear polynomial degree but increases the polynomial degree of the temperature from quadratic to cubic. Using the original, four times globally refined mesh, we then get the following output:

Number of active cells: 256 (on 5 levels)
Number of degrees of freedom: 4,868 (2,178+289+2,401)

*** Timestep 0:  t=0 seconds
   Solving temperature system... 0 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 30+0 iterations.

[... ...]

*** Timestep 1621:  t=0.5 seconds
   Solving temperature system... 0 iterations.
   Solving Stokes system... 1+0 iterations.

   Postprocessing:
     RMS, max velocity:                  42.9 m/s, 69.5 m/s
     Temperature min/avg/max:            0 K, 0.5 K, 1 K
     Heat fluxes through boundary parts: -0.004602 W, 0.004602 W, -4.849 W, 4.849 W

Termination requested by criterion: end time

+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |      53.6s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble Stokes system          |      1622 |      4.04s |       7.5% |
| Assemble temperature system     |      1622 |      24.4s |        46% |
| Build Stokes preconditioner     |         1 |    0.0121s |         0% |
| Build temperature preconditioner|      1622 |      8.05s |        15% |
| Solve Stokes system             |      1622 |      8.92s |        17% |
| Solve temperature system        |      1622 |      1.67s |       3.1% |
| Initialization                  |         1 |    0.0327s |         0% |
| Postprocessing                  |      1622 |      4.27s |         8% |
| Setup dof systems               |         1 |   0.00418s |         0% |
| Setup initial conditions        |         1 |   0.00236s |         0% |
+---------------------------------+-----------+------------+------------+

The heat flux through the top and bottom boundaries is now computed as 4.878. Using the five times globally refined mesh, it is 4.8837 (an error of 0.015%). This is 6 times more accurate than the once more globally refined mesh with the original quadratic elements, at a cost significantly smaller. Furthermore, we can of course combine this with the mesh that is gradually refined as simulation time progresses, and we then get a heat flux that is equal to 4.884446, also only 0.01% away from the accepted value!

As a final remark, to test our hypothesis that it was indeed the temperature polynomial degree that was the limiting factor, we can increase the Stokes polynomial degree to 3 while leaving the temperature polynomial degree at 2. A quick computation shows that in that case we get a heat flux of 4.747 – almost the same value as we got initially with the lower order Stokes element. In other words, at least for this testcase, it really was the temperature variable that limits the accuracy.

5.2.2 Convection in a 3d box

The world is not two-dimensional. While the previous section introduced a number of the knobs one can play with with ASPECT, things only really become interesting once one goes to 3d. The setup from the previous section is easily adjusted to this and in the following, let us walk through some of the changes we have to consider when going from 2d to 3d. The full input file that contains these modifications and that was used for the simulations we will show subsequently can be found at cookbooks/convection-_box-_3d.prm.

The first set of changes has to do with the geometry: it is three-dimensional, and we will have to address the fact that a box in 3d has 6 sides, not the 4 we had previously. The documentation of the “box” geometry (see Section ??) states that these sides are numbered as follows: “in 3d, boundary indicators 0 through 5 indicate left, right, front, back, bottom and top boundaries.” Recalling that we want tangential flow all around and want to fix the temperature to known values at the bottom and top, the following will make sense:

set ??                              = 3

subsection ??
  set ?? = box

  subsection ??
    set ?? = 1
    set ?? = 1
    set ?? = 1
  end
end

subsection ??
  set ??   = bottom, top
  set ?? = box

  subsection ??
    set ?? = 1
    set ??    = 0
  end
end

subsection ??
  set ?? = left, right, front, back, bottom, top
end

The next step is to describe the initial conditions. As before, we will use an unstably layered medium but the perturbation now needs to be both in $x$ - and $y$ -direction

subsection ??
  set ?? = function

  subsection ??
    set ??      = x,y,z
    set ??  = p=0.01, L=1, pi=3.1415926536, k=1
    set ?? = (1.0-z) - p*cos(k*pi*x/L)*sin(pi*z)*y^3
  end
end

The third issue we need to address is that we can likely not afford a mesh as fine as in 2d. We choose a mesh that is refined 3 times globally at the beginning, then 3 times adaptively, and is then adapted every 15 time steps. We also allow one additional mesh refinement in the first time step following $t = 0.003$ once the initial instability has given way to a more stable pattern:

subsection ??
  set ??                = 3
  set ??              = 3
  set ??       = 15

  set ??              = 0.003
end

Finally, as we have seen in the previous section, a computation with $R a = 1 0^{4}$ does not lead to a simulation that is exactly exciting. Let us choose $R a = 1 0^{6}$ instead (the mesh chosen above with up to 7 refinement levels after some time is fine enough to resolve this). We can achieve this in the same way as in the previous section by choosing $α = 1 0^{- 10}$ and setting

subsection ??
  set ?? = vertical

  subsection ??
    set ?? = 1e16   # = Ra / Thermal expansion coefficient
  end
end

This has some interesting implications. First, a higher Rayleigh number makes time scales correspondingly smaller; where we generated graphical output only once every 0.01 time units before, we now need to choose the corresponding increment smaller by a factor of 100:

subsection ??
  set ?? = velocity statistics, temperature statistics, ...
                       ...heat flux statistics, visualization

  subsection ??
    set ?? = 0.0001
  end
end

Secondly, a simulation like this – in 3d, with a significant number of cells, and for a significant number of time steps – will likely take a good amount of time. The computations for which we show results below was run using 64 processors by running it using the command mpirun -n 64 ./aspect convection-box-3d.prm. If the machine should crash during such a run, a significant amount of compute time would be lost if we had to run everything from the start. However, we can avoid this by periodically checkpointing the state of the computation:

subsection ??
set ?? = 50
end

If the computation does crash (or if a computation runs out of the time limit imposed by a scheduling system), then it can be restarted from such checkpointing files (see the parameter Resume computation in Section ??).

Running with this input file requires a bit of patience²¹ since the number of degrees of freedom is just so large: it starts with a bit over 330,000…

Running with 64 MPI tasks.
Number of active cells: 512 (on 4 levels)
Number of degrees of freedom: 20,381 (14,739+729+4,913)

*** Timestep 0:  t=0 seconds
   Solving temperature system... 0 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 18 iterations.

Number of active cells: 1,576 (on 5 levels)
Number of degrees of freedom: 63,391 (45,909+2,179+15,303)

*** Timestep 0:  t=0 seconds
   Solving temperature system... 0 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 19 iterations.

Number of active cells: 3,249 (on 5 levels)
Number of degrees of freedom: 122,066 (88,500+4,066+29,500)

*** Timestep 0:  t=0 seconds
   Solving temperature system... 0 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 20 iterations.

Number of active cells: 8,968 (on 5 levels)
Number of degrees of freedom: 331,696 (240,624+10,864+80,208)

*** Timestep 0:  t=0 seconds
   Solving temperature system... 0 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 21 iterations.
[...]

…but then increases quickly to around 2 million as the solution develops some structure and, after time $t = 0.003$ where we allow for an additional refinement, increases to over 10 million where it then hovers between 8 and 14 million with a maximum of 15,147,534. Clearly, even on a reasonably quick machine, this will take some time: running this on a machine bought in 2011, doing the 10,000 time steps to get to $t = 0.0219$ takes approximately 484,000 seconds (about five and a half days).

The structure or the solution is easiest to grasp by looking at isosurfaces of the temperature. This is shown in Fig. 13 and you can find a movie of the motion that ensues from the heating at the bottom at http://www.youtube.com/watch?v=_bKqU_P4j48. The simulation uses adaptively changing meshes that are fine in rising plumes and sinking blobs and are coarse where nothing much happens. This is most easily seen in the movie at http://www.youtube.com/watch?v=CzCKYyR-_cmg. Fig. 14 shows some of these meshes, though still pictures do not do the evolving nature of the mesh much justice. The effect of increasing the Rayleigh number is apparent when comparing the size of features with, for example, the picture at the right of Fig. 7. In contrast to that picture, the simulation is also clearly non-stationary.

Figure 13: Convection in a 3d box: Temperature isocontours and some velocity vectors at the first time step after times

t = 0.001, 0.004, 0.006

(top row, left to right) an

t = 0.01, 0.013, 0.018

(bottom row).

Figure 14: Convection in a 3d box: Meshes and large-scale velocity field for the third, fourth and sixth of the snapshots shown in Fig. 13.

As before, we could analyze all sorts of data from the statistics file but we will leave this to those interested in specific data. Rather, Fig. 15 only shows the upward heat flux through the bottom and top boundaries of the domain as a function of time.²² The figure reinforces a pattern that can also be seen by watching the movie of the temperature field referenced above, namely that the simulation can be subdivided into three distinct phases. The first phase corresponds to the initial overturning of the unstable layering of the temperature field and is associated with a large spike in heat flux as well as large velocities (not shown here). The second phase, until approximately $t = 0.01$ corresponds to a relative lull: some plumes rise up, but not very fast because the medium is now stably layered but not fully mixed. This can be seen in the relatively low heat fluxes, but also in the fact that there are almost horizontal temperature isosurfaces in the second of the pictures in Fig. 13. After that, the general structure of the temperature field is that the interior of the domain is well mixed with a mostly constant average temperature and thin thermal boundary layers at the top and bottom from which plumes rise and sink. In this regime, the average heat flux is larger but also more variable depending on the number of plumes currently active. Many other analyses would be possible by using what is in the statistics file or by enabling additional postprocessors.

Figure 15: Convection in a 3d box: Upward heat flux through the bottom and top boundaries as a function of time.

5.2.3 Convection in a box with prescribed, variable velocity boundary conditions

A similarly simple setup to the ones considered in the previous subsections is to equip the model we had with a different set of boundary conditions. There, we used slip boundary conditions, i.e., the fluid can flow tangentially along the four sides of our box but this tangential velocity is unspecified. On the other hand, in many situations, one would like to actually prescribe the tangential flow velocity as well. A typical application would be to use boundary conditions at the top that describe experimentally determined velocities of plates. This cookbook shows a simple version of something like this. To make it slightly more interesting, we choose a $2 \times 1$ domain in 2d.

Like for many other things, ASPECT has a set of plugins for prescribed velocity boundary values (see Sections ?? and 6.3.6). These plugins allow one to write sophisticated models for the boundary velocity on parts or all of the boundary, but there is also one simple implementation that just takes a formula for the components of the velocity.

To illustrate this, let us consider the cookbooks/platelike-_boundary.prm input file. It essentially extends the input file considered in the previous example. The part of this file that we are particularly interested in in the current context is the selection of the kind of velocity boundary conditions on the four sides of the box geometry, which we do using a section like this:

subsection ??
  set ?? = left, right, bottom
  set ?? = top: function

  subsection ??
    set ??      = x,z,t
    set ??  = pi=3.1415926
    set ?? = if(x>1+sin(0.5*pi*t), 1, -1); 0
  end
end

We use tangential flow at boundaries named left, right and bottom. Additionally, we specify a comma separated list (here with only a single element) of pairs consisting of the name of a boundary and the name of a prescribed velocity boundary model. Here, we use the function model on the top boundary, which allows us to provide a function-like notation for the components of the velocity vector at the boundary.

The second part we need is that we actually describe the function that sets the velocity. We do this in the subsection Function. The first of these parameters gives names to the components of the position vector (here, we are in 2d and we use $x$ and $z$ as spatial variable names) and the time. We could have left this entry at its default, x,y,t, but since we often think in terms of “depth” as the vertical direction, let us use z for the second coordinate. In the second parameter we define symbolic constants that can be used in the formula for the velocity that is specified in the last parameter. This formula needs to have as many components as there are space dimensions, separated by semicolons. As stated, this means that we prescribe the (horizontal) $x$ -velocity and set the vertical velocity to zero. The horizontal component is here either $1$ or $- 1$ , depending on whether we are to the right or the left of the point $1 + sin (π t ∕ 2)$ that is moving back and forth with time once every four time units. The if statement understood by the parser we use for these formulas has the syntax if(condition, value-if-true, value-if-false).

Note: While you can enter most any expression into the parser for these velocity boundary conditions, not all make sense. In particular, if you use an incompressible medium like we do here, then you need to make sure that either the flow you prescribe is indeed tangential, or that at least the flow into and out of the boundary this function applies to is balanced so that in sum the amount of material in the domain stays constant.

It is in general not possible for ASPECT to verify that a given input is sensible. However, you will quickly find out if it isn’t: The linear solver for the Stokes equations will simply not converge. For example, if your function expression in the input file above read
if(x>1+sin(0.5*pi*t), 1, -1); 1
then at the time of writing this you would get the following error message:
*** Timestep 0: t=0 seconds
Solving temperature system... 0 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system...

…some timing output …

––––––––––––––––––––––––––
Exception on processing:
Iterative method reported convergence failure in step 9539 with residual 6.0552
Aborting!
––––––––––––––––––––––––––

The reason is, of course, that there is no incompressible (divergence free) flow field that allows for a constant vertical outflow component along the top boundary without corresponding inflow anywhere else.

The remainder of the setup is described in the following, complete input file:

############### Global parameters

set ??                              = 2
set ??                             = 0
set ??                               = 20
set ?? = false
set ??                       = output-platelike-boundary

############### Parameters describing the model
# Let us here choose again a box domain of size 2x1
# where we fix the temperature at the bottom and top,
# allow free slip along the bottom, left and right,
# and prescribe the velocity along the top using the
# ‘function’ description.

subsection ??
  set ?? = box

  subsection ??
    set ?? = 2
    set ?? = 1
  end
end

# We then set the temperature to one at the bottom and zero
# at the top:
subsection ??
  set ?? = bottom, top
  set ?? = box

  subsection ??
    set ?? = 1
    set ??    = 0
  end
end

# The velocity along the top boundary models a spreading
# center that is moving left and right:
subsection ??
  set ?? = left, right, bottom
  set ?? = top: function

  subsection ??
    set ??      = x,z,t
    set ??  = pi=3.1415926
    set ?? = if(x>1+sin(0.5*pi*t), 1, -1); 0
  end
end

# We then choose a vertical gravity model and describe the
# initial temperature with a vertical gradient. The default
# strength for gravity is one. The material model is the
# same as before.
subsection ??
  set ?? = vertical
end

subsection ??
  set ?? = function

  subsection ??
    set ??      = x,z
    set ?? = (1-z)
  end
end

subsection ??
  set ?? = simple

  subsection ??
    set ??          = 1e-6
    set ?? = 1e-4
    set ??                     = 1
  end
end

# The final part of this input file describes how many times the
# mesh is refined and what to do with the solution once computed
subsection ??
  set ??        = 0
  set ??          = 5
  set ?? = 0
end

subsection ??
  set ?? = visualization, temperature statistics, heat flux statistics

  subsection ??
    set ?? = 0.1
  end
end

This model description yields a setup with a Rayleigh number of 200 (taking into account that the domain has size 2). It would, thus, be dominated by heat conduction rather than convection if the prescribed velocity boundary conditions did not provide a stirring action. Visualizing the results of this simulation²³ yields images like the ones shown in Fig. 16.

Figure 16: Variable velocity boundary conditions: Temperature and velocity fields at the initial time (top left) and at various other points in time during the simulation.

5.2.4 Using passive and active compositional fields

One frequently wants to track where material goes, either because one simply wants to see where stuff ends up (e.g., to determine if a particular model yields mixing between the lower and upper mantle) or because the material model in fact depends not only pressure, temperature and location but also on the mass fractions of certain chemical or other species. We will refer to the first case as passive and the latter as active to indicate the role of the additional quantities whose distribution we want to track. We refer to the whole process as compositional since we consider quantities that have the flavor of something that denotes the composition of the material at any given point.

There are basically two ways to achieve this: one can advect a set of particles (“tracers”) along with the velocity field, or one can advect along a field. In the first case, where the closest particle came from indicates the value of the concentration at any given position. In the latter case, the concentration(s) at any given position is simply given by the value of the field(s) at this location.

ASPECT implements both strategies, at least to a certain degree. In this cookbook, we will follow the route of advected fields.

The passive case. We will consider the exact same situation as in the previous section but we will ask where the material that started in the bottom 20% of the domain ends up, as well as the material that started in the top 20%. For the moment, let us assume that there is no material between the materials at the bottom, the top, and the middle. The way to describe this situation is to simply add the following block of definitions to the parameter file (you can find the full parameter file in cookbooks/composition-_passive.prm:

# This is the new part: We declare that there will
# be two compositional fields that will be
# advected along. Their initial conditions are given by
# a function that is one for the lowermost 0.2 height
# units of the domain and zero otherwise in the first case,
# and one in the top most 0.2 height units in the latter.
subsection ??
  set ?? = 2
end

subsection ??
  set ?? = function

  subsection ??
    set ??      = x,y
    set ?? = if(y<0.2, 1, 0) ; if(y>0.8, 1, 0)
  end
end

Running this simulation yields results such as the ones shown in Fig. 17 where we show the values of the functions $c_{1} (x, t)$ and $c_{2} (x, t)$ at various times in the simulation. Because these fields were one only inside the lowermost and uppermost parts of the domain, zero everywhere else, and because they have simply been advected along with the flow field, the places where they are larger than one half indicate where material has been transported to so far.²⁴

Figure 17: Passive compositional fields: The figures show, at different times in the simulation, the velocity field along with those locations where the first compositional field is larger than 0.5 (in red, indicating the locations where material from the bottom of the domain has gone) as well as where the second compositional field is larger than 0.5 (in blue, indicating material from the top of the domain. The results were obtained with two more global refinement steps compared to the cookbooks/composition-_passive.prm input file.

Figure 18: Passive compositional fields: A later image of the simulation corresponding to the sequence shown in Fig. 17 (left) and zoom-in on the center, also showing the mesh (right).

Fig. 17 shows one aspect of compositional fields that occasionally makes them difficult to use for very long time computations. The simulation shown here runs for 20 time units, where every cycle of the spreading center at the top moving left and right takes 4 time units, for a total of 5 such cycles. While this is certainly no short-term simulation, it is obviously visible in the figure that the interface between the materials has diffused over time. Fig. 18 shows a zoom into the center of the domain at the final time of the simulation. The figure only shows values that are larger than 0.5, and it looks like the transition from red or blue to the edge of the shown region is no wider than 3 cells. This means that the computation is not overly diffusive but it is nevertheless true that this method has difficulty following long and thin filaments.²⁵ This is an area in which ASPECT may see improvements in the future.

Figure 19: Passive compositional fields: Minimum and maximum of the first compositional variable over time, as well as the mass

m_{1} (t) = \int_{Ω} c_{1} (x, t)

stored in this variable.

A different way of looking at the quality of compositional fields as opposed to particles is to ask whether they conserve mass. In the current context, the mass contained in the $i$ th compositional field is $m_{i} (t) = \int_{Ω} c_{i} (x, t)$ . This can easily be achieve in the following way, by adding the composition statistics postprocessor:

subsection ??
set ?? = visualization, temperature statistics, composition statistics
end

While the scheme we use to advect the compositional fields is not strictly conservative, it is almost perfectly so in practice. For example, in the computations shown in this section (using two additional global mesh refinements over the settings in the parameter file cookbooks/composition-_passive.prm), Fig. 19 shows the maximal and minimal values of the first compositional fields over time, along with the mass $m_{1} (t)$ (these are all tabulated in columns of the statistics file, see Sections 4.1 and 4.4.2). While the maximum and minimum fluctuate slightly due to the instability of the finite element method in resolving discontinuous functions, the mass appears stable at a value of 0.403646 (the exact value, namely the area that was initially filled by each material, is 0.4; the difference is a result of the fact that we can’t exactly represent the step function on our mesh with the finite element space). In fact, the maximal difference in this value between time steps 1 and 500 is only $1.1 \cdot 1 0^{- 6}$ . In other words, these numbers show that the compositional field approach is almost exactly mass conservative.

The active case. The next step, of course, is to make the flow actually depend on the composition. After all, compositional fields are not only intended to indicate where material come from, but also to indicate the properties of this material. In general, the way to achieve this is to write material models where the density, viscosity, and other parameters depend on the composition, taking into account what the compositional fields actually denote (e.g., if they simply indicate the origin of material, or the concentration of things like olivine, perovskite, …). The construction of material models is discussed in much greater detail in Section 6.3.1; we do not want to revisit this issue here and instead choose – once again – the simplest material model that is implemented in ASPECT: the simple model.

The place where we are going to hook in a compositional dependence is the density. In the simple model, the density is fundamentally described by a material that expands linearly with the temperature; for small density variations, this corresponds to a density model of the form $ρ (T) = ρ_{0} (1 - α (T - T_{0}))$ . This is, by virtue of its simplicity, the most often considered density model. But the simple model also has a hook to make the density depend on the first compositional field $c_{1} (x, t)$ , yielding a dependence of the form $ρ (T) = ρ_{0} (1 - α (T - T_{0})) + γ c_{1}$ . Here, let us choose $ρ_{0} = 1, α = 0.01, T_{0} = 0, γ = 100$ . The rest of our model setup will be as in the passive case above. Because the temperature will be between zero and one, the temperature induced density variations will be restricted to 0.01, whereas the density variation by origin of the material is 100. This should make sure that dense material remains at the bottom despite the fact that it is hotter than the surrounding material.²⁶

This setup of the problem can be described using an input file that is almost completely unchanged from the passive case. The only difference is the use of the following section (the complete input file can be found in cookbooks/composition-_active.prm:

subsection ??
  set ?? = simple

  subsection ??
    set ??                           = 1e-6
    set ??                  = 0.01
    set ??                                      = 1
    set ??                              = 1
    set ??                          = 0
    set ?? = 0.1
  end
end

To debug the model, we will also want to visualize the density in our graphical output files. This is done using the following addition to the postprocessing section, using the density visualization plugin:

subsection ??
  set ?? = visualization, temperature statistics, composition statistics

  subsection ??
    set ?? = density
    set ?? = 0.1
  end
end

Figure 20: Active compositional fields: Compositional field 1 at the time

t = 0, 10, 20

. Compared to the results shown in Fig. 17 it is clear that the heavy material stays at the bottom of the domain now. The effect of the density on the velocity field is also clearly visible by noting that at all three times the spreading center at the top boundary is in exactly the same position; this would result in exactly the same velocity field if the density and temperature were constant.

Figure 21: Active compositional fields: Temperature fields at

t = 0, 2, 4, 8, 12, 20

. The black line is the isocontour line

c_{1} (x, t) = 0.5

delineating the position of the dense material at the bottom.

Results of this model are visualized in Fig.s 20 and 21. What is visible is that over the course of the simulation, the material that starts at the bottom of the domain remains there. This can only happen if the circulation is significantly affected by the high density material once the interface starts to become non-horizontal, and this is indeed visible in the velocity vectors. As a second consequence, if the material at the bottom does not move away, then there needs to be a different way for the heat provided at the bottom to get through the bottom layer: either there must be a secondary convection system in the bottom layer, or heat is simply conducted. The pictures in the figure seem to suggest that the latter is the case.

It is easy, using the outline above, to play with the various factors that drive this system, namely:

The magnitude of the velocity prescribed at the top.
The magnitude of the velocities induced by thermal buoyancy, as resulting from the magnitude of gravity and the thermal expansion coefficient.
The magnitude of the velocities induced by compositional variability, as described by the coefficient $γ$ and the magnitude of gravity.

Using the coefficients involved in these considerations, it is trivially possible to map out the parameter space to find which of these effects is dominant. As mentioned in discussing the values in the input file, what is important is the relative size of these parameters, not the fact that currently the density in the material at the bottom is 100 times larger than in the rest of the domain, an effect that from a physical perspective clearly makes no sense at all.

The active case with reactions. This section was contributed by Juliane Dannberg and René Gaßmöller.

In addition, there are setups where one wants the compositional fields to interact with each other. One example would be material upwelling at a mid-ocean ridge and changing the composition to that of oceanic crust when it reaches a certain depth. In this cookbook, we will describe how this kind of behaviour can be achieved by using the composition reaction function of the material model.

We will consider the exact same setup as in the previous paragraphs, except for the initial conditions and properties of the two compositional fields. There is one material that initially fills the bottom half of the domain and is less dense than the material above. In addition, there is another material that only gets created when the first material reaches the uppermost 20% of the domain, and that has a higher density. This should cause the first material to move upwards, get partially converted to the second material, which then sinks down again. This means we want to change the initial conditions for the compositional fields:

subsection ??
  set ?? = function

  subsection ??
    set ??      = x,z
    set ?? = if(z<0.5, 1, 0); 0
  end
end

Moreover, instead of the simple material model, we will use the composition reaction material model, which basically behaves in the same way, but can handle two active compositional field and a reaction between those two fields. In the input file, the user defines a depth and above this reaction depth the first compositional fields is converted to the second field. This can be done by changing the following section (the complete input file can be found in cookbooks/composition-_reaction.prm).

subsection ??
  set ?? = composition reaction

  subsection ??
    set ??                           = 1e-6
    set ??                  = 0.01
    set ??                                      = 1
    set ?? = -5
    set ?? = 5
    set ??                                 = 0.2
  end
end

Figure 22: Reaction between compositional fields: Temperature fields at

t = 0, 2, 4, 8, 12, 20

. The black line is the isocontour line

c_{1} (x, t) = 0.5

delineating the position of the material starting at the bottom and the white line is the isocontour line

c_{2} (x, t) = 0.5

delineating the position of the material that is created by the reaction.

Results of this model are visualized in Fig 22. What is visible is that over the course of the simulation, the material that starts at the bottom of the domain ascends, reaches the reaction depth and gets converted to the second material, which starts to sink down.

5.2.5 Using particles

Using compositional fields to trace where material has come from or is going to has many advantages from a computational point of view. For example, the numerical methods to advect along fields are well developed and we can do so at a cost that is equivalent to one temperature solve for each of the compositional fields. Unless you have many compositional fields, this cost is therefore relatively small compared to the overall cost of a time step. Another advantage is that the value of a compositional field is well defined at every point within the domain. On the other hand, compositional fields over time diffuse initially sharp interfaces, as we have seen in the images of the previous section.

At the same time, the geodynamics community has a history of using particles for this purpose. Historically, this may have been because it is conceptually simpler to advect along individual particles rather than whole fields, since it only requires an ODE integrator rather than the stabilization techniques necessary to advect fields. They also provide the appearance of no diffusion, though this is arguable. Leaving aside the debate whether fields or particles are the way to go, ASPECT supports both: using fields and using particles.

In order to advect particles along with the flow field, one just needs to add the particles postprocessor to the list of postprocessors and specify a few parameters. We do so in the cookbooks/composition-_passive-_particles.prm input file, which is otherwise just a minor variation of the cookbooks/composition-_passive.prm case discussed in the previous Section 5.2.4. In particular, the postprocess section now looks like this:

subsection ??
  set ?? = visualization, particles

  subsection ??
    set ?? = 0.1
  end

  subsection ??
    set ??        = 1000
    set ?? = 0.1
    set ??       = vtu
  end
end

The 1000 particles we are asking here are initially uniformly distributed throughout the domain and are, at the end of each time step, advected along with the velocity field just computed. (There are a number of options to decide which method to use for advecting particles, see Section ??.)

If you run this cookbook, information about all particles will be written into the output directory selected in the input file (as discussed in Section 4.1). In the current case, in addition to the files already discussed there, a directory listing at the end of a run will show several particle related files:

aspect> ls -l output/
total 932
-rw-rw-r-- 1 bangerth bangerth  11134 Dec 11 10:08 depth_average.gnuplot
-rw-rw-r-- 1 bangerth bangerth  11294 Dec 11 10:08 log.txt
-rw-rw-r-- 1 bangerth bangerth 326074 Dec 11 10:07 parameters.prm
-rw-rw-r-- 1 bangerth bangerth 577138 Dec 11 10:07 parameters.tex
drwxrwxr-x 2 bangerth bangerth   4096 Dec 11 18:40 particles
-rw-rw-r-- 1 bangerth bangerth    335 Dec 11 18:40 particles.pvd
-rw-rw-r-- 1 bangerth bangerth    168 Dec 11 18:40 particles.visit
drwxr-xr-x 2 bangerth bangerth   4096 Dec 11 10:08 solution
-rw-rw-r-- 1 bangerth bangerth    484 Dec 11 10:08 solution.pvd
-rw-rw-r-- 1 bangerth bangerth    451 Dec 11 10:08 solution.visit
-rw-rw-r-- 1 bangerth bangerth   8267 Dec 11 10:08 statistics

Here, the particles.pvd and particles.visit files contain a list of all visualization files from all processors and time steps. These files can be loaded in much the same way as the solution.pvd and solution.visit files that were discussed in Section 4.4. The actual data files – possibly a large number, but not of much immediate interest to users – are located in the particles subdirectory.

Coming back to the example at hand, we can visualize the particles that were advected along by opening both the field-based output files and the ones that correspond to particles (for example, output/solution.visit and output/particles.visit) and using a pseudo-color plot for the particles, selecting the “id” of particles to color each particle. By going to, for example, the output from the 72nd visualization output, this then results in a plot like the one shown in Fig. 23.

Figure 23: Passively advected quantities visualized through both a compositional field and a set of 1,000 particles, at

t = 7.2

The particles shown here are not too impressive in still pictures since they are colorized by their particle number, which does not carry any particular meaning other than the fact that it enumerates the particles.²⁷ The particle “id” can, however, be useful when viewing an animation of time steps. There, the different colors of particles allows the eye to follow the motion of a single particle. This is especially true if, after some time, particles have become well mixed by the flow field and adjacent particles no longer have similar colors. In any case, viewing such animations makes it rather intuitive to understand a flow field, but it can of course not be reproduced in a static medium such as this manual.

In any case, we will see in the next section how to attach more interesting information to particles, and how to visualize these.

Using particle properties. The particles in the above example only fulfil the purpose of visualizing the convection pattern. A more meaningful use for particles is to attach “properties” to them. A property consists of one or more numbers (or vectors or tensors) that may either be set at the beginning of the model run and stay constant, or are updated during the model runtime. These properties can then be used for many applications, e.g., storing an initial property (like the position, or initial composition), evaluating a property at a defined particle path (like the pressure-temperature evolution of a certain piece of rock), or by integrating a quantity along a particle path (like the integrated strain a certain domain has experienced). We illustrate these properties in the cookbook cookbooks/composition-_passive-_particles-_properties.prm, in which we add the following lines to the Particles subsection (we also increase the number of particles compared to the previous section to make the visualizations below more obvious):

subsection ??
  subsection ??
    set ??        = 50000

    set ?? = function, initial composition, initial position, pT path

    subsection ??
      set ??      = x,y
      set ?? = if(y<0.2, 1, 0)
    end
  end
end

These commands make sure that every particle will carry four different properties (function, pT path, initial position and initial composition), some of which may be scalars and others can have multiple components. (A full list of particle properties that can currently be selected can be found in Section ??, and new particle properties can be added as plugins as described in Section 6.2.) The properties selected above do the following:

initial position: This particle property simply stores the initial coordinates of the particle and then never changes them. This can be useful to compare the final position of a particle with its initial position and therefore determine how far certain domains traveled during the model runtime. Alternatively, one may want to simply visualize the norm of this vector property (i.e., the norm of the initial position, which is of course equal to the distance from the origin at which a particle started): in mantle simulations in spherical coordinates, the radius is indicative of which part of the mantle a particle comes from, and can therefore be used to visualize where material gets transported over the course of a simulation.
initial composition: This property uses the same method to initialize particle properties as is used to initialize the corresponding compositional fields. Using this, it stores the compositional field initialization values at the location where the particle started, and again never changes them. This is useful in the same context as shown for the field-based example in Section 5.2.4 where we would like to figure where materials ends up. In this case, one would set the initial composition to an indicator function for certain parts of the domain, and then set the initial composition property for the particles to match this composition. Letting the particles advect and at a later time visualizing this particle property will then show where particles came from. In cases where compositional variables undergo changes, e.g., by describing phase changes or chemical reactions, the “initial composition” property can also be useful to compare the final composition of a particle with its initial composition and therefore determine which regions underwent reactions such as those described in Section 5.2.4, and where the material that underwent this reaction got transported to.
function: This particle property can be used to assign to each particle values that are described based on a function of space. It provides an alternative way to set initial values if you don’t want to first set a compositional field’s initial values based on a function, and then copy these values via the “initial composition” property to particles. In the example above, we use the same function as for the compositional initial composition of field number one in Section 5.2.4. Therefore, this property should behave identical to the compositional field (except that the compositional field may have a reaction term that this particle property does not), although the two are of course advected using very different methods. This allows to compare the error in particle position to the numerical diffusion of the compositional field.
pT path: This property is interesting in that the particle property’s values always exactly mirror the pressure and temperature at the particle’s current location. This does not seem to be very useful since the information is already available. However, because each particle has a unique id, one can select a particular particle and output its properties (including pressure and temperature based on the pT path property) at all time steps. This allows for the creation of a pressure-temperature curve of a certain piece of rock. This property is interesting in many lithosphere and crustal scale models, because it is determining the metamorphic reactions that happen during deformation processes (e.g., in a subduction zone).

Figure 24: Passively advected particle properties visualized. Top row: Composition

C_{1}

and particle property “initial

C_{1}

”. The blue line in both figures is the 0.5-isocontour for the

C_{1}

field. Bottom row: Norm of the “initial position” of particles at

t = 0

and

t = 20

The results of all of these properties can of course be visualized. Fig. 24 shows some of the pictures one can create with particles. The top row shows both the composition field $C_{1}$ (along with the mesh on which it is defined) and the corresponding “initial $C_{1}$ ” particle property, at $t = 7.2$ . Because the compositional field does not undergo any reactions, it should of course simply be the initial composition advected along with the flow field, and therefore equal the values of the corresponding particle property. However, field-based compositions suffer from diffusion. On the other hand, the amount of diffusion can easily be decreased by mesh refinement.

The bottom of the figure shows the norm of the “initial position” property at the initial time and at time $t = 20$ . These images therefore show how far from the origin each of the particles shown was at the initial time.

Using active particles. In the examples above, particle properties passively track distinct model properties. These particle properties, however, may also be used to actively influence the model as it runs. For instance, a composition-dependent material model may use particles’ initial composition rather than an advected compositional field. To make this work – i.e., to get information from particles that are located at unpredictable locations, to the quadrature points at which material models and other parts of the code need to evaluate these properties – we need to somehow get the values from particles back to fields that can then be evaluated at any point where this is necessary. A slightly modified version of the active-composition cookbook (cookbooks/composition-_active.prm) illustrates how to use ‘active particles’ in this manner.

This cookbook, cookbooks/composition-_active-_particles.prm, modifies two sections of the input file. First, particles are added under the Postprocess section:

subsection ??
  subsection ??
    set ??         = 100000
    set ??  = 0
    set ??        = vtu
    set ?? = velocity, initial composition
    set ??      = cell average
    set ??   = random uniform
  end
end

Here, each particle will carry the velocity and initial composition properties. In order to use the particle initial composition value to modify the flow through the material model, we now modify the Composition section:

subsection ??
  set ??            = 2
  set ??             = lower, upper
  set ?? = particles, particles
  set ??  = lower:initial lower, upper:initial upper
end

What this does is the following: It says that there will be two compositional fields, called lower and upper (because we will use them to indicate material that comes from either the lower or upper part of the domain). Next, the Compositional field methods states that each of these fields will be computed by interpolation from the particles (if we had left this parameter at its default value, field, for each field, then it would have solved an advection PDE in each time step, as we have done in all previous examples).

In this case, we specify that both of the compositional fields are in fact interpolated from particle properties in each time step. How this is done is described in the fourth line. To understand it, it is important to realize that particles and fields have matching names: We have named the fields lower and upper, whereas the properties that result from the initial composition entry in the particles section are called initial lower and initial upper, since they inherit the names of the fields.

The syntax for interpolation from particles to fields then states that the lower field will be set to the interpolated value of the initial lower particle property at the end of each time step, and similarly for the upper field. In turn, the initial composition particle property was using the same method that one would have used for the compositional field initialization if these fields were actually advected along in each time step.

In this model the given global refinement level (5), associated number of cells (1024) and 100,000 total particles produces an average particle-per-cell count slightly below 100. While on the high end compared to most geodynamic studies using active particles, increasing the number of particles per cell further may alter the solution. As with the numerical resolution, any study using active particles should systematically vary the number of particles per cell in order to determine this parameter’s influence on the simulation.

Note: ASPECT’s particle implementation is in a preliminary state. While the accuracy and scalability of the implementation is benchmarked, other limitations remain. This in particular means that it is not optimized for performance, and more than a few thousand particles per process can slow down a model significantly. Moreover, models with a highly adaptive mesh and many particles do encounter a significant slowdown, because ASPECT only considers the number of degrees of freedom for load balancing across processes and not the number of particles. Therefore processes that compute the solution for coarse-grid regions have to process many more particles than other processes. Additionally, the checkpoint/restart functionality for particles is only implemented in models with a constant number of processes before and after the checkpoint and when the selected particle properties do not change. These limitations might be removed over time, but for current models the user should be aware of them.

5.2.6 Using a free surface

This section was contributed by Ian Rose.

Free surfaces are numerically challenging but can be useful for self consistently tracking dynamic topography and may be quite important as a boundary condition for tectonic processes like subduction. The parameter file cookbooks/free-_surface.prm provides a simple example of how to set up a model with a free surface, as well as demonstrates some of the challenges associated with doing so.

ASPECT supports models with a free surface using an Arbitrary Lagrangian-Eulerian framework (see Section 2.12). Most of this is done internally, so you do not need to worry about the details to run this cookbook. Here we demonstrate the evolution of surface topography that results when a blob of hot material rises in the mantle, pushing up the free surface as it does. Usually the amplitude of free surface topography will be small enough that it is difficult to see with the naked eye in visualizations, but the topography postprocessor can help by outputting the maximum and minimum topography on the free surface at every time step.

The bulk of the parameter file for this cookbook is similar to previous ones in this manual. We use initial temperature conditions that set up a hot blob of rock in the center of the domain.

The main addition is the Free surface subsection. In this subsection you need to give ASPECT a comma separated list of the free surface boundary indicators. In this case, we are dealing with the ‘top’ boundary of a box in 2D. There is another significant parameter that needs to be set here: the value for the stabilization parameter “theta”. If this parameter is zero, then there is no stabilization, and you are likely to see instabilities develop in the free surface. If this parameter is one then it will do a good job of stabilizing the free surface, but it may overly damp its motions. The default value is 0.5.

Also worth mentioning is the change to the CFL number. Stability concerns typically mean that when making a model with a free surface you will want to take smaller time steps. In general just how much smaller will depend on the problem at hand as well as the desired accuracy.

Following are the sections in the input file specific to this testcase. The full parameter file may be found at cookbooks/free-_surface.prm.

set ??                             = 0.1

subsection ??
  set ?? = function
  subsection ??
    set ??      = x,y
    set ?? = if( sqrt( (x-250.e3)^2 + (y-100.e3)^2 ) < 25.e3, 200.0, 0.0)
  end
end

subsection ??
  set ??        = top
  set ?? = 0.5
end

subsection ??
  set ??   = left, right, bottom, top
end

subsection ??
  set ?? = left, right, bottom
end

subsection ??
  set ?? = visualization,topography,velocity statistics,
  subsection ??
    set ?? = 1.e6
  end
end

Running this input file will produce results like those in Figure 25. The model starts with a single hot blob of rock which rises in the domain. As it rises, it pushes up the free surface in the middle, creating a topographic high there. This is similar to the kind of dynamic topography that you might see above a mantle plume on Earth. As the blob rises and diffuses, it loses the buoyancy to push up the boundary, and the surface begins to relax.

After running the cookbook, you may modify it in a number of ways:

Add a more complicated initial temperature field to see how that affects topography.
Add a high-viscosity lithosphere to the top using a compositional field to tamp down on topography.
Explore different values for the stabilization theta and the CFL number to understand the nature of when and why stabilization is necessary.
Try a model in a different geometry, such as spherical shells.

Figure 25: Evolution of surface topography due to a rising blob. On the left is a snapshot of the model setup. The right shows the value of the highest topography in the domain over 18 Myr of model time. The topography peaks at 165 meters after 5.2 Myr. This cookbook may be run with the cookbooks/free-_surface.prm input file.

5.2.7 Using a free surface in a model with a crust

This section was contributed by William Durkin.

This cookbook is a modification of the previous example that explores changes in the way topography develops when a highly viscous crust is added. In this cookbook, we use a material model in which the material changes from low viscosity mantle to high viscosity crust at $z = z_{j} = jump height$ , i.e., the piecewise viscosity function is defined as

\begin{array}{l} η (z) = \{\begin{matrix} η_{U} & for z > z_{j}, \\ η_{L} & for z \leq z_{j} . \end{matrix} \end{array}

where $η_{U}$ and $η_{L}$ are the viscosities of the upper and lower layers, respectively. This viscosity model can be implemented by creating a plugin that is a small modification of the simpler material model (from which it is otherwise simply copied). We call this material model “SimplerWithCrust”. In particular, what is necessary is an evaluation function that looks like this:

    template <int dim>
    void
    SimplerWithCrust<dim>::
    evaluate(const typename Interface<dim>::MaterialModelInputs &in,
              typename Interface<dim>::MaterialModelOutputs &out ) const
    {
      for (unsigned int i=0; i<in.position.size(); ++i)
        {
          const double z = in.position[i][1];

          if (z>jump_height)
            out.viscosities[i] = eta_U;
          else
            out.viscosities[i] = eta_L;

          out.densities[i] = reference_rho*(1.0-thermal_alpha*(in.temperature[i]-reference_T));
          out.thermal_expansion_coefficients[i] = thermal_alpha;
          out.specific_heat[i] = reference_specific_heat;
          out.thermal_conductivities[i] = k_value;
          out.compressibilities[i] = 0.0;
        }
    }

Additional changes make the new parameters Jump height, Lower viscosity, and Upper viscosity available to the input parameter file, and corresponding variables available in the class and used in the code snippet above. The entire code can be found in cookbooks/free-_surface-_with-_crust/plugin/simpler-_with-_crust.cc. Refer to Section 6.1 for more information about writing and running plugins.

The following changes are necessary compared to the input file from the cookbook shown in Section 5.2.6 to include a crust:

Load the plugin implementing the new material model:

set ?? = ./plugin/libsimpler-with-crust.so
Declare values for the new parameters:

subsection ??
  set ?? = simpler with crust
  subsection ??
    set ??             = 3300
    set ??       = 1250
    set ??         = 0.0
    set ??          = 1.0 # low thermal conductivity for a sharp blob
    set ?? = 4e-5

    # Parameters added for this cookbook:
    # The box is 200km high and has its origin set at the bottom left corner.
    # Setting the jump height to 170km creates a 30km thick crust
    set ??               = 1e20
    set ??               = 1e23
    set ??                   = 170e3
  end
end

Note that the height of the interface at 170km is interpreted in the coordinate system in which the box geometry of this cookbook lives. The box has dimensions $500 km \times 200 km$ , so an interface height of 170km implies a depth of 30km.

The entire script is located in cookbooks/free-_surface-_with-_crust/free-_surface-_with-_crust.prm.

Running this input file yields a crust that is 30km thick and 1000 times as viscous as the lower layer. Figure 26 shows that adding a crust to the model causes the maximum topography to both decrease and occur at a later time. Heat flows through the system primarily by advection until the temperature anomaly reaches the base of the crustal layer (approximately at the time for which Fig 26 shows the temperature profile). The crust’s high viscosity reduces the temperature anomaly’s velocity substantially, causing it to affect the surface topography at a later time. Just as the cookbook shown in Section 5.2.6, the topography returns to zero after some time.

Figure 26: Adding a viscous crust to a model with surface topography. The thermal anomaly spreads horizontally as it collides with the highly viscous crust (left). The addition of a crustal layer both dampens and delays the appearance of the topographic maximum and minimum (right).

5.2.8 Averaging material properties

The original motivation for the functionality discussed here, as well as the setup of the input file, were provided by Cedric Thieulot.

Geophysical models are often characterized by abrupt and large jumps in material properties, in particular in the viscosity. An example is a subducting, cold slab surrounded by the hot mantle: Here, the strong temperature-dependence of the viscosity will lead to a sudden jump in the viscosity between mantle and slab. The length scale over which this jump happens will be a few or a few tens of kilometers. Such length scales cannot be adequately resolved in three-dimensional computations with typical meshes for global computations.

Having large viscosity variations in models poses a variety of problems to numerical computations. First, you will find that they lead to very long compute times because our solvers and preconditioners break down. This may be acceptable if it would at least lead to accurate solutions, but large viscosity gradients lead also to large pressure gradients, and this in turn leads to over- and undershoots in the numerical approximation of the gradient. We will demonstrate both of these issues experimentally below.

One of the solution to such problems is the realization that one can mitigate some of the effects by averaging material properties on each cell somehow (see, for example, [SBE $^{+}$ 08, DK08, DMGT11, Thi15, TMK14]). Before going into detail, it is important to realize that if we choose material properties not per quadrature point when doing the integrals for forming the finite element matrix, but per cell, then we will lose accuracy in the solution in those cases where the solution is smooth. More specifically, we will likely lose one or more orders of convergence. In other words, it would be a bad idea to do this averaging unconditionally. On the other hand, if the solution has essentially discontinuous gradients and kinks in the velocity field, then at least at these locations we cannot expect a particularly high convergence order anyway, and the averaging will not hurt very much either. In cases where features of the solution that are due to strongly varying viscosities or other parameters, dominate, we may then as well do the averaging per cell.

To support such cases, ASPECT supports an operation where we evaluate the material model at every quadrature point, given the temperature, pressure, strain rate, and compositions at this point, and then either (i) use these values, (ii) replace the values by their arithmetic average $\bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_{i}$ , (iii) replace the values by their harmonic average $\bar{x} = {(\frac{1}{N} \sum_{i = 1}^{N} \frac{1}{x_{i}})}^{- 1}$ , (iv) replace the values by their geometric average $\bar{x} = {(\prod_{i = 1}^{N} \frac{1}{x_{i}})}^{- 1 ∕ N}$ , or (v) replace the values by the largest value over all quadrature points on this cell. Option (vi) is to project the values from the quadrature points to a bi- (in 2d) or trilinear (in 3d) $Q_{1}$ finite element space on every cell, and then evaluate this finite element representation again at the quadrature points. Unlike the other five operations, the values we get at the quadrature points are not all the same here.

We do this operation for all quantities that the material model computes, i.e., in particular, the viscosity, the density, the compressibility, and the various thermal and thermodynamic properties. In the first 4 cases, the operation guarantees that the resulting material properties are bounded below and above by the minimum and maximum of the original data set. In the last case, the situation is a bit more complicated: The nodal values of the $Q_{1}$ projection are not necessarily bounded by the minimal or maximal original values at the quadrature points, and then neither are the output values after re-interpolation to the quadrature points. Consequently, after projection, we limit the nodal values of the projection to the minimal and maximal original values, and only then interpolate back to the quadrature points.

We demonstrate the effect of all of this with the “sinker” benchmark. This benchmark is defined by a high-viscosity, heavy sphere at the center of a two-dimensional box. This is achieved by defining a compositional field that is one inside and zero outside the sphere, and assigning a compositional dependence to the viscosity and density. We run only a single time step for this benchmark. This is all modeled in the following input file that can also be found in cookbooks/sinker-_with-_averaging/sinker-_with-_averaging.prm:

set ??                              = 2
set ??                             = 0
set ??                               = 0
set ??                       = output_sinker-with-averaging

set ??                 = volume

subsection ??
  set ?? = box
  subsection ??
    set ??  = 1.0000
    set ??  = 1.0000
  end
end

subsection ??
  set ??       = left, right, bottom, top
end

subsection ??
  set ?? = simple

  subsection ??
    set ??             = 1
    set ??                     = 1
    set ?? = 0
    set ?? = 1e6
    set ?? = 10
  end

  set ?? = none
end

subsection ??
  set ?? = vertical
  subsection ??
    set ?? = 1
  end
end

############### Parameters describing the temperature field
# Note: The temperature plays no role in this model

subsection ??
  set ?? = box
end

subsection ??
  set ?? = function
  subsection ??
    set ?? = 0
  end
end

############### Parameters describing the compositional field
# Note: The compositional field is what drives the flow
# in this example

subsection ??
  set ?? = 1
end

subsection ??
  set ?? = function
  subsection ??
    set ??      = x,y
    set ?? = if( (sqrt((x-0.5)^2+(y-0.5)^2)>0.22) , 0 , 1 )
  end
end

############### Parameters describing the discretization

subsection ??
  set ??          = 6
  set ??        = 0
end

############### Parameters describing what to do with the solution

subsection ??
  set ?? = visualization, velocity statistics, composition statistics
  subsection ??
    set ??                 = vtu
    set ?? = 0
    set ?? = density, viscosity
  end
end

The type of averaging on each cell is chosen using this part of the input file:

subsection ??
set ?? = harmonic average
end

For the various different averaging options, and for different levels of mesh refinement, Fig. 27 shows pressure plots that illustrate the problem with oscillations of the discrete pressure. The important part of these plots is not that the solution looks discontinuous – in fact, the exact solution is discontinuous at the edge of the circle²⁸ – but the spikes that go far above and below the “cliff” in the pressure along the edge of the circle. Without averaging, these spikes are obviously orders of magnitude larger than the actual jump height. The spikes do not disappear under mesh refinement nor averaging, but they become far less pronounced with averaging. The results shown in the figure do not really allow to draw conclusions as to which averaging approach is the best; a discussion of this question can also be found in [SBE $^{+}$ 08, DK08, DMGT11, TMK14]).

A very pleasant side effect of averaging is that not only does the solution become better, but it also becomes cheaper to compute. Table 2 shows the number of outer GMRES iterations when solving the Stokes equations (1)–(2).²⁹ The implication of these results is that the averaging gives us a solution that not only reduces the degree of pressure over- and undershoots, but is also significantly faster to compute: for example, the total run time for 8 global refinement steps is reduced from 5,250s for no averaging to 358s for harmonic averaging.


$[- 45.2, 45.2]$	$[- 2.67, 2.67]$	$[- 3.58, 3.58]$	$[- 3.57, 3.57]$	$[- 1.80, 1.80]$	$[- 2.77, 2.77]$


$[- 44.5, 44.5]$	$[- 5.18, 5.18]$	$[- 5.09, 5.09]$	$[- 5.18, 5.18]$	$[- 5.20, 5.20]$	$[- 7.99, 7.99]$

Figure 27: Visualization of the pressure field for the “sinker” problem. Left to right: No averaging, arithmetic averaging, harmonic averaging, geometric averaging, pick largest, project to

Q_{1}

. Top: 7 global refinement steps. Bottom: 8 global refinement steps. The minimal and maximal pressure values are indicated below every picture. This range is symmetric because we enforce that the average of the pressure equals zero. The color scale is adjusted to show only values between

p = - 3

and

p = 3


# of global	no averaging	arithmetic	harmonic	geometric	pick	project
refinement steps		averaging	averaging	averaging	largest	to $Q_{1}$

4	30+64	30+13	30+10	30+12	30+13	30+15
5	30+87	30+14	30+13	30+14	30+14	30+16
6	30+171	30+14	30+15	30+14	30+15	30+17
7	30+143	30+27	30+28	30+26	30+26	30+28
8	30+188	30+27	30+26	30+27	30+28	30+28

Table 2: Number of outer GMRES iterations to solve the Stokes equations for various numbers of global mesh refinement steps and for different material averaging operations. The GMRES solver first tries to run 30 iterations with a cheaper preconditioner before switching to a more expensive preconditioner (see Section ??).

Such improvements carry over to more complex and realistic models. For example, in a simulation of flow under the East African Rift by Sarah Stamps, using approximately 17 million unknowns and run on 64 processors, the number of outer and inner iterations is reduced from 169 and 114,482 without averaging to 77 and 23,180 with harmonic averaging, respectively. This translates into a reduction of run-time from 145 hours to 17 hours. Assessing the accuracy of the answers is of course more complicated in such cases because we do not know the exact solution. However, the results without and with averaging do not differ in any significant way.

A final comment is in order. First, one may think that the results should be better in cases of discontinuous pressures if the numerical approximation actually allowed for discontinuous pressures. This is in fact possible: We can use a finite element in which the pressure space contains piecewise constants (see Section ??). To do so, one simply needs to add the following piece to the input file:

subsection ??
set ?? = true
end

Disappointingly, however, this makes no real difference: the pressure oscillations are no better (maybe even worse) than for the standard Stokes element we use, as shown in Fig. 28 and Table 3. Furthermore, as shown in Table 4, the iteration numbers are also largely unaffected if any kind of averaging is used – though they are far worse using the locally conservative discretization if no averaging has been selected. On the positive side, the visualization of the discontinuous pressure finite element solution makes it much easier to see that the true pressure is in fact discontinuous along the edge of the circle.

Figure 28: Visualization of the pressure field for the “sinker” problem. Like Fig. 27 but using the locally conservative, enriched Stokes element. Pressure values are shown in Table 3.


# of global	no averaging	arithmetic	harmonic	geometric	pick	project
refinement steps		averaging	averaging	averaging	largest	to $Q_{1}$

4	66.32	2.66	2.893	1.869	3.412	3.073
5	81.06	3.537	4.131	3.997	3.885	3.991
6	75.98	4.596	4.184	4.618	4.568	5.093
7	84.36	4.677	5.286	4.362	4.635	5.145
8	83.96	5.701	5.664	4.686	5.524	6.42

Table 3: Maximal pressure values for the “sinker” benchmark, using the locally conservative, enriched Stokes element. The corresponding pressure solutions are shown in Fig. 28.


# of global	no averaging	arithmetic	harmonic	geometric	pick	project
refinement steps		averaging	averaging	averaging	largest	to $Q_{1}$

4	30+376	30+16	30+12	30+14	30+14	30+17
5	30+484	30+16	30+14	30+14	30+14	30+16
6	30+583	30+16	30+17	30+14	30+17	30+17
7	30+1319	30+27	30+28	30+26	30+28	30+28
8	30+1507	30+28	30+27	30+28	30+28	30+29

Table 4: Like Table 2, but using the locally conservative, enriched Stokes element.

5.2.9 Prescribed internal velocity constraints

This section was contributed by Jonathan Perry-Houts

In cases where it is desirable to investigate the behavior of one part of the model domain but the controlling physics of another part is difficult to capture, such as corner flow in subduction zones, it may be useful to force the desired behavior in some parts of the model domain and solve for the resulting flow everywhere else. This is possible through the use of ASPECT’s “signal” mechanism, as documented in Section 6.5.

Internally, ASPECT adds “constraints” to the finite element system for boundary conditions and hanging nodes. These are places in the finite element system where certain solution variables are required to match some prescribed value. Although it is somewhat mathematically inadmissible to prescribe constraints on nodes inside the model domain, $Ω$ , it is nevertheless possible so long as the prescribed velocity field fits in to the finite element’s solution space, and satisfies the other constraints (i.e., is divergence free).

Using ASPECT’s signals mechanism, we write a shared library which provides a “slot” that listens for the signal which is triggered after the regular model constraints are set, but before they are “distributed.”

As an example of this functionality, below is a plugin which allows the user to prescribe internal velocities with functions in a parameter file:

#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/parsed_function.h>
#include <deal.II/fe/fe_values.h>
#include <aspect/global.h>
#include <aspect/simulator_signals.h>

namespace aspect
{
  using namespace dealii;

  // Global variables (to be set by parameters)
  bool prescribe_internal_velocities;

  // Because we do not initially know what dimension we’re in, we need
  // function parser objects for both 2d and 3d.
  Functions::ParsedFunction<2> prescribed_velocity_indicator_function_2d (2);
  Functions::ParsedFunction<3> prescribed_velocity_indicator_function_3d (3);
  Functions::ParsedFunction<2> prescribed_velocity_function_2d (2);
  Functions::ParsedFunction<3> prescribed_velocity_function_3d (3);

  /**
   * Declare additional parameters.
   */
  void declare_parameters(const unsigned int dim,
                          ParameterHandler &prm)
  {
    prm.declare_entry ("Prescribe internal velocities", "false",
                       Patterns::Bool (),
                       "Whether or not to use any prescribed internal velocities. "
                       "Locations in which to prescribe velocities are defined "
                       "in section ‘‘Prescribed velocities/Indicator function’’ "
                       "and the velocities are defined in section ‘‘Prescribed "
                       "velocities/Velocity function’’. Indicators are evaluated "
                       "at the center of each cell, and all DOFs associated with "
                       "the specified velocity component at the indicated cells "
                       "are constrained."
                      );

    prm.enter_subsection ("Prescribed velocities");
    {
      prm.enter_subsection ("Indicator function");
      {
        if (dim == 2)
          Functions::ParsedFunction<2>::declare_parameters (prm, 2);
        else
          Functions::ParsedFunction<3>::declare_parameters (prm, 3);
      }
      prm.leave_subsection ();

      prm.enter_subsection ("Velocity function");
      {
        if (dim == 2)
          Functions::ParsedFunction<2>::declare_parameters (prm, 2);
        else
          Functions::ParsedFunction<3>::declare_parameters (prm, 3);
      }
      prm.leave_subsection ();
    }
    prm.leave_subsection ();
  }

  template <int dim>
  void parse_parameters(const Parameters<dim> &,
                        ParameterHandler &prm)
  {
    prescribe_internal_velocities = prm.get_bool ("Prescribe internal velocities");
    prm.enter_subsection ("Prescribed velocities");
    {
      prm.enter_subsection("Indicator function");
      {
        try
          {
            if (dim == 2)
              prescribed_velocity_indicator_function_2d.parse_parameters (prm);
            else
              prescribed_velocity_indicator_function_3d.parse_parameters (prm);
          }
        catch (...)
          {
            std::cerr << "ERROR: FunctionParser failed to parse\n"
                      << "\t’Prescribed velocities.Indicator function’\n"
                      << "with expression\n"
                      << "\t’" << prm.get("Function expression") << "’";
            throw;
          }
      }
      prm.leave_subsection();

      prm.enter_subsection("Velocity function");
      {
        try
          {
            if (dim == 2)
              prescribed_velocity_function_2d.parse_parameters (prm);
            else
              prescribed_velocity_function_3d.parse_parameters (prm);
          }
        catch (...)
          {
            std::cerr << "ERROR: FunctionParser failed to parse\n"
                      << "\t’Prescribed velocities.Velocity function’\n"
                      << "with expression\n"
                      << "\t’" << prm.get("Function expression") << "’";
            throw;
          }
      }
      prm.leave_subsection();
    }
    prm.leave_subsection ();
  }

  /**
   * This function retrieves the unit support points (in the unit cell) for the current element.
   * The DGP element used when ’set Use locally conservative discretization = true’ does not
   * have support points. If these elements are in use, a fictitious support point at the cell
   * center is returned for each shape function that corresponds to the pressure variable,
   * whereas the support points for the velocity are correct; the fictitious points don’t matter
   * because we only use this function when interpolating the velocity variable, and ignore the
   * evaluation at the pressure support points.
   */
  template <int dim>
  std::vector< Point<dim> >
  get_unit_support_points_for_velocity(const SimulatorAccess<dim> &simulator_access)
  {
    std::vector< Point<dim> > unit_support_points;
    if ( !simulator_access.get_parameters().use_locally_conservative_discretization )
      {
        return simulator_access.get_fe().get_unit_support_points();
      }
    else
      {
        //special case for discontinuous pressure elements, which lack unit support points
        std::vector< Point<dim> > unit_support_points;
        const unsigned int dofs_per_cell = simulator_access.get_fe().dofs_per_cell;
        for (unsigned int dof=0; dof < dofs_per_cell; ++dof)
          {
            // base will hold element, base_index holds node/shape function within that element
            const unsigned int
            base       = simulator_access.get_fe().system_to_base_index(dof).first.first,
            base_index = simulator_access.get_fe().system_to_base_index(dof).second;
            // get the unit support points for the relevant element
            std::vector< Point<dim> > my_support_points = simulator_access.get_fe().base_element(base).get_unit_support_points();
            if ( my_support_points.size() == 0 )
              {
                //manufacture a support point, arbitrarily at cell center
                if ( dim==2 )
                  unit_support_points.push_back( Point<dim> (0.5,0.5) );
                if ( dim==3 )
                  unit_support_points.push_back( Point<dim> (0.5,0.5,0.5) );
              }
            else
              {
                unit_support_points.push_back(my_support_points[base_index]);
              }
          }
        return unit_support_points;
      }
  }

  namespace
  {
    template <int dim>
    const Point<2> &as_2d(const Point<dim> &p);

    template <>
    const Point<2> &as_2d(const Point<2> &p)
    {
      return p;
    }

    template <int dim>
    const Point<3> &as_3d(const Point<dim> &p);

    template <>
    const Point<3> &as_3d(const Point<3> &p)
    {
      return p;
    }

  }

  /**
   * This function is called by a signal which is triggered after the other constraints
   * have been calculated. This enables us to define additional constraints in the mass
   * matrix on any arbitrary degree of freedom in the model space.
   */
  template <int dim>
  void constrain_internal_velocities (const SimulatorAccess<dim> &simulator_access,
                                      ConstraintMatrix &current_constraints)
  {
    if (prescribe_internal_velocities)
      {
        const std::vector< Point<dim> > points = get_unit_support_points_for_velocity(simulator_access);
        const Quadrature<dim> quadrature (points);
        FEValues<dim> fe_values (simulator_access.get_fe(), quadrature, update_q_points);
        typename DoFHandler<dim>::active_cell_iterator cell;

        // Loop over all cells
        for (cell = simulator_access.get_dof_handler().begin_active();
             cell != simulator_access.get_dof_handler().end();
             ++cell)
          if (! cell->is_artificial())
            {
              fe_values.reinit (cell);
              std::vector<unsigned int> local_dof_indices(simulator_access.get_fe().dofs_per_cell);
              cell->get_dof_indices (local_dof_indices);

              for (unsigned int q=0; q<quadrature.size(); q++)
                // If it’s okay to constrain this DOF
                if (current_constraints.can_store_line(local_dof_indices[q]) &&
                    !current_constraints.is_constrained(local_dof_indices[q]))
                  {
                    // Get the velocity component index
                    const unsigned int c_idx =
                      simulator_access.get_fe().system_to_component_index(q).first;

                    // If we’re on one of the velocity DOFs
                    if ((c_idx >=
                         simulator_access.introspection().component_indices.velocities[0])
                        &&
                        (c_idx <=
                         simulator_access.introspection().component_indices.velocities[dim-1]))
                      {
                        // Which velocity component is this DOF associated with?
                        const unsigned int component_direction
                          = (c_idx
                             - simulator_access.introspection().component_indices.velocities[0]);

                        // we get time passed as seconds (always) but may want
                        // to reinterpret it in years
                        double time = simulator_access.get_time();
                        if (simulator_access.convert_output_to_years())
                          time /= year_in_seconds;

                        prescribed_velocity_indicator_function_2d.set_time (time);
                        prescribed_velocity_indicator_function_3d.set_time (time);
                        prescribed_velocity_function_2d.set_time (time);
                        prescribed_velocity_function_3d.set_time (time);

                        const Point<dim> p = fe_values.quadrature_point(q);

                        // Because we defined and parsed our parameter
                        // file differently for 2d and 3d we need to
                        // be sure to query the correct object for
                        // function values. The function parser
                        // objects expect points of a certain
                        // dimension, but Point p will be compiled for
                        // both 2d and 3d, so we need to do some trickery
                        // to make this compile.
                        double indicator, u_i;
                        if (dim == 2)
                          {
                            indicator = prescribed_velocity_indicator_function_2d.value
                                        (as_2d(p),
                                         component_direction);
                            u_i       = prescribed_velocity_function_2d.value
                                        (as_2d(p),
                                         component_direction);
                          }
                        else
                          {
                            indicator = prescribed_velocity_indicator_function_3d.value
                                        (as_3d(p),
                                         component_direction);
                            u_i       = prescribed_velocity_function_3d.value
                                        (as_3d(p),
                                         component_direction);
                          }

                        if (indicator > 0.5)
                          {
                            // Add a constraint of the form dof[q] = u_i
                            // to the list of constraints.
                            current_constraints.add_line (local_dof_indices[q]);
                            current_constraints.set_inhomogeneity (local_dof_indices[q], u_i);
                          }
                      }
                  }
            }
      }
  }

  // Connect declare_parameters and parse_parameters to appropriate signals.
  void parameter_connector ()
  {
    SimulatorSignals<2>::declare_additional_parameters.connect (&declare_parameters);
    SimulatorSignals<3>::declare_additional_parameters.connect (&declare_parameters);

    SimulatorSignals<2>::parse_additional_parameters.connect (&parse_parameters<2>);
    SimulatorSignals<3>::parse_additional_parameters.connect (&parse_parameters<3>);
  }

  // Connect constraints function to correct signal.
  template <int dim>
  void signal_connector (SimulatorSignals<dim> &signals)
  {
    signals.post_constraints_creation.connect (&constrain_internal_velocities<dim>);
  }

  // Tell Aspect to send signals to the connector functions
  ASPECT_REGISTER_SIGNALS_PARAMETER_CONNECTOR(parameter_connector)
  ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, signal_connector<3>)
}

The above plugin can be compiled with cmake . && make in the cookbooks/prescribed_velocity directory. It can be loaded in a parameter file as an “Additional shared library.” By setting parameters like those shown below, it is possible to produce many interesting flow fields such as the ones visualized in (Figure 29).

# Load the signal library.
set ?? = ./libprescribed_velocity.so

## Turn prescribed velocities on
set ?? = true

subsection ??
  subsection ??
    set ?? = x,y,t
    # Return where to prescribe u_x; u_y; u_z
    # (last one only used if dimension = 3)
    # 1 if velocity should be prescribed, 0 otherwise
    set ?? = if((x-.5)^2+(y-.5)^2<.125,1,0); \
                                if((x-.5)^2+(y-.5)^2<.125,1,0)
  end
  subsection ??
    set ?? = x,y,t
    # Return u_x; u_y; u_z (u_z only used if in 3d)
    set ?? = 1;-1
  end
end

(a)

(b)

Figure 29: Examples of flows with prescribed internal velocities, as described in Section 5.2.9.

5.2.10 Artificial viscosity smoothing

This section was contributed by Ryan Grove

Standard finite element discretizations of advection-diffusion equations introduce unphysical oscillations around steep gradients. Therefore, stabilization must be added to the discrete formulation to obtain correct solutions. In ASPECT, we use the Entropy Viscosity scheme developed by Guermond et al. in the paper [Jea11]. In this scheme, an artificial viscosity is calculated on every cell and used to try to combat these oscillations that cause unwanted overshoot and undershoot. More information about how ASPECT does this is located at https://dealii.org/developer/doxygen/deal.II/step_31.html.

Instead of just looking at an individual cell’s artificial viscosity, improvements in the minimizing of the oscillations can be made by smoothing. Smoothing is the act of finding the maximum artificial viscosity taken over a cell $T$ and the neighboring cells across the faces of $T$ , i.e.,

\bar{v_{h}} (T) = max_{K \in N (T)} v_{h} (K)

where $N (T)$ is the set containing $T$ and the neighbors across the faces of $T$ .

This feature can be turned on by setting the ?? flag inside the ?? subsection inside the ?? subsection in your parameter file.

To show how this can be used in practice, let us consider the simple convection in a quarter of a 2d annulus cookbook in Section 5.3.1, a radial compositional field was added to help show the advantages of using the artificial viscosity smoothing feature.

By applying the following changes shown below to the parameters of the already existing file

cookbooks/shell_simple_2d.prm,

subsection ??
  set ?? = 2

  subsection ??
    set ?? = true
  end
end

subsection ??
  set ?? = 1
end

subsection ??
  set ?? = function

  subsection ??
    set ??      = x,y
    set ?? = if(sqrt(x*x+y*y)<4000000,1,0)
  end
end

it is possible to produce pictures of the simple convection in a quarter of a 2d annulus such as the ones visualized in Figure 30.

(a)

(b)

Figure 30: Artificial viscosity smoothing: Example of the output of two similar runs. The run on the left has the artificial viscosity smoothing turned on and the run on the right does not, as described in Section 5.2.10.

5.2.11 Tracking finite strain

This section was contributed by Juliane Dannberg and Rene Gassmöller

Note: In this section, following [BKEO03, DT98], we denote the velocity gradient tensor as

G

, where

G = \nabla u^{T}

, and

u

is the velocity. Note that this is different from the definition of the strain rate

ϵ (u)

, which only contains the symmetric part of

G

. We then denote the deformation gradient (or deformation) tensor by

F

, where

F

is the tensor that deforms an initial state

x

into an deformed state

r = F x

In many geophysical settings, material properties, and in particular the rheology, do not only depend on the current temperature, pressure and strain rate, but also on the history of the system. This can be incorporated in ASPECT models by tracking history variables through compositional fields. In this cookbook, we will show how to do this by tracking the strain that idealized little grains of finite size accumulate over time at every (Lagrangian) point in the model.

Here, we use a material model plugin that defines the compositional fields as the components of the deformation gradient tensor $F_{i j}$ , and modifies the right-hand side of the corresponding advection equations to accumulate strain over time. This is done by adjusting the out.reaction_terms variable:

for (unsigned int q=0; q < in.position.size(); ++q)
  {
    // Convert the compositional fields into the tensor quantity they represent.
    Tensor<2,dim> strain;
    for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
      strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = in.composition[q][i];

    // Compute the strain accumulated in this timestep.
    const Tensor<2,dim> strain_increment = this->get_timestep() * (velocity_gradients[q] * strain);

    // Output the strain increment component-wise to its respective compositional field’s reaction terms.
    for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
      out.reaction_terms[q][i] = strain_increment[Tensor<2,dim>::unrolled_to_component_indices(i)];
  }

Let us denote the accumulated deformation at time step $n$ as $F^{n}$ . We can calculate its time derivative as the product of two tensors, namely the current velocity gradient $G_{i j} = \frac{\partial u_{i}}{\partial x_{j}}$ and the deformation gradient $F^{n - 1}$ accumulated up to the previous time step, in other words $\frac{\partial F}{\partial t} = G F$ , and $F^{0} = I$ , with $I$ being the identity tensor. While we refer to other studies [MJ83, DT98, BKEO03] for a derivation of this relationship, we can give an intuitive example for the necessity to apply the velocity gradient to the already accumulated deformation, instead of simply integrating the velocity gradient over time. Consider a simple one-dimensional “grain” of length $1.0$ , in which case the deformation tensor only has one component, the compression in $x$ -direction. If one embeds this grain into a convergent flow field for a compressible medium where the dimensionless velocity gradient is $- 0.5$ (e.g. a velocity of zero at its left end at $x = 0.0$ , and a velocity of $- 0.5$ at its right end at $x = 1.0$ ), simply integrating the velocity gradient would suggest that the grain reaches a length of zero after two units of time, and would then “flip” its orientation, which is clearly non-physical. What happens instead can be seen by solving the equation of motion for the right end of the grain $\frac{d x}{d t} = v = - 0.5 x$ . Solving this equation for $x$ leads to $x (t) = e^{- 0.5 t}$ . This is therefore also the solution for $F$ since $F x$ transforms the initial position of $x (t = 0) = 1.0$ into the deformed position of $x (t = 1) = e^{- 0.5}$ , which is the definition of $F$ .

In more general cases a visualization of $F$ is not intuitive, because it contains rotational components that represent a rigid body rotation without deformation. Following [BKEO03] we can polar-decompose the tensor into a positive-definite and symmetric left stretching tensor $L$ , and an orthogonal rotation tensor $Q$ , as $F = L Q$ , therefore $L^{2} = L L^{T} = F F^{T}$ . The left stretching tensor $L$ (or finite strain tensor) then describes the deformation we are interested in, and its eigenvalues $λ_{i}$ and eigenvectors $e_{i}$ describe the length and orientation of the half-axes of the finite strain ellipsoid. Moreover, we will represent the amount of relative stretching at every point by the ratio $ln (λ_{1} ∕ λ_{2})$ , called the natural strain [Rib92].

The full plugin implementing the integration of $F$ can be found in cookbooks/finite_strain/finite_strain.cc and can be compiled with cmake . && make in the cookbooks/finite_strain directory. It can be loaded in a parameter file as an “Additional shared library”, and selected as material model. As it is derived from the “simple” material model, all input parameters for the material properties are read in from the subsection Simple model.

set ??            = ./libfinite_strain.so

subsection ??
  set ?? = finite strain

  subsection ??
    set ?? = 4.7
    set ?? = 3400
    set ?? = 2e-5
    set ?? = 5e21
    set ?? = 7
    set ?? = 1600
  end
end

Figure 31: Accumulated finite strain in an example convection model, as described in Section 5.2.11 at a time of 67.6 Ma. Top panel: Temperature distribution. Bottom panel: Natural strain distribution. Additional black crosses are the scaled eigenvectors of the stretching tensor

L

, showing the direction of stretching and compression.

The plugin was tested against analytical solutions for the deformation gradient tensor in simple and pure shear as described in benchmarks/finite_strain/pure_shear.prm and benchmarks/finite_strain/simple_shear.prm.

We will demonstrate its use at the example of a 2D Cartesian convection model (Figure 31): Heating from the bottom leads to the ascent of plumes from the boundary layer (top panel), and the amount of stretching is visible in the distribution of natural strain (color in lower panel). Additionally, the black crosses show the direction of stretching and compression (the eigenvectors of $L$ ). Material moves to the sides at the top of the plume head, so that it is shortened in vertical direction (short vertical lines) and stretched in horizontal direction (long horizontal lines). The sides of the plume head show the opposite effect. Shear occurs mostly at the edges of the plume head, in the plume tail, and in the bottom boundary layer (black areas in the natural strain distribution).

The example used here shows how history variables can be integrated up over the model evolution. While we do not use these variables actively in the computation (in our example, there is no influence of the accumulated strain on the rheology or any other material property), it would be trivial to extend this material model in a way that material properties depend on the integrated strain: Because the values of the compositional fields are part of what the material model gets as inputs, they can easily be used for computing material model outputs such as the viscosity.

Note: In this model we present the use of multiple compositional fields for other purposes than chemical composition. It would have been feasible to run the same model with particles that track the deformation gradient, as additionally implemented and tested in the simple shear and pure shear benchmarks mentioned in this section. Both approaches have specific advantages, and for scientific computations one needs to evaluate the more suitable strategy. Compositional fields cover the whole domain, but are affected by numerical diffusion, effectively reducing the maximum accumulated strain. Particles only provide finite strain values at discrete positions, but can, if this is desired, be used in fewer numbers and only a part of the model domain (and are much faster in this case). If however there needs to be a large number of particles (possibly because they are used for other purposes as well), then they can be much more expensive. Both approaches can be used to actively influence the rheology in the material model.

5.2.12 Reading in compositional initial composition files generated with geomIO

This section was contributed by Juliane Dannberg

Many geophysical setups require initial conditions with several different materials and complex geometries. Hence, sometimes it would be easier to generate the initial geometries of the materials as a drawing instead of by writing code. The MATLAB-based library geomIO (http://geomio.bitbucket.org/) provides a convenient tool to convert a drawing generated with the vector graphics editor Inkscape (https://inkscape.org/en/) to a data file that can be read into ASPECT. Here, we will demonstrate how this can be done for a 2D setup for a model with one compositional field, but geomIO also has the capability to create 3D volumes based on a series of 2D vector drawings using any number of different materials. Similarly, initial conditions defined in this way can also be used with particles instead of compositional fields.

To obtain the developer version of geomIO, you can clone the bitbucket repository by executing the command

git clone https://bitbucket.org/geomio/geomio.git

or you can download geomIO here. You will then need to add the geomIO source folders to your MATLAB path by running the file located in /path/to/geomio/installation/InstallGeomIO.m. An extensive documentation for how to use geomIO can be found here. Among other things, it explains how to generate drawings in Inkscape that can be read in by geomIO, which involves assigning new attributes to paths in Inkscape’s XML editor. In particular, a new property ‘phase’ has to be added to each path, and set to a value corresponding to the index of the material that should be present in this region in the initial condition of the geodynamic model.

Note: geomIO currently only supports the latest stable version of Inkscape (0.91), and other versions might not work with geomIO or cause errors. Moreover, geomIO currently does not support grouping paths (paths can still be combined using Path

\to

Union, Path

\to

Difference or similar commands), and only the outermost closed contour of a path will be considered. This means that, for example, for modeling a spherical annulus, you would have to draw two circles, and assign the inner one the same phase as the background of your drawing.

We will here use a drawing of a jellyfish located in cookbooks/geomio/jellyfish.svg, where different phases have already been assigned to each path (Figure 32).

Figure 32: Vector drawing of a jellyfish.

Note: The page of your drawing in Inkscape should already have the extents (in px) that you later want to use in your model (in m).

After geomIO is initialized in MATLAB, we run geomIO as described in the documentation, loading the default options and then specifying all the option we want to change, such as the path to the input file, or the resolution:

% set options for geomIO
opt                 = geomIO_Options();
opt.inputFileName   = [’/path/to/aspect/doc/manual/cookbooks/geomio/jellyfish.svg’];
opt.DrawCoordRes    = 21; % optionally change resolution with opt.DrawCoordRes = your value;

% run geomIO
[PathCoord]         = run_geomIO(opt,’2D’);

You can view all of the options available by typing opt in MATLAB.

In the next step we create the grid that is used for the coordinates in the ascii data initial conditions file and assign a phase to each grid point:

% define the bounding box for the output mesh
% (this should be the X extent and Y extent in your Aspect model)
xmin = 0; xmax = opt.svg.width;
ymin = 0; ymax = opt.svg.height;

% set the resolution in the output file:
% [Xp,Yp] = ndgrid(xmin:your_steplength_x:xmax,ymin:your_steplength_y:ymax);
[Xp,Yp] = ndgrid(xmin:15:xmax,ymin:15:ymax);
Phase = zeros(size(Xp));

% assign a phase to each grid point according to your drawing
Phase = assignPhase2Markers(PathCoord, opt, Xp, Yp, Phase);

% plot your output
figure(2)
scatter(Xp(:),Yp(:),10,Phase(:),’filled’);
axis equal
axis([xmin xmax ymin ymax])

You can plot the Phase variable in MATLAB to see if the drawing was read in and all phases are assigned correctly (Figure 33).

Figure 33: Plot of the Phase variable in MATLAB.

Finally, we want to write output in a format that can be read in by ASPECT’s ascii data compositional initial conditions plugin. We write the data into the file jelly.txt:

% the headers ASPECT needs for the ascii data plugin
header1 = ’x’;
header2 = ’y’;
header3 = ’phase’;

% create an array in the correct format for the ascii data plugin
Vx = Xp(:);
Vy = Yp(:);
VPhase = Phase(:);
[m,n] = size(Phase);

% write the data into the output file
fid=fopen(’jelly.txt’,’w’);
fprintf(fid, ’# POINTS: %d %d \n’,[m n]);
fprintf(fid, [’# Columns: ’ header1 ’ ’ header2 ’ ’ header3 ’\n’]);
fprintf(fid, ’%f %f %f \n’, [Vx Vy VPhase]’);
fclose(fid);

To read in the file we just created (a copy is located in ASPECT’s data directory), we set up a model with a box geometry with the same extents we specified for the drawing in px and one compositional field. We choose the ascii data compositional initial conditions and specify that we want to read in our jellyfish. The relevant parts of the input file are listed below:

subsection Geometry model
  set Model name = box

# The extents of the box is the same as the width
# and height of the drawing in px
# (an A4 page = 7350x10500 px).
  subsection Box
    set X extent = 7350
    set Y extent = 10500
  end
end

# We need one compositional field that will be assigned
# the values read in from the ascii data plugin.
subsection Compositional fields
  set Number of fields = 1
end

# We use the ascii data plugin to read in the file created with geomIO.
subsection Initial composition model
  set Model name = ascii data

  subsection Ascii data model
    set Data directory       = $ASPECT_SOURCE_DIR/data/initial-composition/ascii-data/test/
    set Data file name       = jelly.txt
  end
end

# We refine the mesh where compositional gradients are
# high, i.e. at the boundary between the different phases
# assigned to the compositional field through the initial
# condition.
subsection Mesh refinement
  set Refinement fraction                      = 0.99
  set Coarsening fraction                      = 0
  set Initial global refinement                = 5
  set Initial adaptive refinement              = 4
  set Time steps between mesh refinement       = 0
  set Strategy                                 = composition
end

If we look at the output in paraview, we can see our jellyfish, with the mesh refined at the boundaries between the different phases (Figure 34).

Figure 34: ASPECT model output of the jellyfish and corresponding mesh in paraview.

For a geophysical setup, the MATLAB code could be extended to write out the phases into several different columns of the ASCII data file (corresponding to different compositional fields). This initial conditions file could then be used in ASPECT with a material model such as the multicomponent model, assigning each phase different material properties.

An animation of a model using the jellyfish as initial condition and assigning it a higher viscosity can be found here: https://www.youtube.com/watch?v=YzNTubNG83Q.

¹⁹This difference is far smaller than the numerical error in the heat flux on the mesh this data is computed on.

²⁰The statistics file gives this value to more digits: 4.89008498. However, these are clearly more digits than the result is accurate.

²¹For computations of this size, one should test a few time steps in debug mode but then, of course, switch to running the actual computation in optimized mode – see Section 4.3.

²²Note that the statistics file actually contains the outward heat flux for each of the six boundaries, which corresponds to the negative of upward flux for the bottom boundary. The figure therefore shows the negative of the values available in the statistics file.

²³In fact, the pictures are generated using a twice more refined mesh to provide adequate resolution. We keep the default setting of five global refinements in the parameter file as documented above to keep compute time reasonable when using the default settings.

²⁴Of course, this interpretation suggests that we could have achieved the same goal by encoding everything into a single function – that would, for example, have had initial values one, zero and minus one in the three parts of the domain we are interested in.

²⁵We note that this is no different for particles where the position of particles has to be integrated over time and is subject to numerical error. In simulations, their location is therefore not the exact one but also subject to a diffusive process resulting from numerical inaccuracies. Furthermore, in long thin filaments, the number of particles per cell often becomes too small and new particles have to be inserted; their properties are then interpolated from the surrounding particles, a process that also incurs a smoothing penalty.

²⁶The actual values do not matter as much here. They are chosen in such a way that the system – previously driven primarily by the velocity boundary conditions at the top – now also feels the impact of the density variations. To have an effect, the buoyancy induced by the density difference between materials must be strong enough to balance or at least approach the forces exerted by whatever is driving the velocity at the top.

²⁷Particles are enumerated in a way so that first the first processor in a parallel computations numbers all of the particles on its first cell, then its second cell, and so on; then the second processor does the same with particles in the order of the cells it owns; etc. Thus, the “id” shown in the picture is not just a random number, but rather shows the order of cells and how they belonged to the processors that participated in the computation at the time when particles were created. After some time, particles may of course have become well mixed. In any case, this ordering is of no real practical use.

²⁸This is also easy to try experimentally – use the input file from above and select 5 global and 10 adaptive refinement steps, with the refinement criteria set to density, then visualize the solution.

²⁹The outer iterations are only part of the problem. As discussed in [KHB12], each GMRES iteration requires solving a linear system with the elliptic operator $- \nabla \cdot 2 η ε (\cdot)$ . For highly heterogeneous models, such as the one discussed in the current section, this may require a lot of Conjugate Gradient iterations. For example, for 8 global refinement steps, the 30+188 outer iterations without averaging shown in Table 2 require a total of 22,096 inner CG iterations for the elliptic block (and a total of 837 for the approximate Schur complement). Using harmonic averaging, the 30+26 outer iterations require only 1258 iterations on the elliptic block (and 84 on the Schur complement). In other words, the number of inner iterations per outer iteration (taking into account the split into “cheap” and “expensive” outer iterations, see [KHB12]) is reduced from 117 to 47 for the elliptic block and from 3.8 to 1.5 for the Schur complement.

[next] [prev] [prev-tail] [front] [up]