numerics/w_method.h

// -*- c++ -*-

// +--------------------------------------------------------------------+
// | This file is part of MathTL - the Mathematical Template Library    |
// |                                                                    |
// | Copyright (c) 2002-2009                                            |
// | Thorsten Raasch, Manuel Werner                                     |
// +--------------------------------------------------------------------+

#ifndef _MATHTL_W_METHOD_H
#define _MATHTL_W_METHOD_H

#include <algebra/triangular_matrix.h>
#include <algebra/vector.h>
#include <numerics/ivp.h>
#include <numerics/one_step_scheme.h>

namespace MathTL
{
  template <class VECTOR>
  class WMethodStageEquationHelper
  {
  public:
    virtual ~WMethodStageEquationHelper() = 0;
    
    virtual void solve_W_stage_equation(const AbstractIVP<VECTOR>* ivp,
                                        const double t,
                                        const VECTOR& v,
                                        const double alpha,
                                        const VECTOR& y,
                                        const double tolerance,
                                        VECTOR& x) const = 0;

    virtual void approximate_ft(const AbstractIVP<VECTOR>* ivp,
                                const double t,
                                const VECTOR& v,
                                const double tolerance,
                                VECTOR& result) const = 0;
  };

  template <class VECTOR>
  class WMethodPreprocessRHSHelper
  {
  public:
    virtual ~WMethodPreprocessRHSHelper() = 0;

    virtual void preprocess_rhs_share(VECTOR& wbeforeandafter,
                                      const double tolerance) const
    { // default behaviour: do nothing
    }
  };
  
  template <class VECTOR>
  class WMethod
    : public OneStepScheme<VECTOR>
  {
  public:
    enum Method {
      GRK4T,     // s=4, p=4
      RODAS3,    // s=4, p=3, index 1, L-stable, stiffly accurate
      ROS2,      // s=2, p=2, L-stable
      ROS3,      // s=3, p=3, L-stable
      ROWDA3,    // s=3, p=3, index 1, L-stable
      ROS3P,     // s=3, p=3, index 1, nonlinear PDEs
      ROS3Pw,    // s=3, p=3, index 1, PDEs
      ROSI2P2,   // s=4, p=3, index 2, L-stable, PDEs, stiffly accurate
      ROSI2PW,   // s=4, p=3, index 2, L-stable, PDEs, stiffly accurate, W-method
      RODAS,     // s=6, p=4
      RODASP     // s=6, p=4, index 1, PDEs, L-stable, stiffly accurate
    };

    WMethod(const Method method,
            const WMethodStageEquationHelper<VECTOR>* stage_equation_solver);
    
    virtual ~WMethod() {}

    void increment(const AbstractIVP<VECTOR>* ivp,
                   const double t_m,
                   const VECTOR& u_m,
                   const double tau,
                   VECTOR& u_mplus1,
                   VECTOR& error_estimate,
                   const double tolerance = 1e-2) const;

    int order() const { return p; }

    void set_preprocessor(WMethodPreprocessRHSHelper<VECTOR>* pp) {
      preprocessor = pp;
    }
    
    static void transform_coefficients(const LowerTriangularMatrix<double>& Alpha,
                                       const LowerTriangularMatrix<double>& Gamma,
                                       const Vector<double>& b,
                                       const Vector<double>& bhat,
                                       LowerTriangularMatrix<double>& A,
                                       LowerTriangularMatrix<double>& C,
                                       Vector<double>& m,
                                       Vector<double>& e,
                                       Vector<double>& alpha_vector,
                                       Vector<double>& gamma_vector);

    static void check_order_conditions(const LowerTriangularMatrix<double>& Alpha,
                                       const LowerTriangularMatrix<double>& Gamma,
                                       const Vector<double>& b,
                                       const Vector<double>& bhat,
                                       const bool wmethod=false);

    LowerTriangularMatrix<double> A;

    LowerTriangularMatrix<double> C;

    Vector<double> m;

    Vector<double> e;

    Vector<double> alpha_vector;

    Vector<double> gamma_vector;


  protected:
    const WMethodStageEquationHelper<VECTOR>* stage_equation_helper;

    WMethodPreprocessRHSHelper<VECTOR>* preprocessor;

    int p;
  };
}

#include <numerics/w_method.cpp>

#endif