TimeIntegrationSchemeFIT.h

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: TimeIntegrationSchemeFIT.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Header file of time integration scheme FIT base class
//
///////////////////////////////////////////////////////////////////////////////
 
// Note: The file is named TimeIntegrationSchemeFIT to parallel the
// TimeIntegrationSchemeGLM file but the class is named
// FractionalInTimeIntegrationScheme so keep with the factory naming
// convention.
 
#ifndef NEKTAR_LIB_UTILITIES_TIME_INTEGRATION_TIME_INTEGRATION_SCHEME_FIT
#define NEKTAR_LIB_UTILITIES_TIME_INTEGRATION_TIME_INTEGRATION_SCHEME_FIT
 
#define LUE LIB_UTILITIES_EXPORT
 
#include <string>
 
#include <LibUtilities/TimeIntegration/TimeIntegrationScheme.h>
 
namespace Nektar
{
namespace LibUtilities
{
 
///////////////////////////////////////////////////////////////////////////////
/// Class for fractional-in-time integration.
class FractionalInTimeIntegrationScheme : public TimeIntegrationScheme
{
public:
    /// Constructor
    FractionalInTimeIntegrationScheme(std::string variant, unsigned int order,
                                      std::vector<NekDouble> freeParams);
 
    /// Destructor
    virtual ~FractionalInTimeIntegrationScheme()
    {
    }
 
    /// Creator
    static TimeIntegrationSchemeSharedPtr create(
        std::string variant, unsigned int order,
        std::vector<NekDouble> freeParams)
    {
        TimeIntegrationSchemeSharedPtr p =
            MemoryManager<FractionalInTimeIntegrationScheme>::AllocateSharedPtr(
                variant, order, freeParams);
 
        return p;
    }
 
    static std::string className;
 
    // Access methods from the base class that are virtual
    LUE virtual std::string GetName() const
    {
        return m_name;
    }
 
    LUE virtual std::string GetVariant() const
    {
        return m_variant;
    }
 
    LUE virtual unsigned int GetOrder() const
    {
        return m_order;
    }
 
    LUE virtual std::vector<NekDouble> GetFreeParams() const
    {
        return m_freeParams;
    }
 
    LUE virtual TimeIntegrationSchemeType GetIntegrationSchemeType() const
    {
        return m_schemeType;
    }
 
    LUE virtual NekDouble GetTimeStability() const
    {
        return 1.0;
    }
 
    LUE unsigned int GetNumIntegrationPhases() const
    {
        return 1;
    }
 
    /**
     * \brief Gets the solution vector of the ODE
     */
    inline const TripleArray &GetSolutionVector() const
    {
        return m_u;
    }
 
    /**
     * \brief Sets the solution vector of the ODE
     */
    inline void SetSolutionVector(const int Offset, const DoubleArray &y)
    {
        m_u[Offset] = y;
    }
 
    // The worker methods from the base class that are virtual
    LUE virtual void InitializeScheme(const NekDouble deltaT,
                                      ConstDoubleArray &y_0,
                                      const NekDouble time,
                                      const TimeIntegrationSchemeOperators &op);
 
    LUE virtual ConstDoubleArray &TimeIntegrate(
        const int timestep, const NekDouble delta_t,
        const TimeIntegrationSchemeOperators &op);
 
    LUE virtual void print(std::ostream &os) const;
    LUE virtual void printFull(std::ostream &os) const;
 
    // Friend classes
    LUE friend std::ostream &operator<<(
        std::ostream &os, const FractionalInTimeIntegrationScheme &rhs);
    LUE friend std::ostream &operator<<(
        std::ostream &os,
        const FractionalInTimeIntegrationSchemeSharedPtr &rhs);
 
protected:
    struct Instance
    {
        int base;
 
        int index;         // Index of this instance
        bool active;       // Used to determine if active
        int activecounter; // counter used to flip active bit
        int activebase;
 
        // Major storage for auxilliary ODE solutions.
        // Storage for values of y currently used to update u
        ComplexTripleArray stage_y;
        std::pair<int, int> stage_ind; // Time-step counters indicating the
                                       // interval ymain is associated with
 
        // Staging allocation
        bool stage_active;
        int stage_ccounter;
        int stage_cbase; // This base is halved after the first cycle
        int stage_fcounter;
        int stage_fbase; // This base is halved after the first cycle
 
        // Ceiling stash allocation
        int cstash_counter; // Counter used to determine
                            // when to stash
        int cstash_base;    // base for counter
        ComplexTripleArray cstash_y;
        std::pair<int, int> cstash_ind; // ind(1) is never used:
                                        // it always matches main.ind(1)
 
