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-2025 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
#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 to keep track of discrete/continous time. It is essential
  /// to have a single Time object when using multiple time-stepping schemes;
  /// e.g., in fluid-structure interaction problems, it is common to use
  /// different schemes for the fluid and solid domains.
  /// Storage is allocated for the current value
  /// of the (continuous) time and a limited history of previous timesteps.
  /// The number of previous timesteps must be equal to the number required
  /// by the "highest order" scheme.
  //====================================================================
  class Time
  {
  private:
    /// Pointer to the value of the continuous time.
    double Continuous_time;
 
    /// Vector that stores the values of the current and previous timesteps.
    Vector<double> Dt;
 
  public:
    /// Constructor: Do not allocate any storage for previous timesteps,
    /// but set the initial value of the time to zero
    Time() : Continuous_time(0.0) {}
 
    /// Constructor: Pass the number of timesteps to be stored
    /// and set the initial value of time to zero.
    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);
    }
 
    /// Broken copy constructor
    Time(const Time&) = delete;
 
    /// Broken assignment operator
    void operator=(const Time&) = delete;
 
    /// Resize the vector holding the number of previous timesteps
    /// and initialise the new values to zero.
    void resize(const unsigned& n_dt)
    {
      Dt.resize(n_dt, 0.0);
    }
 
    /// Set all timesteps to the same value, dt.
    void initialise_dt(const double& dt_)
    {
      unsigned ndt = Dt.size();
      Dt.assign(ndt, dt_);
    }
 
    /// Set the value of the timesteps to be equal to the values passed
    /// in a vector
    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];
      }
    }
 
    /// Destructor: empty
    ~Time() {}
 
    /// Return the current value of the continuous time
    double& time()
    {
      return Continuous_time;
    }
 
    /// Return the number of timesteps stored
    unsigned ndt() const
    {
      return Dt.size();
    }
 
    /// Return the value of the t-th stored timestep (t=0: present; t>0:
    /// previous).
    double& dt(const unsigned& t = 0)
    {
      return Dt[t];
    }
 
    /// Return the value of the t-th stored timestep (t=0: present; t>0:
    /// previous), const version.
    double dt(const unsigned& t = 0) const
    {
      return Dt[t];
    }
 
    /// Return the value of the continuous time at the t-th previous
    /// time level (t=0: current; t>0 previous).
    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;
    }
 
    /// Update all stored values of dt by shifting each value along
    /// the array. This function must be called before starting to solve at a
    /// new time level.
    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];
      }
    }
  };
 
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Base class for time-stepping schemes.
  /// Timestepper provides an approximation of the temporal derivatives of
  /// Data such that the i-th derivative of the j-th value in Data is
  /// represented as
  /// \code
  /// unsigned n_value=data_pt->nvalue();
  /// Vector<double> time_derivative(n_value);
  /// for (unsigned j=0;j<n_value;j++)
  ///  {
  ///  time_derivative[j]=0.0;
  ///  for (unsigned t=0;t<ntstorage();t++)
  ///   {
  ///    time_derivative[j] += data_pt->value(t,j)*weight(i,t);
  ///   }
  ///  }
  /// \endcode
  /// where the time index \c t is such that
  /// - \c t \c = \c 0 refers to the current value of the Data
  /// - \c t \c > \c 0 refers to the `history' values, i.e. the doubles that
  ///   are stored in Data to represent the value's time derivatives. For BDF
  ///   schemes these doubles are the values at previous times;
  ///   for other timestepping schemes they can represent quantities
  ///   such as previous first and second derivatives, etc.
  ///
  /// TimeSteppers can evaluate all derivatives up to their \c order().
  ///
  /// The first \c nprev_values() `history values' represent the
  /// values at the previous timesteps. For BDF schemes,
  /// \c nprev_values() \c = \c  ntstorage()-1 , i.e. all `history values'
  /// represent previous values. For \c Newmark<NSTEPS>, only the first
  /// NSTEPS `history values' represent previous values (NSTEPS=1
  /// gives the normal Newmark scheme).
  //====================================================================
  class TimeStepper
  {
  protected:
    /// Pointer to discrete time storage scheme
    Time* Time_pt;
 
    /// Storage for the weights associated with the timestepper
    DenseMatrix<double> Weight;
 
    /// String that indicates the type of the timestepper
    /// (e.g. "BDF", "Newmark", etc.)
    std::string Type;
 
    /// Boolean variable to indicate whether the timestepping scheme can
    /// be adaptive
    bool Adaptive_Flag;
 
    /// Bool to indicate if the timestepper is steady, i.e. its
    /// time-derivatives evaluate to zero. This status may be achieved
    /// temporarily by calling make_steady(). It can be reset to the
    /// appropriate default by the function undo_make_steady().
    bool Is_steady;
 
    /// Boolean to indicate if the timestepper will output warnings when
    /// setting possibly an incorrect number of initial data values from
    /// function pointers
    bool Shut_up_in_assign_initial_data_values;
 
    /// Flag: is adaptivity done by taking a separate step using an
    /// ExplicitTimeStepper object?
    bool Predict_by_explicit_step;
 
    /// Pointer to explicit time stepper to use as predictor if
    /// Predict_by_explicit_step is set.
    /// ??ds merge the two? predict by explicit if pointer is set?
    ExplicitTimeStepper* Explicit_predictor_pt;
 
    /// Store the time that the predicted values currently stored are at, to
    /// compare for paranoid checks.
    double Predicted_time;
 
    /// The time-index in each Data object where predicted values are
    /// stored. -1 if not set.
    int Predictor_storage_index;
 
  public:
    /// Constructor. Pass the amount of storage required by
    /// timestepper (present value + history values) and the
    /// order of highest time-derivative.
    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;
    }
 
    /// Broken empty constructor
    TimeStepper()
    {
      throw OomphLibError("Don't call default constructor for TimeStepper!",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Broken copy constructor
    TimeStepper(const TimeStepper&) = delete;
 
    /// Broken assignment operator
    void operator=(const TimeStepper&) = delete;
 
    /// virtual destructor
    virtual ~TimeStepper();
 
    /// Highest order derivative that the scheme can compute
    unsigned highest_derivative() const
    {
      return Weight.nrow() - 1;
    }
 
    /// Actual order (accuracy) of the scheme
    virtual unsigned order() const
    {
      return 0;
    }
 
    /// Return current value of continous time
    // (can't have a paranoid test for null pointers because this could be
    // used as a set function)
    double& time()
    {
      return Time_pt->time();
    }
 
    /// Return current value of continous 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();
    }
 
    /// Number of timestep increments that are required by the scheme
    virtual unsigned ndt() const = 0;
 
    /// Number of previous values needed to calculate the value at the
    /// current time. i.e. how many previous values must we loop over to
    /// calculate the values at the evaluation time. For most methods this is
    /// 1, i.e. just use the value at t_{n+1}. See midpoint method for a
    /// counter-example.
    virtual unsigned nprev_values_for_value_at_evaluation_time() const
    {
      return 1;
    }
 
    /// Number of previous values available: 0 for static, 1 for BDF<1>,...
    virtual unsigned nprev_values() const = 0;
 
    /// Function to set the weights for present timestep (don't
    /// need to pass present timestep or previous timesteps as they
    /// are available via Time_pt)
    virtual void set_weights() = 0;
 
    /// Function to make the time stepper temporarily steady. This
    /// is trivially achieved by setting all the weights to zero
    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;
    }
 
    /// Flag to indicate if a timestepper has been made steady (possibly
    /// temporarily to switch off time-dependence)
    bool is_steady() const
    {
      return Is_steady;
    }
 
    /// Flag: is adaptivity done by taking a separate step using an
    /// ExplicitTimeStepper object?
    bool predict_by_explicit_step() const
    {
      return Predict_by_explicit_step;
    }
 
    /// Get the pointer to the explicit timestepper to use as a predictor in
    /// adaptivity if Predict_by_explicit_step is set.
    ExplicitTimeStepper* explicit_predictor_pt()
    {
      return Explicit_predictor_pt;
    }
 
    /// Set the pointer to the explicit timestepper to use as a predictor in
    /// adaptivity if Predict_by_explicit_step is set.
    void set_predictor_pt(ExplicitTimeStepper* _pred_pt)
    {
      Explicit_predictor_pt = _pred_pt;
    }
 
    /// Set the time that the current predictions apply for, only needed for
    /// paranoid checks when doing Predict_by_explicit_step.
    void update_predicted_time(const double& new_time)
    {
      Predicted_time = new_time;
    }
 
    /// Check that the predicted values are the ones we want.
    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
    }
 
    /// Return the time-index in each Data where predicted values are
    /// stored if the timestepper is adaptive.
    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
    }
 
    /// Enable the output of warnings due to possible fct pointer vector
    /// size mismatch in assign_initial_data_values (Default)
    void enable_warning_in_assign_initial_data_values()
    {
      Shut_up_in_assign_initial_data_values = false;
    }
 
    /// Disable the output of warnings due to possible fct pointer vector
    /// size mismatch in assign_initial_data_values
    void disable_warning_in_assign_initial_data_values()
    {
      Shut_up_in_assign_initial_data_values = true;
    }
 
    /// Get a (const) pointer to the weights.
    const DenseMatrix<double>* weights_pt() const
    {
      return &Weight;
    }
 
    /// Reset the is_steady status of a specific TimeStepper to its
    /// default and re-assign the weights.
    virtual void undo_make_steady()
    {
      Is_steady = false;
      set_weights();
    }
 
    /// Return string that indicates the type of the timestepper
    /// (e.g. "BDF", "Newmark", etc.)
    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.
 
    /// Evaluate i-th derivative of all values in Data and return in
    /// Vector deriv[].
    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);
      }
    }
 
    /// Evaluate i-th derivative of j-th value in Data.
    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;
    }
 
    /// Evaluate i-th derivative of all values in Node and return in
    /// Vector deriv[] (this can't be simply combined with time_derivative(..,
    /// Data, ...) because of differences with haning nodes).
    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);
      }
    }
 
    /// Evaluate i-th derivative of j-th value in Node. Note the use of
    /// the node's value() function so that hanging nodes are taken into
    /// account (this is why the functions for Data and Node cannot be
    /// combined through simple polymorphism: value is not virtual).
    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;
    }
 
 
    /// Access function for the pointer to time (const version)
    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;
    }
 
    /// Access function for the pointer to time (can't have a paranoid test
    /// for null pointers because this is used as a set function).
    Time*& time_pt()
    {
      return Time_pt;
    }
 
    /// Access function for j-th weight for the i-th derivative
    virtual double weight(const unsigned& i, const unsigned& j) const
    {
      return Weight(i, j);
    }
 
    /// Return the number of doubles required to represent history
    /// (one for steady)
    unsigned ntstorage() const
    {
      return (Weight.ncol());
    }
 
    /// Initialise the time-history for the Data values
    /// corresponding to an impulsive start.
    virtual void assign_initial_values_impulsive(Data* const& data_pt) = 0;
 
    /// Initialiset the positions for the node corresponding to
    /// an impulsive start
    virtual void assign_initial_positions_impulsive(Node* const& node_pt) = 0;
 
    /// This function advances the Data's time history so that
    /// we can move on to the next timestep
    virtual void shift_time_values(Data* const& data_pt) = 0;
 
    /// This function advances the time history of the positions
    /// at a node. The default should be OK, but would need to be overloaded
    virtual void shift_time_positions(Node* const& node_pt) = 0;
 
    /// Function to indicate whether the scheme is adaptive (false by default)
    bool adaptive_flag() const
    {
      return Adaptive_Flag;
    }
 
    /// Set the weights for the predictor
    /// previous timestep (currently empty -- overwrite for specific scheme)
    virtual void set_predictor_weights() {}
 
    /// Do the predictor step for data stored in a Data object
    /// (currently empty -- overwrite for specific scheme)
    virtual void calculate_predicted_values(Data* const& data_pt) {}
 
    /// Do the predictor step for the positions at a node
    /// (currently empty --- overwrite for a specific scheme)
    virtual void calculate_predicted_positions(Node* const& node_pt) {}
 
    /// Set the weights for the error computation,
    /// (currently empty -- overwrite for specific scheme)
    virtual void set_error_weights() {}
 
    /// Compute the error in the position i at a node
    /// zero here -- overwrite for specific scheme.
    virtual double temporal_error_in_position(Node* const& node_pt,
                                              const unsigned& i)
    {
      return 0.0;
    }
 
    /// Compute the error in the value i in a Data structure
    /// zero here -- overwrite for specific scheme.
    virtual double temporal_error_in_value(Data* const& data_pt,
                                           const unsigned& i)
    {
      return 0.0;
    }
 
    /// Interface for any actions that need to be performed before a time
    /// step.
    virtual void actions_before_timestep(Problem* problem_pt) {}
 
    /// Interface for any actions that need to be performed after a time
    /// step.
    virtual void actions_after_timestep(Problem* problem_pt) {}
  };
 
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Faux time-stepper for steady problems. Allows storage
  /// for NSTEPS previous values.
  //====================================================================
  template<unsigned NSTEPS>
  class Steady : virtual public TimeStepper
  {
  public:
    /// Constructor: Creates storage for NSTEPS previous timesteps
    /// and can evaluate up to 2nd derivatives (though it doesn't
    /// actually do anything -- all time-derivatives
    /// evaluate to zero)
    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;
    }
 
 
    /// Broken copy constructor
    Steady(const Steady&) = delete;
 
    /// Broken assignment operator
    void operator=(const Steady&) = delete;
 
    /// Return the actual order of the scheme. Returning zero here --
    /// doesn't make much sense, though
    unsigned order() const
    {
      return 0;
    }
 
    ///  Initialise the time-history for the Data values,
    /// corresponding to an impulsive start.
    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));
          }
        }
      }
    }
 
    ///  Initialise the time-history for the nodal positions
    /// corresponding to an impulsive start.
    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 for function that returns the (scalar) initial
    /// value at a given value of the continuous time t.
    typedef double (*InitialConditionFctPt)(const double& t);
 
    ///  Initialise the time-history for the Data values,
    /// corresponding to given time history, specified by
    /// Vector of function pointers.
    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));
        }
      }
    }
 
    /// This function updates the Data's time history so that
    /// we can advance to the next timestep. As for BDF schemes,
    /// we simply push the values backwards...
    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));
          }
        }
      }
    }
 
    /// This function advances the time history of the positions
    /// at a node.
    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);
            }
          }
        }
      }
    }
 
    /// Set weights
    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;
    }
 
    /// Number of previous values available.
    unsigned nprev_values() const
    {
      return NSTEPS;
    }
 
    /// Number of timestep increments that need to be stored by the scheme
    unsigned ndt() const
    {
      return NSTEPS;
    }
 
    /// Dummy: Access function for j-th weight for the i-th derivative
    double weight(const unsigned& i, const unsigned& j) const
    {
      if ((i == 0) && (j == 0))
      {
        return One;
      }
      else
      {
        return Zero;
      }
    }
 
  private:
    /// Static variable to hold the value 1.0
    static double One;
 
    /// Static variable to hold the value 0.0
    static double Zero;
 
    // Default Time object
    static Time Dummy_time;
  };
 
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Newmark scheme for second time deriv. Stored data represents
  /// - t=0: value at at present time, Time_pt->time()
  /// - t=1: value at previous time, Time_pt->time()-dt
  /// -  ...
  /// - t=NSTEPS:   value at previous time, Time_pt->time()-NSTEPS*dt
  /// - t=NSTEPS+1: 1st time deriv (= "velocity") at previous time,
  ///               Time_pt->time()-dt
  /// - t=NSTEPS+2: 2nd time deriv (= "acceleration") at previous time,
  ///                Time_pt->time()-dt
  ///
  /// NSTEPS=1 gives normal Newmark.
  //====================================================================
  template<unsigned NSTEPS>
  class Newmark : public TimeStepper
  {
  public:
    /// Constructor: Pass pointer to global time. We set up a
    /// timestepping scheme with NSTEPS+2 doubles to represent the history and
    /// the highest deriv is 2.
    Newmark() : TimeStepper(NSTEPS + 3, 2)
    {
      Type = "Newmark";
 
      // Standard Newmark parameters
      Beta1 = 0.5;
      Beta2 = 0.5;
    }
 
    /// Broken copy constructor
    Newmark(const Newmark&) = delete;
 
    /// Broken assignment operator
    void operator=(const Newmark&) = delete;
 
 
    /// The actual order (accuracy of the scheme)
    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;
    }
 
    /// Initialise the time-history for the values,
    /// corresponding to an impulsive start.
    void assign_initial_values_impulsive(Data* const& data_pt);
 
    /// Initialise the time-history for the values,
    /// corresponding to an impulsive start.
    void assign_initial_positions_impulsive(Node* const& node_pt);
 
    /// Typedef for function that returns the (scalar) initial
    /// value at a given value of the continuous time t.
    typedef double (*InitialConditionFctPt)(const double& t);
 
    ///  Initialise the time-history for the Data values,
    /// so that the Newmark representations for current veloc and
    /// acceleration are exact.
    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 for function that returns the (scalar) initial
    /// value at a given value of the continuous time t and the spatial
    /// coordinate -- appropriate for assignement of initial conditions for
    /// nodes
    typedef double (*NodeInitialConditionFctPt)(const double& t,
                                                const Vector<double>& x);
 
    ///  Initialise the time-history for the nodal values,
    /// so that the Newmark representations for current veloc and
    /// acceleration are exact.
    void assign_initial_data_values(
      Node* const& node_pt,
      Vector<NodeInitialConditionFctPt> initial_value_fct,
      Vector<NodeInitialConditionFctPt> initial_veloc_fct,
      Vector<NodeInitialConditionFctPt> initial_accel_fct);
 
    ///  First step in a two-stage procedure to assign
    /// the history values for the Newmark scheme so
    /// that the veloc and accel that are computed by the scheme
    /// are correct at the current time.
    ///
    /// Call this function for t_deriv=0,1,2,3.
    /// When calling with
    /// - t_deriv=0 :  data_pt(0,*) should contain the values at the
    ///                previous timestep
    /// - t_deriv=1 :  data_pt(0,*) should contain the current values;
    ///                they get stored (temporarily) in data_pt(1,*).
    /// - t_deriv=2 :  data_pt(0,*) should contain the current velocities
    ///                (first time derivs); they get stored (temporarily)
    ///                in data_pt(NSTEPS+1,*).
    /// - t_deriv=3 :  data_pt(0,*) should contain the current accelerations
    ///                (second time derivs); they get stored (temporarily)
    ///                in data_pt(NSTEPS+2,*).
    /// .
    /// Follow this by calls to
    /// \code
    /// assign_initial_data_values_stage2(...)
    /// \endcode
    void assign_initial_data_values_stage1(const unsigned t_deriv,
                                           Data* const& data_pt);
 
    /// Second step in a two-stage procedure to assign
    /// the history values for the Newmark scheme so
    /// that the veloc and accel that are computed by the scheme
    /// are correct at the current time.
    ///
    /// This assigns appropriate values for the "previous
    /// velocities and accelerations" so that their current
    /// values, which were defined in assign_initial_data_values_stage1(...),
    /// are represented exactly by the Newmark scheme.
    void assign_initial_data_values_stage2(Data* const& data_pt);
 
 
    /// This function updates the Data's time history so that
    /// we can advance to the next timestep.
    void shift_time_values(Data* const& data_pt);
 
    /// This function updates a nodal time history so that
    /// we can advance to the next timestep.
    void shift_time_positions(Node* const& node_pt);
 
    /// Set weights
    void set_weights();
 
    /// Number of previous values available.
    unsigned nprev_values() const
    {
      return NSTEPS;
    }
 
    /// Number of timestep increments that need to be stored by the scheme
    unsigned ndt() const
    {
      return NSTEPS;
    }
 
 
  protected:
    /// First Newmark parameter (usually 0.5)
    double Beta1;
 
    /// Second Newmark parameter (usually 0.5)
    double Beta2;
  };
 
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Newmark scheme for second time deriv with first derivatives
  /// calculated using BDF. . Stored data represents
  /// - t=0: value at at present time, Time_pt->time()
  /// - t=1: value at previous time, Time_pt->time()-dt
  /// -  ...
  /// - t=NSTEPS:   value at previous time, Time_pt->time()-NSTEPS*dt
  /// - t=NSTEPS+1: 1st time deriv (= "velocity") at previous time,
  ///               Time_pt->time()-dt
  /// - t=NSTEPS+2: 2nd time deriv (= "acceleration") at previous time,
  ///                Time_pt->time()-dt
  ///
  /// NSTEPS=1 gives normal Newmark.
  //====================================================================
  template<unsigned NSTEPS>
  class NewmarkBDF : public Newmark<NSTEPS>
  {
  public:
    /// Constructor: Pass pointer to global time. We set up a
    /// timestepping scheme with NSTEPS+2 doubles to represent the history and
    /// the highest deriv is 2.
    NewmarkBDF()
    {
      this->Type = "NewmarkBDF";
      Degrade_to_bdf1_for_first_derivs = false;
      Newmark_veloc_weight.resize(NSTEPS + 3);
    }
 
    /// Broken copy constructor
    NewmarkBDF(const NewmarkBDF&) = delete;
 
    /// Broken assignment operator
    void operator=(const NewmarkBDF&) = delete;
 
    /// Set weights
    void set_weights();
 
    /// This function updates the Data's time history so that
    /// we can advance to the next timestep.
    void shift_time_values(Data* const& data_pt);
 
    /// This function updates a nodal time history so that
    /// we can advance to the next timestep.
    void shift_time_positions(Node* const& node_pt);
 
    /// Degrade scheme to first order BDF (for first derivs/veloc);
    /// usually for start-up.
    void enable_degrade_first_derivatives_to_bdf1()
    {
      Degrade_to_bdf1_for_first_derivs = true;
    }
 
 
    /// Disable degradation to first order BDF.
    void disable_degrade_first_derivatives_to_bdf1()
    {
      Degrade_to_bdf1_for_first_derivs = false;
    }
 
  private:
    /// Set original Newmark weights for velocities (needed when
    /// shifting history values -- they're used when updating the
    /// previous accelerations and doing this with bdf can make the
    /// scheme unstable...
    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);
    }
 
    /// Boolean flag to indicate degradation of scheme to first
    /// order BDF (for first derivs/veloc); usually for start-up.
    bool Degrade_to_bdf1_for_first_derivs;
 
    /// Original Newmark weights for velocities (needed when
    /// shifting history values -- they're used when updating the
    /// previous accelerations and doing this with bdf can make the
    /// scheme unstable...
    Vector<double> Newmark_veloc_weight;
  };
 
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //====================================================================
  /// Templated class for BDF-type time-steppers with fixed or
  /// variable timestep.
  /// 1st time derivative recovered directly from the previous function
  /// values. Template parameter represents the number of previous timesteps
  /// stored, so that BDF<1> is the classical first order
  /// backward Euler scheme.
  /// Need to reset weights after every change in timestep.
  //====================================================================
  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:
    /// Private data for the predictor weights
    Vector<double> Predictor_weight;
 
    /// Private data for the error weight
    double Error_weight;
 
  public:
    /// Constructor for the case when we allow adaptive timestepping
    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;
    }
 
 
    /// Broken copy constructor
    BDF(const BDF&) = delete;
 
    /// Broken assignment operator
    void operator=(const BDF&) = delete;
 
    /// Return the actual order of the scheme
    unsigned order() const
    {
      return NSTEPS;
    }
 
    ///  Initialise the time-history for the Data values,
    /// corresponding to an impulsive start.
    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));
          }
        }
      }
    }
 
    ///  Initialise the time-history for the nodal positions
    /// corresponding to an impulsive start.
    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 for function that returns the (scalar) initial
    /// value at a given value of the continuous time t.
    typedef double (*InitialConditionFctPt)(const double& t);
 
    ///  Initialise the time-history for the Data values,
    /// corresponding to given time history, specified by
    /// Vector of function pointers.
    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));
        }
      }
    }
 
    /// This function updates the Data's time history so that
    /// we can advance to the next timestep. For BDF schemes,
    /// we simply push the values backwards...
    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]);
          }
        }
      }
    }
 
    /// This function advances the time history of the positions
    /// at a node.
    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];
            }
          }
        }
      }
    }
 
    /// Set the weights
    void set_weights();
 
    /// Number of previous values available.
    unsigned nprev_values() const
    {
      return NSTEPS;
    }
 
    /// Number of timestep increments that need to be stored by the scheme
    unsigned ndt() const
    {
      return NSTEPS;
    }
 
    /// Function to set the predictor weights
    void set_predictor_weights();
 
    /// Function to calculate predicted positions at a node
    void calculate_predicted_positions(Node* const& node_pt);
 
    /// Function to calculate predicted data values in a Data object
    void calculate_predicted_values(Data* const& data_pt);
 
    /// Function to set the error weights
    void set_error_weights();
 
    /// Compute the error in the position i at a node
    double temporal_error_in_position(Node* const& node_pt, const unsigned& i);
 
    /// Compute the error in the value i in a Data structure
    double temporal_error_in_value(Data* const& data_pt, const unsigned& i);
  };
 
} // namespace oomph
 
#endif