Skip to content

Commit

Permalink
Merge pull request #4115 from mfem/operator-doc
Browse files Browse the repository at this point in the history 
Refactored TimeDependentOperator Documentation [operator-doc]
  • Loading branch information
@tzanio
tzanio committed May 20, 2024
2 parents c444b17 + 3e83791 commit 800b178
Showing 1 changed file with 135 additions and 57 deletions.
192 changes: 135 additions & 57 deletions linalg/operator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -309,48 +309,82 @@ class Operator


/// Base abstract class for first order time dependent operators.
/** Operator of the form: (x,t) -> f(x,t), where k = f(x,t) generally solves the
algebraic equation F(x,k,t) = G(x,t). The functions F and G represent the
_implicit_ and _explicit_ parts of the operator, respectively. For explicit
operators, F(x,k,t) = k, so f(x,t) = G(x,t). */
/** Operator of the form: (u,t) -> k(u,t), where k generally solves the
algebraic equation F(u,k,t) = G(u,t). The functions F and G represent the
_implicit_ and _explicit_ parts of the operator, respectively.
A common use for this class is representing a differential algebraic
equation of the form $ F(y,\frac{dy}{dt},t) = G(y,t) $.
For example, consider an ordinary differential equation of the form
$ M \frac{dy}{dt} = g(y,t) $. There are various ways of expressing this ODE
as a TimeDependentOperator depending on the choices for F and G. Here are
some common choices:
1. F(u,k,t) = k and G(u,t) = inv(M) g(u,t),
2. F(u,k,t) = M k and G(u,t) = g(u,t),
3. F(u,k,t) = M k - g(u,t) and G(u,t) = 0.
Note that depending on the ODE solver, some of the above choices may be
preferable to the others.
*/
class TimeDependentOperator : public Operator
{
public:
/// Enum used to describe the form of the time-dependent operator.
/** The type should be set by classes derived from TimeDependentOperator to
describe the form, in terms of the functions F and G, used by the
specific derived class. This information can be queried by classes or
functions (like time stepping algorithms) to make choices about the
algorithm to use, or to ensure that the TimeDependentOperator uses the
form expected by the class/function.
For example, assume that a derived class is implementing the ODE
$M \frac{dy}{dt} = g(y,t)$ and chooses to define $F(u,k,t) = M k$ and
$G(u,t) = g(u,t)$. Then it cannot use type EXPLICIT, unless $M = I$, or
type HOMOGENEOUS, unless $g(u,t) = 0$. If, on the other hand, the derived
class chooses to define $F(u,k,t) = k$ and $G(u,t) = M^{-1} g(y,t)$, then
the natural choice is to set the type to EXPLICIT, even though setting it
to IMPLICIT is also not wrong -- doing so will simply fail to inform
methods that query this information that it uses a more specific
implementation, EXPLICIT, that may allow the use of algorithms that
support only the EXPLICIT type. */
enum Type
{
EXPLICIT, ///< This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
EXPLICIT, ///< This type assumes F(u,k,t) = k.
IMPLICIT, ///< This is the most general type, no assumptions on F and G.
HOMOGENEOUS ///< This type assumes that G(x,t) = 0.
HOMOGENEOUS ///< This type assumes that G(u,t) = 0.
};

/// Evaluation mode. See SetEvalMode() for details.
enum EvalMode
{
/** Normal evaluation. */
NORMAL,
/** Assuming additive split, f(x,t) = f1(x,t) + f2(x,t), evaluate the
first term, f1. */
/** Assuming additive split, k(u,t) = k1(u,t) + k2(u,t), evaluate the
first term, k1. */
ADDITIVE_TERM_1,
/** Assuming additive split, f(x,t) = f1(x,t) + f2(x,t), evaluate the
second term, f2. */
/** Assuming additive split, k(u,t) = k1(u,t) + k2(u,t), evaluate the
second term, k2. */
ADDITIVE_TERM_2
};

protected:
real_t t; ///< Current time.
Type type; ///< Describes the form of the TimeDependentOperator.
Type type; /**< @brief Describes the form of the TimeDependentOperator, see
the documentation of #Type. */
EvalMode eval_mode; ///< Current evaluation mode.

public:
/** @brief Construct a "square" TimeDependentOperator y = f(x,t), where x and
y have the same dimension @a n. */
/** @brief Construct a "square" TimeDependentOperator (u,t) -> k(u,t), where
u and k have the same dimension @a n. */
explicit TimeDependentOperator(int n = 0, real_t t_ = 0.0,
Type type_ = EXPLICIT)
: Operator(n) { t = t_; type = type_; eval_mode = NORMAL; }

/** @brief Construct a TimeDependentOperator y = f(x,t), where x and y have
dimensions @a w and @a h, respectively. */
TimeDependentOperator(int h, int w, real_t t_ = 0.0, Type type_ = EXPLICIT)
/** @brief Construct a TimeDependentOperator (u,t) -> k(u,t), where u and k
have dimensions @a w and @a h, respectively. */
TimeDependentOperator(int h, int w, double t_ = 0.0, Type type_ = EXPLICIT)
: Operator(h, w) { t = t_; type = type_; eval_mode = NORMAL; }

/// Read the currently set time.
Expand All @@ -373,7 +407,7 @@ class TimeDependentOperator : public Operator
/** The evaluation mode is a switch that allows time-stepping methods to
request evaluation of separate components/terms of the time-dependent
operator. For example, IMEX methods typically assume additive split of
the operator: f(x,t) = f1(x,t) + f2(x,t) and they rely on the ability to
the operator: k(u,t) = k1(u,t) + k2(u,t) and they rely on the ability to
evaluate the two terms separately.
Generally, setting the evaluation mode should affect the behavior of all
Expand All @@ -384,62 +418,104 @@ class TimeDependentOperator : public Operator
{ eval_mode = new_eval_mode; }

/** @brief Perform the action of the explicit part of the operator, G:
@a y = G(@a x, t) where t is the current time.
@a v = G(@a u, t) where t is the current time.
Presently, this method is used by some PETSc ODE solvers, for more
details, see the PETSc Manual. */
virtual void ExplicitMult(const Vector &x, Vector &y) const;
virtual void ExplicitMult(const Vector &u, Vector &v) const;

/** @brief Perform the action of the implicit part of the operator, F:
@a y = F(@a x, @a k, t) where t is the current time.
@a v = F(@a u, @a k, t) where t is the current time.
Presently, this method is used by some PETSc ODE solvers, for more
details, see the PETSc Manual.*/
virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const;

/** @brief Perform the action of the operator: @a y = k = f(@a x, t), where
k solves the algebraic equation F(@a x, k, t) = G(@a x, t) and t is the
current time. */
virtual void Mult(const Vector &x, Vector &y) const;

/** @brief Solve the equation: @a k = f(@a x + @a dt @a k, t), for the
unknown @a k at the current time t.
For general F and G, the equation for @a k becomes:
F(@a x + @a dt @a k, @a k, t) = G(@a x + @a dt @a k, t).
The input vector @a x corresponds to time index (or cycle) n, while the
currently set time, #t, and the result vector @a k correspond to time
index n+1. The time step @a dt corresponds to the time interval between
cycles n and n+1.
This method allows for the abstract implementation of some time
integration methods, including diagonal implicit Runge-Kutta (DIRK)
methods and the backward Euler method in particular.
virtual void ImplicitMult(const Vector &u, const Vector &k, Vector &v) const;

/** @brief Perform the action of the operator (u,t) -> k(u,t) where t is the
current time set by SetTime() and @a k satisfies
F(@a u, @a k, t) = G(@a u, t).
For solving an ordinary differential equation of the form
$ M \frac{dy}{dt} = g(y,t) $, recall that F and G can be defined in
various ways, e.g.:
1. F(u,k,t) = k and G(u,t) = inv(M) g(u,t)
2. F(u,k,t) = M k and G(u,t) = g(u,t)
3. F(u,k,t) = M k - g(u,t) and G(u,t) = 0.
Regardless of the choice of F and G, this function should always compute
@a k = inv(M) g(@a u, t). */
virtual void Mult(const Vector &u, Vector &v) const override;

/** @brief Solve for the unknown @a k, at the current time t, the following
equation:
F(@a u + @a gamma @a k, @a k, t) = G(@a u + @a gamma @a k, t).
For solving an ordinary differential equation of the form
$ M \frac{dy}{dt} = g(y,t) $, recall that F and G can be defined in
various ways, e.g.:
1. F(u,k,t) = k and G(u,t) = inv(M) g(u,t)
2. F(u,k,t) = M k and G(u,t) = g(u,t)
3. F(u,k,t) = M k - g(u,t) and G(u,t) = 0
Regardless of the choice of F and G, this function should solve for @a k
in M @a k = g(@a u + @a gamma @a k, t).
To see how @a k can be useful, consider the backward Euler method defined
by $ y(t + \Delta t) = y(t) + \Delta t k_0 $ where
$ M k_0 = g \big( y(t) + \Delta t k_0, t + \Delta t \big) $. A backward
Euler integrator can use @a k from this function for $k_0$, with the call
using @a u set to $ y(t) $, @a gamma set to $ \Delta t$, and time set to
$t + \Delta t$. See class BackwardEulerSolver.
Generalizing further, consider a diagonally implicit Runge-Kutta (DIRK)
method defined by
$ y(t + \Delta t) = y(t) + \Delta t \sum_{i=1}^s b_i k_i $ where
$ M k_i = g \big( y(t) + \Delta t \sum_{j=1}^i a_{ij} k_j,
t + c_i \Delta t \big) $.
A DIRK integrator can use @a k from this function, with @a u set to
$ y(t) + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j $ and @a gamma set to
$ a_{ii} \Delta t $, for $ k_i $. For example, see class SDIRK33Solver.
If not re-implemented, this method simply generates an error. */
virtual void ImplicitSolve(const real_t dt, const Vector &x, Vector &k);
virtual void ImplicitSolve(const real_t gamma, const Vector &u, Vector &k);

/** @brief Return an Operator representing (dF/dk @a shift + dF/dx) at the
given @a x, @a k, and the currently set time.
/** @brief Return an Operator representing (dF/dk @a shift + dF/du) at the
given @a u, @a k, and the currently set time.
Presently, this method is used by some PETSc ODE solvers, for more
details, see the PETSc Manual. */
virtual Operator& GetImplicitGradient(const Vector &x, const Vector &k,
virtual Operator& GetImplicitGradient(const Vector &u, const Vector &k,
real_t shift) const;

/** @brief Return an Operator representing dG/dx at the given point @a x and
/** @brief Return an Operator representing dG/du at the given point @a u and
the currently set time.
Presently, this method is used by some PETSc ODE solvers, for more
details, see the PETSc Manual. */
virtual Operator& GetExplicitGradient(const Vector &x) const;
virtual Operator& GetExplicitGradient(const Vector &u) const;

/** @brief Setup the ODE linear system $ A(x,t) = (I - gamma J) $ or
$ A = (M - gamma J) $, where $ J(x,t) = \frac{df}{dt(x,t)} $.
/** @brief Setup a linear system as needed by some SUNDIALS ODE solvers.
For solving an ordinary differential equation of the form
$ M \frac{dy}{dt} = g(y,t) $, recall that F and G can be defined as one
of the following:
1. F(u,k,t) = k and G(u,t) = inv(M) g(u,t)
2. F(u,k,t) = M k and G(u,t) = g(u,t)
3. F(u,k,t) = M k - g(u,t) and G(u,t) = 0
This function performs setup to solve $ A x = b $ where A is either
1. A(@a y,t) = I - @a gamma inv(M) J(@a y,t)
2. A(@a y,t) = M - @a gamma J(@a y,t)
3. A(@a y,t) = M - @a gamma J(@a y,t)
with J = dg/dy (or a reasonable approximation thereof).
@param[in] x The state at which $A(x,t)$ should be evaluated.
@param[in] fx The current value of the ODE rhs function, $f(x,t)$.
@param[in] y The state at which A(@a y,t) should be evaluated.
@param[in] v The value of inv(M) g(y,t) for 1 or g(y,t) for 2 & 3.
@param[in] jok Flag indicating if the Jacobian should be updated.
@param[out] jcur Flag to signal if the Jacobian was updated.
@param[in] gamma The scaled time step value.
Expand All @@ -448,10 +524,10 @@ class TimeDependentOperator : public Operator
Presently, this method is used by SUNDIALS ODE solvers, for more
details, see the SUNDIALS User Guides. */
virtual int SUNImplicitSetup(const Vector &x, const Vector &fx,
virtual int SUNImplicitSetup(const Vector &y, const Vector &v,
int jok, int *jcur, real_t gamma);

/** @brief Solve the ODE linear system $ A x = b $ as setup by
/** @brief Solve the ODE linear system A @a x = @a b, where A is defined by
the method SUNImplicitSetup().
@param[in] b The linear system right-hand side.
Expand All @@ -464,16 +540,17 @@ class TimeDependentOperator : public Operator
details, see the SUNDIALS User Guides. */
virtual int SUNImplicitSolve(const Vector &b, Vector &x, real_t tol);

/** @brief Setup the mass matrix in the ODE system $ M y' = f(y,t) $ .
/** @brief Setup the mass matrix in the ODE system
$ M \frac{dy}{dt} = g(y,t) $ .
If not re-implemented, this method simply generates an error.
Presently, this method is used by SUNDIALS ARKStep integrator, for more
details, see the ARKode User Guide. */
virtual int SUNMassSetup();

/** @brief Solve the mass matrix linear system $ M x = b $
as setup by the method SUNMassSetup().
/** @brief Solve the mass matrix linear system M @a x = @a b, where M is
defined by the method SUNMassSetup().
@param[in] b The linear system right-hand side.
@param[in,out] x On input, the initial guess. On output, the solution.
Expand All @@ -485,7 +562,8 @@ class TimeDependentOperator : public Operator
details, see the ARKode User Guide. */
virtual int SUNMassSolve(const Vector &b, Vector &x, real_t tol);

/** @brief Compute the mass matrix-vector product $ v = M x $ .
/** @brief Compute the mass matrix-vector product @a v = M @a x, where M is
defined by the method SUNMassSetup().
@param[in] x The vector to multiply.
@param[out] v The result of the matrix-vector product.
Expand Down

0 comments on commit 800b178

Please sign in to comment.