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

#include <TimeIntegrationScheme.h>

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

Public Types

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

Public Member Functions

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

Member Typedef Documentation

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

Definition at line 367 of file TimeIntegrationScheme.h.

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

Definition at line 369 of file TimeIntegrationScheme.h.

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

Definition at line 365 of file TimeIntegrationScheme.h.

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

Definition at line 368 of file TimeIntegrationScheme.h.

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

Definition at line 371 of file TimeIntegrationScheme.h.

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

Definition at line 372 of file TimeIntegrationScheme.h.

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

Definition at line 370 of file TimeIntegrationScheme.h.

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

Definition at line 366 of file TimeIntegrationScheme.h.

Constructor & Destructor Documentation

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

Definition at line 376 of file TimeIntegrationScheme.h.

377  {
378  }
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, Nektar::LibUtilities::eRungeKutta2_ImprovedEuler, Nektar::LibUtilities::eRungeKutta2_SSP, Nektar::LibUtilities::eRungeKutta3_SSP, Nektar::LibUtilities::eRungeKutta4, 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  case eRungeKutta2:
442  {
443  m_numsteps = 1;
444  m_numstages = 2;
445 
446  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
447  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
448 
449  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
450  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
451  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
452  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
453 
454  m_A[0][1][0] = 0.5;
455  m_B[0][0][1] = 1.0;
456 
457  m_schemeType = eExplicit;
458  m_numMultiStepValues = 1;
459  m_numMultiStepDerivs = 0;
460  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
461  m_timeLevelOffset[0] = 0;
462  }
463  break;
464  case eRungeKutta2_ImprovedEuler:
465  case eRungeKutta2_SSP:
466  {
467  m_numsteps = 1;
468  m_numstages = 2;
469 
470  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
471  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
472 
473  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
474  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
475  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
476  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
477 
478  m_A[0][1][0] = 1.0;
479 
480  m_B[0][0][0] = 0.5;
481  m_B[0][0][1] = 0.5;
482 
483  m_schemeType = eExplicit;
484  m_numMultiStepValues = 1;
485  m_numMultiStepDerivs = 0;
486  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
487  m_timeLevelOffset[0] = 0;
488  }
489  break;
490  case eRungeKutta3_SSP:
491  {
492  m_numsteps = 1;
493  m_numstages = 3;
494 
495  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
496  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
497 
498  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
499  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
500  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
501  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
502 
503  m_A[0][1][0] = 1.0;
504  m_A[0][2][0] = 0.25;
505  m_A[0][2][1] = 0.25;
506 
507  m_B[0][0][0] = 1.0/6.0;
508  m_B[0][0][1] = 1.0/6.0;
509  m_B[0][0][2] = 2.0/3.0;
510 
511  m_schemeType = eExplicit;
512  m_numMultiStepValues = 1;
513  m_numMultiStepDerivs = 0;
514  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
515  m_timeLevelOffset[0] = 0;
516  }
517  break;
518  case eClassicalRungeKutta4:
519  case eRungeKutta4:
520  {
521  m_numsteps = 1;
522  m_numstages = 4;
523 
524  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
525  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
526 
527  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
528  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
529  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
530  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
531 
532  m_A[0][1][0] = 0.5;
533  m_A[0][2][1] = 0.5;
534  m_A[0][3][2] = 1.0;
535 
536  m_B[0][0][0] = 1.0/6.0;
537  m_B[0][0][1] = 1.0/3.0;
538  m_B[0][0][2] = 1.0/3.0;
539  m_B[0][0][3] = 1.0/6.0;
540 
541  m_schemeType = eExplicit;
542  m_numMultiStepValues = 1;
543  m_numMultiStepDerivs = 0;
544  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
545  m_timeLevelOffset[0] = 0;
546  }
547  break;
548  case eDIRKOrder2:
549  {
550  m_numsteps = 1;
551  m_numstages = 2;
552 
553  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
554  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
555 
556  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
557  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
558  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
559  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
560 
561  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
562 
563  m_A[0][0][0] = lambda;
564  m_A[0][1][0] = 1.0 - lambda;
565  m_A[0][1][1] = lambda;
566 
567  m_B[0][0][0] = 1.0 - lambda;
568  m_B[0][0][1] = lambda;
569 
570  m_schemeType = eDiagonallyImplicit;
571  m_numMultiStepValues = 1;
572  m_numMultiStepDerivs = 0;
573  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
574  m_timeLevelOffset[0] = 0;
575  }
576  break;
577  case eDIRKOrder3:
578  {
579  m_numsteps = 1;
580  m_numstages = 3;
581 
582  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
583  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
584 
585  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
586  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
587  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
588  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
589 
590  NekDouble lambda = 0.4358665215;
591 
592  m_A[0][0][0] = lambda;
593  m_A[0][1][0] = 0.5 * (1.0 - lambda);
594  m_A[0][2][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
595  m_A[0][1][1] = lambda;
596  m_A[0][2][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
597  m_A[0][2][2] = lambda;
598 
599  m_B[0][0][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
600  m_B[0][0][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
601  m_B[0][0][2] = lambda;
602 
603  m_schemeType = eDiagonallyImplicit;
604  m_numMultiStepValues = 1;
605  m_numMultiStepDerivs = 0;
606  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
607  m_timeLevelOffset[0] = 0;
608  }
609  break;
610  case eIMEXdirk_2_3_2:
611  {
612  m_numsteps = 1;
613  m_numstages = 3;
614 
615  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
616  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
617 
618  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
619  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
620  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
621  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
622  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
623  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
624 
625  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
626  NekDouble delta = -2.0*sqrt(2.0)/3.0;
627 
628  m_A[0][1][1] = lambda;
629  m_A[0][2][1] = 1.0 - lambda;
630  m_A[0][2][2] = lambda;
631 
632  m_B[0][0][1] = 1.0 - lambda;
633  m_B[0][0][2] = lambda;
634 
635  m_A[1][1][0] = lambda;
636  m_A[1][2][0] = delta;
637  m_A[1][2][1] = 1.0 - delta;
638 
639  m_B[1][0][1] = 1.0 - lambda;
640  m_B[1][0][2] = lambda;
641 
642  m_schemeType = eIMEX;
643  m_numMultiStepValues = 1;
644  m_numMultiStepDerivs = 0;
645  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
646  m_timeLevelOffset[0] = 0;
647  }
648  break;
649  case eIMEXdirk_3_4_3:
650  {
651  m_numsteps = 1;
652  m_numstages = 4;
653 
654  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
655  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
656 
657  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
658  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
659  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
660  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
661  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
662  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
663 
664  NekDouble lambda = 0.4358665215;
665 
666  m_A[0][1][1] = lambda;
667  m_A[0][2][1] = 0.5 * (1.0 - lambda);
668  m_A[0][3][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
669  m_A[0][2][2] = lambda;
670  m_A[0][3][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
671  m_A[0][3][3] = lambda;
672 
673  m_B[0][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
674  m_B[0][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
675  m_B[0][0][3] = lambda;
676 
677  m_A[1][1][0] = 0.4358665215;
678  m_A[1][2][0] = 0.3212788860;
679  m_A[1][2][1] = 0.3966543747;
680  m_A[1][3][0] =-0.105858296;
681  m_A[1][3][1] = 0.5529291479;
682  m_A[1][3][2] = 0.5529291479;
683 
684  m_B[1][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
685  m_B[1][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
686  m_B[1][0][3] = lambda;
687 
688  m_schemeType = eIMEX;
689  m_numMultiStepValues = 1;
690  m_numMultiStepDerivs = 0;
691  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
692  m_timeLevelOffset[0] = 0;
693  }
694  break;
695  case eCNAB:
696  {
697  NekDouble secondth = 1.0/2.0;
698  m_numsteps = 4;
699  m_numstages = 1;
700 
701  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
702  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
703 
704  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,secondth);
705  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
706  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
707  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
708  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,secondth);
709  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
710 
711  m_B[0][0][0] = secondth;
712  m_B[0][1][0] = 1.0;
713  m_B[1][2][0] = 1.0;
714  m_U[0][0] = 2*secondth;
715  m_U[0][2] = 3*secondth;
716  m_U[0][3] = -1*secondth;
717 
718  m_V[0][0] = 2*secondth;
719  m_V[0][1] = secondth;
720  m_V[0][2] = 3*secondth;
721  m_V[0][3] = -1*secondth;
722  m_V[3][2] = 1.0;
723 
724  m_schemeType = eIMEX;
725  m_numMultiStepValues = 1;
726  m_numMultiStepDerivs = 3;
727  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
728  m_timeLevelOffset[0] = 0;
729  m_timeLevelOffset[1] = 0;
730  m_timeLevelOffset[2] = 0;
731  m_timeLevelOffset[3] = 1;
732  }
733  break;
734  case eIMEXGear:
735  {
736  NekDouble twothirdth = 2.0/3.0;
737  m_numsteps = 2;
738  m_numstages = 1;
739 
740  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
741  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
742 
743  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
744  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
745  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
746  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
747  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,twothirdth);
748  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
749 
750  m_B[0][0][0] = twothirdth;
751  m_B[1][0][0] = twothirdth;
752  m_U[0][0] = 2*twothirdth;
753  m_U[0][1] = -0.5*twothirdth;
754 
755  m_V[0][0] = 2*twothirdth;
756  m_V[0][1] = -0.5*twothirdth;
757  m_V[1][0] = 1.0;
758 
759  m_schemeType = eIMEX;
760  m_numMultiStepValues = 2;
761  m_numMultiStepDerivs = 0;
762  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
763  m_timeLevelOffset[0] = 0;
764  m_timeLevelOffset[1] = 1;
765  }
766  break;
767  case eMCNAB:
768  {
769  NekDouble sixthx = 9.0/16.0;
770  m_numsteps = 5;
771  m_numstages = 1;
772 
773  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
774  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
775 
776  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,sixthx);
777  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
778  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
779  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
780  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
781  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
782 
783  m_B[0][0][0] = sixthx;
784  m_B[0][1][0] = 1.0;
785  m_B[1][3][0] = 1.0;
786  m_U[0][0] = 1.0;
787  m_U[0][1] = 6.0/16.0;
788  m_U[0][2] = 1.0/16.0;
789  m_U[0][3] = 1.5;
790  m_U[0][4] = -0.5;
791 
792  m_V[0][0] = 1.0;
793  m_V[0][1] = 6.0/16.0;
794  m_V[0][2] = 1.0/16.0;
795  m_V[0][3] = 1.5;
796  m_V[0][4] = -0.5;
797  m_V[2][1] = 1.0;
798  m_V[4][3] = 1.0;
799 
800  m_schemeType = eIMEX;
801  m_numMultiStepValues = 1;
802  m_numMultiStepDerivs = 4;
803  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
804  m_timeLevelOffset[0] = 0;
805  m_timeLevelOffset[1] = 0;
806  m_timeLevelOffset[2] = 1;
807  m_timeLevelOffset[3] = 0;
808  m_timeLevelOffset[4] = 1;
809  }
810  break;
811  case eIMEXdirk_2_2_2:
812  {
813  m_numsteps = 1;
814  m_numstages = 3;
815 
816  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
817  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
818 
819  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
820  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
821  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
822  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
823  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
824  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
825 
826  NekDouble glambda = 0.788675134594813;
827  NekDouble gdelta = 0.366025403784439;
828 
829  m_A[0][1][1] = glambda;
830  m_A[0][2][1] = 1.0 - glambda;
831  m_A[0][2][2] = glambda;
832 
833  m_B[0][0][1] = 1.0 - glambda;
834  m_B[0][0][2] = glambda;
835 
836  m_A[1][1][0] = glambda;
837  m_A[1][2][0] = gdelta;
838  m_A[1][2][1] = 1.0 - gdelta;
839 
840  m_B[1][0][0] = gdelta;
841  m_B[1][0][1] = 1.0 - gdelta;
842 
843  m_schemeType = eIMEX;
844  m_numMultiStepValues = 1;
845  m_numMultiStepDerivs = 0;
846  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
847  m_timeLevelOffset[0] = 0;
848  }
849  break;
850  case eIMEXdirk_2_3_3:
851  {
852  m_numsteps = 1;
853  m_numstages = 3;
854 
855  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
856  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
857 
858  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
859  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
860  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
861  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
862  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
863  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
864 
865  NekDouble glambda = 0.788675134594813;
866 
867  m_A[0][1][1] = glambda;
868  m_A[0][2][1] = 1.0 - 2.0*glambda;
869  m_A[0][2][2] = glambda;
870 
871  m_B[0][0][1] = 0.5;
872  m_B[0][0][2] = 0.5;
873 
874  m_A[1][1][0] = glambda;
875  m_A[1][2][0] = glambda - 1.0;
876  m_A[1][2][1] = 2.0*(1-glambda);
877 
878  m_B[1][0][1] = 0.5;
879  m_B[1][0][2] = 0.5;
880 
881  m_schemeType = eIMEX;
882  m_numMultiStepValues = 1;
883  m_numMultiStepDerivs = 0;
884  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
885  m_timeLevelOffset[0] = 0;
886  }
887  break;
888  case eIMEXdirk_1_1_1:
889  {
890  m_numsteps = 1;
891  m_numstages = 2;
892 
893  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
894  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
895 
896  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
897  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
898  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
899  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
900  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
901  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
902 
903  m_A[0][1][1] = 1.0;
904 
905  m_B[0][0][1] = 1.0;
906 
907  m_A[1][1][0] = 1.0;
908 
909  m_B[1][0][0] = 1.0;
910 
911  m_schemeType = eIMEX;
912  m_numMultiStepValues = 1;
913  m_numMultiStepDerivs = 0;
914  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
915  m_timeLevelOffset[0] = 0;
916  }
917  break;
918  case eIMEXdirk_1_2_1:
919  {
920  m_numsteps = 1;
921  m_numstages = 2;
922 
923  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
924  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
925 
926  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
927  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
928  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
929  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
930  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
931  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
932 
933 
934  m_A[0][1][1] = 1.0;
935 
936  m_B[0][0][1] = 1.0;
937 
938  m_A[1][1][0] = 1.0;
939 
940  m_B[1][0][1] = 1.0;
941 
942  m_schemeType = eIMEX;
943  m_numMultiStepValues = 1;
944  m_numMultiStepDerivs = 0;
945  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
946  m_timeLevelOffset[0] = 0;
947  }
948  break;
949  case eIMEXdirk_1_2_2:
950  {
951  m_numsteps = 1;
952  m_numstages = 2;
953 
954  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
955  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
956 
957  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
958  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
959  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
960  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
961  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
962  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
963 
964  m_A[0][1][1] = 1.0/2.0;
965 
966  m_B[0][0][1] = 1.0;
967 
968  m_A[1][1][0] = 1.0/2.0;
969 
970  m_B[1][0][1] = 1.0;
971 
972  m_schemeType = eIMEX;
973  m_numMultiStepValues = 1;
974  m_numMultiStepDerivs = 0;
975  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
976  m_timeLevelOffset[0] = 0;
977  }
978  break;
979  case eIMEXdirk_4_4_3:
980  {
981  m_numsteps = 1;
982  m_numstages = 5;
983 
984  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
985  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
986 
987  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
988  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
989  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
990  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
991  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
992  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
993 
994  m_A[0][1][1] = 1.0/2.0;
995  m_A[0][2][1] = 1.0/6.0;
996  m_A[0][2][2] = 1.0/2.0;
997  m_A[0][3][1] = -1.0/2.0;
998  m_A[0][3][2] = 1.0/2.0;
999  m_A[0][3][3] = 1.0/2.0;
1000  m_A[0][4][1] = 3.0/2.0;
1001  m_A[0][4][2] = -3.0/2.0;
1002  m_A[0][4][3] = 1.0/2.0;
1003  m_A[0][4][4] = 1.0/2.0;
1004 
1005  m_B[0][0][1] = 3.0/2.0;
1006  m_B[0][0][2] = -3.0/2.0;
1007  m_B[0][0][3] = 1.0/2.0;
1008  m_B[0][0][4] = 1.0/2.0;
1009 
1010  m_A[1][1][0] = 1.0/2.0;
1011  m_A[1][2][0] = 11.0/18.0;
1012  m_A[1][2][1] = 1.0/18.0;
1013  m_A[1][3][0] = 5.0/6.0;
1014  m_A[1][3][1] = -5.0/6.0;
1015  m_A[1][3][2] = 1.0/2.0;
1016  m_A[1][4][0] = 1.0/4.0;
1017  m_A[1][4][1] = 7.0/4.0;
1018  m_A[1][4][2] = 3.0/4.0;
1019  m_A[1][4][3] = -7.0/4.0;
1020 
1021  m_B[1][0][0] = 1.0/4.0;
1022  m_B[1][0][1] = 7.0/4.0;
1023  m_B[1][0][2] = 3.0/4.0;
1024  m_B[1][0][3] = -7.0/4.0;
1025 
1026  m_schemeType = eIMEX;
1027  m_numMultiStepValues = 1;
1028  m_numMultiStepDerivs = 0;
1029  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1030  m_timeLevelOffset[0] = 0;
1031  }
1032  break;
1033  default:
1034  {
1035  NEKERROR(ErrorUtil::efatal,"Invalid Time Integration Scheme");
1036  }
1037  }
1038 
1039  m_firstStageEqualsOldSolution = CheckIfFirstStageEqualsOldSolution(m_A,m_B,m_U,m_V);
1040  m_lastStageEqualsNewSolution = CheckIfLastStageEqualsNewSolution(m_A,m_B,m_U,m_V);
1041 
1042  ASSERTL1(VerifyIntegrationSchemeType(m_schemeType,m_A,m_B,m_U,m_V),
1043  "Time integration scheme coefficients do not match its type");
1044  }
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
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:536
Nektar::LibUtilities::eClassicalRungeKutta4
Runge-Kutta multi-stage scheme 4th order explicit (old name)
Definition: TimeIntegrationScheme.h:74
Nektar::LibUtilities::eIMEXdirk_1_2_2
Implicit-Explicit Midpoint IMEX DIRK(1,2,2)
Definition: TimeIntegrationScheme.h:93
Nektar::LibUtilities::eRungeKutta4
Classical RungeKutta4 method (new name for eClassicalRungeKutta4)
Definition: TimeIntegrationScheme.h:75
Nektar::LibUtilities::eIMEXdirk_4_4_3
L-stable, four stage, third order IMEX DIRK(4,4,3)
Definition: TimeIntegrationScheme.h:98
Nektar::LibUtilities::eIMEXdirk_1_2_1
Forward-Backward Euler IMEX DIRK(1,2,1)
Definition: TimeIntegrationScheme.h:92
Nektar::LibUtilities::eIMEX
Implicit Explicit General Linear Method.
Definition: TimeIntegrationScheme.h:144
Nektar::LibUtilities::eRungeKutta3_SSP
Nonlinear SSP RungeKutta3 explicit.
Definition: TimeIntegrationScheme.h:76
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:88
Nektar::LibUtilities::eIMEXOrder2
IMEX 2nd order scheme using Backward Different Formula & Extrapolation.
Definition: TimeIntegrationScheme.h:82
Nektar::LibUtilities::eDiagonallyImplicit
Diagonally implicit scheme (e.g. the DIRK schemes)
Definition: TimeIntegrationScheme.h:143
Nektar::LibUtilities::eRungeKutta2
Classical RungeKutta2 method (new name for eMidpoint)
Definition: TimeIntegrationScheme.h:85
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:1687
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:97
Nektar::LibUtilities::eRungeKutta2_SSP
Nonlinear SSP RungeKutta2 explicit (surrogate for eRungeKutta2_ImprovedEuler)
Definition: TimeIntegrationScheme.h:78
Nektar::LibUtilities::eRungeKutta2_ImprovedEuler
Improved RungeKutta2 explicit (old name meaning Heun's method)
Definition: TimeIntegrationScheme.h:77
Nektar::LibUtilities::eIMEXdirk_1_1_1
Forward-Backward Euler IMEX DIRK(1,1,1)
Definition: TimeIntegrationScheme.h:91
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
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:94
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:1048
Nektar::LibUtilities::eDIRKOrder2
Diagonally Implicit Runge Kutta scheme of order 3.
Definition: TimeIntegrationScheme.h:86
Nektar::LibUtilities::eIMEXOrder3
IMEX 3rd order scheme using Backward Different Formula & Extrapolation.
Definition: TimeIntegrationScheme.h:83
Nektar::LibUtilities::eIMEXdirk_2_3_2
L-stable, three stage, third order IMEX DIRK(3,4,3)
Definition: TimeIntegrationScheme.h:95
Nektar::LibUtilities::eDIRKOrder3
Diagonally Implicit Runge Kutta scheme of order 3.
Definition: TimeIntegrationScheme.h:87
Nektar::LibUtilities::eIMEXdirk_2_3_3
L-stable, two stage, third order IMEX DIRK(2,3,3)
Definition: TimeIntegrationScheme.h:96
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:1725
Nektar::LibUtilities::TimeIntegrationScheme::m_A
Array< OneD, Array< TwoD, NekDouble > > m_A
Definition: TimeIntegrationScheme.h:535
Nektar::LibUtilities::eMCNAB
Modified Crank-Nicolson/Adams-Bashforth Order 2 (MCNAB)
Definition: TimeIntegrationScheme.h:90
Nektar::LibUtilities::eIMEXOrder1
IMEX 1st order scheme using Euler Backwards/Euler Forwards.
Definition: TimeIntegrationScheme.h:81
Nektar::LibUtilities::eMidpoint
midpoint method (old name)
Definition: TimeIntegrationScheme.h:84
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrationScheme ( )
inlineprivate

Definition at line 560 of file TimeIntegrationScheme.h.

References ErrorUtil::efatal, and NEKERROR.

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

Definition at line 565 of file TimeIntegrationScheme.h.

References ErrorUtil::efatal, and NEKERROR.

565  :m_schemeKey(NullTimeIntegrationSchemeKey)
566  {
567  NEKERROR(ErrorUtil::efatal,"Copy Constructor for the TimeIntegrationScheme class should not be called");
568  }
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
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 395 of file TimeIntegrationScheme.h.

References m_A.

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

396  {
397  return m_A[0][i][j];
398  }
Nektar::LibUtilities::TimeIntegrationScheme::m_A
Array< OneD, Array< TwoD, NekDouble > > m_A
Definition: TimeIntegrationScheme.h:535
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::A_IMEX ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 415 of file TimeIntegrationScheme.h.

References m_A.

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

416  {
417  return m_A[1][i][j];
418  }
Nektar::LibUtilities::TimeIntegrationScheme::m_A
Array< OneD, Array< TwoD, NekDouble > > m_A
Definition: TimeIntegrationScheme.h:535
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::B ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 400 of file TimeIntegrationScheme.h.

References m_B.

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

401  {
402  return m_B[0][i][j];
403  }
Nektar::LibUtilities::TimeIntegrationScheme::m_B
Array< OneD, Array< TwoD, NekDouble > > m_B
Definition: TimeIntegrationScheme.h:536
NekDouble Nektar::LibUtilities::TimeIntegrationScheme::B_IMEX ( const unsigned int  i,
const unsigned int  j 
) const
inline

Definition at line 420 of file TimeIntegrationScheme.h.

References m_B.

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

421  {
422  return m_B[1][i][j];
423  }
Nektar::LibUtilities::TimeIntegrationScheme::m_B
Array< OneD, Array< TwoD, NekDouble > > m_B
Definition: TimeIntegrationScheme.h:536
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 1687 of file TimeIntegrationScheme.cpp.

References Nektar::NekConstants::kNekZeroTol, m_numstages, and m_numsteps.

Referenced by TimeIntegrationScheme().

1691  {
1692  int i,m;
1693  // First stage equals old solution if:
1694  // 1. the first row of the coefficient matrix A consists of zeros
1695  // 2. U[0][0] is equal to one and all other first row entries of U are zero
1696 
1697  // 1. Check first condition
1698  for(m = 0; m < A.num_elements(); m++)
1699  {
1700  for(i = 0; i < m_numstages; i++)
1701  {
1702  if( fabs(A[m][0][i]) > NekConstants::kNekZeroTol )
1703  {
1704  return false;
1705  }
1706  }
1707  }
1708 
1709  // 2. Check second condition
1710  if( fabs(U[0][0] - 1.0) > NekConstants::kNekZeroTol )
1711  {
1712  return false;
1713  }
1714  for(i = 1; i < m_numsteps; i++)
1715  {
1716  if( fabs(U[0][i]) > NekConstants::kNekZeroTol )
1717  {
1718  return false;
1719  }
1720  }
1721 
1722  return true;
1723  }
Nektar::LibUtilities::TimeIntegrationScheme::U
NekDouble U(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:405
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
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 1725 of file TimeIntegrationScheme.cpp.

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

Referenced by TimeIntegrationScheme().

1729  {
1730  int i,m;
1731  // Last stage equals new solution if:
1732  // 1. the last row of the coefficient matrix A is equal to the first row of matrix B
1733  // 2. the last row of the coefficient matrix U is equal to the first row of matrix V
1734 
1735  // 1. Check first condition
1736  for(m = 0; m < A.num_elements(); m++)
1737  {
1738  for(i = 0; i < m_numstages; i++)
1739  {
1740  if( fabs(A[m][m_numstages-1][i]-B[m][0][i]) > NekConstants::kNekZeroTol )
1741  {
1742  return false;
1743  }
1744  }
1745  }
1746 
1747  // 2. Check second condition
1748  for(i = 0; i < m_numsteps; i++)
1749  {
1750  if( fabs(U[m_numstages-1][i]-V[0][i]) > NekConstants::kNekZeroTol )
1751  {
1752  return false;
1753  }
1754  }
1755 
1756  return true;
1757  }
Nektar::LibUtilities::TimeIntegrationScheme::U
NekDouble U(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:405
Nektar::LibUtilities::TimeIntegrationScheme::B
NekDouble B(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:400
Nektar::LibUtilities::TimeIntegrationScheme::V
NekDouble V(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:410
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
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 1759 of file TimeIntegrationScheme.cpp.

References ASSERTL1, and m_numsteps.

Referenced by TimeIntegrate().

1765  {
1766  // Check if arrays are all of consistent size
1767  ASSERTL1(y_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1768  ASSERTL1(y_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1769 
1770  ASSERTL1(y_old[0]. num_elements()==y_new[0]. num_elements(),"Non-matching number of variables.");
1771  ASSERTL1(y_old[0][0].num_elements()==y_new[0][0].num_elements(),"Non-matching number of coefficients.");
1772 
1773  ASSERTL1(t_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1774  ASSERTL1(t_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1775 
1776  return true;
1777  }
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
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 584 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 385 of file TimeIntegrationScheme.h.

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

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

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

Definition at line 380 of file TimeIntegrationScheme.h.

References m_schemeKey.

Referenced by TimeIntegrate().

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

Definition at line 390 of file TimeIntegrationScheme.h.

References m_schemeType.

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

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

Definition at line 440 of file TimeIntegrationScheme.h.

References m_numMultiStepDerivs.

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

Definition at line 435 of file TimeIntegrationScheme.h.

References m_numMultiStepValues.

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

Definition at line 430 of file TimeIntegrationScheme.h.

References m_numstages.

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

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

Definition at line 425 of file TimeIntegrationScheme.h.

References m_numsteps.

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

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

Definition at line 589 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 505 of file TimeIntegrationScheme.h.

References m_timeLevelOffset.

506  {
507  return m_timeLevelOffset;
508  }
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 1103 of file TimeIntegrationScheme.cpp.

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

1107  {
1108  // create a TimeIntegrationSolution object based upon the
1109  // initial value. Initialise all other multi-step values
1110  // and derivatives to zero
1111  TimeIntegrationSolutionSharedPtr y_out =
1112  MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(m_schemeKey,y_0,time,timestep);
1113 
1114  // calculate the initial derivative, if is part of the
1115  // solution vector of the current scheme
1116  if(m_numMultiStepDerivs)
1117  {
1118  if(m_timeLevelOffset[m_numMultiStepValues] == 0)
1119  {
1120  int i;
1121  int nvar = y_0.num_elements();
1122  int npoints = y_0[0].num_elements();
1123  DoubleArray f_y_0(nvar);
1124  for(i = 0; i < nvar; i++)
1125  {
1126  f_y_0[i] = Array<OneD,NekDouble>(npoints);
1127  }
1128  // calculate the derivative of the initial value
1129  op.DoOdeRhs(y_0,f_y_0,time);
1130 
1131  // multiply by the step size
1132  for(i = 0; i < nvar; i++)
1133  {
1134  Blas::Dscal(npoints,timestep,f_y_0[i].get(),1);
1135  }
1136  y_out->SetDerivative(0,f_y_0,timestep);
1137  }
1138  }
1139 
1140  return y_out;
1141  }
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:368
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 1144 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().

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

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

Definition at line 405 of file TimeIntegrationScheme.h.

References m_U.

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

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

Definition at line 410 of file TimeIntegrationScheme.h.

References m_V.

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

411  {
412  return m_V[i][j];
413  }
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 1048 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().

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

Friends And Related Function Documentation

template<typename >
friend class Nektar::MemoryManager
friend

Time at the different stages.

Definition at line 552 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:354
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 535 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 536 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 547 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Array corresponding to the stage Derivatives.

Definition at line 548 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 517 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate(), and TimeIntegrationScheme().

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

Definition at line 541 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 518 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 543 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 522 of file TimeIntegrationScheme.h.

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

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

Definition at line 520 of file TimeIntegrationScheme.h.

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

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

Definition at line 515 of file TimeIntegrationScheme.h.

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

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

Definition at line 514 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 542 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 512 of file TimeIntegrationScheme.h.

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

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

Definition at line 513 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 550 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 524 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 545 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

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

Definition at line 537 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrationScheme(), and U().

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

Definition at line 538 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 544 of file TimeIntegrationScheme.h.

Referenced by TimeIntegrate().

Generated on Wed Mar 8 2017 12:29:48 for Nektar++ by   doxygen 1.8.8