]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Give it a once over. 196/head
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 13 Oct 2014 23:59:45 +0000 (18:59 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 13 Oct 2014 23:59:45 +0000 (18:59 -0500)
No functional changes, just updates to documentation and style.

examples/step-52/doc/intro.dox
examples/step-52/doc/results.dox
examples/step-52/doc/tooltip
examples/step-52/step-52.cc

index 2998354decded367e30b15da76dd194f086284c6..ba3cc54946780624f57082b54fcbb02687f32700 100644 (file)
@@ -2,7 +2,7 @@
 
 <i>This program was contributed by Bruno Turcksin and Damien Lebrun-Grandie.</i>
 
-<b>Note:</b> In order to run this program, deal.II must be configured to use
+@note In order to run this program, deal.II must be configured to use
 the UMFPACK sparse direct solver. Refer to the <a
 href="../../readme.html#umfpack">ReadMe</a> for instructions how to do this.
 
@@ -45,8 +45,12 @@ $[0,b]\times[0,b]$ of the form
 \phi(x,t) = A\sin(\omega t)(bx-x^2).
 @f}
 By using quadratic finite elements, we can represent this function exactly at
-any particular time, and all the error will be due to the time discretization. We
-impose the following boundary conditions: homogeneous Dirichlet for $x=0$ and
+any particular time, and all the error will be due to the time
+discretization. We do this because it is then easy to observe the order of
+convergence of the various time stepping schemes we will consider, without
+having to separate spatial and temporal errors.
+
+We impose the following boundary conditions: homogeneous Dirichlet for $x=0$ and
 $x=b$ and homogeneous Neumann conditions for $y=0$ and $y=b$. We choose the
 source term so that the corresponding solution is
 in fact of the form stated above:
@@ -54,9 +58,15 @@ in fact of the form stated above:
 S=A\left(\frac{1}{v}\omega \cos(\omega t)(bx -x^2) + \sin(\omega t)
 \left(\Sigma_a (bx-x^2)+2D\right) \right).
 @f}
-Because the solution is a sine, we know that $\phi\left(x,\pi\right) = 0$.
+Because the solution is a sine in time, we know that the exact solution
+satisfies $\phi\left(x,\pi\right) = 0$.
 Therefore, the error at time $t=\pi$ is simply the norm of the numerical
-solution and is particularly easily evaluated.
+solution, i.e., $\|e(\cdot,t=\pi)\|_{L_2} = \|\phi_h(\cdot,t=\pi)\|_{L_2}$,
+and is particularly easily evaluated. In the code, we evaluate the $l_2$ norm
+of the vector of nodal values of $\phi_h$ instead of the $L_2$ norm of the
+associated spatial function, since the former is simpler to compute; however,
+on uniform meshes, the two are just related by a constant and we can
+consequently observe the temporal convergence order with either.
 
 
 <h3>Runge-Kutta methods</h3>
@@ -85,7 +95,7 @@ spatially discretized equation has the form
 M\frac{dU}{dt} = f(t,U),
 @f}
 where $M$ is the mass matrix and $f(t,U)$ is the spatially discretized version
-of the $q(t,u(x,t))$ (where $q$ is typically the place where spatial
+of $q(t,u(x,t))$ (where $q$ is typically the place where spatial
 derivatives appear, but this is not of much concern for the moment given that
 we only consider time derivatives). In other words, this form fits the general
 scheme above if we write 
@@ -113,11 +123,16 @@ Butcher tableaus.)
 At the time of the writing of this tutorial, the methods implemented in
 deal.II can be divided in three categories:
 <ol>
-<li> explicit Runge-Kutta
-<li> embedded (or adaptive) Runge-Kutta
-<li> implicit Runge-Kutta
+<li> Explicit Runge-Kutta; in order for a method to be explicit, it is
+necessary that in the formula above defining $k_i$, $k_i$ does not appear
+on the right hand side. In other words, these methods have to satisfy
+$a_{ii}=0, i=1,\ldots,s$.
+<li> Embedded (or adaptive) Runge-Kutta; we will discuss their properties below.
+<li> Implicit Runge-Kutta; this class of methods require the solution of a
+possibly nonlinear system the stages $k_i$ above, i.e., they have
+$a_{ii}\neq 0$ for at least one of the stages $i=1,\ldots,s$.
 </ol> 
