]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add Runge-Kutta methods for time integration.
authorturcksin <turcksin@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 18 Feb 2014 04:08:08 +0000 (04:08 +0000)
committerturcksin <turcksin@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 18 Feb 2014 04:08:08 +0000 (04:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@32500 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/base/time_stepping.h [new file with mode: 0644]
deal.II/include/deal.II/base/time_stepping.templates.h [new file with mode: 0644]
deal.II/source/base/CMakeLists.txt
deal.II/source/base/time_stepping.cc [new file with mode: 0644]
deal.II/source/base/time_stepping.inst.in [new file with mode: 0644]
tests/base/time_stepping_01.cc [new file with mode: 0644]
tests/base/time_stepping_01.output [new file with mode: 0644]

index cc1f4f4b89734b2751f8df5f1d08f6059e16343a..beb94ec9ee9a57e8b53b6a8686e66f4605f46cfa 100644 (file)
@@ -135,6 +135,12 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>New: There is a new namespace TimeStepping for the algorithms that do time 
+  integrations. In this new namespace, several Runge-Kutta methods have been 
+  implemented: explicit methods, implicit methods, and embedded explicit methods.
+  <br>
+  (Damien Lebrun-Grandie, Bruno Turcksin, 2014/02/17)
+
   <li>New: There is now a class FEEvaluationDGP that implements matrix-free
   evaluation routines by truncated tensor products for FE_DGP elements.
   <br>
diff --git a/deal.II/include/deal.II/base/time_stepping.h b/deal.II/include/deal.II/base/time_stepping.h
new file mode 100644 (file)
index 0000000..1730cc8
--- /dev/null
@@ -0,0 +1,575 @@
+// ---------------------------------------------------------------------
+// $Id: time_stepping.h 32217 2014-01-15 16:34:36Z bangerth $
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__time_stepping_h
+#define __deal2__time_stepping_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/std_cxx1x/function.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * Namespace containing the time stepping methods.
+ *
+ * @author Bruno Turcksin
+ * @date 2014
+ */
+
+namespace TimeStepping
+{
+  /**
+   * Runge-Kutta methods available:
+   *   - Explicit methods:
+   *     - FORWARD_EULER: first order
+   *     - RK_THIRD_ORDER: third order Runge-Kutta
+   *     - RK_CLASSIC_FOURTH_ORDER: classical fourth order Runge-Kutta
+   *   - Implicit methods:
+   *     - BACKWARD_EULER: first order
+   *     - IMPLICIT_MIDPOINT: second order
+   *     - CRANK_NICHOLSON: second order
+   *     - SDIRK_TWO_STAGES: second order
+   *   - Embedded explicit methods:
+   *     - HEUN_EULER: second order
+   *     - BOGACKI_SHAMPINE: third order
+   *     - DOPRI: Dormand-Prince fifth order (method used by ode45 in MATLAB)
+   *     - FEHLBERG: fifth order
+   *     - CASH_KARP: firth order
+   */
+  enum runge_kutta_method { FORWARD_EULER, RK_THIRD_ORDER, RK_CLASSIC_FOURTH_ORDER,
+                            BACKWARD_EULER, IMPLICIT_MIDPOINT, CRANK_NICOLSON,
+                            SDIRK_TWO_STAGES, HEUN_EULER, BOGACKI_SHAMPINE, DOPRI,
+                            FEHLBERG, CASH_KARP
+                          };
+
+
+
+  /**
+   * Reason for exiting evolve_one_time_step when using an embedded
+   * method: DELTA_T (the time step is in the valid range), MIN_DELTA_T (the
+   * time step was increased to the minimum acceptable time step), MAX_DELTA_T
+   * (the time step was reduced to the maximum acceptable time step).
+   */
+  enum embedded_runge_kutta_time_step { DELTA_T, MIN_DELTA_T, MAX_DELTA_T };
+
+
+
+  /**
+   * Abstract class for time stepping methods. These methods assume that the
+   * equation has the form: \f$ \frac{\partial y}{\partial t} = f(t,y) \f$.
+   */
+  template <typename VECTOR>
+  class TimeStepping
+  {
+  public:
+    /**
+     * Purely virtual function. This function is used to advance from time @p
+     * t to p+ @p delta_t. @p F is a vector of functions \f$ f(t,y) \f$ that should be
+     * integrated, the input parameters are the time t and the vector y and the
+     * output is value of f at this point. @p J_inverse is a vector
+     * functions that compute the inverse of the Jacobians associated to the
+     * implicit problems. The input parameters are the
+     * time, \f$ \tau \f$, and a vector. The output is the value of function
+     * at this point. This function returns the time at the end of the
+     * time step.
+     */
+    virtual double evolve_one_time_step(
+      std::vector<std_cxx1x::function<VECTOR (const double, const VECTOR &)> > &F,
+      std::vector<std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> > & J_inverse,
+      double t,
+      double delta_t,
+      VECTOR &y) = 0;
+
+    /**
+     * Empty structure used to store informations.
+     */
+    struct Status {};
+
+    /**
+     * Purely virtual function that return Status.
+     */
+    virtual const Status &get_status() const = 0;
+  };
+
+
+
+  /**
+   * Base class for the Runge-Kutta method
+   *
+   * @author Damien Lebrun-Grandie, Bruno Turcksin
+   * @date 2014
+   */
+  template <typename VECTOR>
+  class RungeKutta : public TimeStepping<VECTOR>
+  {
+  public:
+    /**
+     * Purely virtual method used to initialize the Runge-Kutta method.
+     */
+    virtual void initialize(runge_kutta_method method) = 0;
+    /**
+     * This function is used to advance from time @p
+     * t to p+ @p delta_t. @p F is a vector of functions \f$ f(t,y) \f$ that should be
+     * integrated, the input parameters are the time t and the vector y and the
+     * output is value of f at this point. @p J_inverse is a vector
+     * functions that compute the inverse of the Jacobians associated to the
+     * implicit problems. The input parameters are the
+     * time, \f$ \tau \f$, and a vector. The output is the value of function
+     * at this point. This function returns the time at the end of the
+     * time step. When using Runge-Kutta methods, @p F and @ J_inverse can
+     * only contain one element.
+     */
+    virtual double evolve_one_time_step(
+      std::vector<std_cxx1x::function<VECTOR (const double, const VECTOR &)> > &F,
+      std::vector<std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> > & J_inverse,
+      double t,
+      double delta_t,
+      VECTOR &y);
+
+    /**
+     * Purely virtual function. This function is used to advance from time @p t
+     * to p+ @p delta_t. @p f  is the function \f$ f(t,y) \f$ that should be
+     * integrated, the input parameters are the time t and the vector y and the
+     * output is value of f at this point. @p id_minus_tau_J_inverse is a function
+     * that computes \f$ inv(I-\tau J)\f$ where \f$ I \f$ is the identity matrix,
+     * \f$ \tau \f$ is given, and \f$ J \f$ is the Jacobian \f$ \frac{\partial
+     * J}{\partial y} \f$. The input parameters are the time, \f$ \tau \f$, and
+     * a vector. The output is the value of function at this point.
+     * evolve_one_time_step returns the time at the end of the time step.
+     */
+    virtual double evolve_one_time_step(
+      std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+      std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> id_minus_tau_J_inverse,
+      double t,
+      double delta_t,
+      VECTOR &y) = 0;
+
+  protected:
+    /**
+     * Number of stages of the Runge-Kutta method.
+     */
+    unsigned int n_stages;
+
+    /**
+     * Butcher tableau coefficients.
+     */
+    std::vector<double> b;
+
+    /**
+     * Butcher tableau coefficients.
+     */
+    std::vector<double> c;
+
+    /**
+     * Butcher tableau coefficients.
+     */
+    std::vector<std::vector<double> > a;
+  };
+
+
+
+  /**
+   * ExplicitRungeKutta is derived from RungeKutta and implement the explicit methods.
+   */
+  template <typename VECTOR>
+  class ExplicitRungeKutta : public RungeKutta<VECTOR>
+  {
+  public:
+    /**
+     * Default constructor. initialize(runge_kutta_method) needs to be called
+     * before the object can be used.
+     */
+    ExplicitRungeKutta() {};
+
+    /**
+     * Constructor. This function calls initialize(runge_kutta_method).
+     */
+    ExplicitRungeKutta(runge_kutta_method method);
+
+    /**
+     * Initialize the explicit Runge-Kutta method.
+     */
+    void initialize(runge_kutta_method method);
+
+    /**
+     * This function is used to advance from time @p t to p+ @p delta_t. @p f
+     * is the function \f$ f(t,y) \f$ that should be integrated, the input
+     * parameters are the time t and the vector y and the output is value of
+     * f at this point. @p id_minus_tau_J_inverse is a function that computes
+     * \f$ inv(I-\tau J)\f$ where \f$ I \f$ is the identity matrix, \f$ \tau
+     * \f$ is given, and \f$ J \f$ is the Jacobian \f$ \frac{\partial
+     * J}{\partial y} \f$. The input parameter are the time, \f$ \tau \f$, and
+     * a vector. The output is the value of function at this point.
+     * evolve_one_time_step returns the time at the end of the time step.
+     */
+    double evolve_one_time_step(
+      std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+      std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> id_minus_tau_J_inverse,
+      double t,
+      double delta_t,
+      VECTOR &y);
+
+    /**
+     * This function is used to advance from time @p t to p+ @p delta_t.
+     * This function is similar to the one derived from RungeKutta, but
+     * does not required id_minus_tau_J_inverse because it is not used for
+     * explicit methods. evolve_one_time_step returns the time at the end of the
+     * time step.
+     */
+    double evolve_one_time_step(std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+                                double t,
+                                double delta_t,
+                                VECTOR &y);
+
+    /**
+     * This structure stores the name of the method used.
+     */
+    struct Status : public TimeStepping<VECTOR>::Status
+    {
+      runge_kutta_method method;
+    };
+
+    /**
+     * Return the status of the current object.
+     */
+    const Status &get_status() const;
+
+  private:
+    /**
+     * Compute the different stages needed.
+     */
+    void compute_stages(std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+                        const double t,
+                        const double delta_t,
+                        const VECTOR &y,
+                        std::vector<VECTOR> &f_stages) const;
+
+    /**
+     * Status structure of the object.
+     */
+    Status status;
+  };
+
+
+
+  /**
+   * This class is derived from RungeKutta and implement the implicit
+   * methods. This class works only for Diagonal Implicit Runge-Kutta
+   * (DIRK) methods.
+   */
+  template <typename VECTOR>
+  class ImplicitRungeKutta : public RungeKutta<VECTOR>
+  {
+  public:
+    /**
+     * Default constructor. initialize(runge_kutta_method) and
+     * set_newton_solver_parameters(unsigned int,double)
+     * need to be called before the object can be used.
+     */
+    ImplicitRungeKutta() {};
+
+    /**
+     * Constructor. This function calls initialize(runge_kutta_method) and
+     * initialize the maximum number of iterations and the tolerance of the
+     * Newton solver.
+     */
+    ImplicitRungeKutta(runge_kutta_method method, unsigned int max_it=100, double tolerance=1e-6);
+
+    /**
+     * Initialize the implicit Runge-Kutta method.
+     */
+    void initialize(runge_kutta_method method);
+
+    /**
+     * This function is used to advance from time @p t to p+ @p delta_t. @p f
+     * is the function \f$ f(t,y) \f$ that should be integrated, the input
+     * parameters are the time t and the vector y and the output is value of
+     * f at this point. @p id_minus_tau_J_inverse is a function that computes
+     * \f$ inv(I-\tau J)\f$ where \f$ I \f$ is the identity matrix, \f$ \tau
+     * \f$ is given, and \f$ J \f$ is the Jacobian \f$ \frac{\partial
+     * J}{\partial y} \f$. The input parameters are the time, \f$ \tau \f$, and
+     * a vector. The output is the value of function at this point.
+     * evolve_one_time_step returns the time at the end of the time step.
+     */
+    double evolve_one_time_step(
+      std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+      std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> id_minus_tau_J_inverse,
+      double t,
+      double delta_t,
+      VECTOR &y);
+
+    /**
+     * Set the maximum number of iterations and the tolerance used by the
+     * Newton solver.
+     */
+    void set_newton_solver_parameters(unsigned int max_it, double tolerance);
+
+    /**
+     * Structure that stores the name of the method, the number of
+     * Newton iterations and the norm of the residual when exiting the
+     * Newton solver.
+     */
+    struct Status : public TimeStepping<VECTOR>::Status
+    {
+      runge_kutta_method method;
+      unsigned int       n_iterations;
+      double             norm_residual;
+    };
+
+    /**
+     * Return the status of the current object.
+     */
+    const Status &get_status() const;
+
+  private:
+    /**
+     * Compute the different stages needed.
+     */
+    void compute_stages(
+      std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+      std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> id_minus_tau_J_inverse,
+      double t,
+      double delta_t,
+      VECTOR &y,
+      std::vector<VECTOR> &f_stages);
+
+    /**
+     * Newton solver used for the implicit stages.
+     */
+    void newton_solve(std_cxx1x::function<void (const VECTOR &,VECTOR &)> get_residual,
+                      std_cxx1x::function<VECTOR (const VECTOR &)> id_minus_tau_J_inverse,
+                      VECTOR &y);
+
+    /**
+     * Compute the residual needed by the Newton solver.
+     */
+    void compute_residual(std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+                          double t,
+                          double delta_t,
+                          const VECTOR &old_y,
+                          const VECTOR &y,
+                          VECTOR &tendency,
+                          VECTOR &residual) const;
+
+    /**
+     * When using SDIRK, there is no need to compute the linear combination
+     * of the stages. Thus, when this flag is true, the linear combination
+     * is skipped.
+     */
+    bool skip_linear_combi;
+
+    /**
+     * Maximum number of iterations of the Newton solver.
+     */
+    unsigned int max_it;
+
+    /**
+     * Tolerance of the Newton solver.
+     */
+    double tolerance;
+
+    /**
+     * Status structure of the object.
+     */
+    Status status;
+  };
+
+
+
+  /**
+   * This is class is derived from RungeKutta and implement embedded explicit
+   * methods.
+   */
+  template <typename VECTOR>
+  class EmbeddedExplicitRungeKutta : public RungeKutta<VECTOR>
+  {
+  public:
+    /**
+     * Default constructor. initialize(runge_kutta_method) and
+     * set_time_adaptation_parameters(double, double, double, double, double, double)
+     * need to be called before the object can be used.
+     */
+    EmbeddedExplicitRungeKutta() {};
+
+    /**
+     * Constructor. This function calls initialize(runge_kutta_method) and
+     * initialize the parameters needed for time adaptation.
+     */
+    EmbeddedExplicitRungeKutta(runge_kutta_method method,
+                               double coarsen_param = 1.2,
+                               double refine_param = 0.8,
+                               double min_delta = 1e-14,
+                               double max_delta = 1e100,
+                               double refine_tol = 1e-8,
+                               double coarsen_tol = 1e-12);
+
+    /**
+     * Destructor.
+     */
+    ~EmbeddedExplicitRungeKutta()
+    {
+      free_memory();
+    }
+
+    /**
+     * If necessary, deallocate memory allocated by the object.
+     */
+    void free_memory();
+
+    /**
+     * Initialize the embedded explicit Runge-Kutta method.
+     */
+    void initialize(runge_kutta_method method);
+
+    /**
+     * This function is used to advance from time @p t to p+ @p delta_t. @p f
+     * is the function \f$ f(t,y) \f$ that should be integrated, the input
+     * parameters are the time t and the vector y and the output is value of
+     * f at this point. @p id_minus_tau_J_inverse is a function that computes
+     * \f$ inv(I-\tau J)\f$ where \f$ I \f$ is the identity matrix, \f$ \tau
+     * \f$ is given, and \f$ J \f$ is the Jacobian \f$ \frac{\partial
+     * J}{\partial y} \f$. The input parameters are the time, \f$ \tau \f$, and
+     * a vector. The output is the value of function at this point.
+     * evolve_one_time_step returns the time at the end of the time step.
+     */
+    double evolve_one_time_step(
+      std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+      std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> id_minus_tau_J_inverse,
+      double t,
+      double delta_t,
+      VECTOR &y);
+
+    /**
+     * This function is used to advance from time @p t to p+ @p delta_t.
+     * This function is similar to the one derived from TimeStepping, but
+     * does not required id_minus_tau_J_inverse because it is not used for
+     * explicit methods. evolve_one_time_step returns the time at the end of the
+     * time step.
+     */
+    double evolve_one_time_step(std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+                                double t,
+                                double delta_t,
+                                VECTOR &y);
+
+    /**
+     * Set the parameters necessary for the time adaptation.
+     */
+    void set_time_adaptation_parameters(double coarsen_param,
+                                        double refine_param,
+                                        double min_delta,
+                                        double max_delta,
+                                        double refine_tol,
+                                        double coarsen_tol);
+
+    /**
+     * Structure that stores the name of the method, the reason to exit
+     * evolve_one_time_step, the number of iteration inside n_iterations, a guess
+     * of what the next time step should be, and an estimate of the norm of
+     * the error.
+     */
+    struct Status : public TimeStepping<VECTOR>::Status
+    {
+      runge_kutta_method method;
+      embedded_runge_kutta_time_step exit_delta_t;
+      unsigned int n_iterations;
+      double delta_t_guess;
+      double error_norm;
+    };
+
+    /**
+     * Return the status of the current object.
+     */
+    const Status &get_status() const;
+
+  private:
+    /**
+     * Compute the different stages needed.
+     */
+    void compute_stages(std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+                        const double t,
+                        const double delta_t,
+                        const VECTOR &y,
+                        std::vector<VECTOR> &f_stages);
+
+    /**
+     * This parameter is the factor (>1) by which the time step is
+     * multiplied when the time stepping can be coarsen.
+     */
+    double coarsen_param;
+
+    /**
+     * This parameter is the factor (<1) by which the time step is
+     * multiplied when the time stepping must be refined.
+     */
+    double refine_param;
+
+    /**
+     * Smallest time step allowed.
+     */
+    double min_delta_t;
+
+    /**
+     * Largest time step allowed.
+     */
+    double max_delta_t;
+
+    /**
+     * Refinement tolerance: if the error estimate is larger than
+     * refine_tol, the time step is refined.
+     */
+    double refine_tol;
+
+    /**
+     * Coarsening tolerance: if the error estimate is smaller than
+     * coarse_tol, the time step is coarsen.
+     */
+    double coarsen_tol;
+
+    /**
+     * If the flag is true, the last stage is the same as the first stage
+     * and one evaluation of f can be saved.
+     */
+    bool last_same_as_first;
+
+    /**
+     * Butcher tableau coefficients.
+     */
+    std::vector<double> b1;
+
+    /**
+     * Butcher tableau coefficients.
+     */
+    std::vector<double> b2;
+
+    /**
+     * If the last_same_as_first flag is set to true, the last stage is
+     * saved and reused as the first stage of the next time step.
+     */
+    VECTOR *last_stage;
+
+    /**
+     * Status structure of the object.
+     */
+    Status status;
+  };
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/include/deal.II/base/time_stepping.templates.h b/deal.II/include/deal.II/base/time_stepping.templates.h
new file mode 100644 (file)
index 0000000..6a2c98a
--- /dev/null
@@ -0,0 +1,836 @@
+// ---------------------------------------------------------------------
+// $Id: time_stepping.h 32217 2014-01-15 16:34:36Z bangerth $
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__time_stepping_templates_h
+#define __deal2__time_stepping_templates_h
+
+#include <deal.II/base/std_cxx1x/bind.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/time_stepping.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TimeStepping
+{
+  // ----------------------------------------------------------------------
+  // RungeKutta
+  // ----------------------------------------------------------------------
+
+  template <typename VECTOR>
+  double RungeKutta<VECTOR>::evolve_one_time_step(
+    std::vector<std_cxx1x::function<VECTOR (const double, const VECTOR &)> > &F,
+    std::vector<std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> > & J_inverse,
+    double t,
+    double delta_t,
+    VECTOR &y)
+  {
+    AssertThrow(F.size()==0,
+                ExcMessage("RungeKutta methods cannot handle more that one function to integate."));
+    AssertThrow(J_inverse.size()==0,
+                ExcMessage("RungeKutta methods cannot handle more that one function to integate."));
+
+    return evolve_one_time_step(F[0],J_inverse[0],t,delta_t,y);
+  }
+
+
+
+  // ----------------------------------------------------------------------
+  // ExplicitRungeKutta
+  // ----------------------------------------------------------------------
+
+  template <typename VECTOR>
+  ExplicitRungeKutta<VECTOR>::ExplicitRungeKutta(runge_kutta_method method)
+  {
+    initialize(method);
+  }
+
+
+
+  template <typename VECTOR>
+  void ExplicitRungeKutta<VECTOR>::initialize(runge_kutta_method method)
+  {
+    status.method = method;
+
+    switch (method)
+      {
+      case (FORWARD_EULER) :
+      {
+        this->n_stages = 1;
+        this->b.push_back(1.0);
+        this->c.push_back(0.0);
+
+        break;
+      }
+      case (RK_THIRD_ORDER) :
+      {
+        this->n_stages = 3;
+        this->b.reserve(this->n_stages);
+        this->c.reserve(this->n_stages);
+        this->b.push_back(1.0/6.0);
+        this->b.push_back(2.0/3.0);
+        this->b.push_back(1.0/6.0);
+        this->c.push_back(0.0);
+        this->c.push_back(0.5);
+        this->c.push_back(1.0);
+        std::vector<double> tmp;
+        this->a.push_back(tmp);
+        tmp.resize(1);
+        tmp[0] = 0.5;
+        this->a.push_back(tmp);
+        tmp.resize(2);
+        tmp[0] = -1.0;
+        tmp[1] = 2.0;
+        this->a.push_back(tmp);
+
+        break;
+      }
+      case (RK_CLASSIC_FOURTH_ORDER) :
+      {
+        this->n_stages = 4;
+        this->b.reserve(this->n_stages);
+        this->c.reserve(this->n_stages);
+        std::vector<double> tmp;
+        this->a.push_back(tmp);
+        tmp.resize(1);
+        tmp[0] = 0.5;
+        this->a.push_back(tmp);
+        tmp.resize(2);
+        tmp[0] = 0.0;
+        tmp[1] = 0.5;
+        this->a.push_back(tmp);
+        tmp.resize(3);
+        tmp[1] = 0.0;
+        tmp[2] = 1.0;
+        this->a.push_back(tmp);
+        this->b.push_back(1.0/6.0);
+        this->b.push_back(1.0/3.0);
+        this->b.push_back(1.0/3.0);
+        this->b.push_back(1.0/6.0);
+        this->c.push_back(0.0);
+        this->c.push_back(0.5);
+        this->c.push_back(0.5);
+        this->c.push_back(1.0);
+
+        break;
+      }
+      default :
+      {
+        AssertThrow(false,ExcMessage("Unimplemented explicit Runge-Kutta method."));
+      }
+      }
+  }
+
+
+
+  template <typename VECTOR>
+  double ExplicitRungeKutta<VECTOR>::evolve_one_time_step(
+    std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+    std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> id_minus_tau_J_inverse,
+    double t,
+    double delta_t,
+    VECTOR &y)
+  {
+    return evolve_one_time_step(f,t,delta_t,y);
+  }
+
+
+
+  template <typename VECTOR>
+  double ExplicitRungeKutta<VECTOR>::evolve_one_time_step(
+    std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+    double t,
+    double delta_t,
+    VECTOR &y)
+  {
+    std::vector<VECTOR> f_stages(this->n_stages,y);
+    // Compute the different stages needed.
+    compute_stages(f,t,delta_t,y,f_stages);
+
+    // Linear combinations of the stages.
+    for (unsigned int i=0; i<this->n_stages; ++i)
+      y.sadd(1.,delta_t*this->b[i],f_stages[i]);
+
+    return (t+delta_t);
+  }
+
+
+
+  template <typename VECTOR>
+  const typename ExplicitRungeKutta<VECTOR>::Status &ExplicitRungeKutta<VECTOR>::get_status() const
+  {
+    return status;
+  }
+
+
+
+  template <typename VECTOR>
+  void ExplicitRungeKutta<VECTOR>::compute_stages(
+    std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+    const double t,
+    const double delta_t,
+    const VECTOR &y,
+    std::vector<VECTOR> &f_stages) const
+  {
+    for (unsigned int i=0; i<this->n_stages; ++i)
+      {
+        VECTOR Y(y);
+        for (unsigned int j=0; j<i; ++j)
+          Y.sadd(1.,delta_t*this->a[i][j],f_stages[j]);
+        // Evaluate the function f at the point (t+c[i]*delta_t,Y).
+        f_stages[i] = f(t+this->c[i]*delta_t,Y);
+      }
+  }
+
+
+
+  // ----------------------------------------------------------------------
+  // ImplicitRungeKutta
+  // ----------------------------------------------------------------------
+
+  template <typename VECTOR>
+  ImplicitRungeKutta<VECTOR>::ImplicitRungeKutta(runge_kutta_method method,
+                                                 unsigned int max_it,
+                                                 double tolerance)
+    :
+    RungeKutta<VECTOR> (),
+    skip_linear_combi(false),
+    max_it(max_it),
+    tolerance(tolerance)
+  {
+    initialize(method);
+  }
+
+
+
+  template <typename VECTOR>
+  void ImplicitRungeKutta<VECTOR>::initialize(runge_kutta_method method)
+  {
+    status.method = method;
+
+    switch (method)
+      {
+      case (BACKWARD_EULER) :
+      {
+        this->n_stages = 1;
+        this->a.push_back(std::vector<double>(1, 1.0));
+        this->b.push_back(1.0);
+        this->c.push_back(1.0);
+
+        break;
+      }
+      case (IMPLICIT_MIDPOINT) :
+      {
+        this->a.push_back(std::vector<double>(1, 0.5));
+        this->b.push_back(1.0);
+        this->c.push_back(0.5);
+        this->n_stages = 1;
+
+        break;
+      }
+      case (CRANK_NICOLSON) :
+      {
+        this->n_stages = 2;
+        this->b.reserve(this->n_stages);
+        this->c.reserve(this->n_stages);
+        this->a.push_back(std::vector<double>(1, 0.0));
+        this->a.push_back(std::vector<double>(2, 0.5));
+        this->b.push_back(0.5);
+        this->b.push_back(0.5);
+        this->c.push_back(0.0);
+        this->c.push_back(1.0);
+
+        break;
+      }
+      case (SDIRK_TWO_STAGES) :
+      {
+        this->n_stages = 2;
+        this->b.reserve(this->n_stages);
+        this->c.reserve(this->n_stages);
+        double const gamma = 1.0 - 1.0 / std::sqrt(2.0);
+        this->b.push_back(1.0 - gamma);
+        this->b.push_back(gamma);
+        this->a.push_back(std::vector<double>(1, gamma));
+        this->a.push_back(this->b);
+        this->c.push_back(gamma);
+        this->c.push_back(1.0);
+
+        break;
+      }
+      default :
+      {
+        AssertThrow(false,ExcMessage("Unimplemented implicit Runge-Kutta method."));
+      }
+      }
+  }
+
+
+
+  template <typename VECTOR>
+  double ImplicitRungeKutta<VECTOR>::evolve_one_time_step(
+    std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+    std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> id_minus_tau_J_inverse,
+    double t,
+    double delta_t,
+    VECTOR &y)
+  {
+    VECTOR old_y(y);
+    std::vector<VECTOR> f_stages(this->n_stages,y);
+    // Compute the different stages needed.
+    compute_stages(f,id_minus_tau_J_inverse,t,delta_t,y,f_stages);
+
+    // If necessary, compute the linear combinations of the stages.
+    if (skip_linear_combi==false)
+      {
+        y = old_y;
+        for (unsigned int i=0; i<this->n_stages; ++i)
+          y.sadd(1.,delta_t*this->b[i],f_stages[i]);
+      }
+
+    return (t+delta_t);
+  }
+
+
+
+  template <typename VECTOR>
+  void ImplicitRungeKutta<VECTOR>::set_newton_solver_parameters(unsigned int max_it_, double tolerance_)
+  {
+    max_it = max_it_;
+    tolerance = tolerance_;
+  }
+
+
+
+  template <typename VECTOR>
+  const typename ImplicitRungeKutta<VECTOR>::Status &ImplicitRungeKutta<VECTOR>::get_status() const
+  {
+    return status;
+  }
+
+
+
+  template <typename VECTOR>
+  void ImplicitRungeKutta<VECTOR>::compute_stages(
+    std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+    std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> id_minus_tau_J_inverse,
+    double t,
+    double delta_t,
+    VECTOR &y,
+    std::vector<VECTOR> &f_stages)
+  {
+    VECTOR z(y);
+    for (unsigned int i=0; i<this->n_stages; ++i)
+      {
+        VECTOR old_y(z);
+        for (unsigned int j=0; j<i; ++j)
+          old_y.sadd(1.,delta_t*this->a[i][j],f_stages[j]);
+
+        // Solve the nonlinear system using Newton's method
+        const double new_t = t+this->c[i]*delta_t;
+        const double new_delta_t = this->a[i][i]*delta_t;
+        newton_solve(std_cxx1x::bind(&ImplicitRungeKutta<VECTOR>::compute_residual,this,f,new_t,new_delta_t,
+                                     std_cxx1x::cref(old_y),std_cxx1x::_1,std::ref(f_stages[i]),std_cxx1x::_2),
+                     std_cxx1x::bind(id_minus_tau_J_inverse,new_t,new_delta_t,std_cxx1x::_1),y);
+      }
+  }
+
+
+
+  template <typename VECTOR>
+  void ImplicitRungeKutta<VECTOR>::newton_solve(
+    std_cxx1x::function<void (const VECTOR &,VECTOR &)> get_residual,
+    std_cxx1x::function<VECTOR (const VECTOR &)> id_minus_tau_J_inverse,
+    VECTOR &y)
+  {
+    VECTOR residual(y);
+    get_residual(y,residual);
+    unsigned int i=0;
+    const double initial_residual_norm = residual.l2_norm();
+    double norm_residual = initial_residual_norm;
+    while (i<max_it)
+      {
+        y.sadd(1.0,-1.0,id_minus_tau_J_inverse(residual));
+        get_residual(y,residual);
+        norm_residual = residual.l2_norm();
+        if (norm_residual < tolerance)
+          break;
+        ++i;
+      }
+    status.n_iterations = i+1;
+    status.norm_residual = norm_residual;
+  }
+
+
+
+  template <typename VECTOR>
+  void ImplicitRungeKutta<VECTOR>::compute_residual(
+    std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+    double t,
+    double delta_t,
+    const VECTOR &old_y,
+    const VECTOR &y,
+    VECTOR &tendency,
+    VECTOR &residual) const
+  {
+    // The tendency is stored to save one evaluation of f.
+    tendency = f(t,y);
+    residual = tendency;
+    residual.sadd(delta_t,1.0,old_y);
+    residual.sadd(1.0,-1.,y);
+  }
+
+
+
+  // ----------------------------------------------------------------------
+  // EmbeddedExplicitRungeKutta
+  // ----------------------------------------------------------------------
+
+  template <typename VECTOR>
+  EmbeddedExplicitRungeKutta<VECTOR>::EmbeddedExplicitRungeKutta(runge_kutta_method method,
+      double coarsen_param,
+      double refine_param,
+      double min_delta,
+      double max_delta,
+      double refine_tol,
+      double coarsen_tol)
+    :
+    coarsen_param(coarsen_param),
+    refine_param(refine_param),
+    min_delta_t(min_delta),
+    max_delta_t(max_delta),
+    refine_tol(refine_tol),
+    coarsen_tol(coarsen_tol),
+    last_same_as_first(false),
+    last_stage(NULL)
+  {
+    initialize(method);
+  }
+
+
+
+  template <typename VECTOR>
+  void EmbeddedExplicitRungeKutta<VECTOR>::initialize(runge_kutta_method method)
+  {
+    status.method = method;
+
+    switch (method)
+      {
+      case (HEUN_EULER) :
+      {
+        this->n_stages = 2;
+        this->a.push_back(std::vector<double>());
+        this->a.push_back(std::vector<double>(1, 1.0));
+        this->c.push_back(0.0);
+        this->c.push_back(1.0);
+        b1.push_back(0.5);
+        b1.push_back(0.5);
+        b2.push_back(1.0);
+        b2.push_back(0.0);
+
+        break;
+      }
+      case (BOGACKI_SHAMPINE) :
+      {
+        last_same_as_first = true;
+        this->n_stages = 4;
+        this->c.reserve(this->n_stages);
+        this->b1.reserve(this->n_stages);
+        this->b2.reserve(this->n_stages);
+        std::vector<double> tmp;
+        this->a.push_back(tmp);
+        tmp.resize(1);
+        tmp[0] = 0.5;
+        this->a.push_back(tmp);
+        tmp.resize(2);
+        tmp[0] = 0.0;
+        tmp[1] = 0.75;
+        this->a.push_back(tmp);
+        tmp.resize(3);
+        tmp[0] = 2.0/9.0;
+        tmp[1] = 1.0/3.0;
+        tmp[2] = 4.0/9.0;
+        this->a.push_back(tmp);
+        this->c.push_back(0.0);
+        this->c.push_back(0.5);
+        this->c.push_back(0.75);
+        this->c.push_back(1.0);
+        this->b1.push_back(2.0/9.0);
+        this->b1.push_back(1.0/3.0);
+        this->b1.push_back(4.0/9.0);
+        this->b1.push_back(0.0);
+        this->b2.push_back(7.0/24.0);
+        this->b2.push_back(0.25);
+        this->b2.push_back(1.0/3.0);
+        this->b2.push_back(0.125);
+
+        break;
+      }
+      case (DOPRI) :
+      {
+        last_same_as_first = true;
+        this->n_stages = 7;
+        this->c.reserve(this->n_stages);
+        this->b1.reserve(this->n_stages);
+        this->b2.reserve(this->n_stages);
+        std::vector<double> tmp;
+        this->a.push_back(tmp);
+        tmp.resize(1);
+        tmp[0] = 1./5.;
+        this->a.push_back(tmp);
+        tmp.resize(2);
+        tmp[0] = 3./40.;
+        tmp[1] = 9./40.;
+        this->a.push_back(tmp);
+        tmp.resize(3);
+        tmp[0] = 44./45.;
+        tmp[1] = -56./15.;
+        tmp[2] = 32./9.;
+        this->a.push_back(tmp);
+        tmp.resize(4);
+        tmp[0] = 19372./6561.;
+        tmp[1] = -25360./2187.;
+        tmp[2] = 64448./6561.;
+        tmp[3] = -212./729.;
+        this->a.push_back(tmp);
+        tmp.resize(5);
+        tmp[0] = -9017./3168.;
+        tmp[1] = -355./33.;
+        tmp[2] = 46732./5247.;
+        tmp[3] = 49./176.;
+        tmp[4] = -5103./18656;
+        this->a.push_back(tmp);
+        tmp.resize(6);
+        tmp[0] = 35./384.;
+        tmp[1] = 0.;
+        tmp[2] = 500./1113.;
+        tmp[3] = 125./192.;
+        tmp[4] = -2187./6784.;
+        tmp[5] = 11./84.;
+        this->a.push_back(tmp);
+        this->c.push_back(0.);
+        this->c.push_back(1./5.);
+        this->c.push_back(3./10.);
+        this->c.push_back(4./5.);
+        this->c.push_back(8./9.);
+        this->c.push_back(1.);
+        this->c.push_back(1.);
+        this->b1.push_back(35./384.);
+        this->b1.push_back(0.);
+        this->b1.push_back(500./1113.);
+        this->b1.push_back(125./192.);
+        this->b1.push_back(-2187./6784.);
+        this->b1.push_back(11./84.);
+        this->b1.push_back(0.);
+        this->b2.push_back(5179./57600.);
+        this->b2.push_back(0.);
+        this->b2.push_back(7571./16695.);
+        this->b2.push_back(393./640.);
+        this->b2.push_back(-92097./339200.);
+        this->b2.push_back(187./2100.);
+        this->b2.push_back(1./40.);
+
+        break;
+      }
+      case (FEHLBERG) :
+      {
+        this->n_stages = 6;
+        this->c.reserve(this->n_stages);
+        this->b1.reserve(this->n_stages);
+        this->b2.reserve(this->n_stages);
+        std::vector<double> tmp;
+        this->a.push_back(tmp);
+        tmp.resize(1);
+        tmp[0] = 0.25;
+        this->a.push_back(tmp);
+        tmp.resize(2);
+        tmp[0] = 0.09375;
+        tmp[1] = 0.28125;
+        this->a.push_back(tmp);
+        tmp.resize(3);
+        tmp[0] = 1932.0/2197.0;
+        tmp[1] = -7200.0/2197.0;
+        tmp[2] = 7296.0/2197.0;
+        this->a.push_back(tmp);
+        tmp.resize(4);
+        tmp[0] = 439.0/216.0;
+        tmp[1] = -8.0;
+        tmp[2] = 3680.0/513.0;
+        tmp[3] = -845.0/4104.0;
+        this->a.push_back(tmp);
+        tmp.resize(5);
+        tmp[0] = -8.0/27.0;
+        tmp[1] = 2.0;
+        tmp[2] = -3544.0/2565.0;
+        tmp[3] = 1859.0/4104.0;
+        tmp[4] = -0.275;
+        this->a.push_back(tmp);
+        this->c.push_back(0.0);
+        this->c.push_back(0.25);
+        this->c.push_back(0.375);
+        this->c.push_back(12.0/13.0);
+        this->c.push_back(1.0);
+        this->c.push_back(0.5);
+        this->b1.push_back(16.0/135.0);
+        this->b1.push_back(0.0);
+        this->b1.push_back(6656.0/12825.0);
+        this->b1.push_back(28561.0/56430.0);
+        this->b1.push_back(-0.18);
+        this->b1.push_back(2.0/55.0);
+        this->b2.push_back(25.0/216.0);
+        this->b2.push_back(0.0);
+        this->b2.push_back(1408.0/2565.0);
+        this->b2.push_back(2197.0/4104.0);
+        this->b2.push_back(-0.2);
+        this->b2.push_back(0.0);
+
+        break;
+      }
+      case (CASH_KARP) :
+      {
+        this->n_stages = 6;
+        this->c.reserve(this->n_stages);
+        this->b1.reserve(this->n_stages);
+        this->b2.reserve(this->n_stages);
+        std::vector<double> tmp;
+        this->a.push_back(tmp);
+        tmp.resize(1);
+        tmp[0] = 0.2;
+        this->a.push_back(tmp);
+        tmp.resize(2);
+        tmp[0] = 0.075;
+        tmp[1] = 0.225;
+        this->a.push_back(tmp);
+        tmp.resize(3);
+        tmp[0] = 0.3;
+        tmp[1] = -0.9;
+        tmp[2] = 1.2;
+        this->a.push_back(tmp);
+        tmp.resize(4);
+        tmp[0] = -11.0/54.0;
+        tmp[1] = 2.5;
+        tmp[2] = -70.0/27.0;
+        tmp[3] = 35.0/27.0;
+        this->a.push_back(tmp);
+        tmp.resize(5);
+        tmp[0] = 1631.0/55296.0;
+        tmp[1] = 175.0/512.0;
+        tmp[2] = 575.0/13824.0;
+        tmp[3] = 44275.0/110592.0;
+        tmp[4] = 253.0/4096.0;
+        this->a.push_back(tmp);
+        this->c.push_back(0.0);
+        this->c.push_back(0.2);
+        this->c.push_back(0.3);
+        this->c.push_back(0.6);
+        this->c.push_back(1.0);
+        this->c.push_back(0.875);
+        this->b1.push_back(37.0/378.0);
+        this->b1.push_back(0.0);
+        this->b1.push_back(250.0/621.0);
+        this->b1.push_back(125.0/594.0);
+        this->b1.push_back(0.0);
+        this->b1.push_back(512.0/1771.0);
+        this->b2.push_back(2825.0/27648.0);
+        this->b2.push_back(0.0);
+        this->b2.push_back(18575.0/48384.0);
+        this->b2.push_back(13525.0/55296.0);
+        this->b2.push_back(277.0/14336.0);
+        this->b2.push_back(0.25);
+
+        break;
+      }
+      default :
+      {
+        AssertThrow(false,ExcMessage("Unimplemented Embedded Runge-Kutta method."));
+      }
+      }
+  }
+
+
+
+  template <typename VECTOR>
+  void EmbeddedExplicitRungeKutta<VECTOR>::free_memory()
+  {
+    if (last_stage!=NULL)
+      delete last_stage;
+
+    last_stage = NULL;
+  }
+
+
+
+  template <typename VECTOR>
+  double EmbeddedExplicitRungeKutta<VECTOR>::evolve_one_time_step(
+    std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+    std_cxx1x::function<VECTOR (const double, const double, const VECTOR &)> id_minus_tau_J_inverse,
+    double t,
+    double delta_t,
+    VECTOR &y)
+  {
+    return evolve_one_time_step(f,t,delta_t,y);
+  }
+
+
+
+  template <typename VECTOR>
+  double EmbeddedExplicitRungeKutta<VECTOR>::evolve_one_time_step(
+    std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+    double t, double delta_t, VECTOR &y)
+  {
+    bool done = false;
+    unsigned int count = 0;
+    double error_norm = 0.;
+    VECTOR old_y(y);
+    VECTOR error(y);
+    std::vector<VECTOR> f_stages(this->n_stages,y);
+
+    while (!done)
+      {
+        error = 0.;
+        y = old_y;
+        // Compute the different stages needed.
+        compute_stages(f,t,delta_t,y,f_stages);
+
+        for (unsigned int i=0; i<this->n_stages; ++i)
+          {
+            y.sadd(1.,delta_t*this->b1[i],f_stages[i]);
+            error.sadd(1.,delta_t*(b2[i]-b1[i]),f_stages[i]);
+          }
+
+        error_norm = error.l2_norm();
+        // Check if the norm of error is less than the coarsening tolerance
+        if (error_norm<coarsen_tol)
+          {
+            done = true;
+            // Increase the guessed time step
+            double new_delta_t = delta_t*coarsen_param;
+            // Check that the guessed time step is smaller than the maximum time
+            // step allowed.
+            if (new_delta_t>max_delta_t)
+              {
+                status.exit_delta_t = MAX_DELTA_T;
+                status.delta_t_guess =  max_delta_t;
+              }
+            else
+              {
+                status.exit_delta_t = DELTA_T;
+                status.delta_t_guess = delta_t;
+              }
+          }
+        // Check if the norm of error is less than the refining tolerance
+        else if (error_norm<refine_tol)
+          {
+            done = true;
+            status.exit_delta_t = DELTA_T;
+            status.delta_t_guess = delta_t;
+          }
+        else
+          {
+            // If the time step is already the smallest acceptable, exit.
+            if (delta_t==min_delta_t)
+              {
+                done = true;
+                status.exit_delta_t = MIN_DELTA_T;
+                status.delta_t_guess = delta_t;
+              }
+            // Reduce the time step.
+            else
+              {
+                delta_t *= refine_param;
+                if (delta_t<min_delta_t)
+                  delta_t = min_delta_t;
+              }
+          }
+
+        ++count;
+      }
+
+    // Save the last stage if necessary
+    if (last_same_as_first==true)
+      {
+        if (last_stage==NULL)
+          last_stage = new VECTOR(f_stages.back());
+        else
+          *last_stage = f_stages.back();
+      }
+
+    status.n_iterations = count;
+    status.error_norm = error_norm;
+
+    return (t+delta_t);
+  }
+
+
+
+  template <typename VECTOR>
+  void EmbeddedExplicitRungeKutta<VECTOR>::set_time_adaptation_parameters(double coarsen_param_,
+      double refine_param_,
+      double min_delta_,
+      double max_delta_,
+      double refine_tol_,
+      double coarsen_tol_)
+  {
+    coarsen_param = coarsen_param_;
+    refine_param = refine_param_;
+    min_delta_t = min_delta_;
+    max_delta_t = max_delta_;
+    refine_tol = refine_tol_;
+    coarsen_tol = coarsen_tol_;
+  }
+
+
+
+  template <typename VECTOR>
+  const typename EmbeddedExplicitRungeKutta<VECTOR>::Status &EmbeddedExplicitRungeKutta<VECTOR>::get_status() const
+  {
+    return status;
+  }
+
+
+  template <typename VECTOR>
+  void EmbeddedExplicitRungeKutta<VECTOR>::compute_stages(
+    std_cxx1x::function<VECTOR (const double, const VECTOR &)> f,
+    const double t,
+    const double delta_t,
+    const VECTOR &y,
+    std::vector<VECTOR> &f_stages)
+  {
+    VECTOR Y(y);
+    unsigned int i = 0;
+
+    // If the last stage is the same as the first, we can skip the evaluation
+    // of the first stage.
+    if (last_same_as_first==true)
+      {
+        if (last_stage!=NULL)
+          {
+            f_stages[0] = *last_stage;
+            i = 1;
+          }
+      }
+
+    for (; i<this->n_stages; ++i)
+      {
+        Y = y;
+        for (unsigned int j = 0; j < i; ++j)
+          Y.sadd(1.0,delta_t*this->a[i][j],f_stages[j]);
+        f_stages[i] = f(t+this->c[i]*delta_t,Y);
+      }
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index ab491a5e57e1aa7556f93a1f8e55df86c4747ccc..6745cb7eaf711bb1764bf2a3b38c10da6d2ca820 100644 (file)
@@ -63,11 +63,13 @@ SET(_src
   tensor_product_polynomials_const.cc
   thread_management.cc
   timer.cc
+  time_stepping.cc
   utilities.cc
   )
 
 SET(_inst
   data_out_base.inst.in
+  time_stepping.inst.in
   )
 
 FILE(GLOB _header
diff --git a/deal.II/source/base/time_stepping.cc b/deal.II/source/base/time_stepping.cc
new file mode 100644 (file)
index 0000000..819658c
--- /dev/null
@@ -0,0 +1,35 @@
+// ---------------------------------------------------------------------
+// $Id: time_stepping.cc 30037 2013-07-18 16:55:40Z maier $
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/time_stepping.templates.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/lac/parallel_block_vector.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_block_vector.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+namespace TimeStepping
+{
+#include "time_stepping.inst"
+}
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/source/base/time_stepping.inst.in b/deal.II/source/base/time_stepping.inst.in
new file mode 100644 (file)
index 0000000..23e50fe
--- /dev/null
@@ -0,0 +1,44 @@
+// ---------------------------------------------------------------------
+// $Id: time_stepping.inst.in 30049 2013-07-18 19:42:40Z maier $
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
+{
+  template class ExplicitRungeKutta<V<S> >;
+  template class ImplicitRungeKutta<V<S> >;
+  template class EmbeddedExplicitRungeKutta<V<S> >;
+}
+
+for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
+{
+  template class ExplicitRungeKutta<parallel::distributed::V<S> >;
+  template class ImplicitRungeKutta<parallel::distributed::V<S> >;
+  template class EmbeddedExplicitRungeKutta<parallel::distributed::V<S> >;
+}
+
+for (V : EXTERNAL_SEQUENTIAL_VECTORS)
+{
+  template class ExplicitRungeKutta<V>;
+  template class ImplicitRungeKutta<V>;
+  template class EmbeddedExplicitRungeKutta<V>;
+}
+
+for (V : EXTERNAL_PARALLEL_VECTORS)
+{
+  template class ExplicitRungeKutta<V>;
+  template class ImplicitRungeKutta<V>;
+  template class EmbeddedExplicitRungeKutta<V>;
+}
diff --git a/tests/base/time_stepping_01.cc b/tests/base/time_stepping_01.cc
new file mode 100644 (file)
index 0000000..b689f16
--- /dev/null
@@ -0,0 +1,236 @@
+// ---------------------------------------------------------------------
+// $Id: time_stepping_01.cc 31349 2013-10-20 19:07:06Z maier $
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+#include "../tests.h"
+#include <deal.II/base/time_stepping.h>
+#include <deal.II/lac/vector.h>
+
+// test Runge-Kutta methods
+Vector<double> f1(double const t, Vector<double> const &y) 
+{ 
+  Vector<double> values(y);
+  for (unsigned int i=0; i<values.size(); ++i)
+    values[i] = 1.0;
+
+  return values;
+}
+
+Vector<double> f2(double const t, Vector<double> const &y) 
+{ 
+  Vector<double> values(y);
+  for (unsigned int i=0; i<values.size(); ++i)
+    values[i] = 2.0*t;
+
+  return values; 
+}
+
+Vector<double> f3(double const t, Vector<double> const &y) 
+{ 
+  Vector<double> values(y);
+  for (unsigned int i=0; i<values.size(); ++i)
+    values[i] = 3.0*t*t;
+
+  return values; 
+}
+
+Vector<double> f4(double const t, Vector<double> const &y) 
+{ 
+  Vector<double> values(y);
+  for (unsigned int i=0; i<values.size(); ++i)
+    values[i] = 4.0*t*t*t;
+  
+  return values; 
+}
+
+Vector<double> f5(double const t, Vector<double> const &y) 
+{
+  Vector<double> values(y);
+  for (unsigned int i=0; i<values.size(); ++i)
+    values[i] = 5.0*t*t*t*t;
+  
+  return values; 
+}
+
+Vector<double> id_minus_tau_J_inv1(double const t, double const tau, Vector<double> const &y) 
+{
+  return y;
+}
+
+Vector<double> id_minus_tau_J_inv2(double const t, double const tau, Vector<double> const &y) 
+{
+  return y;
+}
+
+Vector<double> id_minus_tau_J_inv3(double const t, double const tau, Vector<double> const &y) 
+{
+  return y;
+}
+
+Vector<double> id_minus_tau_J_inv4(double const t, double const tau, Vector<double> const &y) 
+{
+  return y;
+}
+
+Vector<double> id_minus_tau_J_inv5(double const t, double const tau, Vector<double> const &y) 
+{
+  return y;
+}
+
+double my1(double const t) 
+{ 
+  return t; 
+}
+
+double my2(double const t) 
+{ 
+  return t*t; 
+}
+
+double my3(double const t) 
+{ 
+  return t*t*t; 
+}
+
+double my4(double const t) 
+{ 
+  return t*t*t*t; 
+}
+
+double my5(double const t) 
+{ 
+  return t*t*t*t*t; 
+}
+
+void test(TimeStepping::RungeKutta<Vector<double> > &solver,
+    std::function<Vector<double> (double const, Vector<double> const &)> f,
+    std::function<Vector<double> (double const, double const, Vector<double> const &)> id_minus_tau_J_inv,
+    std::function<double (double const)> my)
+{
+  unsigned int n_time_steps = 1;
+  unsigned int size = 1;
+  double initial_time = 0.0, final_time = 1.0;
+  double time_step = (final_time-initial_time)/static_cast<double> (n_time_steps);
+  double time = initial_time;
+  Vector<double> solution(size);
+  Vector<double> exact_solution(size);
+  for (unsigned int i=0; i<size; ++i)
+  {
+    solution[i] = my(initial_time);
+    exact_solution[i] = my(final_time);
+  }
+
+  for (unsigned int i=0; i<n_time_steps; ++i)
+    time = solver.evolve_one_time_step(f,id_minus_tau_J_inv,time,time_step,solution);
+
+  Vector<double> error(exact_solution);
+  error.sadd(1.0,-1.0,solution);
+  double error_norm = error.l2_norm();
+  deallog << error_norm <<std::endl;
+}
+
+void test2(TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> > &solver,
+    std::function<Vector<double> (double const, Vector<double> const &)> f,
+    std::function<Vector<double> (double const, double const, Vector<double> const &)> id_minus_tau_J_inv,
+    std::function<double (double const)> my)
+{
+  double initial_time = 0.0, final_time = 1.0;
+  double time_step = 1.0;
+  unsigned int size = 1;
+  unsigned int n_time_steps = 0;
+  double time = initial_time;
+  Vector<double> solution(size);
+  Vector<double> exact_solution(size);
+  for (unsigned int i=0; i<size; ++i)
+  {
+    solution[i] = my(initial_time);
+    exact_solution[i] = my(final_time);
+  }
+
+
+  while (time<final_time)
+  {
+    if (time+time_step>final_time)
+      time_step = final_time - time;
+    time = solver.evolve_one_time_step(f,id_minus_tau_J_inv,time,time_step,solution);
+    time_step = solver.get_status().delta_t_guess;
+  }
+
+  Vector<double> error(exact_solution);
+  error.sadd(1.0,-1.0,solution);
+  double error_norm = error.l2_norm();
+  deallog<<error_norm<<std::endl;
+}
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  deallog<<"Forward Euler"<<std::endl;
+  TimeStepping::ExplicitRungeKutta<Vector<double> > fe(TimeStepping::FORWARD_EULER);
+  test(fe,f1,id_minus_tau_J_inv1,my1);
+
+  deallog<<"Runge-Kutta third order"<<std::endl;
+  TimeStepping::ExplicitRungeKutta<Vector<double> > rk3(TimeStepping::RK_THIRD_ORDER);
+  test(rk3,f3,id_minus_tau_J_inv3,my3);
+
+  deallog<<"Runge-Kutta fourth order"<<std::endl;
+  TimeStepping::ExplicitRungeKutta<Vector<double> > rk4(TimeStepping::RK_CLASSIC_FOURTH_ORDER);
+  test(rk4,f4,id_minus_tau_J_inv4,my4);
+
+  deallog<<"Backward Euler"<<std::endl;
+  TimeStepping::ImplicitRungeKutta<Vector<double> > be(TimeStepping::BACKWARD_EULER);
+  test(be,f1,id_minus_tau_J_inv1,my1);
+
+  deallog<<"Implicit midpoint"<<std::endl;
+  TimeStepping::ImplicitRungeKutta<Vector<double> > im(TimeStepping::IMPLICIT_MIDPOINT);
+  test(im,f2,id_minus_tau_J_inv2,my2);
+
+  deallog<<"Crank-Nicolson"<<std::endl;
+  TimeStepping::ImplicitRungeKutta<Vector<double> > cn(TimeStepping::CRANK_NICOLSON);
+  test(cn,f2,id_minus_tau_J_inv2,my2);
+
+  deallog<<"SDIRK"<<std::endl;
+  TimeStepping::ImplicitRungeKutta<Vector<double> > sdirk(TimeStepping::SDIRK_TWO_STAGES);
+  test(sdirk,f2,id_minus_tau_J_inv2,my2);
+
+  deallog<<"Heun-Euler"<<std::endl;
+  TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> > he(TimeStepping::HEUN_EULER);
+  test2(he,f2,id_minus_tau_J_inv2,my2);
+
+  deallog<<"Bogacki-Shampine"<<std::endl;
+  TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> > bs(TimeStepping::BOGACKI_SHAMPINE);
+  test2(bs,f3,id_minus_tau_J_inv3,my3);
+  bs.free_memory();
+
+  deallog<<"DOPRI"<<std::endl;
+  TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> > dopri(TimeStepping::DOPRI);
+  test2(dopri,f5,id_minus_tau_J_inv5,my5);
+  dopri.free_memory();
+
+  deallog<<"Fehlberg"<<std::endl;
+  TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> > fehlberg(TimeStepping::FEHLBERG);
+  test2(fehlberg,f5,id_minus_tau_J_inv5,my5);
+
+  deallog<<"Cash-Karp"<<std::endl;
+  TimeStepping::EmbeddedExplicitRungeKutta<Vector<double> > ck(TimeStepping::CASH_KARP);
+  test2(ck,f5,id_minus_tau_J_inv5,my5);
+
+  return 0;
+}
diff --git a/tests/base/time_stepping_01.output b/tests/base/time_stepping_01.output
new file mode 100644 (file)
index 0000000..78c9c17
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL::Forward Euler
+DEAL::0
+DEAL::Runge-Kutta third order
+DEAL::0
+DEAL::Runge-Kutta fourth order
+DEAL::0
+DEAL::Backward Euler
+DEAL::0
+DEAL::Implicit midpoint
+DEAL::0
+DEAL::Crank-Nicolson
+DEAL::0
+DEAL::SDIRK
+DEAL::0
+DEAL::Heun-Euler
+DEAL::0
+DEAL::Bogacki-Shampine
+DEAL::0
+DEAL::DOPRI
+DEAL::0
+DEAL::Fehlberg
+DEAL::0
+DEAL::Cash-Karp
+DEAL::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.