        // Ceiling sandbox allocation
        bool csandbox_active; // Flag to determine when
                              // stash 2 is utilized
        int csandbox_counter;
        ComplexTripleArray csandbox_y;
        std::pair<int, int> csandbox_ind;
 
        // Floor stash
        int fstash_base;
        ComplexTripleArray fstash_y;
        std::pair<int, int> fstash_ind;
 
        // Floor sandbox
        bool fsandbox_active;
        int fsandbox_activebase;
        int fsandbox_stashincrement;
        ComplexTripleArray fsandbox_y;
        std::pair<int, int> fsandbox_ind;
 
        // Talbot quadrature rule
        ComplexSingleArray z;
        ComplexSingleArray w;
 
        TripleArray As;
 
        ComplexSingleArray E;
        ComplexDoubleArray Eh;
        ComplexDoubleArray AtEh;
    };
 
    inline unsigned int modIncrement(const unsigned int counter,
                                     const unsigned int base) const;
 
    inline unsigned int computeL(const unsigned int base,
                                 const unsigned int m) const;
 
    inline unsigned int computeQML(const unsigned int base,
                                   const unsigned int m);
 
    inline unsigned int computeTaus(const unsigned int base,
                                    const unsigned int m);
 
    void talbotQuadrature(const unsigned int nQuadPts, const NekDouble mu,
                          const NekDouble nu, const NekDouble sigma,
                          ComplexSingleArray &lamb,
                          ComplexSingleArray &w) const;
 
    void integralClassInitialize(const unsigned int index,
                                 Instance &instance) const;
 
    void updateStage(const unsigned int timeStep, Instance &instance);
 
    void finalIncrement(const unsigned int timeStep,
                        const TimeIntegrationSchemeOperators &op);
 
    void integralContribution(const unsigned int timeStep,
                              const unsigned int tauml,
                              const Instance &instance);
 
    void timeAdvance(const unsigned int timeStep,
                     const TimeIntegrationSchemeOperators &op,
                     Instance &instance, ComplexTripleArray &y);
 
    void advanceSandbox(const unsigned int timeStep,
                        const TimeIntegrationSchemeOperators &op,
                        Instance &instance);
 
    // Variables common to all schemes.
    std::string m_name;
    std::string m_variant;
    unsigned int m_order{0};
    std::vector<NekDouble> m_freeParams;
 
    TimeIntegrationSchemeType m_schemeType{eFractionalInTime};
 
    // Varaibles and methods specific to FIT integration schemes.
    NekDouble m_deltaT{0};
    NekDouble m_T{0};            // Finial time
    unsigned int m_maxTimeSteps; // Number of time steps.
    NekDouble m_alpha{0.3};      // Value for exp integration.
    unsigned int m_base{4};      // "Base" of the algorithm.
    unsigned int m_nQuadPts{20}; // Number of Talbot quadrature rule points
    NekDouble m_sigma{0};
    NekDouble m_mu0{8};
    NekDouble m_nu{0.6};
 
    int m_nvars{0};   // Number of variables in the integration scheme.
    int m_npoints{0}; // Number of points    in the integration scheme.
 
    unsigned int m_Lmax{0}; // Maxium number of integral groups.
    Array<OneD, Instance> m_integral_classes;
 
    // Demarcation integers
    Array<OneD, int> m_qml;
    // Demarcation interval markers
    Array<OneD, int> m_taus;
 
    // Storage of the initial values.
    DoubleArray m_u0;
    // Storage of the next solution from the final increment.
    DoubleArray m_uNext;
    // Storage for the integral contribution.
    ComplexDoubleArray m_uInt;
    // Storage for the exponential factor in the integral contribution.
    ComplexSingleArray m_expFactor;
 
    // Storage of previous states and associated timesteps.
    TripleArray m_u;
 
    // Storage for the stage derivative as the data will be re-used to
    // update the solution.
    DoubleArray m_F;
 
    // J
    SingleArray m_J;
 
    // Ahat array one for each order.
    TripleArray m_Ahats;
 
    // Mulitply the last Ahat array, transposed by J
    SingleArray m_AhattJ;
 
}; // end class FractionalInTimeIntegrator
 
LUE std::ostream &operator<<(std::ostream &os,
                             const FractionalInTimeIntegrationScheme &rhs);
LUE std::ostream &operator<<(
    std::ostream &os, const FractionalInTimeIntegrationSchemeSharedPtr &rhs);
 
} // end namespace LibUtilities
} // end namespace Nektar
 
#endif