Nektar++
TimeIntegrationScheme.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File TimeIntegrationScheme.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: implementation of time integration key class
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 #include <LibUtilities/TimeIntegration/TimeIntegrationScheme.h>
35 #include <LibUtilities/BasicUtils/VmathArray.hpp>
36 #include <LibUtilities/BasicConst/NektarUnivConsts.hpp>
37 #include <iostream>
38 #include <math.h>
39 
40 namespace Nektar
41 {
42  namespace LibUtilities
43  {
44  TimeIntegrationSchemeManagerT &TimeIntegrationSchemeManager(void)
45  {
46  static TimeIntegrationSchemeManagerT instance;
47  instance.RegisterGlobalCreator(TimeIntegrationScheme::Create);
48  return instance;
49  }
50 
51 
52  bool operator==(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs)
53  {
54  return (lhs.m_method == rhs.m_method);
55  }
56 
57  bool operator<(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs)
58  {
59  return (lhs.m_method < rhs.m_method);
60  }
61 
62  bool TimeIntegrationSchemeKey::opLess::operator()(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs) const
63  {
64  return (lhs.m_method < rhs.m_method);
65  }
66 
67  std::ostream& operator<<(std::ostream& os, const TimeIntegrationSchemeKey& rhs)
68  {
69  os << "Time Integration Scheme: " << TimeIntegrationMethodMap[rhs.GetIntegrationMethod()] << std::endl;
70 
71  return os;
72  }
73 
74  TimeIntegrationSolution::TimeIntegrationSolution(const TimeIntegrationSchemeKey &key,
75  const DoubleArray& y,
76  const NekDouble time,
77  const NekDouble timestep):
78  m_scheme(TimeIntegrationSchemeManager()[key]),
79  m_solVector(m_scheme->GetNsteps()),
80  m_t(m_scheme->GetNsteps())
81  {
82  m_solVector[0] = y;
83  m_t[0] = time;
84 
85  int nsteps = m_scheme->GetNsteps();
86  int nvar = y.num_elements();
87  int npoints = y[0].num_elements();
88  int nMultiStepVals = m_scheme->GetNmultiStepValues();
89  const Array<OneD, const unsigned int>& timeLevels = GetTimeLevelOffset();
90  for(int i = 1; i < nsteps; i++)
91  {
92  m_solVector[i] = Array<OneD, Array<OneD, NekDouble> >(nvar);
93  for(int j = 0; j < nvar; j++)
94  {
95  m_solVector[i][j] = Array<OneD,NekDouble>(npoints,0.0);
96  }
97  if(i < nMultiStepVals)
98  {
99  m_t[i] = time - i*timestep*timeLevels[i];
100  }
101  else
102  {
103  m_t[i] = timestep;
104  }
105  }
106  }
107 
108  TimeIntegrationSolution::TimeIntegrationSolution(const TimeIntegrationSchemeKey &key,
109  const TripleArray& y,
110  const Array<OneD, NekDouble>& t):
111  m_scheme(TimeIntegrationSchemeManager()[key]),
112  m_solVector(y),
113  m_t(t)
114  {
115  ASSERTL1(y.num_elements()==m_scheme->GetNsteps(),"Amount of Entries does not match number of (multi-) steps");
116  }
117 
118  TimeIntegrationSolution::TimeIntegrationSolution(const TimeIntegrationSchemeKey &key,
119  unsigned int nvar,
120  unsigned int npoints):
121  m_scheme(TimeIntegrationSchemeManager()[key]),
122  m_solVector(m_scheme->GetNsteps()),
123  m_t(m_scheme->GetNsteps())
124  {
125  for(int i = 0; i < m_scheme->GetNsteps(); i++)
126  {
127  m_solVector[i] = Array<OneD, Array<OneD, NekDouble> >(nvar);
128  for(int j = 0; j < nvar; j++)
129  {
130  m_solVector[i][j] = Array<OneD,NekDouble>(npoints);
131  }
132  }
133  }
134 
135  TimeIntegrationSolution::TimeIntegrationSolution(const TimeIntegrationSchemeKey &key):
136  m_scheme(TimeIntegrationSchemeManager()[key]),
137  m_solVector(m_scheme->GetNsteps()),
138  m_t(m_scheme->GetNsteps())
139  {
140  }
141 
142  TimeIntegrationSchemeSharedPtr TimeIntegrationScheme::Create(const TimeIntegrationSchemeKey &key)
143  {
144  TimeIntegrationSchemeSharedPtr returnval(
145  MemoryManager<TimeIntegrationScheme>::AllocateSharedPtr(key));
146  return returnval;
147  }
148 
149  TimeIntegrationScheme::TimeIntegrationScheme(const TimeIntegrationSchemeKey &key):
150  m_schemeKey(key),
151  m_initialised(false)
152  {
153  switch(key.GetIntegrationMethod())
154  {
155  case eForwardEuler:
156  case eAdamsBashforthOrder1:
157  {
158  m_numsteps = 1;
159  m_numstages = 1;
160 
161  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
162  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
163 
164  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
165  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
166  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
167  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 1.0);
168 
169  m_schemeType = eExplicit;
170  m_numMultiStepValues = 1;
171  m_numMultiStepDerivs = 0;
172  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
173  m_timeLevelOffset[0] = 0;
174  }
175  break;
176  case eAdamsBashforthOrder2:
177  {
178  m_numsteps = 2;
179  m_numstages = 1;
180 
181  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
182  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
183 
184  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
185  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
186  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
187  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
188 
189  m_B[0][0][0] = 3.0/2.0;
190  m_B[0][1][0] = 1.0;
191 
192  m_U[0][0] = 1.0;
193 
194  m_V[0][0] = 1.0;
195  m_V[0][1] = -0.5;
196 
197  m_schemeType = eExplicit;
198  m_numMultiStepValues = 1;
199  m_numMultiStepDerivs = 1;
200  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
201  m_timeLevelOffset[0] = 0;
202  m_timeLevelOffset[1] = 1;
203  }
204  break;
205  case eAdamsBashforthOrder3:
206  {
207  m_numsteps = 4;
208  m_numstages = 1;
209 
210  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
211  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
212 
213  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
214  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
215  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
216  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
217 
218  m_B[0][1][0] = 1.0;
219 
220  m_U[0][0] = 1.0;
221  m_U[0][1] = 23.0/12.0;
222  m_U[0][2] = -4.0/3.0;
223  m_U[0][3] = 5.0/12.0;
224 
225  m_V[0][0] = 1.0;
226  m_V[0][1] = 23.0/12.0;
227  m_V[0][2] = -4.0/3.0;
228  m_V[0][3] = 5.0/12.0;
229  m_V[2][1] = 1.0;
230  m_V[3][2] = 1.0;
231 
232  m_schemeType = eExplicit;
233  m_numMultiStepValues = 1;
234  m_numMultiStepDerivs = 3;
235  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
236  m_timeLevelOffset[0] = 0;
237  m_timeLevelOffset[1] = 1;
238  m_timeLevelOffset[2] = 2;
239  m_timeLevelOffset[3] = 3;
240  }
241  break;
242  case eAdamsBashforthOrder4:
243  {
244  m_numsteps = 5;
245  m_numstages = 1;
246 
247  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
248  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
249 
250  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
251  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
252  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,0.0);
253  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps,0.0);
254 
255  m_B[0][1][0] = 1.0;
256 
257  m_U[0][0] = 1.0;
258  m_U[0][1] = 55.0/24.0;
259  m_U[0][2] = -59.0/24.0;
260  m_U[0][3] = 37.0/24.0;
261  m_U[0][4] = -9.0/24.0;
262 
263  m_V[0][0] = 1.0;
264  m_V[0][1] = 55.0/24.0;
265  m_V[0][2] = -59.0/24.0;
266  m_V[0][3] = 37.0/24.0;
267  m_V[0][4] = -9.0/24.0;
268  m_V[2][1] = 1.0;
269  m_V[3][2] = 1.0;
270  m_V[4][3] = 1.0;
271 
272  m_schemeType = eExplicit;
273  m_numMultiStepValues = 1;
274  m_numMultiStepDerivs = 4;
275  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
276  m_timeLevelOffset[0] = 0;
277  m_timeLevelOffset[1] = 1;
278  m_timeLevelOffset[2] = 2;
279  m_timeLevelOffset[3] = 3;
280  m_timeLevelOffset[4] = 4;
281  }
282  break;
283  case eBackwardEuler:
284  case eBDFImplicitOrder1:
285  case eAdamsMoultonOrder1:
286  {
287  m_numsteps = 1;
288  m_numstages = 1;
289 
290  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
291  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
292 
293  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,1.0);
294  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
295  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
296  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 1.0);
297 
298  m_schemeType = eDiagonallyImplicit;
299  m_numMultiStepValues = 1;
300  m_numMultiStepDerivs = 0;
301  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
302  m_timeLevelOffset[0] = 0; // SJS: Not sure whether this is correct
303  }
304  break;
305  case eIMEXOrder1:
306  {
307  m_numsteps = 2;
308  m_numstages = 1;
309 
310  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
311  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
312 
313  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,1.0);
314  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
315  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
316  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,1.0);
317  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
318  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
319 
320  m_B[0][1][0] = 0.0;
321  m_B[1][0][0] = 0.0;
322  m_V[0][0] = 1.0;
323  m_V[0][1] = 1.0;
324 
325  m_schemeType = eIMEX;
326  m_numMultiStepValues = 1;
327  m_numMultiStepDerivs = 1;
328  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
329  m_timeLevelOffset[0] = 0;
330  m_timeLevelOffset[1] = 0;
331  }
332  break;
333  case eIMEXOrder2:
334  {
335  NekDouble third = 1.0/3.0;
336  m_numsteps = 4;
337  m_numstages = 1;
338 
339  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
340  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
341 
342  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,2*third);
343  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
344  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
345  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
346  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 4*third);
347  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
348 
349  m_B[0][0][0] = 2*third;
350  m_B[1][2][0] = 1.0;
351  m_U[0][1] = -third;
352  m_U[0][3] = -2*third;
353 
354  m_V[0][0] = 4*third;
355  m_V[0][1] = -third;
356  m_V[0][2] = 4*third;
357  m_V[0][3] = -2*third;
358  m_V[1][0] = 1.0;
359  m_V[3][2] = 1.0;
360 
361  m_schemeType = eIMEX;
362  m_numMultiStepValues = 2;
363  m_numMultiStepDerivs = 2;
364  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
365  m_timeLevelOffset[0] = 0;
366  m_timeLevelOffset[1] = 1;
367  m_timeLevelOffset[2] = 0;
368  m_timeLevelOffset[3] = 1;
369  }
370  break;
371  case eIMEXOrder3:
372  {
373  NekDouble eleventh = 1.0/11.0;
374  m_numsteps = 6;
375  m_numstages = 1;
376 
377  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
378  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
379 
380  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,6*eleventh);
381  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
382  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
383  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
384  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 18*eleventh);
385  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
386 
387  m_B[0][0][0] = 6*eleventh;
388  m_B[1][3][0] = 1.0;
389  m_U[0][1] = -9*eleventh;
390  m_U[0][2] = 2*eleventh;
391  m_U[0][4] = -18*eleventh;
392  m_U[0][5] = 6*eleventh;
393 
394  m_V[0][0] = 18*eleventh;
395  m_V[0][1] = -9*eleventh;
396  m_V[0][2] = 2*eleventh;
397  m_V[0][3] = 18*eleventh;
398  m_V[0][4] = -18*eleventh;
399  m_V[0][5] = 6*eleventh;
400  m_V[1][0] = 1.0;
401  m_V[2][1] = 1.0;
402  m_V[4][3] = 1.0;
403  m_V[5][4] = 1.0;
404 
405  m_schemeType = eIMEX;
406  m_numMultiStepValues = 3;
407  m_numMultiStepDerivs = 3;
408  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
409  m_timeLevelOffset[0] = 0;
410  m_timeLevelOffset[1] = 1;
411  m_timeLevelOffset[2] = 2;
412  m_timeLevelOffset[3] = 0;
413  m_timeLevelOffset[4] = 1;
414  m_timeLevelOffset[5] = 2;
415  }
416  break;
417  case eIMEXOrder4:
418  {
419  NekDouble twentyfifth = 1.0/25.0;
420  m_numsteps = 8;
421  m_numstages = 1;
422 
423  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
424  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
425 
426  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,12*twentyfifth);
427  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
428  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
429  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
430  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 48*twentyfifth);
431  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
432 
433  m_B[0][0][0] = 12*twentyfifth;
434  m_B[1][4][0] = 1.0;
435  m_U[0][1] = -36*twentyfifth;
436  m_U[0][2] = 16*twentyfifth;
437  m_U[0][3] = -3*twentyfifth;
438  m_U[0][5] = -72*twentyfifth;
439  m_U[0][7] = -12*twentyfifth;
440 
441  m_V[0][0] = 48*twentyfifth;
442  m_V[0][1] = -36*twentyfifth;
443  m_V[0][2] = 16*twentyfifth;
444  m_V[0][3] = -3*twentyfifth;
445  m_V[0][4] = 48*twentyfifth;
446  m_V[0][5] = -72*twentyfifth;
447  m_V[0][6] = 48*twentyfifth;
448  m_V[0][7] = -12*twentyfifth;
449  m_V[1][0] = 1.0;
450  m_V[2][1] = 1.0;
451  m_V[3][2] = 1.0;
452  m_V[5][4] = 1.0;
453  m_V[6][5] = 1.0;
454  m_V[7][6] = 1.0;
455 
456  m_schemeType = eIMEX;
457  m_numMultiStepValues = 4;
458  m_numMultiStepDerivs = 4;
459  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
460  m_timeLevelOffset[0] = 0;
461  m_timeLevelOffset[1] = 1;
462  m_timeLevelOffset[2] = 2;
463  m_timeLevelOffset[3] = 3;
464  m_timeLevelOffset[4] = 0;
465  m_timeLevelOffset[5] = 1;
466  m_timeLevelOffset[6] = 2;
467  m_timeLevelOffset[7] = 3;
468  }
469  break;
470  case eAdamsMoultonOrder2:
471  {
472  m_numsteps = 2;
473  m_numstages = 1;
474 
475  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
476  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
477 
478  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.5);
479  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
480  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
481  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.0);
482 
483  m_B[0][0][0] = 0.5;
484  m_B[0][1][0] = 1.0;
485 
486  m_U[0][0] = 1.0;
487  m_U[0][1] = 0.5;
488 
489  m_V[0][0] = 1.0;
490  m_V[0][1] = 0.5;
491 
492  m_schemeType = eDiagonallyImplicit;
493  m_numMultiStepValues = 1;
494  m_numMultiStepDerivs = 1;
495  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
496  m_timeLevelOffset[0] = 0;
497  m_timeLevelOffset[1] = 0;
498  }
499  break;
500  case eBDFImplicitOrder2:
501  {
502  NekDouble third = 1.0/3.0;
503  m_numsteps = 2;
504  m_numstages = 1;
505 
506  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
507  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
508 
509  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,2*third);
510  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
511  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
512  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.0);
513 
514  m_B[0][0][0] = 2*third;
515  m_B[0][1][0] = 0.0;
516 
517  m_U[0][0] = 4*third;
518  m_U[0][1] = -third;
519 
520  m_V[0][0] = 4*third;
521  m_V[0][1] = -third;
522  m_V[1][0] = 1.0;
523 
524  m_schemeType = eDiagonallyImplicit;
525  m_numMultiStepValues = 2;
526  m_numMultiStepDerivs = 0;
527  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
528  m_timeLevelOffset[0] = 0;
529  m_timeLevelOffset[1] = 1;
530  }
531  break;
532  case eMidpoint:
533  case eRungeKutta2:
534  {
535  m_numsteps = 1;
536  m_numstages = 2;
537 
538  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
539  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
540 
541  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
542  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
543  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
544  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
545 
546  m_A[0][1][0] = 0.5;
547  m_B[0][0][1] = 1.0;
548 
549  m_schemeType = eExplicit;
550  m_numMultiStepValues = 1;
551  m_numMultiStepDerivs = 0;
552  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
553  m_timeLevelOffset[0] = 0;
554  }
555  break;
556  case eRungeKutta2_ImprovedEuler:
557  case eRungeKutta2_SSP:
558  {
559  m_numsteps = 1;
560  m_numstages = 2;
561 
562  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
563  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
564 
565  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
566  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
567  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
568  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
569 
570  m_A[0][1][0] = 1.0;
571 
572  m_B[0][0][0] = 0.5;
573  m_B[0][0][1] = 0.5;
574 
575  m_schemeType = eExplicit;
576  m_numMultiStepValues = 1;
577  m_numMultiStepDerivs = 0;
578  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
579  m_timeLevelOffset[0] = 0;
580  }
581  break;
582  case eRungeKutta3_SSP:
583  {
584  m_numsteps = 1;
585  m_numstages = 3;
586 
587  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
588  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
589 
590  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
591  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
592  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
593  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
594 
595  m_A[0][1][0] = 1.0;
596  m_A[0][2][0] = 0.25;
597  m_A[0][2][1] = 0.25;
598 
599  m_B[0][0][0] = 1.0/6.0;
600  m_B[0][0][1] = 1.0/6.0;
601  m_B[0][0][2] = 2.0/3.0;
602 
603  m_schemeType = eExplicit;
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 eClassicalRungeKutta4:
611  case eRungeKutta4:
612  {
613  m_numsteps = 1;
614  m_numstages = 4;
615 
616  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
617  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
618 
619  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
620  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
621  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
622  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
623 
624  m_A[0][1][0] = 0.5;
625  m_A[0][2][1] = 0.5;
626  m_A[0][3][2] = 1.0;
627 
628  m_B[0][0][0] = 1.0/6.0;
629  m_B[0][0][1] = 1.0/3.0;
630  m_B[0][0][2] = 1.0/3.0;
631  m_B[0][0][3] = 1.0/6.0;
632 
633  m_schemeType = eExplicit;
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 eRungeKutta5:
641  {
642  m_numsteps = 1;
643  m_numstages = 6;
644 
645  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
646  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
647 
648  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
649  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
650  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
651  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
652 
653  m_A[0][1][0] = 1.0/4.0;
654  m_A[0][2][0] = 1.0/8.0;
655  m_A[0][2][1] = 1.0/8.0;
656  m_A[0][3][2] = 1.0/2.0;
657  m_A[0][4][0] = 3.0/16.0;
658  m_A[0][4][1] = -3.0/8.0;
659  m_A[0][4][2] = 3.0/8.0;
660  m_A[0][4][3] = 9.0/16.0;
661  m_A[0][5][0] = -3.0/7.0;
662  m_A[0][5][1] = 8.0/7.0;
663  m_A[0][5][2] = 6.0/7.0;
664  m_A[0][5][3] = -12.0/7.0;
665  m_A[0][5][4] = 8.0/7.0;
666 
667  m_B[0][0][0] = 7.0/90.0;
668  m_B[0][0][2] = 32.0/90.0;
669  m_B[0][0][3] = 12.0/90.0;
670  m_B[0][0][4] = 32.0/90.0;
671  m_B[0][0][5] = 7.0/90.0;
672 
673  m_schemeType = eExplicit;
674  m_numMultiStepValues = 1;
675  m_numMultiStepDerivs = 0;
676  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
677  m_timeLevelOffset[0] = 0;
678  }
679  break;
680  case eDIRKOrder2:
681  {
682  m_numsteps = 1;
683  m_numstages = 2;
684 
685  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
686  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
687 
688  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
689  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
690  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
691  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
692 
693  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
694 
695  m_A[0][0][0] = lambda;
696  m_A[0][1][0] = 1.0 - lambda;
697  m_A[0][1][1] = lambda;
698 
699  m_B[0][0][0] = 1.0 - lambda;
700  m_B[0][0][1] = lambda;
701 
702  m_schemeType = eDiagonallyImplicit;
703  m_numMultiStepValues = 1;
704  m_numMultiStepDerivs = 0;
705  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
706  m_timeLevelOffset[0] = 0;
707  }
708  break;
709  case eDIRKOrder3:
710  {
711  m_numsteps = 1;
712  m_numstages = 3;
713 
714  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
715  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
716 
717  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
718  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
719  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
720  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
721 
722  NekDouble lambda = 0.4358665215084589;
723 
724  m_A[0][0][0] = lambda;
725  m_A[0][1][0] = 0.5 * (1.0 - lambda);
726  m_A[0][2][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
727  m_A[0][1][1] = lambda;
728  m_A[0][2][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
729  m_A[0][2][2] = lambda;
730 
731  m_B[0][0][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
732  m_B[0][0][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
733  m_B[0][0][2] = lambda;
734 
735  m_schemeType = eDiagonallyImplicit;
736  m_numMultiStepValues = 1;
737  m_numMultiStepDerivs = 0;
738  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
739  m_timeLevelOffset[0] = 0;
740  }
741  break;
742  case eIMEXdirk_2_3_2:
743  {
744  m_numsteps = 1;
745  m_numstages = 3;
746 
747  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
748  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
749 
750  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
751  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
752  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
753  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
754  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
755  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
756 
757  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
758  NekDouble delta = -2.0*sqrt(2.0)/3.0;
759 
760  m_A[0][1][1] = lambda;
761  m_A[0][2][1] = 1.0 - lambda;
762  m_A[0][2][2] = lambda;
763 
764  m_B[0][0][1] = 1.0 - lambda;
765  m_B[0][0][2] = lambda;
766 
767  m_A[1][1][0] = lambda;
768  m_A[1][2][0] = delta;
769  m_A[1][2][1] = 1.0 - delta;
770 
771  m_B[1][0][1] = 1.0 - lambda;
772  m_B[1][0][2] = lambda;
773 
774  m_schemeType = eIMEX;
775  m_numMultiStepValues = 1;
776  m_numMultiStepDerivs = 0;
777  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
778  m_timeLevelOffset[0] = 0;
779  }
780  break;
781  case eIMEXdirk_3_4_3:
782  {
783  m_numsteps = 1;
784  m_numstages = 4;
785 
786  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
787  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
788 
789  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
790  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
791  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
792  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
793  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
794  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
795 
796  NekDouble lambda = 0.4358665215;
797 
798  m_A[0][1][1] = lambda;
799  m_A[0][2][1] = 0.5 * (1.0 - lambda);
800  m_A[0][3][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
801  m_A[0][2][2] = lambda;
802  m_A[0][3][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
803  m_A[0][3][3] = lambda;
804 
805  m_B[0][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
806  m_B[0][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
807  m_B[0][0][3] = lambda;
808 
809  m_A[1][1][0] = 0.4358665215;
810  m_A[1][2][0] = 0.3212788860;
811  m_A[1][2][1] = 0.3966543747;
812  m_A[1][3][0] =-0.105858296;
813  m_A[1][3][1] = 0.5529291479;
814  m_A[1][3][2] = 0.5529291479;
815 
816  m_B[1][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
817  m_B[1][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
818  m_B[1][0][3] = lambda;
819 
820  m_schemeType = eIMEX;
821  m_numMultiStepValues = 1;
822  m_numMultiStepDerivs = 0;
823  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
824  m_timeLevelOffset[0] = 0;
825  }
826  break;
827  case eCNAB:
828  {
829  NekDouble secondth = 1.0/2.0;
830  m_numsteps = 4;
831  m_numstages = 1;
832 
833  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
834  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
835 
836  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,secondth);
837  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
838  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
839  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
840  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,secondth);
841  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
842 
843  m_B[0][0][0] = secondth;
844  m_B[0][1][0] = 1.0;
845  m_B[1][2][0] = 1.0;
846  m_U[0][0] = 2*secondth;
847  m_U[0][2] = 3*secondth;
848  m_U[0][3] = -1*secondth;
849 
850  m_V[0][0] = 2*secondth;
851  m_V[0][1] = secondth;
852  m_V[0][2] = 3*secondth;
853  m_V[0][3] = -1*secondth;
854  m_V[3][2] = 1.0;
855 
856  m_schemeType = eIMEX;
857  m_numMultiStepValues = 1;
858  m_numMultiStepDerivs = 3;
859  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
860  m_timeLevelOffset[0] = 0;
861  m_timeLevelOffset[1] = 0;
862  m_timeLevelOffset[2] = 0;
863  m_timeLevelOffset[3] = 1;
864  }
865  break;
866  case eIMEXGear:
867  {
868  NekDouble twothirdth = 2.0/3.0;
869  m_numsteps = 3;
870  m_numstages = 1;
871 
872  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
873  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
874 
875  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
876  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
877  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
878  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
879  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
880  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
881 
882  m_B[0][0][0] = twothirdth;
883  m_B[1][2][0] = 1.0;
884  m_U[0][0] = 2*twothirdth;
885  m_U[0][1] = -0.5*twothirdth;
886  m_U[0][2] = twothirdth;
887 
888  m_V[0][0] = 2*twothirdth;
889  m_V[0][1] = -0.5*twothirdth;
890  m_V[0][2] = twothirdth;
891  m_V[1][0] = 1.0;
892 
893  m_schemeType = eIMEX;
894  m_numMultiStepValues = 2;
895  m_numMultiStepDerivs = 1;
896  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
897  m_timeLevelOffset[0] = 0;
898  m_timeLevelOffset[1] = 1;
899  m_timeLevelOffset[2] = 0;
900  }
901  break;
902  case eMCNAB:
903  {
904  NekDouble sixthx = 9.0/16.0;
905  m_numsteps = 5;
906  m_numstages = 1;
907 
908  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
909  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
910 
911  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,sixthx);
912  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
913  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
914  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
915  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
916  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
917 
918  m_B[0][0][0] = sixthx;
919  m_B[0][1][0] = 1.0;
920  m_B[1][3][0] = 1.0;
921  m_U[0][0] = 1.0;
922  m_U[0][1] = 6.0/16.0;
923  m_U[0][2] = 1.0/16.0;
924  m_U[0][3] = 1.5;
925  m_U[0][4] = -0.5;
926 
927  m_V[0][0] = 1.0;
928  m_V[0][1] = 6.0/16.0;
929  m_V[0][2] = 1.0/16.0;
930  m_V[0][3] = 1.5;
931  m_V[0][4] = -0.5;
932  m_V[2][1] = 1.0;
933  m_V[4][3] = 1.0;
934 
935  m_schemeType = eIMEX;
936  m_numMultiStepValues = 1;
937  m_numMultiStepDerivs = 4;
938  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
939  m_timeLevelOffset[0] = 0;
940  m_timeLevelOffset[1] = 0;
941  m_timeLevelOffset[2] = 1;
942  m_timeLevelOffset[3] = 0;
943  m_timeLevelOffset[4] = 1;
944  }
945  break;
946  case eIMEXdirk_2_2_2:
947  {
948  m_numsteps = 1;
949  m_numstages = 3;
950 
951  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
952  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
953 
954  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
955  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
956  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
957  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
958  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
959  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
960 
961  NekDouble glambda = 0.2928932188134524756;
962  NekDouble gdelta = -0.7071067811865475244;
963 
964  m_A[0][1][1] = glambda;
965  m_A[0][2][1] = 1.0 - glambda;
966  m_A[0][2][2] = glambda;
967 
968  m_B[0][0][1] = 1.0 - glambda;
969  m_B[0][0][2] = glambda;
970 
971  m_A[1][1][0] = glambda;
972  m_A[1][2][0] = gdelta;
973  m_A[1][2][1] = 1.0 - gdelta;
974 
975  m_B[1][0][0] = gdelta;
976  m_B[1][0][1] = 1.0 - gdelta;
977 
978  m_schemeType = eIMEX;
979  m_numMultiStepValues = 1;
980  m_numMultiStepDerivs = 0;
981  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
982  m_timeLevelOffset[0] = 0;
983  }
984  break;
985  case eIMEXdirk_2_3_3:
986  {
987  m_numsteps = 1;
988  m_numstages = 3;
989 
990  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
991  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
992 
993  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
994  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
995  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
996  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
997  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
998  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
999 
1000  NekDouble glambda = 0.78867513459481288226;
1001 
1002  m_A[0][1][1] = glambda;
1003  m_A[0][2][1] = 1.0 - 2.0*glambda;
1004  m_A[0][2][2] = glambda;
1005 
1006  m_B[0][0][1] = 0.5;
1007  m_B[0][0][2] = 0.5;
1008 
1009  m_A[1][1][0] = glambda;
1010  m_A[1][2][0] = glambda - 1.0;
1011  m_A[1][2][1] = 2.0*(1-glambda);
1012 
1013  m_B[1][0][1] = 0.5;
1014  m_B[1][0][2] = 0.5;
1015 
1016  m_schemeType = eIMEX;
1017  m_numMultiStepValues = 1;
1018  m_numMultiStepDerivs = 0;
1019  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1020  m_timeLevelOffset[0] = 0;
1021  }
1022  break;
1023  case eIMEXdirk_1_1_1:
1024  {
1025  m_numsteps = 1;
1026  m_numstages = 2;
1027 
1028  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
1029  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
1030 
1031  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1032  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1033  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1034  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1035  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
1036  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
1037 
1038  m_A[0][1][1] = 1.0;
1039 
1040  m_B[0][0][1] = 1.0;
1041 
1042  m_A[1][1][0] = 1.0;
1043 
1044  m_B[1][0][0] = 1.0;
1045 
1046  m_schemeType = eIMEX;
1047  m_numMultiStepValues = 1;
1048  m_numMultiStepDerivs = 0;
1049  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1050  m_timeLevelOffset[0] = 0;
1051  }
1052  break;
1053  case eIMEXdirk_1_2_1:
1054  {
1055  m_numsteps = 1;
1056  m_numstages = 2;
1057 
1058  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
1059  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
1060 
1061  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1062  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1063  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1064  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1065  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
1066  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
1067 
1068  m_A[0][1][1] = 1.0;
1069 
1070  m_B[0][0][1] = 1.0;
1071 
1072  m_A[1][1][0] = 1.0;
1073 
1074  m_B[1][0][1] = 1.0;
1075 
1076  m_schemeType = eIMEX;
1077  m_numMultiStepValues = 1;
1078  m_numMultiStepDerivs = 0;
1079  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1080  m_timeLevelOffset[0] = 0;
1081  }
1082  break;
1083  case eIMEXdirk_1_2_2:
1084  {
1085  m_numsteps = 1;
1086  m_numstages = 2;
1087 
1088  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
1089  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
1090 
1091  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1092  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1093  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1094  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1095  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
1096  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
1097 
1098  m_A[0][1][1] = 1.0/2.0;
1099 
1100  m_B[0][0][1] = 1.0;
1101 
1102  m_A[1][1][0] = 1.0/2.0;
1103 
1104  m_B[1][0][1] = 1.0;
1105 
1106  m_schemeType = eIMEX;
1107  m_numMultiStepValues = 1;
1108  m_numMultiStepDerivs = 0;
1109  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1110  m_timeLevelOffset[0] = 0;
1111  }
1112  break;
1113  case eIMEXdirk_4_4_3:
1114  {
1115  m_numsteps = 1;
1116  m_numstages = 5;
1117 
1118  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
1119  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
1120 
1121  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1122  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1123  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
1124  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
1125  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
1126  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
1127 
1128  m_A[0][1][1] = 1.0/2.0;
1129  m_A[0][2][1] = 1.0/6.0;
1130  m_A[0][2][2] = 1.0/2.0;
1131  m_A[0][3][1] = -1.0/2.0;
1132  m_A[0][3][2] = 1.0/2.0;
1133  m_A[0][3][3] = 1.0/2.0;
1134  m_A[0][4][1] = 3.0/2.0;
1135  m_A[0][4][2] = -3.0/2.0;
1136  m_A[0][4][3] = 1.0/2.0;
1137  m_A[0][4][4] = 1.0/2.0;
1138 
1139  m_B[0][0][1] = 3.0/2.0;
1140  m_B[0][0][2] = -3.0/2.0;
1141  m_B[0][0][3] = 1.0/2.0;
1142  m_B[0][0][4] = 1.0/2.0;
1143 
1144  m_A[1][1][0] = 1.0/2.0;
1145  m_A[1][2][0] = 11.0/18.0;
1146  m_A[1][2][1] = 1.0/18.0;
1147  m_A[1][3][0] = 5.0/6.0;
1148  m_A[1][3][1] = -5.0/6.0;
1149  m_A[1][3][2] = 1.0/2.0;
1150  m_A[1][4][0] = 1.0/4.0;
1151  m_A[1][4][1] = 7.0/4.0;
1152  m_A[1][4][2] = 3.0/4.0;
1153  m_A[1][4][3] = -7.0/4.0;
1154 
1155  m_B[1][0][0] = 1.0/4.0;
1156  m_B[1][0][1] = 7.0/4.0;
1157  m_B[1][0][2] = 3.0/4.0;
1158  m_B[1][0][3] = -7.0/4.0;
1159 
1160  m_schemeType = eIMEX;
1161  m_numMultiStepValues = 1;
1162  m_numMultiStepDerivs = 0;
1163  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1164  m_timeLevelOffset[0] = 0;
1165  }
1166  break;
1167  default:
1168  {
1169  NEKERROR(ErrorUtil::efatal,"Invalid Time Integration Scheme");
1170  }
1171  }
1172 
1173  m_firstStageEqualsOldSolution = CheckIfFirstStageEqualsOldSolution(m_A,m_B,m_U,m_V);
1174  m_lastStageEqualsNewSolution = CheckIfLastStageEqualsNewSolution(m_A,m_B,m_U,m_V);
1175 
1176  ASSERTL1(VerifyIntegrationSchemeType(m_schemeType,m_A,m_B,m_U,m_V),
1177  "Time integration scheme coefficients do not match its type");
1178  }
1179 
1180 
1181  bool TimeIntegrationScheme::
1182  VerifyIntegrationSchemeType(TimeIntegrationSchemeType type,
1183  const Array<OneD, const Array<TwoD, NekDouble> >& A,
1184  const Array<OneD, const Array<TwoD, NekDouble> >& B,
1185  const Array<TwoD, const NekDouble>& U,
1186  const Array<TwoD, const NekDouble>& V) const
1187  {
1188  boost::ignore_unused(B, U, V);
1189 
1190  int i;
1191  int j;
1192  int m;
1193  int IMEXdim = A.num_elements();
1194  int dim = A[0].GetRows();
1195 
1196  Array<OneD, TimeIntegrationSchemeType> vertype(IMEXdim,eExplicit);
1197 
1198  for(m = 0; m < IMEXdim; m++)
1199  {
1200  for(i = 0; i < dim; i++)
1201  {
1202  if( fabs(A[m][i][i]) > NekConstants::kNekZeroTol )
1203  {
1204  vertype[m] = eDiagonallyImplicit;
1205  }
1206  }
1207 
1208  for(i = 0; i < dim; i++)
1209  {
1210  for(j = i+1; j < dim; j++)
1211  {
1212  if( fabs(A[m][i][j]) > NekConstants::kNekZeroTol )
1213  {
1214  vertype[m] = eImplicit;
1215  ASSERTL1(false,"Fully Implicit schemes cannnot be handled by the TimeIntegrationScheme class");
1216  }
1217  }
1218  }
1219  }
1220 
1221  if(IMEXdim == 2)
1222  {
1223  ASSERTL1(B.num_elements()==2,"Coefficient Matrix B should have an implicit and explicit part for IMEX schemes");
1224  if((vertype[0] == eDiagonallyImplicit) &&
1225  (vertype[1] == eExplicit))
1226  {
1227  vertype[0] = eIMEX;
1228  }
1229  else
1230  {
1231  ASSERTL1(false,"This is not a proper IMEX scheme");
1232  }
1233  }
1234 
1235  return (vertype[0] == type);
1236  }
1237 
1238  TimeIntegrationSolutionSharedPtr
1239  TimeIntegrationScheme::InitializeScheme(const NekDouble timestep,
1240  ConstDoubleArray &y_0 ,
1241  const NekDouble time ,
1242  const TimeIntegrationSchemeOperators &op)
1243  {
1244  // create a TimeIntegrationSolution object based upon the
1245  // initial value. Initialise all other multi-step values
1246  // and derivatives to zero
1247  TimeIntegrationSolutionSharedPtr y_out =
1248  MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(m_schemeKey,y_0,time,timestep);
1249 
1250  if( GetIntegrationSchemeType() == eExplicit)
1251  {
1252  // ensure initial solution is in correct space
1253  op.DoProjection(y_0,y_out->UpdateSolution(),time);
1254  }
1255 
1256  // calculate the initial derivative, if is part of the
1257  // solution vector of the current scheme
1258  if(m_numMultiStepDerivs)
1259  {
1260  if(m_timeLevelOffset[m_numMultiStepValues] == 0)
1261  {
1262  int i;
1263  int nvar = y_0.num_elements();
1264  int npoints = y_0[0].num_elements();
1265  DoubleArray f_y_0(nvar);
1266  for(i = 0; i < nvar; i++)
1267  {
1268  f_y_0[i] = Array<OneD,NekDouble>(npoints);
1269  }
1270  // calculate the derivative of the initial value
1271  op.DoOdeRhs(y_0,f_y_0,time);
1272 
1273  // multiply by the step size
1274  for(i = 0; i < nvar; i++)
1275  {
1276  Blas::Dscal(npoints,timestep,f_y_0[i].get(),1);
1277  }
1278  y_out->SetDerivative(0,f_y_0,timestep);
1279  }
1280  }
1281 
1282  return y_out;
1283  }
1284 
1285  TimeIntegrationScheme::ConstDoubleArray&
1286  TimeIntegrationScheme::TimeIntegrate(const NekDouble timestep,
1287  TimeIntegrationSolutionSharedPtr &solvector,
1288  const TimeIntegrationSchemeOperators &op)
1289  {
1290  ASSERTL1(!(GetIntegrationSchemeType() == eImplicit),
1291  "Fully Implicit integration scheme cannot be handled by this routine.");
1292 
1293  int nvar = solvector->GetFirstDim ();
1294  int npoints = solvector->GetSecondDim();
1295 
1296  if( (solvector->GetIntegrationScheme()).get() != this )
1297  {
1298  // This branch will be taken when the solution vector
1299  // (solvector) is set up for a different scheme than
1300  // the object this method is called from. (typically
1301  // needed to calculate the first time-levels of a
1302  // multi-step scheme)
1303 
1304  // To do this kind of 'non-matching' integration, we
1305  // perform the following three steps:
1306  //
1307  // 1: copy the required input information from the
1308  // solution vector of the master scheme to the
1309  // input solution vector of the current scheme
1310  //
1311  // 2: time-integrate for one step using the current
1312  // scheme
1313  //
1314  // 3: copy the information contained in the output
1315  // vector of the current scheme to the solution
1316  // vector of the master scheme
1317 
1318  // STEP 1: copy the required input information from
1319  // the solution vector of the master scheme
1320  // to the input solution vector of the
1321  // current scheme
1322 
1323  // 1.1 Determine which information is required for the
1324  // current scheme
1325  int n;
1326  DoubleArray y_n;
1327  NekDouble t_n = 0;
1328  DoubleArray dtFy_n;
1329  unsigned int nCurSchemeVals = m_numMultiStepValues; // number of required values of the current scheme
1330  unsigned int nCurSchemeDers = m_numMultiStepDerivs; // number of required derivs of the current scheme
1331  unsigned int nCurSchemeSteps = m_numsteps; // number of steps in the current scheme
1332  unsigned int nMasterSchemeVals = solvector->GetNvalues(); // number of values of the master scheme
1333  unsigned int nMasterSchemeDers = solvector->GetNderivs(); // number of derivs of the master scheme
1334  // The arrays below contains information to which
1335  // time-level the values and derivatives of the
1336  // schemes belong
1337  const Array<OneD, const unsigned int>& curTimeLevels = m_timeLevelOffset;
1338  const Array<OneD, const unsigned int>& masterTimeLevels = solvector->GetTimeLevelOffset();
1339 
1340  // 1.2 Copy the required information from the master
1341  // solution vector to the input solution vector of
1342  // the current scheme
1343  TimeIntegrationSolutionSharedPtr solvector_in = MemoryManager<TimeIntegrationSolution>::
1344  AllocateSharedPtr(GetIntegrationSchemeKey()); // input solution vector of the current scheme
1345 
1346  for(n = 0; n < nCurSchemeVals; n++)
1347  {
1348  // Get the required value out of the master solution vector
1349  //DoubleArray& y_n = solvector->GetValue ( curTimeLevels[n] );
1350  //NekDouble t_n = solvector->GetValueTime( curTimeLevels[n] );
1351 
1352  y_n = solvector->GetValue ( curTimeLevels[n] );
1353  t_n = solvector->GetValueTime( curTimeLevels[n] );
1354 
1355  // Set the required value in the input solution
1356  // vector of the current scheme
1357  solvector_in->SetValue(curTimeLevels[n],y_n,t_n);
1358  }
1359  for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
1360  {
1361  // Get the required derivative out of the master
1362  // solution vector
1363  //DoubleArray& dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
1364  dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
1365 
1366  // Set the required derivative in the input
1367  // solution vector of the current scheme
1368  solvector_in->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
1369  }
1370 
1371  // STEP 2: time-integrate for one step using the
1372  // current scheme
1373  TimeIntegrationSolutionSharedPtr solvector_out = MemoryManager<TimeIntegrationSolution>:: AllocateSharedPtr(GetIntegrationSchemeKey(),nvar,npoints); // output solution vector of the current scheme
1374 
1375  // integrate
1376  TimeIntegrate(timestep, solvector_in->GetSolutionVector(),
1377  solvector_in->GetTimeVector(),
1378  solvector_out->UpdateSolutionVector(),
1379  solvector_out->UpdateTimeVector(),op);
1380 
1381 
1382  // STEP 3: copy the information contained in the
1383  // output vector of the current scheme to the
1384  // solution vector of the master scheme
1385 
1386  // 3.1 Check whether the current time scheme updates
1387  // the most recent derivative that should be
1388  // updated in the master scheme. If not,
1389  // calculate the derivative. This can be done
1390  // based upon the corresponding value and the
1391  // DoOdeRhs operator.
1392  int j;
1393  bool CalcNewDeriv = false; // flag inidicating whether the new derivative is availble in the output of
1394  // of the current scheme or whether it should be calculated
1395  if( nMasterSchemeDers > 0 )
1396  {
1397  if(nCurSchemeDers == 0)
1398  {
1399  CalcNewDeriv = true;
1400  }
1401  else
1402  {
1403  if( masterTimeLevels[nMasterSchemeVals] < curTimeLevels[nCurSchemeVals] )
1404  {
1405  CalcNewDeriv = true;
1406  }
1407  }
1408  }
1409 
1410  if(CalcNewDeriv)
1411  {
1412  int newDerivTimeLevel = masterTimeLevels[nMasterSchemeVals]; // contains the time level at which
1413  // we want to know the derivative of the
1414  // master scheme
1415  //DoubleArray y_n;
1416  //NekDouble t_n;
1417  // if the time level correspond to 0, calculate the derivative based upon the solution value
1418  // at the new time-level
1419  if (newDerivTimeLevel == 0)
1420  {
1421  y_n = solvector_out->GetValue(0);
1422  t_n = solvector_out->GetValueTime(0);
1423  }
1424  // if the time level correspond to 1, calculate the derivative based upon the solution value
1425  // at the new old-level
1426  else if( newDerivTimeLevel == 1 )
1427  {
1428  y_n = solvector->GetValue(0);
1429  t_n = solvector->GetValueTime(0);
1430  }
1431  else
1432  {
1433  ASSERTL1(false,"Problems with initialising scheme");
1434  }
1435 
1436  DoubleArray f_n(nvar);
1437  for(j = 0; j < nvar; j++)
1438  {
1439  f_n[j] = Array<OneD, NekDouble>(npoints);
1440  }
1441 
1442  // calculate the derivative
1443  op.DoOdeRhs(y_n, f_n, t_n);
1444 
1445  // multiply by dt (as required by the General Linear Method framework)
1446  for(j = 0; j < nvar; j++)
1447  {
1448  Vmath::Smul(npoints,timestep,f_n[j],1,
1449  f_n[j],1);
1450  }
1451 
1452  // Rotate the solution vector
1453  // (i.e. updating without calculating/inserting new values)
1454  solvector->RotateSolutionVector();
1455  // Set the calculated derivative in the master solution vector
1456  solvector->SetDerivative(newDerivTimeLevel,f_n,timestep);
1457  }
1458  else
1459  {
1460  // Rotate the solution vector (i.e. updating
1461  // without calculating/inserting new values)
1462  solvector->RotateSolutionVector();
1463  }
1464 
1465 
1466  // 1.2 Copy the information calculated using the
1467  // current scheme from the output solution vector
1468  // to the master solution vector
1469  for(n = 0; n < nCurSchemeVals; n++)
1470  {
1471  // Get the calculated value out of the output
1472  // solution vector of the current scheme
1473  //DoubleArray& y_n = solvector_out->GetValue ( curTimeLevels[n] );
1474  //NekDouble t_n = solvector_out->GetValueTime( curTimeLevels[n] );
1475  y_n = solvector_out->GetValue ( curTimeLevels[n] );
1476  t_n = solvector_out->GetValueTime( curTimeLevels[n] );
1477 
1478  // Set the calculated value in the master solution vector
1479  solvector->SetValue(curTimeLevels[n],y_n,t_n);
1480  }
1481 
1482  for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
1483  {
1484  // Get the calculated derivative out of the output
1485  // solution vector of the current scheme
1486  // DoubleArray& dtFy_n =
1487  // solvector_out->GetDerivative (curTimeLevels[n]);
1488  dtFy_n = solvector_out->GetDerivative ( curTimeLevels[n] );
1489 
1490  // Set the calculated derivative in the master
1491  // solution vector
1492  solvector->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
1493  }
1494  }
1495  else
1496  {
1497  const TimeIntegrationSchemeKey& key = solvector->GetIntegrationSchemeKey();
1498 
1499  TimeIntegrationSolutionSharedPtr solvector_new = MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(key,nvar,npoints);
1500 
1501  TimeIntegrate(timestep,solvector->GetSolutionVector(),
1502  solvector->GetTimeVector(),
1503  solvector_new->UpdateSolutionVector(),
1504  solvector_new->UpdateTimeVector(),op);
1505 
1506  solvector = solvector_new;
1507  }
1508  return solvector->GetSolution();
1509  }
1510 
1511  void TimeIntegrationScheme::TimeIntegrate(const NekDouble timestep,
1512  ConstTripleArray &y_old ,
1513  ConstSingleArray &t_old ,
1514  TripleArray &y_new ,
1515  SingleArray &t_new ,
1516  const TimeIntegrationSchemeOperators &op)
1517  {
1518  ASSERTL1(CheckTimeIntegrateArguments(timestep,y_old,t_old,y_new,t_new,op), "Arguments not well defined");
1519 
1520  unsigned int i,j,k;
1521  TimeIntegrationSchemeType type = GetIntegrationSchemeType();
1522 
1523  // Check if storage has already been initialised.
1524  // If so, we just zero the temporary storage.
1525  if (m_initialised && m_nvar == GetFirstDim(y_old)
1526  && m_npoints == GetSecondDim(y_old))
1527  {
1528  for(j = 0; j < m_nvar; j++)
1529  {
1530  Vmath::Zero(m_npoints, m_tmp[j], 1);
1531  }
1532  }
1533  else
1534  {
1535  m_nvar = GetFirstDim(y_old);
1536  m_npoints = GetSecondDim(y_old);
1537 
1538  // First, we are going to calculate the various stage
1539  // values and stage derivatives (this is the multi-stage
1540  // part of the method)
1541  // - m_Y corresponds to the stage values
1542  // - m_F corresponds to the stage derivatives
1543  // - m_T corresponds to the time at the different stages
1544  // - m_tmp corresponds to the explicit right hand side of
1545  // each stage equation
1546  // (for explicit schemes, this correspond to m_Y)
1547 
1548  // Allocate memory for the arrays m_Y and m_F and m_tmp The same
1549  // storage will be used for every stage -> m_Y is a
1550  // DoubleArray
1551  m_tmp = DoubleArray(m_nvar);
1552  for(j = 0; j < m_nvar; j++)
1553  {
1554  m_tmp[j] = Array<OneD, NekDouble>(m_npoints,0.0);
1555  }
1556 
1557  // The same storage will be used for every stage -> m_tmp is
1558  // a DoubleArray
1559  if(type == eExplicit)
1560  {
1561  m_Y = m_tmp;
1562  }
1563  else
1564  {
1565  m_Y = DoubleArray(m_nvar);
1566  for(j = 0; j < m_nvar; j++)
1567  {
1568  m_Y[j] = Array<OneD, NekDouble>(m_npoints,0.0);
1569  }
1570  }
1571 
1572  // Different storage for every stage derivative as the data
1573  // will be re-used to update the solution -> m_F is a TripleArray
1574  m_F = TripleArray(m_numstages);
1575  for(i = 0; i < m_numstages; ++i)
1576  {
1577  m_F[i] = DoubleArray(m_nvar);
1578  for(j = 0; j < m_nvar; j++)
1579  {
1580  m_F[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
1581  }
1582  }
1583 
1584  if(type == eIMEX)
1585  {
1586  m_F_IMEX = TripleArray(m_numstages);
1587  for(i = 0; i < m_numstages; ++i)
1588  {
1589  m_F_IMEX[i] = DoubleArray(m_nvar);
1590  for(j = 0; j < m_nvar; j++)
1591  {
1592  m_F_IMEX[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
1593  }
1594  }
1595  }
1596 
1597  // Finally, flag that we have initialised the memory.
1598  m_initialised = true;
1599  }
1600 
1601  // The loop below calculates the stage values and derivatives
1602  for(i = 0; i < m_numstages; i++)
1603  {
1604  if( (i==0) && m_firstStageEqualsOldSolution )
1605  {
1606  for(k = 0; k < m_nvar; k++)
1607  {
1608  Vmath::Vcopy(m_npoints,y_old[0][k],1,m_Y[k],1);
1609  }
1610  m_T = t_old[0];
1611  }
1612  else
1613  {
1614  // The stage values m_Y are a linear combination of:
1615  // 1: the stage derivatives
1616 
1617  if( i != 0 )
1618  {
1619  for(k = 0; k < m_nvar; k++)
1620  {
1621  Vmath::Smul(m_npoints,timestep*A(i,0),m_F[0][k],1,
1622  m_tmp[k],1);
1623 
1624  if(type == eIMEX)
1625  {
1626  Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,0),
1627  m_F_IMEX[0][k],1,
1628  m_tmp[k],1,m_tmp[k],1);
1629  }
1630  }
1631  }
1632  m_T = A(i,0)*timestep;
1633 
1634  for( j = 1; j < i; j++ )
1635  {
1636  for(k = 0; k < m_nvar; k++)
1637  {
1638  Vmath::Svtvp(m_npoints,timestep*A(i,j),m_F[j][k],1,
1639  m_tmp[k],1,m_tmp[k],1);
1640  if(type == eIMEX)
1641  {
1642  Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,j),
1643  m_F_IMEX[j][k],1,
1644  m_tmp[k],1,m_tmp[k],1);
1645  }
1646  }
1647 
1648  m_T += A(i,j)*timestep;
1649  }
1650 
1651  // 2: the imported multi-step solution of the
1652  // previous time level
1653  for(j = 0; j < m_numsteps; j++)
1654  {
1655  for(k = 0; k < m_nvar; k++)
1656  {
1657  Vmath::Svtvp(m_npoints,U(i,j),y_old[j][k],1,
1658  m_tmp[k],1,m_tmp[k],1);
1659  }
1660  m_T += U(i,j)*t_old[j];
1661  }
1662  }
1663 
1664  // Calculate the stage derivative based upon the stage value
1665  if(type == eDiagonallyImplicit)
1666  {
1667  if(m_numstages==1)
1668  {
1669  m_T= t_old[0]+timestep;
1670  }
1671  else
1672  {
1673  m_T= t_old[0];
1674  for(int j=0; j<=i; ++j)
1675  {
1676  m_T += A(i,j)*timestep;
1677  }
1678  }
1679 
1680  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
1681 
1682  for(k = 0; k < m_nvar; k++)
1683  {
1684  Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
1685  Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),m_F[i][k],1,m_F[i][k],1);
1686  }
1687  }
1688  else if(type == eIMEX)
1689  {
1690  if(m_numstages==1)
1691  {
1692  m_T= t_old[0]+timestep;
1693  }
1694  else
1695  {
1696  m_T= t_old[0];
1697  for(int j=0; j<=i; ++j)
1698  {
1699  m_T += A(i,j)*timestep;
1700  }
1701  }
1702 
1703  if(fabs(A(i,i)) > NekConstants::kNekZeroTol)
1704  {
1705  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
1706 
1707  for(k = 0; k < m_nvar; k++)
1708  {
1709  Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
1710  Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),
1711  m_F[i][k],1,m_F[i][k],1);
1712  }
1713  }
1714  op.DoOdeRhs(m_Y, m_F_IMEX[i], m_T);
1715  }
1716  else if( type == eExplicit)
1717  {
1718  // Avoid projecting the same solution twice
1719  if( ! ((i==0) && m_firstStageEqualsOldSolution) )
1720  {
1721  // ensure solution is in correct space
1722  op.DoProjection(m_Y,m_Y,m_T);
1723  }
1724  op.DoOdeRhs(m_Y, m_F[i], m_T);
1725  }
1726  }
1727 
1728  // Next, the solution vector y at the new time level will
1729  // be calculated.
1730  //
1731  // For multi-step methods, this includes updating the
1732  // values of the auxiliary parameters
1733  //
1734  // The loop below calculates the solution at the new time
1735  // level
1736  //
1737  // If last stage equals the new solution, the new solution
1738  // needs not be calculated explicitly but can simply be
1739  // copied. This saves a solve.
1740  int i_start = 0;
1741  if( (m_lastStageEqualsNewSolution) )
1742  {
1743  for(k = 0; k < m_nvar; k++)
1744  {
1745  Vmath::Vcopy(m_npoints,m_Y[k],1,y_new[0][k],1);
1746  }
1747 
1748  if (m_numstages==1 && type == eIMEX)
1749  {
1750  t_new[0] = t_old[0]+timestep;
1751  }
1752  else
1753  {
1754  t_new[0] = B(0,0)*timestep;
1755  for(j = 1; j < m_numstages; j++)
1756  {
1757  t_new[0] += B(0,j)*timestep;
1758  }
1759  for(j = 0; j < m_numsteps; j++)
1760  {
1761  t_new[0] += V(0,j)*t_old[j];
1762  }
1763  }
1764  i_start = 1;
1765  }
1766 
1767  for(i = i_start; i < m_numsteps; i++)
1768  {
1769  // The solution at the new time level is a linear
1770  // combination of:
1771  // 1: the stage derivatives
1772  for(k = 0; k < m_nvar; k++)
1773  {
1774  Vmath::Smul(m_npoints,timestep*B(i,0),m_F[0][k],1,
1775  y_new[i][k],1);
1776 
1777  if(type == eIMEX)
1778  {
1779  Vmath::Svtvp(m_npoints,timestep*B_IMEX(i,0),
1780  m_F_IMEX[0][k],1,y_new[i][k],1,
1781  y_new[i][k],1);
1782  }
1783  }
1784  if(m_numstages != 1 || type != eIMEX)
1785  {
1786  t_new[i] = B(i,0)*timestep;
1787  }
1788 
1789 
1790  for(j = 1; j < m_numstages; j++)
1791  {
1792  for(k = 0; k < m_nvar; k++)
1793  {
1794  Vmath::Svtvp(m_npoints,timestep*B(i,j),m_F[j][k],1,
1795  y_new[i][k],1,y_new[i][k],1);
1796 
1797  if(type == eIMEX)
1798  {
1799  Vmath::Svtvp(m_npoints,timestep*B_IMEX(i,j),
1800  m_F_IMEX[j][k],1,y_new[i][k],1,
1801  y_new[i][k],1);
1802  }
1803  }
1804  if(m_numstages != 1 || type != eIMEX)
1805  {
1806  t_new[i] += B(i,j)*timestep;
1807  }
1808  }
1809 
1810  // 2: the imported multi-step solution of the previous
1811  // time level
1812  for(j = 0; j < m_numsteps; j++)
1813  {
1814  for(k = 0; k < m_nvar; k++)
1815  {
1816  Vmath::Svtvp(m_npoints,V(i,j),y_old[j][k],1,
1817  y_new[i][k],1,y_new[i][k],1);
1818  }
1819  if(m_numstages != 1 || type != eIMEX)
1820  {
1821  t_new[i] += V(i,j)*t_old[j];
1822  }
1823  }
1824  }
1825 
1826  // Ensure that the new solution is projected if necessary
1827  if(type == eExplicit)
1828  {
1829  op.DoProjection(y_new[0],y_new[0],t_new[0]);
1830  }
1831  }
1832 
1833  bool TimeIntegrationScheme::CheckIfFirstStageEqualsOldSolution(const Array<OneD, const Array<TwoD, NekDouble> >& A,
1834  const Array<OneD, const Array<TwoD, NekDouble> >& B,
1835  const Array<TwoD, const NekDouble>& U,
1836  const Array<TwoD, const NekDouble>& V) const
1837  {
1838  boost::ignore_unused(B, V);
1839 
1840  int i,m;
1841  // First stage equals old solution if:
1842  // 1. the first row of the coefficient matrix A consists of zeros
1843  // 2. U[0][0] is equal to one and all other first row entries of U are zero
1844 
1845  // 1. Check first condition
1846  for(m = 0; m < A.num_elements(); m++)
1847  {
1848  for(i = 0; i < m_numstages; i++)
1849  {
1850  if( fabs(A[m][0][i]) > NekConstants::kNekZeroTol )
1851  {
1852  return false;
1853  }
1854  }
1855  }
1856 
1857  // 2. Check second condition
1858  if( fabs(U[0][0] - 1.0) > NekConstants::kNekZeroTol )
1859  {
1860  return false;
1861  }
1862  for(i = 1; i < m_numsteps; i++)
1863  {
1864  if( fabs(U[0][i]) > NekConstants::kNekZeroTol )
1865  {
1866  return false;
1867  }
1868  }
1869 
1870  return true;
1871  }
1872 
1873  bool TimeIntegrationScheme::CheckIfLastStageEqualsNewSolution(const Array<OneD, const Array<TwoD, NekDouble> >& A,
1874  const Array<OneD, const Array<TwoD, NekDouble> >& B,
1875  const Array<TwoD, const NekDouble>& U,
1876  const Array<TwoD, const NekDouble>& V) const
1877  {
1878  int i,m;
1879  // Last stage equals new solution if:
1880  // 1. the last row of the coefficient matrix A is equal to the first row of matrix B
1881  // 2. the last row of the coefficient matrix U is equal to the first row of matrix V
1882 
1883  // 1. Check first condition
1884  for(m = 0; m < A.num_elements(); m++)
1885  {
1886  for(i = 0; i < m_numstages; i++)
1887  {
1888  if( fabs(A[m][m_numstages-1][i]-B[m][0][i]) > NekConstants::kNekZeroTol )
1889  {
1890  return false;
1891  }
1892  }
1893  }
1894 
1895  // 2. Check second condition
1896  for(i = 0; i < m_numsteps; i++)
1897  {
1898  if( fabs(U[m_numstages-1][i]-V[0][i]) > NekConstants::kNekZeroTol )
1899  {
1900  return false;
1901  }
1902  }
1903 
1904  return true;
1905  }
1906 
1907  bool TimeIntegrationScheme::CheckTimeIntegrateArguments(const NekDouble timestep,
1908  ConstTripleArray &y_old ,
1909  ConstSingleArray &t_old ,
1910  TripleArray &y_new ,
1911  SingleArray &t_new ,
1912  const TimeIntegrationSchemeOperators &op) const
1913  {
1914  boost::ignore_unused(timestep, y_old, t_old, y_new, t_new, op);
1915 
1916  // Check if arrays are all of consistent size
1917  ASSERTL1(y_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1918  ASSERTL1(y_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1919 
1920  ASSERTL1(y_old[0]. num_elements()==y_new[0]. num_elements(),"Non-matching number of variables.");
1921  ASSERTL1(y_old[0][0].num_elements()==y_new[0][0].num_elements(),"Non-matching number of coefficients.");
1922 
1923  ASSERTL1(t_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1924  ASSERTL1(t_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1925 
1926  return true;
1927  }
1928 
1929  std::ostream& operator<<(std::ostream& os, const TimeIntegrationSchemeSharedPtr& rhs)
1930  {
1931  return operator<<(os,*rhs);
1932  }
1933 
1934  std::ostream& operator<<(std::ostream& os, const TimeIntegrationScheme& rhs)
1935  {
1936  int i,j;
1937  int r = rhs.GetNsteps();
1938  int s = rhs.GetNstages();
1939  TimeIntegrationSchemeType type = rhs.GetIntegrationSchemeType();
1940 
1941  int oswidth = 9;
1942  int osprecision = 6;
1943 
1944  os << "Time Integration Scheme: " << TimeIntegrationMethodMap[rhs.GetIntegrationMethod()] << std::endl;
1945  os << "- number of steps: " << r << std::endl;
1946  os << "- number of stages: " << s << std::endl;
1947  os << "- type of scheme: " << TimeIntegrationSchemeTypeMap[rhs.GetIntegrationSchemeType()] << std::endl;
1948  os << "General linear method tableau: " << std::endl;
1949 
1950  for(i = 0; i < s; i++)
1951  {
1952  for(j = 0; j < s; j++)
1953  {
1954  os.width(oswidth);
1955  os.precision(osprecision);
1956  os << std::right << rhs.A(i,j) << " ";
1957  }
1958  if(type == eIMEX)
1959  {
1960  os << " '";
1961  for(j = 0; j < s; j++)
1962  {
1963  os.width(oswidth);
1964  os.precision(osprecision);
1965  os << std::right << rhs.A_IMEX(i,j) << " ";
1966  }
1967  }
1968  os << " |";
1969 
1970  for(j = 0; j < r; j++)
1971  {
1972  os.width(oswidth);
1973  os.precision(osprecision);
1974  os << std::right << rhs.U(i,j);
1975  }
1976  os << std::endl;
1977  }
1978  int imexflag = (type == eIMEX)?2:1;
1979  for(int i = 0; i < (r+imexflag*s)*(oswidth+1)+imexflag*2-1; i++)
1980  {
1981  os << "-";
1982  }
1983  os << std::endl;
1984  for(i = 0; i < r; i++)
1985  {
1986  for(j = 0; j < s; j++)
1987  {
1988  os.width(oswidth);
1989  os.precision(osprecision);
1990  os << std::right << rhs.B(i,j) << " ";
1991  }
1992  if(type == eIMEX)
1993  {
1994  os << " '";
1995  for(j = 0; j < s; j++)
1996  {
1997  os.width(oswidth);
1998  os.precision(osprecision);
1999  os << std::right << rhs.B_IMEX(i,j) << " ";
2000  }
2001  }
2002  os << " |";
2003 
2004  for(j = 0; j < r; j++)
2005  {
2006  os.width(oswidth);
2007  os.precision(osprecision);
2008  os << std::right << rhs.V(i,j);
2009  }
2010  os << std::endl;
2011  }
2012  return os;
2013  }
2014  }
2015 }
Nektar::LibUtilities::TimeIntegrationScheme::A
NekDouble A(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:412
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:1873
Nektar::LibUtilities::TimeIntegrationScheme::B
NekDouble B(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:417
Nektar::LibUtilities::TimeIntegrationSchemeSharedPtr
std::shared_ptr< TimeIntegrationScheme > TimeIntegrationSchemeSharedPtr
Definition: TimeIntegrationScheme.h:52
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:209
Nektar::LibUtilities::eBDFImplicitOrder1
BDF multi-step scheme of order 1 (implicit)
Definition: TimeIntegrationScheme.h:74
Nektar::LibUtilities::eAdamsBashforthOrder2
Adams-Bashforth Forward multi-step scheme of order 2.
Definition: TimeIntegrationScheme.h:69
Nektar::LibUtilities::TimeIntegrationScheme::m_B
Array< OneD, Array< TwoD, NekDouble > > m_B
Definition: TimeIntegrationScheme.h:553
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoImplicitSolve
void DoImplicitSolve(InArrayType &inarray, OutArrayType &outarray, const NekDouble time, const NekDouble lambda) const
Definition: TimeIntegrationScheme.h:263
Nektar::LibUtilities::TimeIntegrationScheme::m_nvar
int m_nvar
bool to identify if array has been initialised
Definition: TimeIntegrationScheme.h:559
Nektar::LibUtilities::eClassicalRungeKutta4
Runge-Kutta multi-stage scheme 4th order explicit (old name)
Definition: TimeIntegrationScheme.h:76
Nektar::LibUtilities::eIMEXdirk_1_2_2
Implicit-Explicit Midpoint IMEX DIRK(1,2,2)
Definition: TimeIntegrationScheme.h:97
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:1907
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:83
Nektar::LibUtilities::TimeIntegrationSchemeTypeMap
const char *const TimeIntegrationSchemeTypeMap[]
Definition: TimeIntegrationScheme.h:155
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs
void DoOdeRhs(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
Definition: TimeIntegrationScheme.h:231
Nektar::LibUtilities::TimeIntegrationScheme::m_tmp
DoubleArray m_tmp
Array containing the stage values.
Definition: TimeIntegrationScheme.h:562
Nektar::LibUtilities::TimeIntegrationScheme::V
NekDouble V(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:427
VmathArray.hpp
Nektar::LibUtilities::eRungeKutta4
Classical RungeKutta4 method (new name for eClassicalRungeKutta4)
Definition: TimeIntegrationScheme.h:77
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:488
Nektar::LibUtilities::TimeIntegrationScheme::InitializeScheme
TimeIntegrationSolutionSharedPtr InitializeScheme(const NekDouble timestep, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)
This function initialises the time integration scheme.
Definition: TimeIntegrationScheme.cpp:1239
Nektar::LibUtilities::eIMEXdirk_4_4_3
L-stable, four stage, third order IMEX DIRK(4,4,3)
Definition: TimeIntegrationScheme.h:102
Nektar::LibUtilities::eIMEXdirk_1_2_1
Forward-Backward Euler IMEX DIRK(1,2,1)
Definition: TimeIntegrationScheme.h:96
Nektar::LibUtilities::eIMEX
Implicit Explicit General Linear Method.
Definition: TimeIntegrationScheme.h:151
Nektar::LibUtilities::TimeIntegrationMethodMap
const char *const TimeIntegrationMethodMap[]
Definition: TimeIntegrationScheme.h:106
Nektar::LibUtilities::eRungeKutta3_SSP
Nonlinear SSP RungeKutta3 explicit.
Definition: TimeIntegrationScheme.h:79
Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationMethod
TimeIntegrationMethod GetIntegrationMethod() const
Definition: TimeIntegrationScheme.h:402
Nektar::LibUtilities::TimeIntegrationScheme::m_Y
DoubleArray m_Y
The size of inner data which is stored for reuse.
Definition: TimeIntegrationScheme.h:561
Nektar::LibUtilities::TimeIntegrationSchemeManager
TimeIntegrationSchemeManagerT & TimeIntegrationSchemeManager(void)
Definition: TimeIntegrationScheme.cpp:44
Nektar::LibUtilities::eAdamsMoultonOrder2
Adams-Moulton Forward multi-step scheme of order 2.
Definition: TimeIntegrationScheme.h:73
Nektar::LibUtilities::eAdamsBashforthOrder3
Adams-Bashforth Forward multi-step scheme of order 3.
Definition: TimeIntegrationScheme.h:70
Nektar::lhs
StandardMatrixTag & lhs
Definition: MatrixOperations.cpp:95
Nektar::LibUtilities::TimeIntegrationScheme::A_IMEX
NekDouble A_IMEX(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:432
Nektar::LibUtilities::operator==
bool operator==(const BasisKey &x, const BasisKey &y)
Definition: Foundations/Basis.cpp:871
Nektar::LibUtilities::eCNAB
Crank-Nicolson/Adams-Bashforth Order 2 (CNAB)
Definition: TimeIntegrationScheme.h:92
Nektar::LibUtilities::eIMEXOrder2
IMEX 2nd order scheme using Backward Different Formula & Extrapolation.
Definition: TimeIntegrationScheme.h:85
Nektar::LibUtilities::operator<<
std::ostream & operator<<(std::ostream &os, const BasisKey &rhs)
Definition: Foundations/Basis.cpp:87
Nektar::LibUtilities::TimeIntegrationScheme::m_F
TripleArray m_F
explicit right hand side of each stage equation
Definition: TimeIntegrationScheme.h:564
Nektar::LibUtilities::eDiagonallyImplicit
Diagonally implicit scheme (e.g. the DIRK schemes)
Definition: TimeIntegrationScheme.h:150
Nektar::LibUtilities::eRungeKutta2
Classical RungeKutta2 method (new name for eMidpoint)
Definition: TimeIntegrationScheme.h:89
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:50
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:1182
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:216
Nektar::LibUtilities::eAdamsMoultonOrder1
Adams-Moulton Forward multi-step scheme of order 1.
Definition: TimeIntegrationScheme.h:72
Nektar::LibUtilities::eAdamsBashforthOrder1
Adams-Bashforth Forward multi-step scheme of order 1.
Definition: TimeIntegrationScheme.h:68
Nektar::LibUtilities::TimeIntegrationSolution::TimeIntegrationSolution
TimeIntegrationSolution(const TimeIntegrationSchemeKey &key, const DoubleArray &y, const NekDouble time, const NekDouble timestep)
Definition: TimeIntegrationScheme.cpp:74
Nektar::LibUtilities::TimeIntegrationScheme::U
NekDouble U(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:422
Nektar::LibUtilities::eIMEXdirk_3_4_3
L-stable, three stage, third order IMEX DIRK(3,4,3)
Definition: TimeIntegrationScheme.h:101
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoProjection
void DoProjection(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
Definition: TimeIntegrationScheme.h:255
Nektar::LibUtilities::TimeIntegrationScheme::m_npoints
int m_npoints
The number of variables in integration scheme.
Definition: TimeIntegrationScheme.h:560
Nektar::LibUtilities::eRungeKutta2_SSP
Nonlinear SSP RungeKutta2 explicit (surrogate for eRungeKutta2_ImprovedEuler)
Definition: TimeIntegrationScheme.h:81
Nektar::LibUtilities::TimeIntegrationScheme::TimeIntegrate
ConstDoubleArray & TimeIntegrate(const NekDouble timestep, TimeIntegrationSolutionSharedPtr &y, const TimeIntegrationSchemeOperators &op)
Explicit integration of an ODE.
Definition: TimeIntegrationScheme.cpp:1286
Nektar::LibUtilities::eRungeKutta2_ImprovedEuler
Improved RungeKutta2 explicit (old name meaning Heun&#39;s method)
Definition: TimeIntegrationScheme.h:80
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::LibUtilities::eIMEXdirk_1_1_1
Forward-Backward Euler IMEX DIRK(1,1,1)
Definition: TimeIntegrationScheme.h:95
Nektar::LibUtilities::TimeIntegrationScheme::DoubleArray
Array< OneD, Array< OneD, NekDouble > > DoubleArray
Definition: TimeIntegrationScheme.h:385
TimeIntegrationScheme.h
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::LibUtilities::TimeIntegrationScheme::B_IMEX
NekDouble B_IMEX(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:437
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:1833
Nektar::LibUtilities::operator<
bool operator<(const BasisKey &lhs, const BasisKey &rhs)
Definition: Foundations/Basis.cpp:51
Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeType
TimeIntegrationSchemeType GetIntegrationSchemeType() const
Definition: TimeIntegrationScheme.h:407
Nektar::LibUtilities::TimeIntegrationScheme::Create
static std::shared_ptr< TimeIntegrationScheme > Create(const TimeIntegrationSchemeKey &key)
Definition: TimeIntegrationScheme.cpp:142
Nektar::LibUtilities::TimeIntegrationSchemeKey::opLess::operator()
bool operator()(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs) const
Definition: TimeIntegrationScheme.cpp:62
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:346
NektarUnivConsts.hpp
Blas::Dscal
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:125
Nektar::LibUtilities::eBDFImplicitOrder2
BDF multi-step scheme of order 2 (implicit)
Definition: TimeIntegrationScheme.h:75
Nektar::LibUtilities::eIMEXdirk_2_2_2
L-stable, two stage, second order IMEX DIRK(2,2,2)
Definition: TimeIntegrationScheme.h:98
Nektar::LibUtilities::TimeIntegrationSchemeKey::m_method
TimeIntegrationMethod m_method
integration method
Definition: TimeIntegrationScheme.h:344
Nektar::LibUtilities::eDIRKOrder2
Diagonally Implicit Runge Kutta scheme of order 2.
Definition: TimeIntegrationScheme.h:90
Nektar::LibUtilities::eIMEXOrder3
IMEX 3rd order scheme using Backward Different Formula & Extrapolation.
Definition: TimeIntegrationScheme.h:86
Nektar::LibUtilities::eIMEXdirk_2_3_2
L-stable, three stage, third order IMEX DIRK(3,4,3)
Definition: TimeIntegrationScheme.h:99
Nektar::LibUtilities::eDIRKOrder3
Diagonally Implicit Runge Kutta scheme of order 3.
Definition: TimeIntegrationScheme.h:91
Nektar::LibUtilities::NekManager::RegisterGlobalCreator
bool RegisterGlobalCreator(const CreateFuncType &createFunc)
Register the Global Create Function. The return value is just to facilitate calling statically...
Definition: NekManager.hpp:176
Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeKey
const TimeIntegrationSchemeKey & GetIntegrationSchemeKey() const
Definition: TimeIntegrationScheme.h:397
Nektar::LibUtilities::eIMEXdirk_2_3_3
L-stable, two stage, third order IMEX DIRK(2,3,3)
Definition: TimeIntegrationScheme.h:100
Nektar::LibUtilities::eIMEXOrder4
IMEX 4th order scheme using Backward Different Formula & Extrapolation.
Definition: TimeIntegrationScheme.h:87
Nektar::rhs
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs
Definition: MatrixOperations.cpp:97
Nektar::LibUtilities::TimeIntegrationScheme::m_A
Array< OneD, Array< TwoD, NekDouble > > m_A
Definition: TimeIntegrationScheme.h:552
Nektar::LibUtilities::eMCNAB
Modified Crank-Nicolson/Adams-Bashforth Order 2 (MCNAB)
Definition: TimeIntegrationScheme.h:94
Nektar::LibUtilities::eIMEXOrder1
IMEX 1st order scheme using Euler Backwards/Euler Forwards.
Definition: TimeIntegrationScheme.h:84
Nektar::LibUtilities::eMidpoint
midpoint method (old name)
Definition: TimeIntegrationScheme.h:88
Nektar::LibUtilities::TimeIntegrationSolution::GetTimeLevelOffset
const Array< OneD, const unsigned int > & GetTimeLevelOffset()
Definition: TimeIntegrationScheme.h:738
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376
Nektar::LibUtilities::TimeIntegrationScheme::m_T
NekDouble m_T
Used to store the Explicit stage derivative of IMEX schemes.
Definition: TimeIntegrationScheme.h:567
Nektar::LibUtilities::TimeIntegrationScheme::TripleArray
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > TripleArray
Definition: TimeIntegrationScheme.h:383
Nektar::LibUtilities::TimeIntegrationSolutionSharedPtr
std::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
Definition: TimeIntegrationScheme.h:58
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
Nektar::LibUtilities::eAdamsBashforthOrder4
Adams-Bashforth Forward multi-step scheme of order 4.
Definition: TimeIntegrationScheme.h:71
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::LibUtilities::TimeIntegrationScheme::m_F_IMEX
TripleArray m_F_IMEX
Array corresponding to the stage Derivatives.
Definition: TimeIntegrationScheme.h:565
Generated on Wed Jan 20 2021 03:33:42 for Nektar++ by   doxygen 1.8.13