add ssp rk 3

diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib
index cbcb203815e4890cbc83ce378c110329930c1700..3583362289c0b5f94f89d97e4ec772c04c8f16eb 100644 (file)
--- a/doc/doxygen/references.bib
+++ b/doc/doxygen/references.bib
@@ -969,3 +969,15 @@ year = {2008},
   publisher =    {Microsoft Press},
   year =         2004,
   edition =   {second}}
+
+@article{gottlieb2001strong,
+  title={Strong stability-preserving high-order time discretization methods},
+  author={Gottlieb, Sigal and Shu, Chi-Wang and Tadmor, Eitan},
+  journal={SIAM review},
+  volume={43},
+  number={1},
+  pages={89--112},
+  year={2001},
+  publisher={SIAM}
+}
+

diff --git a/doc/news/changes/minor/20201005JiaqiZhang b/doc/news/changes/minor/20201005JiaqiZhang
new file mode 100644 (file)
index 0000000..eaead32
--- /dev/null
+++ b/doc/news/changes/minor/20201005JiaqiZhang
@@ -0,0 +1,5 @@
+New: SSP_THIRD_ORDER is added to the namespace TimeStepping to 
+implement the explicit third order Strong Stability Preserving (SSP) Runge-Kutta method, 
+which is also called the third order Total Variation Diminishing (TVD) Runge-Kutta method, see @cite gottlieb2001strong.
+<br>
+(Jiaqi Zhang, 2020/10/05)

diff --git a/include/deal.II/base/time_stepping.h b/include/deal.II/base/time_stepping.h
index c0a5096cd4e542fa7a4be6214d0c0c70bb17059c..5b1f418c28531c35f00872d9b06baab32141660a 100644 (file)
--- a/include/deal.II/base/time_stepping.h
+++ b/include/deal.II/base/time_stepping.h
@@ -37,6 +37,7 @@ namespace TimeStepping
    * - Explicit methods (see ExplicitRungeKutta::initialize):
    *   - FORWARD_EULER (first order)
    *   - RK_THIRD_ORDER (third order Runge-Kutta)
+   *   - SSP_THIRD_ORDER (third order SSP Runge-Kutta)
    *   - RK_CLASSIC_FOURTH_ORDER (classical fourth order Runge-Kutta)
    * - Implicit methods (see ImplicitRungeKutta::initialize):
    *   - BACKWARD_EULER (first order)
@@ -61,6 +62,12 @@ namespace TimeStepping
      * Third order Runge-Kutta method.
      */
     RK_THIRD_ORDER,
+    /**
+     * Third order Strong Stability Preserving (SSP) Runge-Kutta method
+     * (SSP time discretizations are also called Total Variation Diminishing
+     * (TVD) methods in the literature, see @cite gottlieb2001strong).
+     */
+    SSP_THIRD_ORDER,
     /**
      * Classical fourth order Runge-Kutta method.
      */
@@ -514,7 +521,7 @@ namespace TimeStepping
 
 
   /**
-   * This is class is derived from RungeKutta and implement embedded explicit
+   * This class is derived from RungeKutta and implements embedded explicit
    * methods.
    */
   template <typename VectorType>

diff --git a/include/deal.II/base/time_stepping.templates.h b/include/deal.II/base/time_stepping.templates.h
index 563ab6154d0d1f537c6e9e76e249e2b8ced567cb..7a363a756225ae9f8208b2a6b676294b36be8b2f 100644 (file)
--- a/include/deal.II/base/time_stepping.templates.h
+++ b/include/deal.II/base/time_stepping.templates.h
@@ -110,6 +110,29 @@ namespace TimeStepping
             tmp[1] = 2.0;
             this->a.push_back(tmp);
 
