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>

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

Public Types

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

Public Member Functions

virtual ~TimeIntegrationScheme ()
 
const TimeIntegrationSchemeKey & GetIntegrationSchemeKey () const
 
TimeIntegrationMethod GetIntegrationMethod () const
 
TimeIntegrationSchemeType GetIntegrationSchemeType () const
 
NekDouble A (const unsigned int i, const unsigned int j) const
 
NekDouble B (const unsigned int i, const unsigned int j) const
 
NekDouble U (const unsigned int i, const unsigned int j) const
 
NekDouble V (const unsigned int i, const unsigned int j) const
 
NekDouble A_IMEX (const unsigned int i, const unsigned int j) const
 
NekDouble B_IMEX (const unsigned int i, const unsigned int j) const
 
unsigned int GetNsteps (void) const
 
unsigned int GetNstages (void) const
 
unsigned int GetNmultiStepValues (void) const
 
unsigned int GetNmultiStepDerivs (void) const
 
TimeIntegrationSolutionSharedPtr InitializeScheme (const NekDouble timestep, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)
 This function initialises the time integration scheme. More...
 
ConstDoubleArray & TimeIntegrate (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, NekDouble > m_U
 
Array< TwoD, NekDouble > m_V
 

Private Member Functions

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

Static Private Member Functions

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

Private Attributes

bool m_initialised
 
int m_nvar
 bool to identify if array has been initialised 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...
 
TimeIntegrationSchemeManagerT & TimeIntegrationSchemeManager (void)
 

Detailed Description

Definition at line 355 of file TimeIntegrationScheme.h.

Member Typedef Documentation

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

Definition at line 361 of file TimeIntegrationScheme.h.

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

Definition at line 363 of file TimeIntegrationScheme.h.

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

Definition at line 359 of file TimeIntegrationScheme.h.

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

Definition at line 362 of file TimeIntegrationScheme.h.

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

Definition at line 365 of file TimeIntegrationScheme.h.

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

Definition at line 366 of file TimeIntegrationScheme.h.

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

Definition at line 364 of file TimeIntegrationScheme.h.

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

Definition at line 360 of file TimeIntegrationScheme.h.

Constructor & Destructor Documentation

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

Definition at line 370 of file TimeIntegrationScheme.h.

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

Definition at line 150 of file TimeIntegrationScheme.cpp.

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

150  :
151  m_schemeKey(key),
152  m_initialised(false)
153  {
154  switch(key.GetIntegrationMethod())
155  {
156  case eForwardEuler:
157  case eAdamsBashforthOrder1:
158  {
159  m_numsteps = 1;
160  m_numstages = 1;
161 
162  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
163  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
164 
165  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
166  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
167  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
168  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 1.0);
169 
170  m_schemeType = eExplicit;
171  m_numMultiStepValues = 1;
172  m_numMultiStepDerivs = 0;
173  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
174  m_timeLevelOffset[0] = 0;
175  }
176  break;
177  case eAdamsBashforthOrder2:
178  {
179  m_numsteps = 2;
180  m_numstages = 1;
181 
182  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
183  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
184 
185  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
186  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages);
187  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
188  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
189 
190  m_B[0][0][0] = 3.0/2.0;
191  m_B[0][1][0] = 1.0;
192 
193  m_U[0][0] = 1.0;
194 
195  m_V[0][0] = 1.0;
196  m_V[0][1] = -0.5;
197 
198  m_schemeType = eExplicit;
199  m_numMultiStepValues = 1;
200  m_numMultiStepDerivs = 1;
201  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
202  m_timeLevelOffset[0] = 0;
203  m_timeLevelOffset[1] = 1;
204  }
205  break;
206  case eAdamsBashforthOrder3:
207  {
208  m_numsteps = 4;
209  m_numstages = 1;
210 
211  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
212  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
213 
214  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
215  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
216  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
217  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
218 
219  m_B[0][1][0] = 1.0;
220 
221  m_U[0][0] = 1.0;
222  m_U[0][1] = 23.0/12.0;
223  m_U[0][2] = -4.0/3.0;
224  m_U[0][3] = 5.0/12.0;
225 
226  m_V[0][0] = 1.0;
227  m_V[0][1] = 23.0/12.0;
228  m_V[0][2] = -4.0/3.0;
229  m_V[0][3] = 5.0/12.0;
230  m_V[2][1] = 1.0;
231  m_V[3][2] = 1.0;
232 
233  m_schemeType = eExplicit;
234  m_numMultiStepValues = 1;
235  m_numMultiStepDerivs = 3;
236  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
237  m_timeLevelOffset[0] = 0;
238  m_timeLevelOffset[1] = 1;
239  m_timeLevelOffset[2] = 2;
240  m_timeLevelOffset[3] = 3;
241  }
242  break;
243  case eBackwardEuler:
244  case eBDFImplicitOrder1:
245  case eAdamsMoultonOrder1:
246  {
247  m_numsteps = 1;
248  m_numstages = 1;
249 
250  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
251  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
252 
253  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,1.0);
254  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
255  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
256  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 1.0);
257 
258  m_schemeType = eDiagonallyImplicit;
259  m_numMultiStepValues = 1;
260  m_numMultiStepDerivs = 0;
261  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
262  m_timeLevelOffset[0] = 0; // SJS: Not sure whether this is correct
263  }
264  break;
265  case eIMEXOrder1:
266  {
267  m_numsteps = 2;
268  m_numstages = 1;
269 
270  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
271  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
272 
273  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,1.0);
274  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
275  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
276  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
277  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
278  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
279 
280  m_B[0][1][0] = 0.0;
281  m_B[1][0][0] = 0.0;
282  m_V[0][0] = 1.0;
283  m_V[0][1] = 1.0;
284 
285  m_schemeType = eIMEX;
286  m_numMultiStepValues = 1;
287  m_numMultiStepDerivs = 1;
288  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
289  m_timeLevelOffset[0] = 0;
290  m_timeLevelOffset[1] = 0;
291  }
292  break;
293  case eIMEXOrder2:
294  {
295  NekDouble third = 1.0/3.0;
296  m_numsteps = 4;
297  m_numstages = 1;
298 
299  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
300  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
301 
302  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,2*third);
303  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
304  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
305  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
306  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 4*third);
307  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
308 
309  m_B[0][0][0] = 2*third;
310  m_B[1][2][0] = 1.0;
311  m_U[0][1] = -third;
312  m_U[0][3] = -2*third;
313 
314  m_V[0][0] = 4*third;
315  m_V[0][1] = -third;
316  m_V[0][2] = 4*third;
317  m_V[0][3] = -2*third;
318  m_V[1][0] = 1.0;
319  m_V[3][2] = 1.0;
320 
321  m_schemeType = eIMEX;
322  m_numMultiStepValues = 2;
323  m_numMultiStepDerivs = 2;
324  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
325  m_timeLevelOffset[0] = 0;
326  m_timeLevelOffset[1] = 1;
327  m_timeLevelOffset[2] = 0;
328  m_timeLevelOffset[3] = 1;
329  }
330  break;
331  case eIMEXOrder3:
332  {
333  NekDouble eleventh = 1.0/11.0;
334  m_numsteps = 6;
335  m_numstages = 1;
336 
337  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
338  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
339 
340  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,6*eleventh);
341  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
342  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
343  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
344  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 18*eleventh);
345  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
346 
347  m_B[0][0][0] = 6*eleventh;
348  m_B[1][3][0] = 1.0;
349  m_U[0][1] = -9*eleventh;
350  m_U[0][2] = 2*eleventh;
351  m_U[0][4] = -18*eleventh;
352  m_U[0][5] = 6*eleventh;
353 
354  m_V[0][0] = 18*eleventh;
355  m_V[0][1] = -9*eleventh;
356  m_V[0][2] = 2*eleventh;
357  m_V[0][3] = 18*eleventh;
358  m_V[0][4] = -18*eleventh;
359  m_V[0][5] = 6*eleventh;
360  m_V[1][0] = 1.0;
361  m_V[2][1] = 1.0;
362  m_V[4][3] = 1.0;
363  m_V[5][4] = 1.0;
364 
365  m_schemeType = eIMEX;
366  m_numMultiStepValues = 3;
367  m_numMultiStepDerivs = 3;
368  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
369  m_timeLevelOffset[0] = 0;
370  m_timeLevelOffset[1] = 1;
371  m_timeLevelOffset[2] = 2;
372  m_timeLevelOffset[3] = 0;
373  m_timeLevelOffset[4] = 1;
374  m_timeLevelOffset[5] = 2;
375  }
376  break;
377  case eAdamsMoultonOrder2:
378  {
379  m_numsteps = 2;
380  m_numstages = 1;
381 
382  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
383  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
384 
385  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.5);
386  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
387  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
388  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.0);
389 
390  m_B[0][0][0] = 0.5;
391  m_B[0][1][0] = 1.0;
392 
393  m_U[0][0] = 1.0;
394  m_U[0][1] = 0.5;
395 
396  m_V[0][0] = 1.0;
397  m_V[0][1] = 0.5;
398 
399 
400  m_schemeType = eDiagonallyImplicit;
401  m_numMultiStepValues = 1;
402  m_numMultiStepDerivs = 1;
403  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
404  m_timeLevelOffset[0] = 0;
405  m_timeLevelOffset[1] = 0;
406  }
407  break;
408  case eBDFImplicitOrder2:
409  {
410  NekDouble third = 1.0/3.0;
411  m_numsteps = 2;
412  m_numstages = 1;
413 
414  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
415  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
416 
417  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,2*third);
418  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
419  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
420  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.0);
421 
422  m_B[0][0][0] = 2*third;
423  m_B[0][1][0] = 0.0;
424 
425  m_U[0][0] = 4*third;
426  m_U[0][1] = -third;
427 
428  m_V[0][0] = 4*third;
429  m_V[0][1] = -third;
430  m_V[1][0] = 1.0;
431 
432  m_schemeType = eDiagonallyImplicit;
433  m_numMultiStepValues = 2;
434  m_numMultiStepDerivs = 0;
435  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
436  m_timeLevelOffset[0] = 0;
437  m_timeLevelOffset[1] = 1;
438  }
439  break;
440  case eMidpoint:
441  {
442  m_numsteps = 1;
443  m_numstages = 2;
444 
445  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
446  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
447 
448  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
449  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
450  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
451  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
452 
453  m_A[0][1][0] = 0.5;
454  m_B[0][0][1] = 1.0;
455 
456  m_schemeType = eExplicit;
457  m_numMultiStepValues = 1;
458  m_numMultiStepDerivs = 0;
459  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
460  m_timeLevelOffset[0] = 0;
461  }
462  break;
463  case eClassicalRungeKutta4:
464  {
465  m_numsteps = 1;
466  m_numstages = 4;
467 
468  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
469  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
470 
471  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
472  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
473  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
474  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
475 
476  m_A[0][1][0] = 0.5;
477  m_A[0][2][1] = 0.5;
478  m_A[0][3][2] = 1.0;
479 
480  m_B[0][0][0] = 1.0/6.0;
481  m_B[0][0][1] = 1.0/3.0;
482  m_B[0][0][2] = 1.0/3.0;
483  m_B[0][0][3] = 1.0/6.0;
484 
485  m_schemeType = eExplicit;
486  m_numMultiStepValues = 1;
487  m_numMultiStepDerivs = 0;
488  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
489  m_timeLevelOffset[0] = 0;
490  }
491  break;
492  case eRungeKutta2_ModifiedEuler:
493  {
494  m_numsteps = 1;
495  m_numstages = 2;
496 
497  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
498  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
499 
500  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
501  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
502  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.5);
503  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.5);
504 
505  m_A[0][1][0] = 0.5;
506 
507  m_B[0][0][1] = 1.0;
508 
509 
510  m_schemeType = eExplicit;
511  m_numMultiStepValues = 1;
512  m_numMultiStepDerivs = 0;
513  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
514  m_timeLevelOffset[0] = 0;
515  }
516  break;
517  case eRungeKutta2_ImprovedEuler:
518  {
519  m_numsteps = 1;
520  m_numstages = 2;
521 
522  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
523  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
524 
525  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
526  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.5);
527  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
528  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
529 
530  m_A[0][1][0] = 1.0;
531 
532  m_schemeType = eExplicit;
533  m_numMultiStepValues = 1;
534  m_numMultiStepDerivs = 0;
535  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
536  m_timeLevelOffset[0] = 0;
537  }
538  break;
539  case eDIRKOrder2:
540  {
541  m_numsteps = 1;
542  m_numstages = 2;
543 
544  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
545  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
546 
547  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
548  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
549  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
550  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
551 
552  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
553 
554  m_A[0][0][0] = lambda;
555  m_A[0][1][0] = 1.0 - lambda;
556  m_A[0][1][1] = lambda;
557 
558  m_B[0][0][0] = 1.0 - lambda;
559  m_B[0][0][1] = lambda;
560 
561  m_schemeType = eDiagonallyImplicit;
562  m_numMultiStepValues = 1;
563  m_numMultiStepDerivs = 0;
564  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
565  m_timeLevelOffset[0] = 0;
566  }
567  break;
568  case eDIRKOrder3:
569  {
570  m_numsteps = 1;
571  m_numstages = 3;
572 
573  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
574  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
575 
576  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
577  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
578  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
579  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
580 
581  NekDouble lambda = 0.4358665215;
582 
583  m_A[0][0][0] = lambda;
584  m_A[0][1][0] = 0.5 * (1.0 - lambda);
585  m_A[0][2][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
586  m_A[0][1][1] = lambda;
587  m_A[0][2][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
588  m_A[0][2][2] = lambda;
589 
590  m_B[0][0][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
591  m_B[0][0][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
592  m_B[0][0][2] = lambda;
593 
594  m_schemeType = eDiagonallyImplicit;
595  m_numMultiStepValues = 1;
596  m_numMultiStepDerivs = 0;
597  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
598  m_timeLevelOffset[0] = 0;
599  }
600  break;
601  case eIMEXdirk_2_3_2:
602  {
603  m_numsteps = 1;
604  m_numstages = 3;
605 
606  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
607  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
608 
609  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
610  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
611  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
612  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
613  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
614  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
615 
616  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
617  NekDouble delta = -2.0*sqrt(2.0)/3.0;
618 
619  m_A[0][1][1] = lambda;
620  m_A[0][2][1] = 1.0 - lambda;
621  m_A[0][2][2] = lambda;
622 
623  m_B[0][0][1] = 1.0 - lambda;
624  m_B[0][0][2] = lambda;
625 
626  m_A[1][1][0] = lambda;
627  m_A[1][2][0] = delta;
628  m_A[1][2][1] = 1.0 - delta;
629 
630  m_B[1][0][1] = 1.0 - lambda;
631  m_B[1][0][2] = lambda;
632 
633  m_schemeType = eIMEX;
634  m_numMultiStepValues = 1;
635  m_numMultiStepDerivs = 0;
636  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
637  m_timeLevelOffset[0] = 0;
638  }
639  break;
640  case eIMEXdirk_3_4_3:
641  {
642  m_numsteps = 1;
643  m_numstages = 4;
644 
645  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
646  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
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_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
651  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
652  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
653  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
654 
655  NekDouble lambda = 0.4358665215;
656 
657  m_A[0][1][1] = lambda;
658  m_A[0][2][1] = 0.5 * (1.0 - lambda);
659  m_A[0][3][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
660  m_A[0][2][2] = lambda;
661  m_A[0][3][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
662  m_A[0][3][3] = lambda;
663 
664  m_B[0][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
665  m_B[0][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
666  m_B[0][0][3] = lambda;
667 
668  m_A[1][1][0] = 0.4358665215;
669  m_A[1][2][0] = 0.3212788860;
670  m_A[1][2][1] = 0.3966543747;
671  m_A[1][3][0] =-0.105858296;
672  m_A[1][3][1] = 0.5529291479;
673  m_A[1][3][2] = 0.5529291479;
674 
675  m_B[1][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
676  m_B[1][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
677  m_B[1][0][3] = lambda;
678 
679  m_schemeType = eIMEX;
680  m_numMultiStepValues = 1;
681  m_numMultiStepDerivs = 0;
682  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
683  m_timeLevelOffset[0] = 0;
684  }
685  break;
686  case eCNAB:
687  {
688  NekDouble secondth = 1.0/2.0;
689  m_numsteps = 4;
690  m_numstages = 1;
691 
692  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
693  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
694 
695  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,secondth);
696  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
697  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
698  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
699  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,secondth);
700  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
701 
702  m_B[0][0][0] = secondth;
703  m_B[0][1][0] = 1.0;
704  m_B[1][2][0] = 1.0;
705  m_U[0][0] = 2*secondth;
706  m_U[0][2] = 3*secondth;
707  m_U[0][3] = -1*secondth;
708 
709  m_V[0][0] = 2*secondth;
710  m_V[0][1] = secondth;
711  m_V[0][2] = 3*secondth;
712  m_V[0][3] = -1*secondth;
713  m_V[3][2] = 1.0;
714 
715  m_schemeType = eIMEX;
716  m_numMultiStepValues = 1;
717  m_numMultiStepDerivs = 3;
718  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
719  m_timeLevelOffset[0] = 0;
720  m_timeLevelOffset[1] = 0;
721  m_timeLevelOffset[2] = 0;
722  m_timeLevelOffset[3] = 1;
723  }
724  break;
725  case eIMEXGear:
726  {
727  NekDouble twothirdth = 2.0/3.0;
728  m_numsteps = 2;
729  m_numstages = 1;
730 
731  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
732  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
733 
734  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
735  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
736  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
737  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
738  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,twothirdth);
739  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
740 
741  m_B[0][0][0] = twothirdth;
742  m_B[1][0][0] = twothirdth;
743  m_U[0][0] = 2*twothirdth;
744  m_U[0][1] = -0.5*twothirdth;
745 
746  m_V[0][0] = 2*twothirdth;
747  m_V[0][1] = -0.5*twothirdth;
748  m_V[1][0] = 1.0;
749 
750  m_schemeType = eIMEX;
751  m_numMultiStepValues = 2;
752  m_numMultiStepDerivs = 0;
753  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
754  m_timeLevelOffset[0] = 0;
755  m_timeLevelOffset[1] = 1;
756  }
757  break;
758  case eMCNAB:
759  {
760  NekDouble sixthx = 9.0/16.0;
761  m_numsteps = 5;
762  m_numstages = 1;
763 
764  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
765  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
766 
767  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,sixthx);
768  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
769  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
770  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
771  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
772  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
773 
774  m_B[0][0][0] = sixthx;
775  m_B[0][1][0] = 1.0;
776  m_B[1][3][0] = 1.0;
777  m_U[0][0] = 1.0;
778  m_U[0][1] = 6.0/16.0;
779  m_U[0][2] = 1.0/16.0;
780  m_U[0][3] = 1.5;
781  m_U[0][4] = -0.5;
782 
783  m_V[0][0] = 1.0;
784  m_V[0][1] = 6.0/16.0;
785  m_V[0][2] = 1.0/16.0;
786  m_V[0][3] = 1.5;
787  m_V[0][4] = -0.5;
788  m_V[2][1] = 1.0;
789  m_V[4][3] = 1.0;
790 
791  m_schemeType = eIMEX;
792  m_numMultiStepValues = 1;
793  m_numMultiStepDerivs = 4;
794  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
795  m_timeLevelOffset[0] = 0;
796  m_timeLevelOffset[1] = 0;
797  m_timeLevelOffset[2] = 1;
798  m_timeLevelOffset[3] = 0;
799  m_timeLevelOffset[4] = 1;
800  }
801  break;
802  case eIMEXdirk_2_2_2:
803  {
804  m_numsteps = 1;
805  m_numstages = 3;
806 
807  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
808  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
809 
810  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
811  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
812  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
813  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
814  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
815  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
816 
817  NekDouble glambda = 0.788675134594813;
818  NekDouble gdelta = 0.366025403784439;
819 
820  m_A[0][1][1] = glambda;
821  m_A[0][2][1] = 1.0 - glambda;
822  m_A[0][2][2] = glambda;
823 
824  m_B[0][0][1] = 1.0 - glambda;
825  m_B[0][0][2] = glambda;
826 
827  m_A[1][1][0] = glambda;
828  m_A[1][2][0] = gdelta;
829  m_A[1][2][1] = 1.0 - gdelta;
830 
831  m_B[1][0][0] = gdelta;
832  m_B[1][0][1] = 1.0 - gdelta;
833 
834  m_schemeType = eIMEX;
835  m_numMultiStepValues = 1;
836  m_numMultiStepDerivs = 0;
837  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
838  m_timeLevelOffset[0] = 0;
839  }
840  break;
841  case eIMEXdirk_2_3_3:
842  {
843  m_numsteps = 1;
844  m_numstages = 3;
845 
846  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
847  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
848 
849  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
850  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
851  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
852  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
853  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
854  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
855 
856  NekDouble glambda = 0.788675134594813;
857 
858  m_A[0][1][1] = glambda;
859  m_A[0][2][1] = 1.0 - 2.0*glambda;
860  m_A[0][2][2] = glambda;
861 
862  m_B[0][0][1] = 0.5;
863  m_B[0][0][2] = 0.5;
864 
865  m_A[1][1][0] = glambda;
866  m_A[1][2][0] = glambda - 1.0;
867  m_A[1][2][1] = 2.0*(1-glambda);
868 
869  m_B[1][0][1] = 0.5;
870  m_B[1][0][2] = 0.5;
871 
872  m_schemeType = eIMEX;
873  m_numMultiStepValues = 1;
874  m_numMultiStepDerivs = 0;
875  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
876  m_timeLevelOffset[0] = 0;
877  }
878  break;
879  case eIMEXdirk_1_1_1:
880  {
881  m_numsteps = 1;
882  m_numstages = 2;
883 
884  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
885  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
886 
887  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
888  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
889  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
890  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
891  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
892  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
893 
894  m_A[0][1][1] = 1.0;
895 
896  m_B[0][0][1] = 1.0;
897 
898  m_A[1][1][0] = 1.0;
899 
900  m_B[1][0][0] = 1.0;
901 
902  m_schemeType = eIMEX;
903  m_numMultiStepValues = 1;
904  m_numMultiStepDerivs = 0;
905  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
906  m_timeLevelOffset[0] = 0;
907  }
908  break;
909  case eIMEXdirk_1_2_1:
910  {
911  m_numsteps = 1;
912  m_numstages = 2;
913 
914  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
915  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
916 
917  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
918  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
919  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
920  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
921  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
922  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
923 
924 
925  m_A[0][1][1] = 1.0;
926 
927  m_B[0][0][1] = 1.0;
928 
929  m_A[1][1][0] = 1.0;
930 
931  m_B[1][0][1] = 1.0;
932 
933  m_schemeType = eIMEX;
934  m_numMultiStepValues = 1;
935  m_numMultiStepDerivs = 0;
936  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
937  m_timeLevelOffset[0] = 0;
938  }
939  break;
940  case eIMEXdirk_1_2_2:
941  {
942  m_numsteps = 1;
943  m_numstages = 2;
944 
945  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
946  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
947 
948  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
949  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
950  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
951  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
952  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
953  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
954 
955  m_A[0][1][1] = 1.0/2.0;
956 
957  m_B[0][0][1] = 1.0;
958 
959  m_A[1][1][0] = 1.0/2.0;
960 
961  m_B[1][0][1] = 1.0;
962 
963  m_schemeType = eIMEX;
964  m_numMultiStepValues = 1;
965  m_numMultiStepDerivs = 0;
966  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
967  m_timeLevelOffset[0] = 0;
968  }
969  break;
970  case eIMEXdirk_4_4_3:
971  {
972  m_numsteps = 1;
973  m_numstages = 5;
974 
975  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
976  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
977 
978  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
979  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
980  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
981  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
982  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
983  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
984 
985  m_A[0][1][1] = 1.0/2.0;
986  m_A[0][2][1] = 1.0/6.0;
987  m_A[0][2][2] = 1.0/2.0;
988  m_A[0][3][1] = -1.0/2.0;
989  m_A[0][3][2] = 1.0/2.0;
990  m_A[0][3][3] = 1.0/2.0;
991  m_A[0][4][1] = 3.0/2.0;
992  m_A[0][4][2] = -3.0/2.0;
993  m_A[0][4][3] = 1.0/2.0;
994  m_A[0][4][4] = 1.0/2.0;
995 
996 
997  m_B[0][0][1] = 3.0/2.0;
998  m_B[0][0][2] = -3.0/2.0;
999  m_B[0][0][3] = 1.0/2.0;
1000  m_B[0][0][4] = 1.0/2.0;
1001 
1002  m_A[1][1][0] = 1.0/2.0;
1003  m_A[1][2][0] = 11.0/18.0;
1004  m_A[1][2][1] = 1.0/18.0;
1005  m_A[1][3][0] = 5.0/6.0;
1006  m_A[1][3][1] = -5.0/6.0;
1007  m_A[1][3][2] = 1.0/2.0;
1008  m_A[1][4][0] = 1.0/4.0;
1009  m_A[1][4][1] = 7.0/4.0;
1010  m_A[1][4][2] = 3.0/4.0;
1011  m_A[1][4][3] = -7.0/4.0;
1012 
1013  m_B[1][0][0] = 1.0/4.0;
1014  m_B[1][0][1] = 7.0/4.0;
1015  m_B[1][0][2] = 3.0/4.0;
1016  m_B[1][0][3] = -7.0/4.0;
1017 
1018  m_schemeType = eIMEX;
1019  m_numMultiStepValues = 1;
1020  m_numMultiStepDerivs = 0;
1021  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1022  m_timeLevelOffset[0] = 0;
1023  }
1024  break;
1025  default:
1026  {
1027  NEKERROR(ErrorUtil::efatal,"Invalid Time Integration Scheme");
1028  }
1029  }
1030 
1031  m_firstStageEqualsOldSolution = CheckIfFirstStageEqualsOldSolution(m_A,m_B,m_U,m_V);
1032  m_lastStageEqualsNewSolution = CheckIfLastStageEqualsNewSolution(m_A,m_B,m_U,m_V);
1033 
1034  ASSERTL1(VerifyIntegrationSchemeType(m_schemeType,m_A,m_B,m_U,m_V),
1035  "Time integration scheme coefficients do not match its type");
1036  }
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:132
Nektar::LibUtilities::eBDFImplicitOrder1
BDF multi-step scheme of order 1 (implicit)
Definition: TimeIntegrationScheme.h:72
Nektar::LibUtilities::eAdamsBashforthOrder2
Adams-Bashforth Forward multi-step scheme of order 2.
Definition: TimeIntegrationScheme.h:68
Nektar::LibUtilities::TimeIntegrationScheme::m_B
Array< OneD, Array< TwoD, NekDouble > > m_B
Definition: TimeIntegrationScheme.h:530
Nektar::LibUtilities::eClassicalRungeKutta4
Runge-Kutta multi-stage scheme 4th order explicit.
Definition: TimeIntegrationScheme.h:74
Nektar::LibUtilities::eIMEXdirk_1_2_2
Implicit-Explicit Midpoint IMEX DIRK(1,2,2)
Definition: TimeIntegrationScheme.h:90
Nektar::LibUtilities::eIMEXdirk_4_4_3
L-stable, four stage, third order IMEX DIRK(4,4,3)
Definition: TimeIntegrationScheme.h:95
Nektar::LibUtilities::eIMEXdirk_1_2_1
Forward-Backward Euler IMEX DIRK(1,2,1)
Definition: TimeIntegrationScheme.h:89
Nektar::LibUtilities::eIMEX
Implicit Explicit General Linear Method.
Definition: TimeIntegrationScheme.h:138
Nektar::LibUtilities::eAdamsMoultonOrder2
Adams-Moulton Forward multi-step scheme of order 2.
Definition: TimeIntegrationScheme.h:71
Nektar::LibUtilities::eAdamsBashforthOrder3
Adams-Bashforth Forward multi-step scheme of order 3.
Definition: TimeIntegrationScheme.h:69
Nektar::LibUtilities::eCNAB
Crank-Nicolson/Adams-Bashforth Order 2 (CNAB)
Definition: TimeIntegrationScheme.h:85
Nektar::LibUtilities::eIMEXOrder2
IMEX 2nd order scheme using Backward Different Formula & Extrapolation.
Definition: TimeIntegrationScheme.h:80
Nektar::LibUtilities::eDiagonallyImplicit
Diagonally implicit scheme (e.g. the DIRK schemes)
Definition: TimeIntegrationScheme.h:137
Nektar::LibUtilities::TimeIntegrationScheme::CheckIfFirstStageEqualsOldSolution
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
Definition: TimeIntegrationScheme.cpp:1679
Nektar::LibUtilities::eAdamsMoultonOrder1
Adams-Moulton Forward multi-step scheme of order 1.
Definition: TimeIntegrationScheme.h:70
Nektar::LibUtilities::eAdamsBashforthOrder1
Adams-Bashforth Forward multi-step scheme of order 1.
Definition: TimeIntegrationScheme.h:67
Nektar::LibUtilities::eIMEXdirk_3_4_3
L-stable, three stage, third order IMEX DIRK(3,4,3)
Definition: TimeIntegrationScheme.h:94
Nektar::LibUtilities::eRungeKutta2_ImprovedEuler
Runge-Kutta multi-stage scheme 2nd order explicit.
Definition: TimeIntegrationScheme.h:76
Nektar::LibUtilities::eIMEXdirk_1_1_1
Forward-Backward Euler IMEX DIRK(1,1,1)
Definition: TimeIntegrationScheme.h:88
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::LibUtilities::eBDFImplicitOrder2
BDF multi-step scheme of order 2 (implicit)
Definition: TimeIntegrationScheme.h:73
Nektar::LibUtilities::eIMEXdirk_2_2_2
L-stable, two stage, second order IMEX DIRK(2,2,2)
Definition: TimeIntegrationScheme.h:91
Nektar::LibUtilities::eRungeKutta2_ModifiedEuler
Runge-Kutta multi-stage scheme 2nd order explicit.
Definition: TimeIntegrationScheme.h:75
Nektar::LibUtilities::TimeIntegrationScheme::VerifyIntegrationSchemeType
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
Definition: TimeIntegrationScheme.cpp:1040
Nektar::LibUtilities::eDIRKOrder2
Diagonally Implicit Runge Kutta scheme of order 3.
Definition: TimeIntegrationScheme.h:83
Nektar::LibUtilities::eIMEXOrder3
IMEX 3rd order scheme using Backward Different Formula & Extrapolation.
Definition: TimeIntegrationScheme.h:81
Nektar::LibUtilities::eIMEXdirk_2_3_2
L-stable, three stage, third order IMEX DIRK(3,4,3)
Definition: TimeIntegrationScheme.h:92
Nektar::LibUtilities::eDIRKOrder3
Diagonally Implicit Runge Kutta scheme of order 3.
Definition: TimeIntegrationScheme.h:84
Nektar::LibUtilities::eIMEXdirk_2_3_3
L-stable, two stage, third order IMEX DIRK(2,3,3)
Definition: TimeIntegrationScheme.h:93
Nektar::LibUtilities::TimeIntegrationScheme::CheckIfLastStageEqualsNewSolution
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
Definition: TimeIntegrationScheme.cpp:1717
Nektar::LibUtilities::TimeIntegrationScheme::m_A
Array< OneD, Array< TwoD, NekDouble > > m_A
Definition: TimeIntegrationScheme.h:529
Nektar::LibUtilities::eMCNAB
Modified Crank-Nicolson/Adams-Bashforth Order 2 (MCNAB)
Definition: TimeIntegrationScheme.h:87
Nektar::LibUtilities::eIMEXOrder1
IMEX 1st order scheme using Euler Backwards/Euler Forwards.
Definition: TimeIntegrationScheme.h:79
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:165
Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrationScheme ( )
inlineprivate

Definition at line 554 of file TimeIntegrationScheme.h.

References ErrorUtil::efatal, and NEKERROR.

554  :m_schemeKey(NullTimeIntegrationSchemeKey)
555  {
556  NEKERROR(ErrorUtil::efatal,"Default Constructor for the TimeIntegrationScheme class should not be called");
557  }
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:132
Nektar::LibUtilities::NullTimeIntegrationSchemeKey
static const TimeIntegrationSchemeKey NullTimeIntegrationSchemeKey(eNoTimeIntegrationMethod)
Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrationScheme ( const TimeIntegrationScheme &  in)
inlineprivate

Definition at line 559 of file TimeIntegrationScheme.h.

References ErrorUtil::efatal, and NEKERROR.

559  :m_schemeKey(NullTimeIntegrationSchemeKey)
560  {
561  NEKERROR(ErrorUtil::efatal,"Copy Constructor for the TimeIntegrationScheme class should not be called");
562  }
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:132
Nektar::LibUtilities::NullTimeIntegrationSchemeKey
static const TimeIntegrationSchemeKey NullTimeIntegrationSchemeKey(eNoTimeIntegrationMethod)

Member Function Documentation

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

Definition at line 389 of file TimeIntegrationScheme.h.

References m_A.

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

390  {
391  return m_A[0][i][j];
392  }
Nektar::LibUtilities::TimeIntegrationScheme::m_A
Array< OneD, Array< TwoD, NekDouble > > m_A
Definition: TimeIntegrationScheme.h:529
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::A_IMEX ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 409 of file TimeIntegrationScheme.h.

References m_A.

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

410  {
411  return m_A[1][i][j];
412  }
Nektar::LibUtilities::TimeIntegrationScheme::m_A
Array< OneD, Array< TwoD, NekDouble > > m_A
Definition: TimeIntegrationScheme.h:529
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::B ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 394 of file TimeIntegrationScheme.h.

References m_B.

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

395  {
396  return m_B[0][i][j];
397  }
Nektar::LibUtilities::TimeIntegrationScheme::m_B
Array< OneD, Array< TwoD, NekDouble > > m_B
Definition: TimeIntegrationScheme.h:530
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::B_IMEX ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 414 of file TimeIntegrationScheme.h.

References m_B.

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

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

Definition at line 1679 of file TimeIntegrationScheme.cpp.

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

Referenced by TimeIntegrationScheme().

1683  {
1684  int i,m;
1685  // First stage equals old solution if:
1686  // 1. the first row of the coefficient matrix A consists of zeros
1687  // 2. U[0][0] is equal to one and all other first row entries of U are zero
1688 
1689  // 1. Check first condition
1690  for(m = 0; m < A.num_elements(); m++)
1691  {
1692  for(i = 0; i < m_numstages; i++)
1693  {
1694  if( fabs(A[m][0][i]) > NekConstants::kNekZeroTol )
1695  {
1696  return false;
1697  }
1698  }
1699  }
1700 
1701  // 2. Check second condition
1702  if( fabs(U[0][0] - 1.0) > NekConstants::kNekZeroTol )
1703  {
1704  return false;
1705  }
1706  for(i = 1; i < m_numsteps; i++)
1707  {
1708  if( fabs(U[0][i]) > NekConstants::kNekZeroTol )
1709  {
1710  return false;
1711  }
1712  }
1713 
1714  return true;
1715  }
Nektar::LibUtilities::TimeIntegrationScheme::U
NekDouble U(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:399
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Nektar::LibUtilities::TimeIntegrationScheme::A
NekDouble A(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:389
bool Nektar::LibUtilities::TimeIntegrationScheme::CheckIfLastStageEqualsNewSolution ( const Array< OneD, const Array< TwoD, NekDouble > > &  A,
const Array< OneD, const Array< TwoD, NekDouble > > &  B,
const Array< TwoD, const NekDouble > &  U,
const Array< TwoD, const NekDouble > &  V 
) const
private

Definition at line 1717 of file TimeIntegrationScheme.cpp.

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

Referenced by TimeIntegrationScheme().

1721  {
1722  int i,m;
1723  // Last stage equals new solution if:
1724  // 1. the last row of the coefficient matrix A is equal to the first row of matrix B
1725  // 2. the last row of the coefficient matrix U is equal to the first row of matrix V
1726 
1727  // 1. Check first condition
1728  for(m = 0; m < A.num_elements(); m++)
1729  {
1730  for(i = 0; i < m_numstages; i++)
1731  {
1732  if( fabs(A[m][m_numstages-1][i]-B[m][0][i]) > NekConstants::kNekZeroTol )
1733  {
1734  return false;
1735  }
1736  }
1737  }
1738 
1739  // 2. Check second condition
1740  for(i = 0; i < m_numsteps; i++)
1741  {
1742  if( fabs(U[m_numstages-1][i]-V[0][i]) > NekConstants::kNekZeroTol )
1743  {
1744  return false;
1745  }
1746  }
1747 
1748  return true;
1749  }
Nektar::LibUtilities::TimeIntegrationScheme::U
NekDouble U(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:399
Nektar::LibUtilities::TimeIntegrationScheme::B
NekDouble B(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:394
Nektar::LibUtilities::TimeIntegrationScheme::V
NekDouble V(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:404
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Nektar::LibUtilities::TimeIntegrationScheme::A
NekDouble A(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:389
bool Nektar::LibUtilities::TimeIntegrationScheme::CheckTimeIntegrateArguments ( const NekDouble  timestep,
ConstTripleArray &  y_old,
ConstSingleArray &  t_old,
TripleArray &  y_new,
SingleArray &  t_new,
const TimeIntegrationSchemeOperators &  op 
) const
private

Definition at line 1751 of file TimeIntegrationScheme.cpp.

References ASSERTL1, and m_numsteps.

Referenced by TimeIntegrate().

1757  {
1758  // Check if arrays are all of consistent size
1759  ASSERTL1(y_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1760  ASSERTL1(y_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1761 
1762  ASSERTL1(y_old[0]. num_elements()==y_new[0]. num_elements(),"Non-matching number of variables.");
1763  ASSERTL1(y_old[0][0].num_elements()==y_new[0][0].num_elements(),"Non-matching number of coefficients.");
1764 
1765  ASSERTL1(t_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1766  ASSERTL1(t_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1767 
1768  return true;
1769  }
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:165
TimeIntegrationSchemeSharedPtr Nektar::LibUtilities::TimeIntegrationScheme::Create ( const TimeIntegrationSchemeKey &  key)
staticprivate

Definition at line 143 of file TimeIntegrationScheme.cpp.

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

144  {
145  TimeIntegrationSchemeSharedPtr returnval(
146  MemoryManager<TimeIntegrationScheme>::AllocateSharedPtr(key));
147  return returnval;
148  }
Nektar::LibUtilities::TimeIntegrationSchemeSharedPtr
boost::shared_ptr< TimeIntegrationScheme > TimeIntegrationSchemeSharedPtr
Definition: TimeIntegrationScheme.h:51
int Nektar::LibUtilities::TimeIntegrationScheme::GetFirstDim ( ConstTripleArray &  y) const
inlineprivate

Definition at line 578 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 379 of file TimeIntegrationScheme.h.

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

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

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

Definition at line 374 of file TimeIntegrationScheme.h.

References m_schemeKey.

Referenced by TimeIntegrate().

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

Definition at line 384 of file TimeIntegrationScheme.h.

References m_schemeType.

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

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

Definition at line 434 of file TimeIntegrationScheme.h.

References m_numMultiStepDerivs.

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

Definition at line 429 of file TimeIntegrationScheme.h.

References m_numMultiStepValues.

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

Definition at line 424 of file TimeIntegrationScheme.h.

References m_numstages.

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

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

Definition at line 419 of file TimeIntegrationScheme.h.

References m_numsteps.

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

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

Definition at line 583 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 499 of file TimeIntegrationScheme.h.

References m_timeLevelOffset.

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

This function initialises the time integration scheme.

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

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

Parameters
timestepThe size of the timestep, i.e. $\Delta t$.
timeon input: the initial time, i.e. $t_0$.
timeon output: the new time-level after initialisation, in general this yields $t = t_0 + (r-1)\Delta t$ where $r$ is the number of steps of the multi-step method.
nstepson output: he number of initialisation steps required. In general this corresponds to $r-1$ where $r$ is the number of steps of the multi-step method.
fan object of the class FuncType, where FuncType should have a method FuncType::ODEforcing to evaluate the right hand side $f(t,\boldsymbol{y})$ of the ODE.
y_0the initial value $\boldsymbol{y}(t_0)$
Returns
An object of the class TimeIntegrationSolution which represents the vectors $\boldsymbol{y}^{[n]}$ and $t^{[n]}$ that can be used to start the actual integration.

Definition at line 1095 of file TimeIntegrationScheme.cpp.

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

1099  {
1100  // create a TimeIntegrationSolution object based upon the
1101  // initial value. Initialise all other multi-step values
1102  // and derivatives to zero
1103  TimeIntegrationSolutionSharedPtr y_out =
1104  MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(m_schemeKey,y_0,time,timestep);
1105 
1106  // calculate the initial derivative, if is part of the
1107  // solution vector of the current scheme
1108  if(m_numMultiStepDerivs)
1109  {
1110  if(m_timeLevelOffset[m_numMultiStepValues] == 0)
1111  {
1112  int i;
1113  int nvar = y_0.num_elements();
1114  int npoints = y_0[0].num_elements();
1115  DoubleArray f_y_0(nvar);
1116  for(i = 0; i < nvar; i++)
1117  {
1118  f_y_0[i] = Array<OneD,NekDouble>(npoints);
1119  }
1120  // calculate the derivative of the initial value
1121  op.DoOdeRhs(y_0,f_y_0,time);
1122 
1123  // multiply by the step size
1124  for(i = 0; i < nvar; i++)
1125  {
1126  Blas::Dscal(npoints,timestep,f_y_0[i].get(),1);
1127  }
1128  y_out->SetDerivative(0,f_y_0,timestep);
1129  }
1130  }
1131 
1132  return y_out;
1133  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::LibUtilities::TimeIntegrationScheme::DoubleArray
Array< OneD, Array< OneD, NekDouble > > DoubleArray
Definition: TimeIntegrationScheme.h:362
Nektar::LibUtilities::TimeIntegrationSolutionSharedPtr
boost::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
Definition: TimeIntegrationScheme.h:57
TimeIntegrationScheme::ConstDoubleArray & Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrate ( const NekDouble  timestep,
TimeIntegrationSolutionSharedPtr &  y,
const TimeIntegrationSchemeOperators &  op 
)

Explicit integration of an ODE.

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

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

Parameters
timestepThe size of the timestep, i.e. $\Delta t$.
fan object of the class FuncType, where FuncType should have a method FuncType::ODEforcing to evaluate the right hand side $f(t,\boldsymbol{y})$ of the ODE.
yon input: the vectors $\boldsymbol{y}^{[n-1]}$ and $t^{[n-1]}$ (which corresponds to the solution at the old time level)
yon output: the vectors $\boldsymbol{y}^{[n]}$ and $t^{[n]}$ (which corresponds to the solution at the old new level)
Returns
The actual solution $\boldsymbol{y}^{n}$ at the new time level (which in fact is also embedded in the argument y).

Definition at line 1136 of file TimeIntegrationScheme.cpp.

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

1139  {
1140  ASSERTL1(!(GetIntegrationSchemeType() == eImplicit),
1141  "Fully Implicit integration scheme cannot be handled by this routine.");
1142 
1143  int nvar = solvector->GetFirstDim ();
1144  int npoints = solvector->GetSecondDim();
1145 
1146  if( (solvector->GetIntegrationScheme()).get() != this )
1147  {
1148  // This branch will be taken when the solution vector
1149  // (solvector) is set up for a different scheme than
1150  // the object this method is called from. (typically
1151  // needed to calculate the first time-levels of a
1152  // multi-step scheme)
1153 
1154  // To do this kind of 'non-matching' integration, we
1155  // perform the following three steps:
1156  //
1157  // 1: copy the required input information from the
1158  // solution vector of the master scheme to the
1159  // input solution vector of the current scheme
1160  //
1161  // 2: time-integrate for one step using the current
1162  // scheme
1163  //
1164  // 3: copy the information contained in the output
1165  // vector of the current scheme to the solution
1166  // vector of the master scheme
1167 
1168  // STEP 1: copy the required input information from
1169  // the solution vector of the master scheme
1170  // to the input solution vector of the
1171  // current scheme
1172 
1173  // 1.1 Determine which information is required for the
1174  // current scheme
1175  int n;
1176  DoubleArray y_n;
1177  NekDouble t_n = 0;
1178  DoubleArray dtFy_n;
1179  unsigned int nCurSchemeVals = m_numMultiStepValues; // number of required values of the current scheme
1180  unsigned int nCurSchemeDers = m_numMultiStepDerivs; // number of required derivs of the current scheme
1181  unsigned int nCurSchemeSteps = m_numsteps; // number of steps in the current scheme
1182  unsigned int nMasterSchemeVals = solvector->GetNvalues(); // number of values of the master scheme
1183  unsigned int nMasterSchemeDers = solvector->GetNderivs(); // number of derivs of the master scheme
1184  // The arrays below contains information to which
1185  // time-level the values and derivatives of the
1186  // schemes belong
1187  const Array<OneD, const unsigned int>& curTimeLevels = m_timeLevelOffset;
1188  const Array<OneD, const unsigned int>& masterTimeLevels = solvector->GetTimeLevelOffset();
1189 
1190  // 1.2 Copy the required information from the master
1191  // solution vector to the input solution vector of
1192  // the current scheme
1193  TimeIntegrationSolutionSharedPtr solvector_in = MemoryManager<TimeIntegrationSolution>::
1194  AllocateSharedPtr(GetIntegrationSchemeKey()); // input solution vector of the current scheme
1195 
1196  for(n = 0; n < nCurSchemeVals; n++)
1197  {
1198  // Get the required value out of the master solution vector
1199  //DoubleArray& y_n = solvector->GetValue ( curTimeLevels[n] );
1200  //NekDouble t_n = solvector->GetValueTime( curTimeLevels[n] );
1201 
1202  y_n = solvector->GetValue ( curTimeLevels[n] );
1203  t_n = solvector->GetValueTime( curTimeLevels[n] );
1204 
1205  // Set the required value in the input solution
1206  // vector of the current scheme
1207  solvector_in->SetValue(curTimeLevels[n],y_n,t_n);
1208  }
1209  for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
1210  {
1211  // Get the required derivative out of the master
1212  // solution vector
1213  //DoubleArray& dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
1214  dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
1215 
1216  // Set the required derivative in the input
1217  // solution vector of the current scheme
1218  solvector_in->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
1219  }
1220 
1221  // STEP 2: time-integrate for one step using the
1222  // current scheme
1223  TimeIntegrationSolutionSharedPtr solvector_out = MemoryManager<TimeIntegrationSolution>:: AllocateSharedPtr(GetIntegrationSchemeKey(),nvar,npoints); // output solution vector of the current scheme
1224 
1225  // integrate
1226  TimeIntegrate(timestep, solvector_in->GetSolutionVector(),
1227  solvector_in->GetTimeVector(),
1228  solvector_out->UpdateSolutionVector(),
1229  solvector_out->UpdateTimeVector(),op);
1230 
1231 
1232  // STEP 3: copy the information contained in the
1233  // output vector of the current scheme to the
1234  // solution vector of the master scheme
1235 
1236  // 3.1 Check whether the current time scheme updates
1237  // the most recent derivative that should be
1238  // updated in the master scheme. If not,
1239  // calculate the derivative. This can be done
1240  // based upon the corresponding value and the
1241  // DoOdeRhs operator.
1242  int j;
1243  bool CalcNewDeriv = false; // flag inidicating whether the new derivative is availble in the output of
1244  // of the current scheme or whether it should be calculated
1245  if( nMasterSchemeDers > 0 )
1246  {
1247  if(nCurSchemeDers == 0)
1248  {
1249  CalcNewDeriv = true;
1250  }
1251  else
1252  {
1253  if( masterTimeLevels[nMasterSchemeVals] < curTimeLevels[nCurSchemeVals] )
1254  {
1255  CalcNewDeriv = true;
1256  }
1257  }
1258  }
1259 
1260  if(CalcNewDeriv)
1261  {
1262  int newDerivTimeLevel = masterTimeLevels[nMasterSchemeVals]; // contains the time level at which
1263  // we want to know the derivative of the
1264  // master scheme
1265  //DoubleArray y_n;
1266  //NekDouble t_n;
1267  // if the time level correspond to 0, calculate the derivative based upon the solution value
1268  // at the new time-level
1269  if (newDerivTimeLevel == 0)
1270  {
1271  y_n = solvector_out->GetValue(0);
1272  t_n = solvector_out->GetValueTime(0);
1273  }
1274  // if the time level correspond to 1, calculate the derivative based upon the solution value
1275  // at the new old-level
1276  else if( newDerivTimeLevel == 1 )
1277  {
1278  y_n = solvector->GetValue(0);
1279  t_n = solvector->GetValueTime(0);
1280  }
1281  else
1282  {
1283  ASSERTL1(false,"Problems with initialising scheme");
1284  }
1285 
1286  DoubleArray f_n(nvar);
1287  for(j = 0; j < nvar; j++)
1288  {
1289  f_n[j] = Array<OneD, NekDouble>(npoints);
1290  }
1291 
1292  // calculate the derivative
1293  op.DoOdeRhs(y_n, f_n, t_n);
1294 
1295  // multiply by dt (as required by the General Linear Method framework)
1296  for(j = 0; j < nvar; j++)
1297  {
1298  Vmath::Smul(npoints,timestep,f_n[j],1,
1299  f_n[j],1);
1300  }
1301 
1302  // Rotate the solution vector
1303  // (i.e. updating without calculating/inserting new values)
1304  solvector->RotateSolutionVector();
1305  // Set the calculated derivative in the master solution vector
1306  solvector->SetDerivative(newDerivTimeLevel,f_n,timestep);
1307  }
1308  else
1309  {
1310  // Rotate the solution vector (i.e. updating
1311  // without calculating/inserting new values)
1312  solvector->RotateSolutionVector();
1313  }
1314 
1315 
1316  // 1.2 Copy the information calculated using the
1317  // current scheme from the output solution vector
1318  // to the master solution vector
1319  for(n = 0; n < nCurSchemeVals; n++)
1320  {
1321  // Get the calculated value out of the output
1322  // solution vector of the current scheme
1323  //DoubleArray& y_n = solvector_out->GetValue ( curTimeLevels[n] );
1324  //NekDouble t_n = solvector_out->GetValueTime( curTimeLevels[n] );
1325  y_n = solvector_out->GetValue ( curTimeLevels[n] );
1326  t_n = solvector_out->GetValueTime( curTimeLevels[n] );
1327 
1328  // Set the calculated value in the master solution vector
1329  solvector->SetValue(curTimeLevels[n],y_n,t_n);
1330  }
1331 
1332  for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
1333  {
1334  // Get the calculated derivative out of the output
1335  // solution vector of the current scheme
1336  // DoubleArray& dtFy_n =
1337  // solvector_out->GetDerivative (curTimeLevels[n]);
1338  dtFy_n = solvector_out->GetDerivative ( curTimeLevels[n] );
1339 
1340  // Set the calculated derivative in the master
1341  // solution vector
1342  solvector->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
1343  }
1344  }
1345  else
1346  {
1347  const TimeIntegrationSchemeKey& key = solvector->GetIntegrationSchemeKey();
1348 
1349  TimeIntegrationSolutionSharedPtr solvector_new = MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(key,nvar,npoints);
1350 
1351  TimeIntegrate(timestep,solvector->GetSolutionVector(),
1352  solvector->GetTimeVector(),
1353  solvector_new->UpdateSolutionVector(),
1354  solvector_new->UpdateTimeVector(),op);
1355 
1356  solvector = solvector_new;
1357  }
1358  return solvector->GetSolution();
1359  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Vmath::Smul
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:199
Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeKey
const TimeIntegrationSchemeKey & GetIntegrationSchemeKey() const
Definition: TimeIntegrationScheme.h:374
Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrate
ConstDoubleArray & TimeIntegrate(const NekDouble timestep, TimeIntegrationSolutionSharedPtr &y, const TimeIntegrationSchemeOperators &op)
Explicit integration of an ODE.
Definition: TimeIntegrationScheme.cpp:1136
Nektar::LibUtilities::TimeIntegrationScheme::DoubleArray
Array< OneD, Array< OneD, NekDouble > > DoubleArray
Definition: TimeIntegrationScheme.h:362
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeType
TimeIntegrationSchemeType GetIntegrationSchemeType() const
Definition: TimeIntegrationScheme.h:384
Nektar::LibUtilities::TimeIntegrationSolutionSharedPtr
boost::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
Definition: TimeIntegrationScheme.h:57
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:165
void Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrate ( const NekDouble  timestep,
ConstTripleArray &  y_old,
ConstSingleArray &  t_old,
TripleArray &  y_new,
SingleArray &  t_new,
const TimeIntegrationSchemeOperators &  op 
)
private

Definition at line 1361 of file TimeIntegrationScheme.cpp.

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

1367  {
1368  ASSERTL1(CheckTimeIntegrateArguments(timestep,y_old,t_old,y_new,t_new,op), "Arguments not well defined");
1369 
1370  unsigned int i,j,k;
1371  TimeIntegrationSchemeType type = GetIntegrationSchemeType();
1372 
1373  // Check if storage has already been initialised.
1374  // If so, we just zero the temporary storage.
1375  if (m_initialised && m_nvar == GetFirstDim(y_old)
1376  && m_npoints == GetSecondDim(y_old))
1377  {
1378  for(j = 0; j < m_nvar; j++)
1379  {
1380  Vmath::Zero(m_npoints, m_tmp[j], 1);
1381  }
1382  }
1383  else
1384  {
1385  m_nvar = GetFirstDim(y_old);
1386  m_npoints = GetSecondDim(y_old);
1387 
1388  // First, we are going to calculate the various stage
1389  // values and stage derivatives (this is the multi-stage
1390  // part of the method)
1391  // - m_Y corresponds to the stage values
1392  // - m_F corresponds to the stage derivatives
1393  // - m_T corresponds to the time at the different stages
1394  // - m_tmp corresponds to the explicit right hand side of
1395  // each stage equation
1396  // (for explicit schemes, this correspond to m_Y)
1397 
1398  // Allocate memory for the arrays m_Y and m_F and m_tmp The same
1399  // storage will be used for every stage -> m_Y is a
1400  // DoubleArray
1401  m_tmp = DoubleArray(m_nvar);
1402  for(j = 0; j < m_nvar; j++)
1403  {
1404  m_tmp[j] = Array<OneD, NekDouble>(m_npoints,0.0);
1405  }
1406 
1407  // The same storage will be used for every stage -> m_tmp is
1408  // a DoubleArray
1409  if(type == eExplicit)
1410  {
1411  m_Y = m_tmp;
1412  }
1413  else
1414  {
1415  m_Y = DoubleArray(m_nvar);
1416  for(j = 0; j < m_nvar; j++)
1417  {
1418  m_Y[j] = Array<OneD, NekDouble>(m_npoints,0.0);
1419  }
1420  }
1421 
1422  // Different storage for every stage derivative as the data
1423  // will be re-used to update the solution -> m_F is a TripleArray
1424  m_F = TripleArray(m_numstages);
1425  for(i = 0; i < m_numstages; ++i)
1426  {
1427  m_F[i] = DoubleArray(m_nvar);
1428  for(j = 0; j < m_nvar; j++)
1429  {
1430  m_F[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
1431  }
1432  }
1433 
1434  if(type == eIMEX)
1435  {
1436  m_F_IMEX = TripleArray(m_numstages);
1437  for(i = 0; i < m_numstages; ++i)
1438  {
1439  m_F_IMEX[i] = DoubleArray(m_nvar);
1440  for(j = 0; j < m_nvar; j++)
1441  {
1442  m_F_IMEX[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
1443  }
1444  }
1445  }
1446 
1447  // Finally, flag that we have initialised the memory.
1448  m_initialised = true;
1449  }
1450 
1451  // The loop below calculates the stage values and derivatives
1452  for(i = 0; i < m_numstages; i++)
1453  {
1454  if( (i==0) && m_firstStageEqualsOldSolution )
1455  {
1456  for(k = 0; k < m_nvar; k++)
1457  {
1458  Vmath::Vcopy(m_npoints,y_old[0][k],1,m_Y[k],1);
1459  }
1460  m_T = t_old[0];
1461  }
1462  else
1463  {
1464  // The stage values m_Y are a linear combination of:
1465  // 1: the stage derivatives
1466 
1467  if( i != 0 )
1468  {
1469  for(k = 0; k < m_nvar; k++)
1470  {
1471  Vmath::Smul(m_npoints,timestep*A(i,0),m_F[0][k],1,
1472  m_tmp[k],1);
1473 
1474  if(type == eIMEX)
1475  {
1476  Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,0),
1477  m_F_IMEX[0][k],1,
1478  m_tmp[k],1,m_tmp[k],1);
1479  }
1480  }
1481  }
1482  m_T = A(i,0)*timestep;
1483 
1484  for( j = 1; j < i; j++ )
1485  {
1486  for(k = 0; k < m_nvar; k++)
1487  {
1488  Vmath::Svtvp(m_npoints,timestep*A(i,j),m_F[j][k],1,
1489  m_tmp[k],1,m_tmp[k],1);
1490  if(type == eIMEX)
1491  {
1492  Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,j),
1493  m_F_IMEX[j][k],1,
1494  m_tmp[k],1,m_tmp[k],1);
1495  }
1496  }
1497 
1498  m_T += A(i,j)*timestep;
1499  }
1500 
1501  // 2: the imported multi-step solution of the
1502  // previous time level
1503  for(j = 0; j < m_numsteps; j++)
1504  {
1505  for(k = 0; k < m_nvar; k++)
1506  {
1507  Vmath::Svtvp(m_npoints,U(i,j),y_old[j][k],1,
1508  m_tmp[k],1,m_tmp[k],1);
1509  }
1510  m_T += U(i,j)*t_old[j];
1511  }
1512  }
1513 
1514  // Calculate the stage derivative based upon the stage value
1515  if(type == eDiagonallyImplicit)
1516  {
1517  if(m_numstages==1)
1518  {
1519  m_T= t_old[0]+timestep;
1520  }
1521  else
1522  {
1523  m_T= t_old[0];
1524  for(int j=0; j<=i; ++j)
1525  {
1526  m_T += A(i,j)*timestep;
1527  }
1528  }
1529 
1530  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
1531 
1532  for(k = 0; k < m_nvar; k++)
1533  {
1534  Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
1535  Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),m_F[i][k],1,m_F[i][k],1);
1536  }
1537  }
1538  else if(type == eIMEX)
1539  {
1540  if(m_numstages==1)
1541  {
1542  m_T= t_old[0]+timestep;
1543  }
1544  else
1545  {
1546  m_T= t_old[0];
1547  for(int j=0; j<=i; ++j)
1548  {
1549  m_T += A(i,j)*timestep;
1550  }
1551  }
1552 
1553  if(fabs(A(i,i)) > NekConstants::kNekZeroTol)
1554  {
1555  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
1556 
1557  for(k = 0; k < m_nvar; k++)
1558  {
1559  Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
1560  Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),
1561  m_F[i][k],1,m_F[i][k],1);
1562  }
1563  }
1564  op.DoOdeRhs(m_Y, m_F_IMEX[i], m_T);
1565  }
1566  else if( type == eExplicit)
1567  {
1568  // ensure solution is in correct space
1569  op.DoProjection(m_Y,m_Y,m_T);
1570  op.DoOdeRhs(m_Y, m_F[i], m_T);
1571  }
1572  }
1573 
1574  // Next, the solution vector y at the new time level will
1575  // be calculated.
1576  //
1577  // For multi-step methods, this includes updating the
1578  // values of the auxiliary parameters
1579  //
1580  // The loop below calculates the solution at the new time
1581  // level
1582  //
1583  // If last stage equals the new solution, the new solution
1584  // needs not be calculated explicitly but can simply be
1585  // copied. This saves a solve.
1586  int i_start = 0;
1587  if( (m_lastStageEqualsNewSolution) )
1588  {
1589  for(k = 0; k < m_nvar; k++)
1590  {
1591  Vmath::Vcopy(m_npoints,m_Y[k],1,y_new[0][k],1);
1592  }
1593 
1594  if (m_numstages==1 && type == eIMEX)
1595  {
1596  t_new[0] = t_old[0]+timestep;
1597  }
1598  else
1599  {
1600  t_new[0] = B(0,0)*timestep;
1601  for(j = 1; j < m_numstages; j++)
1602  {
1603  t_new[0] += B(0,j)*timestep;
1604  }
1605  for(j = 0; j < m_numsteps; j++)
1606  {
1607  t_new[0] += V(0,j)*t_old[j];
1608  }
1609  i_start = 1;
1610  }
1611  }
1612 
1613  for(i = i_start; i < m_numsteps; i++)
1614  {
1615  // The solution at the new time level is a linear
1616  // combination of:
1617  // 1: the stage derivatives
1618  for(k = 0; k < m_nvar; k++)
1619  {
1620  Vmath::Smul(m_npoints,timestep*B(i,0),m_F[0][k],1,
1621  y_new[i][k],1);
1622 
1623  if(type == eIMEX)
1624  {
1625  Vmath::Svtvp(m_npoints,timestep*B_IMEX(i,0),
1626  m_F_IMEX[0][k],1,y_new[i][k],1,
1627  y_new[i][k],1);
1628  }
1629  }
1630  if(m_numstages != 1 || type != eIMEX)
1631  {
1632  t_new[i] = B(i,0)*timestep;
1633  }
1634 
1635 
1636  for(j = 1; j < m_numstages; j++)
1637  {
1638  for(k = 0; k < m_nvar; k++)
1639  {
1640  Vmath::Svtvp(m_npoints,timestep*B(i,j),m_F[j][k],1,
1641  y_new[i][k],1,y_new[i][k],1);
1642 
1643  if(type == eIMEX)
1644  {
1645  Vmath::Svtvp(m_npoints,timestep*B_IMEX(i,j),
1646  m_F_IMEX[j][k],1,y_new[i][k],1,
1647  y_new[i][k],1);
1648  }
1649  }
1650  if(m_numstages != 1 || type != eIMEX)
1651  {
1652  t_new[i] += B(i,j)*timestep;
1653  }
1654  }
1655 
1656  // 2: the imported multi-step solution of the previous
1657  // time level
1658  for(j = 0; j < m_numsteps; j++)
1659  {
1660  for(k = 0; k < m_nvar; k++)
1661  {
1662  Vmath::Svtvp(m_npoints,V(i,j),y_old[j][k],1,
1663  y_new[i][k],1,y_new[i][k],1);
1664  }
1665  if(m_numstages != 1 || type != eIMEX)
1666  {
1667  t_new[i] += V(i,j)*t_old[j];
1668  }
1669  }
1670  }
1671 
1672  // Ensure that the new solution is projected if necessary
1673  if(type == eExplicit)
1674  {
1675  op.DoProjection(y_new[0],y_new[0],t_new[0]);
1676  }
1677  }
Nektar::LibUtilities::TimeIntegrationScheme::B_IMEX
NekDouble B_IMEX(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:414
Nektar::LibUtilities::TimeIntegrationScheme::U
NekDouble U(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:399
Nektar::LibUtilities::TimeIntegrationScheme::B
NekDouble B(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:394
Nektar::LibUtilities::TimeIntegrationScheme::m_nvar
int m_nvar
bool to identify if array has been initialised
Definition: TimeIntegrationScheme.h:536
Nektar::LibUtilities::TimeIntegrationScheme::m_tmp
DoubleArray m_tmp
Array containing the stage values.
Definition: TimeIntegrationScheme.h:539
Nektar::LibUtilities::TimeIntegrationScheme::V
NekDouble V(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:404
Vmath::Svtvp
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:471
Nektar::LibUtilities::eIMEX
Implicit Explicit General Linear Method.
Definition: TimeIntegrationScheme.h:138
Nektar::LibUtilities::TimeIntegrationScheme::m_Y
DoubleArray m_Y
The size of inner data which is stored for reuse.
Definition: TimeIntegrationScheme.h:538
Nektar::LibUtilities::TimeIntegrationScheme::GetSecondDim
int GetSecondDim(ConstTripleArray &y) const
Definition: TimeIntegrationScheme.h:583
Nektar::LibUtilities::TimeIntegrationScheme::m_F
TripleArray m_F
explicit right hand side of each stage equation
Definition: TimeIntegrationScheme.h:541
Nektar::LibUtilities::eDiagonallyImplicit
Diagonally implicit scheme (e.g. the DIRK schemes)
Definition: TimeIntegrationScheme.h:137
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Vmath::Smul
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:199
Nektar::LibUtilities::TimeIntegrationScheme::A_IMEX
NekDouble A_IMEX(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:409
Nektar::LibUtilities::TimeIntegrationScheme::m_npoints
int m_npoints
The number of variables in integration scheme.
Definition: TimeIntegrationScheme.h:537
Nektar::LibUtilities::TimeIntegrationScheme::CheckTimeIntegrateArguments
bool CheckTimeIntegrateArguments(const NekDouble timestep, ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new, SingleArray &t_new, const TimeIntegrationSchemeOperators &op) const
Definition: TimeIntegrationScheme.cpp:1751
Nektar::LibUtilities::TimeIntegrationScheme::DoubleArray
Array< OneD, Array< OneD, NekDouble > > DoubleArray
Definition: TimeIntegrationScheme.h:362
Nektar::LibUtilities::TimeIntegrationScheme::A
NekDouble A(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:389
Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeType
TimeIntegrationSchemeType GetIntegrationSchemeType() const
Definition: TimeIntegrationScheme.h:384
Vmath::Vsub
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:329
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::LibUtilities::TimeIntegrationScheme::m_T
NekDouble m_T
Used to store the Explicit stage derivative of IMEX schemes.
Definition: TimeIntegrationScheme.h:544
Nektar::LibUtilities::TimeIntegrationScheme::TripleArray
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > TripleArray
Definition: TimeIntegrationScheme.h:360
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:165
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1016
Nektar::LibUtilities::TimeIntegrationScheme::m_F_IMEX
TripleArray m_F_IMEX
Array corresponding to the stage Derivatives.
Definition: TimeIntegrationScheme.h:542
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::U ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 399 of file TimeIntegrationScheme.h.

References m_U.

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

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

Definition at line 404 of file TimeIntegrationScheme.h.

References m_V.

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

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

Definition at line 1040 of file TimeIntegrationScheme.cpp.

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

Referenced by TimeIntegrationScheme().

1045  {
1046  int i;
1047  int j;
1048  int m;
1049  int IMEXdim = A.num_elements();
1050  int dim = A[0].GetRows();
1051 
1052  Array<OneD, TimeIntegrationSchemeType> vertype(IMEXdim,eExplicit);
1053 
1054  for(m = 0; m < IMEXdim; m++)
1055  {
1056  for(i = 0; i < dim; i++)
1057  {
1058  if( fabs(A[m][i][i]) > NekConstants::kNekZeroTol )
1059  {
1060  vertype[m] = eDiagonallyImplicit;
1061  }
1062  }
1063 
1064  for(i = 0; i < dim; i++)
1065  {
1066  for(j = i+1; j < dim; j++)
1067  {
1068  if( fabs(A[m][i][j]) > NekConstants::kNekZeroTol )
1069  {
1070  vertype[m] = eImplicit;
1071  ASSERTL1(false,"Fully Impplicit schemes cannnot be handled by the TimeIntegrationScheme class");
1072  }
1073  }
1074  }
1075  }
1076 
1077  if(IMEXdim == 2)
1078  {
1079  ASSERTL1(B.num_elements()==2,"Coefficient Matrix B should have an implicit and explicit part for IMEX schemes");
1080  if((vertype[0] == eDiagonallyImplicit) &&
1081  (vertype[1] == eExplicit))
1082  {
1083  vertype[0] = eIMEX;
1084  }
1085  else
1086  {
1087  ASSERTL1(false,"This is not a proper IMEX scheme");
1088  }
1089  }
1090 
1091  return (vertype[0] == type);
1092  }
Nektar::LibUtilities::TimeIntegrationScheme::B
NekDouble B(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:394
Nektar::LibUtilities::eIMEX
Implicit Explicit General Linear Method.
Definition: TimeIntegrationScheme.h:138
Nektar::LibUtilities::eDiagonallyImplicit
Diagonally implicit scheme (e.g. the DIRK schemes)
Definition: TimeIntegrationScheme.h:137
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Nektar::LibUtilities::TimeIntegrationScheme::A
NekDouble A(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:389
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:165

Friends And Related Function Documentation

template<typename >
friend class Nektar::MemoryManager
friend

Time at the different stages.

Definition at line 546 of file TimeIntegrationScheme.h.

TimeIntegrationSchemeManagerT& TimeIntegrationSchemeManager ( void  )
friend

Definition at line 45 of file TimeIntegrationScheme.cpp.

46  {
47  TimeIntegrationSchemeManagerT& m = Loki::SingletonHolder<TimeIntegrationSchemeManagerT>::Instance();
48  m.RegisterGlobalCreator(TimeIntegrationScheme::Create);
49  return m;
50  }
Nektar::LibUtilities::TimeIntegrationSchemeManagerT
NekManager< TimeIntegrationSchemeKey, TimeIntegrationScheme, TimeIntegrationSchemeKey::opLess > TimeIntegrationSchemeManagerT
Definition: TimeIntegrationScheme.h:348
Nektar::LibUtilities::TimeIntegrationScheme::Create
static boost::shared_ptr< TimeIntegrationScheme > Create(const TimeIntegrationSchemeKey &key)
Definition: TimeIntegrationScheme.cpp:143

Member Data Documentation

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

Definition at line 529 of file TimeIntegrationScheme.h.

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

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

Definition at line 530 of file TimeIntegrationScheme.h.

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

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

explicit right hand side of each stage equation

Definition at line 541 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Array corresponding to the stage Derivatives.

Definition at line 542 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 511 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate(), and TimeIntegrationScheme().

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

Definition at line 535 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 512 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate(), and TimeIntegrationScheme().

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

The number of variables in integration scheme.

Definition at line 537 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 516 of file TimeIntegrationScheme.h.

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

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

Definition at line 514 of file TimeIntegrationScheme.h.

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

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

Definition at line 509 of file TimeIntegrationScheme.h.

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

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

Definition at line 508 of file TimeIntegrationScheme.h.

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

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

bool to identify if array has been initialised

Definition at line 536 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 506 of file TimeIntegrationScheme.h.

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

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

Definition at line 507 of file TimeIntegrationScheme.h.

Referenced by GetIntegrationSchemeType(), and TimeIntegrationScheme().

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

Used to store the Explicit stage derivative of IMEX schemes.

Definition at line 544 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 518 of file TimeIntegrationScheme.h.

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

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

Array containing the stage values.

Definition at line 539 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 531 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme(), and U().

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

Definition at line 532 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme(), and V().

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

The size of inner data which is stored for reuse.

Definition at line 538 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

Generated on Tue Dec 8 2015 23:41:37 for Nektar++ by   doxygen 1.8.9.1