explicit_timesteppers.cc

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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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//====================================================================
//Non-inline member functions for the explicit timesteppers
#include "explicit_timesteppers.h"
#include "timesteppers.h"

namespace oomph
{
 //===================================================================
 ///Dummy value of time always set to zero
 //==================================================================
 double ExplicitTimeSteppableObject::Dummy_time_value = 0.0;

 //==================================================================
 ///\short A single virtual function that returns the residuals
 ///vector multiplied by the inverse mass matrix
 //=================================================================
 void ExplicitTimeSteppableObject::
 get_dvaluesdt(DoubleVector &minv_res)
 {
  std::ostringstream error_stream;
  error_stream 
   << "Empty default function called.\n"
   << "The function must return the solution x of the linear system\n"
   << "                    M x = R\n"
   << "in order for the object to be used by an ExplicitTimeStepper.\n"
   << "NOTE: It is the responsibility of the object to set the size \n"
   << "      of the vector x\n";
  
  
  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
 }
 
 //=======================================================================
 /// Function that should get the values of the dofs in the object
 //=======================================================================
 void ExplicitTimeSteppableObject::get_dofs(DoubleVector &dofs) const
 {
  std::ostringstream error_stream;
  error_stream 
   << "Empty default function called.\n"
   << "The function must return the current values of the degrees of \n"
   << "freedom in the object.\n"
   << "Note: It is the responsibility of the object to set the size\n" 
   << "of the vector\n";
  
  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
  
 }

 //=======================================================================
 /// Function that should get the values of the dofs in the object
 //=======================================================================
 void ExplicitTimeSteppableObject::get_dofs(const unsigned& t, 
                                            DoubleVector &dofs) const
 {
  std::ostringstream error_stream;
  error_stream 
   << "Empty default function called.\n"
   << "The function must return the t'th history values of the degrees of \n"
   << "freedom in the object.\n"
   << "Note: It is the responsibility of the object to set the size\n" 
   << "of the vector\n";
  
  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
  
 }
 
 //=====================================================================
 /// Function that sets the values of the dofs in the object
 //====================================================================
 void ExplicitTimeSteppableObject::set_dofs(const DoubleVector &dofs)
 {
  std::ostringstream error_stream;
  error_stream 
   << "Empty default function called.\n"
   << "The function must set the current values of the degrees of \n"
   << "freedom in the object.\n";
  
  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
  
 }
 

 //====================================================================
 ///Function that adds the lambda multiplied by the increment_dofs
 ///vector to the dofs in the object
 //====================================================================
 void ExplicitTimeSteppableObject::add_to_dofs(
  const double &lambda,
  const DoubleVector &increment_dofs)
 {
  std::ostringstream error_stream;
  error_stream 
   << "Empty default function called.\n"
   << "The function must add lambda multiplied by the contents of the\n"
   << "input vector to the degrees of freedom in the object.\n"
   << "Note: It is the responsibility of the object to ensure that the\n"
   << "      the degrees of freedom are in the same order as those \n"
   << "      returned by get_dvaluesdt()\n";
  
  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  
 }

 //==================================================================
 ///\short Virtual function that should be overloaded to return access
 ///to the local time in the object
 //=================================================================
 double& ExplicitTimeSteppableObject::time()
 {
  std::ostringstream error_stream;
  error_stream 
   << "Empty default function called.\n"
   << "The function must return a reference to the local time in the object\n";
  
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);

  return Dummy_time_value;
 }

 /// \short Virtual function that should be overloaded to return a pointer to a
 /// Time object.
 Time* ExplicitTimeSteppableObject::time_pt() const
 {
  std::ostringstream error_stream;
  error_stream 
   << "Empty default function called.\n"
   << "The function must return a pointer to an oomph-lib Time object.\n";
  throw OomphLibError(error_stream.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);

  return 0;
 }


 //================================================================
 /// Euler timestepping x^{t+1} = x^{t} + dt M^{-1} L(x^{t})
 //=================================================================
 void Euler::timestep(ExplicitTimeSteppableObject* const &object_pt,
                      const double &dt)
 {
  object_pt->actions_before_explicit_timestep();
  object_pt->actions_before_explicit_stage();

  //Vector that will hold the inverse mass matrix multiplied by the
  //residuals
  DoubleVector minv_res;
  //Get M^{-1} R
  object_pt->get_dvaluesdt(minv_res);
  
  //Add the result to the unknowns
  object_pt->add_to_dofs(dt,minv_res);
  //Increment the time by the appropriate amount
  object_pt->time() += dt;

  //Call any actions required after the change in the unknowns
  object_pt->actions_after_explicit_stage();
  object_pt->actions_after_explicit_timestep();
 }