+            break;
+          }
+        case (SSP_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(1.0 / 6.0);
+            this->b.push_back(2.0 / 3.0);
+            this->c.push_back(0.0);
+            this->c.push_back(1.0);
+            this->c.push_back(0.5);
+            std::vector<double> tmp;
+            this->a.push_back(tmp);
+            tmp.resize(1);
+            tmp[0] = 1.0;
+            this->a.push_back(tmp);
+            tmp.resize(2);
+            tmp[0] = 1.0 / 4.0;
+            tmp[1] = 1.0 / 4.0;
+            this->a.push_back(tmp);
+
             break;
           }
         case (RK_CLASSIC_FOURTH_ORDER):

diff --git a/tests/base/time_stepping_01.cc b/tests/base/time_stepping_01.cc
index 3f48f8e1c399a802bcff88ac232553f2ecec7042..a7d7b51dba985394ece5d1ffeb72b991ee9480ce 100644 (file)
--- a/tests/base/time_stepping_01.cc
+++ b/tests/base/time_stepping_01.cc
@@ -13,14 +13,14 @@
 //
 // ---------------------------------------------------------------------
 
-
+// test Runge-Kutta methods in TimeStepping with a) a polynomial with expected
+// error 0 and b) convergence order for y=exp(t^2)
 #include <deal.II/base/time_stepping.h>
 
 #include <deal.II/lac/vector.h>
 
 #include "../tests.h"
 
-// test Runge-Kutta methods
 Vector<double>
 f1(double const t, Vector<double> const &y)
 {
@@ -71,6 +71,17 @@ f5(double const t, Vector<double> const &y)
   return values;
 }
 
+
+Vector<double>
+my_rhs_function(double const t, Vector<double> const &y)
+{
+  Vector<double> values(y);
+  for (unsigned int i = 0; i < values.size(); ++i)
+    values[i] = y[i] * 2.0 * t;
+
+  return values;
+}
+
 Vector<double>
 id_minus_tau_J_inv1(double const t, double const tau, Vector<double> const &y)
 {
@@ -131,6 +142,12 @@ my5(double const t)
   return t * t * t * t * t;
 }
 
