API reference

TSRK(subtype = "5dp", petsc_options = String[]; comm = MPI.COMM_SELF, dm = nothing)

Explicit Runge-Kutta from PETSc's TSRK. subtype is a PETSc TSRKType without its prefix, such as "3bs", "5dp", "5f" or "5bs".

Adapts on its embedded error estimate, so reltol and abstol apply. PETSc gives "1fe", "2b", "3" and "4" no such estimate, so those step at the dt you give and warn if you pass a tolerance. Being explicit it never forms a Jacobian and ignores an ODEFunction's jac, and it cannot carry a mass matrix.

petsc_options are command-line style tokens passed to PETSc for this solve, for example ["-ts_adapt_type", "none"]. They are parsed after the options this package sets, so they win.

A comm other than MPI.COMM_SELF runs the solve distributed over it, with u0 holding this rank's rows; see the MPI section of the documentation. A dm, a DMDA from PETSc.jl, runs it on the DM's communicator and gives f the state ghosted from the DM, as that section describes.

TSRosW(subtype = "ra34pw2", petsc_options = String[]; autodiff = AutoForwardDiff(), comm = MPI.COMM_SELF, dm = nothing)

Rosenbrock-W from PETSc's TSROSW. subtype is a PETSc TSRosWType without its prefix, such as "2m", "ra34pw2" or "r34prw".

Adapts on its embedded error estimate, except for "theta1" and "theta2", which PETSc gives none, so they step at the dt you give and warn if you pass a tolerance. Linearly implicit, so it uses an ODEFunction's jac, and it accepts a mass matrix.

Without a jac the Jacobian comes from autodiff: ForwardDiff by default, colouring a sparse jac_prototype, or AutoFiniteDiff() to have PETSc difference the step's own equations, colouring a sparse prototype too. For a complex state f has to be holomorphic, and one that is not is refused under ForwardDiff.

PETSc's implementation assumes a right-hand side that does not depend on t. When it does, "2p", "2m", "ra3pw", "ra34pw2", "r34prw" and "assp3p3s1c" keep their order, and every other type, including all the fourth-order ones, converges at first order with or without a jac. Carrying t as an extra state whose derivative is 1 restores their order.

On a linear right-hand side that does not depend on t, "ra3pw"'s embedded error estimate is zero with an exact Jacobian and far too small with a differenced one, so an adaptive solve reports success with an error well above the tolerance. Step it with a fixed dt on such problems.

"assp3p3s1c" needs a Jacobian, which AutoFiniteDiff() does not give it, and cannot take a mass matrix, which PETSc leaves out of its explicit first stage. "lassp3p4s2c", "llssp3p4s2c" and "ark3" are refused. They end on an explicit stage: without a Jacobian PETSc stops and asks for one, and with one it does not restore its Jacobian lag after that stage, so an adaptive solve fails within its first two steps and a fixed-step solve diverges.

A comm other than MPI.COMM_SELF runs the solve distributed over it, as for TSRK. There autodiff defaults to AutoFiniteDiff(), and a jac fills this rank's rows of a sparse jac_prototype whose columns are global; see the MPI section of the documentation. With a dm the Jacobian is the DM's own matrix, with the pattern of its stencil, which a jac fills through PETSc's matrix API from the ghosted u, or which PETSc colours and differences f into when there is none, with autodiff then at its default there, AutoFiniteDiff().

TSImplicit(subtype = "beuler"; order = nothing, autodiff = AutoForwardDiff(), comm = MPI.COMM_SELF, dm = nothing)
TSImplicit(subtype, theta; ...)
TSImplicit(subtype, [theta,] petsc_options; ...)

Fully implicit methods from PETSc: "beuler", "cn", "theta" and "bdf". theta sets the parameter of the theta method, where 0.5 is Crank-Nicolson and 1.0 is backward Euler.

order sets the BDF order, 1 through 6. PETSc's own default is 2, which on a stiff problem can cost an order of magnitude in steps against a higher-order method, so raise it when comparing against one.

Only "bdf" carries an embedded error estimate and adapts; the others step at the dt you give and warn if you pass a tolerance. All of them use an ODEFunction's jac and accept a mass matrix, which makes a singular mass matrix an index-1 differential-algebraic problem.

Without a jac the Jacobian comes from autodiff: ForwardDiff by default, colouring a sparse jac_prototype, or AutoFiniteDiff() to have PETSc difference the step's own equations, colouring a sparse prototype too. For a complex state f has to be holomorphic, since PETSc's Newton iteration takes a complex Jacobian, and one that is not is refused under ForwardDiff.

A comm other than MPI.COMM_SELF runs the solve distributed over it, as for TSRK. There autodiff defaults to AutoFiniteDiff(), and a jac fills this rank's rows of a sparse jac_prototype whose columns are global; see the MPI section of the documentation. With a dm the Jacobian is the DM's own matrix, with the pattern of its stencil, which a jac fills through PETSc's matrix API from the ghosted u, or which PETSc colours and differences f into when there is none, with autodiff then at its default there, AutoFiniteDiff().

TSIRK(nstages = 3, petsc_options = String[]; autodiff = AutoForwardDiff(), comm = MPI.COMM_SELF, dm = nothing)

Gauss-Legendre implicit Runge-Kutta from PETSc's TSIRK, of order 2 * nstages: one stage is the implicit midpoint rule at order 2, two stages give order 4 and three give order 6. Measured at each of those.

Fixed step: PETSc gives this family no embedded error estimate, so it steps at the dt you give and warns if you pass a tolerance.

Needs a Jacobian, from the ODEFunction's jac or from autodiff, and refuses AutoFiniteDiff() rather than letting PETSc fail, since it solves all stages as one coupled system whose matrix it cannot build from finite differences. That coupled matrix is a Kronecker product with the Jacobian, which has no LU factorisation, so this algorithm defaults to -pc_type pbjacobi unless your own petsc_options set a -pc_type. So does "irk" picked by TSGeneric or by -ts_type.

A wrong Jacobian is not caught here. Where the other implicit families fail to converge, this one reports success and returns a wrong answer, so check a hand-written jac against a solve without one before trusting it.

A mass matrix is rejected. PETSc's coupled-stage matrix assumes dF/du̇ = I, and with a non-identity mass matrix the answer drifts further from the true one as dt shrinks instead of failing, which is worse than an error.

A comm other than MPI.COMM_SELF runs the solve distributed over it, as for TSRK. There it needs a jac, filling this rank's rows of a sparse jac_prototype whose columns are global, and each rank has to hold PETSc's own share of the state, which splits it evenly with the first ranks taking one row more; see the MPI section of the documentation. A dm is refused.

TSARKIMEX(subtype = "3", petsc_options = String[]; autodiff = AutoForwardDiff(), comm = MPI.COMM_SELF, dm = nothing)

Additive Runge-Kutta IMEX from PETSc's TSARKIMEX. subtype is a PETSc TSARKIMEXType without its prefix, such as "2e", "3", "4" or "5".

Takes a SplitODEProblem whose f1 is integrated implicitly and whose f2 is integrated explicitly, and uses f1's Jacobian when the problem carries one. Without one, f1's Jacobian comes from autodiff, as for TSImplicit. A plain ODEProblem is treated as fully implicit with the explicit part left at zero, which is PETSc's own default. Adapts on its embedded error estimate, except for "prssp2", "ars443" and "bpr3", which PETSc gives none, so they step at the dt you give and warn if you pass a tolerance.

"ars122" needs a SplitODEProblem: it has an explicit first stage and is not stiffly accurate, so PETSc cannot evaluate its first-stage slope when the whole problem is implicit.

"bpr3" is refused on a SplitODEProblem, where it converges at first order. On a plain ODEProblem PETSc does not use its explicit tableau, and it keeps order 3.

A comm other than MPI.COMM_SELF runs the solve distributed over it, as for TSRK. There autodiff defaults to AutoFiniteDiff(), and a jac fills this rank's rows of a sparse jac_prototype whose columns are global; see the MPI section of the documentation. With a dm the Jacobian is the DM's own matrix, with the pattern of its stencil, which a jac fills through PETSc's matrix API from the ghosted u, or which PETSc colours and differences f into when there is none, with autodiff then at its default there, AutoFiniteDiff().

TSDAE(subtype = "bdf", petsc_options = String[]; order = nothing, autodiff = AutoForwardDiff(), comm = MPI.COMM_SELF, dm = nothing)

Fully implicit methods applied to a DAEProblem, whose residual G(t, u, u') = 0 is exactly the form PETSc's IFunction takes. subtype is "beuler", "cn", "theta" or "bdf", as for TSImplicit.

A DAEFunction's jac(J, du, u, p, gamma, t) is gamma * dG/du' + dG/du, which is what PETSc's IJacobian wants whole, so it is passed straight through and gamma is PETSc's shift. Without one the Jacobian comes from autodiff, as for TSImplicit.

order sets the BDF order, 1 through 6, and carries the same warning as TSImplicit: PETSc's own default is 2.

Only "bdf" adapts; the others step at the dt you give. PETSc derives the initial derivative itself, so du0 is used only to check and solve for a consistent start, as initializealg asks; see the DAE initialization section of the documentation.

A comm other than MPI.COMM_SELF runs the solve distributed over it, as for TSRK. There autodiff defaults to AutoFiniteDiff(), and a jac fills this rank's rows of a sparse jac_prototype whose columns are global; see the MPI section of the documentation. With a dm the Jacobian is the DM's own matrix, with the pattern of its stencil, which a jac fills through PETSc's matrix API from the ghosted u, or which PETSc colours and differences f into when there is none, with autodiff then at its default there, AutoFiniteDiff().

TSMPRK(slow, subtype = "p2", petsc_options = String[]; comm = MPI.COMM_SELF, dm = nothing)
TSMPRK(slow, medium, subtype = "2a23", petsc_options = String[]; comm = MPI.COMM_SELF, dm = nothing)

PETSc's multirate partitioned Runge-Kutta. slow lists the indices of the state that are integrated with the outer step; everything else is advanced on a smaller one. Both parts come from the problem's own f, which is evaluated in full and then read row by row, so nothing beyond the index list is asked of the caller.

Passing a medium set as well splits the state three ways, which is what PETSc's "2a23" and "2a33" want. Whatever neither set names is the fast part.

Explicit, so a mass matrix and an analytic Jacobian are both refused. Without a medium set subtype is one of "2a22", "2a32", "p2" or "p3"; with one it is "2a23" or "2a33".

The fast part is stepped on a shorter clock, so the step this tolerates is larger than an explicit single-rate method's: on u' = [-u1, -100u2] with only the first component slow, "p3" holds to dt = 0.0425 against TSRK("5dp")'s 0.03. The ratio is fixed by the tableau rather than by the stiffness, so a separation much wider than that is not something these methods can absorb.

A comm other than MPI.COMM_SELF runs the solve distributed over it, as for TSRK. There slow and medium index this rank's rows and may be empty on some ranks, as long as some rank names a slow row and, for "2a23" and "2a33", a medium one. A dm is refused.

TSBasicSymplectic(subtype = "velverlet", petsc_options = String[]; comm = MPI.COMM_SELF)

PETSc's basic symplectic integrators, TSBASICSYMPLECTIC, for a DynamicalODEProblem or a SecondOrderODEProblem. subtype is "sieuler" (or "1"), semi-implicit Euler of order 1, "velverlet" (or "2"), velocity Verlet of order 2, "3", Ruth's method of order 3, or "4", Forest and Ruth's method of order 4.

PETSc steps the two parts in turn, the velocity v with f1(dv, v, u, p, t) at the current positions and the position u with f2(du, v, u, p, t) at the current velocities, so, as for OrdinaryDiffEq's symplectic methods, f1 must not depend on v nor f2 on u. On such a separable problem the energy error stays bounded over long times instead of drifting. f1 is given the time of the positions it is evaluated at, so a force that depends on t keeps the method's order.

Fixed step: these methods have no error estimate, so they step at the dt you give and warn if you pass a tolerance. Being explicit they ignore a jac and refuse a mass matrix.

Distributed partitioned states are not supported yet, so a comm other than MPI.COMM_SELF is refused.

TSAlpha2(petsc_options = String[]; radius = nothing, autodiff = AutoForwardDiff(), comm = MPI.COMM_SELF)

PETSc's generalized-alpha method for second-order systems, TSALPHA2, for a SecondOrderODEProblemu'' = f(u', u, p, t). Implicit and of order 2. radius is the spectral radius at an infinite step, from 0 to 1: PETSc's default, 1, damps nothing, and a smaller one damps the highest frequencies more.

