below, things here will not be that simply scaled: $\eta$ will be around
$10^{21}$, velocities on the order of $10^{-8}$, pressure around $10^8$, and
the diameter of the domain is $10^7$. In other words, the order of magnitude
-for the first equation is going to be
+for the first equation is going to be
$\eta\text{div}\varepsilon(\mathbf u) \approx 10^{21} \frac{10^{-8}}{(10^7)^2}
\approx 10^{-1}$, whereas the second equation will be around
$\text{div}{\mathbf u}\approx \frac{10^{-8}}{10^7} \approx 10^{-15}$. Well, so
free. As an auxiliary problem, it turns out that it is difficult to find a
tolerance that always works; in practice, one often ends up with a tolerance
that requires 30 or 40 iterations for most time steps, and 10,000 for some
-others.
+others.
So what's a numerical analyst to do in a case like this? The answer is to
start at the root and first make sure that everything is mathematically
operator at the top right). This, however, can be cured by introducing a
scaled pressure $\hat p = \frac{L}{\eta}p$, and we get the scaled equations
@f{eqnarray*}
- -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) +
+ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) +
\nabla \left(\frac{\eta}{L} \hat p\right) &=&
-\rho \; \beta \; T \mathbf{g},
\\
Schur complement $S=B A^{-1} B^T$ by approximating
$-\text{div}(-\eta\Delta)^{-1}\nabla \approx \frac 1{\eta} \mathbf{1}$. Now,
however, we have re-scaled the $B$ and $B^T$ operators. So $S$ should now
-approximate
+approximate
$-\frac{\eta}{L}\text{div}(-\eta\Delta)^{-1}\nabla \frac{\eta}{L} \approx
\left(\frac{\eta}{L}\right)^2 \frac 1{\eta} \mathbf{1}$.
This is exactly the operator we use to approximate $S$.
As a reminder, let us again state the equations we want to solve are these:
@f{eqnarray*}
- -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) +
+ -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) +
\nabla \left( \frac{\eta}{L} \hat p\right) &=&
-\rho \; \beta \; T \mathbf{g},
\\
All of these pieces of equation data are defined in the program in the
<code>EquationData</code> namespace. When run, the program produces
-maximal velocities around 10 centimeters per year, approximately the
-physically correct order of magnitude. We will set the end time to 30
-million years.
+long-term maximal velocities around 10-40 centimeters per year (see
+the results section below), approximately the physically correct order
+of magnitude. We will set the end time to 1 billion years.
<h1>Results</h1>
When run, the program simulates convection in 3d in much the same way
-as @ref step_31 "step-31" did, though with an entirely different testcase.
+as @ref step_31 "step-31" did, though with an entirely different testcase.
<h3>Comparison of results with step-31</h3>
If we run the program as shown above, the output will look roughly like this,
producing the final part of the output after some 2 days when run on 10
-processors:
+processors:
<code>
<pre>
</code>
As can be seen here, we spend most of the compute time in assembling linear
-systems, refining the mesh, and in particular in solving the Stokes and
-temperature linear systems.
+systems, refining the mesh, and in particular in solving the Stokes and
+temperature linear systems.
The 50% spent on solving the linear systems are affected in large part
because the Brazos cluster has a relatively slow ethernet interconnect. A
with estimating the error by making sure that each processor only estimates
the error on those cells it owns.
+
+The program writes output every 25th time step, but we won't show all
+2100 or so images this produces. Rather, let us only show the output
+from every 2500th time step here, even though this does, of course,
+not do full justice to the dynamics that are going on:
+<table>
+ <tr>
+ <td>
+ @image html step-32.2d.temperature.0000.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.0100.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.0200.png
+ <td>
+ </tr>
+
+ <tr>
+ <td>
+ @image html step-32.2d.temperature.0300.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.0400.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.0500.png
+ <td>
+ </tr>
+
+ <tr>
+ <td>
+ @image html step-32.2d.temperature.0600.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.0700.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.0800.png
+ <td>
+ </tr>
+
+ <tr>
+ <td>
+ @image html step-32.2d.temperature.0900.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.1000.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.1100.png
+ <td>
+ </tr>
+
+ <tr>
+ <td>
+ @image html step-32.2d.temperature.1200.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.1300.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.1400.png
+ <td>
+ </tr>
+
+ <tr>
+ <td>
+ @image html step-32.2d.temperature.1500.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.1600.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.1700.png
+ <td>
+ </tr>
+
+ <tr>
+ <td>
+ @image html step-32.2d.temperature.1800.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.1900.png
+ <td>
+
+ <td>
+ @image html step-32.2d.temperature.2000.png
+ <td>
+ </tr>
+
+
+ <tr>
+ <td>
+ @image html step-32.2d.temperature.2100.png
+ <td>
+
+ <td>
+ @image html step-32.2d.grid.2100.png
+ <td>
+
+ <td>
+ @image html step-32.2d.partition.2100.png
+ <td>
+ </tr>
+</table>
+
+The last two images show the grid as well as the partitioning of the
+mesh for the last timestep shown into the 10 subdomains used for this
+computation. The full dynamics are really only visible by looking at
+an animation. <a
+href="http:/www.math.tamu.edu/~bangerth/images/pictures/convection-outward/step-32.2d.convection.gif">At
+this site</a> is such an animation. Beware that this animation is
+about 20MB large, though it is well worth watching due to its almost
+artistic quality.
+
+If you watch the movie, you'll see that the convection pattern goes
+through several stages: First, it gets rid of the instable temperature
+layering with the hot material overlaid by the dense cold
+material. After this great driver is removed and we have a sort of
+stable situation, a few blobs start to separate from the hot boundary
+layer at the inner ring and rise up, with a few cold fingers also
+dropping down from the outer ring. During this phase, the solution
+remains mostly symmetric, reflecting the 12-fold symmetry of the
+original mesh. In a final phase, the fluid enters vigorous chaotic
+stirring in which all symmetries are lost. This is a pattern that then
+continues to dominate flow.
+
+These different phases can also be identified if we look at the
+maximal velocity as a function of time in the simulation:
+
+@image html step-32.2d.t_vs_vmax.png
+
+Here, the velocity (shown in centimeters per year) becomes very large,
+to the order of several meters per year) at the beginning when the
+temperature layering is instable. It then calms down to relatively
+small values before picking up again in the chaotic stirring
+regime. There, it remains in the range of 10-40 centimeters per year,
+quite within the physically expected region.