From c0f510fee8b06199909e10cfebbc80468d86ad9d Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 08:56:44 +0200
Subject: [PATCH 01/10] WIP: more TS tests and examples

---
 examples/ex51.jl     | 186 +++++++++++++++++++++++++++++++++++++++++++
 src/PETSc.jl         |   1 +
 src/ts.jl            | 106 ++++++++++++++++++++++++
 test/examples.jl     |   2 +-
 test/mpi_examples.jl |   2 +-
 5 files changed, 295 insertions(+), 2 deletions(-)
 create mode 100644 examples/ex51.jl
 create mode 100644 src/ts.jl

diff --git a/examples/ex51.jl b/examples/ex51.jl
new file mode 100644
index 00000000..f6e9adce
--- /dev/null
+++ b/examples/ex51.jl
@@ -0,0 +1,186 @@
+using PETSc, MPI, Printf
+
+mutable struct Ex51Context{PetscLib <: PETSc.LibPETSc.PetscLibType}
+    petsclib::PetscLib
+end
+
+function ex51_rhs_callback(
+    ::PETSc.LibPETSc.CTS,
+    t::Float64,
+    u_ptr::PETSc.LibPETSc.CVec,
+    f_ptr::PETSc.LibPETSc.CVec,
+    ctx_ptr::Ptr{Cvoid},
+)::PETSc.LibPETSc.PetscErrorCode
+    ctx = unsafe_pointer_to_objref(ctx_ptr)::Ex51Context
+    petsclib = ctx.petsclib
+    age = petsclib.age
+    u = PETSc.LibPETSc.PetscVec(u_ptr, petsclib, age)
+    f = PETSc.LibPETSc.PetscVec(f_ptr, petsclib, age)
+
+    PETSc.withlocalarray!(
+        (u, f);
+        read = (true, false),
+        write = (false, true),
+    ) do u_array, f_array
+        f_array[1] = cos(t)
+        f_array[2] = sin(u_array[2])
+    end
+
+    return PETSc.LibPETSc.PetscErrorCode(0)
+end
+
+const EX51_RHS_FUNCTION_PTR = @cfunction(
+    ex51_rhs_callback,
+    PETSc.LibPETSc.PetscErrorCode,
+    (
+        PETSc.LibPETSc.CTS,
+        Float64,
+        PETSc.LibPETSc.CVec,
+        PETSc.LibPETSc.CVec,
+        Ptr{Cvoid},
+    ),
+)
+
+function set_initial_condition!(u::PETSc.LibPETSc.PetscVec)
+    PETSc.withlocalarray!(u; read = false, write = true) do u_array
+        u_array[1] = 0.0
+        u_array[2] = 1.0
+    end
+    return nothing
+end
+
+function exact_solution!(u::PETSc.LibPETSc.PetscVec, t::Real)
+    PETSc.withlocalarray!(u; read = false, write = true) do u_array
+        u_array[1] = sin(t)
+        u_array[2] = 2 * atan(exp(t) * tan(0.5))
+    end
+    return nothing
+end
+
+"""
+    solve_ex51(; kwargs...)
+
+Solve the small ODE from PETSc TS tutorial `ex51.c` using PETSc.jl's low-level
+TS wrappers.
+
+Keyword arguments:
+
+- `petsclib`: PETSc library instance. Defaults to `Float64`.
+- `final_time`: final integration time. Default `1.0`.
+- `dt`: initial time step. Default `0.25`.
+- `save_trajectory`: whether to call `TSSetSaveTrajectory`. Default `true`.
+- `options`: additional PETSc command-line options. By default the script uses
+  `ARGS`, so invocations like `julia --project=. examples/ex51.jl -ts_type rk
+  -ts_rk_type 5dp` work.
+- `verbose`: print a short run summary.
+
+Returns a named tuple with the final time, error norm, and numerical solution.
+"""
+function solve_ex51(;
+    petsclib = PETSc.getlib(PetscScalar = Float64),
+    final_time::Real = 1.0,
+    dt::Real = 0.25,
+    save_trajectory::Bool = true,
+    options = nothing,
+    verbose::Bool = true,
+)
+    options === nothing && (options = copy(String.(ARGS)))
+
+    comm = MPI.COMM_WORLD
+    PetscInt = petsclib.PetscInt
+    PetscScalar = petsclib.PetscScalar
+    did_initialize = !PETSc.initialized(petsclib)
+
+    if did_initialize
+        PETSc.initialize(petsclib; options = copy(options))
+    end
+
+    MPI.Comm_size(comm) == 1 || error("This example only supports sequential runs.")
+
+    ts = PETSc.LibPETSc.TS(petsclib)
+    u = PETSc.LibPETSc.PetscVec(petsclib)
+    u_exact = PETSc.LibPETSc.PetscVec(petsclib)
+    current_time = petsclib.PetscReal(NaN)
+    error_norm = petsclib.PetscReal(NaN)
+    solution = PetscScalar[]
+    ctx = Ex51Context(petsclib)
+
+    try
+        ts = PETSc.LibPETSc.TSCreate(petsclib, comm)
+        PETSc.LibPETSc.TSSetType(petsclib, ts, "rosw")
+        PETSc.LibPETSc.TSSetProblemType(petsclib, ts, PETSc.LibPETSc.TS_NONLINEAR)
+
+        u = PETSc.LibPETSc.VecCreate(petsclib, comm)
+        PETSc.LibPETSc.VecSetSizes(petsclib, u, PetscInt(2), PetscInt(2))
+        PETSc.LibPETSc.VecSetUp(petsclib, u)
+        set_initial_condition!(u)
+        PETSc.LibPETSc.TSSetSolution(petsclib, ts, u)
+
+        callback_ctx_ptr = pointer_from_objref(ctx)
+        GC.@preserve ctx begin
+            PETSc.LibPETSc.TSSetRHSFunction(
+                petsclib,
+                ts,
+                nothing,
+                EX51_RHS_FUNCTION_PTR,
+                callback_ctx_ptr,
+            )
+
+            save_trajectory && PETSc.LibPETSc.TSSetSaveTrajectory(petsclib, ts)
+            PETSc.LibPETSc.TSSetMaxTime(petsclib, ts, petsclib.PetscReal(final_time))
+            PETSc.LibPETSc.TSSetExactFinalTime(
+                petsclib,
+                ts,
+                PETSc.LibPETSc.TS_EXACTFINALTIME_STEPOVER,
+            )
+            PETSc.LibPETSc.TSSetTimeStep(petsclib, ts, petsclib.PetscReal(dt))
+
+            adapt = PETSc.LibPETSc.TSGetAdapt(petsclib, ts)
+            PETSc.LibPETSc.TSAdaptSetType(petsclib, adapt, "none")
+
+            PETSc.LibPETSc.TSSetFromOptions(petsclib, ts)
+            PETSc.LibPETSc.TSSolve(petsclib, ts, u)
+        end
+
+        current_time = PETSc.LibPETSc.TSGetTime(petsclib, ts)
+        u_exact = PETSc.LibPETSc.VecDuplicate(petsclib, u)
+        exact_solution!(u_exact, current_time)
+
+        PETSc.LibPETSc.VecAYPX(petsclib, u_exact, PetscScalar(-1), u)
+        error_norm = PETSc.LibPETSc.VecNorm(petsclib, u_exact, PETSc.NORM_2)
+        solution = copy(u[:])
+
+        if verbose
+            if abs(current_time - final_time) > 100 * eps(petsclib.PetscReal)
+                @printf(
+                    "Note: prescribed final time %.16g differs from actual final time %.16g\n",
+                    final_time,
+                    current_time,
+                )
+            end
+
+            @printf("TS type: %s\n", PETSc.LibPETSc.TSGetType(petsclib, ts))
+            @printf("Final time: %.16g\n", current_time)
+            @printf("Error at final time: %.2E\n", error_norm)
+        end
+
+        return (final_time = current_time, error = error_norm, solution = solution)
+    finally
+        if u_exact.ptr != C_NULL
+            PETSc.destroy(u_exact)
+        end
+        if u.ptr != C_NULL
+            PETSc.destroy(u)
+        end
+        if ts.ptr != C_NULL
+            PETSc.LibPETSc.TSDestroy(petsclib, ts)
+        end
+        if did_initialize
+            PETSc.finalize(petsclib)
+        end
+    end
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    solve_ex51()
+end
diff --git a/src/PETSc.jl b/src/PETSc.jl
index 7a13b0ce..53a8d1f9 100644
--- a/src/PETSc.jl
+++ b/src/PETSc.jl
@@ -30,6 +30,7 @@ include("init.jl")
 include("vec.jl")       
 include("mat.jl")          
 include("options.jl")