//====================================================================
//Broken default timestep function for RungeKutta schemes
//====================================================================
template<unsigned ORDER>
void RungeKutta<ORDER>::timestep(ExplicitTimeSteppableObject* const &object_pt,
                                 const double &dt)
{
 std::ostringstream error_stream;
 error_stream << "Timestep not implemented for order " << ORDER << "\n";
 throw OomphLibError(error_stream.str(),
                     OOMPH_CURRENT_FUNCTION,
                     OOMPH_EXCEPTION_LOCATION);
}

//===================================================================
///Explicit specialisation for fourth-order RK scheme
//==================================================================
template<>
void RungeKutta<4>::timestep(ExplicitTimeSteppableObject* const &object_pt,
                             const double &dt)
{
 object_pt->actions_before_explicit_timestep();

 // Stage 1
 // ============================================================
 object_pt->actions_before_explicit_stage();

 //Store the initial values and initial time
 DoubleVector u;
 object_pt->get_dofs(u);

 //Now get the first unknowns
 DoubleVector k1;
 object_pt->get_dvaluesdt(k1);
 
 //Add to the residuals
 object_pt->add_to_dofs(0.5*dt,k1);
 //Increment the time
 object_pt->time() += 0.5*dt;
 object_pt->actions_after_explicit_stage();


 // Stage 2
 // ============================================================
 object_pt->actions_before_explicit_stage();

 //Get the next unknowns
 DoubleVector k2;
 object_pt->get_dvaluesdt(k2);
 
 //Now reset the residuals
 object_pt->set_dofs(u);
 object_pt->add_to_dofs(0.5*dt,k2);
 //Time remains the same
 object_pt->actions_after_explicit_stage();

 // Stage 3
 // ============================================================
 object_pt->actions_before_explicit_stage();

 //Get the next unknowns
 DoubleVector k3;
 object_pt->get_dvaluesdt(k3);
 
 //Now reset the residuals
 object_pt->set_dofs(u);
 object_pt->add_to_dofs(dt,k3);
 //Increment the time (now at initial_time + dt)
 object_pt->time() += 0.5*dt;
 object_pt->actions_after_explicit_stage();

 // Stage 4
 // ============================================================
 object_pt->actions_before_explicit_stage();

 //Get the final unknowns
 DoubleVector k4;
 object_pt->get_dvaluesdt(k4);
 
 //Set the final values of the unknowns
 object_pt->set_dofs(u);
 object_pt->add_to_dofs((dt/6.0),k1);
 object_pt->add_to_dofs((dt/3.0),k2);
 object_pt->add_to_dofs((dt/3.0),k3);
 object_pt->add_to_dofs((dt/6.0),k4);
 object_pt->actions_after_explicit_stage();

 object_pt->actions_after_explicit_timestep();
}

 //===================================================================
 ///Explicit specialisation for second-order RK scheme
 //==================================================================
 template<>
 void RungeKutta<2>::timestep(ExplicitTimeSteppableObject* const &object_pt,
                              const double &dt)
 {
  object_pt->actions_before_explicit_timestep();

  // Stage 1
  // ============================================================
  object_pt->actions_before_explicit_stage();

  // Store the initial values
  DoubleVector u;
  object_pt->get_dofs(u);

  // Get f1 (time derivative at t0, y0) and add to dofs
  DoubleVector f1;
  object_pt->get_dvaluesdt(f1);
  object_pt->add_to_dofs(dt, f1);

  // Advance time to t1 = t0 + dt
  object_pt->time() += dt;

  object_pt->actions_after_explicit_stage();


  // Stage 2
  // ============================================================
  object_pt->actions_before_explicit_stage();

  // get f2 (with t=t1, y = y0 + h f0)
  DoubleVector f2;
  object_pt->get_dvaluesdt(f2);

  // Final answer is starting dofs + h/2 * (f1 + f2)
  object_pt->set_dofs(u);
  object_pt->add_to_dofs(0.5*dt, f1);
  object_pt->add_to_dofs(0.5*dt, f2);

  object_pt->actions_after_explicit_stage();

  // Done, do actions
  object_pt->actions_after_explicit_timestep();
 }


