Nektar++
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>

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 NekDoubleConstSingleArray
 
typedef Array< OneD, NekDoubleSingleArray
 
typedef std::function< void(ConstDoubleArray &, DoubleArray &, const NekDouble) > FunctorType1
 
typedef std::function< void(ConstDoubleArray &, DoubleArray &, const NekDouble, const NekDouble) > FunctorType2
 

Public Member Functions

virtual ~TimeIntegrationScheme ()
 
const TimeIntegrationSchemeKeyGetIntegrationSchemeKey () 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. More...
 
ConstDoubleArrayTimeIntegrate (const NekDouble timestep, TimeIntegrationSolutionSharedPtr &y, const TimeIntegrationSchemeOperators &op)
 Explicit integration of an ODE. More...
 
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, NekDoublem_U
 
Array< TwoD, NekDoublem_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 std::shared_ptr< TimeIntegrationSchemeCreate (const TimeIntegrationSchemeKey &key)
 

Private Attributes

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

Friends

template<typename >
class Nektar::MemoryManager
 Time at the different stages. More...
 
TimeIntegrationSchemeManagerTTimeIntegrationSchemeManager (void)
 

Detailed Description

Definition at line 378 of file TimeIntegrationScheme.h.

Member Typedef Documentation

◆ ConstDoubleArray

Definition at line 384 of file TimeIntegrationScheme.h.

◆ ConstSingleArray

Definition at line 386 of file TimeIntegrationScheme.h.

◆ ConstTripleArray

Definition at line 382 of file TimeIntegrationScheme.h.

◆ DoubleArray

Definition at line 385 of file TimeIntegrationScheme.h.

◆ FunctorType1

Definition at line 388 of file TimeIntegrationScheme.h.

◆ FunctorType2

Definition at line 389 of file TimeIntegrationScheme.h.

◆ SingleArray

Definition at line 387 of file TimeIntegrationScheme.h.

◆ TripleArray

Definition at line 383 of file TimeIntegrationScheme.h.

Constructor & Destructor Documentation

◆ ~TimeIntegrationScheme()

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

Definition at line 393 of file TimeIntegrationScheme.h.

394  {
395  }

◆ TimeIntegrationScheme() [1/3]

Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrationScheme ( const TimeIntegrationSchemeKey key)
private

Definition at line 149 of file TimeIntegrationScheme.cpp.

References ASSERTL1, CheckIfFirstStageEqualsOldSolution(), CheckIfLastStageEqualsNewSolution(), Nektar::LibUtilities::eAdamsBashforthOrder1, Nektar::LibUtilities::eAdamsBashforthOrder2, Nektar::LibUtilities::eAdamsBashforthOrder3, Nektar::LibUtilities::eAdamsBashforthOrder4, 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, Nektar::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::eIMEXOrder4, Nektar::LibUtilities::eMCNAB, Nektar::LibUtilities::eMidpoint, Nektar::LibUtilities::eRungeKutta2, Nektar::LibUtilities::eRungeKutta2_ImprovedEuler, Nektar::LibUtilities::eRungeKutta2_SSP, Nektar::LibUtilities::eRungeKutta3_SSP, Nektar::LibUtilities::eRungeKutta4, Nektar::LibUtilities::eRungeKutta5, 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().

