]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add low-storage rk and update step-67 11282/head
authorJiaqi Zhang <jiaqi2@clemson.edu>
Mon, 30 Nov 2020 18:14:07 +0000 (13:14 -0500)
committerJiaqi Zhang <jiaqi2@clemson.edu>
Fri, 11 Dec 2020 16:52:14 +0000 (16:52 +0000)
doc/news/changes/minor/20201118JiaqiZhang [new file with mode: 0644]
examples/step-67/step-67.cc
include/deal.II/base/time_stepping.h
include/deal.II/base/time_stepping.templates.h
source/base/time_stepping.inst.in
tests/base/time_stepping_01.cc
tests/base/time_stepping_01.output
tests/base/time_stepping_03.cc [new file with mode: 0644]
tests/base/time_stepping_03.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20201118JiaqiZhang b/doc/news/changes/minor/20201118JiaqiZhang
new file mode 100644 (file)
index 0000000..0054c73
--- /dev/null
@@ -0,0 +1,4 @@
+New: A new class LowStorageRungeKutta is added to the namespace TimeStepping to 
+implement the explicit low-storage Runge-Kutta methods, see @cite KennedyCarpenterLewis2000 and step-67.
+<br>
+(Jiaqi Zhang, 2020/11/18)
index d1685c5bea1009eb56d482b06d13dbffdb8398f6..9186562b23c2d3487c9f7fcb9a4a84e4e7eb6b8a 100644 (file)
@@ -23,6 +23,7 @@
 #include <deal.II/base/function.h>
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/timer.h>
+#include <deal.II/base/time_stepping.h>
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/vectorization.h>
 
@@ -250,6 +251,7 @@ namespace Euler_DG
   public:
     LowStorageRungeKuttaIntegrator(const LowStorageRungeKuttaScheme scheme)
     {
+      TimeStepping::runge_kutta_method lsrk;
       // First comes the three-stage scheme of order three by Kennedy et al.
       // (2000). While its stability region is significantly smaller than for
       // the other schemes, it only involves three stages, so it is very
@@ -258,9 +260,7 @@ namespace Euler_DG
         {
           case stage_3_order_3:
             {
-              bi = {{0.245170287303492, 0.184896052186740, 0.569933660509768}};
-              ai = {{0.755726351946097, 0.386954477304099}};
-
+              lsrk = TimeStepping::LOW_STORAGE_RK_STAGE3_ORDER3;
               break;
             }
 
@@ -268,16 +268,7 @@ namespace Euler_DG
             // defined in the paper by Kennedy et al. (2000).
           case stage_5_order_4:
             {
-              bi = {{1153189308089. / 22510343858157.,
-                     1772645290293. / 4653164025191.,
-                     -1672844663538. / 4480602732383.,
-                     2114624349019. / 3568978502595.,
-                     5198255086312. / 14908931495163.}};
-              ai = {{970286171893. / 4311952581923.,
-                     6584761158862. / 12103376702013.,
-                     2251764453980. / 15575788980749.,
-                     26877169314380. / 34165994151039.}};
-
+              lsrk = TimeStepping::LOW_STORAGE_RK_STAGE5_ORDER4;
               break;
             }
 
@@ -296,20 +287,7 @@ namespace Euler_DG
             // flux.
           case stage_7_order_4:
             {
-              bi = {{0.0941840925477795334,
-                     0.149683694803496998,
-                     0.285204742060440058,
-                     -0.122201846148053668,
-                     0.0605151571191401122,
-                     0.345986987898399296,
-                     0.186627171718797670}};
-              ai = {{0.241566650129646868 + bi[0],
-                     0.0423866513027719953 + bi[1],
-                     0.215602732678803776 + bi[2],
-                     0.232328007537583987 + bi[3],
-                     0.256223412574146438 + bi[4],
-                     0.0978694102142697230 + bi[5]}};
-
+              lsrk = TimeStepping::LOW_STORAGE_RK_STAGE7_ORDER4;
               break;
             }
 
@@ -320,30 +298,17 @@ namespace Euler_DG
             // stage is less than for the fourth order schemes.
           case stage_9_order_5:
             {
-              bi = {{2274579626619. / 23610510767302.,
-                     693987741272. / 12394497460941.,
-                     -347131529483. / 15096185902911.,
-                     1144057200723. / 32081666971178.,
-                     1562491064753. / 11797114684756.,
-                     13113619727965. / 44346030145118.,
-                     393957816125. / 7825732611452.,
-                     720647959663. / 6565743875477.,
-                     3559252274877. / 14424734981077.}};
-              ai = {{1107026461565. / 5417078080134.,
-                     38141181049399. / 41724347789894.,
-                     493273079041. / 11940823631197.,
-                     1851571280403. / 6147804934346.,
-                     11782306865191. / 62590030070788.,
-                     9452544825720. / 13648368537481.,
-                     4435885630781. / 26285702406235.,
-                     2357909744247. / 11371140753790.}};
-
+              lsrk = TimeStepping::LOW_STORAGE_RK_STAGE9_ORDER5;
               break;
             }
 
           default:
             AssertThrow(false, ExcNotImplemented());
         }
+      TimeStepping::LowStorageRungeKutta<
+        LinearAlgebra::distributed::Vector<Number>>
+        rk_integrator(lsrk);
+      rk_integrator.get_coefficients(ai, bi, ci);
     }
 
     unsigned int n_stages() const
