TimeIntegrationSchemeFIT.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: TimeIntegrationSchemeFIT.cpp
//
// 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: implementation of time integration scheme FIT 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.
 
#include <LibUtilities/TimeIntegration/TimeIntegrationSchemeFIT.h>
 
namespace Nektar
{
namespace LibUtilities
{
/**
 * @class FractionalInTimeIntegrationScheme
 *
 * A fast convolution algorithm for computing solutions to (Caputo)
 * time-fractional differential equations. This is an explicit solver
 * that expresses the solution as an integral over a Talbot curve,
 * which is discretized with quadrature. First-order quadrature is
 * currently implemented (Soon be expanded to forth order).
 */
FractionalInTimeIntegrationScheme::FractionalInTimeIntegrationScheme(
    std::string variant, unsigned int order, std::vector<NekDouble> freeParams)
    : TimeIntegrationScheme(variant, order, freeParams),
      m_name("FractionalInTime")
{
    m_variant    = variant;
    m_order      = order;
    m_freeParams = freeParams;
 
    // Currently up to 4th order is implemented.
    ASSERTL1(0 < order && order <= 4,
             "FractionalInTime Time integration scheme bad order: " +
                 std::to_string(order));
 
    ASSERTL1(freeParams.size() == 0 ||     // Use defaults
                 freeParams.size() == 1 || // Alpha
                 freeParams.size() == 2 || // Base
                 freeParams.size() == 6,   // Talbot quadrature rule
             "FractionalInTime Time integration scheme invalid number "
             "of free parameters, expected zero, one <alpha>, "
             "two <alpha, base>, or "
             "six <alpha, base, nQuadPts, sigma, mu0, nu> received  " +
                 std::to_string(freeParams.size()));
 
    if (freeParams.size() >= 1)
    {
        m_alpha = freeParams[0]; // Value for exp integration.
    }
 
    if (freeParams.size() >= 2)
    {
        m_base = freeParams[1]; // "Base" of the algorithm.
    }
 
    if (freeParams.size() == 6)
    {
        m_nQuadPts = freeParams[2]; // Number of Talbot quadrature rule points
        m_sigma    = freeParams[3];
        m_mu0      = freeParams[4];
        m_nu       = freeParams[5];
    }
}
 
/**
 * @brief Worker method to initialize the integration scheme.
 */
void FractionalInTimeIntegrationScheme::v_InitializeScheme(
    const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time,
    const TimeIntegrationSchemeOperators &op)
 
{
    boost::ignore_unused(op);
 
    m_nvars   = y_0.size();
    m_npoints = y_0[0].size();
 
    m_deltaT = deltaT;
 
    m_T            = time; // Finial time;
    m_maxTimeSteps = m_T / m_deltaT;
 
    // The +2 below is a buffer, and keeps +2 extra rectangle groups
    // in case T needs to be increased later.
    m_Lmax = computeL(m_base, m_maxTimeSteps) + 2;
 
    // Demarcation integers - one array that is re-used
    m_qml = Array<OneD, int>(m_Lmax - 1, 0);
    // Demarcation interval markers - one array that is re-used
    m_taus = Array<OneD, int>(m_Lmax + 1, 0);
 
    // Storage of the initial values.
    m_u0 = y_0;
 
    // Storage for the exponential factor in the integral
    // contribution. One array that is re-used
    m_expFactor = ComplexSingleArray(m_nQuadPts, 0.0);
 
    // Storage of previous states and associated timesteps.
    m_u = TripleArray(m_order + 1);
 
    for (unsigned int m = 0; m <= m_order; ++m)
    {
        m_u[m] = DoubleArray(m_nvars);
 
        for (unsigned int i = 0; i < m_nvars; ++i)
        {
            m_u[m][i] = SingleArray(m_npoints, 0.0);
 
            for (unsigned int j = 0; j < m_npoints; ++j)
            {
                // Store the initial values as the first previous state.
                if (m == 0)
                    m_u[m][i][j] = m_u0[i][j];
                else
                    m_u[m][i][j] = 0;
            }
        }
    }
 
    // Storage for the stage derivative as the data will be re-used to
    // update the solution.
    m_F = DoubleArray(m_nvars);
    // Storage of the next solution from the final increment.
    m_uNext = DoubleArray(m_nvars);
    // Storage for the integral contribution.
    m_uInt = ComplexDoubleArray(m_nvars);
 
    for (unsigned int i = 0; i < m_nvars; ++i)
    {
        m_F[i]     = SingleArray(m_npoints, 0.0);
        m_uNext[i] = SingleArray(m_npoints, 0.0);
        m_uInt[i]  = ComplexSingleArray(m_npoints, 0.0);
    }
 
    // J
    m_J = SingleArray(m_order, 0.0);
 
    m_J[0] = pow(m_deltaT, m_alpha) / tgamma(m_alpha + 1.);
 
    for (unsigned int m = 1, m_1 = 0; m < m_order; ++m, ++m_1)
    {
        m_J[m] = m_J[m_1] * NekDouble(m) / (NekDouble(m) + m_alpha);
    }
 
    // Ahat array, one for each order.
    // These are elements in a multi-step exponential integrator tableau
    m_Ahats = TripleArray(m_order + 1);
 
    for (unsigned int m = 1; m <= m_order; ++m)
    {
        m_Ahats[m] = DoubleArray(m);
 
        for (unsigned int n = 0; n < m; ++n)
        {
            m_Ahats[m][n] = SingleArray(m, 0.0);
        }
 
        switch (m)
        {
            case 1:
                m_Ahats[m][0][0] = 1.;
                break;
 
            case 2:
                m_Ahats[m][0][0] = 1.;
                m_Ahats[m][0][1] = 0.;
                m_Ahats[m][1][0] = 1.;
                m_Ahats[m][1][1] = -1.;
                break;
 
            case 3:
                m_Ahats[m][0][0] = 1.;
                m_Ahats[m][0][1] = 0.;
                m_Ahats[m][0][2] = 0;
                m_Ahats[m][1][0] = 3. / 2.;
                m_Ahats[m][1][1] = -2.;
                m_Ahats[m][1][2] = 1. / 2.;
                m_Ahats[m][2][0] = 1. / 2.;
                m_Ahats[m][2][1] = -1.;
                m_Ahats[m][2][2] = 1. / 2.;
                break;
 
            case 4:
                m_Ahats[m][0][0] = 1.;
                m_Ahats[m][0][1] = 0.;
                m_Ahats[m][0][2] = 0.;
                m_Ahats[m][0][3] = 0.;
 
                m_Ahats[m][1][0] = 11. / 6.;
                m_Ahats[m][1][1] = -3;
                m_Ahats[m][1][2] = 3. / 2.;
                m_Ahats[m][1][3] = -1. / 3.;
 
                m_Ahats[m][2][0] = 1.;
                m_Ahats[m][2][1] = -5. / 2.;
                m_Ahats[m][2][2] = 2.;
                m_Ahats[m][2][3] = -1. / 2.;
 
                m_Ahats[m][3][0] = 1. / 6.;
                m_Ahats[m][3][1] = -1. / 2.;
                m_Ahats[m][3][2] = 1. / 2.;
                m_Ahats[m][3][3] = -1. / 6.;
                break;
 
            default:
 
                m_Ahats[m][0][0] = 1;
 
                for (unsigned int j = 2; j <= m; ++j)
                {
                    for (unsigned int i = 0; i < m; ++i)
                    {
                        m_Ahats[m][j - 1][i] = pow((1 - j), i);
                    }
                }
 
                ASSERTL1(false, "No matrix inverse.");
 
                // Future code: m_Ahats[m] = inv(m_Ahats[m]);
 
                break;
        }
    }
 
    // Mulitply the last Ahat array, transposed, by J
    m_AhattJ = SingleArray(m_order, 0.0);
 
    for (unsigned int i = 0; i < m_order; ++i)
    {
        for (unsigned int j = 0; j < m_order; ++j)
        {
            m_AhattJ[i] += m_Ahats[m_order][j][i] * m_J[j];
        }
    }
 
    m_integral_classes = Array<OneD, Instance>(m_Lmax);
 
    for (int l = 0; l < m_Lmax; ++l)
    {
        integralClassInitialize(l + 1, m_integral_classes[l]);
    }
}
 
/**
 * @brief Worker method that performs the time integration.
 */
ConstDoubleArray &FractionalInTimeIntegrationScheme::v_TimeIntegrate(
    const int timestep, const NekDouble delta_t,
    const TimeIntegrationSchemeOperators &op)
{
    boost::ignore_unused(delta_t);
 
    ASSERTL1(delta_t == m_deltaT,
             "Delta T has changed which is not permitted.");
 
    // The Fractional in Time works via the logical? time step value.
    int timeStep = timestep + 1;
 
    // Update the storage and counters for integral classes.  Performs
    // staging for updating u.
    for (int l = 0; l < m_Lmax; ++l)
    {
        updateStage(timeStep, m_integral_classes[l]);
    }
 
    // Compute u update to time timeStep * m_deltaT.  Stored in
    // m_uNext.
    finalIncrement(timeStep, op);
 
    // Contributions to the current integral
    int L = computeTaus(m_base, timeStep);
 
    for (int l = 0; l < L; ++l)
    {
        // Integral contribution over [taus(i+1) taus(i)]. Stored in
        // m_uInt.
        integralContribution(timeStep, m_taus[l], m_integral_classes[l]);
 
        for (int i = 0; i < m_nvars; ++i)
        {
            for (int j = 0; j < m_npoints; ++j)
            {
                m_uNext[i][j] += m_uInt[i][j].real();
            }
        }
    }
 
    // Shuffle the previous solutions back one in the history.
    for (int m = m_order; m > 0; --m)
    {
        for (int i = 0; i < m_nvars; ++i)
        {
            for (int j = 0; j < m_npoints; ++j)
            {
                m_u[m][i][j] = m_u[m - 1][i][j];
            }
        }
    }
 
    // Get the current solution.
    for (int i = 0; i < m_nvars; ++i)
    {
        for (int j = 0; j < m_npoints; ++j)
        {
            m_u[0][i][j] = m_uNext[i][j] + m_u0[i][j];
 
            m_uNext[i][j] = 0; // Zero out for the next itereation.
        }
    }
 
    // Dump the current solution.
    // std::cout << "timeStep  " << timeStep << std::endl;
    // for( int j=0; j<m_npoints; ++j )
    // {
    //     for( int i=0; i<m_nvars; ++i )
    //     {
    //         for( int m=0; m<m_order+1; ++m )
    //         {
    //             std::cout << m_u[m][i][j] << "  ";
    //         }
    //         std::cout << std::endl;
    //     }
    //     std::cout << std::endl;
    // }
    // std::cout << std::endl;
 
    // Update the storage and counters for integral classes to
    // time timeStep * m_deltaT. Also time-steps the sandboxes and stashes.
    for (int i = 0; i < m_Lmax; ++i)
    {
        advanceSandbox(timeStep, op, m_integral_classes[i]);
    }
 
    return m_u[0];
}
 
/**
 * @brief Method that increments the counter then performs mod
 * calculation.
 */
unsigned int FractionalInTimeIntegrationScheme::modIncrement(
    const int unsigned counter, const int unsigned base) const
{
    return (counter + 1) % base;
}
 
/**
 * @brief Method to compute the smallest integer L such that base < 2
 * * base^l.
 */
unsigned int FractionalInTimeIntegrationScheme::computeL(
    const unsigned int base, const unsigned int l) const
{
    unsigned int L = ceil(log(l / 2.0) / log(base));
 
    if (l % (unsigned int)(2 * pow(base, L)) == 0)
    {
        ++L;
    }
 
    return L;
}
 
/**
 * @brief Method to compute the demarcation integers q_{m, ell}.
 *
 *  Returns a length-(L-1) vector qml such that h*taus are interval
 *  boundaries for a partition of [0, m h]. The value of h is not
 *  needed to compute this vector.
 */
unsigned int FractionalInTimeIntegrationScheme::computeQML(
    const unsigned int base, const unsigned int m)
{
    int L = computeL(base, m);
 
    // m_qml is set in InitializeScheme to be the largest length expected.
    // qml = Array<OneD, int>( L-1, 0 );
 
    for (unsigned int i = 0; i < L - 1; ++i)
    {
        m_qml[i] = floor(m / pow(base, i + 1)) - 1;
    }
 
    return L;
}
 
/**
 * @brief Method to compute the demarcation interval marker tau_{m, ell}.
 *
 * Returns a length-(L+1) vector tauml such that h*taus are interval
 * boundaries for a partition of [0, m h]. The value of h is not
 * needed to compute this vector.
 */
unsigned int FractionalInTimeIntegrationScheme::computeTaus(
    const unsigned int base, const unsigned int m)
{
    if (m == 1)
    {
        m_taus[0] = 0;
 
        return 0;
    }
    else
    {
        unsigned int L = computeQML(base, m);
 
        // m_taus is set in InitializeScheme to be the largest length
        // expected.
 
        m_taus[0] = m - 1;
 
        for (unsigned int i = 1; i < L; ++i)
        {
            m_taus[i] = m_qml[i - 1] * pow(base, i);
        }
 
        m_taus[L] = 0;
 
        return L;
    }
}
 
/**
 * @brief Method to compute the quadrature rule over Tablot contour
 *
 * Returns a quadrature rule over the Tablot contour defined by the
 * parameterization.
 *
 * gamma(th) = sigma + mu * ( th*cot(th) + i*nu*th ),  -pi < th < pi
 *
 * An N-point rule is returned, equidistant in the parameter theta. The
 * returned quadrature rule approximes an integral over the contour.
 */
void FractionalInTimeIntegrationScheme::talbotQuadrature(
    const unsigned int nQuadPts, const NekDouble mu, const NekDouble nu,
    const NekDouble sigma, ComplexSingleArray &lamb,
    ComplexSingleArray &w) const
{
    lamb = ComplexSingleArray(nQuadPts, 0.0);
    w    = ComplexSingleArray(nQuadPts, 0.0);
 
    for (unsigned int q = 0; q < nQuadPts; ++q)
    {
        NekDouble th =
            (NekDouble(q) + 0.5) / NekDouble(nQuadPts) * 2.0 * M_PI - M_PI;
 
        lamb[q] = sigma + mu * th * std::complex<NekDouble>(1. / tan(th), nu);
 
        w[q] = std::complex<NekDouble>(0, -1. / NekDouble(nQuadPts)) * mu *
               std::complex<NekDouble>(1. / tan(th) - th / (sin(th) * sin(th)),
                                       nu);
    }
 
    // Special case for th = 0 which happens when there is an odd
    // number of quadrature points.
    if (nQuadPts % 2 == 1)
    {
        unsigned int q = (nQuadPts + 1) / 2;
 
        lamb[q] = std::complex<NekDouble>(sigma + mu, 0);
 
        w[q] = std::complex<NekDouble>(nu * mu / nQuadPts, 0);
    }
}
 
/**
 * @brief Method to initialize the integral class
 */
void FractionalInTimeIntegrationScheme::integralClassInitialize(
    const unsigned int index, Instance &instance) const
{
    /**
     * /brief
     *
     * This object stores information for performing integration over
     * an interval [a, b]. (Defined by taus in the parent calling
     * function.)
     *
     * The "main" object stores information about [a,b]. In
     * particular, main.ind identifies [a,b] via multiples of h.
     *
     * Periodically the values of [a,b] need to be incremented. The
     * necessary background storage to accomplish this increment
     * depends whether a or b is being incremented.
     *
     * The objects with "f" ("Floor") modifiers are associated with
     * increments of the interval floor a.
     *
     * The objects with "c" ("Ceiling") modifiers are associated with
     * increments of the interval ceiling b.
     *
     * Items on the "stage" are stored for use in computing u at the
     * current time.  Items in the "stash" are stored for use for
     * future staging. Items in the "sandbox" are being actively
     * updated at the current time for future stashing. Only items in
     * the sandbox are time-stepped. the stage and stash locations are
     * for storage only.
     *
     * This is the same for all integral classes, so there's probably
     * a better way to engineer this. And technically, all that's
     * needed is the array K(instance.z) anyway.
     */
 
    instance.base          = m_base;
    instance.index         = index; // Index of this instance
    instance.active        = false; // Used to determine if active
    instance.activecounter = 0;     // Counter used to flip active bit
    instance.activebase    = 2. * pow(m_base, (index - 1));
 
    // Storage for values of y currently used to update u
    instance.stage_y    = ComplexTripleArray(m_nvars);
    instance.cstash_y   = ComplexTripleArray(m_nvars);
    instance.csandbox_y = ComplexTripleArray(m_nvars);
    instance.fstash_y   = ComplexTripleArray(m_nvars);
    instance.fsandbox_y = ComplexTripleArray(m_nvars);
 
    for (unsigned int q = 0; q < m_nvars; ++q)
    {
        instance.stage_y[q]    = ComplexDoubleArray(m_npoints);
        instance.cstash_y[q]   = ComplexDoubleArray(m_npoints);
        instance.csandbox_y[q] = ComplexDoubleArray(m_npoints);
        instance.fstash_y[q]   = ComplexDoubleArray(m_npoints);
        instance.fsandbox_y[q] = ComplexDoubleArray(m_npoints);
 
        for (unsigned int i = 0; i < m_npoints; ++i)
        {
            instance.stage_y[q][i]    = ComplexSingleArray(m_nQuadPts, 0.0);
            instance.cstash_y[q][i]   = ComplexSingleArray(m_nQuadPts, 0.0);
            instance.csandbox_y[q][i] = ComplexSingleArray(m_nQuadPts, 0.0);
            instance.fstash_y[q][i]   = ComplexSingleArray(m_nQuadPts, 0.0);
            instance.fsandbox_y[q][i] = ComplexSingleArray(m_nQuadPts, 0.0);
        }
    }
 
    // Major storage for auxilliary ODE solutions.
    instance.stage_ind = std::pair<int, int>(0, 0); // Time-step counters
                                                    // indicating the interval
                                                    // ymain is associated with
 
    // Staging allocation
    instance.stage_active   = false;
    instance.stage_ccounter = 0;
    instance.stage_cbase    = pow(m_base, index - 1); // This base is halved
                                                      // after the first cycle
    instance.stage_fcounter = 0;
    instance.stage_fbase    = pow(m_base, index); // This base is halved
                                                  // after the first cycle
 
    // Ceiling stash allocation
    instance.cstash_counter = 0; // Counter used to determine
                                 // when to stash
 
    instance.cstash_base = pow(m_base, index - 1);   // base for counter ind(1)
    instance.cstash_ind = std::pair<int, int>(0, 0); // is never used: it always
                                                     // matches main.ind(1)
 
    // Ceiling sandbox allocation
    instance.csandbox_active = false; // Flag to determine when stash 2
                                      // is utilized
    instance.csandbox_counter = 0;
    instance.csandbox_ind     = std::pair<int, int>(0, 0);
 
    // Floor stash
    instance.fstash_base = 2 * pow(m_base, index);
    instance.fstash_ind  = std::pair<int, int>(0, 0);
 
    // Floor sandbox
    instance.fsandbox_active         = false;
    instance.fsandbox_activebase     = pow(m_base, index);
    instance.fsandbox_stashincrement = (m_base - 1) * pow(m_base, index - 1);
    instance.fsandbox_ind            = std::pair<int, int>(0, 0);
 
    // Defining parameters of the Talbot contour quadrature rule
    NekDouble Tl =
        m_deltaT * (2. * pow(m_base, index) - 1. - pow(m_base, index - 1));
    NekDouble mu = m_mu0 / Tl;
 
    // Talbot quadrature rule
    talbotQuadrature(m_nQuadPts, mu, m_nu, m_sigma, instance.z, instance.w);
 
    /**
     * /brief
     *
     * With sigma == 0, the dependence of z and w on index is just a
     * multiplicative scaling factor (mu). So technically we'll only
     * need one instance of this N-point rule and can scale it
     * accordingly inside each integral_class instance. Not sure if
     * this optimization is worth it. Cumulative memory savings would
     * only be about 4*N*Lmax floats.
 
     * Below: precomputation for time integration of auxiliary
     * variables.  Everything below here is independent of the
     * instance index index. Therefore, we could actually just
     * generate and store one copy of this stuff and use it
     * everywhere.
     */
 
    // 'As' array - one for each order.
    TripleArray &As = instance.As;
 
    As = TripleArray(m_order + 2);
 
    for (unsigned int m = 1; m <= m_order + 1; ++m)
    {
        As[m] = DoubleArray(m);
 
        for (unsigned int n = 0; n < m; ++n)
        {
            As[m][n] = SingleArray(m, 0.0);
        }
 
        switch (m)
        {
            case 1:
                As[m][0][0] = 1.;
                break;
 
            case 2:
                As[m][0][0] = 0.;
                As[m][0][1] = 1.;
                As[m][1][0] = 1.;
                As[m][1][1] = -1.;
                break;
 
            case 3:
                As[m][0][0] = 0.;
                As[m][0][1] = 1.;
                As[m][0][2] = 0;
                As[m][1][0] = 1. / 2.;
                As[m][1][1] = 0.;
                As[m][1][2] = -1. / 2.;
                As[m][2][0] = 1. / 2.;
                As[m][2][1] = -1.;
                As[m][2][2] = 1. / 2.;
                break;
 
            case 4:
                As[m][0][0] = 0.;
                As[m][0][1] = 1.;
                As[m][0][2] = 0.;
                As[m][0][3] = 0.;
 
                As[m][1][0] = 1. / 3.;
                As[m][1][1] = 1. / 2.;
                As[m][1][2] = -1.;
                As[m][1][3] = 1. / 6.;
 
                As[m][2][0] = 1. / 2.;
                As[m][2][1] = -1.;
                As[m][2][2] = 1. / 2.;
                As[m][2][3] = 0.;
 
                As[m][3][0] = 1. / 6.;
                As[m][3][1] = -1. / 2.;
                As[m][3][2] = 1. / 2.;
                As[m][3][3] = -1. / 6.;
                break;
 
            case 5:
                As[m][0][0] = 0.;
                As[m][0][1] = 1.;
                As[m][0][2] = 0.;
                As[m][0][3] = 0.;
                As[m][0][4] = 0.;
 
                As[m][1][0] = 1. / 4.;
                As[m][1][1] = 5. / 6.;
                As[m][1][2] = -3. / 2.;
                As[m][1][3] = 1. / 2.;
                As[m][1][4] = -1. / 12.;
 
                As[m][2][0] = 11. / 24.;
                As[m][2][1] = -5. / 6.;
                As[m][2][2] = 1. / 4.;
                As[m][2][3] = 1. / 6.;
                As[m][2][4] = -1. / 24.;
 
                As[m][3][0] = 1. / 4.;
                As[m][3][1] = -5. / 6.;
                As[m][3][2] = 1.;
                As[m][3][3] = -1. / 2.;
                As[m][3][4] = 1. / 12.;
 
                As[m][4][0] = 1. / 24.;
                As[m][4][1] = -1. / 6.;
                As[m][4][2] = 1. / 4.;
                As[m][4][3] = -1. / 6.;
                As[m][4][4] = 1. / 24.;
                break;
 
                // The default is a general formula, but the matrix inversion
                // involved is ill-conditioned, so the special cases above are
                // epxlicitly given to combat roundoff error in the most-used
                // scenarios.
            default:
                ASSERTL1(false, "No matrix inverse.");
 
                // Ainv = zeros(counter);
                // Ainv(1,:) = 1;
                // Ainv(2,1) = 1;
 
                // for( unsigned int j = 3; j<=counter; ++j )
                // {
                //          Ainv(j,:) = pow((j-2)., (0:counter));
                //          Ainv(j,:) = Ainv(j,:). * pow((-1)., (0:counter));
                // }
 
                // As[m] = inv(Ainv);
                break;
        }
    }
 
    // Initialize the exponenetial integrators.
    instance.E = ComplexSingleArray(m_nQuadPts, 0.0);
 
    for (unsigned int q = 0; q < m_nQuadPts; ++q)
    {
        instance.E[q] = exp(instance.z[q] * m_deltaT);
    }
 
    instance.Eh = ComplexDoubleArray(m_order + 1);
 
    for (unsigned int m = 0; m < m_order + 1; ++m)
    {
        instance.Eh[m] = ComplexSingleArray(m_nQuadPts, 0.0);
 
        for (unsigned int q = 0; q < m_nQuadPts; ++q)
        {
            if (m == 0)
                instance.Eh[0][q] =
                    1. / instance.z[q] * (exp(instance.z[q] * m_deltaT) - 1.0);
            else
                instance.Eh[m][q] = -1. / instance.z[q] +
                                    NekDouble(m) / (instance.z[q] * m_deltaT) *
                                        instance.Eh[m - 1][q];
        }
    }
 
    // 'AtEh' is set for the primary order. If a lower order method is
    // needed for initializing it will be changed in time_advance then
    // restored.
    instance.AtEh = ComplexDoubleArray(m_order + 1);
 
    for (unsigned int m = 0; m <= m_order; ++m)
    {
        instance.AtEh[m] = ComplexSingleArray(m_nQuadPts, 0.0);
 
        for (unsigned int q = 0; q < m_nQuadPts; ++q)
        {
            for (unsigned int i = 0; i <= m_order; ++i)
            {
                instance.AtEh[m][q] +=
                    instance.As[m_order + 1][m][i] * instance.Eh[i][q];
            }
        }
    }
}
 
/**
 * @brief Method to rearrange of staging/stashing for current time
 *
 * (1) activates ceiling staging
 * (2) moves ceiling stash ---> stage
 * (3) moves floor stash --> stage (+ updates all ceiling data)
 */
void FractionalInTimeIntegrationScheme::updateStage(const unsigned int timeStep,
                                                    Instance &instance)
{
    // Counter to flip active bit
    if (!instance.active)
    {
        instance.active = (timeStep % instance.activebase == 0);
    }
 
    // Determine if staging is necessary
    if (instance.active)
    {
        // Floor staging superscedes ceiling staging
        if (timeStep % instance.fstash_base == 0)
        {
            // Here a swap of the contents can be done because values
            // will copied into the stash and the f sandbox values will
            // cleared next.
            std::swap(instance.stage_y, instance.fstash_y);
            instance.stage_ind = instance.fstash_ind;
 
            std::swap(instance.csandbox_y, instance.fsandbox_y);
            instance.csandbox_ind = instance.fsandbox_ind;
 
            // After floor staging happens once, new base is base^index
            instance.fstash_base = pow(instance.base, instance.index);
 
            // Restart floor sandbox
            instance.fsandbox_ind    = std::pair<int, int>(0, 0);
            instance.fsandbox_active = false;
 
            // Clear the floor sandbox values.
            for (unsigned int i = 0; i < m_nvars; ++i)
            {
                for (unsigned int j = 0; j < m_npoints; ++j)
                {
                    for (unsigned int q = 0; q < m_nQuadPts; ++q)
                    {
                        instance.fsandbox_y[i][j][q] = 0;
                    }
                }
            }
        }
 
        // Check for ceiling staging
        else if (timeStep % instance.stage_cbase == 0)
        {
            instance.stage_ind = instance.cstash_ind;
 
            // A swap of the contents can be done because values will
            // copied into the stash.
            std::swap(instance.stage_y, instance.cstash_y);
        }
    }
}
 
/**
 * @brief Method to approximate the integral
 *
 *   \int_{(m-1) h}^{m h} k(m*h -s) f(u, s) dx{s}
 *
 * Using a time-stepping scheme of a particular order. Here, k depends
 * on alpha, the derivative order.
 */
void FractionalInTimeIntegrationScheme::finalIncrement(
    const unsigned int timeStep, const TimeIntegrationSchemeOperators &op)
{
    // Note: m_uNext is initialized to zero and then reset to zero
    // after it is used to update the current solution in TimeIntegrate.
    for (unsigned int m = 0; m < m_order; ++m)
    {
        op.DoOdeRhs(m_u[m], m_F, m_deltaT * (timeStep - m));
 
        for (unsigned int i = 0; i < m_nvars; ++i)
        {
            for (unsigned int j = 0; j < m_npoints; ++j)
            {
                m_uNext[i][j] += m_F[i][j] * m_AhattJ[m];
            }
        }
    }
}
 
/**
 * @brief Method to get the integral contribution over [taus(i+1)
 * taus(i)]. Stored in m_uInt.
 */
void FractionalInTimeIntegrationScheme::integralContribution(
    const unsigned int timeStep, const unsigned int tauml,
    const Instance &instance)
{
    // Assume y has already been updated to time level m
    for (unsigned int q = 0; q < m_nQuadPts; ++q)
    {
        m_expFactor[q] =
            exp(instance.z[q] * m_deltaT * NekDouble(timeStep - tauml)) *
            pow(instance.z[q], -m_alpha) * instance.w[q];
    }
 
    for (unsigned int i = 0; i < m_nvars; ++i)
    {
        for (unsigned int j = 0; j < m_npoints; ++j)
        {
            m_uInt[i][j] = 0;
 
            for (unsigned int q = 0; q < m_nQuadPts; ++q)
            {
                m_uInt[i][j] += instance.stage_y[i][j][q] * m_expFactor[q];
            }
 
            if (m_uInt[i][j].real() < 1e8)
            {
                m_uInt[i][j] = m_uInt[i][j].real();
            }
        }
    }
}
 
/**
 * @brief Method to get the solution to y' = z*y + f(u), using an
 * exponential integrator with implicit order (m_order + 1)
 * interpolation of the f(u) term.
 */
void FractionalInTimeIntegrationScheme::timeAdvance(
    const unsigned int timeStep, const TimeIntegrationSchemeOperators &op,
    Instance &instance, ComplexTripleArray &y)
{
    int order;
 
    // Try automated high-order method.
    if (timeStep <= m_order)
    {
        // Not enough history. For now, demote to lower-order method.
        // TODO: use multi-stage method.
        order = timeStep;
 
        // Prep for the time step.
        for (unsigned int m = 0; m <= order; ++m)
        {
            for (unsigned int q = 0; q < m_nQuadPts; ++q)
            {
                instance.AtEh[m][q] = 0;
 
                for (unsigned int i = 0; i <= order; ++i)
                {
                    instance.AtEh[m][q] +=
                        instance.As[order + 1][m][i] * instance.Eh[i][q];
                }
            }
        }
    }
    else
    {
        order = m_order;
    }
 
    // y = y * instance.E + F * instance.AtEh;
    for (unsigned int m = 0; m <= order; ++m)
    {
        op.DoOdeRhs(m_u[m], m_F, m_deltaT * (timeStep - m));
 
        for (unsigned int i = 0; i < m_nvars; ++i)
        {
            for (unsigned int j = 0; j < m_npoints; ++j)
            {
                for (unsigned int q = 0; q < m_nQuadPts; ++q)
                {
                    // y * instance.E
                    if (m == 0)
                        y[i][j][q] *= instance.E[q];
 
                    // F * instance.AtEh
                    y[i][j][q] += m_F[i][j] * instance.AtEh[m][q];
                }
            }
        }
    }
}
 
/**
 * @brief Method to update sandboxes to the current time.
 *
 * (1) advances ceiling sandbox
 * (2) moves ceiling sandbox ---> stash
 * (3) activates floor sandboxing
 * (4) advances floor sandbox
 * (5) moves floor sandbox ---> stash
 */
void FractionalInTimeIntegrationScheme::advanceSandbox(
    const unsigned int timeStep, const TimeIntegrationSchemeOperators &op,
    Instance &instance)
{
    // (1)
    // update(instance.csandbox.y)
    timeAdvance(timeStep, op, instance, instance.csandbox_y);
    instance.csandbox_ind.second = timeStep;
 
    // (2)
    // Determine if ceiling stashing is necessary
    instance.cstash_counter =
        modIncrement(instance.cstash_counter, instance.cstash_base);
 
    if (timeStep % instance.cstash_base == 0)
    {
        // Then need to stash
        // instance.cstash_y   = instance.csandbox_y;
        instance.cstash_ind = instance.csandbox_ind;
 
        // Stash the c sandbox value. This step has to be a deep copy
        // because the values in the sandbox are still needed for the
        // time advance.
        for (unsigned int i = 0; i < m_nvars; ++i)
        {
            for (unsigned int j = 0; j < m_npoints; ++j)
            {
                for (unsigned int q = 0; q < m_nQuadPts; ++q)
                {
                    instance.cstash_y[i][j][q] = instance.csandbox_y[i][j][q];
                }
            }
        }
    }
 
    if (instance.fsandbox_active)
    {
        // (4)
        timeAdvance(timeStep, op, instance, instance.fsandbox_y);
 
        instance.fsandbox_ind.second = timeStep;
 
        // (5) Move floor sandbox to stash
        if ((instance.fsandbox_ind.second - instance.fsandbox_ind.first) %
                instance.fsandbox_stashincrement ==
            0)
        {
            // instance.fstash_y   = instance.fsandbox_y;
            instance.fstash_ind = instance.fsandbox_ind;
 
            // Stash the f sandbox values. This step has to be a deep
            // copy because the values in the sandbox are still needed
            // for the time advance.
            for (unsigned int i = 0; i < m_nvars; ++i)
            {
                for (unsigned int j = 0; j < m_npoints; ++j)
                {
                    for (unsigned int q = 0; q < m_nQuadPts; ++q)
                    {
                        instance.fstash_y[i][j][q] =
                            instance.fsandbox_y[i][j][q];
                    }
                }
            }
        }
    }
    else // Determine if advancing floor sandbox is necessary at next time
    {
        // (3)
        if (timeStep % instance.fsandbox_activebase == 0)
        {
            instance.fsandbox_active = true;
            instance.fsandbox_ind    = std::pair<int, int>(timeStep, timeStep);
        }
    }
}
 
/**
 * @brief Worker method to print details on the integration scheme
 */
void FractionalInTimeIntegrationScheme::v_print(std::ostream &os) const
{
    os << "Time Integration Scheme: " << GetFullName() << std::endl
       << "       Alpha " << m_alpha << std::endl
       << "       Base " << m_base << std::endl
       << "       Number of instances " << m_Lmax << std::endl
       << "       Number of quadature points " << m_nQuadPts << std::endl
       << "             Talbot Parameter: sigma " << m_sigma << std::endl
       << "             Talbot Parameter: mu0 " << m_mu0 << std::endl
       << "             Talbot Parameter: nu " << m_nu << std::endl;
}
 
void FractionalInTimeIntegrationScheme::v_printFull(std::ostream &os) const
{
    os << "Time Integration Scheme: " << GetFullName() << std::endl
       << "       Alpha " << m_alpha << std::endl
       << "       Base " << m_base << std::endl
       << "       Number of instances " << m_Lmax << std::endl
       << "       Number of quadature points " << m_nQuadPts << std::endl
       << "             Talbot Parameter: sigma " << m_sigma << std::endl
       << "             Talbot Parameter: mu0 " << m_mu0 << std::endl
       << "             Talbot Parameter: nu " << m_nu << std::endl;
}
 
// Friend Operators
std::ostream &operator<<(std::ostream &os,
                         const FractionalInTimeIntegrationScheme &rhs)
{
    rhs.print(os);
 
    return os;
}
 
std::ostream &operator<<(std::ostream &os,
                         const FractionalInTimeIntegrationSchemeSharedPtr &rhs)
{
    os << *rhs.get();
 
    return os;
}
 
} // end namespace LibUtilities
} // namespace Nektar