<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>
--- /dev/null
+// ---------------------------------------------------------------------
+// $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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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
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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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>;
+}
--- /dev/null
+// ---------------------------------------------------------------------
+// $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;
+}
--- /dev/null
+
+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
projects / dealii-svn.git / commitdiff