From 64546e32bdaa0b3c11934d50a5e38bcbe93e4104 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 4 Jun 2013 06:11:19 +0000 Subject: [PATCH] Increase size of stabilization parameter to get stable solutions at large times. git-svn-id: https://svn.dealii.org/trunk@29745 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 +++++++ deal.II/examples/step-31/doc/results.dox | 18 +++++++++++------- deal.II/examples/step-31/step-31.cc | 10 +++++----- 3 files changed, 23 insertions(+), 12 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index f8e3cf0f7d..6aceb090d3 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -137,6 +137,13 @@ this function.

Specific improvements

    +
  1. Fixed: The stabilization parameter for the artificial diffusion +in the step-31 tutorial program has been increased slightly to avoid +instabilities at later times (t > 60). +
    +(Martin Kronbichler, 2013/06/04) +
  2. +
  3. Fixed: If an exception was generated on a task created by Threads::new_task, the program would terminate with a segmentation fault, leaving little trace of what had happened. This is now handled diff --git a/deal.II/examples/step-31/doc/results.dox b/deal.II/examples/step-31/doc/results.dox index 9310d49090..e933bb368f 100644 --- a/deal.II/examples/step-31/doc/results.dox +++ b/deal.II/examples/step-31/doc/results.dox @@ -506,13 +506,16 @@ 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 +for the $Q_1$ element, and then 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$). +is $h_K$, but for $Q_2$ elements is $h_K/2$). It turns out that $\beta$ +needs to be sligthly larger for obtaining stable results also late in +the simulation at times larger than 60, so we actually choose it as +$\beta = 0.034$ in the code.
    Results for 3d
    @@ -520,7 +523,7 @@ is $h_K$, but for $Q_2$ elements is $h_K/2$). One can repeat these experiments in 3d and find the optimal time step for each value of $\beta$ and find the best value of $\beta$. What one finds is that for the same $\beta$ already used in 2d, the time steps -needs to be a bit small, by around a factor of 1.2 or so. This is +needs to be a bit smaller, by around a factor of 1.2 or so. This is easily explained: the time step restriction is $k=\min_K \frac{ch_K}{\|\mathbf{u}\|_{L^\infty(K)}}$ where $h_K$ is the diameter of the cell. However, what is really needed is the @@ -537,12 +540,13 @@ can take.
    Conclusions
    -Concluding, from the simple computations above, $\beta=0.03$ appears to be a +Concluding, from the simple computations above, $\beta=0.034$ appears to be a good choice for the stabilization parameter in 2d, and $\beta=0.05$ in 3d. In -a dimension independent way, we can model this as $\beta=0.015d$. If one does +a dimension independent way, we can model this as $\beta=0.017d$. If one does longer computations (several thousand time steps) on finer meshes, one -realizes that that's not quite small enough and that for stability one will -have to reduce the above values a bit more (by about a factor of $\frac 78$). +realizes that the time step size is not quite small enough and that for +stability one will have to reduce the above values a bit more (by about a +factor of $\frac 78$). As a consequence, a formula that reconciles 2d, 3d, and variable polynomial degree and takes all factors in account reads as follows: diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 7d9a55ac91..a61da97508 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -3,7 +3,7 @@ /* $Id$ */ /* */ -/* Copyright (C) 2007-2012 by the deal.II authors */ +/* Copyright (C) 2007-2013 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -750,7 +750,7 @@ namespace Step31 // introduction. // // There are some universal constants worth mentioning here. First, we need - // to fix $\beta$; we choose $\beta=0.015\cdot dim$, a choice discussed in + // to fix $\beta$; we choose $\beta=0.017\cdot dim$, a choice discussed in // detail in the results section of this tutorial program. The second is the // exponent $\alpha$; $\alpha=1$ appears to work fine for the current // program, even though some additional benefit might be expected from @@ -779,7 +779,7 @@ namespace Step31 const double global_T_variation, const double cell_diameter) const { - const double beta = 0.015 * dim; + const double beta = 0.017 * dim; const double alpha = 1; if (global_u_infty == 0) @@ -1828,12 +1828,12 @@ namespace Step31 const double maximal_velocity = get_maximal_velocity(); if (maximal_velocity >= 0.01) - time_step = 1./(1.6*dim*std::sqrt(1.*dim)) / + time_step = 1./(1.7*dim*std::sqrt(1.*dim)) / temperature_degree * GridTools::minimal_cell_diameter(triangulation) / maximal_velocity; else - time_step = 1./(1.6*dim*std::sqrt(1.*dim)) / + time_step = 1./(1.7*dim*std::sqrt(1.*dim)) / temperature_degree * GridTools::minimal_cell_diameter(triangulation) / .01; -- 2.39.5