From c53b8bf9e77e7c667ae61f0cf779b2aeaf970249 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 18 Apr 2026 18:14:14 +0200 Subject: [PATCH 1/4] WIP: add IMEX example ex16.jl --- examples/ex16.jl | 526 +++++++++++++++++++++++++++++++++++++++++++++++ src/ts.jl | 143 +++++++++++++ test/runtests.jl | 1 + test/ts_ex16.jl | 101 +++++++++ 4 files changed, 771 insertions(+) create mode 100644 examples/ex16.jl create mode 100644 test/ts_ex16.jl diff --git a/examples/ex16.jl b/examples/ex16.jl new file mode 100644 index 00000000..a64e98c7 --- /dev/null +++ b/examples/ex16.jl @@ -0,0 +1,526 @@ +using PETSc, MPI, Printf + +# Van der Pol example adapted from PETSc TS tutorial `ex16.c`. +# +# The second-order ODE +# +# y'' - mu * ((1 - y^2) * y' - y) = 0 +# +# is rewritten as the first-order system +# +# u1_t = u2 +# u2_t = mu * ((1 - u1^2) * u2 - u1) +# +# and then split into explicit and implicit pieces for IMEX timestepping: +# +# G(u, t) = [u2, 0] +# F(u_t, u, t) = [u1_t, u2_t - mu * ((1 - u1^2) * u2 - u1)] +# +# with the special case `imex = false` folding `u2` back into the implicit +# residual, matching the upstream example. + +mutable struct Ex16Context{PetscLib <: PETSc.LibPETSc.PetscLibType} + petsclib::PetscLib + mu::Float64 + imex::Bool + next_output::Float64 +end + +const EX16_REGISTERED_AGES = IdDict{Any, Int}() + +function ex16_option_bool(parsed_options::NamedTuple, key::Symbol, default::Bool) + haskey(parsed_options, key) || return default + + value = getfield(parsed_options, key) + value === nothing && return true + value isa Bool && return value + + value_string = lowercase(String(value)) + if value_string in ("1", "true", "yes", "on") + return true + elseif value_string in ("0", "false", "no", "off") + return false + end + + throw(ArgumentError("Option -$(key) expects a boolean value, got $(repr(value)).")) +end + +function ex16_solver_options(parsed_options::NamedTuple) + filtered_pairs = Pair{Symbol, Union{String, Nothing}}[] + for (key, value) in pairs(parsed_options) + key in (:mu, :imex, :monitor) && continue + push!(filtered_pairs, key => value) + end + return (; filtered_pairs...) +end + +function ex16_register_myark2!(petsclib) + get(EX16_REGISTERED_AGES, petsclib, -1) == petsclib.age && return nothing + + PetscReal = petsclib.PetscReal + A = PetscReal[ + 0.0, + 0.0, + 0.0, + 0.41421356237309504880, + 0.0, + 0.0, + 0.75, + 0.25, + 0.0, + ] + At = PetscReal[ + 0.0, + 0.0, + 0.0, + 0.12132034355964257320, + 0.29289321881345247560, + 0.0, + 0.20710678118654752440, + 0.50000000000000000000, + 0.29289321881345247560, + ] + + PETSc.LibPETSc.TSARKIMEXRegister( + petsclib, + "myark2", + 2, + 3, + At, + nothing, + nothing, + A, + nothing, + nothing, + nothing, + nothing, + 0, + nothing, + nothing, + ) + EX16_REGISTERED_AGES[petsclib] = petsclib.age + return nothing +end + +""" + ex16_rhs! + +Right-hand side callback for the explicit part of the van der Pol split. +""" +function ex16_rhs!( + ::PETSc.LibPETSc.CTS, + ::Float64, + x_ptr::PETSc.LibPETSc.CVec, + f_ptr::PETSc.LibPETSc.CVec, + ctx_ptr::Ptr{Cvoid}, +)::PETSc.LibPETSc.PetscErrorCode + ctx = unsafe_pointer_to_objref(ctx_ptr)::Ex16Context + petsclib = ctx.petsclib + x = PETSc.VecPtr(petsclib, x_ptr, false) + f = PETSc.VecPtr(petsclib, f_ptr, false) + + PETSc.withlocalarray!( + (x, f); + read = (true, false), + write = (false, true), + ) do x_array, f_array + f_array[1] = ctx.imex ? x_array[2] : 0.0 + f_array[2] = 0.0 + end + + return PETSc.LibPETSc.PetscErrorCode(0) +end + +const EX16_RHS_FUNCTION_PTR = @cfunction( + ex16_rhs!, + PETSc.LibPETSc.PetscErrorCode, + ( + PETSc.LibPETSc.CTS, + Float64, + PETSc.LibPETSc.CVec, + PETSc.LibPETSc.CVec, + Ptr{Cvoid}, + ), +) + +""" + ex16_ifunction! + +Implicit residual callback for the van der Pol problem. +""" +function ex16_ifunction!( + ::PETSc.LibPETSc.CTS, + ::Float64, + x_ptr::PETSc.LibPETSc.CVec, + xdot_ptr::PETSc.LibPETSc.CVec, + f_ptr::PETSc.LibPETSc.CVec, + ctx_ptr::Ptr{Cvoid}, +)::PETSc.LibPETSc.PetscErrorCode + ctx = unsafe_pointer_to_objref(ctx_ptr)::Ex16Context + petsclib = ctx.petsclib + x = PETSc.VecPtr(petsclib, x_ptr, false) + xdot = PETSc.VecPtr(petsclib, xdot_ptr, false) + f = PETSc.VecPtr(petsclib, f_ptr, false) + + PETSc.withlocalarray!( + (x, xdot, f); + read = (true, true, false), + write = (false, false, true), + ) do x_array, xdot_array, f_array + f_array[1] = xdot_array[1] + (ctx.imex ? 0.0 : x_array[2]) + f_array[2] = + xdot_array[2] - ctx.mu * ((1.0 - x_array[1] * x_array[1]) * x_array[2] - x_array[1]) + end + + return PETSc.LibPETSc.PetscErrorCode(0) +end + +const EX16_IFUNCTION_PTR = @cfunction( + ex16_ifunction!, + PETSc.LibPETSc.PetscErrorCode, + ( + PETSc.LibPETSc.CTS, + Float64, + PETSc.LibPETSc.CVec, + PETSc.LibPETSc.CVec, + PETSc.LibPETSc.CVec, + Ptr{Cvoid}, + ), +) + +function ex16_fill_ijacobian!( + A::PETSc.LibPETSc.PetscMat, + B::PETSc.LibPETSc.PetscMat, + a, + x1, + x2, + ctx::Ex16Context, +) + PetscScalar = ctx.petsclib.PetscScalar + mu = ctx.mu + + row1 = PetscScalar.([a, ctx.imex ? 0.0 : 1.0]) + row2 = PetscScalar.([mu * (2.0 * x1 * x2 + 1.0), a - mu * (1.0 - x1 * x1)]) + + A[1, [1, 2]] = row1 + A[2, [1, 2]] = row2 + PETSc.assemble!(A) + + if B.ptr != A.ptr + B[1, [1, 2]] = row1 + B[2, [1, 2]] = row2 + PETSc.assemble!(B) + end + + return nothing +end + +""" + ex16_ijacobian! + +Implicit Jacobian callback corresponding to `ex16.c`. +""" +function ex16_ijacobian!( + ::PETSc.LibPETSc.CTS, + ::Float64, + x_ptr::PETSc.LibPETSc.CVec, + ::PETSc.LibPETSc.CVec, + a::Float64, + A_ptr::PETSc.LibPETSc.CMat, + B_ptr::PETSc.LibPETSc.CMat, + ctx_ptr::Ptr{Cvoid}, +)::PETSc.LibPETSc.PetscErrorCode + ctx = unsafe_pointer_to_objref(ctx_ptr)::Ex16Context + petsclib = ctx.petsclib + x = PETSc.VecPtr(petsclib, x_ptr, false) + A = PETSc.LibPETSc.PetscMat(A_ptr, petsclib) + B = PETSc.LibPETSc.PetscMat(B_ptr, petsclib) + + x1, x2 = PETSc.withlocalarray!(x; read = true, write = false) do x_array + (x_array[1], x_array[2]) + end + + ex16_fill_ijacobian!(A, B, a, x1, x2, ctx) + return PETSc.LibPETSc.PetscErrorCode(0) +end + +const EX16_IJACOBIAN_PTR = @cfunction( + ex16_ijacobian!, + PETSc.LibPETSc.PetscErrorCode, + ( + PETSc.LibPETSc.CTS, + Float64, + PETSc.LibPETSc.CVec, + PETSc.LibPETSc.CVec, + Float64, + PETSc.LibPETSc.CMat, + PETSc.LibPETSc.CMat, + Ptr{Cvoid}, + ), +) + +""" + ex16_monitor! + +Optional monitor callback that interpolates the TS solution at multiples of +`0.1`, mirroring the upstream example. +""" +function ex16_monitor!( + ts_ptr::PETSc.LibPETSc.CTS, + step::PETSc.LibPETSc.PetscInt, + t::Float64, + x_ptr::PETSc.LibPETSc.CVec, + ctx_ptr::Ptr{Cvoid}, +)::PETSc.LibPETSc.PetscErrorCode + ctx = unsafe_pointer_to_objref(ctx_ptr)::Ex16Context + petsclib = ctx.petsclib + ts = PETSc.LibPETSc.TS(ts_ptr, petsclib) + x = PETSc.VecPtr(petsclib, x_ptr, false) + dt = PETSc.LibPETSc.TSGetTimeStep(petsclib, ts) + tfinal = PETSc.LibPETSc.TSGetMaxTime(petsclib, ts) + + while ctx.next_output <= t && ctx.next_output <= tfinal + interpolated_x = PETSc.VecSeq(petsclib, length(x)) + try + PETSc.LibPETSc.TSInterpolate( + petsclib, + ts, + petsclib.PetscReal(ctx.next_output), + interpolated_x, + ) + PETSc.withlocalarray!(interpolated_x; read = true, write = false) do x_array + @printf( + "[%.1f] %d TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", + ctx.next_output, + step, + t, + dt, + x_array[1], + x_array[2], + ) + end + finally + PETSc.destroy(interpolated_x) + end + + ctx.next_output += 0.1 + end + + return PETSc.LibPETSc.PetscErrorCode(0) +end + +const EX16_MONITOR_PTR = @cfunction( + ex16_monitor!, + PETSc.LibPETSc.PetscErrorCode, + ( + PETSc.LibPETSc.CTS, + PETSc.LibPETSc.PetscInt, + Float64, + PETSc.LibPETSc.CVec, + Ptr{Cvoid}, + ), +) + +function ex16_create_jacobian_template(petsclib) + PetscScalar = petsclib.PetscScalar + jac = PETSc.MatSeqAIJ(petsclib, 2, 2, petsclib.PetscInt(2)) + jac[1, [1, 2]] = PetscScalar.([1.0, 1.0]) + jac[2, [1, 2]] = PetscScalar.([1.0, 1.0]) + PETSc.assemble!(jac) + PETSc.LibPETSc.MatZeroEntries(petsclib, jac) + return jac +end + +function ex16_initial_condition!(u::PETSc.LibPETSc.PetscVec, mu::Real) + PETSc.withlocalarray!(u; read = false, write = true) do u_array + u_array[1] = 2.0 + u_array[2] = -2.0 / 3.0 + 10.0 / (81.0 * mu) - 292.0 / (2187.0 * mu * mu) + end + return nothing +end + +""" + solve_ex16(; kwargs...) + +Solve the van der Pol IMEX example from PETSc TS tutorial `ex16.c`. + +Keyword arguments: + +- `petsclib`: PETSc library instance. Defaults to `Float64`. +- `mu`: van der Pol stiffness parameter. Default `1000.0`. +- `imex`: whether to keep the `u2` term in the explicit split. Default `true`. +- `monitor`: register the interpolation monitor from the upstream example. + Default `false`. +- `final_time`: final integration time. Default `0.5`. +- `dt`: initial time step. Default `0.01`. +- `finalize_petsc`: whether to finalize PETSc at the end when this function had + to initialize it. Default `false`. +- `options`: additional PETSc command-line options for this solve. By default + the script uses `ARGS`. +- `verbose`: print a short run summary after the solve. + +Returns a named tuple with the final time, step count, solution, and effective +problem parameters. +""" +function solve_ex16(; + petsclib = PETSc.getlib(PetscScalar = Float64), + mu::Real = 1000.0, + imex::Bool = true, + monitor::Bool = false, + final_time::Real = 0.5, + dt::Real = 0.01, + finalize_petsc::Bool = false, + options = nothing, + verbose::Bool = true, +) + options === nothing && (options = copy(String.(ARGS))) + parsed_options = PETSc.parse_options(options) + effective_mu = PETSc.typedget(parsed_options, :mu, Float64(mu)) + effective_imex = ex16_option_bool(parsed_options, :imex, imex) + effective_monitor = ex16_option_bool(parsed_options, :monitor, monitor) + solver_options = ex16_solver_options(parsed_options) + + comm = MPI.COMM_WORLD + PetscScalar = petsclib.PetscScalar + + did_initialize = !PETSc.initialized(petsclib) + if did_initialize + PETSc.initialize(petsclib) + end + + MPI.Comm_size(comm) == 1 || error("This example only supports sequential runs.") + + ts = PETSc.LibPETSc.TS(petsclib) + u = PETSc.LibPETSc.PetscVec(petsclib) + jac = PETSc.LibPETSc.PetscMat(petsclib) + current_time = petsclib.PetscReal(NaN) + steps = petsclib.PetscInt(-1) + solution = PetscScalar[] + ctx = Ex16Context(petsclib, effective_mu, effective_imex, 0.0) + petsc_options = PETSc.Options(petsclib; solver_options...) + pushed_options = false + + try + ex16_register_myark2!(petsclib) + + ts = PETSc.LibPETSc.TSCreate(petsclib, comm) + PETSc.LibPETSc.TSSetType(petsclib, ts, "beuler") + PETSc.LibPETSc.TSSetProblemType(petsclib, ts, PETSc.LibPETSc.TS_NONLINEAR) + + u = PETSc.VecSeq(petsclib, 2) + ex16_initial_condition!(u, effective_mu) + PETSc.LibPETSc.TSSetSolution(petsclib, ts, u) + + jac = ex16_create_jacobian_template(petsclib) + + callback_ctx_ptr = pointer_from_objref(ctx) + GC.@preserve ctx begin + PETSc.LibPETSc.TSSetRHSFunction( + petsclib, + ts, + nothing, + EX16_RHS_FUNCTION_PTR, + callback_ctx_ptr, + ) + PETSc.LibPETSc.TSSetIFunction( + petsclib, + ts, + nothing, + EX16_IFUNCTION_PTR, + callback_ctx_ptr, + ) + PETSc.LibPETSc.TSSetIJacobian( + petsclib, + ts, + jac, + jac, + EX16_IJACOBIAN_PTR, + callback_ctx_ptr, + ) + 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)) + + if effective_monitor + PETSc.LibPETSc.TSMonitorSet( + petsclib, + ts, + EX16_MONITOR_PTR, + callback_ctx_ptr, + ) + end + + push!(petsc_options) + pushed_options = true + PETSc.LibPETSc.TSSetFromOptions(petsclib, ts) + pop!(petsc_options) + pushed_options = false + + PETSc.LibPETSc.TSSolve(petsclib, ts, u) + end + + current_time = PETSc.LibPETSc.TSGetSolveTime(petsclib, ts) + steps = PETSc.LibPETSc.TSGetStepNumber(petsclib, ts) + 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 + + ts_type = PETSc.LibPETSc.TSGetType(petsclib, ts) + @printf("TS type: %s\n", ts_type) + if ts_type == "arkimex" + @printf("ARKIMEX type: %s\n", PETSc.LibPETSc.TSARKIMEXGetType(petsclib, ts)) + end + @printf("mu: %.16g\n", effective_mu) + @printf("IMEX split: %s\n", effective_imex ? "true" : "false") + @printf("Final time: %.16g\n", current_time) + @printf("Steps: %d\n", steps) + @printf( + "Final solution: [% .16e, % .16e]\n", + solution[1], + solution[2], + ) + end + + return ( + final_time = current_time, + steps = steps, + solution = solution, + mu = effective_mu, + imex = effective_imex, + ) + finally + if pushed_options + pop!(petsc_options) + end + if petsc_options.ptr != C_NULL + PETSc.destroy(petsc_options) + end + if jac.ptr != C_NULL + PETSc.destroy(jac) + end + if u.ptr != C_NULL + PETSc.destroy(u) + end + if ts.ptr != C_NULL + PETSc.LibPETSc.TSDestroy(petsclib, ts) + end + if did_initialize && finalize_petsc + PETSc.finalize(petsclib) + end + end +end + +if !isinteractive() && abspath(PROGRAM_FILE) == @__FILE__ + solve_ex16(finalize_petsc = true) +end diff --git a/src/ts.jl b/src/ts.jl index 8a743e32..e21d5f12 100644 --- a/src/ts.jl +++ b/src/ts.jl @@ -228,3 +228,146 @@ function LibPETSc.TSAdaptSetType( LibPETSc.TSAdaptSetType(petsclib, adapt, ptr) return nothing end + +""" + TSMonitorSet(petsclib, ts, monitor::Ptr{Cvoid}, ctx = C_NULL, mdestroy = C_NULL) + +Convenience overload for low-level TS monitor callbacks created with +`@cfunction`. +""" +function LibPETSc.TSMonitorSet( + petsclib::LibPETSc.PetscLibType, + ts::LibPETSc.TS, + monitor::Ptr{Cvoid}, + ctx::Ptr{Cvoid} = C_NULL, + mdestroy::Ptr{Cvoid} = C_NULL, +) end + +LibPETSc.@for_petsc function LibPETSc.TSMonitorSet( + petsclib::$UnionPetscLib, + ts::LibPETSc.TS, + monitor::Ptr{Cvoid}, + ctx::Ptr{Cvoid} = C_NULL, + mdestroy::Ptr{Cvoid} = C_NULL, +) + typed_destroy = Ptr{LibPETSc.PetscCtxDestroyFn}(mdestroy) + LibPETSc.@chk ccall( + (:TSMonitorSet, $petsc_library), + LibPETSc.PetscErrorCode, + (LibPETSc.CTS, LibPETSc.external, Ptr{Cvoid}, Ptr{LibPETSc.PetscCtxDestroyFn}), + ts, + monitor, + ctx, + typed_destroy, + ) + return nothing +end + +""" + TSARKIMEXRegister( + petsclib, + name::String, + order, + s, + At, + bt, + ct, + A, + b, + c, + bembedt, + bembed, + pinterp, + binterpt, + binterp, + ) + +Julia-friendly overload for registering a custom `TSARKIMEX` tableau. Optional +PETSc arrays may be passed as `nothing`, which is translated to `NULL`. +""" +function LibPETSc.TSARKIMEXRegister( + petsclib::LibPETSc.PetscLibType, + name::String, + order::Integer, + s::Integer, + At::AbstractVector, + bt::Union{Nothing, AbstractVector}, + ct::Union{Nothing, AbstractVector}, + A::AbstractVector, + b::Union{Nothing, AbstractVector}, + c::Union{Nothing, AbstractVector}, + bembedt::Union{Nothing, AbstractVector}, + bembed::Union{Nothing, AbstractVector}, + pinterp::Integer, + binterpt::Union{Nothing, AbstractVector}, + binterp::Union{Nothing, AbstractVector}, +) end + +LibPETSc.@for_petsc function LibPETSc.TSARKIMEXRegister( + petsclib::$UnionPetscLib, + name::String, + order::Integer, + s::Integer, + At::AbstractVector, + bt::Union{Nothing, AbstractVector}, + ct::Union{Nothing, AbstractVector}, + A::AbstractVector, + b::Union{Nothing, AbstractVector}, + c::Union{Nothing, AbstractVector}, + bembedt::Union{Nothing, AbstractVector}, + bembed::Union{Nothing, AbstractVector}, + pinterp::Integer, + binterpt::Union{Nothing, AbstractVector}, + binterp::Union{Nothing, AbstractVector}, +) + At_vals = $PetscReal.(At) + A_vals = $PetscReal.(A) + bt_vals = bt === nothing ? Ptr{$PetscReal}(C_NULL) : $PetscReal.(bt) + ct_vals = ct === nothing ? Ptr{$PetscReal}(C_NULL) : $PetscReal.(ct) + b_vals = b === nothing ? Ptr{$PetscReal}(C_NULL) : $PetscReal.(b) + c_vals = c === nothing ? Ptr{$PetscReal}(C_NULL) : $PetscReal.(c) + bembedt_vals = + bembedt === nothing ? Ptr{$PetscReal}(C_NULL) : $PetscReal.(bembedt) + bembed_vals = + bembed === nothing ? Ptr{$PetscReal}(C_NULL) : $PetscReal.(bembed) + binterpt_vals = + binterpt === nothing ? Ptr{$PetscReal}(C_NULL) : $PetscReal.(binterpt) + binterp_vals = + binterp === nothing ? Ptr{$PetscReal}(C_NULL) : $PetscReal.(binterp) + + LibPETSc.@chk ccall( + (:TSARKIMEXRegister, $petsc_library), + LibPETSc.PetscErrorCode, + ( + LibPETSc.TSARKIMEXType, + $PetscInt, + $PetscInt, + Ptr{$PetscReal}, + Ptr{$PetscReal}, + Ptr{$PetscReal}, + Ptr{$PetscReal}, + Ptr{$PetscReal}, + Ptr{$PetscReal}, + Ptr{$PetscReal}, + Ptr{$PetscReal}, + $PetscInt, + Ptr{$PetscReal}, + Ptr{$PetscReal}, + ), + name, + $PetscInt(order), + $PetscInt(s), + At_vals, + bt_vals, + ct_vals, + A_vals, + b_vals, + c_vals, + bembedt_vals, + bembed_vals, + $PetscInt(pinterp), + binterpt_vals, + binterp_vals, + ) + return nothing +end diff --git a/test/runtests.jl b/test/runtests.jl index 6175a626..19f69e55 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,6 +37,7 @@ 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("ts_ex51_implicit.jl") # Regression test for repeated implicit Gauss solves +include("ts_ex16.jl") # Regression test for the van der Pol IMEX example 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 diff --git a/test/ts_ex16.jl b/test/ts_ex16.jl new file mode 100644 index 00000000..8f391dca --- /dev/null +++ b/test/ts_ex16.jl @@ -0,0 +1,101 @@ +using Test +using PETSc + +include(joinpath(dirname(@__DIR__), "examples", "ex16.jl")) + +function capture_stdout_to_string(f::Function) + path, io = mktemp() + close(io) + try + open(path, "w") do out + redirect_stdout(out) do + f() + end + end + return read(path, String) + finally + isfile(path) && rm(path) + end +end + +@testset "TS ex16 example" begin + petsclib = PETSc.getlib(PetscScalar = Float64) + PETSc.initialize(petsclib) + try + result_default_1 = solve_ex16(; + petsclib, + options = String[], + verbose = false, + ) + result_default_2 = solve_ex16(; + petsclib, + options = String[], + verbose = false, + ) + result_no_imex = solve_ex16(; + petsclib, + imex = false, + options = String[], + verbose = false, + ) + result_mu_100 = solve_ex16(; + petsclib, + mu = 100.0, + options = String[], + verbose = false, + ) + result_myark2 = solve_ex16(; + petsclib, + options = [ + "-ts_type", + "arkimex", + "-ts_arkimex_type", + "myark2", + "-ts_adapt_type", + "none", + ], + verbose = false, + ) + monitor_output = capture_stdout_to_string() do + solve_ex16(; + petsclib, + monitor = true, + options = String[], + verbose = false, + ) + end + + @test result_default_1.final_time ≈ 0.5 atol = 100 * eps(Float64) + @test result_default_1.steps > 0 + @test length(result_default_1.solution) == 2 + @test all(isfinite, result_default_1.solution) + + @test result_default_2.final_time ≈ 0.5 atol = 100 * eps(Float64) + @test result_default_2.steps == result_default_1.steps + @test result_default_2.solution ≈ result_default_1.solution rtol = 1e-12 + + @test result_no_imex.final_time ≈ 0.5 atol = 100 * eps(Float64) + @test result_no_imex.steps > 0 + @test !result_no_imex.imex + @test length(result_no_imex.solution) == 2 + @test all(isfinite, result_no_imex.solution) + + @test result_mu_100.final_time ≈ 0.5 atol = 100 * eps(Float64) + @test result_mu_100.steps > 0 + @test result_mu_100.mu ≈ 100.0 + @test length(result_mu_100.solution) == 2 + + @test result_myark2.final_time ≈ 0.5 atol = 100 * eps(Float64) + @test result_myark2.steps > 0 + @test length(result_myark2.solution) == 2 + @test all(isfinite, result_myark2.solution) + + @test occursin("[0.0]", monitor_output) + @test occursin("[0.5]", monitor_output) + @test length(split(chomp(monitor_output), "\n")) >= 6 + finally + if PETSc.initialized(petsclib) && !PETSc.finalized(petsclib) + PETSc.finalize(petsclib) + end + end +end From c5df47a784014ddbcb9d3182249f3dc2853d7df6 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 19 Apr 2026 11:06:00 +0200 Subject: [PATCH 2/4] add some comments and simplify the implementaiton --- examples/ex16.jl | 148 ++++++++++++++++++++++++++++++++--------------- test/ts_ex16.jl | 42 ++++++++++---- 2 files changed, 131 insertions(+), 59 deletions(-) diff --git a/examples/ex16.jl b/examples/ex16.jl index a64e98c7..e6eb8625 100644 --- a/examples/ex16.jl +++ b/examples/ex16.jl @@ -1,6 +1,7 @@ using PETSc, MPI, Printf -# Van der Pol example adapted from PETSc TS tutorial `ex16.c`. +# Van der Pol example adapted from PETSc TS tutorial `ex16.c`, see +# https://petsc.org/main/src/ts/tutorials/ex16.c.html. # # The second-order ODE # @@ -19,32 +20,29 @@ using PETSc, MPI, Printf # with the special case `imex = false` folding `u2` back into the implicit # residual, matching the upstream example. -mutable struct Ex16Context{PetscLib <: PETSc.LibPETSc.PetscLibType} +mutable struct Ex16Context{ + PetscLib <: PETSc.LibPETSc.PetscLibType, + PetscReal <: Real, +} petsclib::PetscLib - mu::Float64 + mu::PetscReal imex::Bool - next_output::Float64 + next_output::PetscReal end +# In Julia it is common to load this file once and call `solve_ex16(...)` +# repeatedly from the same interactive session. The custom `myark2` ARKIMEX +# method therefore needs a small "register once per PETSc lifetime" guard so we +# do not try to re-register it on every solve. The upstream PETSc C example +# does not need this bookkeeping because it is a one-shot program: it starts, +# registers the method, runs once, and exits. const EX16_REGISTERED_AGES = IdDict{Any, Int}() -function ex16_option_bool(parsed_options::NamedTuple, key::Symbol, default::Bool) - haskey(parsed_options, key) || return default - - value = getfield(parsed_options, key) - value === nothing && return true - value isa Bool && return value - - value_string = lowercase(String(value)) - if value_string in ("1", "true", "yes", "on") - return true - elseif value_string in ("0", "false", "no", "off") - return false - end - - throw(ArgumentError("Option -$(key) expects a boolean value, got $(repr(value)).")) -end - +# Split parsed command-line options into: +# - example-specific options handled directly in this file (`mu`, `imex`, +# `monitor`) +# - all remaining PETSc solver options, which we still want to forward to +# `TSSetFromOptions` function ex16_solver_options(parsed_options::NamedTuple) filtered_pairs = Pair{Symbol, Union{String, Nothing}}[] for (key, value) in pairs(parsed_options) @@ -54,6 +52,39 @@ function ex16_solver_options(parsed_options::NamedTuple) return (; filtered_pairs...) end +# Read the example-specific runtime options (`mu`, `imex`, `monitor`) using +# PETSc's own option-query routines so this part stays close to the upstream C +# example. This complements `ex16_solver_options(...)`: that function filters +# these keys out before the remaining options are passed on to +# `TSSetFromOptions`. +function ex16_runtime_options( + petsclib, + parsed_options::NamedTuple; + mu::Real, + imex::Bool, + monitor::Bool, +) + PetscReal = petsclib.PetscReal + query_options = PETSc.Options(petsclib; parsed_options...) + + try + mu_value, mu_set = + PETSc.LibPETSc.PetscOptionsGetReal(petsclib, query_options, "", "-mu") + imex_value, imex_set = + PETSc.LibPETSc.PetscOptionsGetBool(petsclib, query_options, "", "-imex") + monitor_value, monitor_set = + PETSc.LibPETSc.PetscOptionsGetBool(petsclib, query_options, "", "-monitor") + + return ( + mu = Bool(mu_set) ? PetscReal(mu_value) : PetscReal(mu), + imex = Bool(imex_set) ? Bool(imex_value) : imex, + monitor = Bool(monitor_set) ? Bool(monitor_value) : monitor, + ) + finally + PETSc.destroy(query_options) + end +end + function ex16_register_myark2!(petsclib) get(EX16_REGISTERED_AGES, petsclib, -1) == petsclib.age && return nothing @@ -109,7 +140,7 @@ Right-hand side callback for the explicit part of the van der Pol split. """ function ex16_rhs!( ::PETSc.LibPETSc.CTS, - ::Float64, + ::PETSc.LibPETSc.PetscReal, x_ptr::PETSc.LibPETSc.CVec, f_ptr::PETSc.LibPETSc.CVec, ctx_ptr::Ptr{Cvoid}, @@ -136,7 +167,7 @@ const EX16_RHS_FUNCTION_PTR = @cfunction( PETSc.LibPETSc.PetscErrorCode, ( PETSc.LibPETSc.CTS, - Float64, + PETSc.LibPETSc.PetscReal, PETSc.LibPETSc.CVec, PETSc.LibPETSc.CVec, Ptr{Cvoid}, @@ -150,7 +181,7 @@ Implicit residual callback for the van der Pol problem. """ function ex16_ifunction!( ::PETSc.LibPETSc.CTS, - ::Float64, + ::PETSc.LibPETSc.PetscReal, x_ptr::PETSc.LibPETSc.CVec, xdot_ptr::PETSc.LibPETSc.CVec, f_ptr::PETSc.LibPETSc.CVec, @@ -180,7 +211,7 @@ const EX16_IFUNCTION_PTR = @cfunction( PETSc.LibPETSc.PetscErrorCode, ( PETSc.LibPETSc.CTS, - Float64, + PETSc.LibPETSc.PetscReal, PETSc.LibPETSc.CVec, PETSc.LibPETSc.CVec, PETSc.LibPETSc.CVec, @@ -188,6 +219,16 @@ const EX16_IFUNCTION_PTR = @cfunction( ), ) +# PETSc's IMEX TS interface does not ask for the Jacobian of the full ODE +# right-hand side `u_t = f(u, t)`. It asks for the shifted Jacobian +# corresponding to the implicit residual supplied via `TSSetIFunction`, with +# the time-derivative contribution folded in through the scalar `a` passed to +# `IJacobian`. In this example +# `G(u, t) = [u2, 0]` +# `F(u_t, u, t) = [u1_t + (imex ? 0 : u2), +# u2_t - mu * ((1 - u1^2) * u2 - u1)]` +# so the 2x2 entries assembled below match that residual, not the explicit +# `RHSFunction` alone. function ex16_fill_ijacobian!( A::PETSc.LibPETSc.PetscMat, B::PETSc.LibPETSc.PetscMat, @@ -199,16 +240,17 @@ function ex16_fill_ijacobian!( PetscScalar = ctx.petsclib.PetscScalar mu = ctx.mu - row1 = PetscScalar.([a, ctx.imex ? 0.0 : 1.0]) - row2 = PetscScalar.([mu * (2.0 * x1 * x2 + 1.0), a - mu * (1.0 - x1 * x1)]) - - A[1, [1, 2]] = row1 - A[2, [1, 2]] = row2 + A[1, 1] = PetscScalar(a) + A[1, 2] = PetscScalar(ctx.imex ? 0.0 : 1.0) + A[2, 1] = PetscScalar(mu * (2.0 * x1 * x2 + 1.0)) + A[2, 2] = PetscScalar(a - mu * (1.0 - x1 * x1)) PETSc.assemble!(A) if B.ptr != A.ptr - B[1, [1, 2]] = row1 - B[2, [1, 2]] = row2 + B[1, 1] = PetscScalar(a) + B[1, 2] = PetscScalar(ctx.imex ? 0.0 : 1.0) + B[2, 1] = PetscScalar(mu * (2.0 * x1 * x2 + 1.0)) + B[2, 2] = PetscScalar(a - mu * (1.0 - x1 * x1)) PETSc.assemble!(B) end @@ -222,10 +264,10 @@ Implicit Jacobian callback corresponding to `ex16.c`. """ function ex16_ijacobian!( ::PETSc.LibPETSc.CTS, - ::Float64, + ::PETSc.LibPETSc.PetscReal, x_ptr::PETSc.LibPETSc.CVec, ::PETSc.LibPETSc.CVec, - a::Float64, + a::PETSc.LibPETSc.PetscReal, A_ptr::PETSc.LibPETSc.CMat, B_ptr::PETSc.LibPETSc.CMat, ctx_ptr::Ptr{Cvoid}, @@ -249,10 +291,10 @@ const EX16_IJACOBIAN_PTR = @cfunction( PETSc.LibPETSc.PetscErrorCode, ( PETSc.LibPETSc.CTS, - Float64, + PETSc.LibPETSc.PetscReal, PETSc.LibPETSc.CVec, PETSc.LibPETSc.CVec, - Float64, + PETSc.LibPETSc.PetscReal, PETSc.LibPETSc.CMat, PETSc.LibPETSc.CMat, Ptr{Cvoid}, @@ -268,7 +310,7 @@ Optional monitor callback that interpolates the TS solution at multiples of function ex16_monitor!( ts_ptr::PETSc.LibPETSc.CTS, step::PETSc.LibPETSc.PetscInt, - t::Float64, + t::PETSc.LibPETSc.PetscReal, x_ptr::PETSc.LibPETSc.CVec, ctx_ptr::Ptr{Cvoid}, )::PETSc.LibPETSc.PetscErrorCode @@ -315,7 +357,7 @@ const EX16_MONITOR_PTR = @cfunction( ( PETSc.LibPETSc.CTS, PETSc.LibPETSc.PetscInt, - Float64, + PETSc.LibPETSc.PetscReal, PETSc.LibPETSc.CVec, Ptr{Cvoid}, ), @@ -375,9 +417,6 @@ function solve_ex16(; ) options === nothing && (options = copy(String.(ARGS))) parsed_options = PETSc.parse_options(options) - effective_mu = PETSc.typedget(parsed_options, :mu, Float64(mu)) - effective_imex = ex16_option_bool(parsed_options, :imex, imex) - effective_monitor = ex16_option_bool(parsed_options, :monitor, monitor) solver_options = ex16_solver_options(parsed_options) comm = MPI.COMM_WORLD @@ -390,13 +429,26 @@ function solve_ex16(; MPI.Comm_size(comm) == 1 || error("This example only supports sequential runs.") + runtime_options = ex16_runtime_options( + petsclib, + parsed_options; + mu, + imex, + monitor, + ) + ts = PETSc.LibPETSc.TS(petsclib) u = PETSc.LibPETSc.PetscVec(petsclib) jac = PETSc.LibPETSc.PetscMat(petsclib) current_time = petsclib.PetscReal(NaN) steps = petsclib.PetscInt(-1) solution = PetscScalar[] - ctx = Ex16Context(petsclib, effective_mu, effective_imex, 0.0) + ctx = Ex16Context( + petsclib, + runtime_options.mu, + runtime_options.imex, + petsclib.PetscReal(0), + ) petsc_options = PETSc.Options(petsclib; solver_options...) pushed_options = false @@ -408,7 +460,7 @@ function solve_ex16(; PETSc.LibPETSc.TSSetProblemType(petsclib, ts, PETSc.LibPETSc.TS_NONLINEAR) u = PETSc.VecSeq(petsclib, 2) - ex16_initial_condition!(u, effective_mu) + ex16_initial_condition!(u, runtime_options.mu) PETSc.LibPETSc.TSSetSolution(petsclib, ts, u) jac = ex16_create_jacobian_template(petsclib) @@ -445,7 +497,7 @@ function solve_ex16(; ) PETSc.LibPETSc.TSSetTimeStep(petsclib, ts, petsclib.PetscReal(dt)) - if effective_monitor + if runtime_options.monitor PETSc.LibPETSc.TSMonitorSet( petsclib, ts, @@ -481,8 +533,8 @@ function solve_ex16(; if ts_type == "arkimex" @printf("ARKIMEX type: %s\n", PETSc.LibPETSc.TSARKIMEXGetType(petsclib, ts)) end - @printf("mu: %.16g\n", effective_mu) - @printf("IMEX split: %s\n", effective_imex ? "true" : "false") + @printf("mu: %.16g\n", runtime_options.mu) + @printf("IMEX split: %s\n", runtime_options.imex ? "true" : "false") @printf("Final time: %.16g\n", current_time) @printf("Steps: %d\n", steps) @printf( @@ -496,8 +548,8 @@ function solve_ex16(; final_time = current_time, steps = steps, solution = solution, - mu = effective_mu, - imex = effective_imex, + mu = runtime_options.mu, + imex = runtime_options.imex, ) finally if pushed_options diff --git a/test/ts_ex16.jl b/test/ts_ex16.jl index 8f391dca..e95515cb 100644 --- a/test/ts_ex16.jl +++ b/test/ts_ex16.jl @@ -22,6 +22,21 @@ end petsclib = PETSc.getlib(PetscScalar = Float64) PETSc.initialize(petsclib) try + function solve_arkimex(method::String) + solve_ex16(; + petsclib, + options = [ + "-ts_type", + "arkimex", + "-ts_arkimex_type", + method, + "-ts_adapt_type", + "none", + ], + verbose = false, + ) + end + result_default_1 = solve_ex16(; petsclib, options = String[], @@ -44,17 +59,13 @@ end options = String[], verbose = false, ) - result_myark2 = solve_ex16(; - petsclib, - options = [ - "-ts_type", - "arkimex", - "-ts_arkimex_type", - "myark2", - "-ts_adapt_type", - "none", - ], - verbose = false, + result_myark2 = solve_arkimex("myark2") + extra_arkimex_results = Dict( + "ars122" => solve_arkimex("ars122"), + "ars443" => solve_arkimex("ars443"), + "3" => solve_arkimex("3"), + "4" => solve_arkimex("4"), + "5" => solve_arkimex("5"), ) monitor_output = capture_stdout_to_string() do solve_ex16(; @@ -90,6 +101,15 @@ end @test length(result_myark2.solution) == 2 @test all(isfinite, result_myark2.solution) + for (method, result) in extra_arkimex_results + @testset "ARKIMEX $method" begin + @test result.final_time ≈ 0.5 atol = 100 * eps(Float64) + @test result.steps > 0 + @test length(result.solution) == 2 + @test all(isfinite, result.solution) + end + end + @test occursin("[0.0]", monitor_output) @test occursin("[0.5]", monitor_output) @test length(split(chomp(monitor_output), "\n")) >= 6 From efa007588fabe2fe4f5006afd7dd0a22af16331c Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 19 Apr 2026 11:24:42 +0200 Subject: [PATCH 3/4] more comments --- examples/ex16.jl | 28 +++++++++------------------- examples/ex51.jl | 3 ++- examples/ex51_implicit.jl | 3 ++- src/ts.jl | 5 +++++ 4 files changed, 18 insertions(+), 21 deletions(-) diff --git a/examples/ex16.jl b/examples/ex16.jl index e6eb8625..a878fd95 100644 --- a/examples/ex16.jl +++ b/examples/ex16.jl @@ -88,28 +88,18 @@ end function ex16_register_myark2!(petsclib) get(EX16_REGISTERED_AGES, petsclib, -1) == petsclib.age && return nothing + # For the stage tables `At` and `A`, PETSc expects flat vectors in row-major + # order, matching the layout used by C arrays. PetscReal = petsclib.PetscReal A = PetscReal[ - 0.0, - 0.0, - 0.0, - 0.41421356237309504880, - 0.0, - 0.0, - 0.75, - 0.25, - 0.0, + 0.0, 0.0, 0.0, + 0.41421356237309504880, 0.0, 0.0, + 0.75, 0.25, 0.0, ] At = PetscReal[ - 0.0, - 0.0, - 0.0, - 0.12132034355964257320, - 0.29289321881345247560, - 0.0, - 0.20710678118654752440, - 0.50000000000000000000, - 0.29289321881345247560, + 0.0, 0.0, 0.0, + 0.12132034355964257320, 0.29289321881345247560, 0.0, + 0.20710678118654752440, 0.5, 0.29289321881345247560, ] PETSc.LibPETSc.TSARKIMEXRegister( @@ -322,7 +312,7 @@ function ex16_monitor!( tfinal = PETSc.LibPETSc.TSGetMaxTime(petsclib, ts) while ctx.next_output <= t && ctx.next_output <= tfinal - interpolated_x = PETSc.VecSeq(petsclib, length(x)) + interpolated_x = similar(x) try PETSc.LibPETSc.TSInterpolate( petsclib, diff --git a/examples/ex51.jl b/examples/ex51.jl index 1d605349..c908d932 100644 --- a/examples/ex51.jl +++ b/examples/ex51.jl @@ -1,6 +1,7 @@ using PETSc, MPI, Printf -# Small ODE to test TS accuracy. +# Small ODE to test TS accuracy, see +# https://petsc.org/main/src/ts/tutorials/ex51.c.html. # # The ODE # u1_t = cos(t), diff --git a/examples/ex51_implicit.jl b/examples/ex51_implicit.jl index f123ff7c..d5441f87 100644 --- a/examples/ex51_implicit.jl +++ b/examples/ex51_implicit.jl @@ -1,6 +1,7 @@ using PETSc, MPI, Printf -# Small ODE to test implicit TS accuracy with Gauss/IRK schemes. +# Small ODE to test implicit TS accuracy with Gauss/IRK schemes; adapted from +# https://petsc.org/main/src/ts/tutorials/ex51.c.html. # # The ODE # u1_t = cos(t), diff --git a/src/ts.jl b/src/ts.jl index e21d5f12..e972f021 100644 --- a/src/ts.jl +++ b/src/ts.jl @@ -284,6 +284,11 @@ end Julia-friendly overload for registering a custom `TSARKIMEX` tableau. Optional PETSc arrays may be passed as `nothing`, which is translated to `NULL`. + +For the stage tables `At` and `A`, PETSc expects flat vectors in row-major +order, matching the layout used by C arrays. If you start from a Julia matrix, +do not pass `vec(A)` directly since Julia stores matrices column-major; flatten +row-by-row instead, for example with `vec(permutedims(A))`. """ function LibPETSc.TSARKIMEXRegister( petsclib::LibPETSc.PetscLibType, From bf0ec979acc3f9b15638931d55a166c1f500ae4b Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 19 Apr 2026 11:29:55 +0200 Subject: [PATCH 4/4] fix --- examples/ex16.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/ex16.jl b/examples/ex16.jl index a878fd95..6f5fc12e 100644 --- a/examples/ex16.jl +++ b/examples/ex16.jl @@ -312,7 +312,9 @@ function ex16_monitor!( tfinal = PETSc.LibPETSc.TSGetMaxTime(petsclib, ts) while ctx.next_output <= t && ctx.next_output <= tfinal - interpolated_x = similar(x) + # `x` is a borrowed `VecPtr` from PETSc's callback interface, so create + # an owning work vector explicitly instead of calling `similar(x)`. + interpolated_x = PETSc.VecSeq(petsclib, length(x)) try PETSc.LibPETSc.TSInterpolate( petsclib,