TS

The TS (Time Stepping) module integrates ordinary differential equations and differential algebraic equations in time, and is the layer PETSc builds its implicit, explicit and IMEX methods on.

The interface described here is the high-level one. It owns the @cfunction trampolines, keeps your callbacks rooted for as long as the stepper lives, and frees the object for you on a sequential communicator. The low-level TS bindings remain available and unchanged.

Overview

A problem is given to TS in one of two forms, or in both at once:

Explicit, $du/dt = G(t, u)$, through PETSc.set_rhs_function! and optionally PETSc.set_rhs_jacobian!.
Implicit, $F(t, u, du/dt) = 0$, through PETSc.set_ifunction! and PETSc.set_ijacobian!. This form also covers a DAE and a problem with a mass matrix.

Setting both is how an IMEX method splits a problem into the stiff part it treats implicitly and the rest it treats explicitly.

Creating a stepper

ts = PETSc.TS(petsclib, MPI.COMM_WORLD)

# or with options, which are applied when you call solve!
ts = PETSc.TS(petsclib, MPI.COMM_WORLD;
    ts_type = "bdf",
    ts_adapt_type = "basic",
)

Options given here are held and applied inside PETSc.solve! rather than at construction, so a DM and the callbacks attached afterwards are in place before PETSc reads them. Until then the object has no type, and PETSc.type_name answers nothing.

On a communicator of size 1 the garbage collector calls PETSc.destroy!. On a larger one, destruction is yours to do, since collection is asynchronous and the call is collective.

An explicit problem

Solving $du/dt = -u$ from $u(0) = 1$:

ts = PETSc.TS(petsclib, MPI.COMM_SELF)
PETSc.set_type!(ts, :rk)

u = PETSc.PetscVec(petsclib, 1)
u[1] = 1.0
PETSc.assemble!(u)

PETSc.set_rhs_function!(ts) do F, ts, t, u
    PETSc.with_local_array!((u, F); read = (true, false), write = (false, true)) do ua, Fa
        Fa[1] = -ua[1]
    end
end

PETSc.set_time!(ts, 0.0)
PETSc.set_timestep!(ts, 0.01)
PETSc.set_max_time!(ts, 1.0)

PETSc.solve!(u, ts)

The callback comes first in the argument list, so do block syntax works. Pass it as a plain argument, set_rhs_function!(f!, ts), when it is already a value: the subject-first order v0.4 also accepted is gone, and has no shim.

An implicit problem

The same equation written as $F(t, u, u_t) = u_t + u = 0$, stepped with backward Euler. The Jacobian asked for is $\partial F/\partial u + \sigma\, \partial F/\partial u_t$, where PETSc supplies the shift $\sigma$:

ts = PETSc.TS(petsclib, MPI.COMM_SELF)
PETSc.set_type!(ts, :beuler)

J = PETSc.PetscMat(petsclib, 1, 1, petsclib.PetscInt(1))

PETSc.set_ifunction!(ts) do F, ts, t, u, u_t
    PETSc.with_local_array!(
        (u, u_t, F);
        read = (true, true, false),
        write = (false, false, true),
    ) do ua, uta, Fa
        Fa[1] = uta[1] + ua[1]
    end
end

PETSc.set_ijacobian!(ts, J) do A, P, ts, t, u, u_t, shift
    A[1, 1] = shift + 1
    PETSc.assemble!(A)
end

A Jacobian callback is always handed both the Jacobian A and the preconditioning matrix P. They are the same object unless you passed two, which is worth testing with A.ptr == P.ptr before filling P a second time.

Callback signatures

Setter	Callback
`PETSc.set_rhs_function!`	`f!(F, ts, t, u)`
`PETSc.set_rhs_jacobian!`	`updateJ!(A, P, ts, t, u)`
`PETSc.set_ifunction!`	`f!(F, ts, t, u, u_t)`
`PETSc.set_ijacobian!`	`updateJ!(A, P, ts, t, u, u_t, shift)`
`PETSc.set_monitor!`	`f(ts, step, t, u)`
`PETSc.set_pre_step!`	`f!(ts)`
`PETSc.set_post_step!`	`f!(ts)`

A callback's return value is ignored; it reports a failure by throwing. Until v0.6, a nonzero Integer returned by one of the first five still fails the call, with a warning.

An exception thrown in a callback never crosses into PETSc's C code. PETSc is told the callback failed, unwinds, and PETSc.solve! or PETSc.step! then rethrows the original exception.