+double
+my_exact_solution(double const t)
+{
+  return std::exp(t * t);
+}
+
 void
 test(TimeStepping::RungeKutta<Vector<double>> &                          solver,
      std::function<Vector<double>(double const, Vector<double> const &)> f,
@@ -200,72 +217,191 @@ test2(TimeStepping::EmbeddedExplicitRungeKutta<Vector<double>> &solver,
   deallog << error_norm << std::endl;
 }
 
+
+void
+test_convergence(
+  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)
+{
+  std::vector<double> errors;
+  double              initial_time = 0.0, final_time = 1.0;
+  unsigned int        size = 1;
+  Vector<double>      solution(size);
+  Vector<double>      exact_solution(size);
+  for (unsigned int i = 0; i < size; ++i)
+    {
+      exact_solution[i] = my(final_time);
+    }
+
+  deallog << "convergence rate" << std::endl;
+  for (unsigned int cycle = 0; cycle < 10; ++cycle)
+    {
+      unsigned int n_time_steps = std::pow(2., static_cast<double>(cycle));
+      double       time_step =
+        (final_time - initial_time) / static_cast<double>(n_time_steps);
+      double time = initial_time;
+      for (unsigned int i = 0; i < size; ++i)
+        solution[i] = my(initial_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();
+      errors.push_back(error_norm);
+      if (cycle > 0)
+        deallog << std::log(std::fabs(errors[cycle - 1] / errors[cycle])) /
+                     std::log(2.)
+                << std::endl;
+    }
+}
+
 int
 main()
 {
   initlog();
-
-  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);
+  {
+    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 << "Strong Stability Preserving Runge-Kutta third order"
+            << std::endl;
+    TimeStepping::ExplicitRungeKutta<Vector<double>> ssp_rk3(
+      TimeStepping::SSP_THIRD_ORDER);
+    test(ssp_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);
+  }
+
+  {
+    deallog << "Forward Euler first order" << std::endl;
+    TimeStepping::ExplicitRungeKutta<Vector<double>> rk1(
+      TimeStepping::FORWARD_EULER);
+    test_convergence(rk1,
+                     my_rhs_function,
+                     id_minus_tau_J_inv1,
+                     my_exact_solution);
+
+    deallog << "Runge-Kutta third order" << std::endl;
+    TimeStepping::ExplicitRungeKutta<Vector<double>> rk3(
+      TimeStepping::RK_THIRD_ORDER);
+    test_convergence(rk3,
+                     my_rhs_function,
+                     id_minus_tau_J_inv3,
+                     my_exact_solution);
+
+    deallog << "Strong Stability Preserving Runge-Kutta third order"
+            << std::endl;
+    TimeStepping::ExplicitRungeKutta<Vector<double>> ssp_rk3(
+      TimeStepping::SSP_THIRD_ORDER);
+    test_convergence(ssp_rk3,
+                     my_rhs_function,
+                     id_minus_tau_J_inv3,
+                     my_exact_solution);
+
+    deallog << "Runge-Kutta fourth order" << std::endl;
+    TimeStepping::ExplicitRungeKutta<Vector<double>> rk4(
+      TimeStepping::RK_CLASSIC_FOURTH_ORDER);
+    test_convergence(rk4,
+                     my_rhs_function,
+                     id_minus_tau_J_inv4,
+                     my_exact_solution);
+
+    deallog << "Backward Euler first order" << std::endl;
+    TimeStepping::ImplicitRungeKutta<Vector<double>> be(
+      TimeStepping::BACKWARD_EULER);
+    test_convergence(be,
+                     my_rhs_function,
+                     id_minus_tau_J_inv1,
+                     my_exact_solution);
+
+    deallog << "Implicit midpoint second order" << std::endl;
+    TimeStepping::ImplicitRungeKutta<Vector<double>> im(
+      TimeStepping::IMPLICIT_MIDPOINT);
+    test_convergence(im,
+                     my_rhs_function,
+                     id_minus_tau_J_inv2,
+                     my_exact_solution);
+
+    deallog << "Crank-Nicolson second order" << std::endl;
+    TimeStepping::ImplicitRungeKutta<Vector<double>> cn(
+      TimeStepping::CRANK_NICOLSON);
+    test_convergence(cn,
+                     my_rhs_function,
+                     id_minus_tau_J_inv2,
+                     my_exact_solution);
+
+    deallog << "SDIRK second order" << std::endl;
+    TimeStepping::ImplicitRungeKutta<Vector<double>> sdirk(
+      TimeStepping::SDIRK_TWO_STAGES);
+    test_convergence(sdirk,
+                     my_rhs_function,
+                     id_minus_tau_J_inv2,
+                     my_exact_solution);
+  }
 
   return 0;
 }

diff --git a/tests/base/time_stepping_01.output b/tests/base/time_stepping_01.output
index 78c9c17007bac49c8e514184e10d1ec78cac3510..0af4592f05265cbc7621dec08f905eef82b1e87f 100644 (file)
--- a/tests/base/time_stepping_01.output
+++ b/tests/base/time_stepping_01.output
@@ -3,6 +3,8 @@ DEAL::Forward Euler
 DEAL::0
 DEAL::Runge-Kutta third order
 DEAL::0
+DEAL::Strong Stability Preserving Runge-Kutta third order
+DEAL::0
 DEAL::Runge-Kutta fourth order
 DEAL::0
 DEAL::Backward Euler
@@ -23,3 +25,91 @@ DEAL::Fehlberg
 DEAL::0
 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::Backward Euler first order
+DEAL::convergence rate
+DEAL::94.3469
+DEAL::6.54344
+DEAL::1.54696
+DEAL::1.20584
+DEAL::1.09169
+DEAL::1.04347
+DEAL::1.02119
+DEAL::1.01048
+DEAL::1.00519
+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::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::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