From: Jiaqi Zhang Date: Mon, 5 Oct 2020 15:06:00 +0000 (+0000) Subject: add ssp rk 3 X-Git-Tag: v9.3.0-rc1~1031^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11009%2Fhead;p=dealii.git add ssp rk 3 --- diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index cbcb203815..3583362289 100644 --- 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 index 0000000000..eaead32090 --- /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. +
+(Jiaqi Zhang, 2020/10/05) diff --git a/include/deal.II/base/time_stepping.h b/include/deal.II/base/time_stepping.h index c0a5096cd4..5b1f418c28 100644 --- 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 diff --git a/include/deal.II/base/time_stepping.templates.h b/include/deal.II/base/time_stepping.templates.h index 563ab6154d..7a363a7562 100644 --- 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 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 3f48f8e1c3..a7d7b51dba 100644 --- 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 #include #include "../tests.h" -// test Runge-Kutta methods Vector f1(double const t, Vector const &y) { @@ -71,6 +71,17 @@ f5(double const t, Vector const &y) return values; } + +Vector +my_rhs_function(double const t, Vector const &y) +{ + Vector values(y); + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = y[i] * 2.0 * t; + + return values; +} + Vector id_minus_tau_J_inv1(double const t, double const tau, Vector 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> & solver, std::function(double const, Vector const &)> f, @@ -200,72 +217,191 @@ test2(TimeStepping::EmbeddedExplicitRungeKutta> &solver, deallog << error_norm << std::endl; } + +void +test_convergence( + 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) +{ + std::vector errors; + double initial_time = 0.0, final_time = 1.0; + unsigned int size = 1; + Vector solution(size); + Vector 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(cycle)); + double time_step = + (final_time - initial_time) / static_cast(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 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> fe( - TimeStepping::FORWARD_EULER); - test(fe, f1, id_minus_tau_J_inv1, my1); - - deallog << "Runge-Kutta third order" << std::endl; - TimeStepping::ExplicitRungeKutta> rk3( - TimeStepping::RK_THIRD_ORDER); - test(rk3, f3, id_minus_tau_J_inv3, my3); - - deallog << "Runge-Kutta fourth order" << std::endl; - TimeStepping::ExplicitRungeKutta> rk4( - TimeStepping::RK_CLASSIC_FOURTH_ORDER); - test(rk4, f4, id_minus_tau_J_inv4, my4); - - deallog << "Backward Euler" << std::endl; - TimeStepping::ImplicitRungeKutta> be( - TimeStepping::BACKWARD_EULER); - test(be, f1, id_minus_tau_J_inv1, my1); - - deallog << "Implicit midpoint" << std::endl; - TimeStepping::ImplicitRungeKutta> im( - TimeStepping::IMPLICIT_MIDPOINT); - test(im, f2, id_minus_tau_J_inv2, my2); - - deallog << "Crank-Nicolson" << std::endl; - TimeStepping::ImplicitRungeKutta> cn( - TimeStepping::CRANK_NICOLSON); - test(cn, f2, id_minus_tau_J_inv2, my2); - - deallog << "SDIRK" << std::endl; - TimeStepping::ImplicitRungeKutta> sdirk( - TimeStepping::SDIRK_TWO_STAGES); - test(sdirk, f2, id_minus_tau_J_inv2, my2); - - deallog << "Heun-Euler" << std::endl; - TimeStepping::EmbeddedExplicitRungeKutta> he( - TimeStepping::HEUN_EULER); - test2(he, f2, id_minus_tau_J_inv2, my2); - - deallog << "Bogacki-Shampine" << std::endl; - TimeStepping::EmbeddedExplicitRungeKutta> bs( - TimeStepping::BOGACKI_SHAMPINE); - test2(bs, f3, id_minus_tau_J_inv3, my3); - bs.free_memory(); - - deallog << "DOPRI" << std::endl; - TimeStepping::EmbeddedExplicitRungeKutta> dopri( - TimeStepping::DOPRI); - test2(dopri, f5, id_minus_tau_J_inv5, my5); - dopri.free_memory(); - - deallog << "Fehlberg" << std::endl; - TimeStepping::EmbeddedExplicitRungeKutta> fehlberg( - TimeStepping::FEHLBERG); - test2(fehlberg, f5, id_minus_tau_J_inv5, my5); - - deallog << "Cash-Karp" << std::endl; - TimeStepping::EmbeddedExplicitRungeKutta> ck( - TimeStepping::CASH_KARP); - test2(ck, f5, id_minus_tau_J_inv5, my5); + { + deallog << "Forward Euler" << std::endl; + TimeStepping::ExplicitRungeKutta> fe( + TimeStepping::FORWARD_EULER); + test(fe, f1, id_minus_tau_J_inv1, my1); + + deallog << "Runge-Kutta third order" << std::endl; + TimeStepping::ExplicitRungeKutta> 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> 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> rk4( + TimeStepping::RK_CLASSIC_FOURTH_ORDER); + test(rk4, f4, id_minus_tau_J_inv4, my4); + + deallog << "Backward Euler" << std::endl; + TimeStepping::ImplicitRungeKutta> be( + TimeStepping::BACKWARD_EULER); + test(be, f1, id_minus_tau_J_inv1, my1); + + deallog << "Implicit midpoint" << std::endl; + TimeStepping::ImplicitRungeKutta> im( + TimeStepping::IMPLICIT_MIDPOINT); + test(im, f2, id_minus_tau_J_inv2, my2); + + deallog << "Crank-Nicolson" << std::endl; + TimeStepping::ImplicitRungeKutta> cn( + TimeStepping::CRANK_NICOLSON); + test(cn, f2, id_minus_tau_J_inv2, my2); + + deallog << "SDIRK" << std::endl; + TimeStepping::ImplicitRungeKutta> sdirk( + TimeStepping::SDIRK_TWO_STAGES); + test(sdirk, f2, id_minus_tau_J_inv2, my2); + + deallog << "Heun-Euler" << std::endl; + TimeStepping::EmbeddedExplicitRungeKutta> he( + TimeStepping::HEUN_EULER); + test2(he, f2, id_minus_tau_J_inv2, my2); + + deallog << "Bogacki-Shampine" << std::endl; + TimeStepping::EmbeddedExplicitRungeKutta> bs( + TimeStepping::BOGACKI_SHAMPINE); + test2(bs, f3, id_minus_tau_J_inv3, my3); + bs.free_memory(); + + deallog << "DOPRI" << std::endl; + TimeStepping::EmbeddedExplicitRungeKutta> dopri( + TimeStepping::DOPRI); + test2(dopri, f5, id_minus_tau_J_inv5, my5); + dopri.free_memory(); + + deallog << "Fehlberg" << std::endl; + TimeStepping::EmbeddedExplicitRungeKutta> fehlberg( + TimeStepping::FEHLBERG); + test2(fehlberg, f5, id_minus_tau_J_inv5, my5); + + deallog << "Cash-Karp" << std::endl; + TimeStepping::EmbeddedExplicitRungeKutta> ck( + TimeStepping::CASH_KARP); + test2(ck, f5, id_minus_tau_J_inv5, my5); + } + + { + deallog << "Forward Euler first order" << std::endl; + TimeStepping::ExplicitRungeKutta> 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> 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> 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> 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> 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> 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> 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> 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 78c9c17007..0af4592f05 100644 --- 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