Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Types | Public Member Functions | Protected Attributes | Private Member Functions | Static Private Member Functions | Private Attributes | Friends | List of all members
Nektar::LibUtilities::TimeIntegrationScheme Class Reference

#include <TimeIntegrationScheme.h>

Collaboration diagram for Nektar::LibUtilities::TimeIntegrationScheme:
Collaboration graph
[legend]

Public Types

typedef const Array< OneD,
const Array< OneD, Array< OneD,
NekDouble > > > 		ConstTripleArray
typedef Array< OneD, Array
< OneD, Array< OneD, NekDouble > > > 		TripleArray
typedef const Array< OneD,
const Array< OneD, NekDouble > > 		ConstDoubleArray
typedef Array< OneD, Array
< OneD, NekDouble > > 		DoubleArray
typedef const Array< OneD,
const NekDouble > 		ConstSingleArray
typedef Array< OneD, NekDouble > SingleArray
typedef boost::function< void(ConstDoubleArray
&, DoubleArray &, const
NekDouble) > 		FunctorType1
typedef boost::function< void(ConstDoubleArray
&, DoubleArray &, const
NekDouble, const NekDouble) > 		FunctorType2

Public Member Functions

virtual ~TimeIntegrationScheme ()
const TimeIntegrationSchemeKey & GetIntegrationSchemeKey () const
TimeIntegrationMethod GetIntegrationMethod () const
TimeIntegrationSchemeType GetIntegrationSchemeType () const
NekDouble A (const unsigned int i, const unsigned int j) const
NekDouble B (const unsigned int i, const unsigned int j) const
NekDouble U (const unsigned int i, const unsigned int j) const
NekDouble V (const unsigned int i, const unsigned int j) const
NekDouble A_IMEX (const unsigned int i, const unsigned int j) const
NekDouble B_IMEX (const unsigned int i, const unsigned int j) const
unsigned int GetNsteps (void) const
unsigned int GetNstages (void) const
unsigned int GetNmultiStepValues (void) const
unsigned int GetNmultiStepDerivs (void) const
TimeIntegrationSolutionSharedPtr InitializeScheme (const NekDouble timestep, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)
 This function initialises the time integration scheme.
ConstDoubleArray & TimeIntegrate (const NekDouble timestep, TimeIntegrationSolutionSharedPtr &y, const TimeIntegrationSchemeOperators &op)
 Explicit integration of an ODE.
const Array< OneD, const
unsigned int > & 		GetTimeLevelOffset ()

Protected Attributes

TimeIntegrationSchemeKey m_schemeKey
TimeIntegrationSchemeType m_schemeType
unsigned int m_numsteps
unsigned int m_numstages
bool m_firstStageEqualsOldSolution
bool m_lastStageEqualsNewSolution
unsigned int m_numMultiStepValues
unsigned int m_numMultiStepDerivs
Array< OneD, unsigned int > 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

Private Member Functions

 TimeIntegrationScheme (const TimeIntegrationSchemeKey &key)
 TimeIntegrationScheme ()
 TimeIntegrationScheme (const TimeIntegrationScheme &in)
bool VerifyIntegrationSchemeType (TimeIntegrationSchemeType type, const Array< OneD, const Array< TwoD, NekDouble > > &A, const Array< OneD, const Array< TwoD, NekDouble > > &B, const Array< TwoD, const NekDouble > &U, const Array< TwoD, const NekDouble > &V) const
void TimeIntegrate (const NekDouble timestep, ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new, SingleArray &t_new, const TimeIntegrationSchemeOperators &op)
int GetFirstDim (ConstTripleArray &y) const
int GetSecondDim (ConstTripleArray &y) const
bool CheckTimeIntegrateArguments (const NekDouble timestep, ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new, SingleArray &t_new, const TimeIntegrationSchemeOperators &op) const
bool CheckIfFirstStageEqualsOldSolution (const Array< OneD, const Array< TwoD, NekDouble > > &A, const Array< OneD, const Array< TwoD, NekDouble > > &B, const Array< TwoD, const NekDouble > &U, const Array< TwoD, const NekDouble > &V) const
bool CheckIfLastStageEqualsNewSolution (const Array< OneD, const Array< TwoD, NekDouble > > &A, const Array< OneD, const Array< TwoD, NekDouble > > &B, const Array< TwoD, const NekDouble > &U, const Array< TwoD, const NekDouble > &V) const

Static Private Member Functions

static boost::shared_ptr
< TimeIntegrationScheme > 		Create (const TimeIntegrationSchemeKey &key)

Private Attributes

bool m_initialised
int m_nvar
 bool to identify if array has been initialised
int m_npoints
 The number of variables in integration scheme.
DoubleArray m_Y
 The size of inner data which is stored for reuse.
DoubleArray m_tmp
 Array containing the stage values.
TripleArray m_F
 explicit right hand side of each stage equation
TripleArray m_F_IMEX
 Array corresponding to the stage Derivatives.
NekDouble m_T
 Used to store the Explicit stage derivative of IMEX schemes.

Friends

class Nektar::MemoryManager
 Time at the different stages.
TimeIntegrationSchemeManagerT & TimeIntegrationSchemeManager (void)

Detailed Description

Definition at line 355 of file TimeIntegrationScheme.h.

Member Typedef Documentation

typedef const Array<OneD, const Array<OneD, NekDouble> > Nektar::LibUtilities::TimeIntegrationScheme::ConstDoubleArray

Definition at line 361 of file TimeIntegrationScheme.h.

typedef const Array<OneD, const NekDouble > Nektar::LibUtilities::TimeIntegrationScheme::ConstSingleArray

Definition at line 363 of file TimeIntegrationScheme.h.

typedef const Array<OneD, const Array<OneD, Array<OneD, NekDouble> > > Nektar::LibUtilities::TimeIntegrationScheme::ConstTripleArray

Definition at line 359 of file TimeIntegrationScheme.h.

typedef Array<OneD, Array<OneD, NekDouble> > Nektar::LibUtilities::TimeIntegrationScheme::DoubleArray

Definition at line 362 of file TimeIntegrationScheme.h.

typedef boost::function< void (ConstDoubleArray&, DoubleArray&, const NekDouble) > Nektar::LibUtilities::TimeIntegrationScheme::FunctorType1

Definition at line 365 of file TimeIntegrationScheme.h.

typedef boost::function< void (ConstDoubleArray&, DoubleArray&, const NekDouble, const NekDouble) > Nektar::LibUtilities::TimeIntegrationScheme::FunctorType2

Definition at line 366 of file TimeIntegrationScheme.h.

typedef Array<OneD, NekDouble > Nektar::LibUtilities::TimeIntegrationScheme::SingleArray

Definition at line 364 of file TimeIntegrationScheme.h.

typedef Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::LibUtilities::TimeIntegrationScheme::TripleArray

Definition at line 360 of file TimeIntegrationScheme.h.

Constructor & Destructor Documentation

