OGS 6.2.2-87-g988ee9c30.dirty.20200123122242
TimeDiscretization.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <vector>
14 
15 #include "MathLib/LinAlg/LinAlg.h"
17 #include "Types.h"
18 
19 namespace NumLib
20 {
23 
27 {
28 public:
41  virtual void pushMatrices() const = 0;
42 
43  virtual ~InternalMatrixStorage() = default;
44 };
45 
115 {
116 public:
117  TimeDiscretization() = default;
118 
120  virtual void setInitialState(const double t0, GlobalVector const& x0) = 0;
121 
128  virtual double getRelativeChangeFromPreviousTimestep(
129  GlobalVector const& x, MathLib::VecNormType norm_type) = 0;
130 
142  virtual void pushState(const double t, GlobalVector const& x,
143  InternalMatrixStorage const& strg) = 0;
144 
155  virtual void popState(GlobalVector& x) = 0;
156 
163  virtual void nextTimestep(const double t, const double delta_t) = 0;
164 
167  virtual double getCurrentTime() const = 0;
168 
171  virtual double getCurrentTimeIncrement() const = 0;
172 
175  void getXdot(GlobalVector const& x_at_new_timestep, GlobalVector& xdot) const
176  {
177  namespace LinAlg = MathLib::LinAlg;
178 
179  auto const dxdot_dx = getNewXWeight();
180 
181  // xdot = dxdot_dx * x_at_new_timestep - x_old
182  getWeightedOldX(xdot);
183  LinAlg::axpby(xdot, dxdot_dx, -1.0, x_at_new_timestep);
184  }
185 
187  virtual double getNewXWeight() const = 0;
188 
190  virtual void getWeightedOldX(GlobalVector& y) const = 0; // = x_old
191 
192  virtual ~TimeDiscretization() = default;
193 
200 
206  virtual bool isLinearTimeDisc() const { return false; }
211  virtual double getDxDx() const { return 1.0; }
217  virtual GlobalVector const& getCurrentX(GlobalVector const& x_at_new_timestep) const
218  {
219  return x_at_new_timestep;
220  }
221 
228  virtual bool needsPreload() const { return false; }
230 
231 protected:
232  std::unique_ptr<GlobalVector> _dx;
233 
246  double computeRelativeChangeFromPreviousTimestep(
247  GlobalVector const& x,
248  GlobalVector const& x_old,
249  MathLib::VecNormType norm_type);
250 };
251 
253 class BackwardEuler final : public TimeDiscretization
254 {
255 public:
257  : _x_old(NumLib::GlobalVectorProvider::provider.getVector())
258  {
259  }
260 
261  ~BackwardEuler() override
262  {
264  }
265 
266  void setInitialState(const double t0, GlobalVector const& x0) override
267  {
268  _t = t0;
269  MathLib::LinAlg::copy(x0, _x_old);
270  }
271 
272  double getRelativeChangeFromPreviousTimestep(
273  GlobalVector const& x, MathLib::VecNormType norm_type) override;
274 
275  void pushState(const double /*t*/, GlobalVector const& x,
276  InternalMatrixStorage const& /*strg*/) override
277  {
278  MathLib::LinAlg::copy(x, _x_old);
279  }
280 
281  void popState(GlobalVector& x) override
282  {
283  MathLib::LinAlg::copy(_x_old, x);
284  }
285 
286  void nextTimestep(const double t, const double delta_t) override
287  {
288  _t = t;
289  _delta_t = delta_t;
290  }
291 
292  double getCurrentTime() const override { return _t; }
293  double getCurrentTimeIncrement() const override { return _delta_t; }
294  double getNewXWeight() const override { return 1.0 / _delta_t; }
295  void getWeightedOldX(GlobalVector& y) const override
296  {
297  namespace LinAlg = MathLib::LinAlg;
298 
299  // y = x_old / delta_t
300  LinAlg::copy(_x_old, y);
301  LinAlg::scale(y, 1.0 / _delta_t);
302  }
303 
304 private:
305  double _t = std::numeric_limits<double>::quiet_NaN();
306  double _delta_t =
307  std::numeric_limits<double>::quiet_NaN();
308  GlobalVector& _x_old;
309 };
310 
312 class ForwardEuler final : public TimeDiscretization
313 {
314 public:
315  ForwardEuler() : _x_old(NumLib::GlobalVectorProvider::provider.getVector())
316  {
317  }
318 
319  ~ForwardEuler() override
320  {
322  }
323 
324  void setInitialState(const double t0, GlobalVector const& x0) override
325  {
326  _t = t0;
327  _t_old = t0;
328  MathLib::LinAlg::copy(x0, _x_old);
329  }
330 
331  double getRelativeChangeFromPreviousTimestep(
332  GlobalVector const& x, MathLib::VecNormType norm_type) override;
333 
334  void pushState(const double /*t*/, GlobalVector const& x,
335  InternalMatrixStorage const& /*strg*/) override
336  {
337  MathLib::LinAlg::copy(x, _x_old);
338  }
339 
340  void popState(GlobalVector& x) override
341  {
342  MathLib::LinAlg::copy(_x_old, x);
343  }
344 
345  void nextTimestep(const double t, const double delta_t) override
346  {
347  _t_old = _t;
348  _t = t;
349  _delta_t = delta_t;
350  }
351 
352  double getCurrentTime() const override
353  {
354  return _t_old; // forward Euler does assembly at the preceding timestep
355  }
356 
357  double getCurrentTimeIncrement() const override { return _delta_t; }
358 
359  GlobalVector const& getCurrentX(
360  const GlobalVector& /*x_at_new_timestep*/) const override
361  {
362  return _x_old;
363  }
364 
365  double getNewXWeight() const override { return 1.0 / _delta_t; }
366  void getWeightedOldX(GlobalVector& y) const override
367  {
368  namespace LinAlg = MathLib::LinAlg;
369 
370  // y = x_old / delta_t
371  LinAlg::copy(_x_old, y);
372  LinAlg::scale(y, 1.0 / _delta_t);
373  }
374 
375  bool isLinearTimeDisc() const override { return true; }
376  double getDxDx() const override { return 0.0; }
378  GlobalVector const& getXOld() const { return _x_old; }
379 private:
380  double _t = std::numeric_limits<double>::quiet_NaN();
381  double _t_old =
382  std::numeric_limits<double>::quiet_NaN();
383  double _delta_t =
385  std::numeric_limits<double>::quiet_NaN();
386  GlobalVector& _x_old;
387 };
388 
390 class CrankNicolson final : public TimeDiscretization
391 {
392 public:
401  explicit CrankNicolson(const double theta)
402  : _theta(theta),
403  _x_old(NumLib::GlobalVectorProvider::provider.getVector())
404  {
405  }
406 
407  ~CrankNicolson() override
408  {
410  }
411 
412  void setInitialState(const double t0, GlobalVector const& x0) override
413  {
414  _t = t0;
415  MathLib::LinAlg::copy(x0, _x_old);
416  }
417 
418  double getRelativeChangeFromPreviousTimestep(
419  GlobalVector const& x, MathLib::VecNormType norm_type) override;
420 
421  void pushState(const double /*t*/, GlobalVector const& x,
422  InternalMatrixStorage const& strg) override
423  {
424  MathLib::LinAlg::copy(x, _x_old);
425  strg.pushMatrices();
426  }
427 
428  void popState(GlobalVector& x) override
429  {
430  MathLib::LinAlg::copy(_x_old, x);
431  }
432 
433  void nextTimestep(const double t, const double delta_t) override
434  {
435  _t = t;
436  _delta_t = delta_t;
437  }
438 
439  double getCurrentTime() const override { return _t; }
440  double getCurrentTimeIncrement() const override { return _delta_t; }
441  double getNewXWeight() const override { return 1.0 / _delta_t; }
442  void getWeightedOldX(GlobalVector& y) const override
443  {
444  namespace LinAlg = MathLib::LinAlg;
445 
446  // y = x_old / delta_t
447  LinAlg::copy(_x_old, y);
448  LinAlg::scale(y, 1.0 / _delta_t);
449  }
450 
451  bool needsPreload() const override { return true; }
453  double getTheta() const { return _theta; }
455  GlobalVector const& getXOld() const { return _x_old; }
456 private:
457  const double _theta;
458  double _t = std::numeric_limits<double>::quiet_NaN();
459  double _delta_t =
460  std::numeric_limits<double>::quiet_NaN();
461  GlobalVector& _x_old;
462 };
463 
466 {
467 public:
480  explicit BackwardDifferentiationFormula(const unsigned num_steps)
481  : _num_steps(num_steps)
482  {
483  assert(1 <= num_steps && num_steps <= 6);
484  _xs_old.reserve(num_steps);
485  }
486 
488  {
489  for (auto* x : _xs_old)
490  {
492  }
493  }
494 
495  void setInitialState(const double t0, GlobalVector const& x0) override
496  {
497  _t = t0;
498  _xs_old.push_back(
500  }
501 
502  double getRelativeChangeFromPreviousTimestep(
503  GlobalVector const& x, MathLib::VecNormType norm_type) override;
504 
505  void pushState(const double /*t*/, GlobalVector const& x,
506  InternalMatrixStorage const& /*strg*/) override;
507 
508  void popState(GlobalVector& x) override
509  {
510  MathLib::LinAlg::copy(*_xs_old[_xs_old.size()-1], x);
511  }
512 
513  void nextTimestep(const double t, const double delta_t) override
514  {
515  _t = t;
516  _delta_t = delta_t;
517  }
518 
519  double getCurrentTime() const override { return _t; }
520  double getCurrentTimeIncrement() const override { return _delta_t; }
521 
522  double getNewXWeight() const override;
523 
524  void getWeightedOldX(GlobalVector& y) const override;
525 
526 private:
527  std::size_t eff_num_steps() const { return _xs_old.size(); }
528  const unsigned _num_steps;
529  double _t = std::numeric_limits<double>::quiet_NaN();
530  double _delta_t =
531  std::numeric_limits<double>::quiet_NaN();
532 
533  std::vector<GlobalVector*> _xs_old;
534  unsigned _offset = 0;
535 };
536 
538 } // namespace NumLib
Generalized Crank-Nicolson scheme.
double getCurrentTime() const override
double getCurrentTimeIncrement() const override
double getNewXWeight() const override
Returns .
double getNewXWeight() const override
Returns .
double getNewXWeight() const override
Returns .
double getCurrentTime() const override
virtual bool needsPreload() const
const double _theta
the implicitness parameter
Backward differentiation formula.
double getDxDx() const override
bool isLinearTimeDisc() const override
void popState(GlobalVector &x) override
void nextTimestep(const double t, const double delta_t) override
double getCurrentTime() const override
VecNormType
Norm type. Not declared as class type in order to use the members as integers.
Definition: LinAlgEnums.h:21
GlobalVector const & getXOld() const
Returns the solution from the preceding timestep.
virtual void pushMatrices() const =0
double getCurrentTimeIncrement() const override
bool needsPreload() const override
void pushState(const double, GlobalVector const &x, InternalMatrixStorage const &) override
void setInitialState(const double t0, GlobalVector const &x0) override
Sets the initial condition.
void getWeightedOldX(GlobalVector &y) const override
Returns .
GlobalVector & _x_old
the solution from the preceding timestep
void nextTimestep(const double t, const double delta_t) override
Backward Euler scheme.
double getCurrentTimeIncrement() const override
void pushState(const double, GlobalVector const &x, InternalMatrixStorage const &) override
static NUMLIB_EXPORT VectorProvider & provider
std::vector< GlobalVector * > _xs_old
solutions from the preceding timesteps
double getTheta() const
Returns .
double getCurrentTimeIncrement() const override
void axpby(MatrixOrVector &y, double const a, double const b, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:65
void popState(GlobalVector &x) override
GlobalVector const & getXOld() const
Returns the solution from the preceding timestep.
void setInitialState(const double t0, GlobalVector const &x0) override
Sets the initial condition.
GlobalVector & _x_old
the solution from the preceding timestep
Forward Euler scheme.
CrankNicolson(const double theta)
void nextTimestep(const double t, const double delta_t) override
void setInitialState(const double t0, GlobalVector const &x0) override
Sets the initial condition.
GlobalVector const & getCurrentX(const GlobalVector &) const override
const unsigned _num_steps
The order of the BDF method.
virtual ~InternalMatrixStorage()=default
BackwardDifferentiationFormula(const unsigned num_steps)
void pushState(const double, GlobalVector const &x, InternalMatrixStorage const &strg) override
std::unique_ptr< GlobalVector > _dx
Used to store .
void getWeightedOldX(GlobalVector &y) const override
Returns .
GlobalVector & _x_old
the solution from the preceding timestep
virtual void releaseVector(GlobalVector const &x)=0
void popState(GlobalVector &x) override
virtual bool isLinearTimeDisc() const
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:37
void nextTimestep(const double t, const double delta_t) override
virtual double getDxDx() const
virtual GlobalVector const & getCurrentX(GlobalVector const &x_at_new_timestep) const
void getWeightedOldX(GlobalVector &y) const override
Returns .
void scale(MatrixOrVector &x, double const a)
Scales x by a.
Definition: LinAlg.h:44
void setInitialState(const double t0, GlobalVector const &x0) override
Sets the initial condition.
void getXdot(GlobalVector const &x_at_new_timestep, GlobalVector &xdot) const
void popState(GlobalVector &x) override