-Many well known time stepping schemes that one does not typically associated
+Many well known time stepping schemes that one does not typically associate
 with the names Runge or Kutta can in fact be written in a way so that they,
 too, can be expressed in these categories. They oftentimes represent the
 lowest-order members of these families.
@@ -159,13 +174,13 @@ form $\alpha y = f(t,y)$
 for $y$ in each (sub-)timestep. Internally, this is
 done using a Newton-type method and, consequently, they require that the user
 provide functions that can evaluate $M^{-1}f(t,y)$ and
-$\left(I-\Delta t M^{-1} \frac{\partial f}{\partial y}\right)^{-1}$ or equivalently 
-$\left(M - \Delta t \frac{\partial f}{\partial y}\right)^{-1} M$. 
+$\left(I-\tau M^{-1} \frac{\partial f}{\partial y}\right)^{-1}$ or equivalently 
+$\left(M - \tau \frac{\partial f}{\partial y}\right)^{-1} M$. 
 
 The particular form of this operator results from the fact that each Newton
 step requires the solution of an equation of the form
 @f{align*}
-  \left(M - \Delta t \frac{\partial f}{\partial y}\right) \Delta y
+  \left(M - \tau \frac{\partial f}{\partial y}\right) \Delta y
   = -M h(t,y)
 @f}
 for some (given) $h(t,y)$. Implicit methods are 
@@ -174,12 +189,17 @@ course affect the <i>accuracy</i> of the solution, even if the numerical
 solution remains stable and bounded.
 
 Methods in this class include backward Euler, implicit midpoint,
-Crank-Nicolson, and a two stage SDIRK. 
+Crank-Nicolson, and the two stage SDIRK method (short for "singly diagonally
+implicit Runge-Kutta", a term coined to indicate that the diagonal elements
+$a_{ii}$ defining the time stepping method are all equal; this property
+allows for the Newton matrix $I-\tau M^{-1}\frac{\partial f}{\partial y}$ to
+be re-used between stages because $\tau$ is the same every time).
 
 
 <h3>Spatially discrete formulation</h3>
 
-By expanding the solution as always using shape functions $\psi_j$ and writing
+By expanding the solution of our model problem
+as always using shape functions $\psi_j$ and writing
 @f{eqnarray*}
 \phi_h(x,t) = \sum_j U_j(t) \psi_j(x),
 @f}
@@ -197,13 +217,16 @@ where
 @f}
 See also step-24 and step-26 to understand how we arrive here.
 %Boundary terms are not necessary due to the chosen boundary conditions for
-the current problem. To use the Runge-Kutta methods, we can then recast this
+the current problem. To use the Runge-Kutta methods, we recast this
 as follows:
 @f{eqnarray*}
-f(y) = -{\cal D}y - {\cal A}y + {\cal S}
+f(y) = -{\cal D}y - {\cal A}y + {\cal S}.
 @f}
 In the code, we will need to be able to evaluate this function $f(U)$ along
-with its derivative.
+with its derivative,
+@f{eqnarray*}
+\frac{\partial f}{\partial y} = -{\cal D} - {\cal A}.
+@f}
 
 
 <h3>Notes on the testcase</h3>
@@ -219,7 +242,7 @@ important to remember the following issues:
 <ol>
 <li> You will need to project the solution to the new mesh when the mesh is changed. Of course,
      the mesh 
-     used should be the same at the beginning and at the end of the time step,
+     used should be the same from the beginning to the end of each time step,
      a question that arises because Runge-Kutta methods use multiple
      evaluations of the equations within each time step.
 <li> You will need to update the mass matrix and its inverse every time the
index 72a5cb05181bdc4f59c3bf274713e01bab5fad20..4c2d9ea0b0d90ce9e98889f434a5d69841035e5c 100644 (file)
@@ -1,10 +1,10 @@
 <h1>Results</h1>
 
