From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Mon, 17 Aug 2009 20:39:30 +0000 (+0000)
Subject: Better document what we were doing here.
X-Git-Tag: v8.0.0~7276
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2b487132c5a21921e838e835595ed5c3cf3b4cc3;p=dealii.git

Better document what we were doing here.


git-svn-id: https://svn.dealii.org/trunk@19292 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc
index f48ee54347..e33facb240 100644
--- a/deal.II/examples/step-31/step-31.cc
+++ b/deal.II/examples/step-31/step-31.cc
@@ -2293,6 +2293,47 @@ void BoussinesqFlowProblem<dim>::solve ()
 				   // time step is discussed in the results
 				   // section of this program.
 				   //
+				   // There is a snatch here. The formula
+				   // contains a division by the maximum value
+				   // of the velocity. However, at the start
+				   // of the computation, we have a constant
+				   // temperature field (we start with a
+				   // constant temperature, and it will be
+				   // non-constant only after the first time
+				   // step during which the source
+				   // acts). Constant temperature means that
+				   // no buoyancy acts, and so the velocity is
+				   // zero. Dividing by it will not likely
+				   // lead to anything good.
+				   //
+				   // To avoid the resulting infinite time
+				   // step, we ask whether the maximal
+				   // velocity is very small (in particular
+				   // smaller than the values we encounter
+				   // during any of the following time steps)
+				   // and if so rather than dividing by zero
+				   // we just divide by a small value,
+				   // resulting in a large but finite time
+				   // step.
+  old_time_step = time_step;
+  const double maximal_velocity = get_maximal_velocity();
+
+  if (maximal_velocity >= 0.01)
+    time_step = 1./(1.6*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)) /
+		temperature_degree *
+		GridTools::minimal_cell_diameter(triangulation) /
+		.01;
+
+  std::cout << "   " << "Time step: " << time_step
+	    << std::endl;
+  
+  temperature_solution = old_temperature_solution;
+
 				   // Next we set up the temperature system
 				   // and the right hand side using the
 				   // function
@@ -2326,19 +2367,6 @@ void BoussinesqFlowProblem<dim>::solve ()
 				   // Finally, we solve, distribute the
 				   // hanging node constraints and write out
 				   // the number of iterations.
-  old_time_step = time_step;
-  const double maximal_velocity = get_maximal_velocity();
-  time_step = 1./(1.6*dim*std::sqrt(1.*dim)) /
-	      temperature_degree *
-	      GridTools::minimal_cell_diameter(triangulation) /
-              std::max (maximal_velocity, .01);
-
-  std::cout << "   " << "Time step: " << time_step
-	    << std::endl;
-  
-  temperature_solution = old_temperature_solution;
-
-
   assemble_temperature_system (maximal_velocity);
   {