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 f48ee54347ed1855d9124f6002d802c6eba03323..e33facb240e632bef16a461fbeef1f6e1b5be32f 100644 (file)
--- 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);
   {