Class for fractional-in-time integration. More...

#include <TimeIntegrationSchemeFIT.h>

Inheritance diagram for Nektar::LibUtilities::FractionalInTimeIntegrationScheme:

Classes
struct	Instance

Public Member Functions
	FractionalInTimeIntegrationScheme (std::string variant, size_t order, std::vector< NekDouble > freeParams)
	Constructor. More...

	~FractionalInTimeIntegrationScheme () override
	Destructor. More...

Public Member Functions inherited from Nektar::LibUtilities::TimeIntegrationScheme
LUE std::string	GetFullName () const

LUE std::string	GetName () const

LUE std::string	GetVariant () const

LUE size_t	GetOrder () const

LUE std::vector< NekDouble >	GetFreeParams ()

LUE TimeIntegrationSchemeType	GetIntegrationSchemeType ()

LUE NekDouble	GetTimeStability () const

LUE size_t	GetNumIntegrationPhases ()

LUE const TripleArray &	GetSolutionVector () const
	Gets the solution vector of the ODE. More...

LUE TripleArray &	UpdateSolutionVector ()

LUE void	SetSolutionVector (const size_t Offset, const DoubleArray &y)
	Sets the solution vector of the ODE. More...

LUE void	InitializeScheme (const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)
	Explicit integration of an ODE. More...

LUE ConstDoubleArray &	TimeIntegrate (const size_t timestep, const NekDouble delta_t)

LUE void	print (std::ostream &os) const

LUE void	printFull (std::ostream &os) const

Static Public Member Functions
static TimeIntegrationSchemeSharedPtr	create (std::string variant, size_t order, std::vector< NekDouble > freeParams)
	Creator. More...

Static Public Attributes
static std::string	className

Protected Member Functions
LUE std::string	v_GetName () const override

LUE std::string	v_GetVariant () const override

LUE size_t	v_GetOrder () const override

LUE std::vector< NekDouble >	v_GetFreeParams () const override

LUE TimeIntegrationSchemeType	v_GetIntegrationSchemeType () const override

LUE NekDouble	v_GetTimeStability () const override

LUE size_t	v_GetNumIntegrationPhases () const override

const TripleArray &	v_GetSolutionVector () const override
	Gets the solution vector of the ODE. More...

TripleArray &	v_UpdateSolutionVector () override

void	v_SetSolutionVector (const size_t Offset, const DoubleArray &y) override
	Sets the solution vector of the ODE. More...

LUE void	v_InitializeScheme (const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op) override
	Worker method to initialize the integration scheme. More...

LUE ConstDoubleArray &	v_TimeIntegrate (const size_t timestep, const NekDouble delta_t) override
	Worker method that performs the time integration. More...

LUE void	v_print (std::ostream &os) const override
	Worker method to print details on the integration scheme. More...

LUE void	v_printFull (std::ostream &os) const override

size_t	modIncrement (const size_t counter, const size_t base) const
	Method that increments the counter then performs mod calculation. More...

size_t	computeL (const size_t base, const size_t m) const
	Method to compute the smallest integer L such that base < 2. More...

size_t	computeQML (const size_t base, const size_t m)
	Method to compute the demarcation integers q_{m, ell}. More...

size_t	computeTaus (const size_t base, const size_t m)
	Method to compute the demarcation interval marker tau_{m, ell}. More...

void	talbotQuadrature (const size_t nQuadPts, const NekDouble mu, const NekDouble nu, const NekDouble sigma, ComplexSingleArray &lamb, ComplexSingleArray &w) const
	Method to compute the quadrature rule over Tablot contour. More...

void	integralClassInitialize (const size_t index, Instance &instance) const
	Method to initialize the integral class. More...

void	updateStage (const size_t timeStep, Instance &instance)
	Method to rearrange of staging/stashing for current time. More...

void	finalIncrement (const size_t timeStep)
	Method to approximate the integral. More...

void	integralContribution (const size_t timeStep, const size_t tauml, const Instance &instance)
	Method to get the integral contribution over [taus(i+1) taus(i)]. Stored in m_uInt. More...

void	timeAdvance (const size_t timeStep, Instance &instance, ComplexTripleArray &y)
	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. More...

void	advanceSandbox (const size_t timeStep, Instance &instance)
	Method to update sandboxes to the current time. More...

Protected Member Functions inherited from Nektar::LibUtilities::TimeIntegrationScheme
virtual LUE std::string	v_GetFullName () const

