`ConstrainedSystems.ConstrainedODEFunction`

— Type` ConstrainedODEFunction(r1,r2,B1,B2[,L][,C])`

This specifies the functions and operators that comprise an ODE problem with the form

$\dfrac{dy}{dt} = Ly - B_1 z + r_1(y,t)$

$B_2 y + C z = r_2(x,t)$

where $y$ is the state, $z$ is a constraint force, and $x$ is an auxiliary state describing the constraints.

The optional linear operator `L`

defaults to zeros. The `B1`

and `B2`

functions must be of the respective in-place forms `B1(dy,z,x,p)`

(to compute the action of `B1`

on `z`

) and `B2(dz,y,x,p)`

(to compute the action of `B2`

on `y`

). The function `r1`

must of the in-place form `r1(dy,y,x,p,t)`

, and `r2`

must be in the in-place form `r2(dz,x,p,t)`

. The `C`

function can be omitted, but if it is included, then it must be of the form `C(dz,z,x,p)`

(to compute the action of `C`

on `z`

). Alternatively, one can supply out-of-place forms, respectively, as `B1(z,x,p)`

, `B2(y,x,p)`

, `C(z,x,p)`

, `r1(y,x,p,t)`

and `r2(x,p,t)`

.

An optional keyword argument `param_update_func`

can be used to set a function that updates problem parameters with the current solution. This function must take the in-place form `f(q,u,p,t)`

or out of place form `f(u,p,t)`

to create some `q`

based on `u`

, where `y = state(u)`

, `z = constraint(u)`

and `x = aux_state(u)`

. (Note that `q`

might enter the function simply as `p`

, to be mutated.) This function can be used to update `B1`

, `B2`

, and `C`

, for example.

We can also include another (unconstrained) set of equations to the set above in order to update `x`

:

$\dfrac{dx}{dt} = r_{1x}(u,p,t)$

In this case, the right-hand side has access to the entire `u`

vector. We would pass the pair of `r1`

functions as an `ArrayPartition`

.

`ConstrainedSystems.SaddleSystem`

— Method`SaddleSystem`

Construct a saddle-point system operator from the constituent operator blocks. The resulting object can be used with `*`

and `\`

to multiply and solve. The saddle-point problem has the form

$\begin{bmatrix}A & B_1^T \\ B_2 & C \end{bmatrix} \begin{pmatrix} u \\ f \end{pmatrix} = \begin{pmatrix} r_1 \\ r_2 \end{pmatrix}$

**Constructors**

`SaddleSystem(A::AbstractMatrix,B₂::AbstractMatrix,B₁ᵀ::AbstractMatrix,C::AbstractMatrix[,eltype=Float64])`

. Blocks are given as matrices. Must have consistent sizes to stack appropriately. If this is called with `SaddleSystem(A,B₂,B₁ᵀ)`

, it sets `C`

to zero automatically.

`SaddleSystem(A,B₂,B₁ᵀ,C,u,f[,eltype=Float64])`

. Operators `A`

, `B₂`

, `B₁ᵀ`

, `C`

are given in various forms, including matrices, functions, and function-like objects. `u`

and `f`

are examples of the data types in the corresponding solution and right-hand side vectors. Guidelines:

- The entries
`A`

and`B₂`

must be able to act upon`u`

(either by multiplication or as a function) and`B₁ᵀ`

and`C`

must be able to act on`f`

(also, either by multiplication or as a function). `A`

and`B₁ᵀ`

should return data of type`u`

, and`B₂`

and`C`

should return data of type`f`

.`A`

must be invertible and be outfitted with operators ``and`

ldiv!`.- Both
`u`

and`f`

must be subtypes of`AbstractArray`

: they must be equipped with`size`

and`vec`

functions and with a constructor of the form`T(data)`

where`T`

is the data type of`u`

or`f`

and`data`

is the wrapped data array.

If called as `SaddleSystem(A,B₂,B₁ᵀ,u,f)`

, the `C`

block is omitted and assumed to be zero.

If called with `SaddleSystem(A,u)`

, this is equivalent to calling `SaddleSystem(A,nothing,nothing,u,[])`

, then this reverts to the unconstrained system described by operator `A`

.

The list of vectors `u`

and `f`

in any of these constructors can be bundled together as a `SaddleVector`

, e.g. `SaddleSystem(A,B₂,B₁ᵀ,SaddleVector(u,f))`

.

An optional keyword argument `solver=`

can be used to specify the type of solution for the Schur complement system. By default, this is set to `Direct`

, and the Schur complement matrix is formed, factorized, and stored. This can be changed to a variety of iterative solvers, e.g. `BiCGStabl`

, `CG`

, `GMRES`

, in which case an iterative solver from `IterativeSolvers.jl`

is used.

`ConstrainedSystems.SaddleVector`

— Type`SaddleVector(u,f)`

Construct a vector of a state part `u`

and constraint part `f`

of a saddle-point vector, to be associated with a `SaddleSystem`

.

`Base.eltype`

— Method`Base.eltype(::SaddleSystem)`

Report the element type of a `SaddleSystem`

.

`Base.size`

— Method`Base.size(::SaddleSystem)`

Report the size of a `SaddleSystem`

.

`ConstrainedSystems.aux_state`

— Method`aux_state(x)`

Provide the auxiliary state part of the given vector `x`

`ConstrainedSystems.constraint`

— Method`constraint(x::SaddleVector)`

Provide the constraint part of the given saddle vector `x`

`ConstrainedSystems.r1vector`

— Method`r1vector([;state_r1=][,aux_r1=])`

Build a vector of the `r1`

functions for the state ODEs and auxiliary state ODEs.

`ConstrainedSystems.solvector`

— Method`solvector([;state=][,constraint=][,aux_state=])`

Build a solution vector for a constrained system. This takes three optional keyword arguments: `state`

, `constraint`

, and `aux_state`

. If only a state is supplied, then the constraint is set to an empty vector and the system is assumed to correspond to an unconstrained system. (`aux_state`

is ignored in this situation.)

`ConstrainedSystems.state`

— Method`state(x::SaddleVector)`

Provide the state part of the given saddle vector `x`