virtual Nektar::LibUtilities::TimeIntegrationScheme::~TimeIntegrationScheme ( )
inlinevirtual

Definition at line 370 of file TimeIntegrationScheme.h.

{
}
Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrationScheme ( const TimeIntegrationSchemeKey &  key)
private

Definition at line 150 of file TimeIntegrationScheme.cpp.

References ASSERTL1, CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), Nektar::LibUtilities::eAdamsBashforthOrder1, Nektar::LibUtilities::eAdamsBashforthOrder2, Nektar::LibUtilities::eAdamsBashforthOrder3, Nektar::LibUtilities::eAdamsMoultonOrder1, Nektar::LibUtilities::eAdamsMoultonOrder2, Nektar::LibUtilities::eBackwardEuler, Nektar::LibUtilities::eBDFImplicitOrder1, Nektar::LibUtilities::eBDFImplicitOrder2, Nektar::LibUtilities::eClassicalRungeKutta4, Nektar::LibUtilities::eCNAB, Nektar::LibUtilities::eDiagonallyImplicit, Nektar::LibUtilities::eDIRKOrder2, Nektar::LibUtilities::eDIRKOrder3, Nektar::LibUtilities::eExplicit, ErrorUtil::efatal, Nektar::LibUtilities::eForwardEuler, Nektar::LibUtilities::eIMEX, Nektar::LibUtilities::eIMEXdirk_1_1_1, Nektar::LibUtilities::eIMEXdirk_1_2_1, Nektar::LibUtilities::eIMEXdirk_1_2_2, Nektar::LibUtilities::eIMEXdirk_2_2_2, Nektar::LibUtilities::eIMEXdirk_2_3_2, Nektar::LibUtilities::eIMEXdirk_2_3_3, Nektar::LibUtilities::eIMEXdirk_3_4_3, Nektar::LibUtilities::eIMEXdirk_4_4_3, Nektar::LibUtilities::eIMEXGear, Nektar::LibUtilities::eIMEXOrder1, Nektar::LibUtilities::eIMEXOrder2, Nektar::LibUtilities::eIMEXOrder3, Nektar::LibUtilities::eMCNAB, Nektar::LibUtilities::eMidpoint, Nektar::LibUtilities::eRungeKutta2_ImprovedEuler, Nektar::LibUtilities::eRungeKutta2_ModifiedEuler, Nektar::LibUtilities::TimeIntegrationSchemeKey::GetIntegrationMethod(), m_A, m_B, m_firstStageEqualsOldSolution, m_lastStageEqualsNewSolution, m_numMultiStepDerivs, m_numMultiStepValues, m_numstages, m_numsteps, m_schemeType, m_timeLevelOffset, m_U, m_V, NEKERROR, and VerifyIntegrationSchemeType().