Precompilation

Each setter builds its @cfunction trampoline when you call it, so nothing here is affected by precompilation, and there is nothing you need to do about it.

It is worth knowing where the hazard is, though, because the low-level examples show the other pattern. A @cfunction yields a pointer valid only for the session that evaluated it, so this at the top level of a package:

const MY_RHS_PTR = @cfunction(my_rhs!, PetscErrorCode, (CTS, PetscReal, CVec, CVec, Ptr{Cvoid}))

captures a pointer during precompilation and hands PETSc a stale one in every later session. In a script such as examples/ex16.jl it is fine, since the file is evaluated afresh each run. Inside a package it is not: build the pointer inside the function that registers it, or assign it from __init__.

Passing your own data

Anything stored with PETSc.set_user_ctx! is handed back as a trailing argument, to whichever callbacks have a method that accepts one:

PETSc.set_user_ctx!(ts, (; viscosity = 1e-3))

PETSc.set_rhs_function!(ts) do F, ts, t, u, ctx
    # ctx.viscosity is available here
end

The object is held on the Julia side by ts, so it stays alive without any pinning of your own.

Watching a solve

PETSc.set_monitor!(ts) do ts, step, t, u
    @printf("step %3d  t = %.4f\n", step, t)
end

The monitor runs once before the first step and once after each accepted one. It does not replace the monitors PETSc installs from the options database, such as -ts_monitor.

Long runs

A monitor only observes. To change the run between steps, for example to update a history variable once a step is accepted, use a post-step hook; a pre-step hook runs before each step is attempted:

PETSc.set_post_step!(ts) do ts
    u = PETSc.solution(ts)       # the accepted solution, which may be changed here
    # ...
end

Both run inside PETSc.solve!, not PETSc.step!.

What stops a run early, and how:

PETSc.set_max_snes_failures!(ts, 5)       # failed nonlinear solves over the run (default 1)
PETSc.set_max_step_rejections!(ts, 20)    # rejected attempts at one step (default 10)
PETSc.set_error_if_step_fails!(ts, false) # return instead of throwing; then read converged_reason

LibPETSc.PETSC_UNLIMITED removes either limit. A domain error reported with PETSc.set_function_domain_error! rejects the step rather than counting as a failed solve.

To continue from a checkpoint, set the clock and the step count before solving:

PETSc.set_time!(ts, t_checkpoint)
PETSc.set_step_number!(ts, step_checkpoint)

For a DAE, PETSc.set_equation_type! tells the integrator the form of the equations, for example LibPETSc.TS_EQ_DAE_IMPLICIT_INDEX1.

Solving, and reading the result

PETSc.solve!(u, ts)          # integrate, starting from u
PETSc.solve!(ts)             # or from the vector already set on ts

PETSc.converged_reason(ts)   # why it stopped
PETSc.solve_time(ts)         # the time actually reached
PETSc.prev_time(ts)          # the time at the start of the last step
PETSc.step_number(ts)        # steps taken
PETSc.snes_iterations(ts)    # nonlinear iterations, summed over the steps

The final step lands exactly on PETSc.set_max_time! by default. This differs from PETSc, whose own default, TS_EXACTFINALTIME_UNSPECIFIED, integrates past the requested time without saying so; pass exact_final_time to the constructor or call PETSc.set_exact_final_time! to choose otherwise.

For finer control, PETSc.step! takes a single step and PETSc.interpolate! samples the solution inside the step just taken.

Cleaning up

PETSc.destroy!(ts)

Safe to call more than once, and a no-op on a handle left over from a previous initialize/finalize cycle.

Functions

PETSc.LibPETSc.TS — Method

TS(petsclib, comm::MPI.Comm; prefix = "", options...)

Create a PETSc time stepper on the communicator comm.

Options are stored and applied in solve! rather than here, so that a DM and callbacks attached after construction are visible to TSSetFromOptions. They are applied on every solve!, and override a setting made in code for the same option. step! does not apply them.

exact_final_time defaults to TS_EXACTFINALTIME_MATCHSTEP, so the last step lands on the time set by set_max_time!. PETSc's own default is TS_EXACTFINALTIME_UNSPECIFIED, which integrates to the wrong time without reporting an error.

If comm has size 1 the garbage collector calls destroy!. Otherwise destruction is the caller's responsibility.

External Links

PETSc Manual: TS/TSCreate