--- /dev/null
+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)
#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>
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
{
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;
}
// 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;
}
// 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;
}
// 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
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 ?
vec_ki,
solution,
vec_ri);
- sum_previous_bi += bi[stage - 1];
}
}
private:
std::vector<double> bi;
std::vector<double> ai;
+ std::vector<double> ci;
};
* - 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)
* 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.
*/
+ /**
+ * 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.
+ // ----------------------------------------------------------------------
+ // 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 ¤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
// ----------------------------------------------------------------------
{
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>>;
}
{
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>>;
}
{
template class RungeKutta<V>;
template class ExplicitRungeKutta<V>;
+ template class LowStorageRungeKutta<V>;
template class ImplicitRungeKutta<V>;
template class EmbeddedExplicitRungeKutta<V>;
}
// ---------------------------------------------------------------------
// 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>
double
my_exact_solution(double const t)
{
- return std::exp(t * t);
+ return 0.1 * std::exp(t * t);
}
void
}
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 =
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;
main()
{
initlog();
+ // deallog.precision(4);
{
deallog << "Forward Euler" << std::endl;
TimeStepping::ExplicitRungeKutta<Vector<double>> fe(
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);
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);
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
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff