]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Start a tutorial using the Runge-Kutta methods implemented in TimeStepping.h
authorturcksin <turcksin@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 2 May 2014 20:25:11 +0000 (20:25 +0000)
committerturcksin <turcksin@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 2 May 2014 20:25:11 +0000 (20:25 +0000)
git-svn-id: https://svn.dealii.org/trunk@32878 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-52/CMakeLists.txt [new file with mode: 0644]
deal.II/examples/step-52/doc/builds-on [new file with mode: 0644]
deal.II/examples/step-52/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-52/doc/kind [new file with mode: 0644]
deal.II/examples/step-52/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-52/doc/tooltip [new file with mode: 0644]
deal.II/examples/step-52/step-52.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-52/CMakeLists.txt b/deal.II/examples/step-52/CMakeLists.txt
new file mode 100644 (file)
index 0000000..ec29e87
--- /dev/null
@@ -0,0 +1,31 @@
+##
+#  CMake script for the step-52 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-52")
+
+# Declare all source files the target consists of:
+SET(TARGET_SRC
+  ${TARGET}.cc
+  # You can specify additional files here!
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
+
+FIND_PACKAGE(deal.II 8.0 QUIET
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/deal.II/examples/step-52/doc/builds-on b/deal.II/examples/step-52/doc/builds-on
new file mode 100644 (file)
index 0000000..48a0f73
--- /dev/null
@@ -0,0 +1 @@
+step-4
diff --git a/deal.II/examples/step-52/doc/intro.dox b/deal.II/examples/step-52/doc/intro.dox
new file mode 100644 (file)
index 0000000..4ae88a5
--- /dev/null
@@ -0,0 +1,102 @@
+<br>
+
+<i>This program was contributed by Bruno Turcksin and Damien Lebrun-Grandie.</i>
+
+<a name="Intro"></a>
+<h1>Introducion</h1>
+
+This program shows how to use Runge-Kutta methods to solve a time-dependent
+problem.
+
+<h3>Problem statement<h3>
+
+In this example, we solve the energy-integrated time-dependent diffusion
+approximation of the neutron transport equation (see step-28 for the
+time-independent multigroup diffusion). We assume that the medium is not
+fissible and therefore, the neutron flux satisfies the following equation:
+@f{eqnarray*}
+\frac{1}{v}\frac{\partial \phi(x,t)}{\partial t} = \nabla D(x) \nabla \phi(x,t)
+- \Sigma_a(x) \phi(x,t) + S(x,t)
+@f}
+augmented by appropriate boundary conditions. Here, $v$ is the velocity of
+neutrons, $D$ is the diffusion coefficient, $\Sigma_a$ is the <i>absorption
+cross section</i>, and $S$ is a source. Because we are only interested in the
+time dependence, we assume that $D$ and $\Sigma_a$ are constant. In this
+example, we are only interested in the error in time and thus, we are looking
+for a solution of the form:
+@f{eqnarray*}
+\phi(x,t) = A\sin(\omega t)(bx-x^2).
+@f}
+By using quadratic finite elements, we will not have any spatial error. We
+impose the following boundary conditions: homogeneous Dirichlet fo $x=0$ and
+$x=b$ and homogeneous Neumann conditions for $y=0$ and $y=b$. The source is
+given by:
+@f{eqnarray*}
+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
+\f$\phi\left(x,\frac{\pi}{\omega}\right) = 0$. Therefore, we can easily
+compute the error at this time since it is simply the norm of the solution
+found.
+
+<h3>Runge-Kutta</h3>
+
+The Runke-Kutta methods implemented in deal.II assume that the equation to be
+solved can be written as:
+@f{eqnarray*}
+\frac{dy}{dt} = f(t,y).
+@f}
+When using finite elements, the previous equation becomes:
+@f{eqnarray*}
+M\frac{dy}{dt} = f(t,y),
+@f}
+where $M$ is the mass matrix. Therefore, we have:
+@f{eqnarray*}
+\frac{dy}{dt} = M^{-1}f(t,y).
+@f}
+Runke-Kutta methods can be written as:
+@f{eqnarray*}
+y_{n+1} = y_n + \sum_{i=1}^s b_i k_i
+@f}
+where
+@f{eqnarray*}
+k_i = h M^{-1} f(t_n+c_ih,y_n+\sum_{j=1}^sa_{ij}k_j)
+@f}
+with $a_{ij}$, $b_i$, and $c_i$ are known coefficient and $h$ is the time step
+used. The methods currently 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
+</ol> 
+
+<h4>Explicit Runge-Kutta</h4>
+These methods that include for forward Euler, third order Runge-Kutta, and
+fourth order Runge-Kutta, require a function to evaluate $M^{-1}f(t,y). These
+methods become unstable when the time step chosen is too large.
+
+<h4>Embedded Runge-Kutta</h4>
+These methods include Heun-Euler, Bogacki-Shampine, Dormand-Prince (ode45 in
+Matlab), Fehlberg, and Cash-Karp. These methods use a low order method to
+estimate the error and decide if the time step needs to be refined or it can be
+coarsen. Only embedded explicit methods have been implemented so far.
+
+<h4>Implicit Runge-Kutta</h4>
+These methods include backward Euler, implicit midpoint, Crank-Nicolson, and the
+two stages SDIRK. These methods require to evaluate $M^{-1}f(t,y)$ and
+$\left(I-\Delta t M^{-1} \frac{\partial f}{\partial Y}\right) = \left(M - \Delta
+t \frac{\partial f}{\partial y}\right)^{-1} M$. These methods are always stable.
+
+<h3>Remarks</h3>
+To simplify the problem, we solve the domain in two dimensional and the mesh is
+uniform (there is no need to adapt the mesh since we use quadratic finite
+elements and the exact solution is quadratic). Going from a two dimensional
+domain to a three dimensional domain is not very challenging. However if the
+mesh must be adapted, we cannot forget to:
+<ol>
+<li> project the solution to the new mesh when the mesh is changed. The mesh
+used should be the same at the beginning and at the end of the time step.
+<li> update the mass matrix and its inverse.
+</ol>
diff --git a/deal.II/examples/step-52/doc/kind b/deal.II/examples/step-52/doc/kind
new file mode 100644 (file)
index 0000000..86a44aa
--- /dev/null
@@ -0,0 +1 @@
+time dependent
diff --git a/deal.II/examples/step-52/doc/results.dox b/deal.II/examples/step-52/doc/results.dox
new file mode 100644 (file)
index 0000000..4a70bbe
--- /dev/null
@@ -0,0 +1,30 @@
+<h1>Results</h1>
+
+The output of this program consist of the console output and solutions given in
+vtu format.
+
+The console output is:
+@code
+Forward Euler error: 1.00883
+Third order Runge-Kutta error: 0.000227982
+Fourth order Runge-Kutta error: 1.90541e-06
+Backward Euler error: 1.03428
+Implicit Midpoint error: 0.00862702
+Crank-Nicolson error: 0.00862675
+SDIRK error: 0.0042349
+Heun-Euler error: 0.0073012
+Number of steps done: 284
+Bogacki-Shampine error: 0.000207511
+Number of steps done: 200
+Dopri error: 4.01775e-09
+Number of steps done: 200
+Fehlberg error: 9.89504e-09
+Number of steps done: 200
+Cash-Karp error: 2.5579e-10
+Number of steps done: 200
+@endcode
+
+Like expected the high-order methods give a more accurate solutions. We see that
+the Heun-Euler method adapted the number of time steps in order to satisfy the
+tolerance. The others embedded methods did not need to change the number of time
+steps.
diff --git a/deal.II/examples/step-52/doc/tooltip b/deal.II/examples/step-52/doc/tooltip
new file mode 100644 (file)
index 0000000..d35ea62
--- /dev/null
@@ -0,0 +1 @@
+Time-dependent diffusion equation. Neutron transport
diff --git a/deal.II/examples/step-52/step-52.cc b/deal.II/examples/step-52/step-52.cc
new file mode 100644 (file)
index 0000000..1031e49
--- /dev/null
@@ -0,0 +1,629 @@
+/* ---------------------------------------------------------------------
+ * $Id: step-52.cc 30526 2013-08-29 20:06:27Z felix.gruber $
+ *
+ * Copyright (C) 2014 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Authors: Damien Lebrun-Grandie, Bruno Turcksin, 2014
+ */
+
+// @sect3{Include files}
+
+// The first task as usal is to include the functionality of these well-known
+// deal.II library files and some C++ header files.
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <fstream>
+#include <iostream>
+#include <cmath>
+#include <map>
+
+// This is the only include file that is new: It includes all the Runge-Kutta
+// methods.
+#include <deal.II/base/time_stepping.h>
+
+
+// The next step is like in all previous tutorial programs: We put everything
+// into a namespace of its own and then import the deal.II classes and functions
+// into it.
+namespace Step52
+{
+  using namespace dealii;
+
+  // @sect3{Diffusion}
+
+  // Now, here comes the declaration of the main class.
+  class Diffusion
+  {
+    public:
+      Diffusion();
+
+      void run();
+
+    private:
+      // Create the sparsity_pattern and initialize system_matrix.
+      void setup_system();
+
+      // Assemble the part of the matrix of the system that does not depend on
+      // the time.
+      void assemble_system();
+
+      // Compute the intensity of the source at the given point.
+      double get_source(double time,const Point<2> &point) const;
+      
+      // Evaluate the diffusion equation \f$M^{-1}(f(t,y))\f$
+      Vector<double> evaluate_diffusion(const double time, const Vector<double> &y) const;
+
+      // Evaluate \f$\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1} = 
+      // \left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M \f$
+      Vector<double> id_minus_tau_J_inverse(const double time, const double tau, 
+          const Vector<double> &y);
+      
+      // Output the results as vtu
+      void output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const;
+
+      // Driver for the explicit methods
+      void explicit_method(TimeStepping::runge_kutta_method method,
+                           const unsigned int n_time_steps,
+                           const double       initial_time,
+                           const double       final_time);
+
+      // Driver for the implicit methods
+      void implicit_method(TimeStepping::runge_kutta_method method,
+                           const unsigned int n_time_steps,
+                           const double       initial_time,
+                           const double       final_time);
+
+      // Driver for the embedded explicit methods. Returns the number of steps
+      // executed.
+      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                 fe_degree;
+
+      double                       diffusion_coefficient;
+      double                       absorption_xs;
+      
+      Triangulation<2>             triangulation;
+      
+      FE_Q<2>                      fe;
+      
+      DoFHandler<2>                dof_handler;
+      
+      ConstraintMatrix             constraint_matrix;
+      
+      SparsityPattern              sparsity_pattern;
+      
+      SparseMatrix<double>         system_matrix;
+      SparseMatrix<double>         mass_matrix;
+      SparseMatrix<double>         mass_minus_tau_Jacobian;
+      
+      SparseDirectUMFPACK          inverse_mass_matrix;
+      
+      Vector<double>               solution;
+  };
+
+
+
+  // We choose quadratic finite elements so that there are no spatial error.
+  Diffusion::Diffusion()
+    :
+      fe_degree(2),
+      diffusion_coefficient(1./30.),
+      absorption_xs(1.),
+      fe(fe_degree),
+      dof_handler(triangulation)
+  {}
+
+
+
+  void Diffusion::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    // Create the constraint matrix.
+    VectorTools::interpolate_boundary_values(dof_handler,1,ZeroFunction<2>(),constraint_matrix);
+    constraint_matrix.close();
+
+    // Create the sparsity_pattern.
+    CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler,c_sparsity,constraint_matrix);
+    sparsity_pattern.copy_from(c_sparsity);
+
+    system_matrix.reinit(sparsity_pattern);
+    mass_matrix.reinit(sparsity_pattern);
+    mass_minus_tau_Jacobian.reinit(sparsity_pattern);
+    solution.reinit(dof_handler.n_dofs());
+  }
+
+
+
+  void Diffusion::assemble_system()
+  {
+    system_matrix = 0.;
+    mass_matrix = 0.;
+
+    const QGauss<2> quadrature_formula(fe_degree+1);
+
+    FEValues<2> fe_values(fe, quadrature_formula,
+        update_values | update_gradients | update_JxW_values);
+
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+    FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    typename DoFHandler<2>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+           
+    // Compute \f$-\int D \nabla b \cdot \nabla b - \int \Sigma_a b b\f$ and \f$\int b b\f$
+    for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0.;
+      cell_mass_matrix = 0.;
+
+      fe_values.reinit (cell);
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+        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)) *
+                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);
+
+      constraint_matrix.distribute_local_to_global(cell_matrix,local_dof_indices,system_matrix);
+      constraint_matrix.distribute_local_to_global(cell_mass_matrix,local_dof_indices,mass_matrix);
+    }
+
+    // Compute the inverse of the mass matrix.
+    inverse_mass_matrix.initialize(mass_matrix);
+  }
+
+
+
+  double Diffusion::get_source(double time,const Point<2> &point) const
+  {
+    const double pi = 3.14159265358979323846;
+    const double intensity = 10.;
+    const double frequency = 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;
+  }
+
+
+
+  Vector<double> Diffusion::evaluate_diffusion(const double time, const Vector<double> &y) const
+  {
+    Vector<double> tmp(dof_handler.n_dofs());
+    tmp = 0.;
+    // Compute system_matrix*y
+    system_matrix.vmult(tmp,y);
+
+
+    // Compute the source term
+    const QGauss<2> quadrature_formula(fe_degree+1);
+
+    FEValues<2> fe_values(fe, quadrature_formula,
+        update_values | update_quadrature_points | update_JxW_values);
+
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    Vector<double>  cell_source(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    typename DoFHandler<2>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+           
+    for (; cell!=endc; ++cell)
+    {
+      cell_source = 0.;
+
+      fe_values.reinit (cell);
+
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+      {
+        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) * 
+            fe_values.JxW(q_point);
+      }
+
+      cell->get_dof_indices(local_dof_indices);
+
+      // Add the source term to the tmp vector.
+      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);
+
+    return value;
+  }
+
+
+
+  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);
+    inverse_mass_minus_tau_Jacobian.initialize(mass_minus_tau_Jacobian);
+    mass_matrix.vmult(tmp,y);
+    inverse_mass_minus_tau_Jacobian.vmult(result,tmp);
+
+    return result;
+  }
+
+
+
+  void Diffusion::output_results(unsigned int time_step,TimeStepping::runge_kutta_method method) const
+  {
+    std::string method_name;
+
+    switch (method)
+    {
+      case TimeStepping::FORWARD_EULER :
+        {
+          method_name = "forward_euler";
+          break;
+        }
+      case TimeStepping::RK_THIRD_ORDER :
+        {
+          method_name = "rk3";
+          break;
+        }
+      case TimeStepping::RK_CLASSIC_FOURTH_ORDER :
+        {
+          method_name = "rk4";
+          break;
+        }
+      case TimeStepping::BACKWARD_EULER :
+        {
+          method_name = "backward_euler";
+          break;
+        }
+      case TimeStepping::IMPLICIT_MIDPOINT :
+        {
+          method_name = "implicit_midpoint";
+          break;
+        }
+      case TimeStepping::SDIRK_TWO_STAGES :
+        {
+          method_name = "sdirk";
+          break;
+        }
+      case TimeStepping::HEUN_EULER :
+        {
+          method_name = "heun_euler";
+          break;
+        }
+      case TimeStepping::BOGACKI_SHAMPINE :
+        {
+          method_name = "bocacki_shampine";
+          break;
+        }
+      case TimeStepping::DOPRI :
+        {
+          method_name = "dopri";
+          break;
+        }
+      case TimeStepping::FEHLBERG :
+        {
+          method_name = "fehlberg";
+          break;
+        }
+      case TimeStepping::CASH_KARP :
+        {
+          method_name = "cash_karp";
+          break;
+        }
+      default :
+        {
+          break;
+        }
+    }
+
+    DataOut<2> data_out;
+
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution, "flux");
+
+    data_out.build_patches();
+
+    const std::string filename = "solution-" + method_name + "-" +
+                                 Utilities::int_to_string (time_step, 3) +
+                                 ".vtu";
+    std::ofstream output(filename.c_str());
+    data_out.write_vtu(output);
+  }
+
+
+
+  void Diffusion::explicit_method(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;
+    solution = 0.;
+
+    TimeStepping::ExplicitRungeKutta<Vector<double> > explicit_runge_kutta(method);
+    output_results(0,method);
+    for (unsigned int i=0; i<n_time_steps; ++i)
+    {
+      // Because we use a member function, we need to bind this to the
+      // function.
+      time = explicit_runge_kutta.evolve_one_time_step(
+          std_cxx1x::bind(&Diffusion::evaluate_diffusion,this,std_cxx1x::_1,std_cxx1x::_2),
+          time,time_step,solution);
+
+      // We output the results every 10 time steps.
+      if ((i+1)%10==0)
+        output_results(i+1,method);
+    }
+  }
+
+
+
+  void Diffusion::implicit_method(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;
+    solution = 0.;
+
+    TimeStepping::ImplicitRungeKutta<Vector<double> > implicit_runge_kutta(method);
+    output_results(0,method);
+    for (unsigned int i=0; i<n_time_steps; ++i)
+    {
+      // Because we use a member function, we need to bind this to the
+      // function.
+      time = implicit_runge_kutta.evolve_one_time_step(
+          std_cxx1x::bind(&Diffusion::evaluate_diffusion,this,std_cxx1x::_1,std_cxx1x::_2),
+          std_cxx1x::bind(&Diffusion::id_minus_tau_J_inverse,this,std_cxx1x::_1,std_cxx1x::_2,
+            std_cxx1x::_3),
+          time,time_step,solution);
+
+      // We output the results every 10 time steps.
+      if ((i+1)%10==0)
+        output_results(i+1,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)
+  {
+    double time_step = (final_time-initial_time)/static_cast<double> (n_time_steps);
+    double time = initial_time;
+    // Factor multiplying the current time step when the error is below the
+    // threshold.
+    const double coarsen_param = 1.2;
+    // Factor multiplying the current time step when the error is above the
+    // threshold.
+    const double refine_param = 0.8;
+    // Smallest time step acceptable.
+    const double min_delta = 1e-8;
+    // Largest time step acceptable.
+    const double max_delta = 10*time_step;
+    // Threshold above which the time step is refined.
+    const double refine_tol = 1e-1;
+    // Threshold below which the time step is coarsen.
+    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);
+    unsigned int n_steps=0;
+    while (time<final_time)
+    {
+      // Choose the last time step to reach final_time.
+      if (time+time_step>final_time)
+        time_step = final_time-time;
+
+      // Because we use a member function, we need to bind this to the
+      // function.
+      time = embedded_explicit_runge_kutta.evolve_one_time_step(
+          std_cxx1x::bind(&Diffusion::evaluate_diffusion,this,std_cxx1x::_1,std_cxx1x::_2),
+          time,time_step,solution);
+
+      // We output the results every 10 time steps.
+      if ((n_steps+1)%10==0)
+        output_results(n_steps+1,method);
+
+      // Update the time step
+      time_step = embedded_explicit_runge_kutta.get_status().delta_t_guess;
+      ++n_steps;
+    }
+
+    return n_steps;
+  }
+
+
+
+  void Diffusion::run()
+  {
+    // Create the grid (a square [0,5]x[0,5]) and refine the mesh four times.
+    // The final gird has 16 times 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();
+           
+    for (; cell!=endc; ++cell)
+      for (unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
+        if (cell->face(f)->at_boundary())
+        {
+          if ((cell->face(f)->center()[0]==0.) || (cell->face(f)->center()[0]==5.))
+            cell->face(f)->set_boundary_indicator(1);  
+          else
+            cell->face(f)->set_boundary_indicator(0);  
+        }
+
+    setup_system();
+
+    assemble_system();
+
+    unsigned int n_steps = 0;
+    const unsigned int n_time_steps = 200;
+    const double initial_time = 0.;
+    const double final_time = 10.;
+
+    // Use forward Euler
+    explicit_method(TimeStepping::FORWARD_EULER,n_time_steps,initial_time,final_time);
+    std::cout<<"Forward Euler error: "<<solution.l2_norm()<<std::endl;
+    // Use third order Runge-Kutta
+    explicit_method(TimeStepping::RK_THIRD_ORDER,n_time_steps,initial_time,final_time);
+    std::cout<<"Third order Runge-Kutta error: "<<solution.l2_norm()<<std::endl;
+    // Use fourth order Runge-Kutta
+    explicit_method(TimeStepping::RK_CLASSIC_FOURTH_ORDER,n_time_steps,initial_time,final_time);
+    std::cout<<"Fourth order Runge-Kutta error: "<<solution.l2_norm()<<std::endl;
+
+
+    // Use backward Euler
+    implicit_method(TimeStepping::BACKWARD_EULER,n_time_steps,initial_time,final_time);
+    std::cout<<"Backward Euler error: "<<solution.l2_norm()<<std::endl;
+    // Use implicit midpoint
+    implicit_method(TimeStepping::IMPLICIT_MIDPOINT,n_time_steps,initial_time,final_time);
+    std::cout<<"Implicit Midpoint error: "<<solution.l2_norm()<<std::endl;
+    // Use Crank-NICOLSON
+    implicit_method(TimeStepping::CRANK_NICOLSON,n_time_steps,initial_time,final_time);
+    std::cout<<"Crank-Nicolson error: "<<solution.l2_norm()<<std::endl;
+    // Use two stages SDIRK
+    implicit_method(TimeStepping::SDIRK_TWO_STAGES,n_time_steps,initial_time,final_time);
+    std::cout<<"SDIRK error: "<<solution.l2_norm()<<std::endl;
+
+    
+    // Use Heun-Euler
+    n_steps = embedded_explicit_method(TimeStepping::HEUN_EULER,n_time_steps,initial_time,final_time);
+    std::cout<<"Heun-Euler error: "<<solution.l2_norm()<<std::endl;
+    std::cout<<"Number of steps done: "<<n_steps<<std::endl;
+    // Use Bogacki-Shampine
+    n_steps = embedded_explicit_method(TimeStepping::BOGACKI_SHAMPINE,n_time_steps,initial_time,final_time);
+    std::cout<<"Bogacki-Shampine error: "<<solution.l2_norm()<<std::endl;
+    std::cout<<"Number of steps done: "<<n_steps<<std::endl;
+    // Use Dopri
+    n_steps = embedded_explicit_method(TimeStepping::DOPRI,n_time_steps,initial_time,final_time);
+    std::cout<<"Dopri error: "<<solution.l2_norm()<<std::endl;
+    std::cout<<"Number of steps done: "<<n_steps<<std::endl;
+    // Use Fehlberg
+    n_steps = embedded_explicit_method(TimeStepping::FEHLBERG,n_time_steps,initial_time,final_time);
+    std::cout<<"Fehlberg error: "<<solution.l2_norm()<<std::endl;
+    std::cout<<"Number of steps done: "<<n_steps<<std::endl;
+    // Use Cash-Karp
+    n_steps = embedded_explicit_method(TimeStepping::CASH_KARP,n_time_steps,initial_time,final_time);
+    std::cout<<"Cash-Karp error: "<<solution.l2_norm()<<std::endl;
+    std::cout<<"Number of steps done: "<<n_steps<<std::endl;
+  }
+}
+
+
+
+// The following <code>main</code> function is similar to previous examples as
+// well, and need not be commented on.
+int main ()
+{
+  try
+  {
+    Step52::Diffusion diffusion;
+    diffusion.run();
+  }
+  catch (std::exception &exc)
+  {
+    std::cerr << std::endl << std::endl
+      << "----------------------------------------------------"
+      << std::endl;
+    std::cerr << "Exception on processing: " << std::endl
+      << exc.what() << std::endl
+      << "Aborting!" << std::endl
+      << "----------------------------------------------------"
+      << std::endl;
+    return 1;
+  }
+  catch (...)
+  {
+    std::cerr << std::endl << std::endl
+      << "----------------------------------------------------"
+      << std::endl;
+    std::cerr << "Unknown exception!" << std::endl
+      << "Aborting!" << std::endl
+      << "----------------------------------------------------"
+      << std::endl;
+    return 1;
+  };
+
+  return 0;
+}

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.