generalised_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//====================================================================
#ifndef OOMPH_GENERALISED_TIMESTEPPERS_HEADER
#define OOMPH_GENERALISED_TIMESTEPPERS_HEADER
 
#include <typeinfo>
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// OOMPH-LIB headers
#include "timesteppers.h"
 
namespace oomph
{
  //======================================================================
  /// Generalised timestepper that can serve a variety of purposes
  /// in continuation, bifurcation detection and periodic-orbit computations.
  /// The key generalisation is that more than one of the entries is actually
  /// a degree of freedom in the problem. These are distinct from our
  /// standard (implict) Timesteppers in which the only dof is the current
  /// value (first entry in the storage scheme). These objects will typically
  /// be used to replace exisiting timesteppers for specific tasks.
  //====================================================================
  class GeneralisedTimeStepper : public TimeStepper
  {
  protected:
    // Set the number of entries that correspond to dof storage
    unsigned Ndof_storage_entries;
 
    // Vector that represents the storage entries
    // Vector<unsigned> Dof_storage_index;
 
    /// Constructor that can only be called by derived objects.
    /// Pass the information directly through to the TimeStepper object
    GeneralisedTimeStepper(const unsigned& n_tstorage,
                           const unsigned& max_deriv,
                           const unsigned& ndof_storage_entries = 1)
      : TimeStepper(n_tstorage, max_deriv),
        Ndof_storage_entries(ndof_storage_entries)
    {
      // Dof_storage_index.resize(1,0.0);
    }
 
    /// Broken empty constructor
    GeneralisedTimeStepper() : TimeStepper() {}
 
    /// Broken copy constructor
    GeneralisedTimeStepper(const GeneralisedTimeStepper&) = delete;
 
    /// Broken assignment operator
    void operator=(const GeneralisedTimeStepper&) = delete;
 
  public:
    /// Return the number of entries that correspond to dof storage
    unsigned ndof_storage_entries() const
    {
      return Ndof_storage_entries;
    }
 
    /// Return the i-th storage index
    // unsigned dof_storage_index(const unsigned &i) {return
    // Dof_storage_index[i];}
  };
 
  //========================================================================
  /// GeneralisedTimestepper used to store the arclength derivatives
  /// and pervious solutions required in continuation problems. The data
  /// is stored as auxilliary data in the (fake) TimeStepper so that
  /// spatial adaptivity will be handled automatically through our standard
  /// mechanisms. The adopted storage scheme is that the continuation
  /// derivatives will be stored at the first auxilliary value and
  /// the previous value will be the second auixilliary value
  //====================================================================
  class ContinuationStorageScheme : public GeneralisedTimeStepper
  {
    // Store the offset for the derivatives
    unsigned Dof_derivative_offset;
 
    // Store the offset for the current values of the dofs
    unsigned Dof_current_offset;
 
  public:
    /// Constructor for the case when we allow adaptive continuation
    /// It can evaulate up to second derivatives, but doesn't do anything
    /// the time-derivatives evaluate to zero.
    ContinuationStorageScheme()
      : GeneralisedTimeStepper(3, 2),
        Dof_derivative_offset(1),
        Dof_current_offset(2)
    {
      Type = "ContinuationStorageScheme";
      Is_steady = true;
    }
 
 
    /// Broken copy constructor
    ContinuationStorageScheme(const ContinuationStorageScheme&) = delete;
 
    /// Modify the scheme based on the underlying timestepper
    void modify_storage(GeneralisedTimeStepper* const& time_stepper_pt)
    {
      // Get the number of "dofs" in the existing timestepper
      this->Ndof_storage_entries = time_stepper_pt->ndof_storage_entries();
      // Get the current amount of storage
      unsigned n_tstorage = time_stepper_pt->ntstorage();
 
      // Find the offsets which is always relative to the number of dofs stored
      // in the existing timestepper
      Dof_derivative_offset = n_tstorage;
      Dof_current_offset = n_tstorage + this->Ndof_storage_entries;
 
      // Set the new amount of storage twice the dofs to store parameter
      // derivatives and initial values plus the original storage
      unsigned n_new_tstorage = 2 * this->Ndof_storage_entries + n_tstorage;
 
      // Resize the weights accordingly
      Weight.resize(3, n_new_tstorage, 0.0);
      // Set the weight for the zero-th derivative which is always 1.0
      Weight(0, 0) = 1.0;
    }
 
 
    /// Broken assignment operator
    void operator=(const ContinuationStorageScheme&) = delete;
 
    /// Return the actual order of the scheme. It's a steady
    /// scheme so it's zero, but that doesn't really make sense.
    unsigned order() const
    {
      return 0;
    }
 
    /// This is a steady scheme, so you can't do this
    void undo_make_steady()
    {
      Is_steady = true;
    }
 
