OGS 6.2.1-76-gbb689931b
NumLib::TimeDiscretization Class Referenceabstract

Detailed Description

Interface of time discretization schemes for first-order ODEs.

The purpose of TimeDiscretization instances is to store the solution history of an ODE, i. e., to keep the solution at as many timestamps as is required by the respective time discretization scheme. Furthermore, TimeDiscretization instances compute the discretized approximation of the time derivative $ \partial x/\partial t $.

The method documentation of this class uses quantities introduced in the following section.
Currently this interface does not yet support adaptive timestepping. While implementing that will lead to no changes for single-step methods, for multi-step methods this interface will have to be extended.

Discretizing first-order ODEs

A first-order (implicit) ODE has the general form

\[ F(\dot x, x, t) \stackrel{!}{=} 0. \]

In order to solve it numerically a certain time discretization scheme, such as the forward or backward Euler methods, is used. The discretized ODE is then given by

\[ F(\hat x, x_C, t_C) \stackrel{!}{=} 0. \]

This interface has been designed with first-order implicit quasi-linear ODEs in mind. They can be expressed as follows and are given here only to serve as an example.

\[ M(x,t)\cdot \dot x + K(x,t) \cdot x - b(x,t) =: r(\dot x, x, t) \stackrel{!}{=} 0. \]

After time discretization this formula becomes:

\[ M(x_C,t_C)\cdot \hat x + K(x_C,t_C) \cdot x_C - b(x_C,t_C) =: r(\hat x, x_C, t_C) \stackrel{!}{=} 0. \]

The meaning of indices for $ x $ and $ t $ is as follows:

  • $ C $ – "Current": The values of $ x $ and $ t $ at which the discretized ODE is being assembled.
  • $ N $ – "New": The values of $ x $ and $ t $ at the new timestep that is being calculated right now by the ODE solver.
  • $ O $ – "Old": The results from the preceding time step (or a linear combination of results of the preceding time steps in the case of multistep methods) weighted by a scalar factor.
  • $ n $ – Numerical index indicating the timestep.

$ \hat x $ is the discrete approximation of $ \dot x := \partial x/\partial t$. It is assumed that $ \hat x $ can be written in the following form:

\[ \hat x = \alpha \cdot x_N - x_O, \]

where $ \alpha := \partial \hat x / \partial x_N $ is a scalar.

For different time discretization schemes $ x_C $, $ t_C $, $ x_N $, $ x_O $ and $ \alpha $ take different values. Those for the time implemented schemes are given in the table below.

Scheme $ x_C $ $ t_C $ $ \alpha $ $ x_N $ $ x_O $
Forward Euler $ x_n $ $ t_n $ $ 1/\Delta t $ $ x_{n+1} $ $ x_n / \Delta t $
Backward Euler $ x_{n+1} $ $ t_{n+1} $ $ 1/\Delta t $ $ x_{n+1} $ $ x_n / \Delta t $
Crank-Nicolson $ x_{n+1} $ $ t_{n+1} $ $ 1/\Delta t $ $ x_{n+1} $ $ x_n / \Delta t $
BDF(2) $ x_{n+2} $ $ t_{n+1} $ $ 3/(2\Delta t) $ $ x_{n+2} $ $ (2\cdot x_{n+1} - x_n/2)/\Delta t $

The other backward differentiation formulas of orders 1 to 6 are also implemented, but only BDF(2) has bee given here for brevity.

Definition at line 114 of file TimeDiscretization.h.

#include <TimeDiscretization.h>

Inheritance diagram for NumLib::TimeDiscretization:

Public Member Functions

 TimeDiscretization ()=default
virtual void setInitialState (const double t0, GlobalVector const &x0)=0
 Sets the initial condition. More...
virtual double getRelativeChangeFromPreviousTimestep (GlobalVector const &x, MathLib::VecNormType norm_type)=0
virtual void pushState (const double t, GlobalVector const &x, InternalMatrixStorage const &strg)=0
virtual void popState (GlobalVector &x)=0
virtual void nextTimestep (const double t, const double delta_t)=0
virtual double getCurrentTime () const =0
void getXdot (GlobalVector const &x_at_new_timestep, GlobalVector &xdot) const
virtual double getNewXWeight () const =0
 Returns $ \alpha = \partial \hat x / \partial x_N $. More...
virtual void getWeightedOldX (GlobalVector &y) const =0
 Returns $ x_O $. More...
