]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Discuss correct choice of beta and time step.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 15:15:15 +0000 (15:15 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 15:15:15 +0000 (15:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@16624 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-31/doc/results.dox
deal.II/examples/step-31/step-31.cc

index 5210a8fdb7291c871d0dc3bca87555cdd7d725ee..6083f95148fa73eda6de2d74acb70ad0b2a21030 100644 (file)
@@ -35,7 +35,7 @@ temperature. Below we'll discuss a few numerical experiments to choose
 constants. 
 
 
-<h5> Choosing <i>c<sub>k</sub></i> and &beta; </h5>
+<h4> Choosing <i>c<sub>k</sub></i> and &beta; </h4>
 
 These two constants are definitely linked in some way. The reason is easy to
 see: In the case of a pure advection problem, 
@@ -80,7 +80,12 @@ in the domain. What we expect to see is this: If we choose the time step too
 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:
 
@@ -110,6 +115,8 @@ variations, but if we choose
 $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)}}$ 
@@ -117,7 +124,7 @@ for $\beta=0.5$ (this time because much diffusion requires a small time
 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
@@ -145,6 +152,74 @@ hand, the larger the artificial diffusion, the more the hot spot is
 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)}}$ &mdash; 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>
 
index 075a0574b60d147815cf1c8fcbceb2c3e64b56f2..9582f59b82b5a53fbaf5c26a4a7b687083ead405 100644 (file)
@@ -476,8 +476,6 @@ class BoussinesqFlowProblem
     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;
@@ -491,6 +489,7 @@ class BoussinesqFlowProblem
     BlockVector<double>       stokes_rhs;
 
     
+    const unsigned int        temperature_degree;    
     FE_Q<dim>                 temperature_fe;
     DoFHandler<dim>           temperature_dof_handler;
     ConstraintMatrix          temperature_constraints;
@@ -780,13 +779,16 @@ template <int dim>
 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),
@@ -1607,7 +1609,7 @@ double compute_viscosity(
   const double                        old_time_step
 )
 {
-  const double beta = 0.1;
+  const double beta = 0.03;
   const double alpha = 1;
   
   if (global_u_infty == 0)
@@ -2017,13 +2019,11 @@ void BoussinesqFlowProblem<dim>::solve ()
     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;
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.