qmat.solvers.generic

Submodule implementing generic solvers that can be used to solve (non-linear) ODE systems of the form :

\[\frac{du}{dt} = f(u,t), \quad u(0) = u_0.\]

All those solvers are based on a DiffOp base class, implementing :

the \(f(u,t)\) evaluations,
a solver for \(u-\alpha f(u,t)=rhs\), considering given \(\alpha,t,rhs\).

While the \(f(u,t)\) evaluations must be implemented, a default implementation of the solver for \(u-\alpha f(u,t)=rhs\) is provided in the base DiffOp class.

🛠️ Various specialized DiffOp classes are implemented in the diffops submodule.

The solvers implemented here discretizes a time-step \([t_0, t_0+\Delta{t}]\) into time nodes \([t_0+\Delta{t}\tau_1, ..., t_0+\Delta{t}\tau_M]\) noted \([t_1,\dots,t_M]\), also called stages for RK methods, at which are defined the node solutions \(u_m \simeq u(t_m)\). And usually, the vector containing the node solutions \({\bf u} = [u_1,\dots,u_M]^T\) satisfy a all-at-once system :

\[{\bf u} - \Delta{t}Q {\bf f} = {\bf u}_0,\]

where \({\bf f} = [f(u_1, t_1),\dots,f(u_M,t_M)]^T\) is the vector with the evaluations of each node solutions and \({\bf u}_0\) is a vector containing \(u_0\) in each entry. The CoeffSolver allows to solve any ODE using this coefficient-based approach, either directly if the \(Q\) matrix is lower triangular, or iteratively with SDC-based sweeps if \(Q\) is a dense matrix.

An alternative solver approach relates all the node solutions using a \(\phi\) representation of a time-integrator, i.e each node solution \(u_{m+1}\) satisfies the following relation :

\[u_{m+1} -\phi(u_0, u_1, ..., u_{m}, u_{m+1}) = u_0,\]

where \(\phi\) is solely defined by the chosen time-integrator. The system above can be solved node-by-node in a sequential approach, or iteratively with a SDC-based approach. It is implemented in the abstract PhiSolver class, that needs to be specialized by a child class implementing the \(\phi\) function.

🛠️ Specialized PhiSolver classes are implemented in the integrators submodule.

Submodules

Classes

`DiffOp`	Base class for a differential operator \(f(u, t)\) used in a generic ODE.
`CoeffSolver`	Solve generic (non-linear) ODE system using \(Q\)-coefficients with lower triangular form.
`PhiSolver`	Solve generic (non-linear) ODE system using \(\phi\) representation of time-integration solvers.

Package Contents

class DiffOp(u0)[source]

Base class for a differential operator \(f(u, t)\) used in a generic ODE.

It defines the evaluation of \(f(u, t)\) at given \(u\) and \(t\) with a evalF(u, t, out) method, that put the result of the evaluation in the out array.

Additionally, this class defines a default fSolve method that solves :

\[u - \alpha f(u,t) = rhs\]

for given \(\alpha\), \(t\) and \(rhs\). This default method can be overridden by a more efficient specific method for a specific differential operator.

Note

Solutions are stored in N-dimensional numpy.ndarray.

Parameters:: u0 (array-like) – The initial solution associated to the differential operator, to which is extracted the generic shape and datatype of \(u(t)\) solutions.

u0: Initial solution for the differential operator.

property uShape: Shape of a \(u\) solution, stored as numpy array.

property dtype: Datatype of a \(u\) solution, stored as numpy array.

abstractmethod evalF(u: numpy.ndarray, t: float, out: numpy.ndarray)[source]

Evaluate \(f(u,t)\) and store the result into out.

Parameters:

u (np.ndarray) – Input solution for the evaluation.
t (float) – Time for the evaluation.
out (np.ndarray) – Output array in which is stored the evaluation.

fSolve(a: float, rhs: numpy.ndarray, t: float, out: numpy.ndarray)[source]

Solve \(u-\alpha f(u,t)=rhs\) for given \(u,t,rhs\), using out as initial guess and storing the final result into it.

Parameters:

a (float) – The \(\alpha\) coefficient.
rhs (np.ndarray) – The right hand side.
t (float) – Time for the evaluation.
out (np.ndarray) – Input-output array used as initial guess, in which is stored the solution.

classmethod test(t0=0, dt=0.1, eps=0.001, instance=None)[source]

Class method to test the DiffOp implementation.

Parameters:

