Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  {
442  m_numsteps = 1;
443  m_numstages = 2;
444 
445  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
446  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
447 
448  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
449  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
450  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
451  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
452 
453  m_A[0][1][0] = 0.5;
454  m_B[0][0][1] = 1.0;
455 
456  m_schemeType = eExplicit;
457  m_numMultiStepValues = 1;
458  m_numMultiStepDerivs = 0;
459  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
460  m_timeLevelOffset[0] = 0;
461  }
462  break;
463  case eClassicalRungeKutta4:
464  {
465  m_numsteps = 1;
466  m_numstages = 4;
467 
468  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
469  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
470 
471  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
472  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
473  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
474  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
475 
476  m_A[0][1][0] = 0.5;
477  m_A[0][2][1] = 0.5;
478  m_A[0][3][2] = 1.0;
479 
480  m_B[0][0][0] = 1.0/6.0;
481  m_B[0][0][1] = 1.0/3.0;
482  m_B[0][0][2] = 1.0/3.0;
483  m_B[0][0][3] = 1.0/6.0;
484 
485  m_schemeType = eExplicit;
486  m_numMultiStepValues = 1;
487  m_numMultiStepDerivs = 0;
488  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
489  m_timeLevelOffset[0] = 0;
490  }
491  break;
492  case eRungeKutta2_ModifiedEuler:
493  {
494  m_numsteps = 1;
495  m_numstages = 2;
496 
497  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
498  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
499 
500  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
501  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
502  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.5);
503  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 0.5);
504 
505  m_A[0][1][0] = 0.5;
506 
507  m_B[0][0][1] = 1.0;
508 
509 
510  m_schemeType = eExplicit;
511  m_numMultiStepValues = 1;
512  m_numMultiStepDerivs = 0;
513  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
514  m_timeLevelOffset[0] = 0;
515  }
516  break;
517  case eRungeKutta2_ImprovedEuler:
518  {
519  m_numsteps = 1;
520  m_numstages = 2;
521 
522  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
523  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
524 
525  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
526  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.5);
527  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
528  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
529 
530  m_A[0][1][0] = 1.0;
531 
532  m_schemeType = eExplicit;
533  m_numMultiStepValues = 1;
534  m_numMultiStepDerivs = 0;
535  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
536  m_timeLevelOffset[0] = 0;
537  }
538  break;
539  case eDIRKOrder2:
540  {
541  m_numsteps = 1;
542  m_numstages = 2;
543 
544  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
545  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
546 
547  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
548  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
549  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
550  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
551 
552  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
553 
554  m_A[0][0][0] = lambda;
555  m_A[0][1][0] = 1.0 - lambda;
556  m_A[0][1][1] = lambda;
557 
558  m_B[0][0][0] = 1.0 - lambda;
559  m_B[0][0][1] = lambda;
560 
561  m_schemeType = eDiagonallyImplicit;
562  m_numMultiStepValues = 1;
563  m_numMultiStepDerivs = 0;
564  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
565  m_timeLevelOffset[0] = 0;
566  }
567  break;
568  case eDIRKOrder3:
569  {
570  m_numsteps = 1;
571  m_numstages = 3;
572 
573  m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
574  m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
575 
576  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
577  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
578  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
579  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
580 
581  NekDouble lambda = 0.4358665215;
582 
583  m_A[0][0][0] = lambda;
584  m_A[0][1][0] = 0.5 * (1.0 - lambda);
585  m_A[0][2][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
586  m_A[0][1][1] = lambda;
587  m_A[0][2][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
588  m_A[0][2][2] = lambda;
589 
590  m_B[0][0][0] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
591  m_B[0][0][1] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
592  m_B[0][0][2] = lambda;
593 
594  m_schemeType = eDiagonallyImplicit;
595  m_numMultiStepValues = 1;
596  m_numMultiStepDerivs = 0;
597  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
598  m_timeLevelOffset[0] = 0;
599  }
600  break;
601  case eIMEXdirk_2_3_2:
602  {
603  m_numsteps = 1;
604  m_numstages = 3;
605 
606  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
607  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
608 
609  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
610  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
611  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
612  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
613  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
614  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
615 
616  NekDouble lambda = (2.0-sqrt(2.0))/2.0;
617  NekDouble delta = -2.0*sqrt(2.0)/3.0;
618 
619  m_A[0][1][1] = lambda;
620  m_A[0][2][1] = 1.0 - lambda;
621  m_A[0][2][2] = lambda;
622 
623  m_B[0][0][1] = 1.0 - lambda;
624  m_B[0][0][2] = lambda;
625 
626  m_A[1][1][0] = lambda;
627  m_A[1][2][0] = delta;
628  m_A[1][2][1] = 1.0 - delta;
629 
630  m_B[1][0][1] = 1.0 - lambda;
631  m_B[1][0][2] = lambda;
632 
633  m_schemeType = eIMEX;
634  m_numMultiStepValues = 1;
635  m_numMultiStepDerivs = 0;
636  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
637  m_timeLevelOffset[0] = 0;
638  }
639  break;
640  case eIMEXdirk_3_4_3:
641  {
642  m_numsteps = 1;
643  m_numstages = 4;
644 
645  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
646  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
647 
648  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
649  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
650  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
651  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
652  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
653  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
654 
655  NekDouble lambda = 0.4358665215;
656 
657  m_A[0][1][1] = lambda;
658  m_A[0][2][1] = 0.5 * (1.0 - lambda);
659  m_A[0][3][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
660  m_A[0][2][2] = lambda;
661  m_A[0][3][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
662  m_A[0][3][3] = lambda;
663 
664  m_B[0][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
665  m_B[0][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
666  m_B[0][0][3] = lambda;
667 
668  m_A[1][1][0] = 0.4358665215;
669  m_A[1][2][0] = 0.3212788860;
670  m_A[1][2][1] = 0.3966543747;
671  m_A[1][3][0] =-0.105858296;
672  m_A[1][3][1] = 0.5529291479;
673  m_A[1][3][2] = 0.5529291479;
674 
675  m_B[1][0][1] = 0.25 * (-6.0*lambda*lambda + 16.0*lambda - 1.0);
676  m_B[1][0][2] = 0.25 * ( 6.0*lambda*lambda - 20.0*lambda + 5.0);
677  m_B[1][0][3] = lambda;
678 
679  m_schemeType = eIMEX;
680  m_numMultiStepValues = 1;
681  m_numMultiStepDerivs = 0;
682  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
683  m_timeLevelOffset[0] = 0;
684  }
685  break;
686  case eCNAB:
687  {
688  NekDouble secondth = 1.0/2.0;
689  m_numsteps = 4;
690  m_numstages = 1;
691 
692  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
693  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
694 
695  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,secondth);
696  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
697  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
698  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
699  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,secondth);
700  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
701 
702  m_B[0][0][0] = secondth;
703  m_B[0][1][0] = 1.0;
704  m_B[1][2][0] = 1.0;
705  m_U[0][0] = 2*secondth;
706  m_U[0][2] = 3*secondth;
707  m_U[0][3] = -1*secondth;
708 
709  m_V[0][0] = 2*secondth;
710  m_V[0][1] = secondth;
711  m_V[0][2] = 3*secondth;
712  m_V[0][3] = -1*secondth;
713  m_V[3][2] = 1.0;
714 
715  m_schemeType = eIMEX;
716  m_numMultiStepValues = 1;
717  m_numMultiStepDerivs = 3;
718  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
719  m_timeLevelOffset[0] = 0;
720  m_timeLevelOffset[1] = 0;
721  m_timeLevelOffset[2] = 0;
722  m_timeLevelOffset[3] = 1;
723  }
724  break;
725  case eIMEXGear:
726  {
727  NekDouble twothirdth = 2.0/3.0;
728  m_numsteps = 2;
729  m_numstages = 1;
730 
731  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
732  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
733 
734  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
735  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
736  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
737  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
738  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,twothirdth);
739  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
740 
741  m_B[0][0][0] = twothirdth;
742  m_B[1][0][0] = twothirdth;
743  m_U[0][0] = 2*twothirdth;
744  m_U[0][1] = -0.5*twothirdth;
745 
746  m_V[0][0] = 2*twothirdth;
747  m_V[0][1] = -0.5*twothirdth;
748  m_V[1][0] = 1.0;
749 
750  m_schemeType = eIMEX;
751  m_numMultiStepValues = 2;
752  m_numMultiStepDerivs = 0;
753  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
754  m_timeLevelOffset[0] = 0;
755  m_timeLevelOffset[1] = 1;
756  }
757  break;
758  case eMCNAB:
759  {
760  NekDouble sixthx = 9.0/16.0;
761  m_numsteps = 5;
762  m_numstages = 1;
763 
764  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
765  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
766 
767  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,sixthx);
768  m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
769  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
770  m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
771  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 0.0);
772  m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
773 
774  m_B[0][0][0] = sixthx;
775  m_B[0][1][0] = 1.0;
776  m_B[1][3][0] = 1.0;
777  m_U[0][0] = 1.0;
778  m_U[0][1] = 6.0/16.0;
779  m_U[0][2] = 1.0/16.0;
780  m_U[0][3] = 1.5;
781  m_U[0][4] = -0.5;
782 
783  m_V[0][0] = 1.0;
784  m_V[0][1] = 6.0/16.0;
785  m_V[0][2] = 1.0/16.0;
786  m_V[0][3] = 1.5;
787  m_V[0][4] = -0.5;
788  m_V[2][1] = 1.0;
789  m_V[4][3] = 1.0;
790 
791  m_schemeType = eIMEX;
792  m_numMultiStepValues = 1;
793  m_numMultiStepDerivs = 4;
794  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
795  m_timeLevelOffset[0] = 0;
796  m_timeLevelOffset[1] = 0;
797  m_timeLevelOffset[2] = 1;
798  m_timeLevelOffset[3] = 0;
799  m_timeLevelOffset[4] = 1;
800  }
801  break;
802  case eIMEXdirk_2_2_2:
803  {
804  m_numsteps = 1;
805  m_numstages = 3;
806 
807  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
808  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
809 
810  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
811  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
812  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
813  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
814  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
815  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
816 
817  NekDouble glambda = 0.788675134594813;
818  NekDouble gdelta = 0.366025403784439;
819 
820  m_A[0][1][1] = glambda;
821  m_A[0][2][1] = 1.0 - glambda;
822  m_A[0][2][2] = glambda;
823 
824  m_B[0][0][1] = 1.0 - glambda;
825  m_B[0][0][2] = glambda;
826 
827  m_A[1][1][0] = glambda;
828  m_A[1][2][0] = gdelta;
829  m_A[1][2][1] = 1.0 - gdelta;
830 
831  m_B[1][0][0] = gdelta;
832  m_B[1][0][1] = 1.0 - gdelta;
833 
834  m_schemeType = eIMEX;
835  m_numMultiStepValues = 1;
836  m_numMultiStepDerivs = 0;
837  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
838  m_timeLevelOffset[0] = 0;
839  }
840  break;
841  case eIMEXdirk_2_3_3:
842  {
843  m_numsteps = 1;
844  m_numstages = 3;
845 
846  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
847  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
848 
849  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
850  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
851  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
852  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
853  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
854  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
855 
856  NekDouble glambda = 0.788675134594813;
857 
858  m_A[0][1][1] = glambda;
859  m_A[0][2][1] = 1.0 - 2.0*glambda;
860  m_A[0][2][2] = glambda;
861 
862  m_B[0][0][1] = 0.5;
863  m_B[0][0][2] = 0.5;
864 
865  m_A[1][1][0] = glambda;
866  m_A[1][2][0] = glambda - 1.0;
867  m_A[1][2][1] = 2.0*(1-glambda);
868 
869  m_B[1][0][1] = 0.5;
870  m_B[1][0][2] = 0.5;
871 
872  m_schemeType = eIMEX;
873  m_numMultiStepValues = 1;
874  m_numMultiStepDerivs = 0;
875  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
876  m_timeLevelOffset[0] = 0;
877  }
878  break;
879  case eIMEXdirk_1_1_1:
880  {
881  m_numsteps = 1;
882  m_numstages = 2;
883 
884  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
885  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
886 
887  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
888  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
889  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
890  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
891  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
892  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
893 
894  m_A[0][1][1] = 1.0;
895 
896  m_B[0][0][1] = 1.0;
897 
898  m_A[1][1][0] = 1.0;
899 
900  m_B[1][0][0] = 1.0;
901 
902  m_schemeType = eIMEX;
903  m_numMultiStepValues = 1;
904  m_numMultiStepDerivs = 0;
905  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
906  m_timeLevelOffset[0] = 0;
907  }
908  break;
909  case eIMEXdirk_1_2_1:
910  {
911  m_numsteps = 1;
912  m_numstages = 2;
913 
914  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
915  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
916 
917  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
918  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
919  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
920  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
921  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
922  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
923 
924 
925  m_A[0][1][1] = 1.0;
926 
927  m_B[0][0][1] = 1.0;
928 
929  m_A[1][1][0] = 1.0;
930 
931  m_B[1][0][1] = 1.0;
932 
933  m_schemeType = eIMEX;
934  m_numMultiStepValues = 1;
935  m_numMultiStepDerivs = 0;
936  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
937  m_timeLevelOffset[0] = 0;
938  }
939  break;
940  case eIMEXdirk_1_2_2:
941  {
942  m_numsteps = 1;
943  m_numstages = 2;
944 
945  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
946  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
947 
948  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
949  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
950  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
951  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
952  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
953  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
954 
955  m_A[0][1][1] = 1.0/2.0;
956 
957  m_B[0][0][1] = 1.0;
958 
959  m_A[1][1][0] = 1.0/2.0;
960 
961  m_B[1][0][1] = 1.0;
962 
963  m_schemeType = eIMEX;
964  m_numMultiStepValues = 1;
965  m_numMultiStepDerivs = 0;
966  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
967  m_timeLevelOffset[0] = 0;
968  }
969  break;
970  case eIMEXdirk_4_4_3:
971  {
972  m_numsteps = 1;
973  m_numstages = 5;
974 
975  m_A = Array<OneD, Array<TwoD,NekDouble> >(2);
976  m_B = Array<OneD, Array<TwoD,NekDouble> >(2);
977 
978  m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
979  m_B[0] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
980  m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
981  m_B[1] = Array<TwoD,NekDouble>(m_numsteps, m_numstages,0.0);
982  m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps, 1.0);
983  m_V = Array<TwoD,NekDouble>(m_numsteps, m_numsteps, 1.0);
984 
985  m_A[0][1][1] = 1.0/2.0;
986  m_A[0][2][1] = 1.0/6.0;
987  m_A[0][2][2] = 1.0/2.0;
988  m_A[0][3][1] = -1.0/2.0;
989  m_A[0][3][2] = 1.0/2.0;
990  m_A[0][3][3] = 1.0/2.0;
991  m_A[0][4][1] = 3.0/2.0;
992  m_A[0][4][2] = -3.0/2.0;
993  m_A[0][4][3] = 1.0/2.0;
994  m_A[0][4][4] = 1.0/2.0;
995 
996 
997  m_B[0][0][1] = 3.0/2.0;
998  m_B[0][0][2] = -3.0/2.0;
999  m_B[0][0][3] = 1.0/2.0;
1000  m_B[0][0][4] = 1.0/2.0;
1001 
1002  m_A[1][1][0] = 1.0/2.0;
1003  m_A[1][2][0] = 11.0/18.0;
1004  m_A[1][2][1] = 1.0/18.0;
1005  m_A[1][3][0] = 5.0/6.0;
1006  m_A[1][3][1] = -5.0/6.0;
1007  m_A[1][3][2] = 1.0/2.0;
1008  m_A[1][4][0] = 1.0/4.0;
1009  m_A[1][4][1] = 7.0/4.0;
1010  m_A[1][4][2] = 3.0/4.0;
1011  m_A[1][4][3] = -7.0/4.0;
1012 
1013  m_B[1][0][0] = 1.0/4.0;
1014  m_B[1][0][1] = 7.0/4.0;
1015  m_B[1][0][2] = 3.0/4.0;
1016  m_B[1][0][3] = -7.0/4.0;
1017 
1018  m_schemeType = eIMEX;
1019  m_numMultiStepValues = 1;
1020  m_numMultiStepDerivs = 0;
1021  m_timeLevelOffset = Array<OneD,unsigned int>(m_numsteps);
1022  m_timeLevelOffset[0] = 0;
1023  }
1024  break;
1025  default:
1026  {
1027  NEKERROR(ErrorUtil::efatal,"Invalid Time Integration Scheme");
1028  }
1029  }
1030 
1031  m_firstStageEqualsOldSolution = CheckIfFirstStageEqualsOldSolution(m_A,m_B,m_U,m_V);
1032  m_lastStageEqualsNewSolution = CheckIfLastStageEqualsNewSolution(m_A,m_B,m_U,m_V);
1033 
1034  ASSERTL1(VerifyIntegrationSchemeType(m_schemeType,m_A,m_B,m_U,m_V),
1035  "Time integration scheme coefficients do not match its type");
1036  }
1037 
1038 
1039  bool TimeIntegrationScheme::
1040  VerifyIntegrationSchemeType(TimeIntegrationSchemeType type,
1041  const Array<OneD, const Array<TwoD, NekDouble> >& A,
1042  const Array<OneD, const Array<TwoD, NekDouble> >& B,
1043  const Array<TwoD, const NekDouble>& U,
1044  const Array<TwoD, const NekDouble>& V) const
1045  {
1046  int i;
1047  int j;
1048  int m;
1049  int IMEXdim = A.num_elements();
1050  int dim = A[0].GetRows();
1051 
1052  Array<OneD, TimeIntegrationSchemeType> vertype(IMEXdim,eExplicit);
1053 
1054  for(m = 0; m < IMEXdim; m++)
1055  {
1056  for(i = 0; i < dim; i++)
1057  {
1058  if( fabs(A[m][i][i]) > NekConstants::kNekZeroTol )
1059  {
1060  vertype[m] = eDiagonallyImplicit;
1061  }
1062  }
1063 
1064  for(i = 0; i < dim; i++)
1065  {
1066  for(j = i+1; j < dim; j++)
1067  {
1068  if( fabs(A[m][i][j]) > NekConstants::kNekZeroTol )
1069  {
1070  vertype[m] = eImplicit;
1071  ASSERTL1(false,"Fully Impplicit schemes cannnot be handled by the TimeIntegrationScheme class");
1072  }
1073  }
1074  }
1075  }
1076 
1077  if(IMEXdim == 2)
1078  {
1079  ASSERTL1(B.num_elements()==2,"Coefficient Matrix B should have an implicit and explicit part for IMEX schemes");
1080  if((vertype[0] == eDiagonallyImplicit) &&
1081  (vertype[1] == eExplicit))
1082  {
1083  vertype[0] = eIMEX;
1084  }
1085  else
1086  {
1087  ASSERTL1(false,"This is not a proper IMEX scheme");
1088  }
1089  }
1090 
1091  return (vertype[0] == type);
1092  }
1093 
1094  TimeIntegrationSolutionSharedPtr
1095  TimeIntegrationScheme::InitializeScheme(const NekDouble timestep,
1096  ConstDoubleArray &y_0 ,
1097  const NekDouble time ,
1098  const TimeIntegrationSchemeOperators &op)
1099  {
1100  // create a TimeIntegrationSolution object based upon the
1101  // initial value. Initialise all other multi-step values
1102  // and derivatives to zero
1103  TimeIntegrationSolutionSharedPtr y_out =
1104  MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(m_schemeKey,y_0,time,timestep);
1105 
1106  // calculate the initial derivative, if is part of the
1107  // solution vector of the current scheme
1108  if(m_numMultiStepDerivs)
1109  {
1110  if(m_timeLevelOffset[m_numMultiStepValues] == 0)
1111  {
1112  int i;
1113  int nvar = y_0.num_elements();
1114  int npoints = y_0[0].num_elements();
1115  DoubleArray f_y_0(nvar);
1116  for(i = 0; i < nvar; i++)
1117  {
1118  f_y_0[i] = Array<OneD,NekDouble>(npoints);
1119  }
1120  // calculate the derivative of the initial value
1121  op.DoOdeRhs(y_0,f_y_0,time);
1122 
1123  // multiply by the step size
1124  for(i = 0; i < nvar; i++)
1125  {
1126  Blas::Dscal(npoints,timestep,f_y_0[i].get(),1);
1127  }
1128  y_out->SetDerivative(0,f_y_0,timestep);
1129  }
1130  }
1131 
1132  return y_out;
1133  }
1134 
1135  TimeIntegrationScheme::ConstDoubleArray&
1136  TimeIntegrationScheme::TimeIntegrate(const NekDouble timestep,
1137  TimeIntegrationSolutionSharedPtr &solvector,
1138  const TimeIntegrationSchemeOperators &op)
1139  {
1140  ASSERTL1(!(GetIntegrationSchemeType() == eImplicit),
1141  "Fully Implicit integration scheme cannot be handled by this routine.");
1142 
1143  int nvar = solvector->GetFirstDim ();
1144  int npoints = solvector->GetSecondDim();
1145 
1146  if( (solvector->GetIntegrationScheme()).get() != this )
1147  {
1148  // This branch will be taken when the solution vector
1149  // (solvector) is set up for a different scheme than
1150  // the object this method is called from. (typically
1151  // needed to calculate the first time-levels of a
1152  // multi-step scheme)
1153 
1154  // To do this kind of 'non-matching' integration, we
1155  // perform the following three steps:
1156  //
1157  // 1: copy the required input information from the
1158  // solution vector of the master scheme to the
1159  // input solution vector of the current scheme
1160  //
1161  // 2: time-integrate for one step using the current
1162  // scheme
1163  //
1164  // 3: copy the information contained in the output
1165  // vector of the current scheme to the solution
1166  // vector of the master scheme
1167 
1168  // STEP 1: copy the required input information from
1169  // the solution vector of the master scheme
1170  // to the input solution vector of the
1171  // current scheme
1172 
1173  // 1.1 Determine which information is required for the
1174  // current scheme
1175  int n;
1176  DoubleArray y_n;
1177  NekDouble t_n = 0;
1178  DoubleArray dtFy_n;
1179  unsigned int nCurSchemeVals = m_numMultiStepValues; // number of required values of the current scheme
1180  unsigned int nCurSchemeDers = m_numMultiStepDerivs; // number of required derivs of the current scheme
1181  unsigned int nCurSchemeSteps = m_numsteps; // number of steps in the current scheme
1182  unsigned int nMasterSchemeVals = solvector->GetNvalues(); // number of values of the master scheme
1183  unsigned int nMasterSchemeDers = solvector->GetNderivs(); // number of derivs of the master scheme
1184  // The arrays below contains information to which
1185  // time-level the values and derivatives of the
1186  // schemes belong
1187  const Array<OneD, const unsigned int>& curTimeLevels = m_timeLevelOffset;
1188  const Array<OneD, const unsigned int>& masterTimeLevels = solvector->GetTimeLevelOffset();
1189 
1190  // 1.2 Copy the required information from the master
1191  // solution vector to the input solution vector of
1192  // the current scheme
1193  TimeIntegrationSolutionSharedPtr solvector_in = MemoryManager<TimeIntegrationSolution>::
1194  AllocateSharedPtr(GetIntegrationSchemeKey()); // input solution vector of the current scheme
1195 
1196  for(n = 0; n < nCurSchemeVals; n++)
1197  {
1198  // Get the required value out of the master solution vector
1199  //DoubleArray& y_n = solvector->GetValue ( curTimeLevels[n] );
1200  //NekDouble t_n = solvector->GetValueTime( curTimeLevels[n] );
1201 
1202  y_n = solvector->GetValue ( curTimeLevels[n] );
1203  t_n = solvector->GetValueTime( curTimeLevels[n] );
1204 
1205  // Set the required value in the input solution
1206  // vector of the current scheme
1207  solvector_in->SetValue(curTimeLevels[n],y_n,t_n);
1208  }
1209  for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
1210  {
1211  // Get the required derivative out of the master
1212  // solution vector
1213  //DoubleArray& dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
1214  dtFy_n = solvector->GetDerivative ( curTimeLevels[n] );
1215 
1216  // Set the required derivative in the input
1217  // solution vector of the current scheme
1218  solvector_in->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
1219  }
1220 
1221  // STEP 2: time-integrate for one step using the
1222  // current scheme
1223  TimeIntegrationSolutionSharedPtr solvector_out = MemoryManager<TimeIntegrationSolution>:: AllocateSharedPtr(GetIntegrationSchemeKey(),nvar,npoints); // output solution vector of the current scheme
1224 
1225  // integrate
1226  TimeIntegrate(timestep, solvector_in->GetSolutionVector(),
1227  solvector_in->GetTimeVector(),
1228  solvector_out->UpdateSolutionVector(),
1229  solvector_out->UpdateTimeVector(),op);
1230 
1231 
1232  // STEP 3: copy the information contained in the
1233  // output vector of the current scheme to the
1234  // solution vector of the master scheme
1235 
1236  // 3.1 Check whether the current time scheme updates
1237  // the most recent derivative that should be
1238  // updated in the master scheme. If not,
1239  // calculate the derivative. This can be done
1240  // based upon the corresponding value and the
1241  // DoOdeRhs operator.
1242  int j;
1243  bool CalcNewDeriv = false; // flag inidicating whether the new derivative is availble in the output of
1244  // of the current scheme or whether it should be calculated
1245  if( nMasterSchemeDers > 0 )
1246  {
1247  if(nCurSchemeDers == 0)
1248  {
1249  CalcNewDeriv = true;
1250  }
1251  else
1252  {
1253  if( masterTimeLevels[nMasterSchemeVals] < curTimeLevels[nCurSchemeVals] )
1254  {
1255  CalcNewDeriv = true;
1256  }
1257  }
1258  }
1259 
1260  if(CalcNewDeriv)
1261  {
1262  int newDerivTimeLevel = masterTimeLevels[nMasterSchemeVals]; // contains the time level at which
1263  // we want to know the derivative of the
1264  // master scheme
1265  //DoubleArray y_n;
1266  //NekDouble t_n;
1267  // if the time level correspond to 0, calculate the derivative based upon the solution value
1268  // at the new time-level
1269  if (newDerivTimeLevel == 0)
1270  {
1271  y_n = solvector_out->GetValue(0);
1272  t_n = solvector_out->GetValueTime(0);
1273  }
1274  // if the time level correspond to 1, calculate the derivative based upon the solution value
1275  // at the new old-level
1276  else if( newDerivTimeLevel == 1 )
1277  {
1278  y_n = solvector->GetValue(0);
1279  t_n = solvector->GetValueTime(0);
1280  }
1281  else
1282  {
1283  ASSERTL1(false,"Problems with initialising scheme");
1284  }
1285 
1286  DoubleArray f_n(nvar);
1287  for(j = 0; j < nvar; j++)
1288  {
1289  f_n[j] = Array<OneD, NekDouble>(npoints);
1290  }
1291 
1292  // calculate the derivative
1293  op.DoOdeRhs(y_n, f_n, t_n);
1294 
1295  // multiply by dt (as required by the General Linear Method framework)
1296  for(j = 0; j < nvar; j++)
1297  {
1298  Vmath::Smul(npoints,timestep,f_n[j],1,
1299  f_n[j],1);
1300  }
1301 
1302  // Rotate the solution vector
1303  // (i.e. updating without calculating/inserting new values)
1304  solvector->RotateSolutionVector();
1305  // Set the calculated derivative in the master solution vector
1306  solvector->SetDerivative(newDerivTimeLevel,f_n,timestep);
1307  }
1308  else
1309  {
1310  // Rotate the solution vector (i.e. updating
1311  // without calculating/inserting new values)
1312  solvector->RotateSolutionVector();
1313  }
1314 
1315 
1316  // 1.2 Copy the information calculated using the
1317  // current scheme from the output solution vector
1318  // to the master solution vector
1319  for(n = 0; n < nCurSchemeVals; n++)
1320  {
1321  // Get the calculated value out of the output
1322  // solution vector of the current scheme
1323  //DoubleArray& y_n = solvector_out->GetValue ( curTimeLevels[n] );
1324  //NekDouble t_n = solvector_out->GetValueTime( curTimeLevels[n] );
1325  y_n = solvector_out->GetValue ( curTimeLevels[n] );
1326  t_n = solvector_out->GetValueTime( curTimeLevels[n] );
1327 
1328  // Set the calculated value in the master solution vector
1329  solvector->SetValue(curTimeLevels[n],y_n,t_n);
1330  }
1331 
1332  for(n = nCurSchemeVals; n < nCurSchemeSteps; n++)
1333  {
1334  // Get the calculated derivative out of the output
1335  // solution vector of the current scheme
1336  // DoubleArray& dtFy_n =
1337  // solvector_out->GetDerivative (curTimeLevels[n]);
1338  dtFy_n = solvector_out->GetDerivative ( curTimeLevels[n] );
1339 
1340  // Set the calculated derivative in the master
1341  // solution vector
1342  solvector->SetDerivative(curTimeLevels[n],dtFy_n,timestep);
1343  }
1344  }
1345  else
1346  {
1347  const TimeIntegrationSchemeKey& key = solvector->GetIntegrationSchemeKey();
1348 
1349  TimeIntegrationSolutionSharedPtr solvector_new = MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(key,nvar,npoints);
1350 
1351  TimeIntegrate(timestep,solvector->GetSolutionVector(),
1352  solvector->GetTimeVector(),
1353  solvector_new->UpdateSolutionVector(),
1354  solvector_new->UpdateTimeVector(),op);
1355 
1356  solvector = solvector_new;
1357  }
1358  return solvector->GetSolution();
1359  }
1360 
1361  void TimeIntegrationScheme::TimeIntegrate(const NekDouble timestep,
1362  ConstTripleArray &y_old ,
1363  ConstSingleArray &t_old ,
1364  TripleArray &y_new ,
1365  SingleArray &t_new ,
1366  const TimeIntegrationSchemeOperators &op)
1367  {
1368  ASSERTL1(CheckTimeIntegrateArguments(timestep,y_old,t_old,y_new,t_new,op), "Arguments not well defined");
1369 
1370  unsigned int i,j,k;
1371  TimeIntegrationSchemeType type = GetIntegrationSchemeType();
1372 
1373  // Check if storage has already been initialised.
1374  // If so, we just zero the temporary storage.
1375  if (m_initialised && m_nvar == GetFirstDim(y_old)
1376  && m_npoints == GetSecondDim(y_old))
1377  {
1378  for(j = 0; j < m_nvar; j++)
1379  {
1380  Vmath::Zero(m_npoints, m_tmp[j], 1);
1381  }
1382  }
1383  else
1384  {
1385  m_nvar = GetFirstDim(y_old);
1386  m_npoints = GetSecondDim(y_old);
1387 
1388  // First, we are going to calculate the various stage
1389  // values and stage derivatives (this is the multi-stage
1390  // part of the method)
1391  // - m_Y corresponds to the stage values
1392  // - m_F corresponds to the stage derivatives
1393  // - m_T corresponds to the time at the different stages
1394  // - m_tmp corresponds to the explicit right hand side of
1395  // each stage equation
1396  // (for explicit schemes, this correspond to m_Y)
1397 
1398  // Allocate memory for the arrays m_Y and m_F and m_tmp The same
1399  // storage will be used for every stage -> m_Y is a
1400  // DoubleArray
1401  m_tmp = DoubleArray(m_nvar);
1402  for(j = 0; j < m_nvar; j++)
1403  {
1404  m_tmp[j] = Array<OneD, NekDouble>(m_npoints,0.0);
1405  }
1406 
1407  // The same storage will be used for every stage -> m_tmp is
1408  // a DoubleArray
1409  if(type == eExplicit)
1410  {
1411  m_Y = m_tmp;
1412  }
1413  else
1414  {
1415  m_Y = DoubleArray(m_nvar);
1416  for(j = 0; j < m_nvar; j++)
1417  {
1418  m_Y[j] = Array<OneD, NekDouble>(m_npoints,0.0);
1419  }
1420  }
1421 
1422  // Different storage for every stage derivative as the data
1423  // will be re-used to update the solution -> m_F is a TripleArray
1424  m_F = TripleArray(m_numstages);
1425  for(i = 0; i < m_numstages; ++i)
1426  {
1427  m_F[i] = DoubleArray(m_nvar);
1428  for(j = 0; j < m_nvar; j++)
1429  {
1430  m_F[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
1431  }
1432  }
1433 
1434  if(type == eIMEX)
1435  {
1436  m_F_IMEX = TripleArray(m_numstages);
1437  for(i = 0; i < m_numstages; ++i)
1438  {
1439  m_F_IMEX[i] = DoubleArray(m_nvar);
1440  for(j = 0; j < m_nvar; j++)
1441  {
1442  m_F_IMEX[i][j] = Array<OneD, NekDouble>(m_npoints,0.0);
1443  }
1444  }
1445  }
1446 
1447  // Finally, flag that we have initialised the memory.
1448  m_initialised = true;
1449  }
1450 
1451  // The loop below calculates the stage values and derivatives
1452  for(i = 0; i < m_numstages; i++)
1453  {
1454  if( (i==0) && m_firstStageEqualsOldSolution )
1455  {
1456  for(k = 0; k < m_nvar; k++)
1457  {
1458  Vmath::Vcopy(m_npoints,y_old[0][k],1,m_Y[k],1);
1459  }
1460  m_T = t_old[0];
1461  }
1462  else
1463  {
1464  // The stage values m_Y are a linear combination of:
1465  // 1: the stage derivatives
1466 
1467  if( i != 0 )
1468  {
1469  for(k = 0; k < m_nvar; k++)
1470  {
1471  Vmath::Smul(m_npoints,timestep*A(i,0),m_F[0][k],1,
1472  m_tmp[k],1);
1473 
1474  if(type == eIMEX)
1475  {
1476  Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,0),
1477  m_F_IMEX[0][k],1,
1478  m_tmp[k],1,m_tmp[k],1);
1479  }
1480  }
1481  }
1482  m_T = A(i,0)*timestep;
1483 
1484  for( j = 1; j < i; j++ )
1485  {
1486  for(k = 0; k < m_nvar; k++)
1487  {
1488  Vmath::Svtvp(m_npoints,timestep*A(i,j),m_F[j][k],1,
1489  m_tmp[k],1,m_tmp[k],1);
1490  if(type == eIMEX)
1491  {
1492  Vmath::Svtvp(m_npoints,timestep*A_IMEX(i,j),
1493  m_F_IMEX[j][k],1,
1494  m_tmp[k],1,m_tmp[k],1);
1495  }
1496  }
1497 
1498  m_T += A(i,j)*timestep;
1499  }
1500 
1501  // 2: the imported multi-step solution of the
1502  // previous time level
1503  for(j = 0; j < m_numsteps; j++)
1504  {
1505  for(k = 0; k < m_nvar; k++)
1506  {
1507  Vmath::Svtvp(m_npoints,U(i,j),y_old[j][k],1,
1508  m_tmp[k],1,m_tmp[k],1);
1509  }
1510  m_T += U(i,j)*t_old[j];
1511  }
1512  }
1513 
1514  // Calculate the stage derivative based upon the stage value
1515  if(type == eDiagonallyImplicit)
1516  {
1517  if(m_numstages==1)
1518  {
1519  m_T= t_old[0]+timestep;
1520  }
1521  else
1522  {
1523  m_T= t_old[0];
1524  for(int j=0; j<=i; ++j)
1525  {
1526  m_T += A(i,j)*timestep;
1527  }
1528  }
1529 
1530  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
1531 
1532  for(k = 0; k < m_nvar; k++)
1533  {
1534  Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
1535  Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),m_F[i][k],1,m_F[i][k],1);
1536  }
1537  }
1538  else if(type == eIMEX)
1539  {
1540  if(m_numstages==1)
1541  {
1542  m_T= t_old[0]+timestep;
1543  }
1544  else
1545  {
1546  m_T= t_old[0];
1547  for(int j=0; j<=i; ++j)
1548  {
1549  m_T += A(i,j)*timestep;
1550  }
1551  }
1552 
1553  if(fabs(A(i,i)) > NekConstants::kNekZeroTol)
1554  {
1555  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(i,i)*timestep);
1556 
1557  for(k = 0; k < m_nvar; k++)
1558  {
1559  Vmath::Vsub(m_npoints,m_Y[k],1,m_tmp[k],1,m_F[i][k],1);
1560  Vmath::Smul(m_npoints,1.0/(A(i,i)*timestep),
1561  m_F[i][k],1,m_F[i][k],1);
1562  }
1563  }
1564  op.DoOdeRhs(m_Y, m_F_IMEX[i], m_T);
1565  }
1566  else if( type == eExplicit)
1567  {
1568  // ensure solution is in correct space
1569  op.DoProjection(m_Y,m_Y,m_T);
1570  op.DoOdeRhs(m_Y, m_F[i], m_T);
1571  }
1572  }
1573 
1574  // Next, the solution vector y at the new time level will
1575  // be calculated.
1576  //
1577  // For multi-step methods, this includes updating the
1578  // values of the auxiliary parameters
1579  //
1580  // The loop below calculates the solution at the new time
1581  // level
1582  //
1583  // If last stage equals the new solution, the new solution
1584  // needs not be calculated explicitly but can simply be
1585  // copied. This saves a solve.
1586  int i_start = 0;
1587  if( (m_lastStageEqualsNewSolution) )
1588  {
1589  for(k = 0; k < m_nvar; k++)
1590  {
1591  Vmath::Vcopy(m_npoints,m_Y[k],1,y_new[0][k],1);
1592  }
1593 
1594  if (m_numstages==1 && type == eIMEX)
1595  {
1596  t_new[0] = t_old[0]+timestep;
1597  }
1598  else
1599  {
1600  t_new[0] = B(0,0)*timestep;
1601  for(j = 1; j < m_numstages; j++)
1602  {
1603  t_new[0] += B(0,j)*timestep;
1604  }
1605  for(j = 0; j < m_numsteps; j++)
1606  {
1607  t_new[0] += V(0,j)*t_old[j];
1608  }
1609  i_start = 1;
1610  }
1611  }
1612 
1613  for(i = i_start; i < m_numsteps; i++)
1614  {
1615  // The solution at the new time level is a linear
1616  // combination of:
1617  // 1: the stage derivatives
1618  for(k = 0; k < m_nvar; k++)
1619  {
1620  Vmath::Smul(m_npoints,timestep*B(i,0),m_F[0][k],1,
1621  y_new[i][k],1);
1622 
1623  if(type == eIMEX)
1624  {
1625  Vmath::Svtvp(m_npoints,timestep*B_IMEX(i,0),
1626  m_F_IMEX[0][k],1,y_new[i][k],1,
1627  y_new[i][k],1);
1628  }
1629  }
1630  if(m_numstages != 1 || type != eIMEX)
1631  {
1632  t_new[i] = B(i,0)*timestep;
1633  }
1634 
1635 
1636  for(j = 1; j < m_numstages; j++)
1637  {
1638  for(k = 0; k < m_nvar; k++)
1639  {
1640  Vmath::Svtvp(m_npoints,timestep*B(i,j),m_F[j][k],1,
1641  y_new[i][k],1,y_new[i][k],1);
1642 
1643  if(type == eIMEX)
1644  {
1645  Vmath::Svtvp(m_npoints,timestep*B_IMEX(i,j),
1646  m_F_IMEX[j][k],1,y_new[i][k],1,
1647  y_new[i][k],1);
1648  }
1649  }
1650  if(m_numstages != 1 || type != eIMEX)
1651  {
1652  t_new[i] += B(i,j)*timestep;
1653  }
1654  }
1655 
1656  // 2: the imported multi-step solution of the previous
1657  // time level
1658  for(j = 0; j < m_numsteps; j++)
1659  {
1660  for(k = 0; k < m_nvar; k++)
1661  {
1662  Vmath::Svtvp(m_npoints,V(i,j),y_old[j][k],1,
1663  y_new[i][k],1,y_new[i][k],1);
1664  }
1665  if(m_numstages != 1 || type != eIMEX)
1666  {
1667  t_new[i] += V(i,j)*t_old[j];
1668  }
1669  }
1670  }
1671 
1672  // Ensure that the new solution is projected if necessary
1673  if(type == eExplicit)
1674  {
1675  op.DoProjection(y_new[0],y_new[0],t_new[0]);
1676  }
1677  }
1678 
1679  bool TimeIntegrationScheme::CheckIfFirstStageEqualsOldSolution(const Array<OneD, const Array<TwoD, NekDouble> >& A,
1680  const Array<OneD, const Array<TwoD, NekDouble> >& B,
1681  const Array<TwoD, const NekDouble>& U,
1682  const Array<TwoD, const NekDouble>& V) const
1683  {
1684  int i,m;
1685  // First stage equals old solution if:
1686  // 1. the first row of the coefficient matrix A consists of zeros
1687  // 2. U[0][0] is equal to one and all other first row entries of U are zero
1688 
1689  // 1. Check first condition
1690  for(m = 0; m < A.num_elements(); m++)
1691  {
1692  for(i = 0; i < m_numstages; i++)
1693  {
1694  if( fabs(A[m][0][i]) > NekConstants::kNekZeroTol )
1695  {
1696  return false;
1697  }
1698  }
1699  }
1700 
1701  // 2. Check second condition
1702  if( fabs(U[0][0] - 1.0) > NekConstants::kNekZeroTol )
1703  {
1704  return false;
1705  }
1706  for(i = 1; i < m_numsteps; i++)
1707  {
1708  if( fabs(U[0][i]) > NekConstants::kNekZeroTol )
1709  {
1710  return false;
1711  }
1712  }
1713 
1714  return true;
1715  }
1716 
1717  bool TimeIntegrationScheme::CheckIfLastStageEqualsNewSolution(const Array<OneD, const Array<TwoD, NekDouble> >& A,
1718  const Array<OneD, const Array<TwoD, NekDouble> >& B,
1719  const Array<TwoD, const NekDouble>& U,
1720  const Array<TwoD, const NekDouble>& V) const
1721  {
1722  int i,m;
1723  // Last stage equals new solution if:
1724  // 1. the last row of the coefficient matrix A is equal to the first row of matrix B
1725  // 2. the last row of the coefficient matrix U is equal to the first row of matrix V
1726 
1727  // 1. Check first condition
1728  for(m = 0; m < A.num_elements(); m++)
1729  {
1730  for(i = 0; i < m_numstages; i++)
1731  {
1732  if( fabs(A[m][m_numstages-1][i]-B[m][0][i]) > NekConstants::kNekZeroTol )
1733  {
1734  return false;
1735  }
1736  }
1737  }
1738 
1739  // 2. Check second condition
1740  for(i = 0; i < m_numsteps; i++)
1741  {
1742  if( fabs(U[m_numstages-1][i]-V[0][i]) > NekConstants::kNekZeroTol )
1743  {
1744  return false;
1745  }
1746  }
1747 
1748  return true;
1749  }
1750 
1751  bool TimeIntegrationScheme::CheckTimeIntegrateArguments(const NekDouble timestep,
1752  ConstTripleArray &y_old ,
1753  ConstSingleArray &t_old ,
1754  TripleArray &y_new ,
1755  SingleArray &t_new ,
1756  const TimeIntegrationSchemeOperators &op) const
1757  {
1758  // Check if arrays are all of consistent size
1759  ASSERTL1(y_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1760  ASSERTL1(y_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1761 
1762  ASSERTL1(y_old[0]. num_elements()==y_new[0]. num_elements(),"Non-matching number of variables.");
1763  ASSERTL1(y_old[0][0].num_elements()==y_new[0][0].num_elements(),"Non-matching number of coefficients.");
1764 
1765  ASSERTL1(t_old.num_elements()==m_numsteps,"Non-matching number of steps.");
1766  ASSERTL1(t_new.num_elements()==m_numsteps,"Non-matching number of steps.");
1767 
1768  return true;
1769  }
1770 
1771  std::ostream& operator<<(std::ostream& os, const TimeIntegrationSchemeSharedPtr& rhs)
1772  {
1773  return operator<<(os,*rhs);
1774  }
1775 
1776  std::ostream& operator<<(std::ostream& os, const TimeIntegrationScheme& rhs)
1777  {
1778  int i,j;
1779  int r = rhs.GetNsteps();
1780  int s = rhs.GetNstages();
1781  TimeIntegrationSchemeType type = rhs.GetIntegrationSchemeType();
1782 
1783  int oswidth = 9;
1784  int osprecision = 6;
1785 
1786  os << "Time Integration Scheme: " << TimeIntegrationMethodMap[rhs.GetIntegrationMethod()] << endl;
1787  os << "- number of steps: " << r << endl;
1788  os << "- number of stages: " << s << endl;
1789  os << "- type of scheme: " << TimeIntegrationSchemeTypeMap[rhs.GetIntegrationSchemeType()] << endl;
1790  os << "General linear method tableau: " << endl;
1791 
1792  for(i = 0; i < s; i++)
1793  {
1794  for(j = 0; j < s; j++)
1795  {
1796  os.width(oswidth);
1797  os.precision(osprecision);
1798  os << right << rhs.A(i,j) << " ";
1799  }
1800  if(type == eIMEX)
1801  {
1802  os << " '";
1803  for(j = 0; j < s; j++)
1804  {
1805  os.width(oswidth);
1806  os.precision(osprecision);
1807  os << right << rhs.A_IMEX(i,j) << " ";
1808  }
1809  }
1810  os << " |";
1811 
1812  for(j = 0; j < r; j++)
1813  {
1814  os.width(oswidth);
1815  os.precision(osprecision);
1816  os << right << rhs.U(i,j);
1817  }
1818  os << endl;
1819  }
1820  int imexflag = (type == eIMEX)?2:1;
1821  for(int i = 0; i < (r+imexflag*s)*(oswidth+1)+imexflag*2-1; i++)
1822  {
1823  os << "-";
1824  }
1825  os << endl;
1826  for(i = 0; i < r; i++)
1827  {
1828  for(j = 0; j < s; j++)
1829  {
1830  os.width(oswidth);
1831  os.precision(osprecision);
1832  os << right << rhs.B(i,j) << " ";
1833  }
1834  if(type == eIMEX)
1835  {
1836  os << " '";
1837  for(j = 0; j < s; j++)
1838  {
1839  os.width(oswidth);
1840  os.precision(osprecision);
1841  os << right << rhs.B_IMEX(i,j) << " ";
1842  }
1843  }
1844  os << " |";
1845 
1846  for(j = 0; j < r; j++)
1847  {
1848  os.width(oswidth);
1849  os.precision(osprecision);
1850  os << right << rhs.V(i,j);
1851  }
1852  os << endl;
1853  }
1854  return os;
1855  }
1856  }
1857 }
1858 
Generated on Wed Mar 4 2015 16:33:18 for Nektar++ by   doxygen 1.8.1.2