Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TimeIntegrationScheme.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: TimeIntegrationScheme.h
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: Header file of time integration scheme base class
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #ifndef NEKTAR_LIB_UTILITIES_FOUNDATIONS_TIMEINTEGRATIONSCHEME_H
37 #define NEKTAR_LIB_UTILITIES_FOUNDATIONS_TIMEINTEGRATIONSCHEME_H
38 
44 
45 namespace Nektar
46 {
47  namespace LibUtilities
48  {
49  // Forward declaration of all classes in this file
50  class TimeIntegrationScheme;
52 
53  // typedefs
54  typedef boost::shared_ptr<TimeIntegrationScheme> TimeIntegrationSchemeSharedPtr;
55  typedef std::vector<TimeIntegrationSchemeSharedPtr> TimeIntegrationSchemeVector;
57  typedef boost::shared_ptr<TimeIntegrationSolution> TimeIntegrationSolutionSharedPtr;
58  typedef std::vector<TimeIntegrationSolutionSharedPtr> TimeIntegrationSolutionVector;
60 
61  // =========================================================================
62  // ==== ENUM LIST OF ALL SUPPORTED INTEGRATION SCHEMES
63  // =========================================================================
65  {
67  eAdamsBashforthOrder1, //!< Adams-Bashforth Forward multi-step scheme of order 1
68  eAdamsBashforthOrder2, //!< Adams-Bashforth Forward multi-step scheme of order 2
69  eAdamsBashforthOrder3, //!< Adams-Bashforth Forward multi-step scheme of order 3
70  eAdamsMoultonOrder1, //!< Adams-Moulton Forward multi-step scheme of order 1
71  eAdamsMoultonOrder2, //!< Adams-Moulton Forward multi-step scheme of order 2
72  eBDFImplicitOrder1, //!< BDF multi-step scheme of order 1 (implicit)
73  eBDFImplicitOrder2, //!< BDF multi-step scheme of order 2 (implicit)
74  eClassicalRungeKutta4, //!< Runge-Kutta multi-stage scheme 4th order explicit (old name)
75  eRungeKutta4, //!< Classical RungeKutta4 method (new name for eClassicalRungeKutta4)
76  eRungeKutta3_SSP, //!< Nonlinear SSP RungeKutta3 explicit
77  eRungeKutta2_ImprovedEuler, //!< Improved RungeKutta2 explicit (old name meaning Heun's method)
78  eRungeKutta2_SSP, //!< Nonlinear SSP RungeKutta2 explicit (surrogate for eRungeKutta2_ImprovedEuler)
79  eForwardEuler, //!< Forward Euler scheme
80  eBackwardEuler, //!< Backward Euler scheme
81  eIMEXOrder1, //!< IMEX 1st order scheme using Euler Backwards/Euler Forwards
82  eIMEXOrder2, //!< IMEX 2nd order scheme using Backward Different Formula & Extrapolation
83  eIMEXOrder3, //!< IMEX 3rd order scheme using Backward Different Formula & Extrapolation
84  eMidpoint, //!< midpoint method (old name)
85  eRungeKutta2, //!< Classical RungeKutta2 method (new name for eMidpoint)
86  eDIRKOrder2, //!< Diagonally Implicit Runge Kutta scheme of order 3
87  eDIRKOrder3, //!< Diagonally Implicit Runge Kutta scheme of order 3
88  eCNAB, //!< Crank-Nicolson/Adams-Bashforth Order 2 (CNAB)
89  eIMEXGear, //!< IMEX Gear Order 2
90  eMCNAB, //!< Modified Crank-Nicolson/Adams-Bashforth Order 2 (MCNAB)
91  eIMEXdirk_1_1_1, //!< Forward-Backward Euler IMEX DIRK(1,1,1)
92  eIMEXdirk_1_2_1, //!< Forward-Backward Euler IMEX DIRK(1,2,1)
93  eIMEXdirk_1_2_2, //!< Implicit-Explicit Midpoint IMEX DIRK(1,2,2)
94  eIMEXdirk_2_2_2, //!< L-stable, two stage, second order IMEX DIRK(2,2,2)
95  eIMEXdirk_2_3_2, //!< L-stable, three stage, third order IMEX DIRK(3,4,3)
96  eIMEXdirk_2_3_3, //!< L-stable, two stage, third order IMEX DIRK(2,3,3)
97  eIMEXdirk_3_4_3, //!< L-stable, three stage, third order IMEX DIRK(3,4,3)
98  eIMEXdirk_4_4_3, //!< L-stable, four stage, third order IMEX DIRK(4,4,3)
99  SIZE_TimeIntegrationMethod //!< Length of enum list
100  };
101 
102  const char* const TimeIntegrationMethodMap[] =
103  {
104  "NoTimeIntegrationMethod",
105  "AdamsBashforthOrder1",
106  "AdamsBashforthOrder2",
107  "AdamsBashforthOrder3",
108  "AdamsMoultonOrder1",
109  "AdamsMoultonOrder2",
110  "BDFImplicitOrder1",
111  "BDFImplicitOrder2",
112  "ClassicalRungeKutta4",
113  "RungeKutta4",
114  "RungeKutta3_SSP",
115  "RungeKutta2_ImprovedEuler",
116  "RungeKutta2_SSP",
117  "ForwardEuler",
118  "BackwardEuler",
119  "IMEXOrder1",
120  "IMEXOrder2",
121  "IMEXOrder3",
122  "Midpoint",
123  "RungeKutta2",
124  "DIRKOrder2",
125  "DIRKOrder3",
126  "CNAB",
127  "IMEXGear",
128  "MCNAB",
129  "IMEXdirk_1_1_1",
130  "IMEXdirk_1_2_1",
131  "IMEXdirk_1_2_2",
132  "IMEXdirk_2_2_2",
133  "IMEXdirk_2_3_2",
134  "IMEXdirk_2_3_3",
135  "IMEXdirk_3_4_3",
136  "IMEXdirk_4_4_3",
137  };
138 
140  {
142  eExplicit, //!< Formally explicit scheme
143  eDiagonallyImplicit, //!< Diagonally implicit scheme (e.g. the DIRK schemes)
144  eIMEX, //!< Implicit Explicit General Linear Method
145  eImplicit, //!< Fully implicit scheme
146  };
147 
148  const char* const TimeIntegrationSchemeTypeMap[] =
149  {
150  "NoTimeIntegrationSchemeType",
151  "Explicit",
152  "DiagonallyImplicit",
153  "IMEX",
154  "Implicit"
155  };
156 
157  // =====================================================================
158 
159  // =====================================================================
160  // ==== DEFINITION OF THE CLASS TimeIntegrationSchemeOperators
161  // =====================================================================
163  {
164  public:
167 
168  typedef boost::function< void (InArrayType&, OutArrayType&, const NekDouble) > FunctorType1;
169  typedef boost::function< void (InArrayType&, OutArrayType&, const NekDouble, const NekDouble) > FunctorType2;
170 
171  typedef const FunctorType1& ConstFunctorType1Ref;
172  typedef const FunctorType2& ConstFunctorType2Ref;
173 
176 
178  m_functors1(4),
179  m_functors2(1)
180  {
181  }
182 
183  template<typename FuncPointerT, typename ObjectPointerT>
184  void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
185  {
186  m_functors1[0] = boost::bind(func, obj, _1, _2, _3);
187  }
188 
189  template<typename FuncPointerT, typename ObjectPointerT>
190  void DefineOdeExplicitRhs(FuncPointerT func, ObjectPointerT obj)
191  {
192  m_functors1[1] = boost::bind(func, obj, _1, _2, _3);
193  }
194 
195  template<typename FuncPointerT, typename ObjectPointerT>
196  void DefineOdeImplicitRhs(FuncPointerT func, ObjectPointerT obj)
197  {
198  m_functors1[2] = boost::bind(func, obj, _1, _2, _3);
199  }
200 
201  template<typename FuncPointerT, typename ObjectPointerT>
202  void DefineProjection(FuncPointerT func, ObjectPointerT obj)
203  {
204  m_functors1[3] = boost::bind(func, obj, _1, _2, _3);
205  }
206 
207  template<typename FuncPointerT, typename ObjectPointerT>
208  void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
209  {
210  m_functors2[0] = boost::bind(func, obj, _1, _2, _3, _4);
211  }
212 
213 
214  inline void DoOdeRhs(InArrayType &inarray,
215  OutArrayType &outarray,
216  const NekDouble time) const
217  {
218  ASSERTL1(!(m_functors1[0].empty()),"OdeRhs should be defined for this time integration scheme");
219  m_functors1[0](inarray,outarray,time);
220  }
221 
222  inline void DoOdeExplicitRhs(InArrayType &inarray,
223  OutArrayType &outarray,
224  const NekDouble time) const
225  {
226  ASSERTL1(!(m_functors1[1].empty()),"OdeExplicitRhs should be defined for this time integration scheme");
227  m_functors1[1](inarray,outarray,time);
228  }
229 
230  inline void DoOdeImplicitRhs(InArrayType &inarray,
231  OutArrayType &outarray,
232  const NekDouble time) const
233  {
234  ASSERTL1(!(m_functors1[2].empty()),"OdeImplictRhs should be defined for this time integration scheme");
235  m_functors1[2](inarray,outarray,time);
236  }
237 
238  inline void DoProjection(InArrayType &inarray,
239  OutArrayType &outarray,
240  const NekDouble time) const
241  {
242  ASSERTL1(!(m_functors1[3].empty()),"Projection operation should be defined for this time integration scheme");
243  m_functors1[3](inarray,outarray,time);
244  }
245 
246  inline void DoImplicitSolve(InArrayType &inarray,
247  OutArrayType &outarray,
248  const NekDouble time,
249  const NekDouble lambda) const
250  {
251  ASSERTL1(!(m_functors2[0].empty()),"ImplicitSolve should be defined for this time integration scheme");
252  m_functors2[0](inarray,outarray,time,lambda);
253  }
254 
255 
256  protected:
257  FunctorType1Array m_functors1;
258  FunctorType2Array m_functors2;
259 
260  private:
261 
262  };
263 
264  // =========================================================================
265  // ==== DEFINITION OF THE CLASS TimeIntegrationSchemeKey
266  // =========================================================================
268  {
269  public:
270 
271  struct opLess
272  {
274  };
275 
277  m_method(method)
278  {
279  }
280 
282  {
283  }
284 
286  {
287  *this = key; // defer to assignment operator
288  }
289 
291  {
292  m_method = key.m_method;
293 
294  return *this;
295  }
296 
298  {
299  return m_method;
300  }
301 
302  inline bool operator==(const TimeIntegrationSchemeKey &key)
303  {
304  return (m_method == key.m_method);
305  }
306 
307  inline bool operator== (const TimeIntegrationSchemeKey *y)
308  {
309  return (*this == *y);
310  }
311 
313  {
314  return (!(*this == y));
315  }
316 
318  {
319  return (!(*this == *y));
320  }
321 
325 
326  protected:
327  TimeIntegrationMethod m_method; //!< integration method
328 
329  private:
330  // This should never be called
332  {
333  NEKERROR(ErrorUtil::efatal,"Default Constructor for the TimeIntegrationSchemeKey class should not be called");
334  }
335  };
336 
337  static const TimeIntegrationSchemeKey NullTimeIntegrationSchemeKey(eNoTimeIntegrationMethod);
338  LIB_UTILITIES_EXPORT bool operator==(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs);
339  LIB_UTILITIES_EXPORT bool operator<(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs);
340  LIB_UTILITIES_EXPORT std::ostream& operator<<(std::ostream& os, const TimeIntegrationSchemeKey& rhs);
341 
342  // =====================================================================
343 
344 
345  // =====================================================================
346  // ==== DEFINITION OF THE FUNCTION TimeIntegrationSolution
347  // =====================================================================
348  // This is the interface the user can use to get hold of the different
349  // time integration schemes. It returns you the NekManager
350  // singleton which manages all the time integration schemes...
351  //
352  typedef NekManager<TimeIntegrationSchemeKey,
353  TimeIntegrationScheme,
355  LIB_UTILITIES_EXPORT TimeIntegrationSchemeManagerT &TimeIntegrationSchemeManager(void);
356  // =====================================================================
357 
358  // =====================================================================
359  // ==== DEFINITION OF THE CLASS TimeIntegrationScheme
360  // =====================================================================
361  class TimeIntegrationScheme
362  {
363  public:
364 
371  typedef boost::function< void (ConstDoubleArray&, DoubleArray&, const NekDouble) > FunctorType1;
372  typedef boost::function< void (ConstDoubleArray&, DoubleArray&, const NekDouble, const NekDouble) > FunctorType2;
373 
374  public:
375 
377  {
378  }
379 
380  inline const TimeIntegrationSchemeKey& GetIntegrationSchemeKey() const
381  {
382  return m_schemeKey;
383  }
384 
386  {
388  }
389 
391  {
392  return m_schemeType;
393  }
394 
395  inline NekDouble A(const unsigned int i, const unsigned int j) const
396  {
397  return m_A[0][i][j];
398  }
399 
400  inline NekDouble B(const unsigned int i, const unsigned int j) const
401  {
402  return m_B[0][i][j];
403  }
404 
405  inline NekDouble U(const unsigned int i, const unsigned int j) const
406  {
407  return m_U[i][j];
408  }
409 
410  inline NekDouble V(const unsigned int i, const unsigned int j) const
411  {
412  return m_V[i][j];
413  }
414 
415  inline NekDouble A_IMEX(const unsigned int i, const unsigned int j) const
416  {
417  return m_A[1][i][j];
418  }
419 
420  inline NekDouble B_IMEX(const unsigned int i, const unsigned int j) const
421  {
422  return m_B[1][i][j];
423  }
424 
425  inline unsigned int GetNsteps(void) const
426  {
427  return m_numsteps;
428  }
429 
430  inline unsigned int GetNstages(void) const
431  {
432  return m_numstages;
433  }
434 
435  inline unsigned int GetNmultiStepValues(void) const
436  {
437  return m_numMultiStepValues;
438  }
439 
440  inline unsigned int GetNmultiStepDerivs(void) const
441  {
442  return m_numMultiStepDerivs;
443  }
444 
445  /**
446  * \brief This function initialises the time integration
447  * scheme
448  *
449  * Given the solution at the initial time level
450  * \f$\boldsymbol{y}(t_0)\f$, this function generates the
451  * vectors \f$\boldsymbol{y}^{[n]}\f$ and \f$t^{[n]}\f$
452  * needed to evaluate the time integration scheme
453  * formulated as a General Linear Method. These vectors
454  * are embedded in an object of the class
455  * TimeIntegrationSolution. This class is the abstraction
456  * of the input and output vectors of the General Linear
457  * Method.
458  *
459  * For single-step methods, this function is trivial as it
460  * just wraps a TimeIntegrationSolution object around the
461  * given initial values and initial time. However, for
462  * multistep methods, actual time stepping is being done
463  * to evaluate the necessary parameters at multiple time
464  * levels needed to start the actual integration.
465  *
466  * \param timestep The size of the timestep, i.e. \f$\Delta t\f$.
467  * \param time on input: the initial time, i.e. \f$t_0\f$.
468  * \param time on output: the new time-level after initialisation, in general this yields
469  * \f$t = t_0 + (r-1)\Delta t\f$ where \f$r\f$ is the number of steps of the multi-step method.
470  * \param nsteps on output: he number of initialisation steps required. In general this corresponds to \f$r-1\f$
471  * where \f$r\f$ is the number of steps of the multi-step method.
472  * \param f an object of the class FuncType, where FuncType should have a method FuncType::ODEforcing
473  * to evaluate the right hand side \f$f(t,\boldsymbol{y})\f$ of the ODE.
474  * \param y_0 the initial value \f$\boldsymbol{y}(t_0)\f$
475  * \return An object of the class TimeIntegrationSolution which represents the vectors \f$\boldsymbol{y}^{[n]}\f$ and \f$t^{[n]}\f$
476  * that can be used to start the actual integration.
477  */
478  LIB_UTILITIES_EXPORT TimeIntegrationSolutionSharedPtr InitializeScheme(const NekDouble timestep,
479  ConstDoubleArray &y_0 ,
480  const NekDouble time ,
481  const TimeIntegrationSchemeOperators &op );
482 
483  /**
484  * \brief Explicit integration of an ODE.
485  *
486  * This function explicitely perfroms a signle integration step of the ODE system:
487  * \f[
488  * \frac{d\boldsymbol{y}}{dt}=\boldsymbol{f}(t,\boldsymbol{y})
489  * \f]
490  *
491  * \param timestep The size of the timestep, i.e. \f$\Delta t\f$.
492  * \param f an object of the class FuncType, where FuncType should have a method FuncType::ODEforcing
493  * to evaluate the right hand side \f$f(t,\boldsymbol{y})\f$ of the ODE.
494  * \param y on input: the vectors \f$\boldsymbol{y}^{[n-1]}\f$ and \f$t^{[n-1]}\f$ (which corresponds to the
495  * solution at the old time level)
496  * \param y on output: the vectors \f$\boldsymbol{y}^{[n]}\f$ and \f$t^{[n]}\f$ (which corresponds to the
497  * solution at the old new level)
498  * \return The actual solution \f$\boldsymbol{y}^{n}\f$ at the new time level
499  * (which in fact is also embedded in the argument y).
500  */
501  LIB_UTILITIES_EXPORT ConstDoubleArray& TimeIntegrate(const NekDouble timestep,
502  TimeIntegrationSolutionSharedPtr &y ,
503  const TimeIntegrationSchemeOperators &op );
504 
506  {
507  return m_timeLevelOffset;
508  }
509 
510 
511  protected:
512  TimeIntegrationSchemeKey m_schemeKey;
514  unsigned int m_numsteps; //< Number of steps in multi-step component.
515  unsigned int m_numstages; //< Number of stages in multi-stage component.
516 
517  bool m_firstStageEqualsOldSolution; //< Optimisation-flag
518  bool m_lastStageEqualsNewSolution; //< Optimisation-flag
519 
520  unsigned int m_numMultiStepValues; // number of entries in input and output vector that correspond
521  // to VALUES at previous time levels
522  unsigned int m_numMultiStepDerivs; // number of entries in input and output vector that correspond
523  // to DERIVATIVES at previous time levels
524  Array<OneD,unsigned int> m_timeLevelOffset; // denotes to which time-level the entries in both
525  // input and output vector correspond, e.g.
526  // INPUT VECTOR --------> m_inputTimeLevelOffset
527  // _ _ _ _
528  // | u^n | | 0 |
529  // | u^{n-1} | | 1 |
530  // | u^{n-2} | -----> | 2 |
531  // | dt f(u^{n-1})| | 1 |
532  // | dt f(u^{n-2})| | 2 |
533  // - - - -
534 
539 
540  private:
541  bool m_initialised; /// bool to identify if array has been initialised
542  int m_nvar; /// The number of variables in integration scheme.
543  int m_npoints; /// The size of inner data which is stored for reuse.
544  DoubleArray m_Y; /// Array containing the stage values
545  DoubleArray m_tmp; /// explicit right hand side of each stage equation
546 
547  TripleArray m_F; /// Array corresponding to the stage Derivatives
548  TripleArray m_F_IMEX; /// Used to store the Explicit stage derivative of IMEX schemes
549 
550  NekDouble m_T; /// Time at the different stages
551 
552  template <typename> friend class Nektar::MemoryManager;
553  LIB_UTILITIES_EXPORT friend TimeIntegrationSchemeManagerT &TimeIntegrationSchemeManager(void);
554 
555  LIB_UTILITIES_EXPORT static boost::shared_ptr<TimeIntegrationScheme> Create(const TimeIntegrationSchemeKey &key);
556 
557  TimeIntegrationScheme(const TimeIntegrationSchemeKey &key);
558 
559  // These should never be called
561  {
562  NEKERROR(ErrorUtil::efatal,"Default Constructor for the TimeIntegrationScheme class should not be called");
563  }
564 
565  TimeIntegrationScheme(const TimeIntegrationScheme &in):m_schemeKey(NullTimeIntegrationSchemeKey)
566  {
567  NEKERROR(ErrorUtil::efatal,"Copy Constructor for the TimeIntegrationScheme class should not be called");
568  }
569 
571  const Array<OneD, const Array<TwoD, NekDouble> >& A,
572  const Array<OneD, const Array<TwoD, NekDouble> >& B,
574  const Array<TwoD, const NekDouble>& V) const;
575 
576  LIB_UTILITIES_EXPORT void TimeIntegrate(const NekDouble timestep,
577  ConstTripleArray &y_old ,
578  ConstSingleArray &t_old ,
579  TripleArray &y_new ,
580  SingleArray &t_new ,
581  const TimeIntegrationSchemeOperators &op );
582 
583 
584  inline int GetFirstDim(ConstTripleArray &y) const
585  {
586  return y[0].num_elements();
587  }
588 
589  inline int GetSecondDim(ConstTripleArray &y) const
590  {
591  return y[0][0].num_elements();
592  }
593 
595  ConstTripleArray &y_old ,
596  ConstSingleArray &t_old ,
597  TripleArray &y_new ,
598  SingleArray &t_new ,
599  const TimeIntegrationSchemeOperators &op) const;
600 
602  const Array<OneD, const Array<TwoD, NekDouble> >& B,
604  const Array<TwoD, const NekDouble>& V) const;
605 
607  const Array<OneD, const Array<TwoD, NekDouble> >& B,
609  const Array<TwoD, const NekDouble>& V) const;
610 
611 
612 
613  };
614 
615 
616  // =========================================================================
617  // ==== DEFINITION OF THE CLASS TimeIntegrationSolution
618  // =========================================================================
620  {
621  public:
624 
625  // Constructor for single step methods
626  LIB_UTILITIES_EXPORT TimeIntegrationSolution(const TimeIntegrationSchemeKey &key,
627  const DoubleArray& y,
628  const NekDouble time,
629  const NekDouble timestep);
630  // Constructor for multi-step methods
631  LIB_UTILITIES_EXPORT TimeIntegrationSolution(const TimeIntegrationSchemeKey &key,
632  const TripleArray& y,
633  const Array<OneD, NekDouble>& t);
634  LIB_UTILITIES_EXPORT TimeIntegrationSolution(const TimeIntegrationSchemeKey &key,
635  unsigned int nvar,
636  unsigned int npoints);
637  LIB_UTILITIES_EXPORT TimeIntegrationSolution(const TimeIntegrationSchemeKey &key);
638 
639  inline const TimeIntegrationSchemeSharedPtr& GetIntegrationScheme() const
640  {
641  return m_scheme;
642  }
643 
644  inline const TimeIntegrationSchemeKey& GetIntegrationSchemeKey() const
645  {
646  return m_scheme->GetIntegrationSchemeKey();
647  }
648 
650  {
651  return m_scheme->GetIntegrationMethod();
652  }
653 
654  inline const TripleArray& GetSolutionVector() const
655  {
656  return m_solVector;
657  }
658 
659  inline TripleArray& UpdateSolutionVector()
660  {
661  return m_solVector;
662  }
663 
664  inline const DoubleArray& GetSolution() const
665  {
666  return m_solVector[0];
667  }
668 
669  inline DoubleArray& UpdateSolution()
670  {
671  return m_solVector[0];
672  }
673 
675  {
676  return m_t;
677  }
678 
680  {
681  return m_t;
682  }
683 
684  inline NekDouble GetTime() const
685  {
686  return m_t[0];
687  }
688 
689  inline int GetNsteps()
690  {
691  return m_scheme->GetNsteps();
692  }
693 
694  inline int GetFirstDim() const
695  {
696  return m_solVector[0].num_elements();
697  }
698 
699  inline int GetSecondDim() const
700  {
701  return m_solVector[0][0].num_elements();
702  }
703 
704  // Return the number of entries in the solution vector that correspond
705  // to (multi-step) values
706  inline unsigned int GetNvalues(void) const
707  {
708  return m_scheme->GetNmultiStepValues();
709  }
710 
711  // Return the number of entries in the solution vector that correspond
712  // to (multi-step) derivatives
713  inline unsigned int GetNderivs(void) const
714  {
715  return m_scheme->GetNmultiStepDerivs();
716  }
717 
718  // Returns an array which indicates to which time-level the entries in the
719  // solution vector correspond
721  {
722  return m_scheme->GetTimeLevelOffset();
723  }
724 
725  // returns the entry in the solution vector which corresponds to
726  // the (multi-step) value at the time-level with specified offset
727  inline DoubleArray& GetValue(const unsigned int timeLevelOffset)
728  {
729  int nMultiStepVals = m_scheme->GetNmultiStepValues();
731 
732  for(int i = 0; i < nMultiStepVals; i++)
733  {
734  if( timeLevelOffset == offsetvec[i] )
735  {
736  return m_solVector[i];
737  }
738  }
739  ASSERTL1(false,"The solution vector of this scheme does not contain a value at the requested time-level");
740  return m_solVector[0];
741  }
742 
743  // returns the entry in the solution vector which corresponds to
744  // the (multi-step) derivative at the time-level with specified offset
745  inline DoubleArray& GetDerivative(const unsigned int timeLevelOffset)
746  {
747  int nMultiStepVals = m_scheme->GetNmultiStepValues();
748  int size = m_scheme->GetNsteps();
750 
751  for(int i = nMultiStepVals; i < size; i++)
752  {
753  if( timeLevelOffset == offsetvec[i] )
754  {
755  return m_solVector[i];
756  }
757  }
758  ASSERTL1(false,"The solution vector of this scheme does not contain a derivative at the requested time-level");
759  return m_solVector[0];
760  }
761 
762  // returns the time associated with the (multi-step) value at the time-level with the
763  // given offset
764  inline NekDouble GetValueTime(const unsigned int timeLevelOffset)
765  {
766  int nMultiStepVals = m_scheme->GetNmultiStepValues();
768 
769  for(int i = 0; i < nMultiStepVals; i++)
770  {
771  if( timeLevelOffset == offsetvec[i] )
772  {
773  return m_t[i];
774  }
775  }
776  ASSERTL1(false,"The solution vector of this scheme does not contain a value at the requested time-level");
777  return m_t[0];
778  }
779 
780  // sets the (multi-step) value and time in the solution
781  // vector which corresponds to
782  // the value at the time-level with specified offset
783  inline void SetValue(const unsigned int timeLevelOffset, const DoubleArray& y, const NekDouble t)
784  {
785  int nMultiStepVals = m_scheme->GetNmultiStepValues();
787 
788  for(int i = 0; i < nMultiStepVals; i++)
789  {
790  if( timeLevelOffset == offsetvec[i] )
791  {
792  m_solVector[i] = y;
793  m_t[i] = t;
794  return;
795  }
796  }
797  }
798 
799  // sets the (multi-step) derivative and time in the
800  // solution vector which corresponds to
801  // the derivative at the time-level with specified offset
802  inline void SetDerivative(const unsigned int timeLevelOffset, const DoubleArray& y, const NekDouble timestep)
803  {
804  int nMultiStepVals = m_scheme->GetNmultiStepValues();
805  int size = m_scheme->GetNsteps();
807 
808  for(int i = nMultiStepVals; i < size; i++)
809  {
810  if( timeLevelOffset == offsetvec[i] )
811  {
812  m_solVector[i] = y;
813  m_t[i] = timestep;
814  return;
815  }
816  }
817  }
818 
819  // sets the soln Vector
820  inline void SetSolVector(const int Offset, const DoubleArray& y)
821  {
822  m_solVector[Offset] = y;
823  }
824 
825  // Rotate the solution vector
826  // (i.e. updating without calculating/inserting new values)
827  inline void RotateSolutionVector()
828  {
829  int i;
830  int nMultiStepVals = m_scheme->GetNmultiStepValues();
831  int size = m_scheme->GetNsteps();
832  for(i = (nMultiStepVals-1); i > 0; i--)
833  {
834  m_solVector[i] = m_solVector[i-1];
835  }
836 
837  for(i = (size-1); i > nMultiStepVals; i--)
838  {
839  m_solVector[i] = m_solVector[i-1];
840  }
841  }
842 
843 
844  private:
845  TimeIntegrationSchemeSharedPtr m_scheme;
846  TripleArray m_solVector;
848  };
849  // =========================================================================
850 
851 
852  LIB_UTILITIES_EXPORT std::ostream& operator<<(std::ostream& os, const TimeIntegrationScheme& rhs);
853  LIB_UTILITIES_EXPORT std::ostream& operator<<(std::ostream& os, const TimeIntegrationSchemeSharedPtr& rhs);
854 
855  // =========================================================================
856 
857  }; // end of namespace
858 } // end of namespace
859 
860 #endif //NEKTAR_LIB_UTILITIES_FOUNDATIONS_TIMEINTEGRATIONSCHEME_H
NekDouble GetValueTime(const unsigned int timeLevelOffset)
TimeIntegrationSchemeKey(const TimeIntegrationMethod &method)
std::vector< TimeIntegrationSchemeSharedPtr > TimeIntegrationSchemeVector
boost::function< void(ConstDoubleArray &, DoubleArray &, const NekDouble, const NekDouble) > FunctorType2
NekDouble B_IMEX(const unsigned int i, const unsigned int j) const
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:185
NekDouble U(const unsigned int i, const unsigned int j) const
NekDouble B(const unsigned int i, const unsigned int j) const
TimeIntegrationSchemeKey & operator=(const TimeIntegrationSchemeKey &key)
BDF multi-step scheme of order 1 (implicit)
Adams-Bashforth Forward multi-step scheme of order 2.
Array< OneD, Array< TwoD, NekDouble > > m_B
static const TimeIntegrationSchemeKey NullTimeIntegrationSchemeKey(eNoTimeIntegrationMethod)
int m_nvar
bool to identify if array has been initialised
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
Runge-Kutta multi-stage scheme 4th order explicit (old name)
Implicit-Explicit Midpoint IMEX DIRK(1,2,2)
void DoImplicitSolve(InArrayType &inarray, OutArrayType &outarray, const NekDouble time, const NekDouble lambda) const
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
const char *const TimeIntegrationSchemeTypeMap[]
DoubleArray m_tmp
Array containing the stage values.
TimeIntegrationMethod GetIntegrationMethod() const
Classical RungeKutta4 method (new name for eClassicalRungeKutta4)
void DefineOdeExplicitRhs(FuncPointerT func, ObjectPointerT obj)
NekDouble V(const unsigned int i, const unsigned int j) const
TimeIntegrationSolutionSharedPtr InitializeScheme(const NekDouble timestep, ConstDoubleArray &y_0, const NekDouble time, const TimeIntegrationSchemeOperators &op)
This function initialises the time integration scheme.
L-stable, four stage, third order IMEX DIRK(4,4,3)
Forward-Backward Euler IMEX DIRK(1,2,1)
void DefineOdeImplicitRhs(FuncPointerT func, ObjectPointerT obj)
bool operator==(const TimeIntegrationSchemeKey &key)
Implicit Explicit General Linear Method.
boost::function< void(InArrayType &, OutArrayType &, const NekDouble, const NekDouble) > FunctorType2
const char *const TimeIntegrationMethodMap[]
void DoOdeRhs(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
NekManager< TimeIntegrationSchemeKey, TimeIntegrationScheme, TimeIntegrationSchemeKey::opLess > TimeIntegrationSchemeManagerT
const TimeIntegrationSchemeSharedPtr & GetIntegrationScheme() const
bool operator!=(const TimeIntegrationSchemeKey &y)
Nonlinear SSP RungeKutta3 explicit.
bool operator()(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs) const
std::vector< TimeIntegrationSolutionSharedPtr >::iterator TimeIntegrationSolutionVectorIter
DoubleArray m_Y
The size of inner data which is stored for reuse.
TimeIntegrationSchemeManagerT & TimeIntegrationSchemeManager(void)
Adams-Moulton Forward multi-step scheme of order 2.
std::vector< TimeIntegrationSchemeSharedPtr >::iterator TimeIntegrationSchemeVectorIter
Adams-Bashforth Forward multi-step scheme of order 3.
StandardMatrixTag & lhs
bool operator==(const BasisKey &x, const BasisKey &y)
Definition: Basis.cpp:703
void SetValue(const unsigned int timeLevelOffset, const DoubleArray &y, const NekDouble t)
Crank-Nicolson/Adams-Bashforth Order 2 (CNAB)
int GetSecondDim(ConstTripleArray &y) const
const Array< OneD, const Array< OneD, Array< OneD, NekDouble > > > ConstTripleArray
DoubleArray & GetDerivative(const unsigned int timeLevelOffset)
boost::function< void(InArrayType &, OutArrayType &, const NekDouble) > FunctorType1
IMEX 2nd order scheme using Backward Different Formula & Extrapolation.
std::ostream & operator<<(std::ostream &os, const BasisKey &rhs)
Definition: Basis.cpp:85
TripleArray m_F
explicit right hand side of each stage equation
Diagonally implicit scheme (e.g. the DIRK schemes)
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
Classical RungeKutta2 method (new name for eMidpoint)
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
const Array< OneD, const Array< OneD, NekDouble > > InArrayType
Adams-Moulton Forward multi-step scheme of order 1.
NekDouble A_IMEX(const unsigned int i, const unsigned int j) const
void SetDerivative(const unsigned int timeLevelOffset, const DoubleArray &y, const NekDouble timestep)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Adams-Bashforth Forward multi-step scheme of order 1.
TimeIntegrationSolution(const TimeIntegrationSchemeKey &key, const DoubleArray &y, const NekDouble time, const NekDouble timestep)
const TimeIntegrationSchemeKey & GetIntegrationSchemeKey() const
#define LIB_UTILITIES_EXPORT
Array< OneD, Array< OneD, NekDouble > > DoubleArray
L-stable, three stage, third order IMEX DIRK(3,4,3)
friend TimeIntegrationSchemeManagerT & TimeIntegrationSchemeManager(void)
TimeIntegrationSchemeKey(const TimeIntegrationSchemeKey &key)
int m_npoints
The number of variables in integration scheme.
Nonlinear SSP RungeKutta2 explicit (surrogate for eRungeKutta2_ImprovedEuler)
ConstDoubleArray & TimeIntegrate(const NekDouble timestep, TimeIntegrationSolutionSharedPtr &y, const TimeIntegrationSchemeOperators &op)
Explicit integration of an ODE.
bool CheckTimeIntegrateArguments(const NekDouble timestep, ConstTripleArray &y_old, ConstSingleArray &t_old, TripleArray &y_new, SingleArray &t_new, const TimeIntegrationSchemeOperators &op) const
Improved RungeKutta2 explicit (old name meaning Heun's method)
Forward-Backward Euler IMEX DIRK(1,1,1)
Array< OneD, Array< OneD, NekDouble > > DoubleArray
NekDouble A(const unsigned int i, const unsigned int j) const
double NekDouble
DoubleArray & GetValue(const unsigned int timeLevelOffset)
void SetSolVector(const int Offset, const DoubleArray &y)
TimeIntegrationMethod GetIntegrationMethod() const
void DoOdeImplicitRhs(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
void DoProjection(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
bool operator<(const BasisKey &lhs, const BasisKey &rhs)
Definition: Basis.cpp:49
TimeIntegrationSchemeType GetIntegrationSchemeType() const
boost::shared_ptr< TimeIntegrationScheme > TimeIntegrationSchemeSharedPtr
static boost::shared_ptr< TimeIntegrationScheme > Create(const TimeIntegrationSchemeKey &key)
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > TripleArray
TimeIntegrationScheme(const TimeIntegrationScheme &in)
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
BDF multi-step scheme of order 2 (implicit)
L-stable, two stage, second order IMEX DIRK(2,2,2)
TimeIntegrationMethod m_method
integration method
void DoOdeExplicitRhs(InArrayType &inarray, OutArrayType &outarray, const NekDouble time) const
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
Array< OneD, Array< OneD, NekDouble > > OutArrayType
Diagonally Implicit Runge Kutta scheme of order 3.
IMEX 3rd order scheme using Backward Different Formula & Extrapolation.
const Array< OneD, const NekDouble > ConstSingleArray
L-stable, three stage, third order IMEX DIRK(3,4,3)
Diagonally Implicit Runge Kutta scheme of order 3.
boost::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
std::vector< TimeIntegrationSolutionSharedPtr > TimeIntegrationSolutionVector
const Array< OneD, const unsigned int > & GetTimeLevelOffset()
L-stable, two stage, third order IMEX DIRK(2,3,3)
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
boost::function< void(ConstDoubleArray &, DoubleArray &, const NekDouble) > FunctorType1
Array< OneD, Array< TwoD, NekDouble > > m_A
const Array< OneD, const NekDouble > & GetTimeVector() const
Modified Crank-Nicolson/Adams-Bashforth Order 2 (MCNAB)
IMEX 1st order scheme using Euler Backwards/Euler Forwards.
midpoint method (old name)
const Array< OneD, const unsigned int > & GetTimeLevelOffset()
NekDouble m_T
Used to store the Explicit stage derivative of IMEX schemes.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > TripleArray
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:218
const TimeIntegrationSchemeKey & GetIntegrationSchemeKey() const
TripleArray m_F_IMEX
Array corresponding to the stage Derivatives.
const Array< OneD, const Array< OneD, NekDouble > > ConstDoubleArray
friend bool operator<(const TimeIntegrationSchemeKey &lhs, const TimeIntegrationSchemeKey &rhs)