timesteppers.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// This header contains classes and function prototypes for the Time,
// TimeStepper and derived objects
 
// Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_TIME_STEPPERS_HEADER
#define OOMPH_TIME_STEPPERS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
// oomph-lib headers
#include "Vector.h"
#include "nodes.h"
#include "matrices.h"
 
namespace oomph
{
  class Problem;
  class ExplicitTimeStepper;
 
  //====================================================================
  //====================================================================
  class Time
  {
  private:
    double Continuous_time;
 
    Vector<double> Dt;
 
  public:
    Time() : Continuous_time(0.0) {}
 
    Time(const unsigned& ndt) : Continuous_time(0.0)
    {
      // Allocate memory for the timestep storage and initialise values to one
      // to avoid division by zero.
      Dt.resize(ndt, 1.0);
    }
 
    Time(const Time&) = delete;
 
    void operator=(const Time&) = delete;
 
    void resize(const unsigned& n_dt)
    {
      Dt.resize(n_dt, 0.0);
    }
 
    void initialise_dt(const double& dt_)
    {
      unsigned ndt = Dt.size();
      Dt.assign(ndt, dt_);
    }
 
    void initialise_dt(const Vector<double>& dt_)
    {
      // Assign the values from the vector dt_, but preserve the size of Dt and
      // any Dt[i] s.t. i > dt_.size() (which is why we can't just use
      // assignment operator).
      unsigned n_dt = dt_.size();
      for (unsigned i = 0; i < n_dt; i++)
      {
        Dt[i] = dt_[i];
      }
    }
 
    ~Time() {}
 
    double& time()
    {
      return Continuous_time;
    }
 
    unsigned ndt() const
    {
      return Dt.size();
    }
 
    double& dt(const unsigned& t = 0)
    {
      return Dt[t];
    }
 
    double dt(const unsigned& t = 0) const
    {
      return Dt[t];
    }
 
    double time(const unsigned& t = 0) const
    {
#ifdef PARANOID
      if (t > ndt())
      {
        using namespace StringConversion;
        std::string error_msg = "Timestep " + to_string(t) + " out of range!";
        throw OomphLibError(
          error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Load the current value of the time
      double time_local = Continuous_time;
      // Loop over the t previous timesteps and subtract each dt
      for (unsigned i = 0; i < t; i++)
      {
        time_local -= Dt[i];
      }
      return time_local;
    }
 
    void shift_dt()
    {
      unsigned n_dt = Dt.size();
      // Return straight away if there are no stored timesteps
      if (n_dt == 0)
      {
        return;
      }
      // Start from the end of the array and shift every entry back by one
      for (unsigned i = (n_dt - 1); i > 0; i--)
      {
        Dt[i] = Dt[i - 1];
      }
    }
  };
 
 
 
 
  //====================================================================
  //====================================================================
  class TimeStepper
  {
  protected:
    Time* Time_pt;
 
    DenseMatrix<double> Weight;
 
    std::string Type;
 
    bool Adaptive_Flag;
 
    bool Is_steady;
 
    bool Shut_up_in_assign_initial_data_values;
 
    bool Predict_by_explicit_step;
 
    ExplicitTimeStepper* Explicit_predictor_pt;
 
    double Predicted_time;
 
    int Predictor_storage_index;
 
  public:
    TimeStepper(const unsigned& tstorage, const unsigned& max_deriv)
      : Time_pt(0),
        Adaptive_Flag(false),
        Is_steady(false),
        Shut_up_in_assign_initial_data_values(false),
        Predict_by_explicit_step(false)
    {
      // Resize Weights matrix and initialise each weight to zero
      Weight.resize(max_deriv + 1, tstorage, 0.0);
 
      // Set the weight for zero-th derivative, which is always 1.0
      Weight(0, 0) = 1.0;
 
      // Set predictor storage index to negative value so that we can catch
      // cases where it has not been set to a correct value by the inheriting
      // constructor.
      Predictor_storage_index = -1;
 
      Explicit_predictor_pt = 0;
    }
 
    TimeStepper()
    {
      throw OomphLibError("Don't call default constructor for TimeStepper!",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    TimeStepper(const TimeStepper&) = delete;
 
    void operator=(const TimeStepper&) = delete;
 
    virtual ~TimeStepper();
 
    unsigned highest_derivative() const
    {
      return Weight.nrow() - 1;
    }
 
    virtual unsigned order() const
    {
      return 0;
    }
 
    // (can't have a paranoid test for null pointers because this could be
    // used as a set function)
    double& time()
    {
      return Time_pt->time();
    }
 
    double time() const
    {
#ifdef PARANOID
      if (Time_pt == 0)
      {
        std::string error_msg = "Time pointer is null!";
        throw OomphLibError(
          error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Time_pt->time();
    }
 
    virtual unsigned ndt() const = 0;
 
    virtual unsigned nprev_values_for_value_at_evaluation_time() const
    {
      return 1;
    }
 
    virtual unsigned nprev_values() const = 0;
 
    virtual void set_weights() = 0;
 
    void make_steady()
    {
      // Zero the matrix
      Weight.initialise(0);
 
      // Weight of current step in calculation of current step is always 1 in
      // steady state. All other entries left at zero.
      Weight(0, 0) = 1.0;
 
      // Update flag
      Is_steady = true;
    }
 
    bool is_steady() const
    {
      return Is_steady;
    }
 
    bool predict_by_explicit_step() const
    {
      return Predict_by_explicit_step;
    }
 
    ExplicitTimeStepper* explicit_predictor_pt()
    {
      return Explicit_predictor_pt;
    }
 
    void set_predictor_pt(ExplicitTimeStepper* _pred_pt)
    {
      Explicit_predictor_pt = _pred_pt;
    }
 
    void update_predicted_time(const double& new_time)
    {
      Predicted_time = new_time;
    }
 
    void check_predicted_values_up_to_date() const
    {
#ifdef PARANOID
      if (std::abs(time_pt()->time() - Predicted_time) > 1e-15)
      {
        throw OomphLibError(
          "Predicted values are not from the correct time step",
          OOMPH_EXCEPTION_LOCATION,
          OOMPH_CURRENT_FUNCTION);
      }
#endif
    }
 
    unsigned predictor_storage_index() const
    {
      // Error if time stepper is non-adaptive (because then the index doesn't
      // exist so using it would give a potentially difficult to find
      // segault).
#ifdef PARANOID
      if (Predictor_storage_index > 0)
      {
        return unsigned(Predictor_storage_index);
      }
      else
      {
        std::string err = "Predictor storage index is negative, this probably";
        err += " means it hasn't been set for this timestepper.";
        throw OomphLibError(
          err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
      }
#else
      return unsigned(Predictor_storage_index);
#endif
    }
 
    void enable_warning_in_assign_initial_data_values()
    {
      Shut_up_in_assign_initial_data_values = false;
    }
 
    void disable_warning_in_assign_initial_data_values()
    {
      Shut_up_in_assign_initial_data_values = true;
    }
 
    const DenseMatrix<double>* weights_pt() const
    {
      return &Weight;
    }
 
    virtual void undo_make_steady()
    {
      Is_steady = false;
      set_weights();
    }
 
    std::string type() const
    {
      return Type;
    }
 
    //??ds I think that the Data and Node cases below couldp robably be
    // handled with a template argument. However they can't be handled by
    // normal polymorphism because data.value is different to node.value and
    // value is not a virtual function.
 
    void time_derivative(const unsigned& i,
                         Data* const& data_pt,
                         Vector<double>& deriv)
    {
      // Number of values stored in the Data object
      unsigned nvalue = data_pt->nvalue();
      deriv.assign(nvalue, 0.0);
 
      // Loop over all values
      for (unsigned j = 0; j < nvalue; j++)
      {
        deriv[j] = time_derivative(i, data_pt, j);
      }
    }
 
    double time_derivative(const unsigned& i,
                           Data* const& data_pt,
                           const unsigned& j)
    {
      double deriv = 0.0;
      unsigned n_tstorage = ntstorage();
      // Loop over all history data and add the appropriate contribution
      // to the derivative
      for (unsigned t = 0; t < n_tstorage; t++)
      {
        deriv += Weight(i, t) * data_pt->value(t, j);
      }
      return deriv;
    }
 
    void time_derivative(const unsigned& i,
                         Node* const& node_pt,
                         Vector<double>& deriv)
    {
      // Number of values stored in the Data object
      unsigned nvalue = node_pt->nvalue();
      deriv.assign(nvalue, 0.0);
 
      // Loop over all values
      for (unsigned j = 0; j < nvalue; j++)
      {
        deriv[j] = time_derivative(i, node_pt, j);
      }
    }
 
    double time_derivative(const unsigned& i,
                           Node* const& node_pt,
                           const unsigned& j)
    {
      double deriv = 0.0;
      unsigned n_tstorage = ntstorage();
      // Loop over all history data and add the appropriate contribution
      // to the derivative
      for (unsigned t = 0; t < n_tstorage; t++)
      {
        deriv += weight(i, t) * node_pt->value(t, j);
      }
      return deriv;
    }
 
 
    Time* const& time_pt() const
    {
#ifdef PARANOID
      if (Time_pt == 0)
      {
        std::string error_msg(
          "Time_pt is null, probably because it is a steady time stepper.");
        throw OomphLibError(
          error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Time_pt;
    }
 
    Time*& time_pt()
    {
      return Time_pt;
    }
 
    virtual double weight(const unsigned& i, const unsigned& j) const
    {
      return Weight(i, j);
    }
 
    unsigned ntstorage() const
    {
      return (Weight.ncol());
    }
 
    virtual void assign_initial_values_impulsive(Data* const& data_pt) = 0;
 
    virtual void assign_initial_positions_impulsive(Node* const& node_pt) = 0;
 
    virtual void shift_time_values(Data* const& data_pt) = 0;
 
    virtual void shift_time_positions(Node* const& node_pt) = 0;
 
    bool adaptive_flag() const
    {
      return Adaptive_Flag;
    }
 
    virtual void set_predictor_weights() {}
 
    virtual void calculate_predicted_values(Data* const& data_pt) {}
 
    virtual void calculate_predicted_positions(Node* const& node_pt) {}
 
    virtual void set_error_weights() {}
 
    virtual double temporal_error_in_position(Node* const& node_pt,
                                              const unsigned& i)
    {
      return 0.0;
    }
 
    virtual double temporal_error_in_value(Data* const& data_pt,
                                           const unsigned& i)
    {
      return 0.0;
    }
 
    virtual void actions_before_timestep(Problem* problem_pt) {}
 
    virtual void actions_after_timestep(Problem* problem_pt) {}
  };
 
 
 
 
  //====================================================================
  //====================================================================
  template<unsigned NSTEPS>
  class Steady : virtual public TimeStepper
  {
  public:
    Steady() : TimeStepper(NSTEPS + 1, 2)
    {
      Type = "Steady";
      Time_pt = &Dummy_time;
 
      // Overwrite default assignment in base constructor -- this TimeStepper
      // actually is steady all the time.
      Is_steady = true;
    }
 
 
    Steady(const Steady&) = delete;
 
    void operator=(const Steady&) = delete;
 
    unsigned order() const
    {
      return 0;
    }
 
    void assign_initial_values_impulsive(Data* const& data_pt)
    {
      // Find number of values stored
      unsigned n_value = data_pt->nvalue();
      // Loop over values
      for (unsigned j = 0; j < n_value; j++)
      {
        // Set previous values to the initial value, if not a copy
        if (data_pt->is_a_copy(j) == false)
        {
          for (unsigned t = 1; t <= NSTEPS; t++)
          {
            data_pt->set_value(t, j, data_pt->value(j));
          }
        }
      }
    }
 
    void assign_initial_positions_impulsive(Node* const& node_pt)
    {
      // Find the number of coordinates
      unsigned n_dim = node_pt->ndim();
      // Find the number of positoin types
      unsigned n_position_type = node_pt->nposition_type();
 
      // Loop over values
      for (unsigned i = 0; i < n_dim; i++)
      {
        // Set previous values to the initial value, if not a copy
        if (node_pt->position_is_a_copy(i) == false)
        {
          for (unsigned k = 0; k < n_position_type; k++)
          {
            for (unsigned t = 1; t <= NSTEPS; t++)
            {
              node_pt->x_gen(t, k, i) = node_pt->x_gen(k, i);
            }
          }
        }
      }
    }
 
 
    typedef double (*InitialConditionFctPt)(const double& t);
 
    void assign_initial_data_values(
      Data* const& data_pt, Vector<InitialConditionFctPt> initial_value_fct)
    {
      // The time history stores the previous function values
      unsigned n_time_value = ntstorage();
 
      // Find number of values stored
      unsigned n_value = data_pt->nvalue();
 
      // Loop over current and stored timesteps
      for (unsigned t = 0; t < n_time_value; t++)
      {
        // Get corresponding continous time
        double time_local = Time_pt->time(t);
 
        // Loop over values
        for (unsigned j = 0; j < n_value; j++)
        {
          data_pt->set_value(t, j, initial_value_fct[j](time_local));
        }
      }
    }
 
    void shift_time_values(Data* const& data_pt)
    {
      // Find number of values stored
      unsigned n_value = data_pt->nvalue();
 
      // Loop over the values
      for (unsigned j = 0; j < n_value; j++)
      {
        // Set previous values to the previous value, if not a copy
        if (data_pt->is_a_copy(j) == false)
        {
          // Loop over times, in reverse order
          for (unsigned t = NSTEPS; t > 0; t--)
          {
            data_pt->set_value(t, j, data_pt->value(t - 1, j));
          }
        }
      }
    }
 
    void shift_time_positions(Node* const& node_pt)
    {
      // Find the number of coordinates
      unsigned n_dim = node_pt->ndim();
      // Find the number of position types
      unsigned n_position_type = node_pt->nposition_type();
 
      // Loop over the positions
      for (unsigned i = 0; i < n_dim; i++)
      {
        // If the position is not a copy
        if (node_pt->position_is_a_copy(i) == false)
        {
          for (unsigned k = 0; k < n_position_type; k++)
          {
            // Loop over stored times, and set values to previous values
            for (unsigned t = NSTEPS; t > 0; t--)
            {
              node_pt->x_gen(t, k, i) = node_pt->x_gen(t - 1, k, i);
            }
          }
        }
      }
    }
 
    void set_weights()
    {
      // Loop over higher derivatives
      unsigned max_deriv = highest_derivative();
      for (unsigned i = 0; i < max_deriv; i++)
      {
        for (unsigned j = 0; j < NSTEPS; j++)
        {
          Weight(i, j) = 0.0;
        }
      }
 
      // Zeroth derivative must return the value itself:
      Weight(0, 0) = 1.0;
    }
 
    unsigned nprev_values() const
    {
      return NSTEPS;
    }
 
    unsigned ndt() const
    {
      return NSTEPS;
    }
 
    double weight(const unsigned& i, const unsigned& j) const
    {
      if ((i == 0) && (j == 0))
      {
        return One;
      }
      else
      {
        return Zero;
      }
    }
 
  private:
    static double One;
 
    static double Zero;
 
    // Default Time object
    static Time Dummy_time;
  };
 
 
 
 
  //====================================================================
  //====================================================================
  template<unsigned NSTEPS>
  class Newmark : public TimeStepper
  {
  public:
    Newmark() : TimeStepper(NSTEPS + 3, 2)
    {
      Type = "Newmark";
 
      // Standard Newmark parameters
      Beta1 = 0.5;
      Beta2 = 0.5;
    }
 
    Newmark(const Newmark&) = delete;
 
    void operator=(const Newmark&) = delete;
 
 
    unsigned order() const
    {
      std::string error_message =
        "Can't remember the order of the Newmark scheme";
      error_message += " -- I think it's 2nd order...\n";
 
      OomphLibWarning(
        error_message, "Newmark::order()", OOMPH_EXCEPTION_LOCATION);
      return 2;
    }
 
    void assign_initial_values_impulsive(Data* const& data_pt);
 
    void assign_initial_positions_impulsive(Node* const& node_pt);
 
    typedef double (*InitialConditionFctPt)(const double& t);
 
    void assign_initial_data_values(
      Data* const& data_pt,
      Vector<InitialConditionFctPt> initial_value_fct,
      Vector<InitialConditionFctPt> initial_veloc_fct,
      Vector<InitialConditionFctPt> initial_accel_fct);
 
 
    typedef double (*NodeInitialConditionFctPt)(const double& t,
                                                const Vector<double>& x);
 
    void assign_initial_data_values(
      Node* const& node_pt,
      Vector<NodeInitialConditionFctPt> initial_value_fct,
      Vector<NodeInitialConditionFctPt> initial_veloc_fct,
      Vector<NodeInitialConditionFctPt> initial_accel_fct);
 
    void assign_initial_data_values_stage1(const unsigned t_deriv,
                                           Data* const& data_pt);
 
    void assign_initial_data_values_stage2(Data* const& data_pt);
 
 
    void shift_time_values(Data* const& data_pt);
 
    void shift_time_positions(Node* const& node_pt);
 
    void set_weights();
 
    unsigned nprev_values() const
    {
      return NSTEPS;
    }
 
    unsigned ndt() const
    {
      return NSTEPS;
    }
 
 
  protected:
    double Beta1;
 
    double Beta2;
  };
 
 
 
 
  //====================================================================
  //====================================================================
  template<unsigned NSTEPS>
  class NewmarkBDF : public Newmark<NSTEPS>
  {
  public:
    NewmarkBDF()
    {
      this->Type = "NewmarkBDF";
      Degrade_to_bdf1_for_first_derivs = false;
      Newmark_veloc_weight.resize(NSTEPS + 3);
    }
 
    NewmarkBDF(const NewmarkBDF&) = delete;
 
    void operator=(const NewmarkBDF&) = delete;
 
    void set_weights();
 
    void shift_time_values(Data* const& data_pt);
 
    void shift_time_positions(Node* const& node_pt);
 
    void enable_degrade_first_derivatives_to_bdf1()
    {
      Degrade_to_bdf1_for_first_derivs = true;
    }
 
 
    void disable_degrade_first_derivatives_to_bdf1()
    {
      Degrade_to_bdf1_for_first_derivs = false;
    }
 
  private:
    void set_newmark_veloc_weights(const double& dt)
    {
      Newmark_veloc_weight[0] = this->Beta1 * dt * this->Weight(2, 0);
      Newmark_veloc_weight[1] = this->Beta1 * dt * this->Weight(2, 1);
      for (unsigned t = 2; t <= NSTEPS; t++)
      {
        Newmark_veloc_weight[t] = 0.0;
      }
      Newmark_veloc_weight[NSTEPS + 1] =
        1.0 + this->Beta1 * dt * this->Weight(2, NSTEPS + 1);
      Newmark_veloc_weight[NSTEPS + 2] =
        dt * (1.0 - this->Beta1) +
        this->Beta1 * dt * this->Weight(2, NSTEPS + 2);
    }
 
    bool Degrade_to_bdf1_for_first_derivs;
 
    Vector<double> Newmark_veloc_weight;
  };
 
 
 
 
  //====================================================================
  //====================================================================
  template<unsigned NSTEPS>
  class BDF : public TimeStepper
  {
    // A BDF<1> time data set consists of:
    // [y_np1, y_n, dy_n, y^P_np1]
    // Or in english:
    // * y-value at time being/just been solved for
    // * y-value at previous time
    // * approximation to y derivative at previous time (also refered to as
    //   velocity in some places, presumably it corresponds to velocity in
    //   solid mechanics).
    // * predicted y-value at time n+1
 
    // A BDF<2> time data set consists of:
    // [y_np1, y_n, y_nm1, dy_n, y^P_np1]
    // i.e. the same thing but with one more previous time value. Also the
    // derivative approximation will be more accurate.
 
    // If the adaptive flag is set to false then the final two pieces of data
    // in each are not stored (derivative and predictor value).
 
  private:
    Vector<double> Predictor_weight;
 
    double Error_weight;
 
  public:
    BDF(const bool& adaptive = false) : TimeStepper(NSTEPS + 1, 1)
    {
      Type = "BDF";
 
      // If it's adaptive, we need to allocate additional space to
      // carry along a prediction and an acceleration
      if (adaptive)
      {
        // Set the adaptive flag to be true
        Adaptive_Flag = true;
 
        // Set the size of the Predictor_Weights Vector N.B. The size is
        // correct for BDF1 and 2, but may be wrong for others.
        Predictor_weight.resize(NSTEPS + 2);
 
        // Resize the weights to the appropriate size
        Weight.resize(2, NSTEPS + 3, 0.0);
 
        // Storing predicted values in slot after other information
        Predictor_storage_index = NSTEPS + 2;
      }
 
      // Set the weight for the zero-th derivative
      Weight(0, 0) = 1.0;
    }
 
 
    BDF(const BDF&) = delete;
 
    void operator=(const BDF&) = delete;
 
    unsigned order() const
    {
      return NSTEPS;
    }
 
    void assign_initial_values_impulsive(Data* const& data_pt)
    {
      // Find number of values stored
      unsigned n_value = data_pt->nvalue();
      // Loop over values
      for (unsigned j = 0; j < n_value; j++)
      {
        // Set previous values to the initial value, if not a copy
        if (data_pt->is_a_copy(j) == false)
        {
          for (unsigned t = 1; t <= NSTEPS; t++)
          {
            data_pt->set_value(t, j, data_pt->value(j));
          }
 
          // If it's adaptive
          if (adaptive_flag())
          {
            // Initial velocity is zero
            data_pt->set_value(NSTEPS + 1, j, 0.0);
            // Initial prediction is the value
            data_pt->set_value(NSTEPS + 2, j, data_pt->value(j));
          }
        }
      }
    }
 
    void assign_initial_positions_impulsive(Node* const& node_pt)
    {
      // Find the dimension of the node
      unsigned n_dim = node_pt->ndim();
      // Find the number of position types at the node
      unsigned n_position_type = node_pt->nposition_type();
 
      // Loop over the position variables
      for (unsigned i = 0; i < n_dim; i++)
      {
        // If the position is not copied
        // We copy entire coordinates at once
        if (node_pt->position_is_a_copy(i) == false)
        {
          // Loop over the position types
          for (unsigned k = 0; k < n_position_type; k++)
          {
            // Set previous values to the initial value, if not a copy
            for (unsigned t = 1; t <= NSTEPS; t++)
            {
              node_pt->x_gen(t, k, i) = node_pt->x_gen(k, i);
            }
 
            // If it's adaptive
            if (adaptive_flag())
            {
              // Initial mesh velocity is zero
              node_pt->x_gen(NSTEPS + 1, k, i) = 0.0;
              // Initial prediction is the value
              node_pt->x_gen(NSTEPS + 2, k, i) = node_pt->x_gen(k, i);
            }
          }
        }
      }
    }
 
 
    typedef double (*InitialConditionFctPt)(const double& t);
 
    void assign_initial_data_values(
      Data* const& data_pt, Vector<InitialConditionFctPt> initial_value_fct)
    {
      // The time history stores the previous function values
      unsigned n_time_value = ntstorage();
 
      // Find number of values stored
      unsigned n_value = data_pt->nvalue();
 
      // Loop over current and stored timesteps
      for (unsigned t = 0; t < n_time_value; t++)
      {
        // Get corresponding continous time
        double time_local = Time_pt->time(t);
 
        // Loop over values
        for (unsigned j = 0; j < n_value; j++)
        {
          data_pt->set_value(t, j, initial_value_fct[j](time_local));
        }
      }
    }
 
    void shift_time_values(Data* const& data_pt)
    {
      // Find number of values stored
      unsigned n_value = data_pt->nvalue();
      // Storage for velocity need to be here to be in scope
      Vector<double> velocity(n_value);
 
      // Find the number of history values that are stored
      const unsigned nt_value = nprev_values();
 
      // If adaptive, find the velocities
      if (adaptive_flag())
      {
        time_derivative(1, data_pt, velocity);
      }
 
      // Loop over the values
      for (unsigned j = 0; j < n_value; j++)
      {
        // Set previous values to the previous value, if not a copy
        if (data_pt->is_a_copy(j) == false)
        {
          // Loop over times, in reverse order
          for (unsigned t = nt_value; t > 0; t--)
          {
            data_pt->set_value(t, j, data_pt->value(t - 1, j));
          }
 
          // If we are using the adaptive scheme
          if (adaptive_flag())
          {
            // Set the velocity
            data_pt->set_value(nt_value + 1, j, velocity[j]);
          }
        }
      }
    }
 
    void shift_time_positions(Node* const& node_pt)
    {
      // Find the number of coordinates
      unsigned n_dim = node_pt->ndim();
      // Find the number of position types
      unsigned n_position_type = node_pt->nposition_type();
 
      // Find number of stored timesteps
      unsigned n_tstorage = ntstorage();
 
      // Storage for the velocity
      double velocity[n_position_type][n_dim];
 
      // If adaptive, find the velocities
      if (adaptive_flag())
      {
        // Loop over the variables
        for (unsigned i = 0; i < n_dim; i++)
        {
          for (unsigned k = 0; k < n_position_type; k++)
          {
            // Initialise velocity to zero
            velocity[k][i] = 0.0;
            // Loop over all history data
            for (unsigned t = 0; t < n_tstorage; t++)
            {
              velocity[k][i] += Weight(1, t) * node_pt->x_gen(t, k, i);
            }
          }
        }
      }
 
      // Loop over the positions
      for (unsigned i = 0; i < n_dim; i++)
      {
        // If the position is not a copy
        if (node_pt->position_is_a_copy(i) == false)
        {
          // Loop over the position types
          for (unsigned k = 0; k < n_position_type; k++)
          {
            // Loop over stored times, and set values to previous values
            for (unsigned t = NSTEPS; t > 0; t--)
            {
              node_pt->x_gen(t, k, i) = node_pt->x_gen(t - 1, k, i);
            }
 
            // If we are using the adaptive scheme, set the velocity
            if (adaptive_flag())
            {
              node_pt->x_gen(NSTEPS + 1, k, i) = velocity[k][i];
            }
          }
        }
      }
    }
 
    void set_weights();
 
    unsigned nprev_values() const
    {
      return NSTEPS;
    }
 
    unsigned ndt() const
    {
      return NSTEPS;
    }
 
    void set_predictor_weights();
 
    void calculate_predicted_positions(Node* const& node_pt);
 
    void calculate_predicted_values(Data* const& data_pt);
 
    void set_error_weights();
 
    double temporal_error_in_position(Node* const& node_pt, const unsigned& i);
 
    double temporal_error_in_value(Data* const& data_pt, const unsigned& i);
  };
 
} // namespace oomph
 
#endif