From: Jiaqi Zhang Date: Mon, 30 Nov 2020 18:14:07 +0000 (-0500) Subject: add low-storage rk and update step-67 X-Git-Tag: v9.3.0-rc1~765^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dcf516a33d9b71fd9dee53cacf1051e8258d823c;p=dealii.git add low-storage rk and update step-67 --- diff --git a/doc/news/changes/minor/20201118JiaqiZhang b/doc/news/changes/minor/20201118JiaqiZhang new file mode 100644 index 0000000000..0054c73138 --- /dev/null +++ b/doc/news/changes/minor/20201118JiaqiZhang @@ -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. +
+(Jiaqi Zhang, 2020/11/18) diff --git a/examples/step-67/step-67.cc b/examples/step-67/step-67.cc index d1685c5bea..9186562b23 100644 --- a/examples/step-67/step-67.cc +++ b/examples/step-67/step-67.cc @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -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> + 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 bi; std::vector ai; + std::vector ci; }; diff --git a/include/deal.II/base/time_stepping.h b/include/deal.II/base/time_stepping.h index 5b1f418c28..8cc69fcb91 100644 --- a/include/deal.II/base/time_stepping.h +++ b/include/deal.II/base/time_stepping.h @@ -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 + class LowStorageRungeKutta : public RungeKutta + { + public: + using RungeKutta::evolve_one_time_step; + + /** + * Default constructor. This constructor creates an object for which + * you will want to call initialize(runge_kutta_method) + * 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 &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 &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 &a, + std::vector &b, + std::vector &c) const; + + /** + * This structure stores the name of the method used. + */ + struct Status : public TimeStepping::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 &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. diff --git a/include/deal.II/base/time_stepping.templates.h b/include/deal.II/base/time_stepping.templates.h index 7a363a7562..1998d005d3 100644 --- a/include/deal.II/base/time_stepping.templates.h +++ b/include/deal.II/base/time_stepping.templates.h @@ -241,6 +241,234 @@ namespace TimeStepping + // ---------------------------------------------------------------------- + // LowStorageRungeKutta + // ---------------------------------------------------------------------- + + template + LowStorageRungeKutta::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::initialize(method); + } + + + + template + void + LowStorageRungeKutta::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 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 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 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 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 + double + LowStorageRungeKutta::evolve_one_time_step( + const std::function &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 + double + LowStorageRungeKutta::evolve_one_time_step( + const std::function &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 + void + LowStorageRungeKutta::get_coefficients( + std::vector &a, + std::vector &b, + std::vector &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 + const typename LowStorageRungeKutta::Status & + LowStorageRungeKutta::get_status() const + { + return status; + } + + template + void + LowStorageRungeKutta::compute_one_stage( + const std::function &f, + const double t, + const double factor_solution, + const double factor_ai, + const VectorType ¤t_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 // ---------------------------------------------------------------------- diff --git a/source/base/time_stepping.inst.in b/source/base/time_stepping.inst.in index 6315226010..e9667475a7 100644 --- a/source/base/time_stepping.inst.in +++ b/source/base/time_stepping.inst.in @@ -18,6 +18,7 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) { template class RungeKutta>; template class ExplicitRungeKutta>; + template class LowStorageRungeKutta>; template class ImplicitRungeKutta>; template class EmbeddedExplicitRungeKutta>; } @@ -26,6 +27,7 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) { template class RungeKutta>; template class ExplicitRungeKutta>; + template class LowStorageRungeKutta>; template class ImplicitRungeKutta>; template class EmbeddedExplicitRungeKutta>; } @@ -34,6 +36,7 @@ for (V : EXTERNAL_PARALLEL_VECTORS) { template class RungeKutta; template class ExplicitRungeKutta; + template class LowStorageRungeKutta; template class ImplicitRungeKutta; template class EmbeddedExplicitRungeKutta; } diff --git a/tests/base/time_stepping_01.cc b/tests/base/time_stepping_01.cc index a7d7b51dba..0fd9bdf480 100644 --- a/tests/base/time_stepping_01.cc +++ b/tests/base/time_stepping_01.cc @@ -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 #include @@ -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(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> 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> 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> 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> 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> lsrk95( + TimeStepping::LOW_STORAGE_RK_STAGE9_ORDER5); + test(lsrk95, f5, id_minus_tau_J_inv5, my5); + deallog << "Backward Euler" << std::endl; TimeStepping::ImplicitRungeKutta> 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> 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> 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> 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> 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> be( TimeStepping::BACKWARD_EULER); diff --git a/tests/base/time_stepping_01.output b/tests/base/time_stepping_01.output index 0af4592f05..fa4944e077 100644 --- a/tests/base/time_stepping_01.output +++ b/tests/base/time_stepping_01.output @@ -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 index 0000000000..d974b0d922 --- /dev/null +++ b/tests/base/time_stepping_03.cc @@ -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 + +#include "../tests.h" + +int +main() +{ + initlog(); + + const std::vector bi = {{1153189308089. / 22510343858157., + 1772645290293. / 4653164025191., + -1672844663538. / 4480602732383., + 2114624349019. / 3568978502595., + 5198255086312. / 14908931495163.}}; + std::vector ai; + ai = {{970286171893. / 4311952581923., + 6584761158862. / 12103376702013., + 2251764453980. / 15575788980749., + 26877169314380. / 34165994151039.}}; + + deallog << "Check low-storage Runge-Kutta coefficients" << std::endl; + TimeStepping::LowStorageRungeKutta> lsrk54( + TimeStepping::LOW_STORAGE_RK_STAGE5_ORDER4); + std::vector 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 index 0000000000..0b32627dc7 --- /dev/null +++ b/tests/base/time_stepping_03.output @@ -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