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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: implementation of time integration key class
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 #include <LibUtilities/TimeIntegration/TimeIntegrationScheme.h>
36 #include <LibUtilities/BasicUtils/VmathArray.hpp>
37 #include <LibUtilities/BasicConst/NektarUnivConsts.hpp>
38 #include <iostream>
39 #include <math.h>
40 
41 namespace Nektar
42 {
43  namespace LibUtilities
44  {
45  TimeIntegrationSchemeManagerT &TimeIntegrationSchemeManager(void)
46  {
47  TimeIntegrationSchemeManagerT& m = Loki::SingletonHolder<TimeIntegrationSchemeManagerT>::Instance();
48  m.RegisterGlobalCreator(TimeIntegrationScheme::Create);
49  return m;
50  }
51 
52 
53  bool operator==(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs)
54  {
55  return (lhs.m_method == rhs.m_method);
56  }
57 
58  bool operator<(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs)
59  {
60  return (lhs.m_method < rhs.m_method);
61  }
62 
63  bool TimeIntegrationSchemeKey::opLess::operator()(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs) const
64  {
65  return (lhs.m_method < rhs.m_method);
66  }
67 
68  std::ostream& operator<<(std::ostream& os, const TimeIntegrationSchemeKey& rhs)
69  {
70  os << "Time Integration Scheme: " << TimeIntegrationMethodMap[rhs.GetIntegrationMethod()] << endl;
71 
72  return os;
73  }
74 
75  TimeIntegrationSolution::TimeIntegrationSolution(const TimeIntegrationSchemeKey &key,
76  const DoubleArray& y,
77  const NekDouble time,
78  const NekDouble timestep):
79  m_scheme(TimeIntegrationSchemeManager()[key]),
80  m_solVector(m_scheme->GetNsteps()),
81  m_t(m_scheme->GetNsteps())
82  {
83  m_solVector[0] = y;
84  m_t[0] = time;
85 
86  int nsteps = m_scheme->GetNsteps();
87  int nvar = y.num_elements();
88  int npoints = y[0].num_elements();
89  int nMultiStepVals = m_scheme->GetNmultiStepValues();
90  const Array<OneD, const unsigned int>& timeLevels = GetTimeLevelOffset();
91  for(int i = 1; i < nsteps; i++)
92  {
93  m_solVector[i] = Array<OneD, Array<OneD, NekDouble> >(nvar);
94  for(int j = 0; j < nvar; j++)
95  {
96  m_solVector[i][j] = Array<OneD,NekDouble>(npoints,0.0);
97  }
98  if(i < nMultiStepVals)
99  {
100  m_t[i] = time - i*timestep*timeLevels[i];
101  }
102  else
103  {
104  m_t[i] = timestep;
105  }
106  }
107  }
108 
109  TimeIntegrationSolution::TimeIntegrationSolution(const TimeIntegrationSchemeKey &key,
110  const TripleArray& y,
111  const Array<OneD, NekDouble>& t):
112  m_scheme(TimeIntegrationSchemeManager()[key]),
113  m_solVector(y),
114  m_t(t)
115  {
116  ASSERTL1(y.num_elements()==m_scheme->GetNsteps(),"Amount of Entries does not match number of (multi-) steps");
117  }
118 
119  TimeIntegrationSolution::TimeIntegrationSolution(const TimeIntegrationSchemeKey &key,
120  unsigned int nvar,
121  unsigned int npoints):
122  m_scheme(TimeIntegrationSchemeManager()[key]),
123  m_solVector(m_scheme->GetNsteps()),
124  m_t(m_scheme->GetNsteps())
125  {
126  for(int i = 0; i < m_scheme->GetNsteps(); i++)
127  {
128  m_solVector[i] = Array<OneD, Array<OneD, NekDouble> >(nvar);
129  for(int j = 0; j < nvar; j++)
130  {
131  m_solVector[i][j] = Array<OneD,NekDouble>(npoints);
132  }
133  }
134  }
135 
136  TimeIntegrationSolution::TimeIntegrationSolution(const TimeIntegrationSchemeKey &key):
137  m_scheme(TimeIntegrationSchemeManager()[key]),
138  m_solVector(m_scheme->GetNsteps()),
139  m_t(m_scheme->GetNsteps())
140  {
141  }
142 
143  TimeIntegrationSchemeSharedPtr TimeIntegrationScheme::Create(const TimeIntegrationSchemeKey &key)
144  {
145  TimeIntegrationSchemeSharedPtr returnval(
146  MemoryManager<TimeIntegrationScheme>::AllocateSharedPtr(key));
147  return returnval;
148  }
149 
150  TimeIntegrationScheme::TimeIntegrationScheme(const TimeIntegrationSchemeKey &key):
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  }
1045 
1046 
1047  bool TimeIntegrationScheme::
1048  VerifyIntegrationSchemeType(TimeIntegrationSchemeType type,
1049  const Array<OneD, const Array<TwoD, NekDouble> >& A,
1050  const Array<OneD, const Array<TwoD, NekDouble> >& B,
1051  const Array<TwoD, const NekDouble>& U,
1052  const Array<TwoD, const NekDouble>& V) const
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  }
1101 
1102  TimeIntegrationSolutionSharedPtr
1103  TimeIntegrationScheme::InitializeScheme(const NekDouble timestep,
1104  ConstDoubleArray &y_0 ,
1105  const NekDouble time ,
1106  const TimeIntegrationSchemeOperators &op)
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  }
1142 
1143  TimeIntegrationScheme::ConstDoubleArray&
1144  TimeIntegrationScheme::TimeIntegrate(const NekDouble timestep,
1145  TimeIntegrationSolutionSharedPtr &solvector,
1146  const TimeIntegrationSchemeOperators &op)
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  }
1368 
1369  void TimeIntegrationScheme::TimeIntegrate(const NekDouble timestep,
1370  ConstTripleArray &y_old ,
1371  ConstSingleArray &t_old ,
1372  TripleArray &y_new ,
1373  SingleArray &t_new ,
1374  const TimeIntegrationSchemeOperators &op)
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  i_start = 1;
1618  }
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  }
1686 
1687  bool TimeIntegrationScheme::CheckIfFirstStageEqualsOldSolution(const Array<OneD, const Array<TwoD, NekDouble> >& A,
1688  const Array<OneD, const Array<TwoD, NekDouble> >& B,
1689  const Array<TwoD, const NekDouble>& U,
1690  const Array<TwoD, const NekDouble>& V) const
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  }
1724 
1725  bool TimeIntegrationScheme::CheckIfLastStageEqualsNewSolution(const Array<OneD, const Array<TwoD, NekDouble> >& A,
1726  const Array<OneD, const Array<TwoD, NekDouble> >& B,
1727  const Array<TwoD, const NekDouble>& U,
1728  const Array<TwoD, const NekDouble>& V) const
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  }
1758 
1759  bool TimeIntegrationScheme::CheckTimeIntegrateArguments(const NekDouble timestep,
1760  ConstTripleArray &y_old ,
1761  ConstSingleArray &t_old ,
1762  TripleArray &y_new ,
1763  SingleArray &t_new ,
1764  const TimeIntegrationSchemeOperators &op) const
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  }
1778 
1779  std::ostream& operator<<(std::ostream& os, const TimeIntegrationSchemeSharedPtr& rhs)
1780  {
1781  return operator<<(os,*rhs);
1782  }
1783 
1784  std::ostream& operator<<(std::ostream& os, const TimeIntegrationScheme& rhs)
1785  {
1786  int i,j;
1787  int r = rhs.GetNsteps();
1788  int s = rhs.GetNstages();
1789  TimeIntegrationSchemeType type = rhs.GetIntegrationSchemeType();
1790 
1791  int oswidth = 9;
1792  int osprecision = 6;
1793 
1794  os << "Time Integration Scheme: " << TimeIntegrationMethodMap[rhs.GetIntegrationMethod()] << endl;
1795  os << "- number of steps: " << r << endl;
1796  os << "- number of stages: " << s << endl;
1797  os << "- type of scheme: " << TimeIntegrationSchemeTypeMap[rhs.GetIntegrationSchemeType()] << endl;
1798  os << "General linear method tableau: " << endl;
1799 
1800  for(i = 0; i < s; i++)
1801  {
1802  for(j = 0; j < s; j++)
1803  {
1804  os.width(oswidth);
1805  os.precision(osprecision);
1806  os << right << rhs.A(i,j) << " ";
1807  }
1808  if(type == eIMEX)
1809  {
1810  os << " '";
1811  for(j = 0; j < s; j++)
1812  {
1813  os.width(oswidth);
1814  os.precision(osprecision);
1815  os << right << rhs.A_IMEX(i,j) << " ";
1816  }
1817  }
1818  os << " |";
1819 
1820  for(j = 0; j < r; j++)
1821  {
1822  os.width(oswidth);
1823  os.precision(osprecision);
1824  os << right << rhs.U(i,j);
1825  }
1826  os << endl;
1827  }
1828  int imexflag = (type == eIMEX)?2:1;
1829  for(int i = 0; i < (r+imexflag*s)*(oswidth+1)+imexflag*2-1; i++)
1830  {
1831  os << "-";
1832  }
1833  os << endl;
1834  for(i = 0; i < r; i++)
1835  {
1836  for(j = 0; j < s; j++)
1837  {
1838  os.width(oswidth);
1839  os.precision(osprecision);
1840  os << right << rhs.B(i,j) << " ";
1841  }
1842  if(type == eIMEX)
1843  {
1844  os << " '";
1845  for(j = 0; j < s; j++)
1846  {
1847  os.width(oswidth);
1848  os.precision(osprecision);
1849  os << right << rhs.B_IMEX(i,j) << " ";
1850  }
1851  }
1852  os << " |";
1853 
1854  for(j = 0; j < r; j++)
1855  {
1856  os.width(oswidth);
1857  os.precision(osprecision);
1858  os << right << rhs.V(i,j);
1859  }
1860  os << endl;
1861  }
1862  return os;
1863  }
1864  }
1865 }
Nektar::LibUtilities::TimeIntegrationScheme::B_IMEX
NekDouble B_IMEX(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:420
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
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::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::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::LibUtilities::TimeIntegrationScheme::m_B
Array< OneD, Array< TwoD, NekDouble > > m_B
Definition: TimeIntegrationScheme.h:536
Nektar::LibUtilities::TimeIntegrationScheme::m_nvar
int m_nvar
bool to identify if array has been initialised
Definition: TimeIntegrationScheme.h:542
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::TimeIntegrationSchemeOperators::DoImplicitSolve
void DoImplicitSolve(InArrayType &inarray, OutArrayType &outarray, const NekDouble time, const NekDouble lambda) const
Definition: TimeIntegrationScheme.h:246
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Nektar::LibUtilities::TimeIntegrationSchemeTypeMap
const char *const TimeIntegrationSchemeTypeMap[]
Definition: TimeIntegrationScheme.h:148
Nektar::LibUtilities::TimeIntegrationScheme::m_tmp
DoubleArray m_tmp
Array containing the stage values.
Definition: TimeIntegrationScheme.h:545
VmathArray.hpp
Nektar::LibUtilities::eRungeKutta4
Classical RungeKutta4 method (new name for eClassicalRungeKutta4)
Definition: TimeIntegrationScheme.h:75
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:471
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:1103
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::TimeIntegrationMethodMap
const char *const TimeIntegrationMethodMap[]
Definition: TimeIntegrationScheme.h:102
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoOdeRhs
void DoOdeRhs(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
Definition: TimeIntegrationScheme.h:214
Nektar::LibUtilities::eRungeKutta3_SSP
Nonlinear SSP RungeKutta3 explicit.
Definition: TimeIntegrationScheme.h:76
Nektar::LibUtilities::TimeIntegrationSchemeKey::opLess::operator()
bool operator()(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs) const
Definition: TimeIntegrationScheme.cpp:63
Nektar::LibUtilities::TimeIntegrationScheme::m_Y
DoubleArray m_Y
The size of inner data which is stored for reuse.
Definition: TimeIntegrationScheme.h:544
Nektar::LibUtilities::TimeIntegrationSchemeManager
TimeIntegrationSchemeManagerT & TimeIntegrationSchemeManager(void)
Definition: TimeIntegrationScheme.cpp:45
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::lhs
StandardMatrixTag & lhs
Definition: MatrixOperations.cpp:96
Nektar::LibUtilities::operator==
bool operator==(const BasisKey &x, const BasisKey &y)
Definition: Basis.cpp:703
Nektar::LibUtilities::eCNAB
Crank-Nicolson/Adams-Bashforth Order 2 (CNAB)
Definition: TimeIntegrationScheme.h:88
Nektar::LibUtilities::TimeIntegrationScheme::GetSecondDim
int GetSecondDim(ConstTripleArray &y) const
Definition: TimeIntegrationScheme.h:589
Nektar::LibUtilities::eIMEXOrder2
IMEX 2nd order scheme using Backward Different Formula & Extrapolation.
Definition: TimeIntegrationScheme.h:82
Nektar::LibUtilities::operator<<
std::ostream & operator<<(std::ostream &os, const BasisKey &rhs)
Definition: Basis.cpp:85
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::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::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::LibUtilities::eAdamsMoultonOrder1
Adams-Moulton Forward multi-step scheme of order 1.
Definition: TimeIntegrationScheme.h:70
Nektar::LibUtilities::TimeIntegrationScheme::A_IMEX
NekDouble A_IMEX(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:415
Nektar::LibUtilities::eAdamsBashforthOrder1
Adams-Bashforth Forward multi-step scheme of order 1.
Definition: TimeIntegrationScheme.h:67
Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeKey
const TimeIntegrationSchemeKey & GetIntegrationSchemeKey() const
Definition: TimeIntegrationScheme.h:380
Nektar::LibUtilities::TimeIntegrationSolution::TimeIntegrationSolution
TimeIntegrationSolution(const TimeIntegrationSchemeKey &key, const DoubleArray &y, const NekDouble time, const NekDouble timestep)
Definition: TimeIntegrationScheme.cpp:75
Nektar::LibUtilities::eIMEXdirk_3_4_3
L-stable, three stage, third order IMEX DIRK(3,4,3)
Definition: TimeIntegrationScheme.h:97
Nektar::LibUtilities::TimeIntegrationScheme::m_npoints
int m_npoints
The number of variables in integration scheme.
Definition: TimeIntegrationScheme.h:543
Nektar::LibUtilities::eRungeKutta2_SSP
Nonlinear SSP RungeKutta2 explicit (surrogate for eRungeKutta2_ImprovedEuler)
Definition: TimeIntegrationScheme.h:78
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::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::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::LibUtilities::TimeIntegrationScheme::DoubleArray
Array< OneD, Array< OneD, NekDouble > > DoubleArray
Definition: TimeIntegrationScheme.h:368
TimeIntegrationScheme.h
Nektar::LibUtilities::TimeIntegrationScheme::A
NekDouble A(const unsigned int i, const unsigned int j) const
Definition: TimeIntegrationScheme.h:395
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationMethod
TimeIntegrationMethod GetIntegrationMethod() const
Definition: TimeIntegrationScheme.h:385
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DoProjection
void DoProjection(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
Definition: TimeIntegrationScheme.h:238
Nektar::LibUtilities::operator<
bool operator<(const BasisKey &lhs, const BasisKey &rhs)
Definition: Basis.cpp:49
Nektar::LibUtilities::TimeIntegrationScheme::GetIntegrationSchemeType
TimeIntegrationSchemeType GetIntegrationSchemeType() const
Definition: TimeIntegrationScheme.h:390
Nektar::LibUtilities::TimeIntegrationSchemeSharedPtr
boost::shared_ptr< TimeIntegrationScheme > TimeIntegrationSchemeSharedPtr
Definition: TimeIntegrationScheme.h:51
Nektar::LibUtilities::TimeIntegrationScheme::Create
static boost::shared_ptr< TimeIntegrationScheme > Create(const TimeIntegrationSchemeKey &key)
Definition: TimeIntegrationScheme.cpp:143
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
NektarUnivConsts.hpp
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::TimeIntegrationSchemeKey::m_method
TimeIntegrationMethod m_method
integration method
Definition: TimeIntegrationScheme.h:327
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::TimeIntegrationSolutionSharedPtr
boost::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
Definition: TimeIntegrationScheme.h:57
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:171
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
Nektar::LibUtilities::TimeIntegrationSolution::GetTimeLevelOffset
const Array< OneD, const unsigned int > & GetTimeLevelOffset()
Definition: TimeIntegrationScheme.h:720
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::LibUtilities::TimeIntegrationScheme::m_T
NekDouble m_T
Used to store the Explicit stage derivative of IMEX schemes.
Definition: TimeIntegrationScheme.h: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:191
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
Nektar::LibUtilities::TimeIntegrationScheme::m_F_IMEX
TripleArray m_F_IMEX
Array corresponding to the stage Derivatives.
Definition: TimeIntegrationScheme.h:548
Generated on Tue Dec 8 2015 23:45:30 for Nektar++ by   doxygen 1.8.9.1