From: bangerth Date: Thu, 21 Aug 2008 15:15:15 +0000 (+0000) Subject: Discuss correct choice of beta and time step. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a3d55800dc54a80a6dd5296359c6e3d49f113a95;p=dealii-svn.git Discuss correct choice of beta and time step. git-svn-id: https://svn.dealii.org/trunk@16624 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/doc/results.dox b/deal.II/examples/step-31/doc/results.dox index 5210a8fdb7..6083f95148 100644 --- a/deal.II/examples/step-31/doc/results.dox +++ b/deal.II/examples/step-31/doc/results.dox @@ -35,7 +35,7 @@ temperature. Below we'll discuss a few numerical experiments to choose constants. -
Choosing ck and β
+

Choosing ck and β

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. + + +
Results for Q1 elements
+ +Here is what we get for $\beta=0.01, \beta=0.1$, and $\beta=0.5$, different choices of $c_k$, and bilinear elements (temperature_degree=1) 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. + + +
Results for Q2 elements
+ +One can repeat the same sequence of experiments for higher order +elements as well. Here are the graphs for bi-quadratic shape functions +(temperature_degree=2) for the temperature, while we +retain the $Q_2/Q_1$ stable Taylor-Hood element for the Stokes system: + + + + + + + + + + +
+ @image html "step-33.timestep.q2.beta=0.01.png" "" width=4cm + + @image html "step-33.timestep.q2.beta=0.03.png" "" width=4cm +
+ @image html "step-33.timestep.q2.beta=0.1.png" "" width=4cm +
+ +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 grid point away (which for $Q_1$ elements +is $h_K$, but for $Q_2$ elements is $h_K/2$). + + +
Conclusions
+ +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 temperature_degree). + +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. +

Possible extensions

diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 075a0574b6..9582f59b82 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -476,8 +476,6 @@ class BoussinesqFlowProblem Triangulation triangulation; const unsigned int stokes_degree; - const unsigned int temperature_degree; - FESystem stokes_fe; DoFHandler stokes_dof_handler; ConstraintMatrix stokes_constraints; @@ -491,6 +489,7 @@ class BoussinesqFlowProblem BlockVector stokes_rhs; + const unsigned int temperature_degree; FE_Q temperature_fe; DoFHandler temperature_dof_handler; ConstraintMatrix temperature_constraints; @@ -780,13 +779,16 @@ template BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) : triangulation (Triangulation::maximum_smoothing), - stokes_degree (degree), - temperature_degree (degree), + + stokes_degree (1), stokes_fe (FE_Q(stokes_degree+1), dim, FE_Q(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::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;