 OGS
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$$.

Note
The method documentation of this class uses quantities introduced in the following section.
Todo:
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$$
Backward Euler $$x_{n+1}$$ $$t_{n+1}$$ $$1/\Delta t$$ $$x_{n+1}$$ $$x_n / \Delta t$$

Definition at line 87 of file TimeDiscretization.h.

#include <TimeDiscretization.h>

Inheritance diagram for NumLib::TimeDiscretization: [legend]

## Public Member Functions

TimeDiscretization ()=default

virtual void setInitialState (const double t0)=0
Sets the initial condition. More...

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

virtual void nextTimestep (const double t, const double delta_t)=0

virtual double getCurrentTime () const =0

virtual double getCurrentTimeIncrement () const =0

void getXdot (GlobalVector const &x_at_new_timestep, GlobalVector const &x_old, GlobalVector &xdot) const

virtual double getNewXWeight () const =0
Returns $$\alpha = \partial \hat x / \partial x_N$$. More...

virtual void getWeightedOldX (GlobalVector &y, GlobalVector const &x_old) const =0
Returns $$x_O$$. More...

virtual ~TimeDiscretization ()=default

## Protected Attributes

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

## ◆ TimeDiscretization()

 NumLib::TimeDiscretization::TimeDiscretization ( )
default

## ◆ ~TimeDiscretization()

 virtual NumLib::TimeDiscretization::~TimeDiscretization ( )
virtualdefault

## ◆ 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}\|$$.

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

Definition at line 18 of file TimeDiscretization.cpp.

22 {
23  if (norm_type == MathLib::VecNormType::INVALID)
24  {
25  OGS_FATAL("An invalid norm type has been passed");
26  }
27
28  if (!_dx)
29  {
31  }
32
33  auto& dx = *_dx;
34  MathLib::LinAlg::copy(x, dx); // copy x to dx.
35
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);
39
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  }
45
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  }
51
52  // Only norm_x is close to zero
53  return norm_dx / std::numeric_limits<double>::epsilon();
54 }

## ◆ 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::BackwardEuler.

## ◆ getCurrentTimeIncrement()

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

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

Implemented in NumLib::BackwardEuler.

## ◆ getNewXWeight()

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

Returns $$\alpha = \partial \hat x / \partial x_N$$.

Implemented in NumLib::BackwardEuler.

Referenced by getXdot().

## ◆ getWeightedOldX()

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

Returns $$x_O$$.

Implemented in NumLib::BackwardEuler.

Referenced by getXdot().

## ◆ getXdot()

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

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

Definition at line 130 of file TimeDiscretization.h.

133  {
134  namespace LinAlg = MathLib::LinAlg;
135
136  auto const dxdot_dx = getNewXWeight();
137
138  // xdot = dxdot_dx * x_at_new_timestep - x_old
139  getWeightedOldX(xdot, x_old);
140  LinAlg::axpby(xdot, dxdot_dx, -1.0, x_at_new_timestep);
141  }

References MathLib::LinAlg::axpby(), getNewXWeight(), and getWeightedOldX().

## ◆ 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.

Warning
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::BackwardEuler.

## ◆ setInitialState()

 virtual void NumLib::TimeDiscretization::setInitialState ( const double t0 )
pure virtual

Sets the initial condition.

Implemented in NumLib::BackwardEuler.

## ◆ _dx

 std::unique_ptr NumLib::TimeDiscretization::_dx
protected

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

Definition at line 153 of file TimeDiscretization.h.

Referenced by computeRelativeChangeFromPreviousTimestep().

The documentation for this class was generated from the following files:
MathLib::VecNormType::INVALID
@ INVALID
MathLib::LinAlg::norm
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition: LinAlg.h:88
NumLib::TimeDiscretization::_dx
std::unique_ptr< GlobalVector > _dx
Used to store .
Definition: TimeDiscretization.h:153
OGS_FATAL
#define OGS_FATAL(...)
Definition: Error.h:25
NumLib::TimeDiscretization::getWeightedOldX
virtual void getWeightedOldX(GlobalVector &y, GlobalVector const &x_old) const =0
Returns .
MathLib::LinAlg::copy
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:37
MathLib::LinAlg::axpy
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
MathLib::LinAlg::axpby
void axpby(MatrixOrVector &y, double const a, double const b, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:65
MathLib::LinAlg
NumLib::TimeDiscretization::getNewXWeight
virtual double getNewXWeight() const =0
Returns .