Adapts on PETSc's error estimate for the method, which PETSc leaves off unless asked, so reltol and abstol apply. PETSc weighs the velocity and the position with one tolerance per component of the position, so a vector tolerance on the state [v; u] takes the smaller of each velocity and position pair. adaptive = false steps at the dt you give.

Each step solves for the acceleration, with the Jacobian of f in u and u'. That comes from a jac, which is the Jacobian of the first-order system [v; u]' = [f(v, u, p, t); v], 2n by 2n in the order of the state, as OrdinaryDiffEq's implicit methods take it, or without one from autodiff, as for TSImplicit. A sparse jac_prototype of that system keeps the n by n matrix PETSc factors sparse: a jac fills the prototype's structure, and without one the Jacobian is coloured, by SparseMatrixColorings under AD or by PETSc under AutoFiniteDiff, as for TSImplicit. A mass matrix is refused. Until SciMLBase's DynamicalODEFunction can be rebuilt with a jac_prototype, solve and init fail on one before reaching this package, so give such a problem to SciMLBase.__solve or SciMLBase.__init.

A reversed tspan is supported. PETSc only steps forward, so the solve runs in s = -t on w(s) = u(-s), whose velocity is -u' and whose acceleration is f(-w', w, p, -s). The states, the saved velocities, jac and the callbacks stay those of the problem as written.

Distributed partitioned states are not supported yet, so a comm other than MPI.COMM_SELF is refused.

TSGeneric(ts_type, petsc_options = String[]; explicit = false, autodiff = AutoForwardDiff(), comm = MPI.COMM_SELF, dm = nothing)

Any other PETSc TSType by name. An implicit one such as "alpha" works with the default; an explicit one such as "euler" or "ssp" needs explicit = true, since PETSc then wants the right-hand side rather than the implicit residual. The constructor refuses a type given the wrong explicit. An explicit type also ignores a jac and rejects a mass matrix. An implicit one without a jac gets its Jacobian from autodiff, as for TSImplicit.

A type another constructor covers, such as "bdf" or "beuler", is treated as that constructor treats it. Whether any other type adapts is not known here, so it needs dt and gets no tolerance warning. Only "euler" and "alpha" have been run through this package's own convergence tests.

"discgrad", "eimex", "mimex" and "mprk" are refused: each is driven through a PETSc setup call this package does not make, and without it they crash or integrate to zero rather than saying anything. "alpha2" and "basicsymplectic" are refused too, since they need a second-order or partitioned problem; use TSAlpha2 or TSBasicSymplectic.

A comm other than MPI.COMM_SELF runs an explicit type as TSRK runs, and an implicit "beuler", "cn", "theta", "bdf", "rosw", "arkimex", "irk", "alpha" or "dirk" as TSImplicit runs, with autodiff defaulting to AutoFiniteDiff(). Other implicit types are refused there: "glle"'s step control follows the round-off of the distributed linear solve, so it takes other steps than a serial solve and ends with another error, larger or smaller. A dm runs an explicit type as it runs TSRK and refuses an implicit one.

PETScIntegrator

The integrator SciMLBase.init returns for a PETSc TS algorithm. Step it with step!, run it to the end with solve!, stop it early with terminate! and restart it with reinit!. Between steps u, uprev, t, tprev and dt are readable, integ(t), integ(t, Val{1}), integ(t; idxs) and integ(out, t) give the state or its slope inside the step just taken, and add_tstop! schedules a time to land on exactly.

These are in the types PETSc steps in. The clock, and so t, dt and the saved times, is Float32 for a Float32 or ComplexF32 state with a Float32 span and Float64 otherwise, whatever the span's type. The state is the problem's own, except that a whole-number one is stepped in Float64 and a single-precision one with a Float64 span in double precision, while the solution it saves stays in single. For a DynamicalODEProblem or SecondOrderODEProblem the state is an ArrayPartition of the velocity and the position, as in OrdinaryDiffEq.

iter counts every step attempted, accepted, rejected or failed, as maxiters does. opts carries the options OrdinaryDiffEq's integrator exposes where PETSc has them. Assigning maxiters, save_everystep, save_start, save_end, unstable_check or isoutofdomain takes effect at once; tstops and saveat are the live queues keyed by tdir * t, with saveat keeping the times already passed; unstable_check and isoutofdomain are nothing unless given; dense, save_idxs, calck, internalnorm and callback are fixed at init.

