46 #include <boost/core/ignore_unused.hpp>
52 namespace LibUtilities
66 this, y_0, time, deltaT);
83 int nvar = y_0.size();
84 int npoints = y_0[0].size();
86 for (i = 0; i < nvar; i++)
94 for (i = 0; i < nvar; i++)
98 y_out->SetDerivative(0, f_y_0, deltaT);
112 int nvar = solvector->GetFirstDim();
113 int npoints = solvector->GetSecondDim();
115 if (solvector->GetIntegrationSchemeData() !=
150 unsigned int nCurSchemeVals =
153 unsigned int nCurSchemeDers =
156 unsigned int nCurSchemeSteps =
161 unsigned int nMasterSchemeVals =
162 solvector->GetNvalues();
163 unsigned int nMasterSchemeDers =
164 solvector->GetNderivs();
171 solvector->GetTimeLevelOffset();
180 for (
int n = 0; n < nCurSchemeVals; n++)
186 y_n = solvector->GetValue(curTimeLevels[n]);
187 t_n = solvector->GetValueTime(curTimeLevels[n]);
191 solvector_in->SetValue(curTimeLevels[n], y_n, t_n);
194 for (
int n = nCurSchemeVals; n < nCurSchemeSteps; n++)
200 dtFy_n = solvector->GetDerivative(curTimeLevels[n]);
204 solvector_in->SetDerivative(curTimeLevels[n], dtFy_n, deltaT);
216 solvector_in->GetTimeVector(),
217 solvector_out->UpdateSolutionVector(),
218 solvector_out->UpdateTimeVector(), op);
230 bool CalcNewDeriv =
false;
236 if (nMasterSchemeDers > 0)
238 if (nCurSchemeDers == 0)
244 if (masterTimeLevels[nMasterSchemeVals] <
245 curTimeLevels[nCurSchemeVals])
254 int newDerivTimeLevel =
255 masterTimeLevels[nMasterSchemeVals];
267 if (newDerivTimeLevel == 0)
269 y_n = solvector_out->GetValue(0);
270 t_n = solvector_out->GetValueTime(0);
275 else if (newDerivTimeLevel == 1)
277 y_n = solvector->GetValue(0);
278 t_n = solvector->GetValueTime(0);
282 ASSERTL1(
false,
"Problems with initialising scheme");
286 for (
int j = 0; j < nvar; j++)
296 for (
int j = 0; j < nvar; j++)
298 Vmath::Smul(npoints, deltaT, f_n[j], 1, f_n[j], 1);
303 solvector->RotateSolutionVector();
307 solvector->SetDerivative(newDerivTimeLevel, f_n, deltaT);
313 solvector->RotateSolutionVector();
319 for (
int n = 0; n < nCurSchemeVals; n++)
327 y_n = solvector_out->GetValue(curTimeLevels[n]);
328 t_n = solvector_out->GetValueTime(curTimeLevels[n]);
331 solvector->SetValue(curTimeLevels[n], y_n, t_n);
334 for (
int n = nCurSchemeVals; n < nCurSchemeSteps; n++)
340 dtFy_n = solvector_out->GetDerivative(curTimeLevels[n]);
344 solvector->SetDerivative(curTimeLevels[n], dtFy_n, deltaT);
351 this, nvar, npoints);
354 solvector->GetTimeVector(),
355 solvector_new->UpdateSolutionVector(),
356 solvector_new->UpdateTimeVector(), op);
358 solvector = solvector_new;
361 return solvector->GetSolution();
373 "Arguments not well defined");
382 for (
int j = 0; j <
m_nvars; j++)
405 for (
int j = 0; j <
m_nvars; j++)
419 for (
int j = 0; j <
m_nvars; j++)
431 for (
int j = 0; j <
m_nvars; j++)
443 for (
int j = 0; j <
m_nvars; j++)
462 ->InitializeSecondaryData(
this, deltaT);
474 for (
int k = 0; k <
m_nvars; k++)
487 for (
int k = 0; k <
m_nvars; k++)
509 m_T =
A(stage, 0) * deltaT;
511 for (
int j = 1; j < stage; j++)
513 for (
int k = 0; k <
m_nvars; k++)
534 m_T +=
A(stage, j) * deltaT;
540 for (
int k = 0; k <
m_nvars; k++)
554 m_T +=
U(stage, j) * t_old[j];
569 m_T = t_old[0] + deltaT;
574 for (
int j = 0; j <= stage; ++j)
576 m_T +=
A(stage, j) * deltaT;
582 for (
int k = 0; k <
m_nvars; ++k)
587 m_F[stage][k], 1,
m_F[stage][k], 1);
591 else if (type ==
eIMEX)
595 m_T = t_old[0] + deltaT;
600 for (
int j = 0; j <= stage; ++j)
602 m_T +=
A(stage, j) * deltaT;
610 for (
int k = 0; k <
m_nvars; k++)
615 m_F[stage][k], 1,
m_F[stage][k], 1);
648 for (
int k = 0; k <
m_nvars; k++)
655 t_new[0] = t_old[0] + deltaT;
659 t_new[0] =
B(0, 0) * deltaT;
663 t_new[0] +=
B(0, j) * deltaT;
668 t_new[0] +=
V(0, j) * t_old[j];
680 for (
int k = 0; k <
m_nvars; k++)
696 1, y_new[i][k], 1, y_new[i][k], 1);
702 t_new[i] =
B(i, 0) * deltaT;
707 for (
int k = 0; k <
m_nvars; k++)
712 m_F[j][k], 1, y_new[i][k], 1, y_new[i][k], 1);
717 y_new[i][k], 1, y_new[i][k], 1);
723 m_F_IMEX[j][k], 1, y_new[i][k], 1, y_new[i][k],
730 t_new[i] +=
B(i, j) * deltaT;
738 for (
int k = 0; k <
m_nvars; k++)
743 y_new[i][k], 1, y_new[i][k], 1);
748 y_new[i][k], 1, y_new[i][k], 1);
754 t_new[i] +=
V(i, j) * t_old[j];
774 for (
int m = 0; m <
m_A.size(); m++)
814 for (
int m = 0; m <
m_A.size(); m++)
844 int IMEXdim =
m_A.size();
845 int dim =
m_A[0].GetRows();
852 for (
int m = 0; m < IMEXdim; m++)
854 for (
int i = 0; i < dim; i++)
862 for (
int i = 0; i < dim; i++)
864 for (
int j = i + 1; j < dim; j++)
869 ASSERTL1(
false,
"Fully Implicit schemes cannnot be handled "
870 "by the TimeIntegrationScheme class");
879 "Coefficient Matrix B should have an "
880 "implicit and explicit part for IMEX schemes");
888 ASSERTL1(
false,
"This is not a proper IMEX scheme");
893 "Time integration scheme coefficients do not match its type");
902 boost::ignore_unused(op);
904 boost::ignore_unused(y_old, t_old, y_new, t_new, op);
912 ASSERTL1(y_old[0].size() == y_new[0].size(),
913 "Non-matching number of variables.");
914 ASSERTL1(y_old[0][0].size() == y_new[0][0].size(),
915 "Non-matching number of coefficients.");
942 <<
"Time Integration Phase : " << rhs.
m_name <<
"\n"
943 <<
"- number of steps: " << r <<
"\n"
944 <<
"- number of stages: " << s <<
"\n"
945 <<
"General linear method tableau:\n";
947 for (
int i = 0; i < s; i++)
949 for (
int j = 0; j < s; j++)
952 os.precision(osprecision);
953 os << std::right << rhs.
A(i, j) <<
" ";
958 for (
int j = 0; j < s; j++)
961 os.precision(osprecision);
962 os << std::right << rhs.
A_IMEX(i, j) <<
" ";
968 for (
int j = 0; j < r; j++)
971 os.precision(osprecision);
972 os << std::right << rhs.
U(i, j);
977 int imexflag = (type ==
eIMEX) ? 2 : 1;
978 for (
int i = 0; i < (r + imexflag * s) * (oswidth + 1) + imexflag * 2 - 1;
985 for (
int i = 0; i < r; i++)
987 for (
int j = 0; j < s; j++)
990 os.precision(osprecision);
991 os << std::right << rhs.
B(i, j) <<
" ";
996 for (
int j = 0; j < s; j++)
999 os.precision(osprecision);
1000 os << std::right << rhs.
B_IMEX(i, j) <<
" ";
1006 for (
int j = 0; j < r; j++)
1009 os.precision(osprecision);
1010 os << std::right << rhs.
V(i, j);
1016 os.precision(osprecision);
1021 os << std::right <<
" value";
1025 os << std::right <<
" derivative";
1033 for (
int k = 0; k < rhs.
m_nvars; k++)
1036 <<
"General linear method exponential tableau for variable " << k
1039 for (
int i = 0; i < s; i++)
1041 for (
int j = 0; j < s; j++)
1044 os.precision(osprecision);
1045 os << std::right << rhs.
m_A_phi[k][i][j] <<
" ";
1050 for (
int j = 0; j < r; j++)
1053 os.precision(osprecision);
1054 os << std::right << rhs.
m_U_phi[k][i][j];
1059 int imexflag = (type ==
eIMEX) ? 2 : 1;
1061 i < (r + imexflag * s) * (oswidth + 1) + imexflag * 2 - 1; i++)
1067 for (
int i = 0; i < r; i++)
1069 for (
int j = 0; j < s; j++)
1072 os.precision(osprecision);
1073 os << std::right << rhs.
m_B_phi[k][i][j] <<
" ";
1078 for (
int j = 0; j < r; j++)
1081 os.precision(osprecision);
1082 os << std::right << rhs.
m_V_phi[k][i][j];
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
void CheckIfLastStageEqualsNewSolution()
TimeIntegrationSchemeType m_schemeType
Array< TwoD, NekDouble > m_U
int GetFirstDim(ConstTripleArray &y) const
void VerifyIntegrationSchemeType()
int GetSecondDim(ConstTripleArray &y) const
NekDouble A(const unsigned int i, const unsigned int j) const
Array< OneD, Array< TwoD, NekDouble > > m_A_phi
unsigned int m_numMultiStepValues
Array< TwoD, NekDouble > m_V
NekDouble m_lastNVars
Last delta T.
TripleArray m_F
Explicit RHS of each stage equation.
bool CheckTimeIntegrateArguments(ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new, SingleArray &t_new, const TimeIntegrationSchemeOperators &op) const
void CheckIfFirstStageEqualsOldSolution()
unsigned int GetNsteps(void) const
Array< OneD, Array< TwoD, NekDouble > > m_B_phi
int m_npoints
The size of inner data which is stored for reuse.
DoubleArray m_tmp
Array containing the stage values.
NekDouble B_IMEX(const unsigned int i, const unsigned int j) const
NekDouble V(const unsigned int i, const unsigned int j) const
Array< OneD, Array< TwoD, NekDouble > > m_B
bool m_lastStageEqualsNewSolution
const TimeIntegrationScheme * m_parent
Parent scheme object.
Array< OneD, Array< TwoD, NekDouble > > m_A
unsigned int GetNmultiStepValues() const
Array< OneD, Array< TwoD, NekDouble > > m_V_phi
NekDouble m_T
Used to store the Explicit stage derivative of IMEX schemes.
const Array< OneD, const unsigned int > & GetTimeLevelOffset() const
bool m_firstStageEqualsOldSolution
Time at the different stages.
unsigned int GetNmultiStepDerivs() const
NekDouble A_IMEX(const unsigned int i, const unsigned int j) const
NekDouble B(const unsigned int i, const unsigned int j) const
LUE ConstDoubleArray & TimeIntegrate(const NekDouble deltaT, TimeIntegrationSolutionGLMSharedPtr &y, const TimeIntegrationSchemeOperators &op)
Explicit integration of an ODE.
LUE TimeIntegrationSolutionGLMSharedPtr InitializeData(const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)
This function initialises the time integration scheme.
int m_nvars
Last number of vars.
NekDouble m_lastDeltaT
bool to identify array initialization
NekDouble U(const unsigned int i, const unsigned int j) const
TripleArray m_F_IMEX
Array corresponding to the stage Derivatives.
Array< OneD, unsigned int > m_timeLevelOffset
Array< OneD, Array< TwoD, NekDouble > > m_U_phi
Base class for GLM time integration schemes.
virtual LUE std::string GetFullName() const
Binds a set of functions for use by time integration schemes.
void DoImplicitSolve(InArrayType &inarray, OutArrayType &outarray, const NekDouble time, const NekDouble lambda) const
void DoProjection(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
void DoOdeRhs(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
AT< AT< NekDouble > > DoubleArray
std::shared_ptr< TimeIntegrationAlgorithmGLM > TimeIntegrationAlgorithmGLMSharedPtr
TimeIntegrationSchemeType
@ eImplicit
Fully implicit scheme.
@ eExplicit
Formally explicit scheme.
@ eExponential
Exponential scheme.
@ eDiagonallyImplicit
Diagonally implicit scheme (e.g. the DIRK schemes)
@ eIMEX
Implicit Explicit General Linear Method.
std::shared_ptr< TimeIntegrationSolutionGLM > TimeIntegrationSolutionGLMSharedPtr
AT< AT< AT< NekDouble > > > TripleArray
std::ostream & operator<<(std::ostream &os, const BasisKey &rhs)
static const NekDouble kNekZeroTol
The above copyright notice and this permission notice shall be included.
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.