:
m_schemeKey(key),
m_initialised(false)
{
switch(key.GetIntegrationMethod())
{
case eForwardEuler:
case eAdamsBashforthOrder1:
{
m_numsteps = 1;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 1.0);
m_schemeType = eExplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eAdamsBashforthOrder2:
{
m_numsteps = 2;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
m_B[0][0][0] = 3.0/2.0;
m_B[0][1][0] = 1.0;
m_U[0][0] = 1.0;
m_V[0][0] = 1.0;
m_V[0][1] = -0.5;
m_schemeType = eExplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 1;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 1;
}
break;
case eAdamsBashforthOrder3:
{
m_numsteps = 4;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
m_B[0][1][0] = 1.0;
m_U[0][0] = 1.0;
m_U[0][1] = 23.0/12.0;
m_U[0][2] = -4.0/3.0;
m_U[0][3] = 5.0/12.0;
m_V[0][0] = 1.0;
m_V[0][1] = 23.0/12.0;
m_V[0][2] = -4.0/3.0;
m_V[0][3] = 5.0/12.0;
m_V[2][1] = 1.0;
m_V[3][2] = 1.0;
m_schemeType = eExplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 3;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 1;
m_timeLevelOffset[2] = 2;
m_timeLevelOffset[3] = 3;
}
break;
case eBackwardEuler:
case eBDFImplicitOrder1:
case eAdamsMoultonOrder1:
{
m_numsteps = 1;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,1.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 1.0);
m_schemeType = eDiagonallyImplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0; // SJS: Not sure whether this is correct
}
break;
case eIMEXOrder1:
{
m_numsteps = 2;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,1.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
m_B[0][1][0] = 0.0;
m_B[1][0][0] = 0.0;
m_V[0][0] = 1.0;
m_V[0][1] = 1.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 1;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 0;
}
break;
case eIMEXOrder2:
{
NekDouble third = 1.0/3.0;
m_numsteps = 4;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,2*third);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 4*third);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
m_B[0][0][0] = 2*third;
m_B[1][2][0] = 1.0;
m_U[0][1] = -third;
m_U[0][3] = -2*third;
m_V[0][0] = 4*third;
m_V[0][1] = -third;
m_V[0][2] = 4*third;
m_V[0][3] = -2*third;
m_V[1][0] = 1.0;
m_V[3][2] = 1.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 2;
m_numMultiStepDerivs = 2;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 1;
m_timeLevelOffset[2] = 0;
m_timeLevelOffset[3] = 1;
}
break;
case eIMEXOrder3:
{
NekDouble eleventh = 1.0/11.0;
m_numsteps = 6;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,6*eleventh);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 18*eleventh);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
m_B[0][0][0] = 6*eleventh;
m_B[1][3][0] = 1.0;
m_U[0][1] = -9*eleventh;
m_U[0][2] = 2*eleventh;
m_U[0][4] = -18*eleventh;
m_U[0][5] = 6*eleventh;
m_V[0][0] = 18*eleventh;
m_V[0][1] = -9*eleventh;
m_V[0][2] = 2*eleventh;
m_V[0][3] = 18*eleventh;
m_V[0][4] = -18*eleventh;
m_V[0][5] = 6*eleventh;
m_V[1][0] = 1.0;
m_V[2][1] = 1.0;
m_V[4][3] = 1.0;
m_V[5][4] = 1.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 3;
m_numMultiStepDerivs = 3;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 1;
m_timeLevelOffset[2] = 2;
m_timeLevelOffset[3] = 0;
m_timeLevelOffset[4] = 1;
m_timeLevelOffset[5] = 2;
}
break;
case eAdamsMoultonOrder2:
{
m_numsteps = 2;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.5);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.0);
m_B[0][0][0] = 0.5;
m_B[0][1][0] = 1.0;
m_U[0][0] = 1.0;
m_U[0][1] = 0.5;
m_V[0][0] = 1.0;
m_V[0][1] = 0.5;
m_schemeType = eDiagonallyImplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 1;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 0;
}
break;
case eBDFImplicitOrder2:
{
NekDouble third = 1.0/3.0;
m_numsteps = 2;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,2*third);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.0);
m_B[0][0][0] = 2*third;
m_B[0][1][0] = 0.0;
m_U[0][0] = 4*third;
m_U[0][1] = -third;
m_V[0][0] = 4*third;
m_V[0][1] = -third;
m_V[1][0] = 1.0;
m_schemeType = eDiagonallyImplicit;
m_numMultiStepValues = 2;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 1;
}
break;
case eMidpoint:
{
m_numsteps = 1;
m_numstages = 2;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
m_A[0][1][0] = 0.5;
m_B[0][0][1] = 1.0;
m_schemeType = eExplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eClassicalRungeKutta4:
{
m_numsteps = 1;
m_numstages = 4;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
m_A[0][1][0] = 0.5;
m_A[0][2][1] = 0.5;
m_A[0][3][2] = 1.0;
m_B[0][0][0] = 1.0/6.0;
m_B[0][0][1] = 1.0/3.0;
m_B[0][0][2] = 1.0/3.0;
m_B[0][0][3] = 1.0/6.0;
m_schemeType = eExplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eRungeKutta2_ModifiedEuler:
{
m_numsteps = 1;
m_numstages = 2;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.5);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.5);
m_A[0][1][0] = 0.5;
m_B[0][0][1] = 1.0;
m_schemeType = eExplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eRungeKutta2_ImprovedEuler:
{
m_numsteps = 1;
m_numstages = 2;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.5);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
m_A[0][1][0] = 1.0;
m_schemeType = eExplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eDIRKOrder2:
{
m_numsteps = 1;
m_numstages = 2;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
NekDouble lambda = (2.0-sqrt(2.0))/2.0;
m_A[0][0][0] = lambda;
m_A[0][1][0] = 1.0 - lambda;
m_A[0][1][1] = lambda;
m_B[0][0][0] = 1.0 - lambda;
m_B[0][0][1] = lambda;
m_schemeType = eDiagonallyImplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eDIRKOrder3:
{
m_numsteps = 1;
m_numstages = 3;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
NekDouble lambda = 0.4358665215;
m_A[0][0][0] = lambda;
m_A[0][1][0] = 0.5 * (1.0 - lambda);
m_A[0][2][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
m_A[0][1][1] = lambda;
m_A[0][2][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
m_A[0][2][2] = lambda;
m_B[0][0][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
m_B[0][0][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
m_B[0][0][2] = lambda;
m_schemeType = eDiagonallyImplicit;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eIMEXdirk_2_3_2:
{
m_numsteps = 1;
m_numstages = 3;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
NekDouble lambda = (2.0-sqrt(2.0))/2.0;
NekDouble delta = -2.0*sqrt(2.0)/3.0;
m_A[0][1][1] = lambda;
m_A[0][2][1] = 1.0 - lambda;
m_A[0][2][2] = lambda;
m_B[0][0][1] = 1.0 - lambda;
m_B[0][0][2] = lambda;
m_A[1][1][0] = lambda;
m_A[1][2][0] = delta;
m_A[1][2][1] = 1.0 - delta;
m_B[1][0][1] = 1.0 - lambda;
m_B[1][0][2] = lambda;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eIMEXdirk_3_4_3:
{
m_numsteps = 1;
m_numstages = 4;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
NekDouble lambda = 0.4358665215;
m_A[0][1][1] = lambda;
m_A[0][2][1] = 0.5 * (1.0 - lambda);
m_A[0][3][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
m_A[0][2][2] = lambda;
m_A[0][3][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
m_A[0][3][3] = lambda;
m_B[0][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
m_B[0][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
m_B[0][0][3] = lambda;
m_A[1][1][0] = 0.4358665215;
m_A[1][2][0] = 0.3212788860;
m_A[1][2][1] = 0.3966543747;
m_A[1][3][0] =-0.105858296;
m_A[1][3][1] = 0.5529291479;
m_A[1][3][2] = 0.5529291479;
m_B[1][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
m_B[1][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
m_B[1][0][3] = lambda;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eCNAB:
{
NekDouble secondth = 1.0/2.0;
m_numsteps = 4;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,secondth);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,secondth);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
m_B[0][0][0] = secondth;
m_B[0][1][0] = 1.0;
m_B[1][2][0] = 1.0;
m_U[0][0] = 2*secondth;
m_U[0][2] = 3*secondth;
m_U[0][3] = -1*secondth;
m_V[0][0] = 2*secondth;
m_V[0][1] = secondth;
m_V[0][2] = 3*secondth;
m_V[0][3] = -1*secondth;
m_V[3][2] = 1.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 3;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 0;
m_timeLevelOffset[2] = 0;
m_timeLevelOffset[3] = 1;
}
break;
case eIMEXGear:
{
NekDouble twothirdth = 2.0/3.0;
m_numsteps = 2;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,twothirdth);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
m_B[0][0][0] = twothirdth;
m_B[1][0][0] = twothirdth;
m_U[0][0] = 2*twothirdth;
m_U[0][1] = -0.5*twothirdth;
m_V[0][0] = 2*twothirdth;
m_V[0][1] = -0.5*twothirdth;
m_V[1][0] = 1.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 2;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 1;
}
break;
case eMCNAB:
{
NekDouble sixthx = 9.0/16.0;
m_numsteps = 5;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,sixthx);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
m_B[0][0][0] = sixthx;
m_B[0][1][0] = 1.0;
m_B[1][3][0] = 1.0;
m_U[0][0] = 1.0;
m_U[0][1] = 6.0/16.0;
m_U[0][2] = 1.0/16.0;
m_U[0][3] = 1.5;
m_U[0][4] = -0.5;
m_V[0][0] = 1.0;
m_V[0][1] = 6.0/16.0;
m_V[0][2] = 1.0/16.0;
m_V[0][3] = 1.5;
m_V[0][4] = -0.5;
m_V[2][1] = 1.0;
m_V[4][3] = 1.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 4;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
m_timeLevelOffset[1] = 0;
m_timeLevelOffset[2] = 1;
m_timeLevelOffset[3] = 0;
m_timeLevelOffset[4] = 1;
}
break;
case eIMEXdirk_2_2_2:
{
m_numsteps = 1;
m_numstages = 3;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
NekDouble glambda = 0.788675134594813;
NekDouble gdelta = 0.366025403784439;
m_A[0][1][1] = glambda;
m_A[0][2][1] = 1.0 - glambda;
m_A[0][2][2] = glambda;
m_B[0][0][1] = 1.0 - glambda;
m_B[0][0][2] = glambda;
m_A[1][1][0] = glambda;
m_A[1][2][0] = gdelta;
m_A[1][2][1] = 1.0 - gdelta;
m_B[1][0][0] = gdelta;
m_B[1][0][1] = 1.0 - gdelta;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eIMEXdirk_2_3_3:
{
m_numsteps = 1;
m_numstages = 3;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
NekDouble glambda = 0.788675134594813;
m_A[0][1][1] = glambda;
m_A[0][2][1] = 1.0 - 2.0*glambda;
m_A[0][2][2] = glambda;
m_B[0][0][1] = 0.5;
m_B[0][0][2] = 0.5;
m_A[1][1][0] = glambda;
m_A[1][2][0] = glambda - 1.0;
m_A[1][2][1] = 2.0*(1-glambda);
m_B[1][0][1] = 0.5;
m_B[1][0][2] = 0.5;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eIMEXdirk_1_1_1:
{
m_numsteps = 1;
m_numstages = 2;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
m_A[0][1][1] = 1.0;
m_B[0][0][1] = 1.0;
m_A[1][1][0] = 1.0;
m_B[1][0][0] = 1.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eIMEXdirk_1_2_1:
{
m_numsteps = 1;
m_numstages = 2;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
m_A[0][1][1] = 1.0;
m_B[0][0][1] = 1.0;
m_A[1][1][0] = 1.0;
m_B[1][0][1] = 1.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eIMEXdirk_1_2_2:
{
m_numsteps = 1;
m_numstages = 2;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
m_A[0][1][1] = 1.0/2.0;
m_B[0][0][1] = 1.0;
m_A[1][1][0] = 1.0/2.0;
m_B[1][0][1] = 1.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
case eIMEXdirk_4_4_3:
{
m_numsteps = 1;
m_numstages = 5;
m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
m_A[0][1][1] = 1.0/2.0;
m_A[0][2][1] = 1.0/6.0;
m_A[0][2][2] = 1.0/2.0;
m_A[0][3][1] = -1.0/2.0;
m_A[0][3][2] = 1.0/2.0;
m_A[0][3][3] = 1.0/2.0;
m_A[0][4][1] = 3.0/2.0;
m_A[0][4][2] = -3.0/2.0;
m_A[0][4][3] = 1.0/2.0;
m_A[0][4][4] = 1.0/2.0;
m_B[0][0][1] = 3.0/2.0;
m_B[0][0][2] = -3.0/2.0;
m_B[0][0][3] = 1.0/2.0;
m_B[0][0][4] = 1.0/2.0;
m_A[1][1][0] = 1.0/2.0;
m_A[1][2][0] = 11.0/18.0;
m_A[1][2][1] = 1.0/18.0;
m_A[1][3][0] = 5.0/6.0;
m_A[1][3][1] = -5.0/6.0;
m_A[1][3][2] = 1.0/2.0;
m_A[1][4][0] = 1.0/4.0;
m_A[1][4][1] = 7.0/4.0;
m_A[1][4][2] = 3.0/4.0;
m_A[1][4][3] = -7.0/4.0;
m_B[1][0][0] = 1.0/4.0;
m_B[1][0][1] = 7.0/4.0;
m_B[1][0][2] = 3.0/4.0;
m_B[1][0][3] = -7.0/4.0;
m_schemeType = eIMEX;
m_numMultiStepValues = 1;
m_numMultiStepDerivs = 0;
m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
m_timeLevelOffset[0] = 0;
}
break;
default:
{
NEKERROR(ErrorUtil::efatal,"Invalid Time Integration Scheme");
}
}
m_firstStageEqualsOldSolution = CheckIfFirstStageEqualsOldSolution(m_A,m_B,m_U,m_V);
m_lastStageEqualsNewSolution = CheckIfLastStageEqualsNewSolution(m_A,m_B,m_U,m_V);
ASSERTL1(VerifyIntegrationSchemeType(m_schemeType,m_A,m_B,m_U,m_V),
"Time integration scheme coefficients do not match its type");
}
Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrationScheme ( )
inlineprivate

Definition at line 554 of file TimeIntegrationScheme.h.

References ErrorUtil::efatal, and NEKERROR.

:m_schemeKey(NullTimeIntegrationSchemeKey)
{
NEKERROR(ErrorUtil::efatal,"Default Constructor for the TimeIntegrationScheme class should not be called");
}
Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrationScheme ( const TimeIntegrationScheme &  in)
inlineprivate

Definition at line 559 of file TimeIntegrationScheme.h.

References ErrorUtil::efatal, and NEKERROR.

:m_schemeKey(NullTimeIntegrationSchemeKey)
{
NEKERROR(ErrorUtil::efatal,"Copy Constructor for the TimeIntegrationScheme class should not be called");
}

Member Function Documentation

NekDouble Nektar::LibUtilities::TimeIntegrationScheme::A ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 389 of file TimeIntegrationScheme.h.

References m_A.

Referenced by CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), Nektar::LibUtilities::operator<<(), TimeIntegrate(), and VerifyIntegrationSchemeType().

{
return m_A[0][i][j];
}
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::A_IMEX ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 409 of file TimeIntegrationScheme.h.

References m_A.

Referenced by Nektar::LibUtilities::operator<<(), and TimeIntegrate().

{
return m_A[1][i][j];
}
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::B ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 394 of file TimeIntegrationScheme.h.

References m_B.

Referenced by CheckIfLastStageEqualsNewSolution(), Nektar::LibUtilities::operator<<(), TimeIntegrate(), and VerifyIntegrationSchemeType().

{
return m_B[0][i][j];
}
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::B_IMEX ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 414 of file TimeIntegrationScheme.h.

References m_B.

Referenced by Nektar::LibUtilities::operator<<(), and TimeIntegrate().

{
return m_B[1][i][j];
}
bool Nektar::LibUtilities::TimeIntegrationScheme::CheckIfFirstStageEqualsOldSolution ( const Array< OneD, const Array< TwoD, NekDouble > > &  A,
const Array< OneD, const Array< TwoD, NekDouble > > &  B,
const Array< TwoD, const NekDouble > &  U,
const Array< TwoD, const NekDouble > &  V 
) const
private

Definition at line 1679 of file TimeIntegrationScheme.cpp.

References A(), Nektar::NekConstants::kNekZeroTol, m_numstages, and m_numsteps.

Referenced by TimeIntegrationScheme().

{
int i,m;
// 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(m = 0; m < A.num_elements(); m++)
{
for(i = 0; i < m_numstages; i++)
{
if( fabs(A[m][0][i]) > NekConstants::kNekZeroTol )
{
return false;
}
}
}
// 2. Check second condition
if( fabs(U[0][0] - 1.0) > NekConstants::kNekZeroTol )
{
return false;
}
for(i = 1; i < m_numsteps; i++)
{
if( fabs(U[0][i]) > NekConstants::kNekZeroTol )
{
return false;
}
}
return true;
}
bool Nektar::LibUtilities::TimeIntegrationScheme::CheckIfLastStageEqualsNewSolution ( const Array< OneD, const Array< TwoD, NekDouble > > &  A,
const Array< OneD, const Array< TwoD, NekDouble > > &  B,
const Array< TwoD, const NekDouble > &  U,
const Array< TwoD, const NekDouble > &  V 
) const
private

Definition at line 1717 of file TimeIntegrationScheme.cpp.

References A(), B(), Nektar::NekConstants::kNekZeroTol, m_numstages, and m_numsteps.

Referenced by TimeIntegrationScheme().

{
int i,m;
// 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(m = 0; m < A.num_elements(); m++)
{
for(i = 0; i < m_numstages; i++)
{
if( fabs(A[m][m_numstages-1][i]-B[m][0][i]) > NekConstants::kNekZeroTol )
{
return false;
}
}
}
// 2. Check second condition
for(i = 0; i < m_numsteps; i++)
{
if( fabs(U[m_numstages-1][i]-V[0][i]) > NekConstants::kNekZeroTol )
{
return false;
}
}
return true;
}
bool Nektar::LibUtilities::TimeIntegrationScheme::CheckTimeIntegrateArguments ( const NekDouble  timestep,
ConstTripleArray &  y_old,
ConstSingleArray &  t_old,
TripleArray &  y_new,
SingleArray &  t_new,
const TimeIntegrationSchemeOperators &  op 
) const
private

Definition at line 1751 of file TimeIntegrationScheme.cpp.

References ASSERTL1, and m_numsteps.

Referenced by TimeIntegrate().

{
// Check if arrays are all of consistent size
ASSERTL1(y_old.num_elements()==m_numsteps,"Non-matching number of steps.");
ASSERTL1(y_new.num_elements()==m_numsteps,"Non-matching number of steps.");
ASSERTL1(y_old[0]. num_elements()==y_new[0]. num_elements(),"Non-matching number of variables.");
ASSERTL1(y_old[0][0].num_elements()==y_new[0][0].num_elements(),"Non-matching number of coefficients.");
ASSERTL1(t_old.num_elements()==m_numsteps,"Non-matching number of steps.");
ASSERTL1(t_new.num_elements()==m_numsteps,"Non-matching number of steps.");
return true;
}
TimeIntegrationSchemeSharedPtr Nektar::LibUtilities::TimeIntegrationScheme::Create ( const TimeIntegrationSchemeKey &  key)
staticprivate

Definition at line 143 of file TimeIntegrationScheme.cpp.

Referenced by Nektar::LibUtilities::TimeIntegrationSchemeManager().

{
TimeIntegrationSchemeSharedPtr returnval(
MemoryManager<TimeIntegrationScheme>::AllocateSharedPtr(key));
return returnval;
}
int Nektar::LibUtilities::TimeIntegrationScheme::GetFirstDim ( ConstTripleArray &  y) const
inlineprivate

Definition at line 578 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

{
return y[0].num_elements();
}
TimeIntegrationMethod Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationMethod ( ) const
inline

Definition at line 379 of file TimeIntegrationScheme.h.

References Nektar::LibUtilities::TimeIntegrationSchemeKey::GetIntegrationMethod(), and m_schemeKey.

Referenced by Nektar::LibUtilities::operator<<().

{
return m_schemeKey.GetIntegrationMethod();
}
const TimeIntegrationSchemeKey& Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeKey ( ) const
inline

Definition at line 374 of file TimeIntegrationScheme.h.

References m_schemeKey.

Referenced by TimeIntegrate().

{
return m_schemeKey;
}
TimeIntegrationSchemeType Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeType ( ) const
inline

Definition at line 384 of file TimeIntegrationScheme.h.

References m_schemeType.

Referenced by Nektar::LibUtilities::operator<<(), and TimeIntegrate().

{
return m_schemeType;
}
unsigned int Nektar::LibUtilities::TimeIntegrationScheme::GetNmultiStepDerivs ( void  ) const
inline

Definition at line 434 of file TimeIntegrationScheme.h.

References m_numMultiStepDerivs.

{
return m_numMultiStepDerivs;
}
unsigned int Nektar::LibUtilities::TimeIntegrationScheme::GetNmultiStepValues ( void  ) const
inline

Definition at line 429 of file TimeIntegrationScheme.h.

References m_numMultiStepValues.

{
return m_numMultiStepValues;
}
unsigned int Nektar::LibUtilities::TimeIntegrationScheme::GetNstages ( void  ) const
inline

Definition at line 424 of file TimeIntegrationScheme.h.

References m_numstages.

Referenced by Nektar::LibUtilities::operator<<().

{
return m_numstages;
}
unsigned int Nektar::LibUtilities::TimeIntegrationScheme::GetNsteps ( void  ) const
inline

Definition at line 419 of file TimeIntegrationScheme.h.

References m_numsteps.

Referenced by Nektar::LibUtilities::operator<<().

{
return m_numsteps;
}
int Nektar::LibUtilities::TimeIntegrationScheme::GetSecondDim ( ConstTripleArray &  y) const
inlineprivate

Definition at line 583 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

{
return y[0][0].num_elements();
}
const Array<OneD, const unsigned int>& Nektar::LibUtilities::TimeIntegrationScheme::GetTimeLevelOffset ( )
inline

Definition at line 499 of file TimeIntegrationScheme.h.

References m_timeLevelOffset.

{
return m_timeLevelOffset;
}
TimeIntegrationSolutionSharedPtr Nektar::LibUtilities::TimeIntegrationScheme::InitializeScheme ( const NekDouble  timestep,
ConstDoubleArray &  y_0,
const NekDouble  time,
const TimeIntegrationSchemeOperators &  op 
)

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 TimeIntegrationSolution. 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 TimeIntegrationSolution 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
timestepThe size of the timestep, i.e. $\Delta t$.
timeon input: the initial time, i.e. $t_0$.
timeon 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.
nstepson 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.
fan 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_0the initial value $\boldsymbol{y}(t_0)$
Returns
An object of the class TimeIntegrationSolution which represents the vectors $\boldsymbol{y}^{[n]}$ and $t^{[n]}$ that can be used to start the actual integration.

Definition at line 1095 of file TimeIntegrationScheme.cpp.

References Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs(), m_numMultiStepDerivs, m_numMultiStepValues, m_schemeKey, and m_timeLevelOffset.

{
// create a TimeIntegrationSolution object based upon the
// initial value. Initialise all other multi-step values
// and derivatives to zero
TimeIntegrationSolutionSharedPtr y_out =
MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(m_schemeKey,y_0,time,timestep);
// calculate the initial derivative, if is part of the
// solution vector of the current scheme
if(m_numMultiStepDerivs)
{
if(m_timeLevelOffset[m_numMultiStepValues] == 0)
{
int i;
int nvar = y_0.num_elements();
int npoints = y_0[0].num_elements();
DoubleArray f_y_0(nvar);
for(i = 0; i < nvar; i++)
{
f_y_0[i] = Array<OneD,NekDouble>(npoints);
}
// calculate the derivative of the initial value
op.DoOdeRhs(y_0,f_y_0,time);
// multiply by the step size
for(i = 0; i < nvar; i++)
{
Blas::Dscal(npoints,timestep,f_y_0[i].get(),1);
}
y_out->SetDerivative(0,f_y_0,timestep);
}
}
return y_out;
}
TimeIntegrationScheme::ConstDoubleArray & Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrate ( const NekDouble  timestep,
TimeIntegrationSolutionSharedPtr &  y,
const TimeIntegrationSchemeOperators &  op 
)

Explicit integration of an ODE.

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

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

Parameters
timestepThe size of the timestep, i.e. $\Delta t$.
fan 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.
yon input: the vectors $\boldsymbol{y}^{[n-1]}$ and $t^{[n-1]}$ (which corresponds to the solution at the old time level)
yon 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 1136 of file TimeIntegrationScheme.cpp.

References ASSERTL1, Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs(), Nektar::LibUtilities::eImplicit, GetIntegrationSchemeKey(), GetIntegrationSchemeType(), m_numMultiStepDerivs, m_numMultiStepValues, m_numsteps, m_timeLevelOffset, and Vmath::Smul().

{
ASSERTL1(!(GetIntegrationSchemeType() == eImplicit),
"Fully Implicit integration scheme cannot be handled by this routine.");
int nvar = solvector->GetFirstDim ();
int npoints = solvector->GetSecondDim();
if( (solvector->GetIntegrationScheme()).get() != 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
int n;
DoubleArray y_n;
NekDouble t_n;
DoubleArray dtFy_n;
unsigned int nCurSchemeVals = m_numMultiStepValues; // number of required values of the current scheme
unsigned int nCurSchemeDers = m_numMultiStepDerivs; // number of required derivs of the current scheme
unsigned int nCurSchemeSteps = m_numsteps; // number of steps in the current scheme
unsigned int nMasterSchemeVals = solvector->GetNvalues(); // number of values of the master scheme
unsigned int nMasterSchemeDers = solvector->GetNderivs(); // number of derivs of the master scheme
// The arrays below contains information to which
// time-level the values and derivatives of the
// schemes belong
const Array<OneD, const unsigned int>& curTimeLevels = m_timeLevelOffset;
const Array<OneD, const unsigned int>& masterTimeLevels = solvector->GetTimeLevelOffset();
// 1.2 Copy the required information from the master
// solution vector to the input solution vector of
// the current scheme
TimeIntegrationSolutionSharedPtr solvector_in = MemoryManager<TimeIntegrationSolution>::
AllocateSharedPtr(GetIntegrationSchemeKey()); // input solution vector of the current scheme
for(n = 0; n < nCurSchemeVals; n++)
{
// Get the required value out of the master solution vector
//DoubleArray& y_n = solvector->GetValue ( curTimeLevels[n] );
//NekDouble t_n = solvector->GetValueTime( curTimeLevels[n] );
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);
}
for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
{
// Get the required derivative out of the master
// solution vector
//DoubleArray& dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
// Set the required derivative in the input
// solution vector of the current scheme
solvector_in->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
}
// STEP 2: time-integrate for one step using the
// current scheme
TimeIntegrationSolutionSharedPtr solvector_out = MemoryManager<TimeIntegrationSolution>:: AllocateSharedPtr(GetIntegrationSchemeKey(),nvar,npoints); // output solution vector of the current scheme
// integrate
TimeIntegrate(timestep, solvector_in->GetSolutionVector(),
solvector_in->GetTimeVector(),
solvector_out->UpdateSolutionVector(),
solvector_out->UpdateTimeVector(),op);
// 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.
int j;
bool CalcNewDeriv = false; // flag inidicating whether the new derivative is availble in the output of
// of the current scheme or whether it should be calculated
if( nMasterSchemeDers > 0 )
{
if(nCurSchemeDers == 0)
{
CalcNewDeriv = true;
}
else
{
if( masterTimeLevels[nMasterSchemeVals] < curTimeLevels[nCurSchemeVals] )
{
CalcNewDeriv = true;
}
}
}
if(CalcNewDeriv)
{
int newDerivTimeLevel = masterTimeLevels[nMasterSchemeVals]; // contains the time level at which
// we want to know the derivative of the
// master scheme
//DoubleArray y_n;
//NekDouble t_n;
// if the time level correspond to 0, calculate the derivative based upon the solution value
// at the new time-level
if (newDerivTimeLevel == 0)
{
y_n = solvector_out->GetValue(0);
t_n = solvector_out->GetValueTime(0);
}
// if the time level correspond to 1, calculate the derivative based upon the solution value
// at the new old-level
else if( newDerivTimeLevel == 1 )
{
y_n = solvector->GetValue(0);
t_n = solvector->GetValueTime(0);
}
else
{
ASSERTL1(false,"Problems with initialising scheme");
}
DoubleArray f_n(nvar);
for(j = 0; j < nvar; j++)
{
f_n[j] = Array<OneD, NekDouble>(npoints);
}
// calculate the derivative
op.DoOdeRhs(y_n, f_n, t_n);
// multiply by dt (as required by the General Linear Method framework)
for(j = 0; j < nvar; j++)
{
Vmath::Smul(npoints,timestep,f_n[j],1,
f_n[j],1);
}
// Rotate the solution vector
// (i.e. updating without calculating/inserting new values)
solvector->RotateSolutionVector();
// Set the calculated derivative in the master solution vector
solvector->SetDerivative(newDerivTimeLevel,f_n,timestep);
}
else
{
// Rotate the solution vector (i.e. updating
// without calculating/inserting new values)
solvector->RotateSolutionVector();
}
// 1.2 Copy the information calculated using the
// current scheme from the output solution vector
// to the master solution vector
for(n = 0; n < nCurSchemeVals; n++)
{
// Get the calculated value out of the output
// solution vector of the current scheme
//DoubleArray& y_n = solvector_out->GetValue ( curTimeLevels[n] );
//NekDouble t_n = solvector_out->GetValueTime( curTimeLevels[n] );
y_n = solvector_out->GetValue ( curTimeLevels[n] );
t_n = solvector_out->GetValueTime( curTimeLevels[n] );
// Set the calculated value in the master solution vector
solvector->SetValue(curTimeLevels[n],y_n,t_n);
}
for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
{
// Get the calculated derivative out of the output
// solution vector of the current scheme
// DoubleArray& dtFy_n =
// solvector_out->GetDerivative (curTimeLevels[n]);
dtFy_n = solvector_out->GetDerivative ( curTimeLevels[n] );
// Set the calculated derivative in the master
// solution vector
solvector->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
}
}
else
{
const TimeIntegrationSchemeKey& key = solvector->GetIntegrationSchemeKey();
TimeIntegrationSolutionSharedPtr solvector_new = MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(key,nvar,npoints);
TimeIntegrate(timestep,solvector->GetSolutionVector(),
solvector->GetTimeVector(),
solvector_new->UpdateSolutionVector(),
solvector_new->UpdateTimeVector(),op);
solvector = solvector_new;
}
return solvector->GetSolution();
}
void Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrate ( const NekDouble  timestep,
ConstTripleArray &  y_old,
ConstSingleArray &  t_old,
TripleArray &  y_new,
SingleArray &  t_new,
const TimeIntegrationSchemeOperators &  op 
)
private

Definition at line 1361 of file TimeIntegrationScheme.cpp.

References A(), A_IMEX(), ASSERTL1, B(), B_IMEX(), CheckTimeIntegrateArguments(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoImplicitSolve(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoProjection(), Nektar::LibUtilities::eDiagonallyImplicit, Nektar::LibUtilities::eExplicit, Nektar::LibUtilities::eIMEX, GetFirstDim(), GetIntegrationSchemeType(), GetSecondDim(), Nektar::NekConstants::kNekZeroTol, m_F, m_F_IMEX, m_firstStageEqualsOldSolution, m_initialised, m_lastStageEqualsNewSolution, m_npoints, m_numstages, m_numsteps, m_nvar, m_T, m_tmp, m_Y, Vmath::Smul(), Vmath::Svtvp(), U(), V(), Vmath::Vcopy(), Vmath::Vsub(), and Vmath::Zero().

{
ASSERTL1(CheckTimeIntegrateArguments(timestep,y_old,t_old,y_new,t_new,op), "Arguments not well defined");
unsigned int i,j,k;
TimeIntegrationSchemeType type = GetIntegrationSchemeType();
// Check if storage has already been initialised.
// If so, we just zero the temporary storage.
if (m_initialised && m_nvar == GetFirstDim(y_old)
&& m_npoints == GetSecondDim(y_old))
{
for(j = 0; j < m_nvar; j++)
{
Vmath::Zero(m_npoints, m_tmp[j], 1);
}
}
else
{
m_nvar = GetFirstDim(y_old);
m_npoints = GetSecondDim(y_old);
// First, we are going to 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_nvar);
for(j = 0; j < m_nvar; 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)
{
m_Y = m_tmp;
}
else
{
m_Y = DoubleArray(m_nvar);
for(j = 0; j < m_nvar; 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(i = 0; i < m_numstages; ++i)
{
m_F[i] = DoubleArray(m_nvar);
for(j = 0; j < m_nvar; j++)
{
m_F[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
}
}
if(type == eIMEX)
{
m_F_IMEX = TripleArray(m_numstages);
for(i = 0; i < m_numstages; ++i)
{
m_F_IMEX[i] = DoubleArray(m_nvar);
for(j = 0; j < m_nvar; j++)
{
m_F_IMEX[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
}
}
}
// Finally, flag that we have initialised the memory.
m_initialised = true;
}
// The loop below calculates the stage values and derivatives
for(i = 0; i < m_numstages; i++)
{
if( (i==0) && m_firstStageEqualsOldSolution )
{
for(k = 0; k < m_nvar; k++)
{
Vmath::Vcopy(m_npoints,y_old[0][k],1,m_Y[k],1);
}
m_T = t_old[0];
}
else
{
// The stage values m_Y are a linear combination of:
// 1: the stage derivatives
if( i != 0 )
{
for(k = 0; k < m_nvar; k++)
{
Vmath::Smul(m_npoints,timestep*A(i,0),m_F[0][k],1,
m_tmp[k],1);
if(type == eIMEX)
{
Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,0),
m_F_IMEX[0][k],1,
m_tmp[k],1,m_tmp[k],1);
}
}
}
m_T = A(i,0)*timestep;
for( j = 1; j < i; j++ )
{
for(k = 0; k < m_nvar; k++)
{
Vmath::Svtvp(m_npoints,timestep*A(i,j),m_F[j][k],1,
m_tmp[k],1,m_tmp[k],1);
if(type == eIMEX)
{
Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,j),
m_F_IMEX[j][k],1,
m_tmp[k],1,m_tmp[k],1);
}
}
m_T += A(i,j)*timestep;
}
// 2: the imported multi-step solution of the
// previous time level
for(j = 0; j < m_numsteps; j++)
{
for(k = 0; k < m_nvar; k++)
{
Vmath::Svtvp(m_npoints,U(i,j),y_old[j][k],1,
m_tmp[k],1,m_tmp[k],1);
}
m_T += U(i,j)*t_old[j];
}
}
// Calculate the stage derivative based upon the stage value
if(type == eDiagonallyImplicit)
{
if(m_numstages==1)
{
m_T= t_old[0]+timestep;
}
else
{
m_T= t_old[0];
for(int j=0; j<=i; ++j)
{
m_T += A(i,j)*timestep;
}
}
op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
for(k = 0; k < m_nvar; k++)
{
Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),m_F[i][k],1,m_F[i][k],1);
}
}
else if(type == eIMEX)
{
if(m_numstages==1)
{
m_T= t_old[0]+timestep;
}
else
{
m_T= t_old[0];
for(int j=0; j<=i; ++j)
{
m_T += A(i,j)*timestep;
}
}
if(fabs(A(i,i)) > NekConstants::kNekZeroTol)
{
op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
for(k = 0; k < m_nvar; k++)
{
Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),
m_F[i][k],1,m_F[i][k],1);
}
}
op.DoOdeRhs(m_Y, m_F_IMEX[i], m_T);
}
else if( type == eExplicit)
{
// ensure solution is in correct space
op.DoProjection(m_Y,m_Y,m_T);
op.DoOdeRhs(m_Y, m_F[i], 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.
int i_start = 0;
if( (m_lastStageEqualsNewSolution) )
{
for(k = 0; k < m_nvar; k++)
{
Vmath::Vcopy(m_npoints,m_Y[k],1,y_new[0][k],1);
}
if (m_numstages==1 && type == eIMEX)
{
t_new[0] = t_old[0]+timestep;
}
else
{
t_new[0] = B(0,0)*timestep;
for(j = 1; j < m_numstages; j++)
{
t_new[0] += B(0,j)*timestep;
}
for(j = 0; j < m_numsteps; j++)
{
t_new[0] += V(0,j)*t_old[j];
}
i_start = 1;
}
}
for(i = i_start; i < m_numsteps; i++)
{
// The solution at the new time level is a linear
// combination of:
// 1: the stage derivatives
for(k = 0; k < m_nvar; k++)
{
Vmath::Smul(m_npoints,timestep*B(i,0),m_F[0][k],1,
y_new[i][k],1);
if(type == eIMEX)
{
Vmath::Svtvp(m_npoints,timestep*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)*timestep;
}
for(j = 1; j < m_numstages; j++)
{
for(k = 0; k < m_nvar; k++)
{
Vmath::Svtvp(m_npoints,timestep*B(i,j),m_F[j][k],1,
y_new[i][k],1,y_new[i][k],1);
if(type == eIMEX)
{
Vmath::Svtvp(m_npoints,timestep*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)*timestep;
}
}
// 2: the imported multi-step solution of the previous
// time level
for(j = 0; j < m_numsteps; j++)
{
for(k = 0; k < m_nvar; k++)
{
Vmath::Svtvp(m_npoints,V(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)
{
op.DoProjection(y_new[0],y_new[0],t_new[0]);
}
}
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::U ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 399 of file TimeIntegrationScheme.h.

References m_U.

Referenced by Nektar::LibUtilities::operator<<(), and TimeIntegrate().

{
return m_U[i][j];
}
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::V ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 404 of file TimeIntegrationScheme.h.

References m_V.

Referenced by Nektar::LibUtilities::operator<<(), and TimeIntegrate().

{
return m_V[i][j];
}
bool Nektar::LibUtilities::TimeIntegrationScheme::VerifyIntegrationSchemeType ( TimeIntegrationSchemeType  type,
const Array< OneD, const Array< TwoD, NekDouble > > &  A,
const Array< OneD, const Array< TwoD, NekDouble > > &  B,
const Array< TwoD, const NekDouble > &  U,
const Array< TwoD, const NekDouble > &  V 
) const
private

Definition at line 1040 of file TimeIntegrationScheme.cpp.

References A(), ASSERTL1, B(), Nektar::LibUtilities::eDiagonallyImplicit, Nektar::LibUtilities::eExplicit, Nektar::LibUtilities::eIMEX, Nektar::LibUtilities::eImplicit, and Nektar::NekConstants::kNekZeroTol.

Referenced by TimeIntegrationScheme().

{
int i;
int j;
int m;
int IMEXdim = A.num_elements();
int dim = A[0].GetRows();
Array<OneD, TimeIntegrationSchemeType> vertype(IMEXdim,eExplicit);
for(m = 0; m < IMEXdim; m++)
{
for(i = 0; i < dim; i++)
{
if( fabs(A[m][i][i]) > NekConstants::kNekZeroTol )
{
vertype[m] = eDiagonallyImplicit;
}
}
for(i = 0; i < dim; i++)
{
for(j = i+1; j < dim; j++)
{
if( fabs(A[m][i][j]) > NekConstants::kNekZeroTol )
{
vertype[m] = eImplicit;
ASSERTL1(false,"Fully Impplicit schemes cannnot be handled by the TimeIntegrationScheme class");
}
}
}
}
if(IMEXdim == 2)
{
ASSERTL1(B.num_elements()==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");
}
}
return (vertype[0] == type);
}

Friends And Related Function Documentation

friend class Nektar::MemoryManager
friend

Time at the different stages.

Definition at line 546 of file TimeIntegrationScheme.h.

TimeIntegrationSchemeManagerT& TimeIntegrationSchemeManager ( void  )
friend

Definition at line 45 of file TimeIntegrationScheme.cpp.

{
TimeIntegrationSchemeManagerT& m = Loki::SingletonHolder<TimeIntegrationSchemeManagerT>::Instance();
m.RegisterGlobalCreator(TimeIntegrationScheme::Create);
return m;
}

Member Data Documentation

Array<OneD, Array<TwoD,NekDouble> > Nektar::LibUtilities::TimeIntegrationScheme::m_A
protected

Definition at line 529 of file TimeIntegrationScheme.h.

Referenced by A(), A_IMEX(), and TimeIntegrationScheme().

Array<OneD, Array<TwoD,NekDouble> > Nektar::LibUtilities::TimeIntegrationScheme::m_B
protected

Definition at line 530 of file TimeIntegrationScheme.h.

Referenced by B(), B_IMEX(), and TimeIntegrationScheme().

TripleArray Nektar::LibUtilities::TimeIntegrationScheme::m_F
private

explicit right hand side of each stage equation

Definition at line 541 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

TripleArray Nektar::LibUtilities::TimeIntegrationScheme::m_F_IMEX
private

Array corresponding to the stage Derivatives.

Definition at line 542 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

bool Nektar::LibUtilities::TimeIntegrationScheme::m_firstStageEqualsOldSolution
protected

Definition at line 511 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate(), and TimeIntegrationScheme().

bool Nektar::LibUtilities::TimeIntegrationScheme::m_initialised
private

Definition at line 535 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

bool Nektar::LibUtilities::TimeIntegrationScheme::m_lastStageEqualsNewSolution
protected

Definition at line 512 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate(), and TimeIntegrationScheme().

int Nektar::LibUtilities::TimeIntegrationScheme::m_npoints
private

The number of variables in integration scheme.

Definition at line 537 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

unsigned int Nektar::LibUtilities::TimeIntegrationScheme::m_numMultiStepDerivs
protected

Definition at line 516 of file TimeIntegrationScheme.h.

Referenced by GetNmultiStepDerivs(), InitializeScheme(), TimeIntegrate(), and TimeIntegrationScheme().

unsigned int Nektar::LibUtilities::TimeIntegrationScheme::m_numMultiStepValues
protected

Definition at line 514 of file TimeIntegrationScheme.h.

Referenced by GetNmultiStepValues(), InitializeScheme(), TimeIntegrate(), and TimeIntegrationScheme().

unsigned int Nektar::LibUtilities::TimeIntegrationScheme::m_numstages
protected

Definition at line 509 of file TimeIntegrationScheme.h.

Referenced by CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), GetNstages(), TimeIntegrate(), and TimeIntegrationScheme().

unsigned int Nektar::LibUtilities::TimeIntegrationScheme::m_numsteps
protected

Definition at line 508 of file TimeIntegrationScheme.h.

Referenced by CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), CheckTimeIntegrateArguments(), GetNsteps(), TimeIntegrate(), and TimeIntegrationScheme().

int Nektar::LibUtilities::TimeIntegrationScheme::m_nvar
private

bool to identify if array has been initialised

Definition at line 536 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

TimeIntegrationSchemeKey Nektar::LibUtilities::TimeIntegrationScheme::m_schemeKey
protected

Definition at line 506 of file TimeIntegrationScheme.h.

Referenced by GetIntegrationMethod(), GetIntegrationSchemeKey(), and InitializeScheme().

TimeIntegrationSchemeType Nektar::LibUtilities::TimeIntegrationScheme::m_schemeType
protected

Definition at line 507 of file TimeIntegrationScheme.h.

Referenced by GetIntegrationSchemeType(), and TimeIntegrationScheme().

NekDouble Nektar::LibUtilities::TimeIntegrationScheme::m_T
private

Used to store the Explicit stage derivative of IMEX schemes.

Definition at line 544 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

Array<OneD,unsigned int> Nektar::LibUtilities::TimeIntegrationScheme::m_timeLevelOffset
protected

Definition at line 518 of file TimeIntegrationScheme.h.

Referenced by GetTimeLevelOffset(), InitializeScheme(), TimeIntegrate(), and TimeIntegrationScheme().

DoubleArray Nektar::LibUtilities::TimeIntegrationScheme::m_tmp
private

Array containing the stage values.

Definition at line 539 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

Array<TwoD,NekDouble> Nektar::LibUtilities::TimeIntegrationScheme::m_U
protected

Definition at line 531 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme(), and U().

Array<TwoD,NekDouble> Nektar::LibUtilities::TimeIntegrationScheme::m_V
protected

Definition at line 532 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme(), and V().

DoubleArray Nektar::LibUtilities::TimeIntegrationScheme::m_Y
private

The size of inner data which is stored for reuse.

Definition at line 538 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

Generated on Sun Mar 15 2015 00:11:55 for Nektar++ by   doxygen 1.8.1.2