//=================================================================
//General constructor for LowOrder RK schemes
//==================================================================
template<unsigned ORDER>
LowStorageRungeKutta<ORDER>::LowStorageRungeKutta()
{
 Type = "LowStorageRungeKutta";
}

//======================================================================
///Broken default timestep for LowStorageRungeKutta
//======================================================================
template<unsigned ORDER>
void LowStorageRungeKutta<ORDER>::timestep(
 ExplicitTimeSteppableObject* const &object_pt,
 const double &dt)
{
 std::ostringstream error_stream;
 error_stream << "Timestep not implemented for order " << ORDER << "\n";
 throw OomphLibError(error_stream.str(),
                     OOMPH_CURRENT_FUNCTION,
                     OOMPH_EXCEPTION_LOCATION);
}

//================================================================
//Specialised constructor for fourth-order RK scheme
//================================================================
template<>
LowStorageRungeKutta<4>::LowStorageRungeKutta()
{
 Type="LowStorageRungeKutta";

 A.resize(5);
 A[0] = 0.0;
 A[1] = - 567301805773.0/1357537059087.0;
 A[2] = - 2404267990393.0/2016746695238.0;
 A[3] = - 3550918686646.0/2091501179385.0;
 A[4] = - 1275806237668.0/842570457699.0;

 B.resize(5);
 B[0] = 1432997174477.0/9575080441755.0;
 B[1] = 5161836677717.0/13612068292357.0;
 B[2] = 1720146321549.0/2090206949498.0;
 B[3] = 3134564353537.0/4481467310338.0;
 B[4] = 2277821191437.0/14882151754819.0;

 C.resize(5);
 C[0] = B[0];
 C[1] = 2526269341429.0/6820363962896.0;
 C[2] = 2006345519317.0/3224310063776.0;
 C[3] = 2802321613138.0/2924317926251.0;
 C[4] = 1.0;
}



//Explicit specialisation for fourth-order RK scheme
template<>
void LowStorageRungeKutta<4>::timestep(
 ExplicitTimeSteppableObject* const &object_pt,
 const double &dt)
{
 object_pt->actions_before_explicit_timestep();

 //Store the initial time
 const double initial_time = object_pt->time();
 //Temporary storage
 DoubleVector k;

 //Temporary storage for the inverse mass matrix multiplied by the residuals
 DoubleVector minv_res;

 //Loop over the number of steps in the scheme
 for(unsigned i=0;i<5;i++)
  {
   object_pt->actions_before_explicit_stage();

   //Get the inverse mass matrix multiplied by the current value
   //of the residuals
   object_pt->get_dvaluesdt(minv_res);
   //Get the values of k
   const unsigned n_dof= minv_res.nrow();

   //First time round resize k and initialise to zero
   if(i==0) {k.build(minv_res.distribution_pt(),0.0);}
   //Now construct the next value of k
   for(unsigned n=0;n<n_dof;n++)
    {
     k[n] *= A[i];
     k[n] += dt*minv_res[n];
    }

   //Add to the residuals
   object_pt->add_to_dofs(B[i],k);
   //Set the new time
   object_pt->time() = initial_time + C[i]*dt;
   
   object_pt->actions_after_explicit_stage();
  }

 object_pt->actions_after_explicit_timestep();
}

// ??ds this could be heavily optimised if needed. Keeping it simple for
// now
void EBDF3::timestep(ExplicitTimeSteppableObject* const &object_pt,
                     const double &dt)
{
 using namespace StringConversion;

 object_pt->actions_before_explicit_timestep();
 object_pt->actions_before_explicit_stage();

 // Storage indicies for the history values that we need
 unsigned tn = 1;
 unsigned tnm1 = tn+1;
 unsigned tnm2 = tnm1+1;

 // Check dts are the same, this will need to be removed if ebdf3 is being
 // used as something other than a predictor... But seeing as it isn't
 // stable that isn't likely.
#ifdef PARANOID
 if(std::abs(dt - object_pt->time_pt()->dt(0)) > 1e-15)
  {
   std::string err = "Inconsistent dts! Predictor is stepping by " + to_string(dt);
   err += " but dt(0) = " + to_string(object_pt->time_pt()->dt(0));
   throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                       OOMPH_CURRENT_FUNCTION);
  }