PETScAdjoint(; petsc_options = String[])

PETSc's own discrete adjoint, TSAdjointSolve, as a sensealg for SciMLSensitivity's adjoint_sensitivities:

using PETScDiffEq, SciMLSensitivity
sol = solve(prob, TSRK("4"); dt = 0.01, adaptive = false)
du0, dp = adjoint_sensitivities(
    sol, TSRK("4"); sensealg = PETScAdjoint(),
    t = ts, dgdu_discrete = dg!, dt = 0.01, adaptive = false,
)

The result is the gradient of the solution PETSc computes at these steps, not of the exact solution, so it agrees with finite differences of the same fixed-step solve. PETSc runs the forward solve again, saving a trajectory, and then the adjoint. Only the problem is taken from sol, so the keywords that set the steps (dt, adaptive, abstol, reltol, dtmin, dtmax, maxiters) must be repeated exactly as solve was given them; otherwise the gradient belongs to a different discretization and nothing says so. Saving keywords are ignored, and callback, tstops, d_discontinuities and save_idxs are refused, whether given here or to the problem.

Supported: an ODEProblem without a mass matrix, in place or out of place, solved with TSRK of any subtype, TSImplicit("beuler"), TSImplicit("cn"), TSImplicit("theta") with its theta, in its midpoint form or with -ts_theta_endpoint, or TSARKIMEX, in either time direction. The ODEFunction's jac and paramjac are used when given, and otherwise built as the forward solve builds a missing jac: with the algorithm's autodiff, and with ForwardDiff for a TSRK. Under AutoFiniteDiff() both have to be given, since PETSc's own differences never reach its adjoint. A sparse jac_prototype is used. paramjac takes the form jac does, paramjac(pJ, u, p, t) in place or paramjac(u, p, t) out of place, with a row per state and a column per entry of p, which therefore has to be a vector of real numbers. A hand-written jac or paramjac goes into the gradient unchecked, so a wrong one gives a wrong gradient without an error; compare it against a gradient computed without it.

TSARKIMEX takes a SplitODEProblem as well, on MPI.COMM_SELF. The implicit part's jac and paramjac are the problem's, which a SplitODEProblem takes from f1, and the explicit part's are those of f2's own ODEFunction; each is built by automatic differentiation when missing. PETSc's ARKIMEX adjoint has no quadrature, so an integral cost is refused, and -ts_arkimex_fully_implicit is refused on a SplitODEProblem, since the adjoint would still take f2 explicitly. On a plain ODEProblem, every type but "1bee", "l2" and "prssp2" has an explicit first stage, which PETSc evaluates at a stale time on the first step after a restart without its adjoint seeing that. Such a type therefore needs tspan to start at 0 and the stages kept in the trajectory, and is refused otherwise; a SplitODEProblem has no such limit.

A DynamicalODEProblem or SecondOrderODEProblem is differentiated as the first-order system on the flat [v; u] that these methods step. Its jac and paramjac, when given, are that system's, taking the state as an ArrayPartition(v, u) as the forward solve's jac does, and are otherwise built from f1 and f2. The cost functions are handed the state as an ArrayPartition and write its derivative into one, and du0 comes back as one. PETSc has no adjoint for TSBasicSymplectic or TSAlpha2, so those are refused.

The adjoint runs in PETSc's double real build. A Float32 problem is solved there in Float64, so jac, paramjac and the cost functions are handed Float64 states, and du0 and dp come back as Float32 where u0 and p are. Its cost times are matched to the steps at single precision, so the times its own solve saved are accepted with dt repeated as that solve was given it. A complex state is refused.

Costs are discrete, integral or both. At each t[i], dgdu_discrete(out, u, p, t, i) writes the cost's derivative with respect to the state and dgdp_discrete(out, u, p, t, i), if given, its direct derivative with respect to p; no_start = true leaves out t[1]. PETSc's adjoint has no derivative of interpolation, so with fixed steps every cost time must be a time the solve steps to, such as tspan[1] plus a multiple of dt. An adaptive solve lands only on the ends of tspan, so its costs are limited to those, and its gradient treats the accepted step sizes as constants rather than differentiating the controller.

