Nektar++
TimeIntegrationAlgorithmGLM.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File TimeIntegrationAlgorithmGLM.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 scheme data class.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 
43 
44 #include <iostream>
45 
46 #include <boost/core/ignore_unused.hpp>
47 
48 #include <cmath>
49 
50 namespace Nektar
51 {
52 namespace LibUtilities
53 {
54 
55 ////////////////////////////////////////////////////////////////////////////////
57  const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time,
59 {
60  // create a TimeIntegrationSolutionGLM object based upon the
61  // initial value. Initialise all other multi-step values
62  // and derivatives to zero
63 
66  this, y_0, time, deltaT);
67 
69  {
70  // ensure initial solution is in correct space
71  op.DoProjection(y_0, y_out->UpdateSolution(), time);
72  }
73 
74  // calculate the initial derivative, if is part of the
75  // solution vector of the current scheme
76  if (GetNmultiStepDerivs() > 0)
77  {
79 
80  if (offsets[GetNmultiStepValues()] == 0)
81  {
82  int i;
83  int nvar = y_0.size();
84  int npoints = y_0[0].size();
85  DoubleArray f_y_0(nvar);
86  for (i = 0; i < nvar; i++)
87  {
88  f_y_0[i] = Array<OneD, NekDouble>(npoints);
89  }
90  // calculate the derivative of the initial value
91  op.DoOdeRhs(y_0, f_y_0, time);
92 
93  // multiply by the step size
94  for (i = 0; i < nvar; i++)
95  {
96  Blas::Dscal(npoints, deltaT, f_y_0[i].get(), 1);
97  }
98  y_out->SetDerivative(0, f_y_0, deltaT);
99  }
100  }
101 
102  return y_out;
103 }
104 
106  const NekDouble deltaT, TimeIntegrationSolutionGLMSharedPtr &solvector,
108 {
109  // ASSERTL1( !(m_parent->GetIntegrationSchemeType() == eImplicit), "Fully
110  // Implicit integration scheme cannot be handled by this routine." );
111 
112  int nvar = solvector->GetFirstDim();
113  int npoints = solvector->GetSecondDim();
114 
115  if (solvector->GetIntegrationSchemeData() !=
116  this) // FIMXE ... verify that this is correct...
117  {
118  // This branch will be taken when the solution vector
119  // (solvector) is set up for a different scheme than
120  // the object this method is called from. (typically
121  // needed to calculate the first time-levels of a
122  // multi-step scheme)
123 
124  // To do this kind of 'non-matching' integration, we
125  // perform the following three steps:
126  //
127  // 1: copy the required input information from the
128  // solution vector of the master scheme to the
129  // input solution vector of the current scheme
130  //
131  // 2: time-integrate for one step using the current
132  // scheme
133  //
134  // 3: copy the information contained in the output
135  // vector of the current scheme to the solution
136  // vector of the master scheme
137 
138  // STEP 1: copy the required input information from
139  // the solution vector of the master scheme
140  // to the input solution vector of the
141  // current scheme
142 
143  // 1.1 Determine which information is required for the
144  // current scheme
145 
146  DoubleArray y_n;
147  NekDouble t_n = 0;
148  DoubleArray dtFy_n;
149 
150  unsigned int nCurSchemeVals =
151  GetNmultiStepValues(); // number of required values of the current
152  // scheme
153  unsigned int nCurSchemeDers =
154  GetNmultiStepDerivs(); // number of required derivs of the current
155  // scheme
156  unsigned int nCurSchemeSteps =
157  GetNsteps(); // number of steps in the current scheme // FIMXE...
158  // what does it mean by "current scheme"... is this now
159  // a SchemeData issue?
160 
161  unsigned int nMasterSchemeVals =
162  solvector->GetNvalues(); // number of values of the master scheme
163  unsigned int nMasterSchemeDers =
164  solvector->GetNderivs(); // number of derivs of the master scheme
165 
166  // The arrays below contains information to which time-level
167  // the values and derivatives of the schemes belong
168  const Array<OneD, const unsigned int> &curTimeLevels =
170  const Array<OneD, const unsigned int> &masterTimeLevels =
171  solvector->GetTimeLevelOffset();
172 
173  // 1.2 Copy the required information from the master
174  // solution vector to the input solution vector of
175  // the current scheme
178  this); // input solution vector of the current scheme
179 
180  for (int n = 0; n < nCurSchemeVals; n++)
181  {
182  // Get the required value out of the master solution vector
183  // DoubleArray& y_n = solvector->GetValue ( curTimeLevels[n] );
184  // NekDouble t_n = solvector->GetValueTime( curTimeLevels[n] );
185 
186  y_n = solvector->GetValue(curTimeLevels[n]);
187  t_n = solvector->GetValueTime(curTimeLevels[n]);
188 
189  // Set the required value in the input solution
190  // vector of the current scheme
191  solvector_in->SetValue(curTimeLevels[n], y_n, t_n);
192  }
193 
194  for (int n = nCurSchemeVals; n < nCurSchemeSteps; n++)
195  {
196  // Get the required derivative out of the master
197  // solution vector
198  // DoubleArray& dtFy_n = solvector->GetDerivative (
199  // curTimeLevels[n] );
200  dtFy_n = solvector->GetDerivative(curTimeLevels[n]);
201 
202  // Set the required derivative in the input
203  // solution vector of the current scheme
204  solvector_in->SetDerivative(curTimeLevels[n], dtFy_n, deltaT);
205  }
206 
207  // STEP 2: time-integrate for one step using the
208  // current scheme
211  this, nvar,
212  npoints); // output solution vector of the current scheme
213 
214  // Integrate one step
215  TimeIntegrate(deltaT, solvector_in->GetSolutionVector(),
216  solvector_in->GetTimeVector(),
217  solvector_out->UpdateSolutionVector(),
218  solvector_out->UpdateTimeVector(), op);
219 
220  // STEP 3: copy the information contained in the
221  // output vector of the current scheme to the
222  // solution vector of the master scheme
223 
224  // 3.1 Check whether the current time scheme updates
225  // the most recent derivative that should be
226  // updated in the master scheme. If not,
227  // calculate the derivative. This can be done
228  // based upon the corresponding value and the
229  // DoOdeRhs operator.
230  bool CalcNewDeriv = false; // This flag inidicates whether
231  // the new derivative is availble
232  // in the output of the current
233  // scheme or whether it should be
234  // calculated.
235 
236  if (nMasterSchemeDers > 0)
237  {
238  if (nCurSchemeDers == 0)
239  {
240  CalcNewDeriv = true;
241  }
242  else
243  {
244  if (masterTimeLevels[nMasterSchemeVals] <
245  curTimeLevels[nCurSchemeVals])
246  {
247  CalcNewDeriv = true;
248  }
249  }
250  }
251 
252  if (CalcNewDeriv)
253  {
254  int newDerivTimeLevel =
255  masterTimeLevels[nMasterSchemeVals]; // Contains the
256  // time level at
257  // which the
258  // derivative of
259  // the master
260  // scheme is
261  // known.
262  // DoubleArray y_n;
263  // NekDouble t_n;
264  // If the time level corresponds to 0, calculate the
265  // derivative based upon the solution value at the new
266  // time-level
267  if (newDerivTimeLevel == 0)
268  {
269  y_n = solvector_out->GetValue(0);
270  t_n = solvector_out->GetValueTime(0);
271  }
272  // If the time level corresponds to 1, calculate the
273  // derivative based upon the solution value at the new
274  // old-level.
275  else if (newDerivTimeLevel == 1)
276  {
277  y_n = solvector->GetValue(0);
278  t_n = solvector->GetValueTime(0);
279  }
280  else
281  {
282  ASSERTL1(false, "Problems with initialising scheme");
283  }
284 
285  DoubleArray f_n(nvar);
286  for (int j = 0; j < nvar; j++)
287  {
288  f_n[j] = Array<OneD, NekDouble>(npoints);
289  }
290 
291  // Calculate the derivative
292  op.DoOdeRhs(y_n, f_n, t_n);
293 
294  // Multiply by dt (as required by the General Linear Method
295  // framework)
296  for (int j = 0; j < nvar; j++)
297  {
298  Vmath::Smul(npoints, deltaT, f_n[j], 1, f_n[j], 1);
299  }
300 
301  // Rotate the solution vector (i.e. updating without
302  // calculating/inserting new values).
303  solvector->RotateSolutionVector();
304 
305  // Set the calculated derivative in the master solution
306  // vector.
307  solvector->SetDerivative(newDerivTimeLevel, f_n, deltaT);
308  }
309  else
310  {
311  // Rotate the solution vector (i.e. updating without
312  // calculating/inserting new values).
313  solvector->RotateSolutionVector();
314  }
315 
316  // 1.2 Copy the information calculated using the
317  // current scheme from the output solution vector
318  // to the master solution vector
319  for (int n = 0; n < nCurSchemeVals; n++)
320  {
321  // Get the calculated value out of the output
322  // solution vector of the current scheme
323  // DoubleArray& y_n = solvector_out->GetValue ( curTimeLevels[n]
324  // );
325  // NekDouble t_n = solvector_out->GetValueTime( curTimeLevels[n]
326  // );
327  y_n = solvector_out->GetValue(curTimeLevels[n]);
328  t_n = solvector_out->GetValueTime(curTimeLevels[n]);
329 
330  // Set the calculated value in the master solution vector
331  solvector->SetValue(curTimeLevels[n], y_n, t_n);
332  }
333 
334  for (int n = nCurSchemeVals; n < nCurSchemeSteps; n++)
335  {
336  // Get the calculated derivative out of the output
337  // solution vector of the current scheme.
338  // DoubleArray& dtFy_n =
339  // solvector_out->GetDerivative (curTimeLevels[n]);
340  dtFy_n = solvector_out->GetDerivative(curTimeLevels[n]);
341 
342  // Set the calculated derivative in the master
343  // solution vector.
344  solvector->SetDerivative(curTimeLevels[n], dtFy_n, deltaT);
345  }
346  }
347  else
348  {
351  this, nvar, npoints);
352 
353  TimeIntegrate(deltaT, solvector->GetSolutionVector(),
354  solvector->GetTimeVector(),
355  solvector_new->UpdateSolutionVector(),
356  solvector_new->UpdateTimeVector(), op);
357 
358  solvector = solvector_new;
359  }
360 
361  return solvector->GetSolution();
362 
363 } // end TimeIntegrate()
364 
365 // Does the actual multi-stage multi-step integration.
367  const NekDouble deltaT, ConstTripleArray &y_old, ConstSingleArray &t_old,
368  TripleArray &y_new, SingleArray &t_new,
370 {
371  ASSERTL1(
372  CheckTimeIntegrateArguments(/*deltaT,*/ y_old, t_old, y_new, t_new, op),
373  "Arguments not well defined");
374 
376 
377  // Check if storage has already been initialised.
378  // If so, we just zero the temporary storage.
379  if (m_initialised && m_nvars == GetFirstDim(y_old) &&
380  m_npoints == GetSecondDim(y_old))
381  {
382  for (int j = 0; j < m_nvars; j++)
383  {
384  Vmath::Zero(m_npoints, m_tmp[j], 1);
385  }
386  }
387  else
388  {
389  m_nvars = GetFirstDim(y_old);
390  m_npoints = GetSecondDim(y_old);
391 
392  // First, calculate the various stage values and stage
393  // derivatives (this is the multi-stage part of the method)
394  // - m_Y corresponds to the stage values
395  // - m_F corresponds to the stage derivatives
396  // - m_T corresponds to the time at the different stages
397  // - m_tmp corresponds to the explicit right hand side of
398  // each stage equation
399  // (for explicit schemes, this correspond to m_Y)
400 
401  // Allocate memory for the arrays m_Y and m_F and m_tmp. The same
402  // storage will be used for every stage -> m_Y is a
403  // DoubleArray
405  for (int j = 0; j < m_nvars; j++)
406  {
408  }
409 
410  // The same storage will be used for every stage -> m_tmp is a
411  // DoubleArray
412  if (type == eExplicit || m_schemeType == eExponential)
413  {
414  m_Y = m_tmp;
415  }
416  else
417  {
419  for (int j = 0; j < m_nvars; j++)
420  {
422  }
423  }
424 
425  // Different storage for every stage derivative as the data
426  // will be re-used to update the solution -> m_F is a TripleArray
428  for (int i = 0; i < m_numstages; ++i)
429  {
430  m_F[i] = DoubleArray(m_nvars);
431  for (int j = 0; j < m_nvars; j++)
432  {
433  m_F[i][j] = Array<OneD, NekDouble>(m_npoints, 0.0);
434  }
435  }
436 
437  if (type == eIMEX)
438  {
440  for (int i = 0; i < m_numstages; ++i)
441  {
443  for (int j = 0; j < m_nvars; j++)
444  {
446  }
447  }
448  }
449 
450  // Finally, flag that the memory has been initialised.
451  m_initialised = true;
452  } // end else
453 
454  // For an exponential integrator if the time increment or the
455  // number of variables has changed then the exponenial matrices
456  // must be recomputed.
457  if (type == eExponential)
458  {
459  if (m_lastDeltaT != deltaT || m_lastNVars != GetFirstDim(y_old))
460  {
462  ->InitializeSecondaryData(this, deltaT);
463 
464  m_lastDeltaT = deltaT;
465  m_lastNVars = GetFirstDim(y_old);
466  }
467  }
468 
469  // The loop below calculates the stage values and derivatives
470  for (int stage = 0; stage < m_numstages; stage++)
471  {
472  if ((stage == 0) && m_firstStageEqualsOldSolution)
473  {
474  for (int k = 0; k < m_nvars; k++)
475  {
476  Vmath::Vcopy(m_npoints, y_old[0][k], 1, m_Y[k], 1);
477  }
478 
479  m_T = t_old[0];
480  }
481  else
482  {
483  // The stage values m_Y are a linear combination of:
484  // 1: The stage derivatives:
485  if (stage != 0)
486  {
487  for (int k = 0; k < m_nvars; k++)
488  {
489  if (type == eExponential)
490  {
491  Vmath::Smul(m_npoints, deltaT * m_A_phi[k][stage][0],
492  m_F[0][k], 1, m_tmp[k], 1);
493  }
494  else
495  {
496  Vmath::Smul(m_npoints, deltaT * A(stage, 0), m_F[0][k],
497  1, m_tmp[k], 1);
498  }
499 
500  if (type == eIMEX)
501  {
502  Vmath::Svtvp(m_npoints, deltaT * A_IMEX(stage, 0),
503  m_F_IMEX[0][k], 1, m_tmp[k], 1, m_tmp[k],
504  1);
505  }
506  }
507  }
508 
509  m_T = A(stage, 0) * deltaT;
510 
511  for (int j = 1; j < stage; j++)
512  {
513  for (int k = 0; k < m_nvars; k++)
514  {
515  if (type == eExponential)
516  {
517  Vmath::Svtvp(m_npoints, deltaT * m_A_phi[k][stage][j],
518  m_F[j][k], 1, m_tmp[k], 1, m_tmp[k], 1);
519  }
520  else
521  {
522  Vmath::Svtvp(m_npoints, deltaT * A(stage, j), m_F[j][k],
523  1, m_tmp[k], 1, m_tmp[k], 1);
524  }
525 
526  if (type == eIMEX)
527  {
528  Vmath::Svtvp(m_npoints, deltaT * A_IMEX(stage, j),
529  m_F_IMEX[j][k], 1, m_tmp[k], 1, m_tmp[k],
530  1);
531  }
532  }
533 
534  m_T += A(stage, j) * deltaT;
535  }
536 
537  // 2: The imported multi-step solution of the previous time level:
538  for (int j = 0; j < m_numsteps; j++)
539  {
540  for (int k = 0; k < m_nvars; k++)
541  {
542  if (type == eExponential)
543  {
544  Vmath::Svtvp(m_npoints, m_U_phi[k][stage][j],
545  y_old[j][k], 1, m_tmp[k], 1, m_tmp[k], 1);
546  }
547  else
548  {
549  Vmath::Svtvp(m_npoints, U(stage, j), y_old[j][k], 1,
550  m_tmp[k], 1, m_tmp[k], 1);
551  }
552  }
553 
554  m_T += U(stage, j) * t_old[j];
555  }
556  } // end else
557 
558  // Calculate the stage derivative based upon the stage value
559  if (type == eDiagonallyImplicit)
560  {
561  if ((stage == 0) && m_firstStageEqualsOldSolution)
562  {
563  op.DoOdeRhs(m_Y, m_F[stage], m_T);
564  }
565  else
566  {
567  if (m_numstages == 1)
568  {
569  m_T = t_old[0] + deltaT;
570  }
571  else
572  {
573  m_T = t_old[0];
574  for (int j = 0; j <= stage; ++j)
575  {
576  m_T += A(stage, j) * deltaT;
577  }
578  }
579 
580  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(stage, stage) * deltaT);
581 
582  for (int k = 0; k < m_nvars; ++k)
583  {
584  Vmath::Vsub(m_npoints, m_Y[k], 1, m_tmp[k], 1,
585  m_F[stage][k], 1);
586  Vmath::Smul(m_npoints, 1.0 / (A(stage, stage) * deltaT),
587  m_F[stage][k], 1, m_F[stage][k], 1);
588  }
589  }
590  }
591  else if (type == eIMEX)
592  {
593  if (m_numstages == 1)
594  {
595  m_T = t_old[0] + deltaT;
596  }
597  else
598  {
599  m_T = t_old[0];
600  for (int j = 0; j <= stage; ++j)
601  {
602  m_T += A(stage, j) * deltaT;
603  }
604  }
605 
606  if (fabs(A(stage, stage)) > NekConstants::kNekZeroTol)
607  {
608  op.DoImplicitSolve(m_tmp, m_Y, m_T, A(stage, stage) * deltaT);
609 
610  for (int k = 0; k < m_nvars; k++)
611  {
612  Vmath::Vsub(m_npoints, m_Y[k], 1, m_tmp[k], 1,
613  m_F[stage][k], 1);
614  Vmath::Smul(m_npoints, 1.0 / (A(stage, stage) * deltaT),
615  m_F[stage][k], 1, m_F[stage][k], 1);
616  }
617  }
618  op.DoOdeRhs(m_Y, m_F_IMEX[stage], m_T);
619  }
620  else if (type == eExplicit || m_schemeType == eExponential)
621  {
622  // Avoid projecting the same solution twice
623  if (!((stage == 0) && m_firstStageEqualsOldSolution))
624  {
625  // ensure solution is in correct space
626  op.DoProjection(m_Y, m_Y, m_T);
627  }
628 
629  op.DoOdeRhs(m_Y, m_F[stage], m_T);
630  }
631  }
632 
633  // Next, the solution vector y at the new time level will
634  // be calculated.
635  //
636  // For multi-step methods, this includes updating the
637  // values of the auxiliary parameters
638  //
639  // The loop below calculates the solution at the new time
640  // level
641  //
642  // If last stage equals the new solution, the new solution
643  // needs not be calculated explicitly but can simply be
644  // copied. This saves a solve.
645  int i_start = 0;
647  {
648  for (int k = 0; k < m_nvars; k++)
649  {
650  Vmath::Vcopy(m_npoints, m_Y[k], 1, y_new[0][k], 1);
651  }
652 
653  if (m_numstages == 1 && type == eIMEX)
654  {
655  t_new[0] = t_old[0] + deltaT;
656  }
657  else
658  {
659  t_new[0] = B(0, 0) * deltaT;
660 
661  for (int j = 1; j < m_numstages; j++)
662  {
663  t_new[0] += B(0, j) * deltaT;
664  }
665 
666  for (int j = 0; j < m_numsteps; j++)
667  {
668  t_new[0] += V(0, j) * t_old[j];
669  }
670  }
671 
672  i_start = 1;
673  }
674 
675  for (int i = i_start; i < m_numsteps; i++)
676  {
677  // The solution at the new time level is a linear
678  // combination of:
679  // 1: the stage derivatives
680  for (int k = 0; k < m_nvars; k++)
681  {
682  if (type == eExponential)
683  {
684  Vmath::Smul(m_npoints, deltaT * m_B_phi[k][i][0], m_F[0][k], 1,
685  y_new[i][k], 1);
686  }
687  else
688  {
689  Vmath::Smul(m_npoints, deltaT * B(i, 0), m_F[0][k], 1,
690  y_new[i][k], 1);
691  }
692 
693  if (type == eIMEX)
694  {
695  Vmath::Svtvp(m_npoints, deltaT * B_IMEX(i, 0), m_F_IMEX[0][k],
696  1, y_new[i][k], 1, y_new[i][k], 1);
697  }
698  }
699 
700  if (m_numstages != 1 || type != eIMEX)
701  {
702  t_new[i] = B(i, 0) * deltaT;
703  }
704 
705  for (int j = 1; j < m_numstages; j++)
706  {
707  for (int k = 0; k < m_nvars; k++)
708  {
709  if (type == eExponential)
710  {
711  Vmath::Svtvp(m_npoints, deltaT * m_B_phi[k][i][j],
712  m_F[j][k], 1, y_new[i][k], 1, y_new[i][k], 1);
713  }
714  else
715  {
716  Vmath::Svtvp(m_npoints, deltaT * B(i, j), m_F[j][k], 1,
717  y_new[i][k], 1, y_new[i][k], 1);
718  }
719 
720  if (type == eIMEX)
721  {
722  Vmath::Svtvp(m_npoints, deltaT * B_IMEX(i, j),
723  m_F_IMEX[j][k], 1, y_new[i][k], 1, y_new[i][k],
724  1);
725  }
726  }
727 
728  if (m_numstages != 1 || type != eIMEX)
729  {
730  t_new[i] += B(i, j) * deltaT;
731  }
732  }
733 
734  // 2: the imported multi-step solution of the previous
735  // time level
736  for (int j = 0; j < m_numsteps; j++)
737  {
738  for (int k = 0; k < m_nvars; k++)
739  {
740  if (type == eExponential)
741  {
742  Vmath::Svtvp(m_npoints, m_V_phi[k][i][j], y_old[j][k], 1,
743  y_new[i][k], 1, y_new[i][k], 1);
744  }
745  else
746  {
747  Vmath::Svtvp(m_npoints, V(i, j), y_old[j][k], 1,
748  y_new[i][k], 1, y_new[i][k], 1);
749  }
750  }
751 
752  if (m_numstages != 1 || type != eIMEX)
753  {
754  t_new[i] += V(i, j) * t_old[j];
755  }
756  }
757  }
758 
759  // Ensure that the new solution is projected if necessary
760  if (type == eExplicit || m_schemeType == eExponential)
761  {
762  op.DoProjection(y_new[0], y_new[0], t_new[0]);
763  }
764 
765 } // end TimeIntegrate()
766 
768 {
769  // First stage equals old solution if:
770  // 1. the first row of the coefficient matrix A consists of zeros
771  // 2. U[0][0] is equal to one and all other first row entries of U are zero
772 
773  // 1. Check first condition
774  for (int m = 0; m < m_A.size(); m++)
775  {
776  for (int i = 0; i < m_numstages; i++)
777  {
778  if (fabs(m_A[m][0][i]) > NekConstants::kNekZeroTol)
779  {
781  return;
782  }
783  }
784  }
785 
786  // 2. Check second condition
787  if (fabs(m_U[0][0] - 1.0) > NekConstants::kNekZeroTol)
788  {
790  return;
791  }
792 
793  for (int i = 1; i < m_numsteps; i++)
794  {
795  if (fabs(m_U[0][i]) > NekConstants::kNekZeroTol)
796  {
798  return;
799  }
800  }
801 
803 }
804 
806 {
807  // Last stage equals new solution if:
808  // 1. the last row of the coefficient matrix A is equal to the first row of
809  // matrix B
810  // 2. the last row of the coefficient matrix U is equal to the first row of
811  // matrix V
812 
813  // 1. Check first condition
814  for (int m = 0; m < m_A.size(); m++)
815  {
816  for (int i = 0; i < m_numstages; i++)
817  {
818  if (fabs(m_A[m][m_numstages - 1][i] - m_B[m][0][i]) >
820  {
822  return;
823  }
824  }
825  }
826 
827  // 2. Check second condition
828  for (int i = 0; i < m_numsteps; i++)
829  {
830  if (fabs(m_U[m_numstages - 1][i] - m_V[0][i]) >
832  {
834  return;
835  }
836  }
837 
839 }
840 
842 {
843 #ifdef DEBUG
844  int IMEXdim = m_A.size();
845  int dim = m_A[0].GetRows();
846 
848 
849  if (m_schemeType == eExponential)
850  vertype[0] = eExponential;
851 
852  for (int m = 0; m < IMEXdim; m++)
853  {
854  for (int i = 0; i < dim; i++)
855  {
856  if (fabs(m_A[m][i][i]) > NekConstants::kNekZeroTol)
857  {
858  vertype[m] = eDiagonallyImplicit;
859  }
860  }
861 
862  for (int i = 0; i < dim; i++)
863  {
864  for (int j = i + 1; j < dim; j++)
865  {
866  if (fabs(m_A[m][i][j]) > NekConstants::kNekZeroTol)
867  {
868  vertype[m] = eImplicit;
869  ASSERTL1(false, "Fully Implicit schemes cannnot be handled "
870  "by the TimeIntegrationScheme class");
871  }
872  }
873  }
874  }
875 
876  if (IMEXdim == 2)
877  {
878  ASSERTL1(m_B.size() == 2,
879  "Coefficient Matrix B should have an "
880  "implicit and explicit part for IMEX schemes");
881 
882  if ((vertype[0] == eDiagonallyImplicit) && (vertype[1] == eExplicit))
883  {
884  vertype[0] = eIMEX;
885  }
886  else
887  {
888  ASSERTL1(false, "This is not a proper IMEX scheme");
889  }
890  }
891 
892  ASSERTL1(vertype[0] == m_schemeType,
893  "Time integration scheme coefficients do not match its type");
894 #endif
895 }
896 
898  ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new,
899  SingleArray &t_new, const TimeIntegrationSchemeOperators &op) const
900 {
901 #ifdef DEBUG
902  boost::ignore_unused(op);
903 #else
904  boost::ignore_unused(y_old, t_old, y_new, t_new, op);
905 #endif
906 
907  // Check if arrays are all of consistent size
908 
909  ASSERTL1(y_old.size() == m_numsteps, "Non-matching number of steps.");
910  ASSERTL1(y_new.size() == m_numsteps, "Non-matching number of steps.");
911 
912  ASSERTL1(y_old[0].size() == y_new[0].size(),
913  "Non-matching number of variables.");
914  ASSERTL1(y_old[0][0].size() == y_new[0][0].size(),
915  "Non-matching number of coefficients.");
916 
917  ASSERTL1(t_old.size() == m_numsteps, "Non-matching number of steps.");
918  ASSERTL1(t_new.size() == m_numsteps, "Non-matching number of steps.");
919 
920  return true;
921 }
922 
923 // Friend Operators
924 std::ostream &operator<<(std::ostream &os,
926 {
927  return operator<<(os, *rhs);
928 }
929 
930 std::ostream &operator<<(std::ostream &os,
931  const TimeIntegrationAlgorithmGLM &rhs)
932 {
933  int r = rhs.m_numsteps;
934  int s = rhs.m_numstages;
936 
937  int oswidth = 9;
938  int osprecision = 6;
939 
940  os << "Time Integration Scheme (Master): " << rhs.m_parent->GetFullName()
941  << "\n"
942  << "Time Integration Phase : " << rhs.m_name << "\n"
943  << "- number of steps: " << r << "\n"
944  << "- number of stages: " << s << "\n"
945  << "General linear method tableau:\n";
946 
947  for (int i = 0; i < s; i++)
948  {
949  for (int j = 0; j < s; j++)
950  {
951  os.width(oswidth);
952  os.precision(osprecision);
953  os << std::right << rhs.A(i, j) << " ";
954  }
955  if (type == eIMEX)
956  {
957  os << " '";
958  for (int j = 0; j < s; j++)
959  {
960  os.width(oswidth);
961  os.precision(osprecision);
962  os << std::right << rhs.A_IMEX(i, j) << " ";
963  }
964  }
965 
966  os << " |";
967 
968  for (int j = 0; j < r; j++)
969  {
970  os.width(oswidth);
971  os.precision(osprecision);
972  os << std::right << rhs.U(i, j);
973  }
974  os << std::endl;
975  }
976 
977  int imexflag = (type == eIMEX) ? 2 : 1;
978  for (int i = 0; i < (r + imexflag * s) * (oswidth + 1) + imexflag * 2 - 1;
979  i++)
980  {
981  os << "-";
982  }
983  os << std::endl;
984 
985  for (int i = 0; i < r; i++)
986  {
987  for (int j = 0; j < s; j++)
988  {
989  os.width(oswidth);
990  os.precision(osprecision);
991  os << std::right << rhs.B(i, j) << " ";
992  }
993  if (type == eIMEX)
994  {
995  os << " '";
996  for (int j = 0; j < s; j++)
997  {
998  os.width(oswidth);
999  os.precision(osprecision);
1000  os << std::right << rhs.B_IMEX(i, j) << " ";
1001  }
1002  }
1003 
1004  os << " |";
1005 
1006  for (int j = 0; j < r; j++)
1007  {
1008  os.width(oswidth);
1009  os.precision(osprecision);
1010  os << std::right << rhs.V(i, j);
1011  }
1012 
1013  os << " |";
1014 
1015  os.width(oswidth);
1016  os.precision(osprecision);
1017  os << std::right << rhs.m_timeLevelOffset[i];
1018 
1019  if (i < rhs.m_numMultiStepValues)
1020  {
1021  os << std::right << " value";
1022  }
1023  else
1024  {
1025  os << std::right << " derivative";
1026  }
1027 
1028  os << std::endl;
1029  }
1030 
1031  if (type == eExponential)
1032  {
1033  for (int k = 0; k < rhs.m_nvars; k++)
1034  {
1035  os << std::endl
1036  << "General linear method exponential tableau for variable " << k
1037  << ":\n";
1038 
1039  for (int i = 0; i < s; i++)
1040  {
1041  for (int j = 0; j < s; j++)
1042  {
1043  os.width(oswidth);
1044  os.precision(osprecision);
1045  os << std::right << rhs.m_A_phi[k][i][j] << " ";
1046  }
1047 
1048  os << " |";
1049 
1050  for (int j = 0; j < r; j++)
1051  {
1052  os.width(oswidth);
1053  os.precision(osprecision);
1054  os << std::right << rhs.m_U_phi[k][i][j];
1055  }
1056  os << std::endl;
1057  }
1058 
1059  int imexflag = (type == eIMEX) ? 2 : 1;
1060  for (int i = 0;
1061  i < (r + imexflag * s) * (oswidth + 1) + imexflag * 2 - 1; i++)
1062  {
1063  os << "-";
1064  }
1065  os << std::endl;
1066 
1067  for (int i = 0; i < r; i++)
1068  {
1069  for (int j = 0; j < s; j++)
1070  {
1071  os.width(oswidth);
1072  os.precision(osprecision);
1073  os << std::right << rhs.m_B_phi[k][i][j] << " ";
1074  }
1075 
1076  os << " |";
1077 
1078  for (int j = 0; j < r; j++)
1079  {
1080  os.width(oswidth);
1081  os.precision(osprecision);
1082  os << std::right << rhs.m_V_phi[k][i][j];
1083  }
1084  os << std::endl;
1085  }
1086  }
1087  }
1088 
1089  return os;
1090 } // end function operator<<
1091 
1092 } // end namespace LibUtilities
1093 } // namespace Nektar
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
NekDouble A(const unsigned int i, const unsigned int j) const
TripleArray m_F
Explicit RHS of each stage equation.
bool CheckTimeIntegrateArguments(ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new, SingleArray &t_new, const TimeIntegrationSchemeOperators &op) const
int m_npoints
The size of inner data which is stored for reuse.
DoubleArray m_tmp
Array containing the stage values.
NekDouble B_IMEX(const unsigned int i, const unsigned int j) const
NekDouble V(const unsigned int i, const unsigned int j) const
const TimeIntegrationScheme * m_parent
Parent scheme object.
NekDouble m_T
Used to store the Explicit stage derivative of IMEX schemes.
const Array< OneD, const unsigned int > & GetTimeLevelOffset() const
NekDouble A_IMEX(const unsigned int i, const unsigned int j) const
NekDouble B(const unsigned int i, const unsigned int j) const
LUE ConstDoubleArray & TimeIntegrate(const NekDouble deltaT, TimeIntegrationSolutionGLMSharedPtr &y, const TimeIntegrationSchemeOperators &op)
Explicit integration of an ODE.
LUE TimeIntegrationSolutionGLMSharedPtr InitializeData(const NekDouble deltaT, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)
This function initialises the time integration scheme.
NekDouble m_lastDeltaT
bool to identify array initialization
NekDouble U(const unsigned int i, const unsigned int j) const
TripleArray m_F_IMEX
Array corresponding to the stage Derivatives.
Base class for GLM time integration schemes.
Binds a set of functions for use by time integration schemes.
void DoImplicitSolve(InArrayType &inarray, OutArrayType &outarray, const NekDouble time, const NekDouble lambda) const
void DoProjection(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
void DoOdeRhs(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:168
AT< AT< NekDouble > > DoubleArray
std::shared_ptr< TimeIntegrationAlgorithmGLM > TimeIntegrationAlgorithmGLMSharedPtr
@ eImplicit
Fully implicit scheme.
@ eExplicit
Formally explicit scheme.
@ eExponential
Exponential scheme.
@ eDiagonallyImplicit
Diagonally implicit scheme (e.g. the DIRK schemes)
@ eIMEX
Implicit Explicit General Linear Method.
std::shared_ptr< TimeIntegrationSolutionGLM > TimeIntegrationSolutionGLMSharedPtr
AT< AT< AT< NekDouble > > > TripleArray
std::ostream & operator<<(std::ostream &os, const BasisKey &rhs)
static const NekDouble kNekZeroTol
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
double NekDouble
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:622
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
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:419