+include("ts.jl")
 include("ksp.jl")
 include("snes.jl")          
 include("dm.jl")          
diff --git a/src/ts.jl b/src/ts.jl
new file mode 100644
index 00000000..16cd0019
--- /dev/null
+++ b/src/ts.jl
@@ -0,0 +1,106 @@
+"""
+    TSSetRHSFunction(petsclib, ts, r, fptr::Ptr{Cvoid}, ctx = C_NULL)
+
+Convenience overload for low-level TS RHS callbacks created with `@cfunction`.
+
+The generated bindings currently accept the PETSc function-wrapper type directly,
+while Julia's `@cfunction` returns a raw pointer. This overload bridges that
+gap so callback-based TS examples can use the low-level interface naturally.
+"""
+function LibPETSc.TSSetRHSFunction(
+    petsclib::LibPETSc.PetscLibType,
+    ts::LibPETSc.TS,
+    r::AbstractPetscVec,
+    fptr::Ptr{Cvoid},
+    ctx::Ptr{Cvoid} = C_NULL,
+) end
+
+LibPETSc.@for_petsc function LibPETSc.TSSetRHSFunction(
+    petsclib::$UnionPetscLib,
+    ts::LibPETSc.TS,
+    r::AbstractPetscVec{$PetscLib},
+    fptr::Ptr{Cvoid},
+    ctx::Ptr{Cvoid} = C_NULL,
+)
+    typed_fptr = Base.unsafe_convert(Ptr{LibPETSc.TSRHSFunctionFn}, fptr)
+    LibPETSc.@chk ccall(
+        (:TSSetRHSFunction, $petsc_library),
+        LibPETSc.PetscErrorCode,
+        (LibPETSc.CTS, LibPETSc.CVec, Ptr{LibPETSc.TSRHSFunctionFn}, Ptr{Cvoid}),
+        ts,
+        r,
+        typed_fptr,
+        ctx,
+    )
+    return nothing
+end
+
+function LibPETSc.TSSetRHSFunction(
+    petsclib::LibPETSc.PetscLibType,
+    ts::LibPETSc.TS,
+    ::Nothing,
+    fptr::Ptr{Cvoid},
+    ctx::Ptr{Cvoid} = C_NULL,
+) end
+
+LibPETSc.@for_petsc function LibPETSc.TSSetRHSFunction(
+    petsclib::$UnionPetscLib,
+    ts::LibPETSc.TS,
+    ::Nothing,
+    fptr::Ptr{Cvoid},
+    ctx::Ptr{Cvoid} = C_NULL,
+)
+    typed_fptr = Base.unsafe_convert(Ptr{LibPETSc.TSRHSFunctionFn}, fptr)
+    LibPETSc.@chk ccall(
+        (:TSSetRHSFunction, $petsc_library),
+        LibPETSc.PetscErrorCode,
+        (LibPETSc.CTS, LibPETSc.CVec, Ptr{LibPETSc.TSRHSFunctionFn}, Ptr{Cvoid}),
+        ts,
+        C_NULL,
+        typed_fptr,
+        ctx,
+    )
+    return nothing
+end
+
+"""
+    adapt = TSGetAdapt(petsclib, ts)
+
+Return the adaptive time-step controller attached to `ts`.
+"""
+function LibPETSc.TSGetAdapt(
+    petsclib::LibPETSc.PetscLibType,
+    ts::LibPETSc.TS,
+) end
+
+LibPETSc.@for_petsc function LibPETSc.TSGetAdapt(
+    petsclib::$UnionPetscLib,
+    ts::LibPETSc.TS,
+)
+    adapt_ref = Ref{LibPETSc.TSAdapt}()
+    LibPETSc.@chk ccall(
+        (:TSGetAdapt, $petsc_library),
+        LibPETSc.PetscErrorCode,
+        (LibPETSc.CTS, Ptr{LibPETSc.TSAdapt}),
+        ts,
+        adapt_ref,
+    )
+    return adapt_ref[]
+end
+
+"""
+    TSAdaptSetType(petsclib, adapt, type::String)
+
+Convenience wrapper for setting the TS adaptivity controller using a Julia
+string such as `"none"` or `"basic"`.
+"""
+function LibPETSc.TSAdaptSetType(
+    petsclib::LibPETSc.PetscLibType,
+    adapt::LibPETSc.TSAdapt,
+    type::String,
+)
+    c_str = Vector{UInt8}(type * "\0")
+    ptr = Base.unsafe_convert(Ptr{Cchar}, pointer(c_str))
+    LibPETSc.TSAdaptSetType(petsclib, adapt, ptr)
+    return nothing
+end
diff --git a/test/examples.jl b/test/examples.jl
index 87368712..666eba5f 100644
--- a/test/examples.jl
+++ b/test/examples.jl
@@ -5,7 +5,7 @@ using MPI
 using .PETScTestUtils: find_sources
 
 @testset "examples" begin
