#include <TimeIntegrationAlgorithmGLM.h>

Public Member Functions
	TimeIntegrationAlgorithmGLM (const TimeIntegrationSchemeGLM *parent)

	~TimeIntegrationAlgorithmGLM ()

void	InitializeScheme (const TimeIntegrationSchemeOperators &op)

LUE TimeIntegrationSolutionGLMSharedPtr	InitializeData (const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time)
	This function initialises the time integration scheme. More...

LUE ConstDoubleArray &	TimeIntegrate (const NekDouble deltaT, TimeIntegrationSolutionGLMSharedPtr &y)
	Explicit integration of an ODE. More...

LUE void	TimeIntegrate (const NekDouble deltaT, ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new, SingleArray &t_new)

LUE void	CheckAndVerify ()

TimeIntegrationSchemeType	GetIntegrationSchemeType () const

size_t	GetNmultiStepValues () const

size_t	GetNmultiStepImplicitDerivs () const

size_t	GetNmultiStepExplicitDerivs () const

const Array< OneD, const size_t > &	GetTimeLevelOffset () const

Public Attributes
const TimeIntegrationSchemeGLM *	m_parent {nullptr}
	Parent scheme object. More...

std::string	m_name

std::string	m_variant

size_t	m_order {0}

std::vector< NekDouble >	m_freeParams

TimeIntegrationSchemeType	m_schemeType {eNoTimeIntegrationSchemeType}

size_t	m_numstages {0}

size_t	m_numsteps {0}

size_t	m_numMultiStepValues {0}

size_t	m_numMultiStepImplicitDerivs {0}

size_t	m_numMultiStepExplicitDerivs {0}

Array< OneD, size_t >	m_timeLevelOffset

Array< OneD, Array< TwoD, NekDouble > >	m_A

Array< OneD, Array< TwoD, NekDouble > >	m_B

Array< TwoD, NekDouble >	m_U

Array< TwoD, NekDouble >	m_V

Array< OneD, std::complex< NekDouble > >	m_L

Array< OneD, Array< TwoD, NekDouble > >	m_A_phi

Array< OneD, Array< TwoD, NekDouble > >	m_B_phi

Array< OneD, Array< TwoD, NekDouble > >	m_U_phi

Array< OneD, Array< TwoD, NekDouble > >	m_V_phi

bool	m_initialised {false}

NekDouble	m_lastDeltaT {0}

NekDouble	m_lastNVars {0}
	Last delta T. More...

size_t	m_nvars
	Last number of vars. More...

size_t	m_npoints
	The number of variables in integration scheme. More...

Private Member Functions
void	CheckIfFirstStageEqualsOldSolution ()
	Optimisation-flag. More...

void	CheckIfLastStageEqualsNewSolution ()

void	VerifyIntegrationSchemeType ()

bool	CheckTimeIntegrateArguments (ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new, SingleArray &t_new) const

NekDouble	A (const size_t i, const size_t j) const

NekDouble	A (const size_t k, const size_t i, const size_t j) const

NekDouble	B (const size_t i, const size_t j) const

NekDouble	B (const size_t k, const size_t i, const size_t j) const

NekDouble	U (const size_t i, const size_t j) const

NekDouble	U (const size_t k, const size_t i, const size_t j) const

NekDouble	V (const size_t i, const size_t j) const

NekDouble	V (const size_t k, const size_t i, const size_t j) const

NekDouble	A_IMEX (const size_t i, const size_t j) const

NekDouble	B_IMEX (const size_t i, const size_t j) const

size_t	GetNsteps (void) const

size_t	GetNstages (void) const

size_t	GetFirstDim (ConstTripleArray &y) const

size_t	GetSecondDim (ConstTripleArray &y) const

Private Attributes
TimeIntegrationSchemeOperators	m_op

DoubleArray	m_Y

DoubleArray	m_tmp
	Array containing the stage values. More...

TripleArray	m_F
	Explicit RHS of each stage equation. More...

TripleArray	m_F_IMEX
	Array corresponding to the stage Derivatives. More...

NekDouble	m_T {0}
	Used to store the Explicit stage derivative of IMEX schemes. More...

bool	m_firstStageEqualsOldSolution {false}
	ime at the different stages More...

bool	m_lastStageEqualsNewSolution {false}
	Optimisation-flag. More...

Friends
LUE friend std::ostream &	operator<< (std::ostream &os, const TimeIntegrationScheme &rhs)
	The size of inner data which is stored for reuse. More...

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

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

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

Detailed Description

Definition at line 50 of file TimeIntegrationAlgorithmGLM.h.

Constructor & Destructor Documentation

◆ TimeIntegrationAlgorithmGLM()

Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::TimeIntegrationAlgorithmGLM ( const TimeIntegrationSchemeGLM * parent )

inline

Definition at line 53 of file TimeIntegrationAlgorithmGLM.h.

        : m_parent(parent)
    {
    }

◆ ~TimeIntegrationAlgorithmGLM()

Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::~TimeIntegrationAlgorithmGLM ( )

inline

Definition at line 58 of file TimeIntegrationAlgorithmGLM.h.

59 {

60 }

Member Function Documentation

◆ A() [1/2]

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::A	(	const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 267 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_A[0][i][j];
    }

References m_A.

Referenced by InitializeData(), and TimeIntegrate().

