From: Bruno Turcksin Date: Tue, 18 Feb 2014 04:08:08 +0000 (+0000) Subject: Add Runge-Kutta methods for time integration. X-Git-Tag: v8.2.0-rc1~815 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ab14bd290208f9d4e9653e802d43808473dc2f3f;p=dealii.git Add Runge-Kutta methods for time integration. git-svn-id: https://svn.dealii.org/trunk@32500 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index cc1f4f4b89..beb94ec9ee 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -135,6 +135,12 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Damien Lebrun-Grandie, Bruno Turcksin, 2014/02/17) +
  2. New: There is now a class FEEvaluationDGP that implements matrix-free evaluation routines by truncated tensor products for FE_DGP elements.
    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 index 0000000000..1730cc81b5 --- /dev/null +++ b/deal.II/include/deal.II/base/time_stepping.h @@ -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 +#include + +#include + +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 + 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 > &F, + std::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 + class RungeKutta : public TimeStepping + { + 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 > &F, + std::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 f, + std_cxx1x::function 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 b; + + /** + * Butcher tableau coefficients. + */ + std::vector c; + + /** + * Butcher tableau coefficients. + */ + std::vector > a; + }; + + + + /** + * ExplicitRungeKutta is derived from RungeKutta and implement the explicit methods. + */ + template + class ExplicitRungeKutta : public RungeKutta + { + 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 f, + std_cxx1x::function 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 f, + double t, + double delta_t, + VECTOR &y); + + /** + * This structure stores the name of the method used. + */ + struct Status : public TimeStepping::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 f, + const double t, + const double delta_t, + const VECTOR &y, + std::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 + class ImplicitRungeKutta : public RungeKutta + { + 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 f, + std_cxx1x::function 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::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 f, + std_cxx1x::function id_minus_tau_J_inverse, + double t, + double delta_t, + VECTOR &y, + std::vector &f_stages); + + /** + * Newton solver used for the implicit stages. + */ + void newton_solve(std_cxx1x::function get_residual, + std_cxx1x::function id_minus_tau_J_inverse, + VECTOR &y); + + /** + * Compute the residual needed by the Newton solver. + */ + void compute_residual(std_cxx1x::function 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 + class EmbeddedExplicitRungeKutta : public RungeKutta + { + 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 f, + std_cxx1x::function 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 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::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 f, + const double t, + const double delta_t, + const VECTOR &y, + std::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 b1; + + /** + * Butcher tableau coefficients. + */ + std::vector 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 index 0000000000..6a2c98a5f2 --- /dev/null +++ b/deal.II/include/deal.II/base/time_stepping.templates.h @@ -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 +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace TimeStepping +{ + // ---------------------------------------------------------------------- + // RungeKutta + // ---------------------------------------------------------------------- + + template + double RungeKutta::evolve_one_time_step( + std::vector > &F, + std::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 + ExplicitRungeKutta::ExplicitRungeKutta(runge_kutta_method method) + { + initialize(method); + } + + + + template + void ExplicitRungeKutta::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 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 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 + double ExplicitRungeKutta::evolve_one_time_step( + std_cxx1x::function f, + std_cxx1x::function id_minus_tau_J_inverse, + double t, + double delta_t, + VECTOR &y) + { + return evolve_one_time_step(f,t,delta_t,y); + } + + + + template + double ExplicitRungeKutta::evolve_one_time_step( + std_cxx1x::function f, + double t, + double delta_t, + VECTOR &y) + { + std::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; in_stages; ++i) + y.sadd(1.,delta_t*this->b[i],f_stages[i]); + + return (t+delta_t); + } + + + + template + const typename ExplicitRungeKutta::Status &ExplicitRungeKutta::get_status() const + { + return status; + } + + + + template + void ExplicitRungeKutta::compute_stages( + std_cxx1x::function f, + const double t, + const double delta_t, + const VECTOR &y, + std::vector &f_stages) const + { + for (unsigned int i=0; in_stages; ++i) + { + VECTOR Y(y); + for (unsigned int j=0; ja[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 + ImplicitRungeKutta::ImplicitRungeKutta(runge_kutta_method method, + unsigned int max_it, + double tolerance) + : + RungeKutta (), + skip_linear_combi(false), + max_it(max_it), + tolerance(tolerance) + { + initialize(method); + } + + + + template + void ImplicitRungeKutta::initialize(runge_kutta_method method) + { + status.method = method; + + switch (method) + { + case (BACKWARD_EULER) : + { + this->n_stages = 1; + this->a.push_back(std::vector(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(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(1, 0.0)); + this->a.push_back(std::vector(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(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 + double ImplicitRungeKutta::evolve_one_time_step( + std_cxx1x::function f, + std_cxx1x::function id_minus_tau_J_inverse, + double t, + double delta_t, + VECTOR &y) + { + VECTOR old_y(y); + std::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; in_stages; ++i) + y.sadd(1.,delta_t*this->b[i],f_stages[i]); + } + + return (t+delta_t); + } + + + + template + void ImplicitRungeKutta::set_newton_solver_parameters(unsigned int max_it_, double tolerance_) + { + max_it = max_it_; + tolerance = tolerance_; + } + + + + template + const typename ImplicitRungeKutta::Status &ImplicitRungeKutta::get_status() const + { + return status; + } + + + + template + void ImplicitRungeKutta::compute_stages( + std_cxx1x::function f, + std_cxx1x::function id_minus_tau_J_inverse, + double t, + double delta_t, + VECTOR &y, + std::vector &f_stages) + { + VECTOR z(y); + for (unsigned int i=0; in_stages; ++i) + { + VECTOR old_y(z); + for (unsigned int j=0; ja[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::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 + void ImplicitRungeKutta::newton_solve( + std_cxx1x::function get_residual, + std_cxx1x::function 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 + void ImplicitRungeKutta::compute_residual( + std_cxx1x::function 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 + EmbeddedExplicitRungeKutta::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 + void EmbeddedExplicitRungeKutta::initialize(runge_kutta_method method) + { + status.method = method; + + switch (method) + { + case (HEUN_EULER) : + { + this->n_stages = 2; + this->a.push_back(std::vector()); + this->a.push_back(std::vector(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 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 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 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 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 + void EmbeddedExplicitRungeKutta::free_memory() + { + if (last_stage!=NULL) + delete last_stage; + + last_stage = NULL; + } + + + + template + double EmbeddedExplicitRungeKutta::evolve_one_time_step( + std_cxx1x::function f, + std_cxx1x::function id_minus_tau_J_inverse, + double t, + double delta_t, + VECTOR &y) + { + return evolve_one_time_step(f,t,delta_t,y); + } + + + + template + double EmbeddedExplicitRungeKutta::evolve_one_time_step( + std_cxx1x::function 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 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; in_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_normmax_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 + void EmbeddedExplicitRungeKutta::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 + const typename EmbeddedExplicitRungeKutta::Status &EmbeddedExplicitRungeKutta::get_status() const + { + return status; + } + + + template + void EmbeddedExplicitRungeKutta::compute_stages( + std_cxx1x::function f, + const double t, + const double delta_t, + const VECTOR &y, + std::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 (; in_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 diff --git a/deal.II/source/base/CMakeLists.txt b/deal.II/source/base/CMakeLists.txt index ab491a5e57..6745cb7eaf 100644 --- a/deal.II/source/base/CMakeLists.txt +++ b/deal.II/source/base/CMakeLists.txt @@ -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 index 0000000000..819658cea1 --- /dev/null +++ b/deal.II/source/base/time_stepping.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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 index 0000000000..23e50fe3b8 --- /dev/null +++ b/deal.II/source/base/time_stepping.inst.in @@ -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 >; + template class ImplicitRungeKutta >; + template class EmbeddedExplicitRungeKutta >; +} + +for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) +{ + template class ExplicitRungeKutta >; + template class ImplicitRungeKutta >; + template class EmbeddedExplicitRungeKutta >; +} + +for (V : EXTERNAL_SEQUENTIAL_VECTORS) +{ + template class ExplicitRungeKutta; + template class ImplicitRungeKutta; + template class EmbeddedExplicitRungeKutta; +} + +for (V : EXTERNAL_PARALLEL_VECTORS) +{ + template class ExplicitRungeKutta; + template class ImplicitRungeKutta; + template class EmbeddedExplicitRungeKutta; +} diff --git a/tests/base/time_stepping_01.cc b/tests/base/time_stepping_01.cc new file mode 100644 index 0000000000..b689f16a33 --- /dev/null +++ b/tests/base/time_stepping_01.cc @@ -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 +#include + +// test Runge-Kutta methods +Vector f1(double const t, Vector const &y) +{ + Vector values(y); + for (unsigned int i=0; i f2(double const t, Vector const &y) +{ + Vector values(y); + for (unsigned int i=0; i f3(double const t, Vector const &y) +{ + Vector values(y); + for (unsigned int i=0; i f4(double const t, Vector const &y) +{ + Vector values(y); + for (unsigned int i=0; i f5(double const t, Vector const &y) +{ + Vector values(y); + for (unsigned int i=0; i id_minus_tau_J_inv1(double const t, double const tau, Vector const &y) +{ + return y; +} + +Vector id_minus_tau_J_inv2(double const t, double const tau, Vector const &y) +{ + return y; +} + +Vector id_minus_tau_J_inv3(double const t, double const tau, Vector const &y) +{ + return y; +} + +Vector id_minus_tau_J_inv4(double const t, double const tau, Vector const &y) +{ + return y; +} + +Vector id_minus_tau_J_inv5(double const t, double const tau, Vector 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 > &solver, + std::function (double const, Vector const &)> f, + std::function (double const, double const, Vector const &)> id_minus_tau_J_inv, + std::function 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 (n_time_steps); + double time = initial_time; + Vector solution(size); + Vector exact_solution(size); + for (unsigned int i=0; i error(exact_solution); + error.sadd(1.0,-1.0,solution); + double error_norm = error.l2_norm(); + deallog << error_norm < > &solver, + std::function (double const, Vector const &)> f, + std::function (double const, double const, Vector const &)> id_minus_tau_J_inv, + std::function 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 solution(size); + Vector exact_solution(size); + for (unsigned int i=0; ifinal_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 error(exact_solution); + error.sadd(1.0,-1.0,solution); + double error_norm = error.l2_norm(); + deallog< > fe(TimeStepping::FORWARD_EULER); + test(fe,f1,id_minus_tau_J_inv1,my1); + + deallog<<"Runge-Kutta third order"< > rk3(TimeStepping::RK_THIRD_ORDER); + test(rk3,f3,id_minus_tau_J_inv3,my3); + + deallog<<"Runge-Kutta fourth order"< > rk4(TimeStepping::RK_CLASSIC_FOURTH_ORDER); + test(rk4,f4,id_minus_tau_J_inv4,my4); + + deallog<<"Backward Euler"< > be(TimeStepping::BACKWARD_EULER); + test(be,f1,id_minus_tau_J_inv1,my1); + + deallog<<"Implicit midpoint"< > im(TimeStepping::IMPLICIT_MIDPOINT); + test(im,f2,id_minus_tau_J_inv2,my2); + + deallog<<"Crank-Nicolson"< > cn(TimeStepping::CRANK_NICOLSON); + test(cn,f2,id_minus_tau_J_inv2,my2); + + deallog<<"SDIRK"< > sdirk(TimeStepping::SDIRK_TWO_STAGES); + test(sdirk,f2,id_minus_tau_J_inv2,my2); + + deallog<<"Heun-Euler"< > he(TimeStepping::HEUN_EULER); + test2(he,f2,id_minus_tau_J_inv2,my2); + + deallog<<"Bogacki-Shampine"< > bs(TimeStepping::BOGACKI_SHAMPINE); + test2(bs,f3,id_minus_tau_J_inv3,my3); + bs.free_memory(); + + deallog<<"DOPRI"< > dopri(TimeStepping::DOPRI); + test2(dopri,f5,id_minus_tau_J_inv5,my5); + dopri.free_memory(); + + deallog<<"Fehlberg"< > fehlberg(TimeStepping::FEHLBERG); + test2(fehlberg,f5,id_minus_tau_J_inv5,my5); + + deallog<<"Cash-Karp"< > 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 index 0000000000..78c9c17007 --- /dev/null +++ b/tests/base/time_stepping_01.output @@ -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