|
| 1 | +# The AbstractSciMLOperator Interface |
| 2 | + |
| 3 | +## Formal Properties of DiffEqOperators |
| 4 | + |
| 5 | +These are the formal properties that an `AbstractSciMLOperator` should obey |
| 6 | +for it to work in the solvers. |
| 7 | + |
| 8 | +## AbstractDiffEqOperator Interface Description |
| 9 | + |
| 10 | +1. Function call and multiplication: `L(du,u,p,t)` for inplace and `du = L(u,p,t)` for |
| 11 | + out-of-place, meaning `L*u` and `mul!`. |
| 12 | +2. If the operator is not a constant, update it with `(u,p,t)`. A mutating form, i.e. |
| 13 | + `update_coefficients!(A,u,p,t)` that changes the internal coefficients, and a |
| 14 | + out-of-place form `B = update_coefficients(A,u,p,t)`. |
| 15 | +3. `isconstant(A)` trait for whether the operator is constant or not. |
| 16 | + |
| 17 | +## AbstractDiffEqLinearOperator Interface Description |
| 18 | + |
| 19 | +1. `AbstractSciMLLinearOperator <: AbstractSciMLOperator` |
| 20 | +2. Can absorb under multiplication by a scalar. In all algorithms things like |
| 21 | + `dt*L` show up all the time, so the linear operator must be able to absorb |
| 22 | + such constants. |
| 23 | +4. `isconstant(A)` trait for whether the operator is constant or not. |
| 24 | +5. Optional: `diagonal`, `symmetric`, etc traits from LinearMaps.jl. |
| 25 | +6. Optional: `exp(A)`. Required for simple exponential integration. |
| 26 | +7. Optional: `expv(A,u,t) = exp(t*A)*u` and `expv!(v,A::AbstractSciMLOperator,u,t)` |
| 27 | + Required for sparse-saving exponential integration. |
| 28 | +8. Optional: factorizations. `ldiv!`, `factorize` et. al. This is only required |
| 29 | + for algorithms which use the factorization of the operator (Crank-Nicolson), |
| 30 | + and only for when the default linear solve is used. |
| 31 | + |
| 32 | +## Note About Affine Operators |
| 33 | + |
| 34 | +Affine operators are operators which have the action `Q*x = A*x + b`. These operators have |
| 35 | +no matrix representation, since if there was it would be a linear operator instead of an |
| 36 | +affine operator. You can only represent an affine operator as a linear operator in a |
| 37 | +dimension of one larger via the operation: `[A b] * [u;1]`, so it would require something modified |
| 38 | +to the input as well. As such, affine operators are a distinct generalization of linear operators. |
| 39 | + |
| 40 | +While it this seems like it might doom the idea of using matrix-free affine operators, it turns out |
| 41 | +that affine operators can be used in all cases where matrix-free linear solvers are used due to |
| 42 | +an easy genearlization of the standard convergence proofs. If Q is the affine operator |
| 43 | +``Q(x) = Ax + b``, then solving ``Qx = c`` is equivalent to solving ``Ax + b = c`` or ``Ax = c-b``. |
| 44 | +If you know do this same "plug-and-chug" handling of the affine operator in into the GMRES/CG/etc. |
| 45 | +convergence proofs, move the affine part to the rhs residual, and show it converges to solving |
| 46 | +``Ax = c-b``, and thus GMRES/CG/etc. solves ``Q(x) = c`` for an affine operator properly. |
| 47 | + |
| 48 | +That same trick then can be used pretty much anywhere you would've had a linear operator to extend |
| 49 | +the proof to affine operators, so then ``exp(A*t)*v`` operations via Krylov methods work for A being |
| 50 | +affine as well, and all sorts of things. Thus affine operators have no matrix representation but they |
| 51 | +are still compatible with essentially any Krylov method which would otherwise be compatible with |
| 52 | +matrix-free representations, hence their support in the SciMLOperators interface. |
0 commit comments