qmat.solvers.generic ==================== .. py:module:: qmat.solvers.generic .. autoapi-nested-parse:: Submodule implementing generic solvers that can be used to solve (non-linear) ODE systems of the form : .. math:: \frac{du}{dt} = f(u,t), \quad u(0) = u_0. All those solvers are based on a :class:`DiffOp` base class, implementing : - the :math:`f(u,t)` evaluations, - a solver for :math:`u-\alpha f(u,t)=rhs`, considering given :math:`\alpha,t,rhs`. While the :math:`f(u,t)` evaluations must be implemented, a default implementation of the solver for :math:`u-\alpha f(u,t)=rhs` is provided in the base :class:`DiffOp` class. 🛠️ Various specialized :class:`DiffOp` classes are implemented in the :class:`diffops` submodule. The solvers implemented here discretizes a time-step :math:`[t_0, t_0+\Delta{t}]` into **time nodes** :math:`[t_0+\Delta{t}\tau_1, ..., t_0+\Delta{t}\tau_M]` noted :math:`[t_1,\dots,t_M]`, also called **stages** for RK methods, at which are defined the **node solutions** :math:`u_m \simeq u(t_m)`. And usually, the vector containing the node solutions :math:`{\bf u} = [u_1,\dots,u_M]^T` satisfy a **all-at-once system** : .. math:: {\bf u} - \Delta{t}Q {\bf f} = {\bf u}_0, where :math:`{\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 :math:`{\bf u}_0` is a vector containing :math:`u_0` in each entry. The :class:`CoeffSolver` allows to solve any ODE using this coefficient-based approach, either directly if the :math:`Q` matrix is lower triangular, or iteratively with SDC-based sweeps if :math:`Q` is a dense matrix. ---- An alternative solver approach relates all the node solutions using a :math:`\phi` **representation** of a time-integrator, *i.e* each node solution :math:`u_{m+1}` satisfies the following relation : .. math:: u_{m+1} -\phi(u_0, u_1, ..., u_{m}, u_{m+1}) = u_0, where :math:`\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 :class:`PhiSolver` class, that needs to be specialized by a child class implementing the :math:`\phi` function. 🛠️ Specialized :class:`PhiSolver` classes are implemented in the :class:`integrators` submodule. Submodules ---------- .. toctree:: :maxdepth: 1 /api/qmat/solvers/generic/diffops/index /api/qmat/solvers/generic/integrators/index Classes ------- .. autoapisummary:: qmat.solvers.generic.DiffOp qmat.solvers.generic.CoeffSolver qmat.solvers.generic.PhiSolver Package Contents ---------------- .. py:class:: DiffOp(u0) Base class for a differential operator :math:`f(u, t)` used in a generic ODE. It defines the evaluation of :math:`f(u, t)` at given :math:`u` and :math:`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 : .. math:: u - \alpha f(u,t) = rhs for given :math:`\alpha`, :math:`t` and :math:`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 :class:`numpy.ndarray`. :Parameters: **u0** (*array-like*) -- The initial solution associated to the differential operator, to which is extracted the generic shape and datatype of :math:`u(t)` solutions. .. py:attribute:: u0 Initial solution for the differential operator. .. py:property:: uShape Shape of a :math:`u` solution, stored as numpy array. .. py:property:: dtype Datatype of a :math:`u` solution, stored as numpy array. .. py:method:: evalF(u: numpy.ndarray, t: float, out: numpy.ndarray) :abstractmethod: Evaluate :math:`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. .. py:method:: fSolve(a: float, rhs: numpy.ndarray, t: float, out: numpy.ndarray) Solve :math:`u-\alpha f(u,t)=rhs` for given :math:`u,t,rhs`, using `out` as initial guess and storing the final result into it. :Parameters: * **a** (*float*) -- The :math:`\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. .. py:method:: test(t0=0, dt=0.1, eps=0.001, instance=None) :classmethod: 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. .. py:class:: CoeffSolver(diffOp: DiffOp, tEnd=1, nSteps=1, t0=0) Solve generic (non-linear) ODE system using :math:`Q`-coefficients with lower triangular form. It can be used to solve generic ODE systems of the form : .. math:: \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. .. py:attribute:: diffOp Differential Operator implementing :math:`f(u,t)`. .. py:attribute:: axpy BLAS-I function executing :math:`y=\alpha x + y` for any solution vectors :math:`x,y`. .. py:attribute:: t0 :value: 0 Initial simulation time. .. py:attribute:: tEnd :value: 1 Final simulation time. .. py:attribute:: nSteps :value: 1 Number of simulation time-steps .. py:attribute:: dt :value: 1.0 Time-step size for the simulation .. py:property:: u0 Initial solution for the problem .. py:property:: uShape Shape of the solution at a given time. .. py:property:: dtype Datatype of the solution at a given time. .. py:property:: times Time values for each time-step .. py:method:: evalF(u: numpy.ndarray, t: float, out: numpy.ndarray) Wrapper for the `DiffOp` function evaluating :math:`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. .. py:method:: fSolve(a: float, rhs: numpy.ndarray, t: float, out: numpy.ndarray) Wrapper for the `DiffOp` function solving :math:`u-\alpha f(u,t) = rhs`. :Parameters: * **a** (*float*) -- The :math:`\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. .. py:method:: lowerTri(Q: numpy.ndarray, strict=False) :staticmethod: 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. :rtype: bool .. py:method:: solve(Q, weights, uNum=None, tInit=0) Solve the ODE considering **lower-triangular** :math:`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) :math:`u_{m}` that is solved using previously computed node solution : .. math:: 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 :math:`t_m = t_0 + \tau_m` and :math:`q_{i,j}` are the coefficients :math:`Q`. Finally, the **step update** is done using all computed node solutions : .. math:: u(t_0+\Delta{t}) \simeq u_0 + \sum_{m=1}^{M} \omega_{m} f(u_m, t_m), where :math:`\omega_{m}` are the weights associated to the :math:`Q`-coefficients. If no weights are provided, then it simply uses the last node solution for the step update : .. math:: u(t_0+\Delta{t}) \simeq u_M :Parameters: * **Q** (*np.2darray-like*) -- The **lower-triangular** :math:`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 :math:`\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. :rtype: np.ndarray .. py:method:: solveSDC(nSweeps, Q, weights, QDelta, uNum=None, tInit=0) Solve the ODE with dense :math:`Q` coefficients using SDC sweeps. Considering a **lower-triangular** approximation :math:`Q_\Delta` of :math:`Q`, it performes for each time-step :math:`K` SDC sweeps : .. math:: \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} where :math:`q^\Delta_{i,j}` and :math:`q_{i,j}` are the coefficients of :math:`Q_\Delta` and :math:`Q`, respectively. It uses a **copy initialization**, that is :math:`u_{m}^0 = u_0`. Finally, the **step update** is done using all computed node solutions : .. math:: u(t_0+\Delta{t}) \simeq u_0 + \sum_{m=1}^{M} \omega_{m} f(u_m, t_m), where :math:`\omega_{m}` are the weights associated to the :math:`Q`-coefficients. If no weights are provided, then it simply uses the last node solution for the step update : .. math:: u(t_0+\Delta{t}) \simeq u_M :Parameters: * **nSweeps** (*int*) -- Number of SDC sweeps :math:`K`. * **Q** (*2D array-like*) -- The dense :math:`Q` matrix. * **weights** (*1D array-like*) -- The associated weights :math:`\omega_{m}` for the step update. * **QDelta** (*2D array-like*) -- The lower-triangular :math:`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. :rtype: np.ndarray .. py:class:: PhiSolver(diffOp: DiffOp, nodes, tEnd=1, nSteps=1, t0=0) .. autoapi-inheritance-diagram:: qmat.solvers.generic.PhiSolver :parts: 1 Solve generic (non-linear) ODE system using :math:`\phi` representation of time-integration solvers. It consider the following ODE : .. math:: \frac{du}{dt} = f(u,t), and compute for each step the solution on **time nodes** :math:`\tau_1, ..., \tau_M` by soving the following system : .. math:: u_{m+1} -\phi(u_0, u_1, ..., u_{m}, u_{m+1}) = u_0. It uses then per default the last node solution :math:`u_{M}` as initial solution for the next step. ⚙️ Requires the implementation of an `evalPhi` method that evaluates the :math:`\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 :math:`\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. .. py:attribute:: nodes Time nodes for each time-step of the time-integrator. .. py:property:: nNodes Number of time-nodes .. py:method:: evalPhi(uVals, fEvals, out, t0=0) :abstractmethod: Evaluate the :math:`\phi` operator on time-node up to :math:`u_{m+1}`. Considering :math:`u_0, u_1, \dots, u_{m+1}`, if evaluates : .. math:: \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 :math:`f(u_0,t_0),f(u_1,\tau_1),...,f(u_{m},\tau_{m})` as arguments, in order to avoid any additional :math:`f(u,t)` evaluations. :Parameters: * **uVals** (list[np.ndarray] of size :math:`m+2`) -- The :math:`m+1` time-node solutions + the initial solution :math:`u_0`. * **fEvals** (list[np.ndarray] of size :math:`m+1` or :math:`m+1`) -- The :math:`f(u,t)` evaluations at each time nodes (+ initial solution), up to time-node :math:`m`. It can eventually contain a pre-computed :math:`f_{m+1}` to spare one :math:`f(u,t)` evaluation. * **out** (*np.ndarray*) -- Array used to store the evaluation. * **t0** (*float, optional*) -- Initial step time. The default is 0. .. py:method:: phiSolve(uPrev, fEvals, out, rhs=0, t0=0) Solve the node update at given time-node :math:`\tau_{m+1}`. Considering :math:`m+1` previous known node solutions :math:`u_0, u_1, ..., u_{m}`, it solves the following system : .. math:: 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 :math:`f` evaluations :math:`f_0, f_1, ..., f_{m}` to avoid supplementar re-computing those. :Parameters: * **uPrev** (list[np.ndarray] of size :math:`m+1`) -- The previous node solutions :math:`u_0, u_1, ..., u_{m}`. * **fEvals** (list[np.ndarray] of size :math:`m+1`) -- Evaluations of previous node solutions :math:`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. .. py:method:: stepUpdate(u0, uNodes, fEvals, out) Update end-step solution to be used as initial guess for next step. .. note:: This method has to ensures that fEvals[0] contains the :math:`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 :math:`u_1,\dots,u_M`. * **fEvals** (*list[np.ndarray]*) -- Precomputed node evaluation :math:`f_1,\dots,f_M`. * **out** (*np.ndarray*) -- Output array to store the result. .. py:method:: solve(uNum=None, tInit=0) Solve using sequential computation of node solutions for each step, using the relation : .. math:: 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 :math:`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. :rtype: np.ndarray .. py:method:: solveSDC(nSweeps, Q=None, weights=None, uNum=None, tInit=0) Solve the ODE with dense :math:`Q` coefficients using SDC sweeps. Considering a **lower-triangular** approximation :math:`Q_\Delta` of :math:`Q`, it performes for each time-step :math:`K` SDC sweeps : .. math:: 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 :math:`\phi_m^k:=\phi(u_0,u_1^k,\dots,u_m^k,f_0,f_1^k,\dots,f_{m-1}^k)` and :math:`q_{i,j}` are the coefficients of the :math:`Q` matrix. It uses a **copy initialization**, that is :math:`u_{m}^0 = u_0`. 💡 If we consider that :math:`\phi_m^{k}` is like a coarse solver applied on iteration :math:`k` and :math:`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 :math:`k`, then the SDC correction above furiously resemble to a **Parareal iteration** 👻 👻 👻 Finally, the **step update** is done using all computed node solutions : .. math:: u(t_0+\Delta{t}) \simeq u_0 + \sum_{m=1}^{M} \omega_{m} f(u_m, t_m), where :math:`\omega_{m}` are the weights associated to the :math:`Q`-coefficients. If weights are not used (`weights=False`), then it simply uses the last node solution for the step update : .. math:: u(t_0+\Delta{t}) \simeq u_M :Parameters: * **nSweeps** (*int*) -- Number of SDC sweeps :math:`K`. * **Q** (*2D array-like, optional*) -- The dense :math:`Q` matrix. If not provided, automatically computed using the :class:`LagrangeApproximation` class and the solver nodes. * **weights** (*1D array-like, optional*) -- The associated weights :math:`\omega_{m}` for the step update. If not provided, automatically computed using the :class:`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. :rtype: np.ndarray