@@ -386,10 +351,10 @@ namespace Euler_DG
                                  vec_ri,
                                  solution,
                                  vec_ri);
-      double sum_previous_bi = 0;
+
       for (unsigned int stage = 1; stage < bi.size(); ++stage)
         {
-          const double c_i = sum_previous_bi + ai[stage - 1];
+          const double c_i = ci[stage];
           pde_operator.perform_stage(current_time + c_i * time_step,
                                      bi[stage] * time_step,
                                      (stage == bi.size() - 1 ?
@@ -399,13 +364,13 @@ namespace Euler_DG
                                      vec_ki,
                                      solution,
                                      vec_ri);
-          sum_previous_bi += bi[stage - 1];
         }
     }
 
   private:
     std::vector<double> bi;
     std::vector<double> ai;
+    std::vector<double> ci;
   };
 
 
index 5b1f418c28531c35f00872d9b06baab32141660a..8cc69fcb91aabeccebf42acc360db1e8725ad623 100644 (file)
@@ -39,6 +39,11 @@ namespace TimeStepping
    *   - RK_THIRD_ORDER (third order Runge-Kutta)
    *   - SSP_THIRD_ORDER (third order SSP Runge-Kutta)
    *   - RK_CLASSIC_FOURTH_ORDER (classical fourth order Runge-Kutta)
+   * - Low-storage (explicit) Runge-Kutta methods
+   *   - LOW_STORAGE_RK_STAGE3_ORDER3 (Three stages and third order)
+   *   - LOW_STORAGE_RK_STAGE5_ORDER4 (Five stages and fourth order)
+   *   - LOW_STORAGE_RK_STAGE7_ORDER4 (Seven stages and fourth order)
+   *   - LOW_STORAGE_RK_STAGE9_ORDER5 (Nine stages and fifth order)
    * - Implicit methods (see ImplicitRungeKutta::initialize):
    *   - BACKWARD_EULER (first order)
    *   - IMPLICIT_MIDPOINT (second order)
@@ -72,6 +77,28 @@ namespace TimeStepping
      * Classical fourth order Runge-Kutta method.
      */
     RK_CLASSIC_FOURTH_ORDER,
+    /**
+     * Three-stage scheme of order three by Kennedy et al.
+     * @cite KennedyCarpenterLewis2000. Its stability region is
+     * significantly smaller than the higher order schemes, but due to three
+     * stages only, it is very competitive in terms of the work per stage.
+     */
+    LOW_STORAGE_RK_STAGE3_ORDER3,
+    /**
+     * Five-stage scheme of order four,
+     * defined in the paper by Kennedy et al. @cite KennedyCarpenterLewis2000.
+     */
+    LOW_STORAGE_RK_STAGE5_ORDER4,
+    /**
+     * Seven-stage scheme of order four defined in the paper by Tselios and
+     * Simos @cite TseliosSimos2007.
+     */
+    LOW_STORAGE_RK_STAGE7_ORDER4,
+    /**
+     * Nine-stage scheme of order five
+     * defined in the paper by Kennedy et al. @cite KennedyCarpenterLewis2000.
+     */
+    LOW_STORAGE_RK_STAGE9_ORDER5,
     /**
      * Backward Euler method, first order.
      */