virtual ~TimeDiscretization ()=default
Extended Interface

These methods are provided primarily to make certain concrete time discretizations with special demands, such as the forward Euler or Crank-Nicolson schemes, possible.

virtual bool isLinearTimeDisc () const
virtual double getDxDx () const
virtual GlobalVector const & getCurrentX (GlobalVector const &x_at_new_timestep) const
virtual bool needsPreload () const

Protected Member Functions

double computeRelativeChangeFromPreviousTimestep (GlobalVector const &x, GlobalVector const &x_old, MathLib::VecNormType norm_type)

Protected Attributes

std::unique_ptr< GlobalVector > _dx
 Used to store $ u_{n+1}-u_{n}$. More...

Constructor & Destructor Documentation

◆ TimeDiscretization()

NumLib::TimeDiscretization::TimeDiscretization ( )

◆ ~TimeDiscretization()

virtual NumLib::TimeDiscretization::~TimeDiscretization ( )

Member Function Documentation

◆ computeRelativeChangeFromPreviousTimestep()

double NumLib::TimeDiscretization::computeRelativeChangeFromPreviousTimestep ( GlobalVector const &  x,
GlobalVector const &  x_old,
MathLib::VecNormType  norm_type 

Compute and return the relative change of solutions between two successive time steps by $ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| $.

xThe current solution
x_oldThe previous solution
norm_typeThe norm type of global vector
$ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| $.
the value of x_old is changed to x - x_old after this computation.

Definition at line 18 of file TimeDiscretization.cpp.

References _dx, MathLib::LinAlg::axpy(), MathLib::LinAlg::copy(), MathLib::INVALID, MathLib::LinAlg::norm(), and OGS_FATAL.

Referenced by NumLib::BackwardEuler::getRelativeChangeFromPreviousTimestep(), NumLib::ForwardEuler::getRelativeChangeFromPreviousTimestep(), NumLib::CrankNicolson::getRelativeChangeFromPreviousTimestep(), and NumLib::BackwardDifferentiationFormula::getRelativeChangeFromPreviousTimestep().

22 {
23  if (norm_type == MathLib::VecNormType::INVALID)
24  {
25  OGS_FATAL("An invalid norm type has been passed");
26  }
28  if (!_dx)
29  {
31  }
33  auto& dx = *_dx;
34  MathLib::LinAlg::copy(x, dx); // copy x to dx.
36  // dx = x - x_old --> x - dx --> dx
37  MathLib::LinAlg::axpy(dx, -1.0, x_old);
38  const double norm_dx = MathLib::LinAlg::norm(dx, norm_type);
40  const double norm_x = MathLib::LinAlg::norm(x, norm_type);
41  if (norm_x > std::numeric_limits<double>::epsilon())
42  {
43  return norm_dx / norm_x;
44  }
46  // Both of norm_x and norm_dx are close to zero
47  if (norm_dx < std::numeric_limits<double>::epsilon())
48  {
49  return 1.0;
50  }
52  // Only norm_x is close to zero
53  return norm_dx / std::numeric_limits<double>::epsilon();
54 }
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition: LinAlg.h:87
std::unique_ptr< GlobalVector > _dx
Used to store .
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:57
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36

◆ getCurrentTime()

virtual double NumLib::TimeDiscretization::getCurrentTime ( ) const
pure virtual

Returns $ t_C $, i.e., the time at which the equation will be assembled.

Implemented in NumLib::BackwardDifferentiationFormula, NumLib::CrankNicolson, NumLib::ForwardEuler, and NumLib::BackwardEuler.

◆ getCurrentX()

virtual GlobalVector const& NumLib::TimeDiscretization::getCurrentX ( GlobalVector const &  x_at_new_timestep) const

Returns $ x_C $, i.e., the state at which the equation will be assembled.

This method is overridden in the ForwardEuler scheme.

Reimplemented in NumLib::ForwardEuler.

Definition at line 213 of file TimeDiscretization.h.

214  {
215  return x_at_new_timestep;
216  }

◆ getDxDx()

virtual double NumLib::TimeDiscretization::getDxDx ( ) const

Returns $ \partial x_C / \partial x_N $.

The ForwardEuler scheme overrides this.

Reimplemented in NumLib::ForwardEuler.

Definition at line 207 of file TimeDiscretization.h.

207 { return 1.0; }

◆ getNewXWeight()

virtual double NumLib::TimeDiscretization::getNewXWeight ( ) const
pure virtual

◆ getRelativeChangeFromPreviousTimestep()

virtual double NumLib::TimeDiscretization::getRelativeChangeFromPreviousTimestep ( GlobalVector const &  x,
MathLib::VecNormType  norm_type 
pure virtual

Get the relative change of solutions between two successive time steps by $ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| $.

xThe solution at the current timestep.
norm_typeThe type of global vector norm.

Implemented in NumLib::BackwardDifferentiationFormula, NumLib::CrankNicolson, NumLib::ForwardEuler, and NumLib::BackwardEuler.

◆ getWeightedOldX()

virtual void NumLib::TimeDiscretization::getWeightedOldX ( GlobalVector &  y) const
pure virtual

◆ getXdot()

void NumLib::TimeDiscretization::getXdot ( GlobalVector const &  x_at_new_timestep,
GlobalVector &  xdot 
) const

Returns $ \hat x $, i.e. the discretized approximation of $ \dot x //! $.

Definition at line 171 of file TimeDiscretization.h.

References MathLib::LinAlg::axpby().

172  {
173  namespace LinAlg = MathLib::LinAlg;
175  auto const dxdot_dx = getNewXWeight();
177  // xdot = dxdot_dx * x_at_new_timestep - x_old
178  getWeightedOldX(xdot);
179  LinAlg::axpby(xdot, dxdot_dx, -1.0, x_at_new_timestep);
180  }
virtual double getNewXWeight() const =0
Returns .
virtual void getWeightedOldX(GlobalVector &y) const =0
Returns .
void axpby(MatrixOrVector &y, double const a, double const b, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:64

◆ isLinearTimeDisc()

virtual bool NumLib::TimeDiscretization::isLinearTimeDisc ( ) const

Tell whether this scheme inherently requires a nonlinear solver or not.

The ForwardEuler scheme is inherently linear in that sense, the others are not.

Reimplemented in NumLib::ForwardEuler.

Definition at line 202 of file TimeDiscretization.h.

202 { return false; }

◆ needsPreload()

virtual bool NumLib::TimeDiscretization::needsPreload ( ) const

Indicate that this scheme needs some additional assembly before the first timestep will be solved.

The CrankNicolson scheme needs such preload.

Reimplemented in NumLib::CrankNicolson.

Definition at line 224 of file TimeDiscretization.h.

224 { return false; }

◆ nextTimestep()

virtual void NumLib::TimeDiscretization::nextTimestep ( const double  t,
const double  delta_t 
pure virtual

Indicate that the computation of a new timestep is being started now.

Currently changing timestep sizes are not supported. Thus, delta_t must not change throughout the entire time integration process! This is not checked by this code!

Implemented in NumLib::BackwardDifferentiationFormula, NumLib::CrankNicolson, NumLib::ForwardEuler, and NumLib::BackwardEuler.

◆ popState()

virtual void NumLib::TimeDiscretization::popState ( GlobalVector &  x)
pure virtual

Restores the given vector x to its old value. Used only for repeating of the time step with new time step size when the current time step is rejected. The restored x is only used as an initial guess for linear solver in the first Picard nonlinear iteration.

xThe solution at the current time step, which is going to be reset to its previous value.

Implemented in NumLib::BackwardDifferentiationFormula, NumLib::CrankNicolson, NumLib::ForwardEuler, and NumLib::BackwardEuler.

◆ pushState()

virtual void NumLib::TimeDiscretization::pushState ( const double  t,
GlobalVector const &  x,
InternalMatrixStorage const &  strg 
pure virtual

Indicate that the current timestep is done and that you will proceed to the next one.

Do not use this method for setting the initial condition, rather use setInitialState()!
tThe current timestep.
xThe solution at the current timestep.
strgTrigger storing some internal state. Currently only used by the CrankNicolson scheme.

Implemented in NumLib::BackwardDifferentiationFormula, NumLib::CrankNicolson, NumLib::ForwardEuler, and NumLib::BackwardEuler.

◆ setInitialState()

virtual void NumLib::TimeDiscretization::setInitialState ( const double  t0,
GlobalVector const &  x0 
pure virtual

Member Data Documentation

◆ _dx

std::unique_ptr<GlobalVector> NumLib::TimeDiscretization::_dx

Used to store $ u_{n+1}-u_{n}$.

Definition at line 228 of file TimeDiscretization.h.

Referenced by computeRelativeChangeFromPreviousTimestep().

The documentation for this class was generated from the following files: