(sec:cookbooks:convection_box_3d)=
# 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, 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](https://github.com/geodynamics/aspect/blob/main/cookbooks/convection_box_3d/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 {ref}`parameters:Geometry_20model`) 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:
``` {literalinclude} start.part.prm
```
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
``` {literalinclude} initial.part.prm
```
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:
``` {literalinclude} amr.part.prm
```
Finally, as we have seen in the previous section, a computation with $Ra=10^4$
does not lead to a simulation that is exactly exciting. Let us choose
$Ra=10^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 $\alpha=10^{-10}$ and setting
``` {literalinclude} gravity.part.prm
```
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:
``` {literalinclude} postprocess.part.prm
```
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:
``` {literalinclude} checkpoint.part.prm
```
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 {ref}`parameters:Global`).
Running with this input file requires a bit of patience[^footnote1] since the number of
degrees of freedom is just so large: it starts with a bit over 330,000…
``` ksh
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 {numref}`fig:box-3d-solution` and you can find a movie of
the motion that ensues from the heating at the bottom at
. 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
. {numref}`fig:box-3d-mesh` 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
{numref}`fig:convection-box-fields`. In contrast to that picture, the
simulation is also clearly non-stationary.
```{figure-md} fig:box-3d-solution
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) and t=0.01, 0.013, 0.018 (bottom row).
```
```{figure-md} fig:box-3d-mesh
Convection in a 3d box: Meshes and large-scale velocity field for the third, fourth and sixth of the snapshots shown in Fig. 6.
```
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, {numref}`fig:box-3d-heat-flux`
only shows the upward heat flux through the bottom and top boundaries of the
domain as a function of time.[^footnote2] 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 {numref}`fig:box-3d-solution`. 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-md} fig:box-3d-heat-flux
Convection in a 3d box: Upward heat flux through the bottom and top boundaries as a function of time.
```
[^footnote1]: 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 {ref}`sec:run-aspect:debug-mode`.
[^footnote2]: 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.