-  examples_dir = joinpath(@__DIR__, "..", "examples")
+  examples_dir = joinpath(dirname(@__DIR__), "examples")
   examples = find_sources(examples_dir)
   filter!(file -> readline(file) != "# EXCLUDE FROM TESTING", examples)
 
diff --git a/test/mpi_examples.jl b/test/mpi_examples.jl
index 18a838ce..7cc0216c 100644
--- a/test/mpi_examples.jl
+++ b/test/mpi_examples.jl
@@ -4,7 +4,7 @@ using MPI
 using .PETScTestUtils: find_sources
 
 @testset "mpi examples" begin
-  examples_dir = joinpath(@__DIR__, "..", "examples")
+  examples_dir = joinpath(dirname(@__DIR__), "examples")
   examples = find_sources(examples_dir)
   filter!(file -> readline(file) == "# INCLUDE IN MPI TEST", examples)
 

From 37c935393fa69b24934c237a33434fddbb2d3c17 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 09:04:56 +0200
Subject: [PATCH 02/10] WIP

---
 .gitignore       |  5 ++++-
 examples/ex51.jl | 29 +++++++++++++++++++++++++++--
 test/runtests.jl |  5 ++---
 3 files changed, 33 insertions(+), 6 deletions(-)

diff --git a/.gitignore b/.gitignore
index 7d095569..9bb70b68 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,4 +8,7 @@ deps/petsc-*.tar.gz
 docs/build/
 docs/site/
 Manifest.toml
-wrapping/.CondaPkg
\ No newline at end of file
+wrapping/.CondaPkg
+**/.vscode
+**/run/
+
diff --git a/examples/ex51.jl b/examples/ex51.jl
index f6e9adce..bc037578 100644
--- a/examples/ex51.jl
+++ b/examples/ex51.jl
@@ -1,10 +1,29 @@
 using PETSc, MPI, Printf
 
+# Small ODE to test TS accuracy.
+#
+# The ODE
+#                 u1_t = cos(t),
+#                 u2_t = sin(u2)
+# with analytical solution
+#                 u1(t) = sin(t),
+#                 u2(t) = 2 * atan(exp(t) * tan(0.5))
+# is used to test the accuracy of TS schemes.
+
+# PETSc callbacks take a raw `void*` context pointer. We keep the active
+# `petsclib` inside a mutable Julia object so we can pass a stable reference
+# through `pointer_from_objref()` and recover it inside the RHS routine.
 mutable struct Ex51Context{PetscLib <: PETSc.LibPETSc.PetscLibType}
     petsclib::PetscLib
 end
 
-function ex51_rhs_callback(
+"""
+    ex51_rhs!
+
+Right-hand side routine passed to `TSSetRHSFunction`, corresponding to
+`RHSFunction` in PETSc's `ex51.c`.
+"""
+function ex51_rhs!(
     ::PETSc.LibPETSc.CTS,
     t::Float64,
     u_ptr::PETSc.LibPETSc.CVec,
@@ -30,7 +49,7 @@ function ex51_rhs_callback(
 end
 
 const EX51_RHS_FUNCTION_PTR = @cfunction(
-    ex51_rhs_callback,
+    ex51_rhs!,
     PETSc.LibPETSc.PetscErrorCode,
     (
         PETSc.LibPETSc.CTS,
@@ -106,10 +125,12 @@ function solve_ex51(;
     ctx = Ex51Context(petsclib)
 
     try
+        # Create timestepping solver context.
         ts = PETSc.LibPETSc.TSCreate(petsclib, comm)
         PETSc.LibPETSc.TSSetType(petsclib, ts, "rosw")
         PETSc.LibPETSc.TSSetProblemType(petsclib, ts, PETSc.LibPETSc.TS_NONLINEAR)
 
+        # Set initial conditions.
         u = PETSc.LibPETSc.VecCreate(petsclib, comm)
         PETSc.LibPETSc.VecSetSizes(petsclib, u, PetscInt(2), PetscInt(2))
         PETSc.LibPETSc.VecSetUp(petsclib, u)
@@ -126,6 +147,8 @@ function solve_ex51(;
                 callback_ctx_ptr,
             )
 
+            # Configure solver options. As in the C example, adaptivity is
+            # forced to take constant time steps.
             save_trajectory && PETSc.LibPETSc.TSSetSaveTrajectory(petsclib, ts)
             PETSc.LibPETSc.TSSetMaxTime(petsclib, ts, petsclib.PetscReal(final_time))
             PETSc.LibPETSc.TSSetExactFinalTime(
@@ -142,6 +165,8 @@ function solve_ex51(;
             PETSc.LibPETSc.TSSolve(petsclib, ts, u)
         end
 
+        # Compute the error against the analytical solution at the achieved
+        # final time, matching the PETSc example's `VecAYPX` + `VecNorm` path.
         current_time = PETSc.LibPETSc.TSGetTime(petsclib, ts)
         u_exact = PETSc.LibPETSc.VecDuplicate(petsclib, u)
         exact_solution!(u_exact, current_time)
diff --git a/test/runtests.jl b/test/runtests.jl
index a6f54027..351a37a9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -30,8 +30,8 @@ include("dmnetwork.jl")     # new test for DMNetwork example
 include("dmshell.jl")       # new test for DMShell example
 include("dmproduct.jl")     # test for DMProduct example
 include("matshell.jl")      # autowrapped!
-include("test_dmstag.jl") 
-include("test_snes.jl")  
+include("test_dmstag.jl")
+include("test_snes.jl")
 include("old_test.jl")
 include("low_level_viewer.jl")  # Low-level viewer convenience functions
 include("low_level_ts.jl")      # Low-level TS functions
@@ -68,4 +68,3 @@ if do_mpi
     end
 end
 
-

From f0133c48303d043b9b9e32548013c613923f31d3 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 09:13:01 +0200
Subject: [PATCH 03/10] comment

---
 examples/ex51.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/examples/ex51.jl b/examples/ex51.jl
index bc037578..98693266 100644
--- a/examples/ex51.jl
+++ b/examples/ex51.jl
@@ -12,7 +12,8 @@ using PETSc, MPI, Printf
 
 # PETSc callbacks take a raw `void*` context pointer. We keep the active
 # `petsclib` inside a mutable Julia object so we can pass a stable reference
-# through `pointer_from_objref()` and recover it inside the RHS routine.
+# through `pointer_from_objref()` and recover it inside the RHS routine
+# to access the library-specific types, functions, and constants.
 mutable struct Ex51Context{PetscLib <: PETSc.LibPETSc.PetscLibType}
     petsclib::PetscLib
 end
@@ -24,7 +25,7 @@ Right-hand side routine passed to `TSSetRHSFunction`, corresponding to
 `RHSFunction` in PETSc's `ex51.c`.
 """
 function ex51_rhs!(
-    ::PETSc.LibPETSc.CTS,
+    ::PETSc.LibPETSc.CTS, # optional argument `TS ts` not needed here
     t::Float64,
     u_ptr::PETSc.LibPETSc.CVec,
     f_ptr::PETSc.LibPETSc.CVec,

From f6ec1e1031702761f3339638ddc6c3f8e98a5bf8 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 09:19:08 +0200
Subject: [PATCH 04/10] wip

---
 examples/ex51.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/examples/ex51.jl b/examples/ex51.jl
index 98693266..3bda986d 100644
--- a/examples/ex51.jl
+++ b/examples/ex51.jl
@@ -117,6 +117,10 @@ function solve_ex51(;
 
     MPI.Comm_size(comm) == 1 || error("This example only supports sequential runs.")
 
+    # Keep these variables concretely typed across the whole function, even
+    # though the actual PETSc objects are created later inside the `try` block.
+    # The null-pointer placeholders are overwritten by `TSCreate`/`VecCreate`,
+    # and the `finally` block checks `ptr != C_NULL` before destroying them.
     ts = PETSc.LibPETSc.TS(petsclib)
     u = PETSc.LibPETSc.PetscVec(petsclib)
     u_exact = PETSc.LibPETSc.PetscVec(petsclib)

From ce1ecaafb95f6548f4a03c8a9f49693a6b173482 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 09:27:38 +0200
Subject: [PATCH 05/10] wip

---
 examples/ex51.jl | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/examples/ex51.jl b/examples/ex51.jl
index 3bda986d..b19b6010 100644
--- a/examples/ex51.jl
+++ b/examples/ex51.jl
@@ -31,11 +31,17 @@ function ex51_rhs!(
     f_ptr::PETSc.LibPETSc.CVec,
     ctx_ptr::Ptr{Cvoid},
 )::PETSc.LibPETSc.PetscErrorCode
+    # This function is called through PETSc's C callback interface, so its argument
+    # types must match the low-level PETSc ABI. In particular, the state and RHS
+    # vectors arrive as raw `CVec` pointers. We immediately wrap them using
+    # PETSc.jl's higher-level `VecPtr(..., own = false)` helper so we can use
+    # Julia-friendly helpers such as `withlocalarray!` without taking ownership
+    # away from PETSc.
     ctx = unsafe_pointer_to_objref(ctx_ptr)::Ex51Context
     petsclib = ctx.petsclib
-    age = petsclib.age
-    u = PETSc.LibPETSc.PetscVec(u_ptr, petsclib, age)
-    f = PETSc.LibPETSc.PetscVec(f_ptr, petsclib, age)
+    # `own = false` since memory is managed by PETSc internally
+    u = PETSc.VecPtr(petsclib, u_ptr, false)
+    f = PETSc.VecPtr(petsclib, f_ptr, false)
 
     PETSc.withlocalarray!(
         (u, f);

From 701f78c5530a92b5a177328de777f5650e1e83d4 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 09:36:01 +0200
Subject: [PATCH 06/10] clean up

---
 examples/ex51.jl | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/examples/ex51.jl b/examples/ex51.jl
index b19b6010..6ff13217 100644
--- a/examples/ex51.jl
+++ b/examples/ex51.jl
@@ -113,7 +113,6 @@ function solve_ex51(;
     options === nothing && (options = copy(String.(ARGS)))
 
     comm = MPI.COMM_WORLD
-    PetscInt = petsclib.PetscInt
     PetscScalar = petsclib.PetscScalar
     did_initialize = !PETSc.initialized(petsclib)
 
@@ -125,7 +124,7 @@ function solve_ex51(;
 
     # Keep these variables concretely typed across the whole function, even
     # though the actual PETSc objects are created later inside the `try` block.
-    # The null-pointer placeholders are overwritten by `TSCreate`/`VecCreate`,
+    # The null-pointer placeholders are overwritten by `TSCreate`/`VecSeq`,
     # and the `finally` block checks `ptr != C_NULL` before destroying them.
     ts = PETSc.LibPETSc.TS(petsclib)
     u = PETSc.LibPETSc.PetscVec(petsclib)
@@ -142,9 +141,7 @@ function solve_ex51(;
         PETSc.LibPETSc.TSSetProblemType(petsclib, ts, PETSc.LibPETSc.TS_NONLINEAR)
 
         # Set initial conditions.
-        u = PETSc.LibPETSc.VecCreate(petsclib, comm)
-        PETSc.LibPETSc.VecSetSizes(petsclib, u, PetscInt(2), PetscInt(2))
-        PETSc.LibPETSc.VecSetUp(petsclib, u)
+        u = PETSc.VecSeq(petsclib, 2)
         set_initial_condition!(u)
         PETSc.LibPETSc.TSSetSolution(petsclib, ts, u)
 
@@ -179,7 +176,7 @@ function solve_ex51(;
         # Compute the error against the analytical solution at the achieved
         # final time, matching the PETSc example's `VecAYPX` + `VecNorm` path.
         current_time = PETSc.LibPETSc.TSGetTime(petsclib, ts)
-        u_exact = PETSc.LibPETSc.VecDuplicate(petsclib, u)
+        u_exact = similar(u)
         exact_solution!(u_exact, current_time)
 
         PETSc.LibPETSc.VecAYPX(petsclib, u_exact, PetscScalar(-1), u)

From b8ae02bb9f9692365a86c605a27790eae0541dd2 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 09:43:29 +0200
Subject: [PATCH 07/10] finalize

---
 examples/ex51.jl | 13 +++++++++----
 test/runtests.jl |  2 +-
 test/ts_ex51.jl  | 24 ++++++++++++++++++++++++
 3 files changed, 34 insertions(+), 5 deletions(-)
 create mode 100644 test/ts_ex51.jl

diff --git a/examples/ex51.jl b/examples/ex51.jl
index 6ff13217..055f2ed9 100644
--- a/examples/ex51.jl
+++ b/examples/ex51.jl
@@ -95,6 +95,9 @@ Keyword arguments:
 - `final_time`: final integration time. Default `1.0`.
 - `dt`: initial time step. Default `0.25`.
 - `save_trajectory`: whether to call `TSSetSaveTrajectory`. Default `true`.
+- `finalize_petsc`: whether to finalize PETSc at the end when this function had
+  to initialize it. Default `false` so repeated calls in one Julia session work
+  reliably.
 - `options`: additional PETSc command-line options. By default the script uses
   `ARGS`, so invocations like `julia --project=. examples/ex51.jl -ts_type rk
   -ts_rk_type 5dp` work.
@@ -107,6 +110,7 @@ function solve_ex51(;
     final_time::Real = 1.0,
     dt::Real = 0.25,
     save_trajectory::Bool = true,
+    finalize_petsc::Bool = false,
     options = nothing,
     verbose::Bool = true,
 )
@@ -114,8 +118,9 @@ function solve_ex51(;
 
     comm = MPI.COMM_WORLD
     PetscScalar = petsclib.PetscScalar
-    did_initialize = !PETSc.initialized(petsclib)
 
+    # `did_initialize`: whether we initialized the library in this call
+    did_initialize = !PETSc.initialized(petsclib)
     if did_initialize
         PETSc.initialize(petsclib; options = copy(options))
     end
@@ -208,12 +213,12 @@ function solve_ex51(;
         if ts.ptr != C_NULL
             PETSc.LibPETSc.TSDestroy(petsclib, ts)
         end
-        if did_initialize
+        if did_initialize && finalize_petsc
             PETSc.finalize(petsclib)
         end
     end
 end
 
-if abspath(PROGRAM_FILE) == @__FILE__
-    solve_ex51()
+if !isinteractive() && abspath(PROGRAM_FILE) == @__FILE__
+    solve_ex51(finalize_petsc = true)
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 351a37a9..2462a15e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -35,6 +35,7 @@ include("test_snes.jl")
 include("old_test.jl")
 include("low_level_viewer.jl")  # Low-level viewer convenience functions
 include("low_level_ts.jl")      # Low-level TS functions
+include("ts_ex51.jl")           # Regression test for repeated ex51 solves
 include("low_level_is.jl")      # Low-level IS functions
 include("low_level_petscsection.jl")  # Low-level PetscSection functions
 include("low_level_tao.jl")     # Low-level Tao functions
@@ -67,4 +68,3 @@ if do_mpi
         end
     end
 end
-
diff --git a/test/ts_ex51.jl b/test/ts_ex51.jl
new file mode 100644
index 00000000..f7e5a5f2
--- /dev/null
+++ b/test/ts_ex51.jl
@@ -0,0 +1,24 @@
+using Test
+using PETSc
+
+include(joinpath(dirname(@__DIR__), "examples", "ex51.jl"))
+
+@testset "TS ex51 example" begin
+    petsclib = PETSc.getlib(PetscScalar = Float64)
+    try
+        result1 = solve_ex51(; petsclib, options = String[], verbose = false)
+        result2 = solve_ex51(; petsclib, options = String[], verbose = false)
+
+        @test result1.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result1.solution) == 2
+        @test result1.error < 0.5
+
+        @test result2.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result2.solution) == 2
+        @test result2.error < 0.5
+    finally
+        if PETSc.initialized(petsclib) && !PETSc.finalized(petsclib)
+            PETSc.finalize(petsclib)
+        end
+    end
+end

From e658482981bed1279f7b290fbda4a0454c1a60fa Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 09:55:02 +0200
Subject: [PATCH 08/10] handle options

---
 examples/ex51.jl | 23 +++++++++++++++++++----
 test/ts_ex51.jl  | 35 +++++++++++++++++++++++++++--------
 2 files changed, 46 insertions(+), 12 deletions(-)

diff --git a/examples/ex51.jl b/examples/ex51.jl
index 055f2ed9..1d605349 100644
--- a/examples/ex51.jl
+++ b/examples/ex51.jl
@@ -98,9 +98,11 @@ Keyword arguments:
 - `finalize_petsc`: whether to finalize PETSc at the end when this function had
   to initialize it. Default `false` so repeated calls in one Julia session work
   reliably.
-- `options`: additional PETSc command-line options. By default the script uses
-  `ARGS`, so invocations like `julia --project=. examples/ex51.jl -ts_type rk
-  -ts_rk_type 5dp` work.
+- `options`: additional PETSc command-line options for this solve. By default
+  the script uses `ARGS`, so invocations like
+  `julia --project=. examples/ex51.jl -ts_type rk -ts_rk_type 5dp` work, and
+  programmatic calls can use values such as
+  `options = ["-ts_type", "rk", "-ts_rk_type", "3bs"]`.
 - `verbose`: print a short run summary.
 
 Returns a named tuple with the final time, error norm, and numerical solution.
@@ -115,6 +117,7 @@ function solve_ex51(;
     verbose::Bool = true,
 )
     options === nothing && (options = copy(String.(ARGS)))
+    parsed_options = PETSc.parse_options(options)
 
     comm = MPI.COMM_WORLD
     PetscScalar = petsclib.PetscScalar
@@ -122,7 +125,7 @@ function solve_ex51(;
     # `did_initialize`: whether we initialized the library in this call
     did_initialize = !PETSc.initialized(petsclib)
     if did_initialize
-        PETSc.initialize(petsclib; options = copy(options))
+        PETSc.initialize(petsclib)
     end
 
     MPI.Comm_size(comm) == 1 || error("This example only supports sequential runs.")
@@ -138,6 +141,8 @@ function solve_ex51(;
     error_norm = petsclib.PetscReal(NaN)
     solution = PetscScalar[]
     ctx = Ex51Context(petsclib)
+    petsc_options = PETSc.Options(petsclib; parsed_options...)
+    pushed_options = false
 
     try
         # Create timestepping solver context.
@@ -174,7 +179,11 @@ function solve_ex51(;
             adapt = PETSc.LibPETSc.TSGetAdapt(petsclib, ts)
             PETSc.LibPETSc.TSAdaptSetType(petsclib, adapt, "none")
 
+            push!(petsc_options)
+            pushed_options = true
             PETSc.LibPETSc.TSSetFromOptions(petsclib, ts)
+            pop!(petsc_options)
+            pushed_options = false
             PETSc.LibPETSc.TSSolve(petsclib, ts, u)
         end
 
@@ -204,6 +213,12 @@ function solve_ex51(;
 
         return (final_time = current_time, error = error_norm, solution = solution)
     finally
+        if pushed_options
+            pop!(petsc_options)
+        end
+        if petsc_options.ptr != C_NULL
+            PETSc.destroy(petsc_options)
+        end
         if u_exact.ptr != C_NULL
             PETSc.destroy(u_exact)
         end
diff --git a/test/ts_ex51.jl b/test/ts_ex51.jl
index f7e5a5f2..48d3d16a 100644
--- a/test/ts_ex51.jl
+++ b/test/ts_ex51.jl
@@ -6,16 +6,35 @@ include(joinpath(dirname(@__DIR__), "examples", "ex51.jl"))
 @testset "TS ex51 example" begin
     petsclib = PETSc.getlib(PetscScalar = Float64)
     try
-        result1 = solve_ex51(; petsclib, options = String[], verbose = false)
-        result2 = solve_ex51(; petsclib, options = String[], verbose = false)
+        result_default_1 = solve_ex51(; petsclib, options = String[], verbose = false)
+        result_3bs = solve_ex51(;
+            petsclib,
+            options = ["-ts_type", "rk", "-ts_rk_type", "3bs"],
+            verbose = false,
+        )
 
-        @test result1.final_time ≈ 1.0 atol = 100 * eps(Float64)
-        @test length(result1.solution) == 2
-        @test result1.error < 0.5
+        result_default_2 = solve_ex51(; petsclib, options = String[], verbose = false)
+        result_5dp = solve_ex51(;
+            petsclib,
+            options = ["-ts_type", "rk", "-ts_rk_type", "5dp"],
+            verbose = false,
+        )
 
-        @test result2.final_time ≈ 1.0 atol = 100 * eps(Float64)
-        @test length(result2.solution) == 2
-        @test result2.error < 0.5
+        @test result_default_1.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result_default_1.solution) == 2
+        @test result_default_1.error < 5.0e-4
+
+        @test result_3bs.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result_3bs.solution) == 2
+        @test result_3bs.error < 5.0e-5
+
+        @test result_default_2.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result_default_2.solution) == 2
+        @test result_default_2.error ≈ result_default_1.error rtol = 1e-12
+
+        @test result_5dp.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result_5dp.solution) == 2
+        @test result_5dp.error < 5.0e-8
     finally
         if PETSc.initialized(petsclib) && !PETSc.finalized(petsclib)
             PETSc.finalize(petsclib)

From 916f1d2dc3667de3d6a34e45458514cfc4854cbe Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 10:27:32 +0200
Subject: [PATCH 09/10] fix save_trajectory issue

---
 src/init.jl     |  2 +-
 test/ts_ex51.jl | 25 +++++++++++++++++++++++--
 2 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/src/init.jl b/src/init.jl
index dd775a25..0e80d379 100644
--- a/src/init.jl
+++ b/src/init.jl
@@ -303,4 +303,4 @@ function check_petsc_wrappers_version(petsclib=nothing)
     end
 
     return (wrappers_version = wrappers_version, installed_version = installed_version, match = match)
-end
\ No newline at end of file
+end
diff --git a/test/ts_ex51.jl b/test/ts_ex51.jl
index 48d3d16a..f9c398a6 100644
--- a/test/ts_ex51.jl
+++ b/test/ts_ex51.jl
@@ -5,18 +5,39 @@ include(joinpath(dirname(@__DIR__), "examples", "ex51.jl"))
 
 @testset "TS ex51 example" begin
     petsclib = PETSc.getlib(PetscScalar = Float64)
+    PETSc.initialize(petsclib)
     try
-        result_default_1 = solve_ex51(; petsclib, options = String[], verbose = false)
+        PETSc.finalize(petsclib)
+        PETSc.initialize(petsclib)
+
+        # `examples.jl` covers the default `save_trajectory = true` path in a
+        # fresh Julia subprocess. Here we disable trajectory saving because
+        # PETSc 3.22 does not reliably re-register the built-in `basic`
+        # trajectory type after earlier TS usage and reinitialization in the
+        # same process.
+        result_default_1 = solve_ex51(;
+            petsclib,
+            options = String[],
+            save_trajectory = false,
+            verbose = false,
+        )
         result_3bs = solve_ex51(;
             petsclib,
             options = ["-ts_type", "rk", "-ts_rk_type", "3bs"],
+            save_trajectory = false,
             verbose = false,
         )
 
-        result_default_2 = solve_ex51(; petsclib, options = String[], verbose = false)
+        result_default_2 = solve_ex51(;
+            petsclib,
+            options = String[],
+            save_trajectory = false,
+            verbose = false,
+        )
         result_5dp = solve_ex51(;
             petsclib,
             options = ["-ts_type", "rk", "-ts_rk_type", "5dp"],
+            save_trajectory = false,
             verbose = false,
         )
 

From 4ef5486dc6cd160fc042414cb8fbe3926bc31efb Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Mon, 6 Apr 2026 10:27:48 +0200
Subject: [PATCH 10/10] clean up

---
 test/ts_ex51.jl | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/test/ts_ex51.jl b/test/ts_ex51.jl
index f9c398a6..28b8ca36 100644
--- a/test/ts_ex51.jl
+++ b/test/ts_ex51.jl
@@ -7,9 +7,6 @@ include(joinpath(dirname(@__DIR__), "examples", "ex51.jl"))
     petsclib = PETSc.getlib(PetscScalar = Float64)
     PETSc.initialize(petsclib)
     try
-        PETSc.finalize(petsclib)
-        PETSc.initialize(petsclib)
-
         # `examples.jl` covers the default `save_trajectory = true` path in a
         # fresh Julia subprocess. Here we disable trajectory saving because
         # PETSc 3.22 does not reliably re-register the built-in `basic`