149  :
150  m_schemeKey(key),
151  m_initialised(false)
152  {
153  switch(key.GetIntegrationMethod())
154  {
155  case eForwardEuler:
157  {
158  m_numsteps = 1;
159  m_numstages = 1;
160 
161  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
162  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
163 
164  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
165  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
166  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
167  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 1.0);
168 
172  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
173  m_timeLevelOffset[0] = 0;
174  }
175  break;
177  {
178  m_numsteps = 2;
179  m_numstages = 1;
180 
181  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
182  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
183 
184  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
185  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
186  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
187  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
188 
189  m_B[0][0][0] = 3.0/2.0;
190  m_B[0][1][0] = 1.0;
191 
192  m_U[0][0] = 1.0;
193 
194  m_V[0][0] = 1.0;
195  m_V[0][1] = -0.5;
196 
200  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
201  m_timeLevelOffset[0] = 0;
202  m_timeLevelOffset[1] = 1;
203  }
204  break;
206  {
207  m_numsteps = 4;
208  m_numstages = 1;
209 
210  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
211  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
212 
213  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
214  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
215  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
216  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
217 
218  m_B[0][1][0] = 1.0;
219 
220  m_U[0][0] = 1.0;
221  m_U[0][1] = 23.0/12.0;
222  m_U[0][2] = -4.0/3.0;
223  m_U[0][3] = 5.0/12.0;
224 
225  m_V[0][0] = 1.0;
226  m_V[0][1] = 23.0/12.0;
227  m_V[0][2] = -4.0/3.0;
228  m_V[0][3] = 5.0/12.0;
229  m_V[2][1] = 1.0;
230  m_V[3][2] = 1.0;
231 
235  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
236  m_timeLevelOffset[0] = 0;
237  m_timeLevelOffset[1] = 1;
238  m_timeLevelOffset[2] = 2;
239  m_timeLevelOffset[3] = 3;
240  }
241  break;
243  {
244  m_numsteps = 5;
245  m_numstages = 1;
246 
247  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
248  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
249 
250  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
251  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
252  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
253  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
254 
255  m_B[0][1][0] = 1.0;
256 
257  m_U[0][0] = 1.0;
258  m_U[0][1] = 55.0/24.0;
259  m_U[0][2] = -59.0/24.0;
260  m_U[0][3] = 37.0/24.0;
261  m_U[0][4] = -9.0/24.0;
262 
263  m_V[0][0] = 1.0;
264  m_V[0][1] = 55.0/24.0;
265  m_V[0][2] = -59.0/24.0;
266  m_V[0][3] = 37.0/24.0;
267  m_V[0][4] = -9.0/24.0;
268  m_V[2][1] = 1.0;
269  m_V[3][2] = 1.0;
270  m_V[4][3] = 1.0;
271 
275  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
276  m_timeLevelOffset[0] = 0;
277  m_timeLevelOffset[1] = 1;
278  m_timeLevelOffset[2] = 2;
279  m_timeLevelOffset[3] = 3;
280  m_timeLevelOffset[4] = 4;
281  }
282  break;
283  case eBackwardEuler:
284  case eBDFImplicitOrder1:
285  case eAdamsMoultonOrder1:
286  {
287  m_numsteps = 1;
288  m_numstages = 1;
289 
290  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
291  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
292 
293  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,1.0);
294  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
295  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
296  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 1.0);
297 
301  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
302  m_timeLevelOffset[0] = 0; // SJS: Not sure whether this is correct
303  }
304  break;
305  case eIMEXOrder1:
306  {
307  m_numsteps = 2;
308  m_numstages = 1;
309 
310  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
311  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
312 
313  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,1.0);
314  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
315  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
316  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
317  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
318  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
319 
320  m_B[0][1][0] = 0.0;
321  m_B[1][0][0] = 0.0;
322  m_V[0][0] = 1.0;
323  m_V[0][1] = 1.0;
324 
328  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
329  m_timeLevelOffset[0] = 0;
330  m_timeLevelOffset[1] = 0;
331  }
332  break;
333  case eIMEXOrder2:
334  {
335  NekDouble third = 1.0/3.0;
336  m_numsteps = 4;
337  m_numstages = 1;
338 
339  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
340  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
341 
342  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,2*third);
343  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
344  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
345  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
346  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 4*third);
347  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
348 
349  m_B[0][0][0] = 2*third;
350  m_B[1][2][0] = 1.0;
351  m_U[0][1] = -third;
352  m_U[0][3] = -2*third;
353 
354  m_V[0][0] = 4*third;
355  m_V[0][1] = -third;
356  m_V[0][2] = 4*third;
357  m_V[0][3] = -2*third;
358  m_V[1][0] = 1.0;
359  m_V[3][2] = 1.0;
360 
364  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
365  m_timeLevelOffset[0] = 0;
366  m_timeLevelOffset[1] = 1;
367  m_timeLevelOffset[2] = 0;
368  m_timeLevelOffset[3] = 1;
369  }
370  break;
371  case eIMEXOrder3:
372  {
373  NekDouble eleventh = 1.0/11.0;
374  m_numsteps = 6;
375  m_numstages = 1;
376 
377  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
378  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
379 
380  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,6*eleventh);
381  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
382  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
383  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
384  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 18*eleventh);
385  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
386 
387  m_B[0][0][0] = 6*eleventh;
388  m_B[1][3][0] = 1.0;
389  m_U[0][1] = -9*eleventh;
390  m_U[0][2] = 2*eleventh;
391  m_U[0][4] = -18*eleventh;
392  m_U[0][5] = 6*eleventh;
393 
394  m_V[0][0] = 18*eleventh;
395  m_V[0][1] = -9*eleventh;
396  m_V[0][2] = 2*eleventh;
397  m_V[0][3] = 18*eleventh;
398  m_V[0][4] = -18*eleventh;
399  m_V[0][5] = 6*eleventh;
400  m_V[1][0] = 1.0;
401  m_V[2][1] = 1.0;
402  m_V[4][3] = 1.0;
403  m_V[5][4] = 1.0;
404 
408  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
409  m_timeLevelOffset[0] = 0;
410  m_timeLevelOffset[1] = 1;
411  m_timeLevelOffset[2] = 2;
412  m_timeLevelOffset[3] = 0;
413  m_timeLevelOffset[4] = 1;
414  m_timeLevelOffset[5] = 2;
415  }
416  break;
417  case eIMEXOrder4:
418  {
419  NekDouble twentyfifth = 1.0/25.0;
420  m_numsteps = 8;
421  m_numstages = 1;
422 
423  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
424  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
425 
426  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,12*twentyfifth);
427  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
428  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
429  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
430  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 48*twentyfifth);
431  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
432 
433  m_B[0][0][0] = 12*twentyfifth;
434  m_B[1][4][0] = 1.0;
435  m_U[0][1] = -36*twentyfifth;
436  m_U[0][2] = 16*twentyfifth;
437  m_U[0][3] = -3*twentyfifth;
438  m_U[0][5] = -72*twentyfifth;
439  m_U[0][7] = -12*twentyfifth;
440 
441  m_V[0][0] = 48*twentyfifth;
442  m_V[0][1] = -36*twentyfifth;
443  m_V[0][2] = 16*twentyfifth;
444  m_V[0][3] = -3*twentyfifth;
445  m_V[0][4] = 48*twentyfifth;
446  m_V[0][5] = -72*twentyfifth;
447  m_V[0][6] = 48*twentyfifth;
448  m_V[0][7] = -12*twentyfifth;
449  m_V[1][0] = 1.0;
450  m_V[2][1] = 1.0;
451  m_V[3][2] = 1.0;
452  m_V[5][4] = 1.0;
453  m_V[6][5] = 1.0;
454  m_V[7][6] = 1.0;
455 
459  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
460  m_timeLevelOffset[0] = 0;
461  m_timeLevelOffset[1] = 1;
462  m_timeLevelOffset[2] = 2;
463  m_timeLevelOffset[3] = 3;
464  m_timeLevelOffset[4] = 0;
465  m_timeLevelOffset[5] = 1;
466  m_timeLevelOffset[6] = 2;
467  m_timeLevelOffset[7] = 3;
468  }
469  break;
470  case eAdamsMoultonOrder2:
471  {
472  m_numsteps = 2;
473  m_numstages = 1;
474 
475  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
476  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
477 
478  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.5);
479  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
480  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
481  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.0);
482 
483  m_B[0][0][0] = 0.5;
484  m_B[0][1][0] = 1.0;
485 
486  m_U[0][0] = 1.0;
487  m_U[0][1] = 0.5;
488 
489  m_V[0][0] = 1.0;
490  m_V[0][1] = 0.5;
491 
495  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
496  m_timeLevelOffset[0] = 0;
497  m_timeLevelOffset[1] = 0;
498  }
499  break;
500  case eBDFImplicitOrder2:
501  {
502  NekDouble third = 1.0/3.0;
503  m_numsteps = 2;
504  m_numstages = 1;
505 
506  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
507  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
508 
509  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,2*third);
510  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
511  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
512  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.0);
513 
514  m_B[0][0][0] = 2*third;
515  m_B[0][1][0] = 0.0;
516 
517  m_U[0][0] = 4*third;
518  m_U[0][1] = -third;
519 
520  m_V[0][0] = 4*third;
521  m_V[0][1] = -third;
522  m_V[1][0] = 1.0;
523 
527  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
528  m_timeLevelOffset[0] = 0;
529  m_timeLevelOffset[1] = 1;
530  }
531  break;
532  case eMidpoint:
533  case eRungeKutta2:
534  {
535  m_numsteps = 1;
536  m_numstages = 2;
537 
538  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
539  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
540 
541  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
542  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
543  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
544  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
545 
546  m_A[0][1][0] = 0.5;
547  m_B[0][0][1] = 1.0;
548 
552  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
553  m_timeLevelOffset[0] = 0;
554  }
555  break;
557  case eRungeKutta2_SSP:
558  {
559  m_numsteps = 1;
560  m_numstages = 2;
561 
562  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
563  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
564 
565  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
566  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
567  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
568  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
569 
570  m_A[0][1][0] = 1.0;
571 
572  m_B[0][0][0] = 0.5;
573  m_B[0][0][1] = 0.5;
574 
578  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
579  m_timeLevelOffset[0] = 0;
580  }
581  break;
582  case eRungeKutta3_SSP:
583  {
584  m_numsteps = 1;
585  m_numstages = 3;
586 
587  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
588  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
589 
590  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
591  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
592  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
593  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
594 
595  m_A[0][1][0] = 1.0;
596  m_A[0][2][0] = 0.25;
597  m_A[0][2][1] = 0.25;
598 
599  m_B[0][0][0] = 1.0/6.0;
600  m_B[0][0][1] = 1.0/6.0;
601  m_B[0][0][2] = 2.0/3.0;
602 
606  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
607  m_timeLevelOffset[0] = 0;
608  }
609  break;
611  case eRungeKutta4:
612  {
613  m_numsteps = 1;
614  m_numstages = 4;
615 
616  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
617  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
618 
619  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
620  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
621  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
622  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
623 
624  m_A[0][1][0] = 0.5;
625  m_A[0][2][1] = 0.5;
626  m_A[0][3][2] = 1.0;
627 
628  m_B[0][0][0] = 1.0/6.0;
629  m_B[0][0][1] = 1.0/3.0;
630  m_B[0][0][2] = 1.0/3.0;
631  m_B[0][0][3] = 1.0/6.0;
632 
636  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
637  m_timeLevelOffset[0] = 0;
638  }
639  break;
640  case eRungeKutta5:
641  {
642  m_numsteps = 1;
643  m_numstages = 6;
644 
645  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
646  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
647 
648  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
649  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
650  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
651  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
652 
653  m_A[0][1][0] = 1.0/4.0;
654  m_A[0][2][0] = 1.0/8.0;
655  m_A[0][2][1] = 1.0/8.0;
656  m_A[0][3][2] = 1.0/2.0;
657  m_A[0][4][0] = 3.0/16.0;
658  m_A[0][4][1] = -3.0/8.0;
659  m_A[0][4][2] = 3.0/8.0;
660  m_A[0][4][3] = 9.0/16.0;
661  m_A[0][5][0] = -3.0/7.0;
662  m_A[0][5][1] = 8.0/7.0;
663  m_A[0][5][2] = 6.0/7.0;
664  m_A[0][5][3] = -12.0/7.0;
665  m_A[0][5][4] = 8.0/7.0;
666 
667  m_B[0][0][0] = 7.0/90.0;
668  m_B[0][0][2] = 32.0/90.0;
669  m_B[0][0][3] = 12.0/90.0;
670  m_B[0][0][4] = 32.0/90.0;
671  m_B[0][0][5] = 7.0/90.0;
672 
676  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
677  m_timeLevelOffset[0] = 0;
678  }
679  break;
680  case eDIRKOrder2:
681  {
682  m_numsteps = 1;
683  m_numstages = 2;
684 
685  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
686  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
687 
688  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
689  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
690  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
691  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
692 
693  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
694 
695  m_A[0][0][0] = lambda;
696  m_A[0][1][0] = 1.0 - lambda;
697  m_A[0][1][1] = lambda;
698 
699  m_B[0][0][0] = 1.0 - lambda;
700  m_B[0][0][1] = lambda;
701 
705  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
706  m_timeLevelOffset[0] = 0;
707  }
708  break;
709  case eDIRKOrder3:
710  {
711  m_numsteps = 1;
712  m_numstages = 3;
713 
714  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
715  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
716 
717  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
718  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
719  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
720  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
721 
722  NekDouble lambda = 0.4358665215084589;
723 
724  m_A[0][0][0] = lambda;
725  m_A[0][1][0] = 0.5 * (1.0 - lambda);
726  m_A[0][2][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
727  m_A[0][1][1] = lambda;
728  m_A[0][2][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
729  m_A[0][2][2] = lambda;
730 
731  m_B[0][0][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
732  m_B[0][0][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
733  m_B[0][0][2] = lambda;
734 
738  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
739  m_timeLevelOffset[0] = 0;
740  }
741  break;
742  case eIMEXdirk_2_3_2:
743  {
744  m_numsteps = 1;
745  m_numstages = 3;
746 
747  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
748  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
749 
750  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
751  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
752  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
753  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
754  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
755  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
756 
757  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
758  NekDouble delta = -2.0*sqrt(2.0)/3.0;
759 
760  m_A[0][1][1] = lambda;
761  m_A[0][2][1] = 1.0 - lambda;
762  m_A[0][2][2] = lambda;
763 
764  m_B[0][0][1] = 1.0 - lambda;
765  m_B[0][0][2] = lambda;
766 
767  m_A[1][1][0] = lambda;
768  m_A[1][2][0] = delta;
769  m_A[1][2][1] = 1.0 - delta;
770 
771  m_B[1][0][1] = 1.0 - lambda;
772  m_B[1][0][2] = lambda;
773 
777  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
778  m_timeLevelOffset[0] = 0;
779  }
780  break;
781  case eIMEXdirk_3_4_3:
782  {
783  m_numsteps = 1;
784  m_numstages = 4;
785 
786  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
787  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
788 
789  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
790  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
791  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
792  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
793  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
794  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
795 
796  NekDouble lambda = 0.4358665215;
797 
798  m_A[0][1][1] = lambda;
799  m_A[0][2][1] = 0.5 * (1.0 - lambda);
800  m_A[0][3][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
801  m_A[0][2][2] = lambda;
802  m_A[0][3][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
803  m_A[0][3][3] = lambda;
804 
805  m_B[0][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
806  m_B[0][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
807  m_B[0][0][3] = lambda;
808 
809  m_A[1][1][0] = 0.4358665215;
810  m_A[1][2][0] = 0.3212788860;
811  m_A[1][2][1] = 0.3966543747;
812  m_A[1][3][0] =-0.105858296;
813  m_A[1][3][1] = 0.5529291479;
814  m_A[1][3][2] = 0.5529291479;
815 
816  m_B[1][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
817  m_B[1][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
818  m_B[1][0][3] = lambda;
819 
823  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
824  m_timeLevelOffset[0] = 0;
825  }
826  break;
827  case eCNAB:
828  {
829  NekDouble secondth = 1.0/2.0;
830  m_numsteps = 4;
831  m_numstages = 1;
832 
833  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
834  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
835 
836  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,secondth);
837  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
838  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
839  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
840  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,secondth);
841  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
842 
843  m_B[0][0][0] = secondth;
844  m_B[0][1][0] = 1.0;
845  m_B[1][2][0] = 1.0;
846  m_U[0][0] = 2*secondth;
847  m_U[0][2] = 3*secondth;
848  m_U[0][3] = -1*secondth;
849 
850  m_V[0][0] = 2*secondth;
851  m_V[0][1] = secondth;
852  m_V[0][2] = 3*secondth;
853  m_V[0][3] = -1*secondth;
854  m_V[3][2] = 1.0;
855 
859  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
860  m_timeLevelOffset[0] = 0;
861  m_timeLevelOffset[1] = 0;
862  m_timeLevelOffset[2] = 0;
863  m_timeLevelOffset[3] = 1;
864  }
865  break;
866  case eIMEXGear:
867  {
868  NekDouble twothirdth = 2.0/3.0;
869  m_numsteps = 3;
870  m_numstages = 1;
871 
872  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
873  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
874 
875  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
876  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
877  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
878  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
879  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
880  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
881 
882  m_B[0][0][0] = twothirdth;
883  m_B[1][2][0] = 1.0;
884  m_U[0][0] = 2*twothirdth;
885  m_U[0][1] = -0.5*twothirdth;
886  m_U[0][2] = twothirdth;
887 
888  m_V[0][0] = 2*twothirdth;
889  m_V[0][1] = -0.5*twothirdth;
890  m_V[0][2] = twothirdth;
891  m_V[1][0] = 1.0;
892 
896  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
897  m_timeLevelOffset[0] = 0;
898  m_timeLevelOffset[1] = 1;
899  m_timeLevelOffset[2] = 0;
900  }
901  break;
902  case eMCNAB:
903  {
904  NekDouble sixthx = 9.0/16.0;
905  m_numsteps = 5;
906  m_numstages = 1;
907 
908  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
909  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
910 
911  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,sixthx);
912  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
913  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
914  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
915  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
916  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
917 
918  m_B[0][0][0] = sixthx;
919  m_B[0][1][0] = 1.0;
920  m_B[1][3][0] = 1.0;
921  m_U[0][0] = 1.0;
922  m_U[0][1] = 6.0/16.0;
923  m_U[0][2] = 1.0/16.0;
924  m_U[0][3] = 1.5;
925  m_U[0][4] = -0.5;
926 
927  m_V[0][0] = 1.0;
928  m_V[0][1] = 6.0/16.0;
929  m_V[0][2] = 1.0/16.0;
930  m_V[0][3] = 1.5;
931  m_V[0][4] = -0.5;
932  m_V[2][1] = 1.0;
933  m_V[4][3] = 1.0;
934 
938  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
939  m_timeLevelOffset[0] = 0;
940  m_timeLevelOffset[1] = 0;
941  m_timeLevelOffset[2] = 1;
942  m_timeLevelOffset[3] = 0;
943  m_timeLevelOffset[4] = 1;
944  }
945  break;
946  case eIMEXdirk_2_2_2:
947  {
948  m_numsteps = 1;
949  m_numstages = 3;
950 
951  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
952  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
953 
954  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
955  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
956  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
957  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
958  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
959  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
960 
961  NekDouble glambda = 0.2928932188134524756;
962  NekDouble gdelta = -0.7071067811865475244;
963 
964  m_A[0][1][1] = glambda;
965  m_A[0][2][1] = 1.0 - glambda;
966  m_A[0][2][2] = glambda;
967 
968  m_B[0][0][1] = 1.0 - glambda;
969  m_B[0][0][2] = glambda;
970 
971  m_A[1][1][0] = glambda;
972  m_A[1][2][0] = gdelta;
973  m_A[1][2][1] = 1.0 - gdelta;
974 
975  m_B[1][0][0] = gdelta;
976  m_B[1][0][1] = 1.0 - gdelta;
977 
981  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
982  m_timeLevelOffset[0] = 0;
983  }
984  break;
985  case eIMEXdirk_2_3_3:
986  {
987  m_numsteps = 1;
988  m_numstages = 3;
989 
990  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
991  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
992 
993  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
994  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
995  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
996  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
997  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
998  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
999 
1000  NekDouble glambda = 0.78867513459481288226;
1001 
1002  m_A[0][1][1] = glambda;
1003  m_A[0][2][1] = 1.0 - 2.0*glambda;
1004  m_A[0][2][2] = glambda;
1005 
1006  m_B[0][0][1] = 0.5;
1007  m_B[0][0][2] = 0.5;
1008 
1009  m_A[1][1][0] = glambda;
1010  m_A[1][2][0] = glambda - 1.0;
1011  m_A[1][2][1] = 2.0*(1-glambda);
1012 
1013  m_B[1][0][1] = 0.5;
1014  m_B[1][0][2] = 0.5;
1015 
1016  m_schemeType = eIMEX;
1019  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1020  m_timeLevelOffset[0] = 0;
1021  }
1022  break;
1023  case eIMEXdirk_1_1_1:
1024  {
1025  m_numsteps = 1;
1026  m_numstages = 2;
1027 
1028  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
1029  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
1030 
1031  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1032  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1033  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1034  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1035  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
1036  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
1037 
1038  m_A[0][1][1] = 1.0;
1039 
1040  m_B[0][0][1] = 1.0;
1041 
1042  m_A[1][1][0] = 1.0;
1043 
1044  m_B[1][0][0] = 1.0;
1045 
1046  m_schemeType = eIMEX;
1049  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1050  m_timeLevelOffset[0] = 0;
1051  }
1052  break;
1053  case eIMEXdirk_1_2_1:
1054  {
1055  m_numsteps = 1;
1056  m_numstages = 2;
1057 
1058  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
1059  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
1060 
1061  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1062  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1063  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1064  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1065  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
1066  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
1067 
1068  m_A[0][1][1] = 1.0;
1069 
1070  m_B[0][0][1] = 1.0;
1071 
1072  m_A[1][1][0] = 1.0;
1073 
1074  m_B[1][0][1] = 1.0;
1075 
1076  m_schemeType = eIMEX;
1079  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1080  m_timeLevelOffset[0] = 0;
1081  }
1082  break;
1083  case eIMEXdirk_1_2_2:
1084  {
1085  m_numsteps = 1;
1086  m_numstages = 2;
1087 
1088  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
1089  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
1090 
1091  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1092  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1093  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1094  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1095  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
1096  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
1097 
1098  m_A[0][1][1] = 1.0/2.0;
1099 
1100  m_B[0][0][1] = 1.0;
1101 
1102  m_A[1][1][0] = 1.0/2.0;
1103 
1104  m_B[1][0][1] = 1.0;
1105 
1106  m_schemeType = eIMEX;
1109  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1110  m_timeLevelOffset[0] = 0;
1111  }
1112  break;
1113  case eIMEXdirk_4_4_3:
1114  {
1115  m_numsteps = 1;
1116  m_numstages = 5;
1117 
1118  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
1119  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
1120 
1121  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1122  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1123  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1124  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1125  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
1126  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
1127 
1128  m_A[0][1][1] = 1.0/2.0;
1129  m_A[0][2][1] = 1.0/6.0;
1130  m_A[0][2][2] = 1.0/2.0;
1131  m_A[0][3][1] = -1.0/2.0;
1132  m_A[0][3][2] = 1.0/2.0;
1133  m_A[0][3][3] = 1.0/2.0;
1134  m_A[0][4][1] = 3.0/2.0;
1135  m_A[0][4][2] = -3.0/2.0;
1136  m_A[0][4][3] = 1.0/2.0;
1137  m_A[0][4][4] = 1.0/2.0;
1138 
1139  m_B[0][0][1] = 3.0/2.0;
1140  m_B[0][0][2] = -3.0/2.0;
1141  m_B[0][0][3] = 1.0/2.0;
1142  m_B[0][0][4] = 1.0/2.0;
1143 
1144  m_A[1][1][0] = 1.0/2.0;
1145  m_A[1][2][0] = 11.0/18.0;
1146  m_A[1][2][1] = 1.0/18.0;
1147  m_A[1][3][0] = 5.0/6.0;
1148  m_A[1][3][1] = -5.0/6.0;
1149  m_A[1][3][2] = 1.0/2.0;
1150  m_A[1][4][0] = 1.0/4.0;
1151  m_A[1][4][1] = 7.0/4.0;
1152  m_A[1][4][2] = 3.0/4.0;
1153  m_A[1][4][3] = -7.0/4.0;
1154 
1155  m_B[1][0][0] = 1.0/4.0;
1156  m_B[1][0][1] = 7.0/4.0;
1157  m_B[1][0][2] = 3.0/4.0;
1158  m_B[1][0][3] = -7.0/4.0;
1159 
1160  m_schemeType = eIMEX;
1163  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1164  m_timeLevelOffset[0] = 0;
1165  }
1166  break;
1167  default:
1168  {
1169  NEKERROR(ErrorUtil::efatal,"Invalid Time Integration Scheme");
1170  }
1171  }
1172 
1175 
1177  "Time integration scheme coefficients do not match its type");
1178  }
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
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:209
BDF multi-step scheme of order 1 (implicit)
Adams-Bashforth Forward multi-step scheme of order 2.
Array< OneD, Array< TwoD, NekDouble > > m_B
Runge-Kutta multi-stage scheme 4th order explicit (old name)
Implicit-Explicit Midpoint IMEX DIRK(1,2,2)
Classical RungeKutta4 method (new name for eClassicalRungeKutta4)
L-stable, four stage, third order IMEX DIRK(4,4,3)
Forward-Backward Euler IMEX DIRK(1,2,1)
Implicit Explicit General Linear Method.
Nonlinear SSP RungeKutta3 explicit.
Adams-Moulton Forward multi-step scheme of order 2.
Adams-Bashforth Forward multi-step scheme of order 3.
Crank-Nicolson/Adams-Bashforth Order 2 (CNAB)
IMEX 2nd order scheme using Backward Different Formula & Extrapolation.
Diagonally implicit scheme (e.g. the DIRK schemes)
Classical RungeKutta2 method (new name for eMidpoint)
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
Adams-Moulton Forward multi-step scheme of order 1.
Adams-Bashforth Forward multi-step scheme of order 1.
L-stable, three stage, third order IMEX DIRK(3,4,3)
Nonlinear SSP RungeKutta2 explicit (surrogate for eRungeKutta2_ImprovedEuler)
Improved RungeKutta2 explicit (old name meaning Heun&#39;s method)
Forward-Backward Euler IMEX DIRK(1,1,1)
double NekDouble
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
BDF multi-step scheme of order 2 (implicit)
L-stable, two stage, second order IMEX DIRK(2,2,2)
Diagonally Implicit Runge Kutta scheme of order 2.
IMEX 3rd order scheme using Backward Different Formula & Extrapolation.
L-stable, three stage, third order IMEX DIRK(3,4,3)
Diagonally Implicit Runge Kutta scheme of order 3.
L-stable, two stage, third order IMEX DIRK(2,3,3)
IMEX 4th order scheme using Backward Different Formula & Extrapolation.
Array< OneD, Array< TwoD, NekDouble > > m_A
Modified Crank-Nicolson/Adams-Bashforth Order 2 (MCNAB)
IMEX 1st order scheme using Euler Backwards/Euler Forwards.
midpoint method (old name)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
Adams-Bashforth Forward multi-step scheme of order 4.

◆ TimeIntegrationScheme() [2/3]

Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrationScheme ( )
inlineprivate

Definition at line 577 of file TimeIntegrationScheme.h.

References Nektar::ErrorUtil::efatal, and NEKERROR.

578  {
579  NEKERROR(ErrorUtil::efatal,"Default Constructor for the TimeIntegrationScheme class should not be called");
580  }
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:209
static const TimeIntegrationSchemeKey NullTimeIntegrationSchemeKey(eNoTimeIntegrationMethod)

◆ TimeIntegrationScheme() [3/3]

Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrationScheme ( const TimeIntegrationScheme in)
inlineprivate

Definition at line 582 of file TimeIntegrationScheme.h.

References Nektar::ErrorUtil::efatal, LIB_UTILITIES_EXPORT, and NEKERROR.

583  {
584  boost::ignore_unused(in);
585  NEKERROR(ErrorUtil::efatal,"Copy Constructor for the TimeIntegrationScheme class should not be called");
586  }
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:209
static const TimeIntegrationSchemeKey NullTimeIntegrationSchemeKey(eNoTimeIntegrationMethod)

Member Function Documentation

◆ A()

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

Definition at line 412 of file TimeIntegrationScheme.h.

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

413  {
414  return m_A[0][i][j];
415  }
Array< OneD, Array< TwoD, NekDouble > > m_A

◆ A_IMEX()

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

Definition at line 432 of file TimeIntegrationScheme.h.

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

433  {
434  return m_A[1][i][j];
435  }
Array< OneD, Array< TwoD, NekDouble > > m_A

◆ B()

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

◆ B_IMEX()

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

Definition at line 437 of file TimeIntegrationScheme.h.

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

438  {
439  return m_B[1][i][j];
440  }
Array< OneD, Array< TwoD, NekDouble > > m_B

◆ CheckIfFirstStageEqualsOldSolution()

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 1833 of file TimeIntegrationScheme.cpp.

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

Referenced by TimeIntegrationScheme().

1837  {
1838  boost::ignore_unused(B, V);
1839 
1840  int i,m;
1841  // First stage equals old solution if:
1842  // 1. the first row of the coefficient matrix A consists of zeros
1843  // 2. U[0][0] is equal to one and all other first row entries of U are zero
1844 
1845  // 1. Check first condition
1846  for(m = 0; m < A.num_elements(); m++)
1847  {
1848  for(i = 0; i < m_numstages; i++)
1849  {
1850  if( fabs(A[m][0][i]) > NekConstants::kNekZeroTol )
1851  {
1852  return false;
1853  }
1854  }
1855  }
1856 
1857  // 2. Check second condition
1858  if( fabs(U[0][0] - 1.0) > NekConstants::kNekZeroTol )
1859  {
1860  return false;
1861  }
1862  for(i = 1; i < m_numsteps; i++)
1863  {
1864  if( fabs(U[0][i]) > NekConstants::kNekZeroTol )
1865  {
1866  return false;
1867  }
1868  }
1869 
1870  return true;
1871  }
NekDouble B(const unsigned int i, const unsigned int j) const
NekDouble V(const unsigned int i, const unsigned int j) const
static const NekDouble kNekZeroTol
NekDouble U(const unsigned int i, const unsigned int j) const

◆ CheckIfLastStageEqualsNewSolution()

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 1873 of file TimeIntegrationScheme.cpp.

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

Referenced by TimeIntegrationScheme().

1877  {
1878  int i,m;
1879  // Last stage equals new solution if:
1880  // 1. the last row of the coefficient matrix A is equal to the first row of matrix B
1881  // 2. the last row of the coefficient matrix U is equal to the first row of matrix V
1882 
1883  // 1. Check first condition
1884  for(m = 0; m < A.num_elements(); m++)
1885  {
1886  for(i = 0; i < m_numstages; i++)
1887  {
1888  if( fabs(A[m][m_numstages-1][i]-B[m][0][i]) > NekConstants::kNekZeroTol )
1889  {
1890  return false;
1891  }
1892  }
1893  }
1894 
1895  // 2. Check second condition
1896  for(i = 0; i < m_numsteps; i++)
1897  {
1898  if( fabs(U[m_numstages-1][i]-V[0][i]) > NekConstants::kNekZeroTol )
1899  {
1900  return false;
1901  }
1902  }
1903 
1904  return true;
1905  }
NekDouble B(const unsigned int i, const unsigned int j) const
NekDouble V(const unsigned int i, const unsigned int j) const
static const NekDouble kNekZeroTol
NekDouble U(const unsigned int i, const unsigned int j) const

◆ CheckTimeIntegrateArguments()

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 1907 of file TimeIntegrationScheme.cpp.

References ASSERTL1, and m_numsteps.

Referenced by TimeIntegrate().

1913  {
1914  boost::ignore_unused(timestep, y_old, t_old, y_new, t_new, op);
1915 
1916  // Check if arrays are all of consistent size
1917  ASSERTL1(y_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1918  ASSERTL1(y_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1919 
1920  ASSERTL1(y_old[0]. num_elements()==y_new[0]. num_elements(),"Non-matching number of variables.");
1921  ASSERTL1(y_old[0][0].num_elements()==y_new[0][0].num_elements(),"Non-matching number of coefficients.");
1922 
1923  ASSERTL1(t_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1924  ASSERTL1(t_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1925 
1926  return true;
1927  }
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250

◆ Create()

TimeIntegrationSchemeSharedPtr Nektar::LibUtilities::TimeIntegrationScheme::Create ( const TimeIntegrationSchemeKey key)
staticprivate

Definition at line 142 of file TimeIntegrationScheme.cpp.

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

143  {
146  return returnval;
147  }
std::shared_ptr< TimeIntegrationScheme > TimeIntegrationSchemeSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.

◆ GetFirstDim()

int Nektar::LibUtilities::TimeIntegrationScheme::GetFirstDim ( ConstTripleArray y) const
inlineprivate

Definition at line 602 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

603  {
604  return y[0].num_elements();
605  }

◆ GetIntegrationMethod()

TimeIntegrationMethod Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationMethod ( ) const
inline

Definition at line 402 of file TimeIntegrationScheme.h.

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

403  {
405  }

◆ GetIntegrationSchemeKey()

const TimeIntegrationSchemeKey& Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeKey ( ) const
inline

Definition at line 397 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

398  {
399  return m_schemeKey;
400  }

◆ GetIntegrationSchemeType()

TimeIntegrationSchemeType Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeType ( ) const
inline

Definition at line 407 of file TimeIntegrationScheme.h.

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

408  {
409  return m_schemeType;
410  }

◆ GetNmultiStepDerivs()

unsigned int Nektar::LibUtilities::TimeIntegrationScheme::GetNmultiStepDerivs ( void  ) const
inline

Definition at line 457 of file TimeIntegrationScheme.h.

References LIB_UTILITIES_EXPORT.

458  {
459  return m_numMultiStepDerivs;
460  }

◆ GetNmultiStepValues()

unsigned int Nektar::LibUtilities::TimeIntegrationScheme::GetNmultiStepValues ( void  ) const
inline

Definition at line 452 of file TimeIntegrationScheme.h.

453  {
454  return m_numMultiStepValues;
455  }

◆ GetNstages()

unsigned int Nektar::LibUtilities::TimeIntegrationScheme::GetNstages ( void  ) const
inline

Definition at line 447 of file TimeIntegrationScheme.h.

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

448  {
449  return m_numstages;
450  }

◆ GetNsteps()

unsigned int Nektar::LibUtilities::TimeIntegrationScheme::GetNsteps ( void  ) const
inline

Definition at line 442 of file TimeIntegrationScheme.h.

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

443  {
444  return m_numsteps;
445  }

◆ GetSecondDim()

int Nektar::LibUtilities::TimeIntegrationScheme::GetSecondDim ( ConstTripleArray y) const
inlineprivate

Definition at line 607 of file TimeIntegrationScheme.h.

References LIB_UTILITIES_EXPORT.

Referenced by TimeIntegrate().

608  {
609  return y[0][0].num_elements();
610  }

◆ GetTimeLevelOffset()

const Array<OneD, const unsigned int>& Nektar::LibUtilities::TimeIntegrationScheme::GetTimeLevelOffset ( )
inline

Definition at line 522 of file TimeIntegrationScheme.h.

523  {
524  return m_timeLevelOffset;
525  }

◆ InitializeScheme()

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 1239 of file TimeIntegrationScheme.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoProjection(), Blas::Dscal(), Nektar::LibUtilities::eExplicit, GetIntegrationSchemeType(), m_numMultiStepDerivs, m_numMultiStepValues, m_schemeKey, and m_timeLevelOffset.

1243  {
1244  // create a TimeIntegrationSolution object based upon the
1245  // initial value. Initialise all other multi-step values
1246  // and derivatives to zero
1249 
1251  {
1252  // ensure initial solution is in correct space
1253  op.DoProjection(y_0,y_out->UpdateSolution(),time);
1254  }
1255 
1256  // calculate the initial derivative, if is part of the
1257  // solution vector of the current scheme
1259  {
1261  {
1262  int i;
1263  int nvar = y_0.num_elements();
1264  int npoints = y_0[0].num_elements();
1265  DoubleArray f_y_0(nvar);
1266  for(i = 0; i < nvar; i++)
1267  {
1268  f_y_0[i] = Array<OneD,NekDouble>(npoints);
1269  }
1270  // calculate the derivative of the initial value
1271  op.DoOdeRhs(y_0,f_y_0,time);
1272 
1273  // multiply by the step size
1274  for(i = 0; i < nvar; i++)
1275  {
1276  Blas::Dscal(npoints,timestep,f_y_0[i].get(),1);
1277  }
1278  y_out->SetDerivative(0,f_y_0,timestep);
1279  }
1280  }
1281 
1282  return y_out;
1283  }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Array< OneD, Array< OneD, NekDouble > > DoubleArray
TimeIntegrationSchemeType GetIntegrationSchemeType() const
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:125
std::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr

◆ TimeIntegrate() [1/2]

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 1286 of file TimeIntegrationScheme.cpp.

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

1289  {
1291  "Fully Implicit integration scheme cannot be handled by this routine.");
1292 
1293  int nvar = solvector->GetFirstDim ();
1294  int npoints = solvector->GetSecondDim();
1295 
1296  if( (solvector->GetIntegrationScheme()).get() != this )
1297  {
1298  // This branch will be taken when the solution vector
1299  // (solvector) is set up for a different scheme than
1300  // the object this method is called from. (typically
1301  // needed to calculate the first time-levels of a
1302  // multi-step scheme)
1303 
1304  // To do this kind of 'non-matching' integration, we
1305  // perform the following three steps:
1306  //
1307  // 1: copy the required input information from the
1308  // solution vector of the master scheme to the
1309  // input solution vector of the current scheme
1310  //
1311  // 2: time-integrate for one step using the current
1312  // scheme
1313  //
1314  // 3: copy the information contained in the output
1315  // vector of the current scheme to the solution
1316  // vector of the master scheme
1317 
1318  // STEP 1: copy the required input information from
1319  // the solution vector of the master scheme
1320  // to the input solution vector of the
1321  // current scheme
1322 
1323  // 1.1 Determine which information is required for the
1324  // current scheme
1325  int n;
1326  DoubleArray y_n;
1327  NekDouble t_n = 0;
1328  DoubleArray dtFy_n;
1329  unsigned int nCurSchemeVals = m_numMultiStepValues; // number of required values of the current scheme
1330  unsigned int nCurSchemeDers = m_numMultiStepDerivs; // number of required derivs of the current scheme
1331  unsigned int nCurSchemeSteps = m_numsteps; // number of steps in the current scheme
1332  unsigned int nMasterSchemeVals = solvector->GetNvalues(); // number of values of the master scheme
1333  unsigned int nMasterSchemeDers = solvector->GetNderivs(); // number of derivs of the master scheme
1334  // The arrays below contains information to which
1335  // time-level the values and derivatives of the
1336  // schemes belong
1337  const Array<OneD, const unsigned int>& curTimeLevels = m_timeLevelOffset;
1338  const Array<OneD, const unsigned int>& masterTimeLevels = solvector->GetTimeLevelOffset();
1339 
1340  // 1.2 Copy the required information from the master
1341  // solution vector to the input solution vector of
1342  // the current scheme
1344  AllocateSharedPtr(GetIntegrationSchemeKey()); // input solution vector of the current scheme
1345 
1346  for(n = 0; n < nCurSchemeVals; n++)
1347  {
1348  // Get the required value out of the master solution vector
1349  //DoubleArray& y_n = solvector->GetValue ( curTimeLevels[n] );
1350  //NekDouble t_n = solvector->GetValueTime( curTimeLevels[n] );
1351 
1352  y_n = solvector->GetValue ( curTimeLevels[n] );
1353  t_n = solvector->GetValueTime( curTimeLevels[n] );
1354 
1355  // Set the required value in the input solution
1356  // vector of the current scheme
1357  solvector_in->SetValue(curTimeLevels[n],y_n,t_n);
1358  }
1359  for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
1360  {
1361  // Get the required derivative out of the master
1362  // solution vector
1363  //DoubleArray& dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
1364  dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
1365 
1366  // Set the required derivative in the input
1367  // solution vector of the current scheme
1368  solvector_in->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
1369  }
1370 
1371  // STEP 2: time-integrate for one step using the
1372  // current scheme
1373  TimeIntegrationSolutionSharedPtr solvector_out = MemoryManager<TimeIntegrationSolution>:: AllocateSharedPtr(GetIntegrationSchemeKey(),nvar,npoints); // output solution vector of the current scheme
1374 
1375  // integrate
1376  TimeIntegrate(timestep, solvector_in->GetSolutionVector(),
1377  solvector_in->GetTimeVector(),
1378  solvector_out->UpdateSolutionVector(),
1379  solvector_out->UpdateTimeVector(),op);
1380 
1381 
1382  // STEP 3: copy the information contained in the
1383  // output vector of the current scheme to the
1384  // solution vector of the master scheme
1385 
1386  // 3.1 Check whether the current time scheme updates
1387  // the most recent derivative that should be
1388  // updated in the master scheme. If not,
1389  // calculate the derivative. This can be done
1390  // based upon the corresponding value and the
1391  // DoOdeRhs operator.
1392  int j;
1393  bool CalcNewDeriv = false; // flag inidicating whether the new derivative is availble in the output of
1394  // of the current scheme or whether it should be calculated
1395  if( nMasterSchemeDers > 0 )
1396  {
1397  if(nCurSchemeDers == 0)
1398  {
1399  CalcNewDeriv = true;
1400  }
1401  else
1402  {
1403  if( masterTimeLevels[nMasterSchemeVals] < curTimeLevels[nCurSchemeVals] )
1404  {
1405  CalcNewDeriv = true;
1406  }
1407  }
1408  }
1409 
1410  if(CalcNewDeriv)
1411  {
1412  int newDerivTimeLevel = masterTimeLevels[nMasterSchemeVals]; // contains the time level at which
1413  // we want to know the derivative of the
1414  // master scheme
1415  //DoubleArray y_n;
1416  //NekDouble t_n;
1417  // if the time level correspond to 0, calculate the derivative based upon the solution value
1418  // at the new time-level
1419  if (newDerivTimeLevel == 0)
1420  {
1421  y_n = solvector_out->GetValue(0);
1422  t_n = solvector_out->GetValueTime(0);
1423  }
1424  // if the time level correspond to 1, calculate the derivative based upon the solution value
1425  // at the new old-level
1426  else if( newDerivTimeLevel == 1 )
1427  {
1428  y_n = solvector->GetValue(0);
1429  t_n = solvector->GetValueTime(0);
1430  }
1431  else
1432  {
1433  ASSERTL1(false,"Problems with initialising scheme");
1434  }
1435 
1436  DoubleArray f_n(nvar);
1437  for(j = 0; j < nvar; j++)
1438  {
1439  f_n[j] = Array<OneD, NekDouble>(npoints);
1440  }
1441 
1442  // calculate the derivative
1443  op.DoOdeRhs(y_n, f_n, t_n);
1444 
1445  // multiply by dt (as required by the General Linear Method framework)
1446  for(j = 0; j < nvar; j++)
1447  {
1448  Vmath::Smul(npoints,timestep,f_n[j],1,
1449  f_n[j],1);
1450  }
1451 
1452  // Rotate the solution vector
1453  // (i.e. updating without calculating/inserting new values)
1454  solvector->RotateSolutionVector();
1455  // Set the calculated derivative in the master solution vector
1456  solvector->SetDerivative(newDerivTimeLevel,f_n,timestep);
1457  }
1458  else
1459  {
1460  // Rotate the solution vector (i.e. updating
1461  // without calculating/inserting new values)
1462  solvector->RotateSolutionVector();
1463  }
1464 
1465 
1466  // 1.2 Copy the information calculated using the
1467  // current scheme from the output solution vector
1468  // to the master solution vector
1469  for(n = 0; n < nCurSchemeVals; n++)
1470  {
1471  // Get the calculated value out of the output
1472  // solution vector of the current scheme
1473  //DoubleArray& y_n = solvector_out->GetValue ( curTimeLevels[n] );
1474  //NekDouble t_n = solvector_out->GetValueTime( curTimeLevels[n] );
1475  y_n = solvector_out->GetValue ( curTimeLevels[n] );
1476  t_n = solvector_out->GetValueTime( curTimeLevels[n] );
1477 
1478  // Set the calculated value in the master solution vector
1479  solvector->SetValue(curTimeLevels[n],y_n,t_n);
1480  }
1481 
1482  for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
1483  {
1484  // Get the calculated derivative out of the output
1485  // solution vector of the current scheme
1486  // DoubleArray& dtFy_n =
1487  // solvector_out->GetDerivative (curTimeLevels[n]);
1488  dtFy_n = solvector_out->GetDerivative ( curTimeLevels[n] );
1489 
1490  // Set the calculated derivative in the master
1491  // solution vector
1492  solvector->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
1493  }
1494  }
1495  else
1496  {
1497  const TimeIntegrationSchemeKey& key = solvector->GetIntegrationSchemeKey();
1498 
1500 
1501  TimeIntegrate(timestep,solvector->GetSolutionVector(),
1502  solvector->GetTimeVector(),
1503  solvector_new->UpdateSolutionVector(),
1504  solvector_new->UpdateTimeVector(),op);
1505 
1506  solvector = solvector_new;
1507  }
1508  return solvector->GetSolution();
1509  }
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
ConstDoubleArray & TimeIntegrate(const NekDouble timestep, TimeIntegrationSolutionSharedPtr &y, const TimeIntegrationSchemeOperators &op)
Explicit integration of an ODE.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Array< OneD, Array< OneD, NekDouble > > DoubleArray
double NekDouble
TimeIntegrationSchemeType GetIntegrationSchemeType() const
const TimeIntegrationSchemeKey & GetIntegrationSchemeKey() const
std::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250

◆ TimeIntegrate() [2/2]

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

1517  {
1518  ASSERTL1(CheckTimeIntegrateArguments(timestep,y_old,t_old,y_new,t_new,op), "Arguments not well defined");
1519 
1520  unsigned int i,j,k;
1522 
1523  // Check if storage has already been initialised.
1524  // If so, we just zero the temporary storage.
1525  if (m_initialised && m_nvar == GetFirstDim(y_old)
1526  && m_npoints == GetSecondDim(y_old))
1527  {
1528  for(j = 0; j < m_nvar; j++)
1529  {
1530  Vmath::Zero(m_npoints, m_tmp[j], 1);
1531  }
1532  }
1533  else
1534  {
1535  m_nvar = GetFirstDim(y_old);
1536  m_npoints = GetSecondDim(y_old);
1537 
1538  // First, we are going to calculate the various stage
1539  // values and stage derivatives (this is the multi-stage
1540  // part of the method)
1541  // - m_Y corresponds to the stage values
1542  // - m_F corresponds to the stage derivatives
1543  // - m_T corresponds to the time at the different stages
1544  // - m_tmp corresponds to the explicit right hand side of
1545  // each stage equation
1546  // (for explicit schemes, this correspond to m_Y)
1547 
1548  // Allocate memory for the arrays m_Y and m_F and m_tmp The same
1549  // storage will be used for every stage -> m_Y is a
1550  // DoubleArray
1552  for(j = 0; j < m_nvar; j++)
1553  {
1554  m_tmp[j] = Array<OneD, NekDouble>(m_npoints,0.0);
1555  }
1556 
1557  // The same storage will be used for every stage -> m_tmp is
1558  // a DoubleArray
1559  if(type == eExplicit)
1560  {
1561  m_Y = m_tmp;
1562  }
1563  else
1564  {
1565  m_Y = DoubleArray(m_nvar);
1566  for(j = 0; j < m_nvar; j++)
1567  {
1568  m_Y[j] = Array<OneD, NekDouble>(m_npoints,0.0);
1569  }
1570  }
1571 
1572  // Different storage for every stage derivative as the data
1573  // will be re-used to update the solution -> m_F is a TripleArray
1575  for(i = 0; i < m_numstages; ++i)
1576  {
1577  m_F[i] = DoubleArray(m_nvar);
1578  for(j = 0; j < m_nvar; j++)
1579  {
1580  m_F[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
1581  }
1582  }
1583 
1584  if(type == eIMEX)
1585  {
1586  m_F_IMEX = TripleArray(m_numstages);
1587  for(i = 0; i < m_numstages; ++i)
1588  {
1589  m_F_IMEX[i] = DoubleArray(m_nvar);
1590  for(j = 0; j < m_nvar; j++)
1591  {
1592  m_F_IMEX[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
1593  }
1594  }
1595  }
1596 
1597  // Finally, flag that we have initialised the memory.
1598  m_initialised = true;
1599  }
1600 
1601  // The loop below calculates the stage values and derivatives
1602  for(i = 0; i < m_numstages; i++)
1603  {
1604  if( (i==0) && m_firstStageEqualsOldSolution )
1605  {
1606  for(k = 0; k < m_nvar; k++)
1607  {
1608  Vmath::Vcopy(m_npoints,y_old[0][k],1,m_Y[k],1);
1609  }
1610  m_T = t_old[0];
1611  }
1612  else
1613  {
1614  // The stage values m_Y are a linear combination of:
1615  // 1: the stage derivatives
1616 
1617  if( i != 0 )
1618  {
1619  for(k = 0; k < m_nvar; k++)
1620  {
1621  Vmath::Smul(m_npoints,timestep*A(i,0),m_F[0][k],1,
1622  m_tmp[k],1);
1623 
1624  if(type == eIMEX)
1625  {
1626  Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,0),
1627  m_F_IMEX[0][k],1,
1628  m_tmp[k],1,m_tmp[k],1);
1629  }
1630  }
1631  }
1632  m_T = A(i,0)*timestep;
1633 
1634  for( j = 1; j < i; j++ )
1635  {
1636  for(k = 0; k < m_nvar; k++)
1637  {
1638  Vmath::Svtvp(m_npoints,timestep*A(i,j),m_F[j][k],1,
1639  m_tmp[k],1,m_tmp[k],1);
1640  if(type == eIMEX)
1641  {
1642  Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,j),
1643  m_F_IMEX[j][k],1,
1644  m_tmp[k],1,m_tmp[k],1);
1645  }
1646  }
1647 
1648  m_T += A(i,j)*timestep;
1649  }
1650 
1651  // 2: the imported multi-step solution of the
1652  // previous time level
1653  for(j = 0; j < m_numsteps; j++)
1654  {
1655  for(k = 0; k < m_nvar; k++)
1656  {
1657  Vmath::Svtvp(m_npoints,U(i,j),y_old[j][k],1,
1658  m_tmp[k],1,m_tmp[k],1);
1659  }
1660  m_T += U(i,j)*t_old[j];
1661  }
1662  }
1663 
1664  // Calculate the stage derivative based upon the stage value
1665  if(type == eDiagonallyImplicit)
1666  {
1667  if(m_numstages==1)
1668  {
1669  m_T= t_old[0]+timestep;
1670  }
1671  else
1672  {
1673  m_T= t_old[0];
1674  for(int j=0; j<=i; ++j)
1675  {
1676  m_T += A(i,j)*timestep;
1677  }
1678  }
1679 
1680  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
1681 
1682  for(k = 0; k < m_nvar; k++)
1683  {
1684  Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
1685  Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),m_F[i][k],1,m_F[i][k],1);
1686  }
1687  }
1688  else if(type == eIMEX)
1689  {
1690  if(m_numstages==1)
1691  {
1692  m_T= t_old[0]+timestep;
1693  }
1694  else
1695  {
1696  m_T= t_old[0];
1697  for(int j=0; j<=i; ++j)
1698  {
1699  m_T += A(i,j)*timestep;
1700  }
1701  }
1702 
1703  if(fabs(A(i,i)) > NekConstants::kNekZeroTol)
1704  {
1705  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
1706 
1707  for(k = 0; k < m_nvar; k++)
1708  {
1709  Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
1710  Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),
1711  m_F[i][k],1,m_F[i][k],1);
1712  }
1713  }
1714  op.DoOdeRhs(m_Y, m_F_IMEX[i], m_T);
1715  }
1716  else if( type == eExplicit)
1717  {
1718  // Avoid projecting the same solution twice
1719  if( ! ((i==0) && m_firstStageEqualsOldSolution) )
1720  {
1721  // ensure solution is in correct space
1722  op.DoProjection(m_Y,m_Y,m_T);
1723  }
1724  op.DoOdeRhs(m_Y, m_F[i], m_T);
1725  }
1726  }
1727 
1728  // Next, the solution vector y at the new time level will
1729  // be calculated.
1730  //
1731  // For multi-step methods, this includes updating the
1732  // values of the auxiliary parameters
1733  //
1734  // The loop below calculates the solution at the new time
1735  // level
1736  //
1737  // If last stage equals the new solution, the new solution
1738  // needs not be calculated explicitly but can simply be
1739  // copied. This saves a solve.
1740  int i_start = 0;
1742  {
1743  for(k = 0; k < m_nvar; k++)
1744  {
1745  Vmath::Vcopy(m_npoints,m_Y[k],1,y_new[0][k],1);
1746  }
1747 
1748  if (m_numstages==1 && type == eIMEX)
1749  {
1750  t_new[0] = t_old[0]+timestep;
1751  }
1752  else
1753  {
1754  t_new[0] = B(0,0)*timestep;
1755  for(j = 1; j < m_numstages; j++)
1756  {
1757  t_new[0] += B(0,j)*timestep;
1758  }
1759  for(j = 0; j < m_numsteps; j++)
1760  {
1761  t_new[0] += V(0,j)*t_old[j];
1762  }
1763  }
1764  i_start = 1;
1765  }
1766 
1767  for(i = i_start; i < m_numsteps; i++)
1768  {
1769  // The solution at the new time level is a linear
1770  // combination of:
1771  // 1: the stage derivatives
1772  for(k = 0; k < m_nvar; k++)
1773  {
1774  Vmath::Smul(m_npoints,timestep*B(i,0),m_F[0][k],1,
1775  y_new[i][k],1);
1776 
1777  if(type == eIMEX)
1778  {
1779  Vmath::Svtvp(m_npoints,timestep*B_IMEX(i,0),
1780  m_F_IMEX[0][k],1,y_new[i][k],1,
1781  y_new[i][k],1);
1782  }
1783  }
1784  if(m_numstages != 1 || type != eIMEX)
1785  {
1786  t_new[i] = B(i,0)*timestep;
1787  }
1788 
1789 
1790  for(j = 1; j < m_numstages; j++)
1791  {
1792  for(k = 0; k < m_nvar; k++)
1793  {
1794  Vmath::Svtvp(m_npoints,timestep*B(i,j),m_F[j][k],1,
1795  y_new[i][k],1,y_new[i][k],1);
1796 
1797  if(type == eIMEX)
1798  {
1799  Vmath::Svtvp(m_npoints,timestep*B_IMEX(i,j),
1800  m_F_IMEX[j][k],1,y_new[i][k],1,
1801  y_new[i][k],1);
1802  }
1803  }
1804  if(m_numstages != 1 || type != eIMEX)
1805  {
1806  t_new[i] += B(i,j)*timestep;
1807  }
1808  }
1809 
1810  // 2: the imported multi-step solution of the previous
1811  // time level
1812  for(j = 0; j < m_numsteps; j++)
1813  {
1814  for(k = 0; k < m_nvar; k++)
1815  {
1816  Vmath::Svtvp(m_npoints,V(i,j),y_old[j][k],1,
1817  y_new[i][k],1,y_new[i][k],1);
1818  }
1819  if(m_numstages != 1 || type != eIMEX)
1820  {
1821  t_new[i] += V(i,j)*t_old[j];
1822  }
1823  }
1824  }
1825 
1826  // Ensure that the new solution is projected if necessary
1827  if(type == eExplicit)
1828  {
1829  op.DoProjection(y_new[0],y_new[0],t_new[0]);
1830  }
1831  }
NekDouble A(const unsigned int i, const unsigned int j) const
NekDouble B(const unsigned int i, const unsigned int j) const
int m_nvar
bool to identify if array has been initialised
bool CheckTimeIntegrateArguments(const NekDouble timestep, ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new, SingleArray &t_new, const TimeIntegrationSchemeOperators &op) const
DoubleArray m_tmp
Array containing the stage values.
NekDouble V(const unsigned int i, const unsigned int j) const
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
Definition: Vmath.cpp:488
Implicit Explicit General Linear Method.
DoubleArray m_Y
The size of inner data which is stored for reuse.
NekDouble A_IMEX(const unsigned int i, const unsigned int j) const
TripleArray m_F
explicit right hand side of each stage equation
Diagonally implicit scheme (e.g. the DIRK schemes)
static const NekDouble kNekZeroTol
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
NekDouble U(const unsigned int i, const unsigned int j) const
int m_npoints
The number of variables in integration scheme.
Array< OneD, Array< OneD, NekDouble > > DoubleArray
NekDouble B_IMEX(const unsigned int i, const unsigned int j) const
TimeIntegrationSchemeType GetIntegrationSchemeType() const
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.
Definition: Vmath.cpp:346
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376
NekDouble m_T
Used to store the Explicit stage derivative of IMEX schemes.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > TripleArray
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
TripleArray m_F_IMEX
Array corresponding to the stage Derivatives.

◆ U()

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

Definition at line 422 of file TimeIntegrationScheme.h.

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

423  {
424  return m_U[i][j];
425  }

◆ V()

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

Definition at line 427 of file TimeIntegrationScheme.h.

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

428  {
429  return m_V[i][j];
430  }

◆ VerifyIntegrationSchemeType()

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 1182 of file TimeIntegrationScheme.cpp.

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

Referenced by TimeIntegrationScheme().

1187  {
1188  boost::ignore_unused(B, U, V);
1189 
1190  int i;
1191  int j;
1192  int m;
1193  int IMEXdim = A.num_elements();
1194  int dim = A[0].GetRows();
1195 
1196  Array<OneD, TimeIntegrationSchemeType> vertype(IMEXdim,eExplicit);
1197 
1198  for(m = 0; m < IMEXdim; m++)
1199  {
1200  for(i = 0; i < dim; i++)
1201  {
1202  if( fabs(A[m][i][i]) > NekConstants::kNekZeroTol )
1203  {
1204  vertype[m] = eDiagonallyImplicit;
1205  }
1206  }
1207 
1208  for(i = 0; i < dim; i++)
1209  {
1210  for(j = i+1; j < dim; j++)
1211  {
1212  if( fabs(A[m][i][j]) > NekConstants::kNekZeroTol )
1213  {
1214  vertype[m] = eImplicit;
1215  ASSERTL1(false,"Fully Implicit schemes cannnot be handled by the TimeIntegrationScheme class");
1216  }
1217  }
1218  }
1219  }
1220 
1221  if(IMEXdim == 2)
1222  {
1223  ASSERTL1(B.num_elements()==2,"Coefficient Matrix B should have an implicit and explicit part for IMEX schemes");
1224  if((vertype[0] == eDiagonallyImplicit) &&
1225  (vertype[1] == eExplicit))
1226  {
1227  vertype[0] = eIMEX;
1228  }
1229  else
1230  {
1231  ASSERTL1(false,"This is not a proper IMEX scheme");
1232  }
1233  }
1234 
1235  return (vertype[0] == type);
1236  }
NekDouble B(const unsigned int i, const unsigned int j) const
NekDouble V(const unsigned int i, const unsigned int j) const
Implicit Explicit General Linear Method.
Diagonally implicit scheme (e.g. the DIRK schemes)
static const NekDouble kNekZeroTol
NekDouble U(const unsigned int i, const unsigned int j) const
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250

Friends And Related Function Documentation

◆ Nektar::MemoryManager

template<typename >
friend class Nektar::MemoryManager
friend

Time at the different stages.

Definition at line 569 of file TimeIntegrationScheme.h.

◆ TimeIntegrationSchemeManager

TimeIntegrationSchemeManagerT& TimeIntegrationSchemeManager ( void  )
friend

Definition at line 44 of file TimeIntegrationScheme.cpp.

45  {
46  static TimeIntegrationSchemeManagerT instance;
47  instance.RegisterGlobalCreator(TimeIntegrationScheme::Create);
48  return instance;
49  }
NekManager< TimeIntegrationSchemeKey, TimeIntegrationScheme, TimeIntegrationSchemeKey::opLess > TimeIntegrationSchemeManagerT
static std::shared_ptr< TimeIntegrationScheme > Create(const TimeIntegrationSchemeKey &key)

Member Data Documentation

◆ m_A

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

Definition at line 552 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme().

◆ m_B

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

Definition at line 553 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme().

◆ m_F

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

explicit right hand side of each stage equation

Definition at line 564 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

◆ m_F_IMEX

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

Array corresponding to the stage Derivatives.

Definition at line 565 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

◆ m_firstStageEqualsOldSolution

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

Definition at line 534 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate(), and TimeIntegrationScheme().

◆ m_initialised

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

Definition at line 558 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

◆ m_lastStageEqualsNewSolution

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

Definition at line 535 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate(), and TimeIntegrationScheme().

◆ m_npoints

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

The number of variables in integration scheme.

Definition at line 560 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

◆ m_numMultiStepDerivs

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

Definition at line 539 of file TimeIntegrationScheme.h.

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

◆ m_numMultiStepValues

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

Definition at line 537 of file TimeIntegrationScheme.h.

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

◆ m_numstages

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

◆ m_numsteps

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

◆ m_nvar

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

bool to identify if array has been initialised

Definition at line 559 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

◆ m_schemeKey

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

Definition at line 529 of file TimeIntegrationScheme.h.

Referenced by InitializeScheme().

◆ m_schemeType

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

Definition at line 530 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme().

◆ m_T

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

Used to store the Explicit stage derivative of IMEX schemes.

Definition at line 567 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

◆ m_timeLevelOffset

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

Definition at line 541 of file TimeIntegrationScheme.h.

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

◆ m_tmp

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

Array containing the stage values.

Definition at line 562 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

◆ m_U

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

Definition at line 554 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme().

◆ m_V

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

Definition at line 555 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme().

◆ m_Y

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

The size of inner data which is stored for reuse.

Definition at line 561 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().