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

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index f8e3cf0f7d728483159b1adaf9a95ab276195071..6aceb090d35a0b828db3dcc3cea395a890b454f4 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -137,6 +137,13 @@ this function.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: The stabilization parameter for the artificial diffusion
+in the step-31 tutorial program has been increased slightly to avoid
+instabilities at later times (<i>t</i> > 60).
+<br>
+(Martin Kronbichler, 2013/06/04)
+</li>
+
 <li> 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 9310d4909052b3e76fd8ce172c9c0f803f1e317a..e933bb368f7984bacf8d5adc5beba3ee36d86919 100644 (file)
--- 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 <i>grid point</i> 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.
 
 
 <h5>Results for 3d</h5>
@@ -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 <i>diameter</i> of the cell. However, what is really needed is the
@@ -537,12 +540,13 @@ can take.
 
 <h5>Conclusions</h5>
 
-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 7d9a55ac919fbe2f6aa0051c629e6b0c7e27c943..a61da9750866364897ba5440d5218b231e88cf05 100644 (file)
--- 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;