virtual LUE std::string	v_GetName () const =0

virtual LUE std::string	v_GetVariant () const =0

virtual LUE size_t	v_GetOrder () const =0

virtual LUE std::vector< NekDouble >	v_GetFreeParams () const =0

virtual LUE TimeIntegrationSchemeType	v_GetIntegrationSchemeType () const =0

virtual LUE NekDouble	v_GetTimeStability () const =0

virtual LUE size_t	v_GetNumIntegrationPhases () const =0

virtual LUE const TripleArray &	v_GetSolutionVector () const =0

virtual LUE TripleArray &	v_UpdateSolutionVector ()=0

virtual LUE void	v_SetSolutionVector (const size_t Offset, const DoubleArray &y)=0

virtual LUE void	v_InitializeScheme (const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)=0

virtual LUE ConstDoubleArray &	v_TimeIntegrate (const size_t timestep, const NekDouble delta_t)=0

virtual LUE void	v_print (std::ostream &os) const =0

virtual LUE void	v_printFull (std::ostream &os) const =0

LUE	TimeIntegrationScheme (std::string variant, size_t order, std::vector< NekDouble > freeParams)

LUE	TimeIntegrationScheme (const TimeIntegrationScheme &in)=delete

virtual	~TimeIntegrationScheme ()=default

Protected Attributes
std::string	m_name

std::string	m_variant

size_t	m_order {0}

std::vector< NekDouble >	m_freeParams

TimeIntegrationSchemeType	m_schemeType {eFractionalInTime}

TimeIntegrationSchemeOperators	m_op

NekDouble	m_deltaT {0}

NekDouble	m_T {0}

size_t	m_maxTimeSteps

NekDouble	m_alpha {0.3}

size_t	m_base {4}

size_t	m_nQuadPts {20}

NekDouble	m_sigma {0}

NekDouble	m_mu0 {8}

NekDouble	m_nu {0.6}

size_t	m_nvars {0}

size_t	m_npoints {0}

size_t	m_Lmax {0}

Array< OneD, Instance >	m_integral_classes

Array< OneD, size_t >	m_qml

Array< OneD, size_t >	m_taus

DoubleArray	m_u0

DoubleArray	m_uNext

ComplexDoubleArray	m_uInt

ComplexSingleArray	m_expFactor

TripleArray	m_u

DoubleArray	m_F

SingleArray	m_J

TripleArray	m_Ahats

SingleArray	m_AhattJ

Friends
LUE friend std::ostream &	operator<< (std::ostream &os, const FractionalInTimeIntegrationScheme &rhs)

LUE friend std::ostream &	operator<< (std::ostream &os, const FractionalInTimeIntegrationSchemeSharedPtr &rhs)

Detailed Description

Class for fractional-in-time integration.

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).

Definition at line 54 of file TimeIntegrationSchemeFIT.h.

Constructor & Destructor Documentation

◆ FractionalInTimeIntegrationScheme()

Nektar::LibUtilities::FractionalInTimeIntegrationScheme::FractionalInTimeIntegrationScheme	(	std::string	variant,
		size_t	order,
		std::vector< NekDouble >	freeParams
	)

Constructor.

Definition at line 53 of file TimeIntegrationSchemeFIT.cpp.

    : TimeIntegrationScheme(variant, order, freeParams),
      m_name("FractionalInTime")
{
    m_variant    = variant;
    m_order      = order;
    m_freeParams = freeParams;
 
    // Currently up to 4th order is implemented.
    ASSERTL1(1 <= 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];
    }
}

References ASSERTL1, m_alpha, m_base, m_freeParams, m_mu0, m_nQuadPts, m_nu, m_order, m_sigma, and m_variant.

◆ ~FractionalInTimeIntegrationScheme()

Nektar::LibUtilities::FractionalInTimeIntegrationScheme::~FractionalInTimeIntegrationScheme ( )

inlineoverride

Destructor.

Definition at line 62 of file TimeIntegrationSchemeFIT.h.

63 {

64 }

Member Function Documentation

◆ advanceSandbox()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::advanceSandbox	(	const size_t	timeStep,
		Instance &	instance
	)

protected

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

Definition at line 965 of file TimeIntegrationSchemeFIT.cpp.

