]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Finish documenting step-26.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 18 Dec 2013 20:18:21 +0000 (20:18 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 18 Dec 2013 20:18:21 +0000 (20:18 +0000)
git-svn-id: https://svn.dealii.org/trunk@32053 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-26/doc/builds-on
deal.II/examples/step-26/doc/kind
deal.II/examples/step-26/doc/results.dox
deal.II/examples/step-26/doc/tooltip
deal.II/examples/step-26/step-26.cc

index 48a0f738761717ac5d7b461e0b04951a515849bc..17402734c787ca2681ca4585d5db3399ad0e9c49 100644 (file)
@@ -1 +1 @@
-step-4
+step-6
index c1d9154931a75b6d96f1ece6725575e19c39d429..86a44aa1efa4783ce24f5c46afffc3c91ce569cc 100644 (file)
@@ -1 +1 @@
-techniques
+time dependent
index 6a401f850da88ef692302e36d44dc482d5a21e7b..798cc8f298e2906a1c15440177df8ecbac71b106 100644 (file)
@@ -38,7 +38,7 @@ Number of active cells: 1803
 Number of degrees of freedom: 2109
 @endcode
 
-Maybe of more interest is a visualization of the solution and the mesh on wich
+Maybe of more interest is a visualization of the solution and the mesh on which
 it was computed:
 
 <img
@@ -87,7 +87,22 @@ Ralf Hartmann, published by the University of Heidelberg, Germany, in 2002.
 
 If you look at the meshes in the movie above, it is clear that they are not
 particularly well suited to the task at hand. In fact, they look rather
-random. The underlying reason is that we only adapt the mesh once every fifth
+random. 
+
+There are two factors at play. First, there are some islands where cells
+have been refined but that are surrounded by non-refined cells (and there
+are probably also a few occasional coarsened islands). These are not terrible,
+as they most of the time do not affect the approximation quality of the mesh,
+but they also don't help because so many of their additional degrees of
+freedom are in fact constrained by hanging node constraints. That said,
+this is easy to fix: the Triangulation class takes an argument to its
+constructor indicating a level of "mesh smoothing". Passing one of many 
+possible flags, this instructs the triangulation to refine some additional 
+cells, or not to refine some cells, so that the resulting mesh does not have 
+these artifacts.
+
+The second problem is more severe: the mesh appears to lag the solution. 
+The underlying reason is that we only adapt the mesh once every fifth
 time step, and only allow for a single refinement in these cases. Whenever a
 source switches on, the solution had been very smooth in this area before and
 the mesh was consequently rather coarse. This implies that the next time step
index 7b3e8b447aac35e188c2a23320b55e56810dfdae..9bd689e03e8b695ecd410d78a5bc683cde3768a0 100644 (file)
@@ -1 +1 @@
-Complicated boundary descriptions.
+The heat equation. Time dependent meshes.
index 28bcb029fa1453ddebe807c79f5489b62c1e51c3..72adfa1b858e35dbed5f2cb7c0a69460910e2c23 100644 (file)
@@ -194,7 +194,12 @@ namespace Step26
 
   // @sect3{The <code>HeatEquation</code> implementation}
   //
-  // The next step then is the implementation of the main class.
+  // It is time now for the implementation of the main class. Let's
+  // start with the constructor which selects a linear element, a time
+  // step constant at 1/500 (remember that one period of the source
+  // on the right hand side was set to 0.2 above, so we resolve each
+  // period with 100 time steps) and chooses the Crank Nicolson method
+  // by setting $\theta=1/2$.
   template<int dim>
   HeatEquation<dim>::HeatEquation ()
     :
@@ -202,11 +207,30 @@ namespace Step26
     dof_handler(triangulation),
     time_step(1. / 500),
     theta(0.5)
-  {
-  }
+  {}
 
 
 
+  // @sect4{<code>HeatEquation::setup_system</code>}
+  //
+  // The next function is the one that sets up the DoFHandler object,
+  // computes the constraints, and sets the linear algebra objects
+  // to their correct sizes. We also compute the mass and Laplace
+  // matrix here by simply calling two functions in the library.
+  //
+  // Note that we compute these matrices taking into account already the
+  // constraints due to hanging nodes. These are all homogenous, i.e.,
+  // they only consist of constraints of the form $U_i = \alpha_{ij} U_j
+  // + \alpha_{ik} U_k$ (whereas inhomogenous constraints would also
+  // have a term not proportional to $U$, i.e., $U_i = \alpha_{ij} U_j
+  // + \alpha_{ik} U_k + c_i$). For this kind of constraint, we can
+  // eliminate hanging nodes independently in the matrix and the
+  // right hand side vectors, but this is not the case for inhomogenous
+  // constraints for which we can eliminate constrained degrees of freedom
+  // only by looking at both the system matrix and corresponding right
+  // right hand side at the same time. This may become a problem when
+  // dealing with non-zero Dirichlet boundary conditions, though we
+  // do not do this here in the current program.
   template<int dim>
   void HeatEquation<dim>::setup_system()
   {
@@ -254,7 +278,10 @@ namespace Step26
   }
 
 
-
+  // @sect4{<code>HeatEquation::solve_time_step</code>}
+  //
+  // The next function is the one that solves the actual linear system
+  // for a single time step. There is nothing surprising here:
   template<int dim>
   void HeatEquation<dim>::solve_time_step()
   {
@@ -264,7 +291,8 @@ namespace Step26
     PreconditionSSOR<> preconditioner;
     preconditioner.initialize(system_matrix, 1.0);
 
-    cg.solve(system_matrix, solution, system_rhs, preconditioner);
+    cg.solve(system_matrix, solution, system_rhs,
+             preconditioner);
 
     constraints.distribute(solution);
 
@@ -274,6 +302,9 @@ namespace Step26
 
 
 
+  // @sect4{<code>HeatEquation::output_results</code>}
+  //
+  // Neither is there anything new in generating graphical output:
   template<int dim>
   void HeatEquation<dim>::output_results() const
   {
@@ -292,40 +323,29 @@ namespace Step26
   }
 
 
-  // @sect4{BoussinesqFlowProblem::refine_mesh}
+  // @sect4{<code>HeatEquation::refine_mesh</code>}
   //
-  // This function takes care of the adaptive mesh refinement. The three tasks
+  // This function is the interesting part of the program. It takes care of
+  // the adaptive mesh refinement. The three tasks
   // this function performs is to first find out which cells to
   // refine/coarsen, then to actually do the refinement and eventually
   // transfer the solution vectors between the two different grids. The first
   // task is simply achieved by using the well-established Kelly error
-  // estimator on the temperature (it is the temperature we're mainly
-  // interested in for this program, and we need to be accurate in regions of
-  // high temperature gradients, also to not have too much numerical
-  // diffusion). The second task is to actually do the remeshing. That
-  // involves only basic functions as well, such as the
+  // estimator on the solution. The second task is to actually do the
+  // remeshing. That involves only basic functions as well, such as the
   // <code>refine_and_coarsen_fixed_fraction</code> that refines those cells
-  // with the largest estimated error that together make up 80 per cent of the
+  // with the largest estimated error that together make up 60 per cent of the
   // error, and coarsens those cells with the smallest error that make up for
-  // a combined 10 per cent of the error.
+  // a combined 40 per cent of the error. Note that for problems such as the
+  // current one where the areas where something is going on are shifting
+  // around, we want to aggressively coarsen so that we can move cells
+  // around to where it is necessary.
   //
-  // If implemented like this, we would get a program that will not make much
-  // progress: Remember that we expect temperature fields that are nearly
-  // discontinuous (the diffusivity $\kappa$ is very small after all) and
-  // consequently we can expect that a freely adapted mesh will refine further
-  // and further into the areas of large gradients. This decrease in mesh size
-  // will then be accompanied by a decrease in time step, requiring an
-  // exceedingly large number of time steps to solve to a given final time. It
-  // will also lead to meshes that are much better at resolving
-  // discontinuities after several mesh refinement cycles than in the
-  // beginning.
-  //
-  // In particular to prevent the decrease in time step size and the
-  // correspondingly large number of time steps, we limit the maximal
-  // refinement depth of the mesh. To this end, after the refinement indicator
-  // has been applied to the cells, we simply loop over all cells on the
-  // finest level and unselect them from refinement if they would result in
-  // too high a mesh level.
+  // As already discussed in the introduction, too small a mesh leads to
+  // too small a time step, whereas too large a mesh leads to too little
+  // resolution. Consequently, after the first two steps, we have two
+  // loops that limit refinement and coarsening to an allowable range of
+  // cells:
   template <int dim>
   void HeatEquation<dim>::refine_mesh (const unsigned int min_grid_level,
                                        const unsigned int max_grid_level)
@@ -341,6 +361,7 @@ namespace Step26
     GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
                                                        estimated_error_per_cell,
                                                        0.6, 0.4);
+
     if (triangulation.n_levels() > max_grid_level)
       for (typename Triangulation<dim>::active_cell_iterator
            cell = triangulation.begin_active(max_grid_level);
@@ -357,53 +378,48 @@ namespace Step26
     // SolutionTransfer class and we have to prepare the solution vectors that
     // should be transferred to the new grid (we will lose the old grid once
     // we have done the refinement so the transfer has to happen concurrently
-    // with refinement). What we definitely need are the current and the old
-    // temperature (BDF-2 time stepping requires two old solutions). Since the
-    // SolutionTransfer objects only support to transfer one object per dof
-    // handler, we need to collect the two temperature solutions in one data
-    // structure. Moreover, we choose to transfer the Stokes solution, too,
-    // since we need the velocity at two previous time steps, of which only
-    // one is calculated on the fly.
+    // with refinement). At the point where we call this function, we will
+    // have just computed the solution, so we no longer need the old_solution
+    // variable (it will be overwritten by the solution just after the mesh
+    // may have been refined, i.e., at the end of the time step; see below).
+    // In other words, we only need the one solution vector, and we copy it
+    // to a temporary object where it is safe from being reset when we further
+    // down below call <code>setup_system()</code>.
     //
-    // Consequently, we initialize two SolutionTransfer objects for the Stokes
-    // and temperature DoFHandler objects, by attaching them to the old dof
-    // handlers. With this at place, we can prepare the triangulation and the
-    // data vectors for refinement (in this order).
-    std::vector<Vector<double> > x_solution (2);
-    x_solution[0] = solution;
-    x_solution[1] = old_solution;
-
+    // Consequently, we initialize a SolutionTransfer object by attaching
+    // it to the old DoF handler. We then prepare the triangulation and the
+    // data vector for refinement (in this order).
     SolutionTransfer<dim> solution_trans(dof_handler);
 
+    Vector<double> previous_solution;
+    previous_solution = solution;
     triangulation.prepare_coarsening_and_refinement();
-    solution_trans.prepare_for_coarsening_and_refinement(x_solution);
+    solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
 
     // Now everything is ready, so do the refinement and recreate the dof
     // structure on the new grid, and initialize the matrix structures and the
-    // new vectors in the <code>setup_dofs</code> function. Next, we actually
-    // perform the interpolation of the solutions between the grids. We create
-    // another copy of temporary vectors for temperature (now corresponding to
-    // the new grid), and let the interpolate function do the job. Then, the
-    // resulting array of vectors is written into the respective vector member
-    // variables. For the Stokes vector, everything is just the same &ndash;
-    // except that we do not need another temporary vector since we just
-    // interpolate a single vector. In the end, we have to tell the program
-    // that the matrices and preconditioners need to be regenerated, since the
-    // mesh has changed.
+    // new vectors in the <code>setup_system</code> function. Next, we actually
+    // perform the interpolation of the solution from old to new grid.
     triangulation.execute_coarsening_and_refinement ();
     setup_system ();
 
-    std::vector<Vector<double> > tmp (2);
-    tmp[0].reinit (solution);
-    tmp[1].reinit (solution);
-    solution_trans.interpolate(x_solution, tmp);
-
-    solution = tmp[0];
-    old_solution = tmp[1];
+    solution_trans.interpolate(previous_solution, solution);
   }
 
 
 
+  // @sect4{<code>HeatEquation::run</code>}
+  //
+  // This is the main driver of the program, where we loop over all
+  // time steps. At the top of the function, we set the number of
+  // initial global mesh refinements and the number of initial cycles of
+  // adaptive mesh refinement by repeating the first time step a few
+  // times. Then we create a mesh, initialize the various objects we will
+  // work with, set a label for where we should start when re-running
+  // the first time step, and interpolate the initial solution onto
+  // out mesh (we choose the zero function here, which of course we could
+  // do in a simpler way by just setting the solution vector to zero). We
+  // also output the initial time step once.
   template<int dim>
   void HeatEquation<dim>::run()
   {
@@ -422,6 +438,10 @@ namespace Step26
 
 start_time_iteration:
 
+    tmp.reinit (solution.size());
+    forcing_terms.reinit (solution.size());
+
+
     VectorTools::interpolate(dof_handler,
                              ZeroFunction<dim>(),
                              old_solution);
@@ -432,6 +452,12 @@ start_time_iteration:
 
     output_results();
 
+    // Then we start the main loop until the computed time exceeds our
+    // end time of 0.5. The first task is to build the right hand
+    // side of the linear system we need to solve in each time step.
+    // Recall that it contains the term $MU^{n-1}-(1-\theta)k_n AU^{n-1}$.
+    // We put these terms into the variable system_rhs, with the
+    // help of a temporary vector:
     while (time <= 0.5)
       {
         time += time_step;
@@ -440,14 +466,18 @@ start_time_iteration:
         std::cout << "Time step " << timestep_number << " at t=" << time
                   << std::endl;
 
-        tmp.reinit (solution.size());
-        forcing_terms.reinit (solution.size());
-
         mass_matrix.vmult(system_rhs, old_solution);
 
         laplace_matrix.vmult(tmp, old_solution);
         system_rhs.add(-(1 - theta) * time_step, tmp);
 
+        // The second piece is to compute the contributions of the source
+        // terms. This corresponds to the term $k_n
+        // \left[ (1-\theta)F^{n-1} + \theta F^n \right]$. The following
+        // code calls VectorTools::create_right_hand_side to compute the
+        // vectors $F$, where we set the time of the right hand side
+        // (source) function before we evaluate it. The result of this
+        // all ends up in the forcing_terms variable:
         RightHandSide<dim> rhs_function;
         rhs_function.set_time(time);
         VectorTools::create_right_hand_side(dof_handler,
@@ -465,8 +495,25 @@ start_time_iteration:
 
         forcing_terms.add(time_step * (1 - theta), tmp);
 
+        // Next, we add the forcing terms to the ones that
+        // come from the time stepping, and also build the matrix
+        // $M+k_n\theta A$ that we have to invert in each time step.
+        // The final piece of these operations is to eliminate
+        // hanging node constrained degrees of freedom from the
+        // linear system:
         system_rhs += forcing_terms;
 
+        system_matrix.copy_from(mass_matrix);
+        system_matrix.add(theta * time_step, laplace_matrix);
+
+        constraints.condense (system_matrix, system_rhs);
+
+        // There is one more operation we need to do before we
+        // can solve it: boundary values. To this end, we create
+        // a boundary value object, set the proper time to the one
+        // of the current time step, and evaluate it as we have
+        // done many times before. The result is used to also
+        // set the correct boundary values in the linear system:
         {
           BoundaryValues<dim> boundary_values_function;
           boundary_values_function.set_time(time);
@@ -477,20 +524,27 @@ start_time_iteration:
                                                    boundary_values_function,
                                                    boundary_values);
 
-          system_matrix.copy_from(mass_matrix);
-          system_matrix.add(theta * time_step, laplace_matrix);
           MatrixTools::apply_boundary_values(boundary_values,
                                              system_matrix,
                                              solution,
                                              system_rhs);
         }
 
-        constraints.condense (system_matrix, system_rhs);
-
+        // With this out of the way, all we have to do is solve the
+        // system, generate graphical data, and...
         solve_time_step();
 
         output_results();
 
+        // ...take care of mesh refinement. Here, what we want to do is
+        // (i) refine the requested number of times at the very beginning
+        // of the solution procedure, after which we jump to the top to
+        // restart the time iteration, (ii) refine every fifth time
+        // step after that.
+        //
+        // The time loop and, indeed, the main part of the program ends
+        // with starting into the next time step by setting old_solution
+        // to the solution we have just computed.
         if ((timestep_number == 1) &&
             (pre_refinement_step < n_adaptive_pre_refinement_steps))
           {
@@ -498,19 +552,32 @@ start_time_iteration:
                          initial_global_refinement + n_adaptive_pre_refinement_steps);
             ++pre_refinement_step;
 
+            tmp.reinit (solution.size());
+            forcing_terms.reinit (solution.size());
+
             std::cout << std::endl;
 
             goto start_time_iteration;
           }
         else if ((timestep_number > 0) && (timestep_number % 5 == 0))
-          refine_mesh (initial_global_refinement,
-                       initial_global_refinement + n_adaptive_pre_refinement_steps);
+          {
+            refine_mesh (initial_global_refinement,
+                         initial_global_refinement + n_adaptive_pre_refinement_steps);
+            tmp.reinit (solution.size());
+            forcing_terms.reinit (solution.size());
+          }
 
         old_solution = solution;
       }
   }
 }
 
+
+// @sect3{The <code>main</code> function}
+//
+// Having made it this far,  there is, again, nothing
+// much to discuss for the main function of this
+// program: it looks like all such functions since step-6.
 int main()
 {
   try

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.