constants.
-<h5> Choosing <i>c<sub>k</sub></i> and β </h5>
+<h4> Choosing <i>c<sub>k</sub></i> and β </h4>
These two constants are definitely linked in some way. The reason is easy to
see: In the case of a pure advection problem,
big (i.e. choose a $c_k$ bigger than theoretically allowed) then we will get
exponential growth of the temperature. If we choose $\beta$ too small, then
the transport stabilization becomes insufficient and the solution will show
-significant oscillations but not exponential growth. Here is what we get for
+significant oscillations but not exponential growth.
+
+
+<h5>Results for Q<sub>1</sub> elements</h5>
+
+Here is what we get for
$\beta=0.01, \beta=0.1$, and $\beta=0.5$, different choices of $c_k$, and
bilinear elements (<code>temperature_degree=1</code>) in 2d:
$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$
or smaller, then the scheme is
stable though a bit wobbly. For more artificial diffusion, we can choose
+$k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$
+or smaller for $\beta=0.03$,
$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$
or smaller for $\beta=0.1$, and again need
$k=\frac 1{15}\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$
step).
So how to choose? If we were simply interested in a large time step, then we
-would go with $\beta=...$ and
+would go with $\beta=0.1$ and
$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$.
On the other hand, we're also interested in accuracy and here it may be of
interest to actually investigate what these curves show. To this end note that
diffused. Note that for this criterion, the time step size does not play a
significant role.
+So to sum up, likely the best choice would appear to be $\beta=0.03$
+and $k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$. The curve is
+a bit wobbly, but overall pictures looks pretty reasonable with the
+exception of some over and undershoots close to the start time due to
+Gibb's phenomenon.
+
+
+<h5>Results for Q<sub>2</sub> elements</h5>
+
+One can repeat the same sequence of experiments for higher order
+elements as well. Here are the graphs for bi-quadratic shape functions
+(<code>temperature_degree=2</code>) for the temperature, while we
+retain the $Q_2/Q_1$ stable Taylor-Hood element for the Stokes system:
+
+<table align="center" border="1" cellspacing="3" cellpadding="3">
+ <tr>
+ <td>
+ @image html "step-33.timestep.q2.beta=0.01.png" "" width=4cm
+ </td>
+ <td>
+ @image html "step-33.timestep.q2.beta=0.03.png" "" width=4cm
+ </td>
+ </tr>
+
+ <tr>
+ <td>
+ @image html "step-33.timestep.q2.beta=0.1.png" "" width=4cm
+ </td>
+ </tr>
+</table>
+
+Again, small values of $\beta$ lead to less diffusion but we have to
+choose the time step very small to keep things under control. Too
+large values of $\beta$ make for more diffusion, but again require
+small time steps. The best value would appear to be $\beta=0.03$, as
+for the $Q_1$ element, and the we have to choose
+$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ — exactly
+half the size for the $Q_1$ element, a fact that may not be surprising
+if we state the CFL condition as the requirement that the time step be
+small enough so that the distance transport advects in each time step
+is no longer than one <i>grid point</i> away (which for $Q_1$ elements
+is $h_K$, but for $Q_2$ elements is $h_K/2$).
+
+
+<h5>Conclusions</h5>
+
+Concluding, $\beta=0.03$ appears to be a good choice for the
+stabilization parameter, and then
+$k=\frac 14 \frac 1{q_T}\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}$
+is an appropriate time step, where $q_T$ is the polynomial degree of
+the temperature shape functions (in the program, this corresponds to
+the variable <code>temperature_degree</code>).
+
+As for the question of whether to use $Q_1$ or $Q_2$ elements for the
+temperature, the following considerations may be useful: First,
+solving the temperature equation is hardly a factor in the overall
+scheme since almost the entire compute time goes into solving the
+Stokes system in each time step. Higher order elements for the
+temperature equation are therefore not a significant drawback. On the
+other hand, if one compares the size of the over- and undershoots the
+solution produces due to the discontinuous source description, one
+notices that for the choice of $\beta$ and $k$ as above, the $Q_1$
+solution dips down to around $-0.47$, whereas the $Q_2$ solution only
+goes to $-0.13$ (remember that the exact solution should never become
+negative at all. This means that the $Q_2$ solution is significantly
+more accurate; the program therefore uses these higher order elements,
+despite the penalty we pay in terms of smaller time steps.
+
<h3> Possible extensions </h3>
Triangulation<dim> triangulation;
const unsigned int stokes_degree;
- const unsigned int temperature_degree;
-
FESystem<dim> stokes_fe;
DoFHandler<dim> stokes_dof_handler;
ConstraintMatrix stokes_constraints;
BlockVector<double> stokes_rhs;
+ const unsigned int temperature_degree;
FE_Q<dim> temperature_fe;
DoFHandler<dim> temperature_dof_handler;
ConstraintMatrix temperature_constraints;
BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (const unsigned int degree)
:
triangulation (Triangulation<dim>::maximum_smoothing),
- stokes_degree (degree),
- temperature_degree (degree),
+
+ stokes_degree (1),
stokes_fe (FE_Q<dim>(stokes_degree+1), dim,
FE_Q<dim>(stokes_degree), 1),
stokes_dof_handler (triangulation),
+
+ temperature_degree (2),
temperature_fe (temperature_degree),
temperature_dof_handler (triangulation),
+
time_step (0),
old_time_step (0),
timestep_number (0),
const double old_time_step
)
{
- const double beta = 0.1;
+ const double beta = 0.03;
const double alpha = 1;
if (global_u_infty == 0)
stokes_constraints.distribute (stokes_solution);
}
- // TODO: determine limit of stability for
- // the time step (whether it needs to be /4
- // or whether we could get away with a
- // bigger time step)
old_time_step = time_step;
- time_step = GridTools::minimal_cell_diameter(triangulation) /
- std::max (get_maximal_velocity(), .01) / 4;
+ time_step = 1./4. /
+ temperature_degree *
+ GridTools::minimal_cell_diameter(triangulation) /
+ std::max (get_maximal_velocity(), .01);
temperature_solution = old_temperature_solution;