◆ A() [2/2]

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::A	(	const size_t	k,
		const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 274 of file TimeIntegrationAlgorithmGLM.h.

    {
        if (m_schemeType == eExponential)
        {
            // For exponential schemes, the parameter A is specific for
            // each variable k
            return m_A_phi[k][i][j];
        }
        else
        {
            return m_A[0][i][j];
        }
    }

References Nektar::LibUtilities::eExponential, m_A, m_A_phi, and m_schemeType.

◆ A_IMEX()

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::A_IMEX	(	const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 351 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_A[1][i][j];
    }

References m_A.

Referenced by TimeIntegrate().

◆ B() [1/2]

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::B	(	const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 288 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_B[0][i][j];
    }

References m_B.

Referenced by TimeIntegrate().

◆ B() [2/2]

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::B	(	const size_t	k,
		const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 295 of file TimeIntegrationAlgorithmGLM.h.

    {
        if (m_schemeType == eExponential)
        {
            // For exponential schemes, the parameter B is specific for
            // each variable k
            return m_B_phi[k][i][j];
        }
        else
        {
            return m_B[0][i][j];
        }
    }

References Nektar::LibUtilities::eExponential, m_B, m_B_phi, and m_schemeType.

◆ B_IMEX()

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::B_IMEX	(	const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 356 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_B[1][i][j];
    }

References m_B.

Referenced by TimeIntegrate().

◆ CheckAndVerify()

LUE void Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::CheckAndVerify ( )

inline

Definition at line 134 of file TimeIntegrationAlgorithmGLM.h.

    {
        CheckIfFirstStageEqualsOldSolution();
        CheckIfLastStageEqualsNewSolution();
        VerifyIntegrationSchemeType();
    };

References CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), and VerifyIntegrationSchemeType().

◆ CheckIfFirstStageEqualsOldSolution()

void Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::CheckIfFirstStageEqualsOldSolution ( )

private

Optimisation-flag.

Definition at line 771 of file TimeIntegrationAlgorithmGLM.cpp.

{
    // First stage equals old solution if:
    // 1. the first row of the coefficient matrix A consists of zeros
    // 2. U[0][0] is equal to one and all other first row entries of U are zero
 
    // 1. Check first condition
    for (size_t m = 0; m < m_A.size(); m++)
    {
        for (size_t i = 0; i < m_numstages; i++)
        {
            if (fabs(m_A[m][0][i]) > NekConstants::kNekZeroTol)
            {
                m_firstStageEqualsOldSolution = false;
                return;
            }
        }
    }
 
    // 2. Check second condition
    if (fabs(m_U[0][0] - 1.0) > NekConstants::kNekZeroTol)
    {
        m_firstStageEqualsOldSolution = false;
        return;
    }
 
    for (size_t i = 1; i < m_numsteps; i++)
    {
        if (fabs(m_U[0][i]) > NekConstants::kNekZeroTol)
        {
            m_firstStageEqualsOldSolution = false;
            return;
        }
    }
 
    m_firstStageEqualsOldSolution = true;
}

References Nektar::NekConstants::kNekZeroTol, m_A, m_firstStageEqualsOldSolution, m_numstages, m_numsteps, and m_U.

Referenced by CheckAndVerify().

◆ CheckIfLastStageEqualsNewSolution()

void Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::CheckIfLastStageEqualsNewSolution ( )

private

Definition at line 809 of file TimeIntegrationAlgorithmGLM.cpp.

{
    // Last stage equals new solution if:
    // 1. the last row of the coefficient matrix A is equal to the first row of
    // matrix B
    // 2. the last row of the coefficient matrix U is equal to the first row of
    // matrix V
 
    // 1. Check first condition
    for (size_t m = 0; m < m_A.size(); m++)
    {
        for (size_t i = 0; i < m_numstages; i++)
        {
            if (fabs(m_A[m][m_numstages - 1][i] - m_B[m][0][i]) >
                NekConstants::kNekZeroTol)
            {
                m_lastStageEqualsNewSolution = false;
                return;
            }
        }
    }
 
    // 2. Check second condition
    for (size_t i = 0; i < m_numsteps; i++)
    {
        if (fabs(m_U[m_numstages - 1][i] - m_V[0][i]) >
            NekConstants::kNekZeroTol)
        {
            m_lastStageEqualsNewSolution = false;
            return;
        }
    }
 
    m_lastStageEqualsNewSolution = true;
}

References Nektar::NekConstants::kNekZeroTol, m_A, m_B, m_lastStageEqualsNewSolution, m_numstages, m_numsteps, m_U, and m_V.

Referenced by CheckAndVerify().

◆ CheckTimeIntegrateArguments()

bool Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::CheckTimeIntegrateArguments	(	ConstTripleArray &	y_old,
		ConstSingleArray &	t_old,
		TripleArray &	y_new,
		SingleArray &	t_new
	)		const

private

Definition at line 903 of file TimeIntegrationAlgorithmGLM.cpp.

{
    // Check if arrays are all of consistent size
 
    ASSERTL1(y_old.size() == m_numsteps, "Non-matching number of steps.");
    ASSERTL1(y_new.size() == m_numsteps, "Non-matching number of steps.");
 
    ASSERTL1(y_old[0].size() == y_new[0].size(),
             "Non-matching number of variables.");
    ASSERTL1(y_old[0][0].size() == y_new[0][0].size(),
             "Non-matching number of coefficients.");
 
    ASSERTL1(t_old.size() == m_numsteps, "Non-matching number of steps.");
    ASSERTL1(t_new.size() == m_numsteps, "Non-matching number of steps.");
 
    return true;
}

References ASSERTL1, and m_numsteps.

Referenced by TimeIntegrate().

◆ GetFirstDim()

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::GetFirstDim ( ConstTripleArray & y ) const

inlineprivate

Definition at line 371 of file TimeIntegrationAlgorithmGLM.h.

    {
        return y[0].size();
    }

Referenced by TimeIntegrate().

◆ GetIntegrationSchemeType()

TimeIntegrationSchemeType Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::GetIntegrationSchemeType ( ) const

inline

Definition at line 141 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_schemeType;
    }

References m_schemeType.

◆ GetNmultiStepExplicitDerivs()

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::GetNmultiStepExplicitDerivs ( ) const

inline

Definition at line 156 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_numMultiStepExplicitDerivs;
    }

References m_numMultiStepExplicitDerivs.

Referenced by Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetNexplicitderivs(), InitializeData(), and TimeIntegrate().

◆ GetNmultiStepImplicitDerivs()

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::GetNmultiStepImplicitDerivs ( ) const

inline

Definition at line 151 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_numMultiStepImplicitDerivs;
    }

References m_numMultiStepImplicitDerivs.