    /// Broken initialisation the time-history for the Data values
    /// corresponding to an impulsive start.
    void assign_initial_values_impulsive(Data* const& data_pt)
    {
      OomphLibWarning(
        "Cannot perform impulsive start for ContinuationStorageScheme",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Broken initialisation of
    /// the positions for the node corresponding to an impulsive start
    void assign_initial_positions_impulsive(Node* const& node_pt)
    {
      OomphLibWarning(
        "Cannot perform impulsive start for ContinuationStorageScheme",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Broken shifting of time values
    void shift_time_values(Data* const& data_pt)
    {
      throw OomphLibError(
        "Cannot shift time values forContinuationStorageScheme",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Broken shifting of time positions
    void shift_time_positions(Node* const& node_pt)
    {
      throw OomphLibError(
        "Cannot shift time positions forContinuationStorageScheme",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Set the weights (Do nothing)
    void set_weights() {}
 
    /// Number of previous values available.
    unsigned nprev_values() const
    {
      return 0;
    }
 
    /// Number of timestep increments that need to be stored by the scheme
    unsigned ndt() const
    {
      return 0;
    }
 
    /// Set consistent values of the derivatives and current value when
    /// the data is pinned. This must be done by the "timestepper" because only
    /// it knows the local storage scheme
    void set_consistent_pinned_values(Data* const& data_pt)
    {
#ifdef PARANOID
      // If the data is not associated with the continuation time stepper then
      // complain
      if (this != data_pt->time_stepper_pt())
      {
        std::ostringstream error_stream;
        error_stream
          << "Data object " << data_pt << " has timestepper of type "
          << typeid(data_pt->time_stepper_pt()).name() << "\n"
          << "Please set the data's timestepper to be a "
          << "ContinuationStorageScheme before calling this function\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Loop over the values in the data object
      const unsigned n_value = data_pt->nvalue();
      for (unsigned i = 0; i < n_value; ++i)
      {
        // Only bother to do anything if the value is pinned and not a copy
        if (data_pt->is_pinned(i) && (data_pt->is_a_copy(i) == false))
        {
          // ASSUMPTION storage is always at the "front"
          for (unsigned t = 0; t < this->Ndof_storage_entries; t++)
          {
            // Set the stored derivatve to be zero
            data_pt->set_value(Dof_derivative_offset + t, i, 0.0);
            // Set the stored current value to be the same as the present value
            data_pt->set_value(Dof_current_offset + t, i, data_pt->value(t, i));
          }
        }
      }
    }
 
    /// Set consistent values of the derivatives and current value when
    /// the Nodes position is pinned. This must be done by the "timestepper"
    /// because only it knows the local storage scheme
    void set_consistent_pinned_positions(Node* const& node_pt)
    {
      // Only need to do anything if this is a solid node
      if (SolidNode* const solid_node_pt = dynamic_cast<SolidNode*>(node_pt))
      {
#ifdef PARANOID
        // If the data is not associated with the continuation time stepper then
        // complain
        if (this != node_pt->position_time_stepper_pt())
        {
          std::ostringstream error_stream;
          error_stream
            << "Node object " << node_pt << " has position timestepper of type "
            << typeid(node_pt->position_time_stepper_pt()).name() << "\n"
            << "Please set the Node's position timestepper to be a "
            << "ContinuationStorageScheme before calling this function\n";
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
 
        // Find the number of coordinates
        const unsigned n_dim = solid_node_pt->ndim();
        // Find the number of position types
        const unsigned n_position_type = solid_node_pt->nposition_type();
 
        // Loop over physical coordinates
        for (unsigned i = 0; i < n_dim; i++)
        {
          // Set the appropriate values if it's not a copy
          if (solid_node_pt->position_is_a_copy(i) == false)
          {
            // Loop over generalised dofs
            for (unsigned k = 0; k < n_position_type; k++)
            {
              // If it's pinned then set the "history" values
              if (solid_node_pt->position_is_pinned(k, i))
              {
                for (unsigned t = 0; t < Ndof_storage_entries; t++)
                {
                  // Set the derivative to 0
                  solid_node_pt->x_gen(Dof_derivative_offset + t, k, i) = 0.0;
                  // Set the stored current value to the present value
                  solid_node_pt->x_gen(Dof_current_offset + t, k, i) =
                    solid_node_pt->x_gen(t, k, i);
                }
              }
            }
          }
        }
      }
    }
 
    // Return the stored derivative offset
    unsigned dof_derivative_offset()
    {
      return Dof_derivative_offset;
    }
 
    // Return the offset for the current values
    unsigned dof_current_offset()
    {
      return Dof_current_offset;
    }
  };
 
} // namespace oomph
 
#endif