-The point of this program is less to show partcular results, but instead to
-show how it is done and that we have already demonstrated using the code
-above. Consequently, the output the program yields is relatively sparse and
-consists only of the console output and the solutions given in
-VTU format for visualization.
+The point of this program is less to show particular results, but instead to
+show how it is done. This we have already demonstrated simply by discussing
+the code above. Consequently, the output the program yields is relatively
+sparse and consists only of the console output and the solutions given in VTU
+format for visualization.
 
 The console output contains both errors and, for some of the methods, the
 number of steps they performed:
@@ -20,7 +20,7 @@ Implicit Midpoint:        error=0.00862702
 Crank-Nicolson:           error=0.00862675
 SDIRK:                    error=0.0042349
 
-Embedded explicit methods:
+Embedded %explicit methods:
 Heun-Euler:               error=0.0073012
                 steps performed=284
 Bogacki-Shampine:         error=0.000207511
index d35ea621a8ae874314c8925818eef4416d9186e3..7754b44e6b16f44c034e8ad506e59ac3e8a062c2 100644 (file)
@@ -1 +1 @@
-Time-dependent diffusion equation. Neutron transport
+Time-dependent diffusion equation. Time stepping methods.
index 2c5cd4ed82f627d486767f14509f96708ff04d19..51c7d9e302c2ba02cf1072edbceb61052d229cd0 100644 (file)
@@ -61,17 +61,18 @@ namespace Step52
 {
   using namespace dealii;
 
-  // @sect3{Diffusion}
-
-  // The next piece is the declaration of the main class. Most of the functions in
-  // this class are not new and have been explained in previous tutorials. The
-  // only interesting functions are <code>evaluate_diffusion</code> and
-  // <code>id_minus_tau_J_inverse</code>. <code>evaluate_diffusion</code> evaluates the
-  // diffusion equation, $M^{-1}(f(t,y))$, at a given time, for a given $\tau$
-  // and a given $y$. <code>id_minus_tau_J_inverse</code> evaluates
-  // $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$ or
-  // equivalently $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$ at
-  // a given time, for a given $\tau$ and $y$. This function is needed when an
+  // @sect3{The <code>Diffusion</code> class}
+
+  // The next piece is the declaration of the main class. Most of the
+  // functions in this class are not new and have been explained in previous
+  // tutorials. The only interesting functions are
+  // <code>evaluate_diffusion()</code> and
+  // <code>id_minus_tau_J_inverse()</code>. <code>evaluate_diffusion()</code>
+  // evaluates the diffusion equation, $M^{-1}(f(t,y))$, at a given time and a
+  // given $y$. <code>id_minus_tau_J_inverse()</code> evaluates $\left(I-\tau
+  // M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$ or equivalently
+  // $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$ at a given
+  // time, for a given $\tau$ and $y$. This function is needed when an
   // implicit method is used.
   class Diffusion
   {
@@ -85,40 +86,45 @@ namespace Step52
 
     void assemble_system();
 
-    double get_source(double time,const Point<2> &point) const;
+    double get_source (const double time,
+                       const Point<2> &point) const;
 
-    Vector<double> evaluate_diffusion(const double time, const Vector<double> &y) const;
+    Vector<double> evaluate_diffusion (const double time,
+                                       const Vector<double> &y) const;
 
-    Vector<double> id_minus_tau_J_inverse(const double time,
-                                          const double tau,
-                                          const Vector<double> &y);
+    Vector<double> id_minus_tau_J_inverse (const double time,
+                                           const double tau,
+                                           const Vector<double> &y);
 
-    void output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const;
+    void output_results (const unsigned int time_step,
+                         TimeStepping::runge_kutta_method method) const;
 
-    // The next three functions are the driver for the explicit methods, the
+    // The next three functions are the drivers for the explicit methods, the
     // implicit methods, and the embedded explicit methods respectively. The
     // driver function for embedded explicit methods returns the number of
-    // steps executed.
-    void explicit_method(TimeStepping::runge_kutta_method method,
-                         const unsigned int n_time_steps,
-                         const double       initial_time,
-                         const double       final_time);
+    // steps executed given that it only takes the number of time steps passed
+    // as an argument as a hint, but internally computed the optimal time step
+    // itself.
+    void explicit_method (const TimeStepping::runge_kutta_method method,
+                          const unsigned int                     n_time_steps,
+                          const double                           initial_time,
+                          const double                           final_time);
 
-    void implicit_method(TimeStepping::runge_kutta_method method,
-                         const unsigned int n_time_steps,
-                         const double       initial_time,
-                         const double       final_time);
+    void implicit_method (const TimeStepping::runge_kutta_method method,
+                          const unsigned int                     n_time_steps,
+                          const double                           initial_time,
+                          const double                           final_time);
 
-    unsigned int embedded_explicit_method(TimeStepping::runge_kutta_method method,
-                                          const unsigned int n_time_steps,
-                                          const double initial_time,
-                                          const double final_time);
+    unsigned int embedded_explicit_method (const TimeStepping::runge_kutta_method method,
+                                           const unsigned int                     n_time_steps,
+                                           const double                     initial_time,
+                                           const double                     final_time);
 
 
     unsigned int                 fe_degree;
 
     double                       diffusion_coefficient;
-    double                       absorption_xs;
+    double                       absorption_cross_section;
 
     Triangulation<2>             triangulation;
 
@@ -146,17 +152,17 @@ namespace Step52
     :
     fe_degree(2),
     diffusion_coefficient(1./30.),
-    absorption_xs(1.),
+    absorption_cross_section(1.),
     fe(fe_degree),
     dof_handler(triangulation)
   {}
 
 
 
-  // @sect5{<code>Diffusion::setup_system</code>}
+  // @sect4{<code>Diffusion::setup_system</code>}
   // Now, we create the constraint matrix and the sparsity pattern. Then, we
   // initialize the matrices and the solution vector.
-  void Diffusion::setup_system()
+  void Diffusion::setup_system ()
   {
     dof_handler.distribute_dofs(fe);
 
@@ -175,12 +181,18 @@ namespace Step52
 
 
 
-  // @sect5{<code>Diffusion::assemble_system</code>}
-  // In this function, we compute
-  // $-\int D \nabla b_i \cdot \nabla b_j d\boldsymbol{r} - \int \Sigma_a b_i b_j d\boldsymbol{r}$
-  // and the mass matrix $\int b_i b_j d\boldsymbol{r}$. The mass matrix is then
-  // inverted using a direct solver.
-  void Diffusion::assemble_system()
+  // @sect4{<code>Diffusion::assemble_system</code>} In this function, we
+  // compute $-\int D \nabla b_i \cdot \nabla b_j d\boldsymbol{r} - \int
+  // \Sigma_a b_i b_j d\boldsymbol{r}$ and the mass matrix $\int b_i b_j
+  // d\boldsymbol{r}$. The mass matrix is then inverted using a direct solver;
+  // the <code>inverse_mass_matrix</code> variable will then store the inverse
+  // of the mass matrix so that $M^{-1}$ can be applied to a vector using the
+  // <code>vmult()</code> function of that object. (Internally, UMFPACK does
+  // not really store the inverse of the matrix, but its LU factors; applying
+  // the inverse matrix is then equivalent to doing one forward and one
+  // backward solves with these two factors, which has the same complexity as
+  // applying an explicit inverse of the matrix).
+  void Diffusion::assemble_system ()
   {
     system_matrix = 0.;
     mass_matrix = 0.;
@@ -214,13 +226,16 @@ namespace Step52
           for (unsigned int i=0; i<dofs_per_cell; ++i)
             for (unsigned int j=0; j<dofs_per_cell; ++j)
               {
-                cell_matrix(i,j) += ((-diffusion_coefficient * fe_values.shape_grad(i,q_point) *
-                                      fe_values.shape_grad(j,q_point) - absorption_xs *
-                                      fe_values.shape_value(i,q_point) * fe_values.shape_value(j,q_point)) *
+                cell_matrix(i,j) += ((-diffusion_coefficient *
+                                      fe_values.shape_grad(i,q_point) *
+                                      fe_values.shape_grad(j,q_point)
+                                      - absorption_cross_section *
+                                      fe_values.shape_value(i,q_point) *
+                                      fe_values.shape_value(j,q_point)) *
                                      fe_values.JxW(q_point));
-                cell_mass_matrix(i,j) +=  fe_values.shape_value(i,q_point) *
-                                          fe_values.shape_value(j,q_point) *
-                                          fe_values.JxW(q_point);
+                cell_mass_matrix(i,j) += fe_values.shape_value(i,q_point) *
+                                         fe_values.shape_value(j,q_point) *
+                                         fe_values.JxW(q_point);
               }
 
         cell->get_dof_indices(local_dof_indices);
@@ -234,32 +249,41 @@ namespace Step52
 
 
 
-  // @sect5{<code>Diffusion::get_source</code>}
+  // @sect4{<code>Diffusion::get_source</code>}
   //
-  // In this function, the source for a given time and a given point is
-  // computed.
-  double Diffusion::get_source(double time,const Point<2> &point) const
+  // In this function, the source term of the equation for a given time and a
+  // given point is computed.
+  double Diffusion::get_source (const double time,
+                                const Point<2> &point) const
   {
-    const double pi = 3.14159265358979323846;
     const double intensity = 10.;
-    const double frequency = pi/10.;
+    const double frequency = numbers::PI/10.;
     const double b = 5.;
     const double x = point(0);
-    double source = 0.;
 
-    source = intensity*(frequency*std::cos(frequency*time)*(b*x-x*x) + std::sin(frequency*time) *
-                        (absorption_xs*(b*x-x*x)+2.*diffusion_coefficient));
-
-    return source;
+    return intensity* (frequency*std::cos(frequency*time)*(b*x-x*x)
+                       +
+                       std::sin(frequency*time) * (absorption_cross_section*(b*x-x*x)
+                                                   +
+                                                   2.*diffusion_coefficient));
   }
 
 
 
-  // @sect5{<code>Diffusion:evaluate_diffusion</code>}
+  // @sect4{<code>Diffusion:evaluate_diffusion</code>}
   //
-  // Now, the weak form of the diffusion equation is evaluated at a given
-  // time $t$ and for a given vector $y$.
-  Vector<double> Diffusion::evaluate_diffusion(const double time, const Vector<double> &y) const
+  // Next, we evaluate the weak form of the diffusion equation at a given time
+  // $t$ and for a given vector $y$. In other words, as outlined in the
+  // introduction, we evaluate $M^{-1}(-{\cal D}y - {\cal A}y + {\cal
+  // S})$. For this, we have to apply the matrix $-{\cal D} - {\cal A}$
+  // (previously computed and stored in the variable
+  // <code>system_matrix</code>) to $y$ and then add the source term which we
+  // integrate as we usually do. (Integrating up the solution could be done
+  // using VectorTools::create_right_hand_side() if you wanted to save a few
+  // lines of code, or wanted to take advantage of doing the integration in
+  // parallel.) The result is then multiplied by $M^{-1}$.
+  Vector<double> Diffusion::evaluate_diffusion (const double time,
+                                                const Vector<double> &y) const
   {
     Vector<double> tmp(dof_handler.n_dofs());
     tmp = 0.;
@@ -267,7 +291,8 @@ namespace Step52
 
     const QGauss<2> quadrature_formula(fe_degree+1);
 
-    FEValues<2> fe_values(fe, quadrature_formula,
+    FEValues<2> fe_values(fe,
+                          quadrature_formula,
                           update_values | update_quadrature_points | update_JxW_values);
 
 
@@ -290,18 +315,21 @@ namespace Step52
 
         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
           {
-            double source = get_source(time,fe_values.quadrature_point(q_point)) ;
+            const double source = get_source(time,
+                                             fe_values.quadrature_point(q_point));
             for (unsigned int i=0; i<dofs_per_cell; ++i)
-              cell_source(i) += source * fe_values.shape_value(i,q_point) *
+              cell_source(i) += source *
+                                fe_values.shape_value(i,q_point) *
                                 fe_values.JxW(q_point);
           }
 
         cell->get_dof_indices(local_dof_indices);
 
-        constraint_matrix.distribute_local_to_global(cell_source,local_dof_indices,tmp);
+        constraint_matrix.distribute_local_to_global(cell_source,
+                                                     local_dof_indices,
+                                                     tmp);
       }
 
-
     Vector<double> value(dof_handler.n_dofs());
     inverse_mass_matrix.vmult(value,tmp);
 
@@ -309,28 +337,30 @@ namespace Step52
   }
 
 
-  // @sect5{<code>Diffusion::id_minus_tau_J_inverse</code>}
+  // @sect4{<code>Diffusion::id_minus_tau_J_inverse</code>}
   //
   // We compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M$. This
   // is done in several steps:
   //   - compute $M-\tau \frac{\partial f}{\partial y}$
-  //   - inverse the matrix to get $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1}$
+  //   - invert the matrix to get $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1}$
   //   - compute $tmp=My$
-  //   - compute $\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp = \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} My$.
-  Vector<double> Diffusion::id_minus_tau_J_inverse(const double time, const double tau,
-                                                   const Vector<double> &y)
+  //   - compute $z=\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp = \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} My$
+  //   - return z.
+  Vector<double> Diffusion::id_minus_tau_J_inverse (const double time,
+                                                    const double tau,
+                                                    const Vector<double> &y)
   {
-    Vector<double> tmp(dof_handler.n_dofs());
-    Vector<double> result(y);
     SparseDirectUMFPACK inverse_mass_minus_tau_Jacobian;
 
     mass_minus_tau_Jacobian.copy_from(mass_matrix);
-    mass_minus_tau_Jacobian.add(-tau,system_matrix);
+    mass_minus_tau_Jacobian.add(-tau, system_matrix);
 
     inverse_mass_minus_tau_Jacobian.initialize(mass_minus_tau_Jacobian);
 
-    mass_matrix.vmult(tmp,y);
+    Vector<double> tmp(dof_handler.n_dofs());
+    mass_matrix.vmult(tmp, y);
 
+    Vector<double> result(y);
     inverse_mass_minus_tau_Jacobian.vmult(result,tmp);
 
     return result;
@@ -338,71 +368,75 @@ namespace Step52
 
 
 
-  // @sect5{<code>Diffusion::output_results</code>}
+  // @sect4{<code>Diffusion::output_results</code>}
   //
-  // We output the solution in vtu files.
-  void Diffusion::output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const
+  // The following function then outputs the solution in vtu files indexed by
+  // the number of the time step and the name of the time stepping method. Of
+  // course, the (exact) result should really be the same for all time
+  // stepping method, but the output here at least allows us to compare them.
+  void Diffusion::output_results (const unsigned int time_step,
+                                  TimeStepping::runge_kutta_method method) const
   {
     std::string method_name;
 
     switch (method)
       {
-      case TimeStepping::FORWARD_EULER :
+      case TimeStepping::FORWARD_EULER:
       {
         method_name = "forward_euler";
         break;
       }
-      case TimeStepping::RK_THIRD_ORDER :
+      case TimeStepping::RK_THIRD_ORDER:
       {
         method_name = "rk3";
         break;
       }
-      case TimeStepping::RK_CLASSIC_FOURTH_ORDER :
+      case TimeStepping::RK_CLASSIC_FOURTH_ORDER:
       {
         method_name = "rk4";
         break;
       }
-      case TimeStepping::BACKWARD_EULER :
+      case TimeStepping::BACKWARD_EULER:
       {
         method_name = "backward_euler";
         break;
       }
-      case TimeStepping::IMPLICIT_MIDPOINT :
+      case TimeStepping::IMPLICIT_MIDPOINT:
       {
         method_name = "implicit_midpoint";
         break;
       }
-      case TimeStepping::SDIRK_TWO_STAGES :
+      case TimeStepping::SDIRK_TWO_STAGES:
       {
         method_name = "sdirk";
         break;
       }
-      case TimeStepping::HEUN_EULER :
+      case TimeStepping::HEUN_EULER:
       {
         method_name = "heun_euler";
         break;
       }
-      case TimeStepping::BOGACKI_SHAMPINE :
+      case TimeStepping::BOGACKI_SHAMPINE:
       {
         method_name = "bocacki_shampine";
         break;
       }
-      case TimeStepping::DOPRI :
+      case TimeStepping::DOPRI:
       {
         method_name = "dopri";
         break;
       }
-      case TimeStepping::FEHLBERG :
+      case TimeStepping::FEHLBERG:
       {
         method_name = "fehlberg";
         break;
       }
-      case TimeStepping::CASH_KARP :
+      case TimeStepping::CASH_KARP:
       {
         method_name = "cash_karp";
         break;
       }
-      default :
+      default:
       {
         break;
       }
@@ -423,17 +457,18 @@ namespace Step52
   }
 
 
-  // @sect5{<code>Diffusion::explicit_method</code>}
-  // This function is the driver for all the explicit method. It calls
-  // <code>evolve_one_time_step</code> which performs one time step.
-  // <code>evolve_one_time_step</code> needs to evaluate $M^{-1}(f(t,y))$,
-  // i.e, it needs <code>evaluate_diffusion</code>. Because
-  // <code>evaluate_diffusion</code> is a member function, it needs to be bind to
-  // $this$. Finally, the solution is output every 10 time steps.
-  void Diffusion::explicit_method(TimeStepping::runge_kutta_method method,
-                                  const unsigned int                n_time_steps,
-                                  const double                      initial_time,
-                                  const double                      final_time)
+  // @sect4{<code>Diffusion::explicit_method</code>}
+  //
+  // This function is the driver for all the explicit methods. It calls
+  // <code>evolve_one_time_step</code> which performs one time step.  For
+  // explicit methods, <code>evolve_one_time_step</code> needs to evaluate
+  // $M^{-1}(f(t,y))$, i.e, it needs <code>evaluate_diffusion</code>. Because
+  // <code>evaluate_diffusion</code> is a member function, it needs to be bound
+  // to <code>this</code>. Finally, the solution is output every 10 time steps.
+  void Diffusion::explicit_method (const TimeStepping::runge_kutta_method method,
+                                   const unsigned int                     n_time_steps,
+                                   const double                           initial_time,
+                                   const double                           final_time)
   {
     const double time_step = (final_time-initial_time)/static_cast<double> (n_time_steps);
     double time = initial_time;
@@ -444,7 +479,10 @@ namespace Step52
     for (unsigned int i=0; i<n_time_steps; ++i)
       {
         time = explicit_runge_kutta.evolve_one_time_step(
-                 std_cxx11::bind(&Diffusion::evaluate_diffusion,this,std_cxx11::_1,std_cxx11::_2),
+                 std_cxx11::bind(&Diffusion::evaluate_diffusion,
+                                 this,
+                                 std_cxx11::_1,
+                                 std_cxx11::_2),
                  time,time_step,solution);
 
         if ((i+1)%10==0)
@@ -454,14 +492,15 @@ namespace Step52
 
 
 
-  // @sect5{<code>Diffusion::implicit_method</code>}
+  // @sect4{<code>Diffusion::implicit_method</code>}
   // This function is equivalent to <code>explicit_method</code> but for implicit
   // methods. When using implicit methods, we need to evaluate $M^{-1}(f(t,y))$
-  // and $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$.
-  void Diffusion::implicit_method(TimeStepping::runge_kutta_method method,
-                                  const unsigned int               n_time_steps,
-                                  const double                     initial_time,
-                                  const double                     final_time)
+  // and $\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}$
+  // for which we use the two member functions previously introduced.
+  void Diffusion::implicit_method (const TimeStepping::runge_kutta_method method,
+                                   const unsigned int                     n_time_steps,
+                                   const double                           initial_time,
+                                   const double                           final_time)
   {
     const double time_step = (final_time-initial_time)/static_cast<double> (n_time_steps);
     double time = initial_time;
@@ -472,8 +511,14 @@ namespace Step52
     for (unsigned int i=0; i<n_time_steps; ++i)
       {
         time = implicit_runge_kutta.evolve_one_time_step(
-                 std_cxx11::bind(&Diffusion::evaluate_diffusion,this,std_cxx11::_1,std_cxx11::_2),
-                 std_cxx11::bind(&Diffusion::id_minus_tau_J_inverse,this,std_cxx11::_1,std_cxx11::_2,
+                 std_cxx11::bind(&Diffusion::evaluate_diffusion,
+                                 this,
+                                 std_cxx11::_1,
+                                 std_cxx11::_2),
+                 std_cxx11::bind(&Diffusion::id_minus_tau_J_inverse,
+                                 this,
+                                 std_cxx11::_1,
+                                 std_cxx11::_2,
                                  std_cxx11::_3),
                  time,time_step,solution);
 
@@ -484,7 +529,7 @@ namespace Step52
 
 
 
-  // @sect5{<code>Diffusion::embedded_explicit_method</code>}
+  // @sect4{<code>Diffusion::embedded_explicit_method</code>}
   // This function is the driver for the embedded explict methods. It requires
   // more parameters:
   //   - coarsen_param: factor multiplying the current time step when the error
@@ -499,10 +544,10 @@ namespace Step52
   // is too large, the time step will be reduced. If the error is below the
   // threshold, a larger time step will be tried for the next time step.
   // </code>delta_t_guess</code> is the guessed time step produced by the embedded method.
-  unsigned int Diffusion::embedded_explicit_method(TimeStepping::runge_kutta_method method,
-                                                   const unsigned int n_time_steps,
-                                                   const double initial_time,
-                                                   const double final_time)
+  unsigned int Diffusion::embedded_explicit_method(const TimeStepping::runge_kutta_method method,
+                                                   const unsigned int                     n_time_steps,
+                                                   const double                           initial_time,
+                                                   const double                           final_time)
   {
     double time_step = (final_time-initial_time)/static_cast<double> (n_time_steps);
     double time = initial_time;
@@ -514,14 +559,21 @@ namespace Step52
     const double coarsen_tol = 1e-5;
     solution = 0.;
 
-    TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> > embedded_explicit_runge_kutta(method,
-        coarsen_param,refine_param,min_delta,max_delta,refine_tol,coarsen_tol);
-    output_results(0,method);
+    TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> >
+    embedded_explicit_runge_kutta(method,
+                                  coarsen_param,
+                                  refine_param,
+                                  min_delta,
+                                  max_delta,
+                                  refine_tol,
+                                  coarsen_tol);
+    output_results (0, method);
+
+    // Now for the time loop. The last time step is chosen such that the final
+    // time is exactly reached.
     unsigned int n_steps=0;
     while (time<final_time)
       {
-        // The last time step is chosend such that the final time is exactly
-        // reached.
         if (time+time_step>final_time)
           time_step = final_time-time;
 
@@ -541,15 +593,17 @@ namespace Step52
 
 
 
-  // @sect5{<code>Diffusion::run</code>}
-  void Diffusion::run()
+  // @sect4{<code>Diffusion::run</code>}
+  //
+  // The following is the main function of the program. At the top, we create
+  // the grid (a [0,5]x[0,5] square) and refine it four times to get a mesh
+  // that has 16 by 16 cells, for a total of 256.  We then set the boundary
+  // indicator to 1 for those parts of the boundary where $x=0$ and $x=5$.
+  void Diffusion::run ()
   {
-    // We create the grid (a [0,5]x[0,5] square) and refine the mesh four times.
-    // The final grid has 16 by 16 cells, for a total of 256.
     GridGenerator::hyper_cube(triangulation, 0., 5.);
     triangulation.refine_global(4);
 
-    // Set the boundary indicator for x=0 and x=5 to 1.
     typename Triangulation<2>::active_cell_iterator
     cell = triangulation.begin_active(),
     endc = triangulation.end();
@@ -564,16 +618,23 @@ namespace Step52
               cell->face(f)->set_boundary_indicator(0);
           }
 
+    // Next, we set up the linear systems and fill them with content so that
+    // they can be used throughout the time stepping process:
     setup_system();
 
     assemble_system();
 
-    unsigned int n_steps = 0;
+    // Finally, we solve the diffusion problem using several of the
+    // Runge-Kutta methods implemented in namespace TimeStepping, each time
+    // outputting the error at the end time. (As explained in the
+    // introduction, since the exact solution is zero at the final time, the
+    // error equals the numerical solution and can be computed by just taking
+    // the $l_2$ norm of the solution vector.)
+    unsigned int       n_steps      = 0;
     const unsigned int n_time_steps = 200;
-    const double initial_time = 0.;
-    const double final_time = 10.;
+    const double       initial_time = 0.;
+    const double       final_time   = 10.;
 
-    // Next, we solve the diffusion problem using different Runge-Kutta methods.
     std::cout << "Explicit methods:" << std::endl;
     explicit_method (TimeStepping::FORWARD_EULER,
                      n_time_steps,

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.