Referenced by Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetExplicitDerivative(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetImplicitDerivative(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetNimplicitderivs(), InitializeData(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::RotateSolutionVector(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::SetExplicitDerivative(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::SetImplicitDerivative(), and TimeIntegrate().

◆ GetNmultiStepValues()

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::GetNmultiStepValues ( ) const

inline

Definition at line 146 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_numMultiStepValues;
    }

References m_numMultiStepValues.

◆ GetNstages()

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::GetNstages ( void ) const

inlineprivate

Definition at line 366 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_numstages;
    }

References m_numstages.

◆ GetNsteps()

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::GetNsteps ( void ) const

inlineprivate

Definition at line 361 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_numsteps;
    }

References m_numsteps.

Referenced by TimeIntegrate().

◆ GetSecondDim()

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::GetSecondDim ( ConstTripleArray & y ) const

inlineprivate

Definition at line 376 of file TimeIntegrationAlgorithmGLM.h.

    {
        return y[0][0].size();
    }

Referenced by TimeIntegrate().

◆ GetTimeLevelOffset()

const Array< OneD, const size_t > & Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::GetTimeLevelOffset ( ) const

inline

Definition at line 161 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_timeLevelOffset;
    }

References m_timeLevelOffset.

Referenced by Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetExplicitDerivative(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetImplicitDerivative(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetTimeLevelOffset(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetValue(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetValueTime(), InitializeData(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::SetExplicitDerivative(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::SetImplicitDerivative(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::SetValue(), and TimeIntegrate().

◆ InitializeData()

TimeIntegrationSolutionGLMSharedPtr Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::InitializeData	(	const NekDouble	deltaT,
		ConstDoubleArray &	y_0,
		const NekDouble	time
	)

This function initialises the time integration scheme.

Given the solution at the initial time level \(\boldsymbol{y}(t_0)\), this function generates the vectors \(\boldsymbol{y}^{[n]}\) and \(t^{[n]}\) needed to evaluate the time integration scheme formulated as a General Linear Method. These vectors are embedded in an object of the class TimeIntegrationSolutionGLM. This class is the abstraction of the input and output vectors of the General Linear Method.

For single-step methods, this function is trivial as it just wraps a TimeIntegrationSolutionGLM object around the given initial values and initial time. However, for multistep methods, actual time stepping is being done to evaluate the necessary parameters at multiple time levels needed to start the actual integration.

Parameters

timestep	The size of the timestep, i.e. \(\Delta t\).
time	on input: the initial time, i.e. \(t_0\).
time	on output: the new time-level after initialisation, in general this yields \(t = t_0 + (r-1)\Delta t\) where \(r\) is the number of steps of the multi-step method.
nsteps	on output: he number of initialisation steps required. In general this corresponds to \(r-1\) where \(r\) is the number of steps of the multi-step method.
f	an object of the class FuncType, where FuncType should have a method FuncType::ODEforcing to evaluate the right hand side \(f(t,\boldsymbol{y})\) of the ODE.
y_0	the initial value \(\boldsymbol{y}(t_0)\)

Returns: An object of the class TimeIntegrationSolutionGLM which represents the vectors \(\boldsymbol{y}^{[n]}\) and \(t^{[n]}\) that can be used to start the actual integration.

Definition at line 55 of file TimeIntegrationAlgorithmGLM.cpp.

{
    // Create a TimeIntegrationSolutionGLM object based upon the initial value
    // y_0. Initialise all other multi-step values and derivatives to zero.
    TimeIntegrationSolutionGLMSharedPtr y_out =
        MemoryManager<TimeIntegrationSolutionGLM>::AllocateSharedPtr(
            this, y_0, time, deltaT);
 
    if (m_schemeType == eExplicit || m_schemeType == eExponential)
    {
        // Ensure initial solution is in correct space.
        m_op.DoProjection(y_out->GetSolution(), y_out->UpdateSolution(), time);
    }
    else if (m_schemeType == eDiagonallyImplicit &&
             fabs(A(0, 0)) < NekConstants::kNekZeroTol)
    {
        // Explicit first-stage when first diagonal coefficient is equal to
        // zero (EDIRK/ESDIRK schemes).
        m_op.DoProjection(y_out->GetSolution(), y_out->UpdateSolution(), time);
    }
 
    // Calculate the initial derivative, if is part of the solution vector of
    // the current scheme.
    if (GetNmultiStepExplicitDerivs() > 0)
    {
        const Array<OneD, const size_t> offsets = GetTimeLevelOffset();
 
        if (offsets[GetNmultiStepValues() + GetNmultiStepImplicitDerivs()] == 0)
        {
            size_t nvar    = y_0.size();
            size_t npoints = y_0[0].size();
            DoubleArray f_y_0(nvar);
            for (size_t i = 0; i < nvar; i++)
            {
                f_y_0[i] = Array<OneD, NekDouble>(npoints);
            }
 
            // Calculate the derivative of the initial value.
            m_op.DoOdeRhs(y_out->GetSolution(), f_y_0, time);
 
            // Multiply by the step size.
            for (size_t i = 0; i < nvar; i++)
            {
                Vmath::Smul(npoints, deltaT, f_y_0[i], 1, f_y_0[i], 1);
            }
            y_out->SetExplicitDerivative(0, f_y_0, deltaT);
        }
    }
 
    return y_out;
}

References A(), Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoProjection(), Nektar::LibUtilities::eDiagonallyImplicit, Nektar::LibUtilities::eExplicit, Nektar::LibUtilities::eExponential, GetNmultiStepExplicitDerivs(), GetNmultiStepImplicitDerivs(), GetNmultiStepValues(), GetTimeLevelOffset(), Nektar::NekConstants::kNekZeroTol, m_op, m_schemeType, and Vmath::Smul().

◆ InitializeScheme()

void Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::InitializeScheme ( const TimeIntegrationSchemeOperators & op )

Definition at line 49 of file TimeIntegrationAlgorithmGLM.cpp.

{
    m_op = op;
}

References m_op.

◆ TimeIntegrate() [1/2]

void Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::TimeIntegrate	(	const NekDouble	deltaT,
		ConstTripleArray &	y_old,
		ConstSingleArray &	t_old,
		TripleArray &	y_new,
		SingleArray &	t_new
	)

Definition at line 430 of file TimeIntegrationAlgorithmGLM.cpp.

{
    ASSERTL1(CheckTimeIntegrateArguments(y_old, t_old, y_new, t_new),
             "Arguments not well defined");
 
    TimeIntegrationSchemeType type = m_schemeType;
 
    // Check if storage has already been initialised. If so, we just zero the
    // temporary storage.
    if (m_initialised && m_nvars == GetFirstDim(y_old) &&
        m_npoints == GetSecondDim(y_old))
    {
        for (size_t j = 0; j < m_nvars; j++)
        {
            Vmath::Zero(m_npoints, m_tmp[j], 1);
        }
    }
    else
    {
        m_nvars   = GetFirstDim(y_old);
        m_npoints = GetSecondDim(y_old);
 
        // First, calculate the various stage values and stage derivatives
        // (this is the multi-stage part of the method)
        // - m_Y   corresponds to the stage values
        // - m_F   corresponds to the stage derivatives
        // - m_T   corresponds to the time at the different stages
        // - m_tmp corresponds to the explicit right hand side of each stage
        //         equation
        //   (for explicit schemes, this correspond to m_Y)
 
        // Allocate memory for the arrays m_Y and m_F and m_tmp. The same
        // storage will be used for every stage -> m_Y is a DoubleArray.
        m_tmp = DoubleArray(m_nvars);
        for (size_t j = 0; j < m_nvars; j++)
        {
            m_tmp[j] = Array<OneD, NekDouble>(m_npoints, 0.0);
        }
 
        // The same storage will be used for every stage -> m_tmp is a
        // DoubleArray.
        if (type == eExplicit || type == eExponential)
        {
            m_Y = m_tmp;
        }
        else
        {
            m_Y = DoubleArray(m_nvars);
            for (size_t j = 0; j < m_nvars; j++)
            {
                m_Y[j] = Array<OneD, NekDouble>(m_npoints, 0.0);
            }
        }
 
        // Different storage for every stage derivative as the data
        // will be re-used to update the solution -> m_F is a TripleArray.
        m_F = TripleArray(m_numstages);
        for (size_t i = 0; i < m_numstages; ++i)
        {
            m_F[i] = DoubleArray(m_nvars);
            for (size_t j = 0; j < m_nvars; j++)
            {
                m_F[i][j] = Array<OneD, NekDouble>(m_npoints, 0.0);
            }
        }
 
        if (type == eIMEX)
        {
            m_F_IMEX = TripleArray(m_numstages);
            for (size_t i = 0; i < m_numstages; ++i)
            {
                m_F_IMEX[i] = DoubleArray(m_nvars);
                for (size_t j = 0; j < m_nvars; j++)
                {
                    m_F_IMEX[i][j] = Array<OneD, NekDouble>(m_npoints, 0.0);
                }
            }
        }
 
        // Finally, flag indicating that the memory has been initialised.
        m_initialised = true;
 
    } // end else
 
    // For an exponential integrator if the time increment or the number of
    // variables has changed then the exponenial matrices must be recomputed.
 
    // NOTE: The eigenvalues are not yet seeded!
    if (type == eExponential)
    {
        if (m_lastDeltaT != deltaT || m_lastNVars != GetFirstDim(y_old))
        {
            m_parent->InitializeSecondaryData(this, deltaT);
 
            m_lastDeltaT = deltaT;
            m_lastNVars  = GetFirstDim(y_old);
        }
    }
 
    // The loop below calculates the stage values and derivatives.
    for (size_t stage = 0; stage < m_numstages; stage++)
    {
        if ((stage == 0) && m_firstStageEqualsOldSolution)
        {
            for (size_t k = 0; k < m_nvars; k++)
            {
                Vmath::Vcopy(m_npoints, y_old[0][k], 1, m_Y[k], 1);
            }
        }
        else
        {
            // The stage values m_Y are a linear combination of:
            // 1: The stage derivatives:
            if (stage != 0)
            {
                for (size_t k = 0; k < m_nvars; k++)
                {
                    Vmath::Smul(m_npoints, deltaT * A(k, stage, 0), m_F[0][k],
                                1, m_tmp[k], 1);
 
                    if (type == eIMEX)
                    {
                        Vmath::Svtvp(m_npoints, deltaT * A_IMEX(stage, 0),
                                     m_F_IMEX[0][k], 1, m_tmp[k], 1, m_tmp[k],
                                     1);
                    }
                }
            }
 
            for (size_t j = 1; j < stage; j++)
            {
                for (size_t k = 0; k < m_nvars; k++)
                {
                    Vmath::Svtvp(m_npoints, deltaT * A(k, stage, j), m_F[j][k],
                                 1, m_tmp[k], 1, m_tmp[k], 1);
 
                    if (type == eIMEX)
                    {
                        Vmath::Svtvp(m_npoints, deltaT * A_IMEX(stage, j),
                                     m_F_IMEX[j][k], 1, m_tmp[k], 1, m_tmp[k],
                                     1);
                    }
                }
            }
 
            // 2: The imported multi-step solution of the previous time level:
            for (size_t j = 0; j < m_numsteps; j++)
            {
                for (size_t k = 0; k < m_nvars; k++)
                {
                    Vmath::Svtvp(m_npoints, U(k, stage, j), y_old[j][k], 1,
                                 m_tmp[k], 1, m_tmp[k], 1);
                }
            }
        } // end else
 
        // Calculate the stage derivative based upon the stage value.
        if (type == eDiagonallyImplicit)
        {
            if (m_numstages == 1)
            {
                m_T = t_old[0] + deltaT;
            }
            else
            {
                m_T = t_old[0];
                for (size_t j = 0; j <= stage; ++j)
                {
                    m_T += A(stage, j) * deltaT;
                }
            }
 
            // Explicit first-stage when first diagonal coefficient is equal to
            // zero (EDIRK/ESDIRK schemes).
            if ((stage == 0) && (fabs(A(0, 0)) < NekConstants::kNekZeroTol))
            {
                m_op.DoOdeRhs(m_Y, m_F[stage], t_old[0]);
            }
            // Otherwise, the stage is updated implicitly.
            else
            {
                m_op.DoImplicitSolve(m_tmp, m_Y, m_T, A(stage, stage) * deltaT);
 
                for (size_t k = 0; k < m_nvars; ++k)
                {
                    Vmath::Vsub(m_npoints, m_Y[k], 1, m_tmp[k], 1,
                                m_F[stage][k], 1);
                    Vmath::Smul(m_npoints, 1.0 / (A(stage, stage) * deltaT),
                                m_F[stage][k], 1, m_F[stage][k], 1);
                }
            }
        }
        else if (type == eIMEX)
        {
            if (m_numstages == 1)
            {
                m_T = t_old[0] + deltaT;
            }
            else
            {
                m_T = t_old[0];
                for (size_t j = 0; j <= stage; ++j)
                {
                    m_T += A(stage, j) * deltaT;
                }
            }
 
            if (fabs(A(stage, stage)) > NekConstants::kNekZeroTol)
            {
                m_op.DoImplicitSolve(m_tmp, m_Y, m_T, A(stage, stage) * deltaT);
 
                for (size_t k = 0; k < m_nvars; k++)
                {
                    Vmath::Vsub(m_npoints, m_Y[k], 1, m_tmp[k], 1,
                                m_F[stage][k], 1);
                    Vmath::Smul(m_npoints, 1.0 / (A(stage, stage) * deltaT),
                                m_F[stage][k], 1, m_F[stage][k], 1);
                }
            }
 
            m_op.DoOdeRhs(m_Y, m_F_IMEX[stage], m_T);
        }
        else if (type == eExplicit || type == eExponential)
        {
            if (m_numstages == 1)
            {
                m_T = t_old[0];
            }
            else
            {
                m_T = t_old[0];
                for (size_t j = 0; j < stage; ++j)
                {
                    m_T += A(stage, j) * deltaT;
                }
            }
 
            // Avoid projecting the same solution twice.
            if (!((stage == 0) && m_firstStageEqualsOldSolution))
            {
                // Ensure solution is in correct space.
                m_op.DoProjection(m_Y, m_Y, m_T);
            }
 
            m_op.DoOdeRhs(m_Y, m_F[stage], m_T);
        }
    }
 
    // Next, the solution vector y at the new time level will be calculated.
    //
    // For multi-step methods, this includes updating the values of the
    // auxiliary parameters.
    //
    // The loop below calculates the solution at the new time level.
    //
    // If last stage equals the new solution, the new solution needs not be
    // calculated explicitly but can simply be copied. This saves a solve.
    size_t i_start = 0;
    if (m_lastStageEqualsNewSolution)
    {
        for (size_t k = 0; k < m_nvars; k++)
        {
            Vmath::Vcopy(m_npoints, m_Y[k], 1, y_new[0][k], 1);
        }
 
        t_new[0] = t_old[0] + deltaT;
 
        i_start = 1;
    }
 
    for (size_t i = i_start; i < m_numsteps; i++)
    {
        // The solution at the new time level is a linear combination of:
        // 1: The stage derivatives:
        for (size_t k = 0; k < m_nvars; k++)
        {
            Vmath::Smul(m_npoints, deltaT * B(k, i, 0), m_F[0][k], 1,
                        y_new[i][k], 1);
 
            if (type == eIMEX)
            {
                Vmath::Svtvp(m_npoints, deltaT * B_IMEX(i, 0), m_F_IMEX[0][k],
                             1, y_new[i][k], 1, y_new[i][k], 1);
            }
        }
 
        if (m_numstages != 1 || type != eIMEX)
        {
            t_new[i] = B(i, 0) * deltaT;
        }
 
        for (size_t j = 1; j < m_numstages; j++)
        {
            for (size_t k = 0; k < m_nvars; k++)
            {
                Vmath::Svtvp(m_npoints, deltaT * B(k, i, j), m_F[j][k], 1,
                             y_new[i][k], 1, y_new[i][k], 1);
 
                if (type == eIMEX)
                {
                    Vmath::Svtvp(m_npoints, deltaT * B_IMEX(i, j),
                                 m_F_IMEX[j][k], 1, y_new[i][k], 1, y_new[i][k],
                                 1);
                }
            }
 
            if (m_numstages != 1 || type != eIMEX)
            {
                t_new[i] += B(i, j) * deltaT;
            }
        }
 
        // 2: The imported multi-step solution of the previous time level:
        for (size_t j = 0; j < m_numsteps; j++)
        {
            for (size_t k = 0; k < m_nvars; k++)
            {
                Vmath::Svtvp(m_npoints, V(k, i, j), y_old[j][k], 1, y_new[i][k],
                             1, y_new[i][k], 1);
            }
 
            if (m_numstages != 1 || type != eIMEX)
            {
                t_new[i] += V(i, j) * t_old[j];
            }
        }
    }
 
    // Ensure that the new solution is projected if necessary.
    if (type == eExplicit || type == eExponential ||
        fabs(m_T - t_new[0]) > NekConstants::kNekZeroTol)
    {
        m_op.DoProjection(y_new[0], y_new[0], t_new[0]);
    }
 
} // end TimeIntegrate()

◆ TimeIntegrate() [2/2]

ConstDoubleArray & Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::TimeIntegrate	(	const NekDouble	deltaT,
		TimeIntegrationSolutionGLMSharedPtr &	y
	)

Explicit integration of an ODE.

This function explicitely perfroms a single integration step of the ODE system:

\[ \frac{d\boldsymbol{y}}{dt}=\boldsymbol{f}(t,\boldsymbol{y}) \]

Parameters

deltaT	The size of the timestep, i.e. \(\Delta t\).
f	an object of the class FuncType, where FuncType should have a method FuncType::ODEforcing to evaluate the right hand side \(f(t,\boldsymbol{y})\) of the ODE.
y	on input: the vectors \(\boldsymbol{y}^{[n-1]}\) and \(t^{[n-1]}\) (which corresponds to the solution at the old time level)
y	on output: the vectors \(\boldsymbol{y}^{[n]}\) and \(t^{[n]}\) (which corresponds to the solution at the old new level)

Returns: The actual solution \(\boldsymbol{y}^{n}\) at the new time level (which in fact is also embedded in the argument y).

Definition at line 108 of file TimeIntegrationAlgorithmGLM.cpp.

{
    size_t nvar    = solvector->GetFirstDim();
    size_t npoints = solvector->GetSecondDim();
 
    if (solvector->GetIntegrationSchemeData() != this)
    {
        // This branch will be taken when the solution vector (solvector) is set
        // up for a different scheme than the object this method is called from.
        // (typically needed to calculate the first time-levels of a multi-step
        // scheme).
 
        // To do this kind of 'non-matching' integration, we perform the
        // following three steps:
        //
        // 1: Copy the required input information from the solution vector of
        //    the master scheme to the input solution vector of the current
        //    scheme.
        //
        // 2: Time-integrate for one step using the current scheme.
        //
        // 3: Copy the information contained in the output vector of the current
        //    scheme to the solution vector of the master scheme.
 
        // STEP 1: Copy the required input information from the solution
        //         vector of the master scheme to the input solution vector
        //         of the current scheme.
 
        // 1.1 Determine which information is required for the current scheme.
 
        // Number of required values of the current scheme.
        size_t nCurSchemeVals = GetNmultiStepValues();
 
        // Number of required implicit derivs of the current scheme.
        size_t nCurSchemeImpDers = GetNmultiStepImplicitDerivs();
 
        // Number of required explicit derivs of the current scheme.
        size_t nCurSchemeExpDers = GetNmultiStepExplicitDerivs();
 
        // Number of steps in the current scheme.
        size_t nCurSchemeSteps = GetNsteps();
 
        // Number of values of the master scheme.
        size_t nMasterSchemeVals = solvector->GetNvalues();
 
        // Number of implicit derivs of the master scheme.
        size_t nMasterSchemeImpDers = solvector->GetNimplicitderivs();
 
        // Number of explicit derivs of the master scheme.
        size_t nMasterSchemeExpDers = solvector->GetNexplicitderivs();
 
        // Array indicating which time-level the values and derivatives of the
        // current schemes belong.
        const Array<OneD, const size_t> &curTimeLevels = GetTimeLevelOffset();
 
        // Array indicating which time-level the values and derivatives of the
        // master schemes belong.
        const Array<OneD, const size_t> &masterTimeLevels =
            solvector->GetTimeLevelOffset();
 
        // 1.2 Copy the required information from the master solution vector to
        //     the input solution vector of the current scheme.
 
        NekDouble t_n = 0.0;
        DoubleArray y_n;
        DoubleArray dtFy_n;
 
        // Input solution vector of the current scheme.
        TimeIntegrationSolutionGLMSharedPtr solvector_in =
            MemoryManager<TimeIntegrationSolutionGLM>::AllocateSharedPtr(this);
 
        // Set solution.
        for (size_t n = 0; n < nCurSchemeVals; n++)
        {
            // Get the required value out of the master solution vector.
            y_n = solvector->GetValue(curTimeLevels[n]);
            t_n = solvector->GetValueTime(curTimeLevels[n]);
 
            // Set the required value in the input solution vector of the
            // current scheme.
            solvector_in->SetValue(curTimeLevels[n], y_n, t_n);
        }
 
        // Set implicit derivatives.
        for (size_t n = nCurSchemeVals; n < nCurSchemeVals + nCurSchemeImpDers;
             n++)
        {
            // Get the required derivative out of the master solution vector.
            dtFy_n = solvector->GetImplicitDerivative(curTimeLevels[n]);
 
            // Set the required derivative in the input solution vector of the
            // current scheme.
            solvector_in->SetImplicitDerivative(curTimeLevels[n], dtFy_n,
                                                deltaT);
        }
 
        // Set explicit derivatives.
        for (size_t n = nCurSchemeVals + nCurSchemeImpDers; n < nCurSchemeSteps;
             n++)
        {
            // Get the required derivative out of the master solution vector.
            dtFy_n = solvector->GetExplicitDerivative(curTimeLevels[n]);
 
            // Set the required derivative in the input solution vector of the
            // current scheme.
            solvector_in->SetExplicitDerivative(curTimeLevels[n], dtFy_n,
                                                deltaT);
        }
 
        // STEP 2: Time-integrate for one step using the current scheme.
 
        // Output solution vector of the current scheme.
        TimeIntegrationSolutionGLMSharedPtr solvector_out =
            MemoryManager<TimeIntegrationSolutionGLM>::AllocateSharedPtr(
                this, nvar, npoints);
 
        // Integrate one step.
        TimeIntegrate(deltaT, solvector_in->GetSolutionVector(),
                      solvector_in->GetTimeVector(),
                      solvector_out->UpdateSolutionVector(),
                      solvector_out->UpdateTimeVector());
 
        // STEP 3: Copy the information contained in the output vector of the
        //         current scheme to the solution vector of the master scheme.
 
        // 3.1 Check whether the current time scheme updates the most recent
        //     derivative that should be updated in the master scheme.  If not,
        //     calculate the derivative. This can be done based upon the
        //     corresponding value and the DoOdeRhs operator.
 
        // This flag indicates whether the new implicit derivative is available
        // in the output of the current scheme or whether it should be
        // calculated.
        bool CalcNewImpDeriv = false;
 
        if (nMasterSchemeImpDers > 0)
        {
            if (nCurSchemeImpDers == 0 || (masterTimeLevels[nMasterSchemeVals] <
                                           curTimeLevels[nCurSchemeVals]))
            {
                CalcNewImpDeriv = true;
            }
        }
 
        // This flag indicates whether the new explicit derivative is available
        // in the output of the current scheme or whether it should be
        // calculated.
        bool CalcNewExpDeriv = false;
 
        if (nMasterSchemeExpDers > 0)
        {
            if (nCurSchemeExpDers == 0 ||
                (masterTimeLevels[nMasterSchemeVals + nMasterSchemeImpDers] <
                 curTimeLevels[nCurSchemeVals + nCurSchemeImpDers]))
            {
                CalcNewExpDeriv = true;
            }
        }
 
        // Indicates the time level at which the implicit derivative of the
        // master scheme has to be updated.
        size_t newImpDerivTimeLevel =
            (masterTimeLevels.size() > nMasterSchemeVals)
                ? masterTimeLevels[nMasterSchemeVals]
                : 0;
 
        DoubleArray f_impn(nvar);
 
        // Calculate implicit derivatives.
        if (CalcNewImpDeriv)
        {
            if (newImpDerivTimeLevel == 0 || newImpDerivTimeLevel == 1)
            {
                y_n = solvector->GetValue(newImpDerivTimeLevel);
                t_n = solvector->GetValueTime(newImpDerivTimeLevel);
            }
            else
            {
                ASSERTL1(false, "Problems with initialising scheme");
            }
 
            for (size_t j = 0; j < nvar; j++)
            {
                f_impn[j] = Array<OneD, NekDouble>(npoints);
            }
 
            // Calculate the derivative of the initial value.
            m_op.DoImplicitSolve(y_n, f_impn, t_n + deltaT, deltaT);
            for (size_t j = 0; j < nvar; j++)
            {
                Vmath::Vsub(m_npoints, f_impn[j], 1, y_n[j], 1, f_impn[j], 1);
            }
        }
 
        // Indiciates the time level at which the explicit derivative of the
        // master scheme has to be updated.
        size_t newExpDerivTimeLevel =
            (masterTimeLevels.size() > nMasterSchemeVals + nMasterSchemeImpDers)
                ? masterTimeLevels[nMasterSchemeVals + nMasterSchemeImpDers]
                : 0;
 
        DoubleArray f_expn(nvar);
 
        // Calculate explicit derivatives.
        if (CalcNewExpDeriv)
        {
            // If the time level corresponds to 0, calculate the derivative
            // based upon the solution value at the new time-level.
            if (newExpDerivTimeLevel == 0)
            {
                y_n = solvector_out->GetValue(0);
                t_n = solvector_out->GetValueTime(0);
            }
            // If the time level corresponds to 1, calculate the derivative
            // based upon the solution value at the old time-level.
            else if (newExpDerivTimeLevel == 1)
            {
                y_n = solvector->GetValue(0);
                t_n = solvector->GetValueTime(0);
            }
            else
            {
                ASSERTL1(false, "Problems with initialising scheme");
            }
 
            for (size_t j = 0; j < nvar; j++)
            {
                f_expn[j] = Array<OneD, NekDouble>(npoints);
            }
 
            // Ensure solution is in correct space.
            if (newExpDerivTimeLevel == 1)
            {
                m_op.DoProjection(y_n, y_n, t_n);
            }
 
            // Calculate the derivative.
            m_op.DoOdeRhs(y_n, f_expn, t_n);
 
            // Multiply by dt (as required by the General Linear Method
            // framework).
            for (size_t j = 0; j < nvar; j++)
            {
                Vmath::Smul(npoints, deltaT, f_expn[j], 1, f_expn[j], 1);
            }
        }
 
        // Rotate the solution vector (i.e. updating without
        // calculating/inserting new values).
        solvector->RotateSolutionVector();
 
        // 3.2 Copy the information calculated using the current scheme from
        //     the output solution vector to the master solution vector.
 
        // Set solution.
        y_n = solvector_out->GetValue(0);
        t_n = solvector_out->GetValueTime(0);
 
        // Ensure solution is in correct space.
        if (newExpDerivTimeLevel == 1)
        {
            m_op.DoProjection(y_n, y_n, t_n);
        }
 
        solvector->SetValue(0, y_n, t_n);
 
        // Set implicit derivative.
        if (CalcNewImpDeriv)
        {
            // Set the calculated derivative in the master solution vector.
            solvector->SetImplicitDerivative(newImpDerivTimeLevel, f_impn,
                                             deltaT);
        }
        else if (nCurSchemeImpDers > 0 && nMasterSchemeImpDers > 0)
        {
            // Get the calculated derivative out of the output solution vector
            // of the current scheme.
            dtFy_n = solvector_out->GetImplicitDerivative(newImpDerivTimeLevel);
 
            // Set the calculated derivative in the master solution vector.
            solvector->SetImplicitDerivative(newImpDerivTimeLevel, dtFy_n,
                                             deltaT);
        }
 
        // Set explicit derivative.
        if (CalcNewExpDeriv)
        {
            // Set the calculated derivative in the master solution vector.
            solvector->SetExplicitDerivative(newExpDerivTimeLevel, f_expn,
                                             deltaT);
        }
        else if (nCurSchemeExpDers > 0 && nMasterSchemeExpDers > 0)
        {
            // Get the calculated derivative out of the output solution vector
            // of the current scheme.
            dtFy_n = solvector_out->GetExplicitDerivative(newExpDerivTimeLevel);
 
            // Set the calculated derivative in the master solution vector.
            solvector->SetExplicitDerivative(newExpDerivTimeLevel, dtFy_n,
                                             deltaT);
        }
    }
    else
    {
        TimeIntegrationSolutionGLMSharedPtr solvector_new =
            MemoryManager<TimeIntegrationSolutionGLM>::AllocateSharedPtr(
                this, nvar, npoints);
 
        TimeIntegrate(deltaT, solvector->GetSolutionVector(),
                      solvector->GetTimeVector(),
                      solvector_new->UpdateSolutionVector(),
                      solvector_new->UpdateTimeVector());
 
        solvector = solvector_new;
    }
 
    return solvector->GetSolution();
 
} // end TimeIntegrate()

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoImplicitSolve(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoProjection(), GetNmultiStepExplicitDerivs(), GetNmultiStepImplicitDerivs(), GetNmultiStepValues(), GetNsteps(), GetTimeLevelOffset(), m_npoints, m_op, Vmath::Smul(), TimeIntegrate(), and Vmath::Vsub().

Referenced by TimeIntegrate().

◆ U() [1/2]

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::U	(	const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 309 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_U[i][j];
    }

References m_U.

Referenced by TimeIntegrate().

◆ U() [2/2]

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::U	(	const size_t	k,
		const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 316 of file TimeIntegrationAlgorithmGLM.h.

    {
        if (m_schemeType == eExponential)
        {
            // For exponential schemes, the parameter U is specific for
            // each variable k
            return m_U_phi[k][i][j];
        }
        else
        {
            return m_U[i][j];
        }
    }

References Nektar::LibUtilities::eExponential, m_schemeType, m_U, and m_U_phi.

◆ V() [1/2]

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::V	(	const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 330 of file TimeIntegrationAlgorithmGLM.h.

    {
        return m_V[i][j];
    }

References m_V.

Referenced by TimeIntegrate().

◆ V() [2/2]

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::V	(	const size_t	k,
		const size_t	i,
		const size_t	j
	)		const

inlineprivate

Definition at line 337 of file TimeIntegrationAlgorithmGLM.h.

    {
        if (m_schemeType == eExponential)
        {
            // For exponential schemes, the parameter V is specific for
            // each variable k
            return m_V_phi[k][i][j];
        }
        else
        {
            return m_V[i][j];
        }
    }

References Nektar::LibUtilities::eExponential, m_schemeType, m_V, and m_V_phi.

◆ VerifyIntegrationSchemeType()

void Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::VerifyIntegrationSchemeType ( )

private

Definition at line 845 of file TimeIntegrationAlgorithmGLM.cpp.

{
#ifdef NEKTAR_DEBUG
    size_t IMEXdim = m_A.size();
    size_t dim     = m_A[0].GetRows();
 
    Array<OneD, TimeIntegrationSchemeType> vertype(IMEXdim, eExplicit);
 
    if (m_schemeType == eExponential)
    {
        vertype[0] = eExponential;
    }
 
    for (size_t m = 0; m < IMEXdim; m++)
    {
        for (size_t i = 0; i < dim; i++)
        {
            if (fabs(m_A[m][i][i]) > NekConstants::kNekZeroTol)
            {
                vertype[m] = eDiagonallyImplicit;
            }
        }
 
        for (size_t i = 0; i < dim; i++)
        {
            for (size_t j = i + 1; j < dim; j++)
            {
                if (fabs(m_A[m][i][j]) > NekConstants::kNekZeroTol)
                {
                    vertype[m] = eImplicit;
                    ASSERTL1(false, "Fully Implicit schemes cannnot be handled "
                                    "by the TimeIntegrationScheme class");
                }
            }
        }
    }
 
    if (IMEXdim == 2)
    {
        ASSERTL1(m_B.size() == 2,
                 "Coefficient Matrix B should have an "
                 "implicit and explicit part for IMEX schemes");
 
        if ((vertype[0] == eDiagonallyImplicit) && (vertype[1] == eExplicit))
        {
            vertype[0] = eIMEX;
        }
        else
        {
            ASSERTL1(false, "This is not a proper IMEX scheme");
        }
    }
 
    ASSERTL1(vertype[0] == m_schemeType,
             "Time integration scheme coefficients do not match its type");
#endif
}

References ASSERTL1, Nektar::LibUtilities::eDiagonallyImplicit, Nektar::LibUtilities::eExplicit, Nektar::LibUtilities::eExponential, Nektar::LibUtilities::eIMEX, Nektar::LibUtilities::eImplicit, Nektar::NekConstants::kNekZeroTol, m_A, m_B, and m_schemeType.

Referenced by CheckAndVerify().

Friends And Related Function Documentation

◆ operator<< [1/4]

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

friend

Definition at line 932 of file TimeIntegrationAlgorithmGLM.cpp.

{
    size_t r                       = rhs.m_numsteps;
    size_t s                       = rhs.m_numstages;
    TimeIntegrationSchemeType type = rhs.m_schemeType;
 
    size_t oswidth     = 9;
    size_t osprecision = 6;
 
    os << "Time Integration Scheme (Master): " << rhs.m_parent->GetFullName()
       << "\n"
       << "Time Integration Phase  : " << rhs.m_name << "\n"
       << "- number of steps:  " << r << "\n"
       << "- number of stages: " << s << "\n"
       << "General linear method tableau:\n";
 
    for (size_t i = 0; i < s; i++)
    {
        for (size_t j = 0; j < s; j++)
        {
            os.width(oswidth);
            os.precision(osprecision);
            os << std::right << rhs.A(i, j) << " ";
        }
        if (type == eIMEX)
        {
            os << " '";
            for (size_t j = 0; j < s; j++)
            {
                os.width(oswidth);
                os.precision(osprecision);
                os << std::right << rhs.A_IMEX(i, j) << " ";
            }
        }
 
        os << " |";
 
        for (size_t j = 0; j < r; j++)
        {
            os.width(oswidth);
            os.precision(osprecision);
            os << std::right << rhs.U(i, j);
        }
        os << std::endl;
    }
 
    size_t imexflag = (type == eIMEX) ? 2 : 1;
    for (size_t i = 0;
         i < (r + imexflag * s) * (oswidth + 1) + imexflag * 2 - 1; i++)
    {
        os << "-";
    }
    os << std::endl;
 
    for (size_t i = 0; i < r; i++)
    {
        for (size_t j = 0; j < s; j++)
        {
            os.width(oswidth);
            os.precision(osprecision);
            os << std::right << rhs.B(i, j) << " ";
        }
        if (type == eIMEX)
        {
            os << " '";
            for (size_t j = 0; j < s; j++)
            {
                os.width(oswidth);
                os.precision(osprecision);
                os << std::right << rhs.B_IMEX(i, j) << " ";
            }
        }
 
        os << " |";
 
        for (size_t j = 0; j < r; j++)
        {
            os.width(oswidth);
            os.precision(osprecision);
            os << std::right << rhs.V(i, j);
        }
 
        os << "  |";
 
        os.width(oswidth);
        os.precision(osprecision);
        os << std::right << rhs.m_timeLevelOffset[i];
 
        if (i < rhs.m_numMultiStepValues)
        {
            os << std::right << " value";
        }
        else
        {
            os << std::right << " derivative";
        }
 
        os << std::endl;
    }
 
    if (type == eExponential)
    {
        for (size_t k = 0; k < rhs.m_nvars; k++)
        {
            os << std::endl
               << "General linear method exponential tableau for variable " << k
               << ":\n";
 
            for (size_t i = 0; i < s; i++)
            {
                for (size_t j = 0; j < s; j++)
                {
                    os.width(oswidth);
                    os.precision(osprecision);
                    os << std::right << rhs.A(k, i, j) << " ";
                }
 
                os << " |";
 
                for (size_t j = 0; j < r; j++)
                {
                    os.width(oswidth);
                    os.precision(osprecision);
                    os << std::right << rhs.B(k, i, j);
                }
                os << std::endl;
            }
 
            size_t imexflag = (type == eIMEX) ? 2 : 1;
            for (size_t i = 0;
                 i < (r + imexflag * s) * (oswidth + 1) + imexflag * 2 - 1; i++)
            {
                os << "-";
            }
            os << std::endl;
 
            for (size_t i = 0; i < r; i++)
            {
                for (size_t j = 0; j < s; j++)
                {
                    os.width(oswidth);
                    os.precision(osprecision);
                    os << std::right << rhs.B(k, i, j) << " ";
                }
 
                os << " |";
 
                for (size_t j = 0; j < r; j++)
                {
                    os.width(oswidth);
                    os.precision(osprecision);
                    os << std::right << rhs.V(k, i, j);
                }
                os << std::endl;
            }
        }
    }
 
    return os;
} // end function operator<<

◆ operator<< [2/4]

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

friend

Definition at line 926 of file TimeIntegrationAlgorithmGLM.cpp.

{
    return operator<<(os, *rhs);
}

◆ operator<< [3/4]

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

friend

The size of inner data which is stored for reuse.

Definition at line 66 of file TimeIntegrationScheme.cpp.

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

◆ operator<< [4/4]

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

friend

Definition at line 73 of file TimeIntegrationScheme.cpp.

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

Member Data Documentation

◆ m_A

Array<OneD, Array<TwoD, NekDouble> > Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_A

Definition at line 207 of file TimeIntegrationAlgorithmGLM.h.

Referenced by A(), A_IMEX(), CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), and VerifyIntegrationSchemeType().

◆ m_A_phi

Array<OneD, Array<TwoD, NekDouble> > Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_A_phi

Definition at line 215 of file TimeIntegrationAlgorithmGLM.h.

Referenced by A(), and Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_B

Array<OneD, Array<TwoD, NekDouble> > Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_B

Definition at line 208 of file TimeIntegrationAlgorithmGLM.h.

Referenced by B(), B_IMEX(), CheckIfLastStageEqualsNewSolution(), and VerifyIntegrationSchemeType().

◆ m_B_phi

Array<OneD, Array<TwoD, NekDouble> > Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_B_phi

Definition at line 216 of file TimeIntegrationAlgorithmGLM.h.

Referenced by B(), and Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_F

TripleArray Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_F

private

Explicit RHS of each stage equation.

Definition at line 248 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

◆ m_F_IMEX

TripleArray Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_F_IMEX

private

Array corresponding to the stage Derivatives.

Definition at line 249 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

◆ m_firstStageEqualsOldSolution

bool Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_firstStageEqualsOldSolution {false}

private

ime at the different stages

Definition at line 254 of file TimeIntegrationAlgorithmGLM.h.

Referenced by CheckIfFirstStageEqualsOldSolution(), and TimeIntegrate().

◆ m_freeParams

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

Definition at line 173 of file TimeIntegrationAlgorithmGLM.h.

◆ m_initialised

bool Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_initialised {false}

Definition at line 222 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

◆ m_L

Array<OneD, std::complex<NekDouble> > Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_L

Definition at line 213 of file TimeIntegrationAlgorithmGLM.h.

Referenced by Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_lastDeltaT

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_lastDeltaT {0}

Definition at line 225 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

◆ m_lastNVars

NekDouble Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_lastNVars {0}

Last delta T.

Definition at line 226 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

◆ m_lastStageEqualsNewSolution

bool Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_lastStageEqualsNewSolution {false}

private

Optimisation-flag.

Definition at line 255 of file TimeIntegrationAlgorithmGLM.h.

Referenced by CheckIfLastStageEqualsNewSolution(), and TimeIntegrate().

◆ m_name

std::string Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_name

Definition at line 170 of file TimeIntegrationAlgorithmGLM.h.

◆ m_npoints

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_npoints

The number of variables in integration scheme.

Definition at line 229 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

◆ m_numMultiStepExplicitDerivs

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_numMultiStepExplicitDerivs {0}

Definition at line 194 of file TimeIntegrationAlgorithmGLM.h.

Referenced by GetNmultiStepExplicitDerivs().

◆ m_numMultiStepImplicitDerivs

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_numMultiStepImplicitDerivs {0}

Definition at line 190 of file TimeIntegrationAlgorithmGLM.h.

Referenced by GetNmultiStepImplicitDerivs().

◆ m_numMultiStepValues

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_numMultiStepValues {0}

Definition at line 186 of file TimeIntegrationAlgorithmGLM.h.

Referenced by GetNmultiStepValues().

◆ m_numstages

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_numstages {0}

Definition at line 179 of file TimeIntegrationAlgorithmGLM.h.

Referenced by CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), GetNstages(), TimeIntegrate(), and Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_numsteps

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_numsteps {0}

Definition at line 182 of file TimeIntegrationAlgorithmGLM.h.

Referenced by CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), CheckTimeIntegrateArguments(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetExplicitDerivative(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::GetNsteps(), GetNsteps(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::RotateSolutionVector(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::SetExplicitDerivative(), TimeIntegrate(), Nektar::LibUtilities::TimeIntegrationSolutionGLM::TimeIntegrationSolutionGLM(), and Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_nvars

size_t Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_nvars

Last number of vars.

Definition at line 228 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate(), and Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_op

TimeIntegrationSchemeOperators Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_op

private

Definition at line 243 of file TimeIntegrationAlgorithmGLM.h.

Referenced by InitializeData(), InitializeScheme(), and TimeIntegrate().

◆ m_order

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

Definition at line 172 of file TimeIntegrationAlgorithmGLM.h.

Referenced by Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_parent

const TimeIntegrationSchemeGLM* Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_parent {nullptr}

Parent scheme object.

Definition at line 168 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

◆ m_schemeType

TimeIntegrationSchemeType Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_schemeType {eNoTimeIntegrationSchemeType}

Definition at line 176 of file TimeIntegrationAlgorithmGLM.h.

Referenced by A(), B(), GetIntegrationSchemeType(), InitializeData(), TimeIntegrate(), U(), V(), and VerifyIntegrationSchemeType().

◆ m_T

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

private

Used to store the Explicit stage derivative of IMEX schemes.

Definition at line 252 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

◆ m_timeLevelOffset

Array<OneD, size_t> Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_timeLevelOffset

Definition at line 206 of file TimeIntegrationAlgorithmGLM.h.

Referenced by GetTimeLevelOffset().

◆ m_tmp

DoubleArray Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_tmp

private

Array containing the stage values.

Definition at line 246 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

◆ m_U

Array<TwoD, NekDouble> Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_U

Definition at line 209 of file TimeIntegrationAlgorithmGLM.h.

Referenced by CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), and U().

◆ m_U_phi

Array<OneD, Array<TwoD, NekDouble> > Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_U_phi

Definition at line 218 of file TimeIntegrationAlgorithmGLM.h.

Referenced by U(), and Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_V

Array<TwoD, NekDouble> Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_V

Definition at line 210 of file TimeIntegrationAlgorithmGLM.h.

Referenced by CheckIfLastStageEqualsNewSolution(), and V().

◆ m_V_phi

Array<OneD, Array<TwoD, NekDouble> > Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_V_phi

Definition at line 219 of file TimeIntegrationAlgorithmGLM.h.

Referenced by V(), and Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_variant

std::string Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_variant

Definition at line 171 of file TimeIntegrationAlgorithmGLM.h.

Referenced by Nektar::LibUtilities::EulerExponentialTimeIntegrationScheme::v_InitializeSecondaryData().

◆ m_Y

DoubleArray Nektar::LibUtilities::TimeIntegrationAlgorithmGLM::m_Y

private

Definition at line 245 of file TimeIntegrationAlgorithmGLM.h.

Referenced by TimeIntegrate().

Public Member Functions

Public Attributes

Private Member Functions

Private Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

◆ TimeIntegrationAlgorithmGLM()

◆ ~TimeIntegrationAlgorithmGLM()

Member Function Documentation

◆ A() [1/2]

◆ A() [2/2]

◆ A_IMEX()

◆ B() [1/2]

◆ B() [2/2]

◆ B_IMEX()

◆ CheckAndVerify()

◆ CheckIfFirstStageEqualsOldSolution()

◆ CheckIfLastStageEqualsNewSolution()

◆ CheckTimeIntegrateArguments()

◆ GetFirstDim()

◆ GetIntegrationSchemeType()

◆ GetNmultiStepExplicitDerivs()

◆ GetNmultiStepImplicitDerivs()

◆ GetNmultiStepValues()

◆ GetNstages()

◆ GetNsteps()

◆ GetSecondDim()

◆ GetTimeLevelOffset()

◆ InitializeData()

◆ InitializeScheme()

◆ TimeIntegrate() [1/2]

◆ TimeIntegrate() [2/2]

◆ U() [1/2]

◆ U() [2/2]

◆ V() [1/2]

◆ V() [2/2]

◆ VerifyIntegrationSchemeType()

Friends And Related Function Documentation

◆ operator<< [1/4]

◆ operator<< [2/4]

◆ operator<< [3/4]

◆ operator<< [4/4]

Member Data Documentation

◆ m_A

◆ m_A_phi

◆ m_B

◆ m_B_phi

◆ m_F

◆ m_F_IMEX

◆ m_firstStageEqualsOldSolution

◆ m_freeParams

◆ m_initialised

◆ m_L

◆ m_lastDeltaT

◆ m_lastNVars

◆ m_lastStageEqualsNewSolution

◆ m_name

◆ m_npoints

◆ m_numMultiStepExplicitDerivs

◆ m_numMultiStepImplicitDerivs

◆ m_numMultiStepValues

◆ m_numstages

◆ m_numsteps

◆ m_nvars

◆ m_op

◆ m_order

◆ m_parent

◆ m_schemeType

◆ m_T

◆ m_timeLevelOffset

◆ m_tmp

◆ m_U

◆ m_U_phi

◆ m_V

◆ m_V_phi

◆ m_variant

◆ m_Y