An integral cost, the integral of g(u, p, t) from tspan[1] to tspan[2], is given as g or as its derivatives dgdu_continuous(out, u, p, t) and dgdp_continuous(out, u, p, t). PETSc sums it with the method's own quadrature, through a quadrature TS whose right-hand side Jacobians are the derivatives of g, so the gradient is that of the sum of dt * b[i] * g over the stages of a TSRK, of dt * g at the end of each backward Euler step, of the trapezoidal sum for Crank-Nicolson, and for the theta method of dt * g at its stage, at t + theta * dt, or in its endpoint form of dt * ((1 - theta) * g(t) + theta * g(t + dt)). Where dgdu_continuous is left out, g is differentiated with respect to u and p as a missing jac is, with ForwardDiff or the algorithm's autodiff; where it is given and dgdp_continuous is not, the direct derivative with respect to p is taken as zero, as SciMLSensitivity's adjoints take it. Under AutoFiniteDiff() the derivatives have to be given. Given together with discrete terms, the two gradients are added. An integral cost is refused on a communicator other than MPI.COMM_SELF.

petsc_options apply to this adjoint's own run and are parsed after the algorithm's. The trajectory is kept in memory, every stage of every step; -ts_trajectory_solution_only 1 keeps only the states and recomputes each step during the adjoint, and -ts_trajectory_type basic writes one file per step to the working directory instead. For an adaptive solve the memory trajectory reserves 8 bytes for each of maxiters steps before starting, 8 MB at the default and 8 GB at maxiters = 10^9. TSImplicit and TSARKIMEX solve transposed linear systems with the same Krylov solver and tolerances as their Newton steps, so with default options the gradient can be off by up to about their relative tolerance of 1e-5 while the forward states are far closer; pass -ksp_type preonly -pc_type lu, or a tighter -ksp_rtol, when that matters. Options PETSc reads only while the adjoint runs, such as -ts_trajectory_view and -ts_adjoint_view_solution, have no effect in petsc_options, though they do when set globally, for example through PETSC_OPTIONS.

With an algorithm whose comm is not MPI.COMM_SELF, the adjoint runs distributed as the solve does, for the same methods on an ODEProblem. The ODEFunction then needs jac, filling this rank's rows of a sparse jac_prototype with global columns, and paramjac when there are parameters, filling this rank's rows; both are collective like f. dgdu_discrete gets this rank's rows of the state and writes their derivative, and dgdp_discrete gives this rank's share of the direct derivative, which the ranks add up. du0 holds this rank's rows, and dp is the whole gradient on every rank. The cost times, no_start, the length of p and whether dgdp_discrete is given must agree across the ranks.

Returns (du0, dp'), where dp is nothing when p is nothing or SciMLBase.NullParameters(). Differentiating solve itself with a reverse-mode AD package is not supported.

PETScDiffEq.reshape_local_array(x, dm)

This rank's part of a vector on the DMDA dm, indexed by grid point in global numbering: a[c, i] on a 1-D grid and a[c, i, j] on a 2-D one, where c runs over the degrees of freedom at each point. x can be the ghosted array f receives or the block of the state the rank owns, such as du, u0 or a saved state, and a shares its memory. This is PETSc.jl's reshape_local_array, which PETSc.jl 0.4 calls reshapelocalarray and pads to three grid axes.

set_stencil_values!(J, rows, cols, vals; add = false)

Write a block of the DM's matrix J, as a jac gets it in a solve with a dm, through PETSc's MatSetValuesStencil. rows and cols are grid indices, or vectors or tuples of them, in the numbering PETScDiffEq.reshape_local_array uses: (c, i) on a 1-D grid, (c, i, j) on a 2-D one and (c, i, j, k) on a 3-D one, as tuples or CartesianIndexes, global and 1-based, with c the degree of freedom at the point. vals[r, s] is the entry at rows[r] and cols[s], a vector or tuple when either side is a single index and a number when both are. The rows are this rank's own and the columns lie within its ghost region; PETSc drops an index past a ghosted edge, which has no global entry, and refuses an entry outside the DM's stencil. add = true adds to the entries instead of setting them, and PETSc wants a PETSc.assemble!(J) between a set and an add. It allocates nothing once warm, so with tuples for cols and vals a jac can call it at every grid point without allocating.