t0 (float, optional) – Evaluation time to test the instance. The default is 0.
dt (float, optional) – Time-step to test the fSolve method. The default is 1e-1.
eps (float, optional) – Perturbation added in the expected solution to test the fSolve method. The default is 1e-3.
instance (DiffOp, optional) – Instance to be tested. If not provided (None), an instance is created using the default constructor.

class CoeffSolver(diffOp: DiffOp, tEnd=1, nSteps=1, t0=0)[source]

Solve generic (non-linear) ODE system using \(Q\)-coefficients with lower triangular form.

It can be used to solve generic ODE systems of the form :

\[\frac{du}{dt} = f(u,t), \quad u(0)=u_0.\]

Parameters:

diffOp (DiffOp) – Differential operator for the ODE.
tEnd (float, optional) – Final simulation time. The default is 1.
nSteps (int, optional) – Number of simulation time-steps. The default is 1.
t0 (float, optional) – Initial simulation time. The default is 0.

diffOp: Differential Operator implementing \(f(u,t)\).

axpy: BLAS-I function executing \(y=\alpha x + y\) for any solution vectors \(x,y\).

t0 = 0: Initial simulation time.

tEnd = 1: Final simulation time.

nSteps = 1: Number of simulation time-steps

dt = 1.0: Time-step size for the simulation

property u0: Initial solution for the problem

property uShape: Shape of the solution at a given time.

property dtype: Datatype of the solution at a given time.

property times: Time values for each time-step

evalF(u: numpy.ndarray, t: float, out: numpy.ndarray)[source]

Wrapper for the DiffOp function evaluating \(f(u,t)\).

Parameters:

u (np.ndarray) – Input solution for the evaluation.
t (float) – Time for the evaluation.
out (np.ndarray) – Output array in which is stored the evaluation.

fSolve(a: float, rhs: numpy.ndarray, t: float, out: numpy.ndarray)[source]

Wrapper for the DiffOp function solving \(u-\alpha f(u,t) = rhs\).

Parameters:

a (float) – The \(\alpha\) coefficient.
rhs (np.ndarray) – The right hand side.
t (float) – Time for the evaluation.
out (np.ndarray) – Input-output array used as initial guess, in which is stored the solution.

static lowerTri(Q: numpy.ndarray, strict=False)[source]

Check if a 2D matrix is lower triangular.

Parameters:

Q (np.ndarray) – Matrix to check.
strict (bool, optional) – Check for strictly lower triangular matrix. The default is False.

Returns:

Is the matrix (strictly) lower triangular or not.

Return type:

bool

solve(Q, weights, uNum=None, tInit=0)[source]

Solve the ODE considering lower-triangular \(Q\) coefficients.

This is equivalent to the classical implementation of a generic Runge-Kutta method using its Butcher table. For each time-step, it defines a node solution (or stage) \(u_{m}\) that is solved using previously computed node solution :

\[u_{m} - \Delta{t}q_{m,m}f(u_m,t_m) = u_0 + \Delta{t}\sum_{j=1}^{m-1}q_{m,j}f(u_j, t_j),\]

where \(t_m = t_0 + \tau_m\) and \(q_{i,j}\) are the coefficients \(Q\). Finally, the step update is done using all computed node solutions :

\[u(t_0+\Delta{t}) \simeq u_0 + \sum_{m=1}^{M} \omega_{m} f(u_m, t_m),\]

where \(\omega_{m}\) are the weights associated to the \(Q\)-coefficients. If no weights are provided, then it simply uses the last node solution for the step update :

\[u(t_0+\Delta{t}) \simeq u_M\]

Parameters:

Q (np.2darray-like) – The lower-triangular \(Q\)-coefficients matrix.
weights (np.1darray-like) – The associated :math:omega_{m}` weights. If not provided, use the last node solution for the update (requires \(\tau_{M} = 1\)).
uNum (np.ndarray, optional) – Array of shape (nSteps+1,*uShape), that can be use to store the result and avoid creating it internally. The default is None.
tInit (float, optional) – Initial time offset to be added to solver’s own t0 for successive solve calls. The default is 0.

Returns:

uNum – Array of shape (nSteps+1,*uShape) that stores the solution at each time-step.

Return type:

np.ndarray

solveSDC(nSweeps, Q, weights, QDelta, uNum=None, tInit=0)[source]

Solve the ODE with dense \(Q\) coefficients using SDC sweeps.

Considering a lower-triangular approximation \(Q_\Delta\) of \(Q\), it performes for each time-step \(K\) SDC sweeps :

\[\begin{split}\begin{align} u_{m}^{k+1} - \Delta{t}q^\Delta_{m,m}f(u_m^{k+1},t_m) =&~ u_0 + \Delta{t}\sum_{j=1}^{M}q_{m,j}f(u_j^k, t_j) \\ &+ \Delta{t}\sum_{j=1}^{m-1}q^\Delta_{m,j}f(u_j^{k+1},t_j) - \Delta{t}\sum_{j=1}^{m}q^\Delta_{m,j}f(u_j^{k},t_j), \end{align}\end{split}\]

where \(q^\Delta_{i,j}\) and \(q_{i,j}\) are the coefficients of \(Q_\Delta\) and \(Q\), respectively. It uses a copy initialization, that is \(u_{m}^0 = u_0\).

Finally, the step update is done using all computed node solutions :

\[u(t_0+\Delta{t}) \simeq u_0 + \sum_{m=1}^{M} \omega_{m} f(u_m, t_m),\]

where \(\omega_{m}\) are the weights associated to the \(Q\)-coefficients. If no weights are provided, then it simply uses the last node solution for the step update :

\[u(t_0+\Delta{t}) \simeq u_M\]

Parameters:

nSweeps (int) – Number of SDC sweeps \(K\).
Q (2D array-like) – The dense \(Q\) matrix.
weights (1D array-like) – The associated weights \(\omega_{m}\) for the step update.
QDelta (2D array-like) – The lower-triangular \(Q_\Delta\) matrix.
uNum (np.ndarray, optional) – Array of shape (nSteps+1,*uShape), that can be use to store the result and avoid creating it internally. The default is None.
tInit (float, optional) – Initial time offset to be added to solver’s own t0 for successive solve calls. The default is 0.

Returns:

uNum – Array of shape (nSteps+1,*uShape) that stores the solution at each time-step.

Return type:

np.ndarray

class PhiSolver(diffOp: DiffOp, nodes, tEnd=1, nSteps=1, t0=0)[source]

Inheritance diagram of qmat.solvers.generic.PhiSolver

Solve generic (non-linear) ODE system using \(\phi\) representation of time-integration solvers.

It consider the following ODE :

\[\frac{du}{dt} = f(u,t),\]

and compute for each step the solution on time nodes \(\tau_1, ..., \tau_M\) by soving the following system :

\[u_{m+1} -\phi(u_0, u_1, ..., u_{m}, u_{m+1}) = u_0.\]

It uses then per default the last node solution \(u_{M}\) as initial solution for the next step.

⚙️ Requires the implementation of an evalPhi method that evaluates the \(\phi\) function. Also, a default phiSolve method is implemented, that solves the system above, and can be overridden for specific time-integrator (in particular for explicit time-integrators). Finally, it implements a default stepUpdate method that setup the next time-step using the last time-node solution.

Parameters:

diffOp (DiffOp) – Differential operator for the ODE.
nodes (1D array-like) – The time nodes \(\tau_1, ..., \tau_M\).
tEnd (float, optional) – Final simulation time. The default is 1.
nSteps (int, optional) – Number of simulation time-steps. The default is 1.
t0 (float, optional) – Initial simulation time. The default is 0.

nodes: Time nodes for each time-step of the time-integrator.

property nNodes: Number of time-nodes

abstractmethod evalPhi(uVals, fEvals, out, t0=0)[source]

Evaluate the \(\phi\) operator on time-node up to \(u_{m+1}\).

Considering \(u_0, u_1, \dots, u_{m+1}\), if evaluates :

\[\phi(u_0, u_1, ..., u_{m}, u_{m+1}),\]

and store its value into the output vector out. It also takes the node evaluation \(f(u_0,t_0),f(u_1,\tau_1),...,f(u_{m},\tau_{m})\) as arguments, in order to avoid any additional \(f(u,t)\) evaluations.

Parameters:

uVals (list[np.ndarray] of size \(m+2\)) – The \(m+1\) time-node solutions + the initial solution \(u_0\).
fEvals (list[np.ndarray] of size \(m+1\) or \(m+1\)) – The \(f(u,t)\) evaluations at each time nodes (+ initial solution), up to time-node \(m\). It can eventually contain a pre-computed \(f_{m+1}\) to spare one \(f(u,t)\) evaluation.
out (np.ndarray) – Array used to store the evaluation.
t0 (float, optional) – Initial step time. The default is 0.

phiSolve(uPrev, fEvals, out, rhs=0, t0=0)[source]

Solve the node update at given time-node \(\tau_{m+1}\).

Considering \(m+1\) previous known node solutions \(u_0, u_1, ..., u_{m}\), it solves the following system :

\[u -\phi(u_0, u_1, ..., u_{m}, u) = rhs,\]

where the value given in out is used as initial guess and to store the computed solution. It also takes as argument the \(f\) evaluations \(f_0, f_1, ..., f_{m}\) to avoid supplementar re-computing those.

Parameters:

uPrev (list[np.ndarray] of size \(m+1\)) – The previous node solutions \(u_0, u_1, ..., u_{m}\).
fEvals (list[np.ndarray] of size \(m+1\)) – Evaluations of previous node solutions \(f_0, f_1, ..., f_{m}\).
out (np.ndarray) – Array with the initial guess, used to store the final solution.
rhs (np.ndarray or float, optional) – Right hand side used to solve the equation above. The default is 0.
t0 (float, optional) – Initial step size. The default is 0.

stepUpdate(u0, uNodes, fEvals, out)[source]

Update end-step solution to be used as initial guess for next step.

Note

This method has to ensures that fEvals[0] contains the \(f(u,t)\) evaluation of the next step initial solution.

Parameters:

u0 (np.ndarray) – Initial solution for the current step.
uNodes (list[np.ndarray]) – Precomputed node solutions \(u_1,\dots,u_M\).
fEvals (list[np.ndarray]) – Precomputed node evaluation \(f_1,\dots,f_M\).
out (np.ndarray) – Output array to store the result.

solve(uNum=None, tInit=0)[source]

Solve using sequential computation of node solutions for each step, using the relation :

\[u_{m+1} -\phi(u_0, u_1, ..., u_{m}, u_{m+1}, f_0, f_1, ..., f_{m}) = u_0.\]

and the step update to compute \(u(t_0+\Delta_t)\) using all computed node solutions.

Parameters:

uNum (np.ndarray, optional) – Array of shape (nSteps+1,*uShape), that can be use to store the result and avoid creating it internally. The default is None.
tInit (float, optional) – Initial time offset to be added to solver’s own t0 for successive solve calls. The default is 0.

Returns:

uNum – Array of shape (nSteps+1,*uShape) that stores the solution at each time-step.

Return type:

np.ndarray

solveSDC(nSweeps, Q=None, weights=None, uNum=None, tInit=0)[source]

Solve the ODE with dense \(Q\) coefficients using SDC sweeps.

Considering a lower-triangular approximation \(Q_\Delta\) of \(Q\), it performes for each time-step \(K\) SDC sweeps :

\[u_{m}^{k+1} - \phi_m^{k+1} = u_0 + \Delta{t}\sum_{j=1}^{M}q_{m,j}f(u_j^k, t_j) - \phi_m^k,\]

where \(\phi_m^k:=\phi(u_0,u_1^k,\dots,u_m^k,f_0,f_1^k,\dots,f_{m-1}^k)\) and \(q_{i,j}\) are the coefficients of the \(Q\) matrix. It uses a copy initialization, that is \(u_{m}^0 = u_0\).

💡 If we consider that \(\phi_m^{k}\) is like a coarse solver applied on iteration \(k\) and \(u_0 + \Delta{t}\sum_{j=1}^{M}q_{m,j}f(u_j^k, t_j)\) is like a fine solver applied to iteration \(k\), then the SDC correction above furiously resemble to a Parareal iteration 👻 👻 👻

Finally, the step update is done using all computed node solutions :

\[u(t_0+\Delta{t}) \simeq u_0 + \sum_{m=1}^{M} \omega_{m} f(u_m, t_m),\]

where \(\omega_{m}\) are the weights associated to the \(Q\)-coefficients. If weights are not used (weights=False), then it simply uses the last node solution for the step update :

\[u(t_0+\Delta{t}) \simeq u_M\]

Parameters:

nSweeps (int) – Number of SDC sweeps \(K\).
Q (2D array-like, optional) – The dense \(Q\) matrix. If not provided, automatically computed using the LagrangeApproximation class and the solver nodes.
weights (1D array-like, optional) – The associated weights \(\omega_{m}\) for the step update. If not provided, automatically computed using the LagrangeApproximation class and the solver nodes.
uNum (np.ndarray, optional) – Array of shape (nSteps+1,*uShape), that can be use to store the result and avoid creating it internally. The default is None.
tInit (float, optional) – Initial time offset to be added to solver’s own t0 for successive solve calls. The default is 0.

Returns:

uNum – Array of shape (nSteps+1,*uShape) that stores the solution at each time-step.

Return type:

np.ndarray