#endif

 // Get older dt values
 double dtnm1 = object_pt->time_pt()->dt(1);
 double dtnm2 = object_pt->time_pt()->dt(2);

 // Calculate weights for these dts
 set_weights(dt, dtnm1, dtnm2);
 
 // Get derivative value at step n (even though this uses values from t=0 =
 // step n+1, it's ok because we haven't changed the values in that slot
 // yet).
 DoubleVector fn;
 object_pt->get_dvaluesdt(fn);
 fn *= Fn_weight;

 // Extract history values and multiply by their weights
 DoubleVector ynp1, yn, ynm1, ynm2;
 object_pt->get_dofs(tn, yn);
 yn *= Yn_weight;
 object_pt->get_dofs(tnm1, ynm1);
 ynm1 *= Ynm1_weight;
 object_pt->get_dofs(tnm2, ynm2);
 ynm2 *= Ynm2_weight;


 // Add everything together
 ynp1 = yn;
 ynp1 += ynm1;
 ynp1 += ynm2;
 ynp1 += fn;

 // Done, update things in the object
 object_pt->set_dofs(ynp1);
 object_pt->time() += dt;

 // Actions functions
 object_pt->actions_after_explicit_stage();
 object_pt->actions_after_explicit_timestep();
}

/// Calculate the weights for this set of step sizes.
void EBDF3::set_weights(const double &dtn, const double &dtnm1, 
                        const double &dtnm2)
 {
  using namespace std;

  // If this is slow we can probably optimise by doing direct
  // multiplication instead of using pow.

  // Generated using sympy from my python ode code.

  double denom = pow(dtnm1,4)*dtnm2 + 2*pow(dtnm1,3)*pow(dtnm2,2) 
   + pow(dtnm1,2)*pow(dtnm2,3);

  Yn_weight = -(2*pow(dtn,3)*dtnm1*dtnm2 + pow(dtn,3)*pow(dtnm2,2) 
               + 3*pow(dtn,2)*pow(dtnm1,2)*dtnm2 
               + 3*pow(dtn,2)*dtnm1*pow(dtnm2,2) + pow(dtn,2)*pow(dtnm2,3) 
               - pow(dtnm1,4)*dtnm2 - 2*pow(dtnm1,3)*pow(dtnm2,2) 
               - pow(dtnm1,2)*pow(dtnm2,3))/denom;

  Ynm1_weight = -(-pow(dtn,3)*pow(dtnm1,2) - 2*pow(dtn,3)*dtnm1*dtnm2 
                 - pow(dtn,3)*pow(dtnm2,2) - pow(dtn,2)*pow(dtnm1,3) 
                 - 3*pow(dtn,2)*pow(dtnm1,2)*dtnm2 
                 - 3*pow(dtn,2)*dtnm1*pow(dtnm2,2) 
                 - pow(dtn,2)*pow(dtnm2,3))/denom;

  Ynm2_weight = -(pow(dtn,3)*pow(dtnm1,2) + pow(dtn,2)*pow(dtnm1,3))/denom;

  Fn_weight = -(-pow(dtn,3)*pow(dtnm1,2)*dtnm2 - pow(dtn,3)*dtnm1*pow(dtnm2,2) 
               - 2*pow(dtn,2)*pow(dtnm1,3)*dtnm2 
               - 3*pow(dtn,2)*pow(dtnm1,2)*pow(dtnm2,2) 
               - pow(dtn,2)*dtnm1*pow(dtnm2,3) - dtn*pow(dtnm1,4)*dtnm2 
               - 2*dtn*pow(dtnm1,3)*pow(dtnm2,2) 
               - dtn*pow(dtnm1,2)*pow(dtnm2,3))/denom;
 }



//Force build of templates
template class RungeKutta<4>;
template class LowStorageRungeKutta<4>;
}