@@ -373,6 +400,126 @@ namespace TimeStepping
 
 
 
+  /**
+   * The LowStorageRungeKutta class is derived from RungeKutta and implements a
+   * specific class of explicit methods. The main advantages of low-storage
+   * methods are the reduced memory consumption and the reduced memory access.
+   */
+  template <typename VectorType>
+  class LowStorageRungeKutta : public RungeKutta<VectorType>
+  {
+  public:
+    using RungeKutta<VectorType>::evolve_one_time_step;
+
+    /**
+     * Default constructor. This constructor creates an object for which
+     * you will want to call <code>initialize(runge_kutta_method)</code>
+     * before it can be used.
+     */
+    LowStorageRungeKutta() = default;
+
+    /**
+     * Constructor. This function calls initialize(runge_kutta_method).
+     */
+    LowStorageRungeKutta(const runge_kutta_method method);
+
+    /**
+     * Initialize the explicit Runge-Kutta method.
+     */
+    void
+    initialize(const runge_kutta_method method) override;
+
+    /**
+     * This function is used to advance from time @p t to t+ @p delta_t. @p f
+     * is the function $ f(t,y) $ 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 $
+     * inv(I-\tau J)$ where $ I $ is the identity matrix, $ \tau $ is given,
+     * and $ J $ is the Jacobian $ \frac{\partial J}{\partial y} $. The input
+     * parameters are the time, $ \tau $, 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(
+      const std::function<VectorType(const double, const VectorType &)> &f,
+      const std::function<
+        VectorType(const double, const double, const VectorType &)>
+        &         id_minus_tau_J_inverse,
+      double      t,
+      double      delta_t,
+      VectorType &y) override;
+
+    /**
+     * This function is used to advance from time @p t to t+ @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. Note that vec_ki holds the evaluation of the differential operator,
+     * and vec_ri holds the right-hand side for the differential operator
+     * application.
+     */
+    double
+    evolve_one_time_step(
+      const std::function<VectorType(const double, const VectorType &)> &f,
+      double                                                             t,
+      double      delta_t,
+      VectorType &solution,
+      VectorType &vec_ri,
+      VectorType &vec_ki);
+
+    /**
+     * Get the coefficients of the scheme.
+     * Note that here vector @p a is not the conventional definition in terms of a
+     * Butcher tableau but merely one of the sub-diagonals. More details can be
+     * found in step-67 and the references therein.
+     */
+    void
+    get_coefficients(std::vector<double> &a,
+                     std::vector<double> &b,
+                     std::vector<double> &c) const;
+
+    /**
+     * This structure stores the name of the method used.
+     */
+    struct Status : public TimeStepping<VectorType>::Status
+    {
+      Status()
+        : method(invalid)
+      {}
+
+      runge_kutta_method method;
+    };
+
+    /**
+     * Return the status of the current object.
+     */
+    const Status &
+    get_status() const override;
+
+  private:
+    /**
+     * Compute  one stage of low storage rk.
+     */
+    void
+    compute_one_stage(
+      const std::function<VectorType(const double, const VectorType &)> &f,
+      const double                                                       t,
+      const double      factor_solution,
+      const double      factor_ai,
+      const VectorType &corrent_ri,
+      VectorType &      vec_ki,
+      VectorType &      solution,
+      VectorType &      next_ri) 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.
index 7a363a756225ae9f8208b2a6b676294b36be8b2f..1998d005d39d0802cb8fbc00e4baba5069cae0e3 100644 (file)
@@ -241,6 +241,234 @@ namespace TimeStepping
 
 
 
+  // ----------------------------------------------------------------------
+  // LowStorageRungeKutta
+  // ----------------------------------------------------------------------
+
+  template <typename VectorType>
+  LowStorageRungeKutta<VectorType>::LowStorageRungeKutta(
+    const runge_kutta_method method)
+  {
+    // virtual functions called in constructors and destructors never use the
+    // override in a derived class
+    // for clarity be explicit on which function is called
+    LowStorageRungeKutta<VectorType>::initialize(method);
+  }
+
+
+
+  template <typename VectorType>
+  void
+  LowStorageRungeKutta<VectorType>::initialize(const runge_kutta_method method)
+  {
+    status.method = method;
+
+    switch (method)
+      {
+        case (LOW_STORAGE_RK_STAGE3_ORDER3):
+          {
+            this->n_stages = 3;
+            this->b.reserve(this->n_stages);
+            this->b.push_back(0.245170287303492);
+            this->b.push_back(0.184896052186740);
+            this->b.push_back(0.569933660509768);
+
+            std::vector<double> tmp;
+            tmp = {{0.755726351946097, 0.386954477304099}};
+            this->a.push_back(tmp);
+            break;
+          }
+        case (LOW_STORAGE_RK_STAGE5_ORDER4):
+          {
+            this->n_stages = 5;
+            this->b        = {{1153189308089. / 22510343858157.,
+                        1772645290293. / 4653164025191.,
+                        -1672844663538. / 4480602732383.,
+                        2114624349019. / 3568978502595.,
+                        5198255086312. / 14908931495163.}};
+            std::vector<double> ai;
+            ai = {{970286171893. / 4311952581923.,
+                   6584761158862. / 12103376702013.,
+                   2251764453980. / 15575788980749.,
+                   26877169314380. / 34165994151039.}};
+            this->a.push_back(ai);
+            break;
+          }
+        case (LOW_STORAGE_RK_STAGE7_ORDER4):
+          {
+            this->n_stages = 7;
+            this->b        = {{0.0941840925477795334,
+                        0.149683694803496998,
+                        0.285204742060440058,
+                        -0.122201846148053668,
+                        0.0605151571191401122,
+                        0.345986987898399296,
+                        0.186627171718797670}};
+            std::vector<double> ai;
+            ai = {{0.241566650129646868 + this->b[0],
+                   0.0423866513027719953 + this->b[1],
+                   0.215602732678803776 + this->b[2],
+                   0.232328007537583987 + this->b[3],
+                   0.256223412574146438 + this->b[4],
+                   0.0978694102142697230 + this->b[5]}};
+            this->a.push_back(ai);
+            break;
+          }
+        case (LOW_STORAGE_RK_STAGE9_ORDER5):
+          {
+            this->n_stages = 9;
+            this->b        = {{2274579626619. / 23610510767302.,
+                        693987741272. / 12394497460941.,
+                        -347131529483. / 15096185902911.,
+                        1144057200723. / 32081666971178.,
+                        1562491064753. / 11797114684756.,
+                        13113619727965. / 44346030145118.,
+                        393957816125. / 7825732611452.,
+                        720647959663. / 6565743875477.,
+                        3559252274877. / 14424734981077.}};
+            std::vector<double> ai;
+            ai = {{1107026461565. / 5417078080134.,
+                   38141181049399. / 41724347789894.,
+                   493273079041. / 11940823631197.,
+                   1851571280403. / 6147804934346.,
+                   11782306865191. / 62590030070788.,
+                   9452544825720. / 13648368537481.,
+                   4435885630781. / 26285702406235.,
+                   2357909744247. / 11371140753790.}};
+            this->a.push_back(ai);
+            break;
+          }
+        default:
+          {
+            AssertThrow(false,
+                        ExcMessage(
+                          "Unimplemented low-storage Runge-Kutta method."));
+          }
+      }
+    // compute ci
+    this->c.reserve(this->n_stages);
+    this->c.push_back(0.);
+    double sum_previous_bi = 0.;
+    for (unsigned int stage = 1; stage < this->n_stages; ++stage)
+      {
+        const double tmp = sum_previous_bi + this->a[0][stage - 1];
+        this->c.push_back(tmp);
+        sum_previous_bi += this->b[stage - 1];
+      }
+  }
+
+
+
+  template <typename VectorType>
+  double
+  LowStorageRungeKutta<VectorType>::evolve_one_time_step(
+    const std::function<VectorType(const double, const VectorType &)> &f,
+    const std::function<
+      VectorType(const double, const double, const VectorType &)>
+      & /*id_minus_tau_J_inverse*/,
+    double      t,
+    double      delta_t,
+    VectorType &y)
+  {
+    // We need two auxiliary vectors, namely the vector ki
+    // to hold the evaluation of the differential operator, and the vector ri
+    // that holds the right-hand side for the differential operator application.
+    VectorType vec_ri;
+    VectorType vec_ki;
+    return evolve_one_time_step(f, t, delta_t, y, vec_ri, vec_ki);
+  }
+
+
+
+  template <typename VectorType>
+  double
+  LowStorageRungeKutta<VectorType>::evolve_one_time_step(
+    const std::function<VectorType(const double, const VectorType &)> &f,
+    double                                                             t,
+    double                                                             delta_t,
+    VectorType &                                                       solution,
+    VectorType &                                                       vec_ri,
+    VectorType &                                                       vec_ki)
+  {
+    compute_one_stage(f,
+                      t,
+                      this->b[0] * delta_t,
+                      this->a[0][0] * delta_t,
+                      solution,
+                      vec_ki,
+                      solution,
+                      vec_ri);
+
+    for (unsigned int stage = 1; stage < this->n_stages; ++stage)
+      {
+        const double c_i = this->c[stage];
+        const double factor_ai =
+          (stage == this->n_stages - 1 ? 0 : this->a[0][stage] * delta_t);
+        compute_one_stage(f,
+                          t + c_i * delta_t,
+                          this->b[stage] * delta_t,
+                          factor_ai,
+                          vec_ri,
+                          vec_ki,
+                          solution,
+                          vec_ri);
+      }
+    return (t + delta_t);
+  }
+
+  template <typename VectorType>
+  void
+  LowStorageRungeKutta<VectorType>::get_coefficients(
+    std::vector<double> &a,
+    std::vector<double> &b,
+    std::vector<double> &c) const
+  {
+    a.resize(this->a[0].size());
+    a = this->a[0];
+
+    b.resize(this->b.size());
+    b = this->b;
+
+    c.resize(this->c.size());
+    c = this->c;
+  }
+
+  template <typename VectorType>
+  const typename LowStorageRungeKutta<VectorType>::Status &
+  LowStorageRungeKutta<VectorType>::get_status() const
+  {
+    return status;
+  }
+
+  template <typename VectorType>
+  void
+  LowStorageRungeKutta<VectorType>::compute_one_stage(
+    const std::function<VectorType(const double, const VectorType &)> &f,
+    const double                                                       t,
+    const double      factor_solution,
+    const double      factor_ai,
+    const VectorType &current_ri,
+    VectorType &      vec_ki,
+    VectorType &      solution,
+    VectorType &      next_ri) const
+  {
+    const double ai = factor_ai;
+    const double bi = factor_solution;
+    vec_ki          = f(t, current_ri);
+
+    if (ai == double())
+      {
+        solution.sadd(1., bi, vec_ki);
+      }
+    else
+      {
+        next_ri = solution;
+        next_ri.sadd(1., ai, vec_ki);
+        solution.add(bi, vec_ki);
+      }
+  }
+
+
   // ----------------------------------------------------------------------
   // ImplicitRungeKutta
   // ----------------------------------------------------------------------
index 6315226010e52a04fea4f7880efc474879fe6857..e9667475a7b18dfa8638cd7c233f4631f30ee8a7 100644 (file)
@@ -18,6 +18,7 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
   {
     template class RungeKutta<V<S>>;
     template class ExplicitRungeKutta<V<S>>;
+    template class LowStorageRungeKutta<V<S>>;
     template class ImplicitRungeKutta<V<S>>;
     template class EmbeddedExplicitRungeKutta<V<S>>;
   }
@@ -26,6 +27,7 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
   {
     template class RungeKutta<LinearAlgebra::distributed::V<S>>;
     template class ExplicitRungeKutta<LinearAlgebra::distributed::V<S>>;
+    template class LowStorageRungeKutta<LinearAlgebra::distributed::V<S>>;
     template class ImplicitRungeKutta<LinearAlgebra::distributed::V<S>>;
     template class EmbeddedExplicitRungeKutta<LinearAlgebra::distributed::V<S>>;
   }
@@ -34,6 +36,7 @@ for (V : EXTERNAL_PARALLEL_VECTORS)
   {
     template class RungeKutta<V>;
     template class ExplicitRungeKutta<V>;
+    template class LowStorageRungeKutta<V>;
     template class ImplicitRungeKutta<V>;
     template class EmbeddedExplicitRungeKutta<V>;
   }
index a7d7b51dba985394ece5d1ffeb72b991ee9480ce..0fd9bdf48017edd3d914c81f09baf002902ac713 100644 (file)
@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 // test Runge-Kutta methods in TimeStepping with a) a polynomial with expected
-// error 0 and b) convergence order for y=exp(t^2)
+// error 0 and b) convergence order for y=0.1*exp(t^2)
 #include <deal.II/base/time_stepping.h>
 
 #include <deal.II/lac/vector.h>
@@ -145,7 +145,7 @@ my5(double const t)
 double
 my_exact_solution(double const t)
 {
-  return std::exp(t * t);
+  return 0.1 * std::exp(t * t);
 }
 
 void
@@ -238,7 +238,7 @@ test_convergence(
     }
 
   deallog << "convergence rate" << std::endl;
-  for (unsigned int cycle = 0; cycle < 10; ++cycle)
+  for (unsigned int cycle = 0; cycle < 8; ++cycle)
     {
       unsigned int n_time_steps = std::pow(2., static_cast<double>(cycle));
       double       time_step =
@@ -255,7 +255,7 @@ test_convergence(
       error.sadd(1.0, -1.0, solution);
       double error_norm = error.l2_norm();
       errors.push_back(error_norm);
-      if (cycle > 0)
+      if (cycle > 1)
         deallog << std::log(std::fabs(errors[cycle - 1] / errors[cycle])) /
                      std::log(2.)
                 << std::endl;
@@ -266,6 +266,7 @@ int
 main()
 {
   initlog();
+  // deallog.precision(4);
   {
     deallog << "Forward Euler" << std::endl;
     TimeStepping::ExplicitRungeKutta<Vector<double>> fe(
@@ -288,6 +289,26 @@ main()
       TimeStepping::RK_CLASSIC_FOURTH_ORDER);
     test(rk4, f4, id_minus_tau_J_inv4, my4);
 
+    deallog << "Low-storage Runge-Kutta stage 3 order 3" << std::endl;
+    TimeStepping::LowStorageRungeKutta<Vector<double>> lsrk33(
+      TimeStepping::LOW_STORAGE_RK_STAGE3_ORDER3);
+    test(lsrk33, f3, id_minus_tau_J_inv3, my3);
+
+    deallog << "Low-storage Runge-Kutta stage 5 order 4" << std::endl;
+    TimeStepping::LowStorageRungeKutta<Vector<double>> lsrk54(
+      TimeStepping::LOW_STORAGE_RK_STAGE5_ORDER4);
+    test(lsrk54, f4, id_minus_tau_J_inv4, my4);
+
+    deallog << "Low-storage Runge-Kutta stage 7 order 4" << std::endl;
+    TimeStepping::LowStorageRungeKutta<Vector<double>> lsrk74(
+      TimeStepping::LOW_STORAGE_RK_STAGE7_ORDER4);
+    test(lsrk74, f4, id_minus_tau_J_inv4, my4);
+
+    deallog << "Low-storage Runge-Kutta stage 9 order 5" << std::endl;
+    TimeStepping::LowStorageRungeKutta<Vector<double>> lsrk95(
+      TimeStepping::LOW_STORAGE_RK_STAGE9_ORDER5);
+    test(lsrk95, f5, id_minus_tau_J_inv5, my5);
+
     deallog << "Backward Euler" << std::endl;
     TimeStepping::ImplicitRungeKutta<Vector<double>> be(
       TimeStepping::BACKWARD_EULER);
@@ -370,6 +391,38 @@ main()
                      id_minus_tau_J_inv4,
                      my_exact_solution);
 
+    deallog << "Low-storage Runge-Kutta stage 3 order 3" << std::endl;
+    TimeStepping::LowStorageRungeKutta<Vector<double>> lsrk33(
+      TimeStepping::LOW_STORAGE_RK_STAGE3_ORDER3);
+    test_convergence(lsrk33,
+                     my_rhs_function,
+                     id_minus_tau_J_inv3,
+                     my_exact_solution);
+
+    deallog << "Low-storage Runge-Kutta stage 5 order 4" << std::endl;
+    TimeStepping::LowStorageRungeKutta<Vector<double>> lsrk54(
+      TimeStepping::LOW_STORAGE_RK_STAGE5_ORDER4);
+    test_convergence(lsrk54,
+                     my_rhs_function,
+                     id_minus_tau_J_inv4,
+                     my_exact_solution);
+
+    deallog << "Low-storage Runge-Kutta stage 7 order 4" << std::endl;
+    TimeStepping::LowStorageRungeKutta<Vector<double>> lsrk74(
+      TimeStepping::LOW_STORAGE_RK_STAGE7_ORDER4);
+    test_convergence(lsrk74,
+                     my_rhs_function,
+                     id_minus_tau_J_inv4,
+                     my_exact_solution);
+
+    deallog << "Low-storage Runge-Kutta stage 9 order 5" << std::endl;
+    TimeStepping::LowStorageRungeKutta<Vector<double>> lsrk95(
+      TimeStepping::LOW_STORAGE_RK_STAGE9_ORDER5);
+    test_convergence(lsrk95,
+                     my_rhs_function,
+                     id_minus_tau_J_inv5,
+                     my_exact_solution);
+
     deallog << "Backward Euler first order" << std::endl;
     TimeStepping::ImplicitRungeKutta<Vector<double>> be(
       TimeStepping::BACKWARD_EULER);
index 0af4592f05265cbc7621dec08f905eef82b1e87f..fa4944e07776e9e036ae4dbe70b0979766ece17d 100644 (file)
@@ -7,6 +7,14 @@ DEAL::Strong Stability Preserving Runge-Kutta third order
 DEAL::0
 DEAL::Runge-Kutta fourth order
 DEAL::0
+DEAL::Low-storage Runge-Kutta stage 3 order 3
+DEAL::0
+DEAL::Low-storage Runge-Kutta stage 5 order 4
+DEAL::0
+DEAL::Low-storage Runge-Kutta stage 7 order 4
+DEAL::0
+DEAL::Low-storage Runge-Kutta stage 9 order 5
+DEAL::0
 DEAL::Backward Euler
 DEAL::0
 DEAL::Implicit midpoint
@@ -27,89 +35,97 @@ DEAL::Cash-Karp
 DEAL::0
 DEAL::Forward Euler first order
 DEAL::convergence rate
-DEAL::0.496119
 DEAL::0.634657
 DEAL::0.763742
 DEAL::0.862050
 DEAL::0.924783
 DEAL::0.960617
 DEAL::0.979834
-DEAL::0.989794
-DEAL::0.994866
 DEAL::Runge-Kutta third order
 DEAL::convergence rate
-DEAL::7.01207
 DEAL::0.0523936
 DEAL::2.00049
 DEAL::2.55207
 DEAL::2.78347
 DEAL::2.89332
 DEAL::2.94704
-DEAL::2.97361
-DEAL::2.98681
 DEAL::Strong Stability Preserving Runge-Kutta third order
 DEAL::convergence rate
-DEAL::2.44463
 DEAL::2.72007
 DEAL::2.86894
 DEAL::2.93941
 DEAL::2.97135
 DEAL::2.98614
 DEAL::2.99319
-DEAL::2.99663
-DEAL::2.99832
 DEAL::Runge-Kutta fourth order
 DEAL::convergence rate
-DEAL::3.32883
 DEAL::3.64708
 DEAL::3.86926
 DEAL::3.95478
 DEAL::3.98399
 DEAL::3.99389
-DEAL::3.99745
-DEAL::3.99887
-DEAL::3.99923
+DEAL::3.99746
+DEAL::Low-storage Runge-Kutta stage 3 order 3
+DEAL::convergence rate
+DEAL::2.65986
+DEAL::2.82914
+DEAL::2.91745
+DEAL::2.95989
+DEAL::2.98029
+DEAL::2.99024
+DEAL::Low-storage Runge-Kutta stage 5 order 4
+DEAL::convergence rate
+DEAL::3.52474
+DEAL::3.77867
+DEAL::3.89561
+DEAL::3.94964
+DEAL::3.97531
+DEAL::3.98778
+DEAL::Low-storage Runge-Kutta stage 7 order 4
+DEAL::convergence rate
+DEAL::-0.579086
+DEAL::3.29557
+DEAL::3.71647
+DEAL::3.87052
+DEAL::3.93798
+DEAL::3.96951
+DEAL::Low-storage Runge-Kutta stage 9 order 5
+DEAL::convergence rate
+DEAL::4.66916
+DEAL::4.85354
+DEAL::4.93262
+DEAL::4.96796
+DEAL::4.98436
+DEAL::4.99866
 DEAL::Backward Euler first order
 DEAL::convergence rate
-DEAL::94.3469
-DEAL::6.54344
+DEAL::6.54345
 DEAL::1.54696
-DEAL::1.20584
-DEAL::1.09169
-DEAL::1.04347
-DEAL::1.02119
-DEAL::1.01048
-DEAL::1.00519
+DEAL::1.20585
+DEAL::1.09170
+DEAL::1.04350
+DEAL::1.02129
 DEAL::Implicit midpoint second order
 DEAL::convergence rate
-DEAL::1.35296
-DEAL::1.96494
-DEAL::1.99730
-DEAL::1.99994
-DEAL::2.00083
-DEAL::2.00283
-DEAL::2.01821
-DEAL::2.00827
-DEAL::2.18675
+DEAL::1.96525
+DEAL::1.99777
+DEAL::2.00269
+DEAL::2.00884
+DEAL::2.04892
+DEAL::2.05461
 DEAL::Crank-Nicolson second order
 DEAL::convergence rate
-DEAL::7.33468
-DEAL::2.34239
-DEAL::2.07462
-DEAL::2.01811
-DEAL::2.00462
-DEAL::2.00144
-DEAL::2.00230
-DEAL::2.00201
-DEAL::2.01221
+DEAL::2.34241
+DEAL::2.07476
+DEAL::2.01846
+DEAL::2.00559
+DEAL::2.00652
+DEAL::2.00632
 DEAL::SDIRK second order
 DEAL::convergence rate
-DEAL::2.55970
-DEAL::2.11577
-DEAL::2.02903
-DEAL::2.00817
-DEAL::2.00306
-DEAL::2.00408
-DEAL::2.00486
-DEAL::2.02211
-DEAL::2.24289
+DEAL::2.11605
+DEAL::2.02959
+DEAL::2.00992
+DEAL::2.01080
+DEAL::2.01420
+DEAL::2.11290
diff --git a/tests/base/time_stepping_03.cc b/tests/base/time_stepping_03.cc
new file mode 100644 (file)
index 0000000..d974b0d
--- /dev/null
@@ -0,0 +1,61 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// test the coefficients of the low-storage Runge-Kutta methods in TimeStepping
+#include <deal.II/base/time_stepping.h>
+
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  const std::vector<double> bi = {{1153189308089. / 22510343858157.,
+                                   1772645290293. / 4653164025191.,
+                                   -1672844663538. / 4480602732383.,
+                                   2114624349019. / 3568978502595.,
+                                   5198255086312. / 14908931495163.}};
+  std::vector<double>       ai;
+  ai = {{970286171893. / 4311952581923.,
+         6584761158862. / 12103376702013.,
+         2251764453980. / 15575788980749.,
+         26877169314380. / 34165994151039.}};
+
+  deallog << "Check low-storage Runge-Kutta coefficients" << std::endl;
+  TimeStepping::LowStorageRungeKutta<Vector<double>> lsrk54(
+    TimeStepping::LOW_STORAGE_RK_STAGE5_ORDER4);
+  std::vector<double> a, b, c;
+  lsrk54.get_coefficients(a, b, c);
+
+  double sum_previous_bi = 0.;
+  for (unsigned int i = 0; i < b.size(); ++i)
+    {
+      if (i != b.size() - 1)
+        Assert(std::fabs(a[i] - ai[i]) < 1e-10, ExcInternalError());
+      Assert(std::fabs(b[i] - bi[i]) < 1e-10, ExcInternalError());
+
+      if (i > 0)
+        {
+          const double ci = sum_previous_bi + ai[i - 1];
+          deallog << " c: " << c[i] - ci << std::endl;
+          sum_previous_bi += bi[i - 1];
+        }
+      else
+        deallog << " c: " << c[0] << std::endl;
+    }
+
+  return 0;
+}
diff --git a/tests/base/time_stepping_03.output b/tests/base/time_stepping_03.output
new file mode 100644 (file)
index 0000000..0b32627
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::Check low-storage Runge-Kutta coefficients
+DEAL:: c: 0
+DEAL:: c: 0
+DEAL:: c: 0
+DEAL:: c: 0
+DEAL:: c: 0
\ No newline at end of file

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.