const unsigned int dofs_per_cell = cd.local_dof_indices.size();
- // TODO EXPLAIN:
-
+ // Interface entries are ignored by the boundary_constraints object
+ // above when filling the mg_matrices[cd.level]. Instead, we copy these
+ // entries into the interface matrix of the current level manually:
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
if (mg_constrained_dofs.is_interface_matrix_entry(
in the pressure variable:
<table align="center" class="doxtable">
<tr>
- <td></td>
- <td colspan="3" align="center">Finite element order</td>
+ <th></th>
+ <th colspan="3" align="center">Finite element order</th>
</tr>
<tr>
- <td>Refinement level</td>
- <td>0</td>
- <td>1</td>
- <td>2</td>
+ <th>Refinement level</th>
+ <th>0</th>
+ <th>1</th>
+ <th>2</th>
</tr>
<tr>
- <td>0</td> <td>1.45344</td> <td>0.0831743</td> <td>0.0235186</td>
+ <th>0</th> <td>1.45344</td> <td>0.0831743</td> <td>0.0235186</td>
</tr>
<tr>
- <td>1</td> <td>0.715099</td> <td>0.0245341</td> <td>0.00293983</td>
+ <th>1</th> <td>0.715099</td> <td>0.0245341</td> <td>0.00293983</td>
</tr>
<tr>
- <td>2</td> <td>0.356383</td> <td>0.0063458</td> <td>0.000367478</td>
+ <th>2</th> <td>0.356383</td> <td>0.0063458</td> <td>0.000367478</td>
</tr>
<tr>
- <td>3</td> <td>0.178055</td> <td>0.00159944</td> <td>4.59349e-05</td>
+ <th>3</th> <td>0.178055</td> <td>0.00159944</td> <td>4.59349e-05</td>
</tr>
<tr>
- <td>4</td> <td>0.0890105</td> <td>0.000400669</td> <td>5.74184e-06</td>
+ <th>4</th> <td>0.0890105</td> <td>0.000400669</td> <td>5.74184e-06</td>
</tr>
<tr>
- <td>5</td> <td>0.0445032</td> <td>0.000100218</td> <td>7.17799e-07</td>
+ <th>5</th> <td>0.0445032</td> <td>0.000100218</td> <td>7.17799e-07</td>
</tr>
<tr>
- <td>6</td> <td>0.0222513</td> <td>2.50576e-05</td> <td>9.0164e-08</td>
+ <th>6</th> <td>0.0222513</td> <td>2.50576e-05</td> <td>9.0164e-08</td>
</tr>
<tr>
- <td></td> <td>$O(h)$</td> <td>$O(h^2)$</td> <td>$O(h^3)$</td>
+ <th></th> <th>$O(h)$</th> <th>$O(h^2)$</th> <th>$O(h^3)$</th>
</tr>
</table>
in the velocity variables:
<table align="center" class="doxtable">
<tr>
- <td></td>
- <td colspan="3" align="center">Finite element order</td>
+ <th></th>
+ <th colspan="3" align="center">Finite element order</th>
</tr>
<tr>
- <td>Refinement level</td>
- <td>0</td>
- <td>1</td>
- <td>2</td>
+ <th>Refinement level</th>
+ <th>0</th>
+ <th>1</th>
+ <th>2</th>
</tr>
<tr>
- <td>0</td> <td>0.367423</td> <td>0.127657</td> <td>5.10388e-14</td>
+ <th>0</th> <td>0.367423</td> <td>0.127657</td> <td>5.10388e-14</td>
</tr>
<tr>
- <td>1</td> <td>0.175891</td> <td>0.0319142</td> <td>9.04414e-15</td>
+ <th>1</th> <td>0.175891</td> <td>0.0319142</td> <td>9.04414e-15</td>
</tr>
<tr>
- <td>2</td> <td>0.0869402</td> <td>0.00797856</td> <td>1.23723e-14</td>
+ <th>2</th> <td>0.0869402</td> <td>0.00797856</td> <td>1.23723e-14</td>
</tr>
<tr>
- <td>3</td> <td>0.0433435</td> <td>0.00199464</td> <td>1.86345e-07</td>
+ <th>3</th> <td>0.0433435</td> <td>0.00199464</td> <td>1.86345e-07</td>
</tr>
<tr>
- <td>4</td> <td>0.0216559</td> <td>0.00049866</td> <td>2.72566e-07</td>
+ <th>4</th> <td>0.0216559</td> <td>0.00049866</td> <td>2.72566e-07</td>
</tr>
<tr>
- <td>5</td> <td>0.010826</td> <td>0.000124664</td> <td>3.57141e-07</td>
+ <th>5</th> <td>0.010826</td> <td>0.000124664</td> <td>3.57141e-07</td>
</tr>
<tr>
- <td>6</td> <td>0.00541274</td> <td>3.1166e-05</td> <td>4.46124e-07</td>
+ <th>6</th> <td>0.00541274</td> <td>3.1166e-05</td> <td>4.46124e-07</td>
</tr>
<tr>
- <td></td> <td>$O(h)$</td> <td>$O(h^2)$</td> <td>$O(h^3)$</td>
+ <th></th> <td>$O(h)$</td> <td>$O(h^2)$</td> <td>$O(h^3)$</td>
</tr>
</table>
The result concerning the convergence order is the same here.
@endcode
For 1 through 9 global refinement steps, we then get the following sequence
of point values:
- <table align="center">
- <tr> <td># of refinements</td> <td>$u_h(\frac 13,\frac13)$</td> </tr>
+ <table align="center" class="doxtable">
+ <tr> <th># of refinements</th> <th>$u_h(\frac 13,\frac13)$</th> </tr>
<tr> <td>1</td> <td>0.166667</td> </tr>
<tr> <td>2</td> <td>0.227381</td> </tr>
<tr> <td>3</td> <td>0.237375</td> </tr>
parameters mean, while the first and third should be obvious. Doing the same
study again where we change the number of global refinement steps, we get
the following result:
- <table align="center">
- <tr> <td># of refinements</td> <td>$\int_\Omega u_h(x)\; dx$</td> </tr>
+ <table align="center" class="doxtable">
+ <tr> <th># of refinements</th> <th>$\int_\Omega u_h(x)\; dx$</th> </tr>
<tr> <td>0</td> <td>0.09375000</td> </tr>
<tr> <td>1</td> <td>0.12790179</td> </tr>
<tr> <td>2</td> <td>0.13733440</td> </tr>
@note If you are interested in using the program as the basis for your own
experiments, you will also want to take a look at its continuation in
step-32. Furthermore, step-32 later was developed into the much larger open
-source code Aspect (see http://aspect.dealii.org/ ) that can solve realistic
+source code ASPECT (see https://aspect.geodynamics.org/ ) that can solve realistic
problems and that you may want to investigate before trying to morph step-31
into something that can solve whatever you want to solve.
step, and that we have to modify the above formula slightly. If
$k_n,k_{n-1}$ are the time steps sizes of the current and previous time
step, then we use the approximations
-
-$\frac{\partial T}{\partial t} \approx
+@f{align*}{
+\frac{\partial T}{\partial t} \approx
\frac 1{k_n}
\left(
\frac{2k_n+k_{n-1}}{k_n+k_{n-1}} T^{n}
\frac{k_n+k_{n-1}}{k_{n-1}}T^{n-1}
+
\frac{k_n^2}{k_{n-1}(k_n+k_{n-1})} T^{n-2}
- \right)$
+ \right)
+ @f}
and
-$T^n \approx
+@f{align*}{
+T^n \approx
T^{n-1} + k_n \frac{\partial T}{\partial t}
\approx
T^{n-1} + k_n
\frac{T^{n-1}-T^{n-2}}{k_{n-1}}
=
- \left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2}$,
+ \left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2},
+@f}
and above equation is generalized as follows:
@f{eqnarray*}
\frac{2k_n+k_{n-1}}{k_n+k_{n-1}} T^n
program, in particular in 3d. This is the topic of the step-32
tutorial program which will implement strategies to solve this problem in
%parallel on a cluster. It is also the basis of the much larger open
-source code Aspect (see http://aspect.dealii.org/ ) that can solve realistic
+source code ASPECT (see https://aspect.geodynamics.org/ ) that can solve realistic
problems and that constitutes the further development of step-32.
Another direction would be to make the fluid flow more realistic. The program
</b>
The continuation of development of this program has led to the much larger open
-source code <i>Aspect</i> (see http://aspect.geodynamics.org/) which is much
+source code <i>ASPECT</i> (see http://aspect.geodynamics.org/) which is much
more flexible in solving many kinds of related problems.
</i>
physics would have to be implemented, and several other aspects are
currently missing from this program as well. We will come back to this
issue in the results section again, but state for now that providing a
-realistic description is a goal of the <i>Aspect</i> code in
+realistic description is a goal of the <i>ASPECT</i> code in
development at the time of writing this.
As a reminder, let us again state the equations we want to solve are these:
That said, both step-31 and the current step-32 have not come about by chance
but are certainly meant as wayposts along the path to a more comprehensive
program that will simulate convection in the earth mantle. We call this code
-<i>Aspect</i> (short for <i>Advanced %Solver for Problems in Earth's
+<i>ASPECT</i> (short for <i>Advanced %Solver for Problems in Earth's
ConvecTion</i>); its development is funded by
the <a href="http://www.geodynamics.org">Computational Infrastructure in
Geodynamics</a> initiative with support from the National Science
-Foundation. More information on <i>Aspect</i> is available at
+Foundation. More information on <i>ASPECT</i> is available at
its <a href="https://aspect.geodynamics.org/">homepage</a>.
requiring more CPU hours than is available on a typical
cluster. Consequently, rather than showing a complete simulation here,
let us simply show a couple of pictures we have obtained using the
-successor to this program, called <i>Aspect</i> (short for <i>Advanced
+successor to this program, called <i>ASPECT</i> (short for <i>Advanced
%Solver for Problems in Earth's ConvecTion</i>), that is being
developed independently of deal.II and that already incorporates some
of the extensions discussed below. The following two pictures show
There are many directions in which this program could be extended. As
mentioned at the end of the introduction, most of these are under active
-development in the <i>Aspect</i> (short for <i>Advanced %Solver for Problems
+development in the <i>ASPECT</i> (short for <i>Advanced %Solver for Problems
in Earth's ConvecTion</i>) code at the time this tutorial program is being
finished. Specifically, the following are certainly topics that one should
address to make the program more useful:
There are many other ways to extend the current program. However, rather than
discussing them here, let us point to the much larger open
-source code Aspect (see http://aspect.dealii.org/ ) that constitutes the
+source code ASPECT (see https://aspect.geodynamics.org/ ) that constitutes the
further development of step-32 and that already includes many such possible
extensions.
On a quarter-circle defined by $\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}$ we are
going to solve the Stokes problem
@f{eqnarray*}
- -\Delta \; \textbf{u} + \nabla p &=& (\exp(-100*\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\
+ -\Delta \; \textbf{u} + \nabla p &=& (\exp(-100\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\
-\textrm{div}\; \textbf{u}&=&0,\\
\textbf{u}|_{\Gamma_1}&=&{\bf 0},
@f}
To illustrate how one describes geometries using charts in deal.II, we will
consider a case that originates in an application of the <a
-href="http://aspect.dealii.org">ASPECT mantle convection code</a>, using a
+href="https://aspect.geodynamics.org">ASPECT mantle convection code</a>, using a
data set provided by D. Sarah Stamps. In the concrete application, we were
interested in describing flow in the Earth mantle under the <a
href="http://en.wikipedia.org/wiki/East_African_rift">East African Rift</a>, a
step-22. Details of this approach can be found under the "Block Schur
complement preconditioner" subsection of the "Possible Extensions"
section of step-22. For the preconditioner of the velocity block, we
-borrow a class from <a href="http://aspect.dealii.org">ASPECT</a>
+borrow a class from <a href="https://aspect.geodynamics.org">ASPECT</a>
called @p BlockSchurPreconditioner that has the option to solve for
the inverse of $A$ or just apply one preconditioner sweep for it
instead, which provides us with an expensive and cheap approach,
// of a single preconditioner application.
// 3. do not use InverseMatrix but explicitly call SolverCG.
// This approach is also used in the ASPECT code
- // (see http://aspect.dealii.org) that solves the Stokes equations in
+ // (see https://aspect.geodynamics.org) that solves the Stokes equations in
// the context of simulating convection in the earth mantle, and which
// has been used to solve problems on many thousands of processors.
//
<table align="center" class="doxtable">
<tr>
<th>degree</th>
- <td>1</td>
- <td>2</td>
- <td>3</td>
- <td>4</td>
- <td>5</td>
- <td>6</td>
- <td>7</td>
- <td>8</td>
- <td>9</td>
- <td>10</td>
- <td>11</td>
- <td>12</td>
+ <th>1</th>
+ <th>2</th>
+ <th>3</th>
+ <th>4</th>
+ <th>5</th>
+ <th>6</th>
+ <th>7</th>
+ <th>8</th>
+ <th>9</th>
+ <th>10</th>
+ <th>11</th>
+ <th>12</th>
</tr>
<tr>
<th>Number of DoFs</th>
<table align="center" class="doxtable">
<tr>
<th>degree</th>
- <td>1</td>
- <td>2</td>
- <td>3</td>
- <td>4</td>
- <td>5</td>
- <td>6</td>
- <td>7</td>
- <td>8</td>
- <td>9</td>
- <td>10</td>
- <td>11</td>
- <td>12</td>
+ <th>1</th>
+ <th>2</th>
+ <th>3</th>
+ <th>4</th>
+ <th>5</th>
+ <th>6</th>
+ <th>7</th>
+ <th>8</th>
+ <th>9</th>
+ <th>10</th>
+ <th>11</th>
+ <th>12</th>
</tr>
<tr>
<th>Number of DoFs</th>
projects / dealii.git / commitdiff