{
    // (1)
    // update(instance.csandbox.y)
    timeAdvance(timeStep, 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 (size_t i = 0; i < m_nvars; ++i)
        {
            for (size_t j = 0; j < m_npoints; ++j)
            {
                for (size_t 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, 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 (size_t i = 0; i < m_nvars; ++i)
            {
                for (size_t j = 0; j < m_npoints; ++j)
                {
                    for (size_t 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<size_t, size_t>(timeStep, timeStep);
        }
    }
}

Referenced by v_TimeIntegrate().

◆ computeL()

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::computeL	(	const size_t	base,
		const size_t	l
	)		const

inlineprotected

Method to compute the smallest integer L such that base < 2.

base^l.

Definition at line 368 of file TimeIntegrationSchemeFIT.cpp.

{
    size_t L = ceil(log(l / 2.0) / log(base));
 
    if (l % (size_t)(2 * pow(base, L)) == 0)
    {
        ++L;
    }
 
    return L;
}

References L, and tinysimd::log().

Referenced by computeQML(), and v_InitializeScheme().

◆ computeQML()

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::computeQML	(	const size_t	base,
		const size_t	m
	)

inlineprotected

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.

Definition at line 388 of file TimeIntegrationSchemeFIT.cpp.

{
    size_t L = computeL(base, m);
 
    // m_qml is set in InitializeScheme to be the largest length expected.
    // qml = Array<OneD, size_t>( L-1, 0 );
 
    for (size_t i = 0; i < L - 1; ++i)
    {
        m_qml[i] = floor(m / pow(base, i + 1)) - 1;
    }
 
    return L;
}

References computeL(), L, and m_qml.

Referenced by computeTaus().

◆ computeTaus()

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::computeTaus	(	const size_t	base,
		const size_t	m
	)

inlineprotected

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.

Definition at line 411 of file TimeIntegrationSchemeFIT.cpp.

{
    if (m == 1)
    {
        m_taus[0] = 0;
 
        return 0;
    }
    else
    {
        size_t L = computeQML(base, m);
 
        // m_taus is set in InitializeScheme to be the largest length
        // expected.
 
        m_taus[0] = m - 1;
 
        for (size_t i = 1; i < L; ++i)
        {
            m_taus[i] = m_qml[i - 1] * pow(base, i);
        }
 
        m_taus[L] = 0;
 
        return L;
    }
}

References computeQML(), L, m_qml, and m_taus.

Referenced by v_TimeIntegrate().

◆ create()

static TimeIntegrationSchemeSharedPtr Nektar::LibUtilities::FractionalInTimeIntegrationScheme::create	(	std::string	variant,
		size_t	order,
		std::vector< NekDouble >	freeParams
	)

inlinestatic

Creator.

Definition at line 67 of file TimeIntegrationSchemeFIT.h.

    {
        TimeIntegrationSchemeSharedPtr p =
            MemoryManager<FractionalInTimeIntegrationScheme>::AllocateSharedPtr(
                variant, order, freeParams);
 
        return p;
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ finalIncrement()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::finalIncrement ( const size_t timeStep )

protected

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.

Definition at line 841 of file TimeIntegrationSchemeFIT.cpp.

{
    // Note: m_uNext is initialized to zero and then reset to zero
    // after it is used to update the current solution in TimeIntegrate.
    for (size_t m = 0; m < m_order; ++m)
    {
        m_op.DoOdeRhs(m_u[m], m_F, m_deltaT * (timeStep - m));
 
        for (size_t i = 0; i < m_nvars; ++i)
        {
            for (size_t j = 0; j < m_npoints; ++j)
            {
                m_uNext[i][j] += m_F[i][j] * m_AhattJ[m];
            }
        }
    }
}

References Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs(), m_AhattJ, m_deltaT, m_F, m_npoints, m_nvars, m_op, m_order, m_u, and m_uNext.

Referenced by v_TimeIntegrate().

◆ integralClassInitialize()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::integralClassInitialize	(	const size_t	index,
		Instance &	instance
	)		const

protected

Method to initialize the integral class.

/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.

/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.

Definition at line 486 of file TimeIntegrationSchemeFIT.cpp.

{
    /**
     * /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 (size_t 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 (size_t 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<size_t, size_t>(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<size_t, size_t>(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<size_t, size_t>(0, 0);
 
    // Floor stash
    instance.fstash_base = 2 * pow(m_base, index);
    instance.fstash_ind  = std::pair<size_t, size_t>(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<size_t, size_t>(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 (size_t m = 1; m <= m_order + 1; ++m)
    {
        As[m] = DoubleArray(m);
 
        for (size_t 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.");
                break;
        }
    }
 
    // Initialize the exponenetial integrators.
    instance.E = ComplexSingleArray(m_nQuadPts, 0.0);
 
    for (size_t q = 0; q < m_nQuadPts; ++q)
    {
        instance.E[q] = exp(instance.z[q] * m_deltaT);
    }
 
    instance.Eh = ComplexDoubleArray(m_order + 1);
 
    for (size_t m = 0; m < m_order + 1; ++m)
    {
        instance.Eh[m] = ComplexSingleArray(m_nQuadPts, 0.0);
 
        for (size_t 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 (size_t m = 0; m <= m_order; ++m)
    {
        instance.AtEh[m] = ComplexSingleArray(m_nQuadPts, 0.0);
 
        for (size_t q = 0; q < m_nQuadPts; ++q)
        {
            for (size_t i = 0; i <= m_order; ++i)
            {
                instance.AtEh[m][q] +=
                    instance.As[m_order + 1][m][i] * instance.Eh[i][q];
            }
        }
    }
}

Referenced by v_InitializeScheme().

◆ integralContribution()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::integralContribution	(	const size_t	timeStep,
		const size_t	tauml,
		const Instance &	instance
	)

protected

Method to get the integral contribution over [taus(i+1) taus(i)]. Stored in m_uInt.

Definition at line 863 of file TimeIntegrationSchemeFIT.cpp.

{
    // Assume y has already been updated to time level m
    for (size_t 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 (size_t i = 0; i < m_nvars; ++i)
    {
        for (size_t j = 0; j < m_npoints; ++j)
        {
            m_uInt[i][j] = 0;
 
            for (size_t 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();
            }
        }
    }
}

References m_alpha, m_deltaT, m_expFactor, m_npoints, m_nQuadPts, m_nvars, m_uInt, Nektar::UnitTests::q(), Nektar::LibUtilities::FractionalInTimeIntegrationScheme::Instance::stage_y, Nektar::LibUtilities::FractionalInTimeIntegrationScheme::Instance::w, and Nektar::LibUtilities::FractionalInTimeIntegrationScheme::Instance::z.

Referenced by v_TimeIntegrate().

◆ modIncrement()

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::modIncrement	(	const size_t	counter,
		const size_t	base
	)		const

inlineprotected

Method that increments the counter then performs mod calculation.

Definition at line 358 of file TimeIntegrationSchemeFIT.cpp.

{
    return (counter + 1) % base;
}

Referenced by advanceSandbox().

◆ talbotQuadrature()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::talbotQuadrature	(	const size_t	nQuadPts,
		const NekDouble	mu,
		const NekDouble	nu,
		const NekDouble	sigma,
		ComplexSingleArray &	lamb,
		ComplexSingleArray &	w
	)		const

protected

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.

Definition at line 451 of file TimeIntegrationSchemeFIT.cpp.

{
    lamb = ComplexSingleArray(nQuadPts, 0.0);
    w    = ComplexSingleArray(nQuadPts, 0.0);
 
    for (size_t 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)
    {
        size_t q = (nQuadPts + 1) / 2;
 
        lamb[q] = std::complex<NekDouble>(sigma + mu, 0);
 
        w[q] = std::complex<NekDouble>(nu * mu / nQuadPts, 0);
    }
}

References Nektar::UnitTests::q(), and Nektar::UnitTests::w().

Referenced by integralClassInitialize().

◆ timeAdvance()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::timeAdvance	(	const size_t	timeStep,
		Instance &	instance,
		ComplexTripleArray &	y
	)

protected

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.

Definition at line 898 of file TimeIntegrationSchemeFIT.cpp.

{
    size_t 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 (size_t m = 0; m <= order; ++m)
        {
            for (size_t q = 0; q < m_nQuadPts; ++q)
            {
                instance.AtEh[m][q] = 0;
 
                for (size_t 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 (size_t m = 0; m <= order; ++m)
    {
        m_op.DoOdeRhs(m_u[m], m_F, m_deltaT * (timeStep - m));
 
        for (size_t i = 0; i < m_nvars; ++i)
        {
            for (size_t j = 0; j < m_npoints; ++j)
            {
                for (size_t 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];
                }
            }
        }
    }
}

References Nektar::LibUtilities::FractionalInTimeIntegrationScheme::Instance::As, Nektar::LibUtilities::FractionalInTimeIntegrationScheme::Instance::AtEh, Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs(), Nektar::LibUtilities::FractionalInTimeIntegrationScheme::Instance::E, Nektar::LibUtilities::FractionalInTimeIntegrationScheme::Instance::Eh, m_deltaT, m_F, m_npoints, m_nQuadPts, m_nvars, m_op, m_order, m_u, and Nektar::UnitTests::q().

Referenced by advanceSandbox().

◆ updateStage()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::updateStage	(	const size_t	timeStep,
		Instance &	instance
	)

protected

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)

Definition at line 777 of file TimeIntegrationSchemeFIT.cpp.

{
    // 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<size_t, size_t>(0, 0);
            instance.fsandbox_active = false;
 
            // Clear the floor sandbox values.
            for (size_t i = 0; i < m_nvars; ++i)
            {
                for (size_t j = 0; j < m_npoints; ++j)
                {
                    for (size_t 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);
        }
    }
}

Referenced by v_TimeIntegrate().

◆ v_GetFreeParams()

LUE std::vector< NekDouble > Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_GetFreeParams ( ) const

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 103 of file TimeIntegrationSchemeFIT.h.

    {
        return m_freeParams;
    }

References m_freeParams.

◆ v_GetIntegrationSchemeType()

LUE TimeIntegrationSchemeType Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_GetIntegrationSchemeType ( ) const

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 108 of file TimeIntegrationSchemeFIT.h.

    {
        return m_schemeType;
    }

References m_schemeType.

◆ v_GetName()

LUE std::string Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_GetName ( ) const

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 88 of file TimeIntegrationSchemeFIT.h.

    {
        return m_name;
    }

References m_name.

◆ v_GetNumIntegrationPhases()

LUE size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_GetNumIntegrationPhases ( ) const

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 118 of file TimeIntegrationSchemeFIT.h.

    {
        return 1;
    }

◆ v_GetOrder()

LUE size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_GetOrder ( ) const

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 98 of file TimeIntegrationSchemeFIT.h.

    {
        return m_order;
    }

References m_order.

◆ v_GetSolutionVector()

const TripleArray & Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_GetSolutionVector ( ) const

inlineoverrideprotectedvirtual

Gets the solution vector of the ODE.

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 126 of file TimeIntegrationSchemeFIT.h.

    {
        return m_u;
    }

References m_u.

◆ v_GetTimeStability()

LUE NekDouble Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_GetTimeStability ( ) const

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 113 of file TimeIntegrationSchemeFIT.h.

    {
        return 1.0;
    }

◆ v_GetVariant()

LUE std::string Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_GetVariant ( ) const

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 93 of file TimeIntegrationSchemeFIT.h.

    {
        return m_variant;
    }

References m_variant.

◆ v_InitializeScheme()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_InitializeScheme	(	const NekDouble	deltaT,
		ConstDoubleArray &	y_0,
		const NekDouble	time,
		const TimeIntegrationSchemeOperators &	op
	)

overrideprotectedvirtual

Worker method to initialize the integration scheme.

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 99 of file TimeIntegrationSchemeFIT.cpp.

{
    m_op      = 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, size_t>(m_Lmax - 1, size_t(0));
 
    // Demarcation interval markers - one array that is re-used
    m_taus = Array<OneD, size_t>(m_Lmax + 1, size_t(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 (size_t m = 0; m <= m_order; ++m)
    {
        m_u[m] = DoubleArray(m_nvars);
 
        for (size_t i = 0; i < m_nvars; ++i)
        {
            m_u[m][i] = SingleArray(m_npoints, 0.0);
 
            for (size_t 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 (size_t 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) / std::tgamma(m_alpha + 1.);
 
    for (size_t 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 (size_t m = 1; m <= m_order; ++m)
    {
        m_Ahats[m] = DoubleArray(m);
 
        for (size_t 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 (size_t j = 2; j <= m; ++j)
                {
                    for (size_t 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 (size_t i = 0; i < m_order; ++i)
    {
        for (size_t 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 (size_t l = 0; l < m_Lmax; ++l)
    {
        integralClassInitialize(l + 1, m_integral_classes[l]);
    }
}

References ASSERTL1, computeL(), integralClassInitialize(), m_Ahats, m_AhattJ, m_alpha, m_base, m_deltaT, m_expFactor, m_F, m_integral_classes, m_J, m_Lmax, m_maxTimeSteps, m_npoints, m_nQuadPts, m_nvars, m_op, m_order, m_qml, m_T, m_taus, m_u, m_u0, m_uInt, and m_uNext.

◆ v_print()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_print ( std::ostream & os ) const

overrideprotectedvirtual

Worker method to print details on the integration scheme.

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 1045 of file TimeIntegrationSchemeFIT.cpp.

{
    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;
}

References Nektar::LibUtilities::TimeIntegrationScheme::GetFullName(), m_alpha, m_base, m_Lmax, m_mu0, m_nQuadPts, m_nu, and m_sigma.

◆ v_printFull()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_printFull ( std::ostream & os ) const

overrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 1057 of file TimeIntegrationSchemeFIT.cpp.

{
    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;
}

References Nektar::LibUtilities::TimeIntegrationScheme::GetFullName(), m_alpha, m_base, m_Lmax, m_mu0, m_nQuadPts, m_nu, and m_sigma.

◆ v_SetSolutionVector()

void Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_SetSolutionVector	(	const size_t	Offset,
		const DoubleArray &	y
	)

inlineoverrideprotectedvirtual

Sets the solution vector of the ODE.

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 138 of file TimeIntegrationSchemeFIT.h.

    {
        m_u[Offset] = y;
    }

References m_u.

◆ v_TimeIntegrate()

ConstDoubleArray & Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_TimeIntegrate	(	const size_t	timestep,
		const NekDouble	delta_t
	)

overrideprotectedvirtual

Worker method that performs the time integration.

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 283 of file TimeIntegrationSchemeFIT.cpp.

{
    ASSERTL1(delta_t == m_deltaT,
             "Delta T has changed which is not permitted.");
 
    // The Fractional in Time works via the logical? time step value.
    size_t timeStep = timestep + 1;
 
    // Update the storage and counters for integral classes.  Performs
    // staging for updating u.
    for (size_t 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);
 
    // Contributions to the current integral
    size_t L = computeTaus(m_base, timeStep);
 
    for (size_t 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 (size_t i = 0; i < m_nvars; ++i)
        {
            for (size_t 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 (size_t m = m_order; m > 0; --m)
    {
        for (size_t i = 0; i < m_nvars; ++i)
        {
            for (size_t j = 0; j < m_npoints; ++j)
            {
                m_u[m][i][j] = m_u[m - 1][i][j];
            }
        }
    }
 
    // Get the current solution.
    for (size_t i = 0; i < m_nvars; ++i)
    {
        for (size_t 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.
        }
    }
 
    // Update the storage and counters for integral classes to
    // time timeStep * m_deltaT. Also time-steps the sandboxes and stashes.
    for (size_t i = 0; i < m_Lmax; ++i)
    {
        advanceSandbox(timeStep, m_integral_classes[i]);
    }
 
    return m_u[0];
}

References advanceSandbox(), ASSERTL1, computeTaus(), finalIncrement(), integralContribution(), L, m_base, m_deltaT, m_integral_classes, m_Lmax, m_npoints, m_nvars, m_order, m_taus, m_u, m_u0, m_uInt, m_uNext, and updateStage().

◆ v_UpdateSolutionVector()

TripleArray & Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_UpdateSolutionVector ( )

inlineoverrideprotectedvirtual

Implements Nektar::LibUtilities::TimeIntegrationScheme.

Definition at line 130 of file TimeIntegrationSchemeFIT.h.

    {
        return m_u;
    }

References m_u.

Friends And Related Function Documentation

◆ operator<< [1/2]

LUE friend std::ostream & operator<<	(	std::ostream &	os,
		const FractionalInTimeIntegrationScheme &	rhs
	)

friend

Definition at line 1070 of file TimeIntegrationSchemeFIT.cpp.

{
    rhs.print(os);
 
    return os;
}

◆ operator<< [2/2]

LUE friend std::ostream & operator<<	(	std::ostream &	os,
		const FractionalInTimeIntegrationSchemeSharedPtr &	rhs
	)

friend

Definition at line 1078 of file TimeIntegrationSchemeFIT.cpp.

{
    os << *rhs.get();
 
    return os;
}

Member Data Documentation

◆ className

std::string Nektar::LibUtilities::FractionalInTimeIntegrationScheme::className

static

Initial value:

=
    GetTimeIntegrationSchemeFactory().RegisterCreatorFunction(
        "FractionalInTime", FractionalInTimeIntegrationScheme::create)

Definition at line 77 of file TimeIntegrationSchemeFIT.h.

◆ m_Ahats

TripleArray Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_Ahats

protected

Definition at line 293 of file TimeIntegrationSchemeFIT.h.

Referenced by v_InitializeScheme().

◆ m_AhattJ

SingleArray Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_AhattJ

protected

Definition at line 296 of file TimeIntegrationSchemeFIT.h.

Referenced by finalIncrement(), and v_InitializeScheme().

◆ m_alpha

NekDouble Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_alpha {0.3}

protected

Definition at line 255 of file TimeIntegrationSchemeFIT.h.

Referenced by FractionalInTimeIntegrationScheme(), integralContribution(), v_InitializeScheme(), v_print(), and v_printFull().

◆ m_base

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_base {4}

protected

Definition at line 256 of file TimeIntegrationSchemeFIT.h.

Referenced by FractionalInTimeIntegrationScheme(), integralClassInitialize(), v_InitializeScheme(), v_print(), v_printFull(), and v_TimeIntegrate().

◆ m_deltaT

NekDouble Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_deltaT {0}

protected

Definition at line 252 of file TimeIntegrationSchemeFIT.h.

Referenced by finalIncrement(), integralClassInitialize(), integralContribution(), timeAdvance(), v_InitializeScheme(), and v_TimeIntegrate().

◆ m_expFactor

ComplexSingleArray Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_expFactor

protected

Definition at line 280 of file TimeIntegrationSchemeFIT.h.

Referenced by integralContribution(), and v_InitializeScheme().

◆ m_F

DoubleArray Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_F

protected

Definition at line 287 of file TimeIntegrationSchemeFIT.h.

Referenced by finalIncrement(), timeAdvance(), and v_InitializeScheme().

◆ m_freeParams

std::vector<NekDouble> Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_freeParams

protected

Definition at line 246 of file TimeIntegrationSchemeFIT.h.

Referenced by FractionalInTimeIntegrationScheme(), and v_GetFreeParams().

◆ m_integral_classes

Array<OneD, Instance> Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_integral_classes

protected

Definition at line 266 of file TimeIntegrationSchemeFIT.h.

Referenced by v_InitializeScheme(), and v_TimeIntegrate().

◆ m_J

SingleArray Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_J

protected

Definition at line 290 of file TimeIntegrationSchemeFIT.h.

Referenced by v_InitializeScheme().

◆ m_Lmax

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_Lmax {0}

protected

Definition at line 265 of file TimeIntegrationSchemeFIT.h.

Referenced by v_InitializeScheme(), v_print(), v_printFull(), and v_TimeIntegrate().

◆ m_maxTimeSteps

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_maxTimeSteps

protected

Definition at line 254 of file TimeIntegrationSchemeFIT.h.

Referenced by v_InitializeScheme().

◆ m_mu0

NekDouble Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_mu0 {8}

protected

Definition at line 259 of file TimeIntegrationSchemeFIT.h.

Referenced by FractionalInTimeIntegrationScheme(), integralClassInitialize(), v_print(), and v_printFull().

◆ m_name

std::string Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_name

protected

Definition at line 243 of file TimeIntegrationSchemeFIT.h.

Referenced by v_GetName().

◆ m_npoints

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_npoints {0}

protected

Definition at line 263 of file TimeIntegrationSchemeFIT.h.

Referenced by advanceSandbox(), finalIncrement(), integralClassInitialize(), integralContribution(), timeAdvance(), updateStage(), v_InitializeScheme(), and v_TimeIntegrate().

◆ m_nQuadPts

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_nQuadPts {20}

protected

Definition at line 257 of file TimeIntegrationSchemeFIT.h.

Referenced by advanceSandbox(), FractionalInTimeIntegrationScheme(), integralClassInitialize(), integralContribution(), timeAdvance(), updateStage(), v_InitializeScheme(), v_print(), and v_printFull().

◆ m_nu

NekDouble Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_nu {0.6}

protected

Definition at line 260 of file TimeIntegrationSchemeFIT.h.

Referenced by FractionalInTimeIntegrationScheme(), integralClassInitialize(), v_print(), and v_printFull().

◆ m_nvars

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_nvars {0}

protected

Definition at line 262 of file TimeIntegrationSchemeFIT.h.

Referenced by advanceSandbox(), finalIncrement(), integralClassInitialize(), integralContribution(), timeAdvance(), updateStage(), v_InitializeScheme(), and v_TimeIntegrate().

◆ m_op

TimeIntegrationSchemeOperators Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_op

protected

Definition at line 249 of file TimeIntegrationSchemeFIT.h.

Referenced by finalIncrement(), timeAdvance(), and v_InitializeScheme().

◆ m_order

size_t Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_order {0}

protected

Definition at line 245 of file TimeIntegrationSchemeFIT.h.

Referenced by finalIncrement(), FractionalInTimeIntegrationScheme(), integralClassInitialize(), timeAdvance(), v_GetOrder(), v_InitializeScheme(), and v_TimeIntegrate().

◆ m_qml

Array<OneD, size_t> Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_qml

protected

Definition at line 269 of file TimeIntegrationSchemeFIT.h.

Referenced by computeQML(), computeTaus(), and v_InitializeScheme().

◆ m_schemeType

TimeIntegrationSchemeType Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_schemeType {eFractionalInTime}

protected

Definition at line 248 of file TimeIntegrationSchemeFIT.h.

Referenced by v_GetIntegrationSchemeType().

◆ m_sigma

NekDouble Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_sigma {0}

protected

Definition at line 258 of file TimeIntegrationSchemeFIT.h.

Referenced by FractionalInTimeIntegrationScheme(), integralClassInitialize(), v_print(), and v_printFull().

◆ m_T

NekDouble Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_T {0}

protected

Definition at line 253 of file TimeIntegrationSchemeFIT.h.

Referenced by v_InitializeScheme().

◆ m_taus

Array<OneD, size_t> Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_taus

protected

Definition at line 271 of file TimeIntegrationSchemeFIT.h.

Referenced by computeTaus(), v_InitializeScheme(), and v_TimeIntegrate().

◆ m_u

TripleArray Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_u

protected

Definition at line 283 of file TimeIntegrationSchemeFIT.h.

Referenced by finalIncrement(), timeAdvance(), v_GetSolutionVector(), v_InitializeScheme(), v_SetSolutionVector(), v_TimeIntegrate(), and v_UpdateSolutionVector().

◆ m_u0

DoubleArray Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_u0

protected

Definition at line 274 of file TimeIntegrationSchemeFIT.h.

Referenced by v_InitializeScheme(), and v_TimeIntegrate().

◆ m_uInt

ComplexDoubleArray Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_uInt

protected

Definition at line 278 of file TimeIntegrationSchemeFIT.h.

Referenced by integralContribution(), v_InitializeScheme(), and v_TimeIntegrate().

◆ m_uNext

DoubleArray Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_uNext

protected

Definition at line 276 of file TimeIntegrationSchemeFIT.h.

Referenced by finalIncrement(), v_InitializeScheme(), and v_TimeIntegrate().

◆ m_variant

std::string Nektar::LibUtilities::FractionalInTimeIntegrationScheme::m_variant

protected

Definition at line 244 of file TimeIntegrationSchemeFIT.h.

Referenced by FractionalInTimeIntegrationScheme(), and v_GetVariant().

Classes

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

◆ FractionalInTimeIntegrationScheme()

◆ ~FractionalInTimeIntegrationScheme()

Member Function Documentation

◆ advanceSandbox()

◆ computeL()

◆ computeQML()

◆ computeTaus()

◆ create()

◆ finalIncrement()

◆ integralClassInitialize()

◆ integralContribution()

◆ modIncrement()

◆ talbotQuadrature()

◆ timeAdvance()

◆ updateStage()

◆ v_GetFreeParams()

◆ v_GetIntegrationSchemeType()

◆ v_GetName()

◆ v_GetNumIntegrationPhases()

◆ v_GetOrder()

◆ v_GetSolutionVector()

◆ v_GetTimeStability()

◆ v_GetVariant()

◆ v_InitializeScheme()

◆ v_print()

◆ v_printFull()

◆ v_SetSolutionVector()

◆ v_TimeIntegrate()

◆ v_UpdateSolutionVector()

Friends And Related Function Documentation

◆ operator<< [1/2]

◆ operator<< [2/2]

Member Data Documentation

◆ className

◆ m_Ahats

◆ m_AhattJ

◆ m_alpha

◆ m_base

◆ m_deltaT

◆ m_expFactor

◆ m_F

◆ m_freeParams

◆ m_integral_classes

◆ m_J

◆ m_Lmax

◆ m_maxTimeSteps

◆ m_mu0

◆ m_name

◆ m_npoints

◆ m_nQuadPts

◆ m_nu

◆ m_nvars

◆ m_op

◆ m_order

◆ m_qml

◆ m_schemeType

◆ m_sigma

◆ m_T

◆ m_taus

◆ m_u

◆ m_u0

◆ m_uInt

◆ m_uNext

◆ m_variant