diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3c658075c..71207cd9b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@
 
 - `local_arrays` and `restore_local_arrays!` warn: use `with_local_arrays!`. An exception between the two calls left both Vecs checked out.
 - `local_coordinate_array(dm)` warns: use `with_local_coordinate_array(f, dm)`. It now returns a copy. It used to leave the coordinate vector checked out until the garbage collector freed the array, and if the DM was destroyed first, that late hand-back crashed Julia.
+- `@simple_fn`, `@residual_fn`, `@jacobian_fn` and `@bd_fn` take the library first, `@residual_fn(petsclib, f0, Nf)`, and use its `PetscInt`, `PetscReal` and `PetscScalar` for the C signature. The library value or its type both work; it is evaluated when the macro expands, so it must be a global defined before the call. Without the library, the macros take `PetscInt`, `PetscScalar` and `PetscReal` from the calling module's globals as before, and warn. That form needed those three consts and allowed only one library per module.
 
 ### Added
 
diff --git a/docs/src/man/dmplex.md b/docs/src/man/dmplex.md
index c881ab571..5b2a00471 100644
--- a/docs/src/man/dmplex.md
+++ b/docs/src/man/dmplex.md
@@ -98,8 +98,11 @@ PETSc.set_name!(petsclib, dm, "stokes")
 ## Pointwise Callbacks
 
 PETSc's FEM assembly calls user-supplied C function pointers at each
-quadrature point.  Three macros generate the required `@cfunction` wrappers
-from plain Julia functions.
+quadrature point.  Four macros (`@simple_fn`, `@residual_fn`, `@jacobian_fn`,
+`@bd_fn`) generate the required `@cfunction` wrappers from plain Julia
+functions.  Their first argument is the library (or its type), which sets the
+C types; it is evaluated when the macro expands, so it must be a global
+defined before the macro call.
 
 ### `@simple_fn` — boundary / exact-solution functions
 
@@ -111,7 +114,7 @@ function bg_vel(t, x, u, ctx)
     u[1] = x[1] * 0.1    # Vx = exx * x
     u[2] = x[2] * (-0.1) # Vy = eyy * y
 end
-const bg_vel_ptr = PETSc.@simple_fn(bg_vel)
+const bg_vel_ptr = PETSc.@simple_fn(petsclib, bg_vel)
 ```
 
 ### `@residual_fn` — f0 / f1 residual and post-processing functions
@@ -127,7 +130,7 @@ function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     f0[1] = 0.0
     for d in 0:dim_-1; f0[1] -= u_x[d*dim_+d+1]; end
 end
-const f0_p_ptr = PETSc.@residual_fn(f0_p, 1)
+const f0_p_ptr = PETSc.@residual_fn(petsclib, f0_p, 1)
 ```
 
 ### `@jacobian_fn` — g0/g1/g2/g3 Jacobian functions
@@ -140,7 +143,7 @@ function g1_pu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g1)
     for d in 0:dim_-1; g1[d*dim_+d+1] = -1.0; end
 end
-const g1_pu_ptr = PETSc.@jacobian_fn(g1_pu, dim_*dim_)
+const g1_pu_ptr = PETSc.@jacobian_fn(petsclib, g1_pu, dim_*dim_)
 ```
 
 !!! note "Automatic differentiation"
diff --git a/examples/dmplex/ex12.jl b/examples/dmplex/ex12.jl
index cde7f0801..eee9e8e0b 100644
--- a/examples/dmplex/ex12.jl
+++ b/examples/dmplex/ex12.jl
@@ -130,8 +130,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -148,13 +146,13 @@ let _b = get(NamedTuple(pairs(opts)), :bc,    nothing); global bc    = _b === no
 function quadratic_u(t, x, u, ctx)
     u[1] = sum(xi^2 for xi in x)
 end
-const quadratic_u_ptr = PETSc.@simple_fn(quadratic_u)
+const quadratic_u_ptr = PETSc.@simple_fn(petsclib, quadratic_u)
 
 # Zero initial guess for the interior DOFs.
 function zero_u(t, x, u, ctx)
     fill!(u, 0)
 end
-const zero_u_ptr = PETSc.@simple_fn(zero_u)
+const zero_u_ptr = PETSc.@simple_fn(petsclib, zero_u)
 
 # ── COEFF_FIELD:  κ(x) = x₂ + 1 ─────────────────────────────────────────────
 #
@@ -166,14 +164,14 @@ function f0_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0)
     f0[1] = (2*dim_+2)*x[2] + 2*dim_
 end
-const f0_field_ptr = PETSc.@residual_fn(f0_field, Nf)
+const f0_field_ptr = PETSc.@residual_fn(petsclib, f0_field, Nf)
 
 function f1_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f1)
     kap = x[2] + 1.0
     @inbounds for d in 1:dim_; f1[d] = kap * u_x[d]; end
 end
-const f1_field_ptr = PETSc.@residual_fn(f1_field, dim_)
+const f1_field_ptr = PETSc.@residual_fn(petsclib, f1_field, dim_)
 
 function g3_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g3)
@@ -181,7 +179,7 @@ function g3_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     fill!(g3, 0)
     @inbounds for d in 1:dim_; g3[(d-1)*dim_ + d] = kap; end
 end
-const g3_field_ptr = PETSc.@jacobian_fn(g3_field, dim_*dim_)
+const g3_field_ptr = PETSc.@jacobian_fn(petsclib, g3_field, dim_*dim_)
 
 # Neumann flux for COEFF_FIELD:  h = κ(x)·∇u·n = (x₂+1)·2(x·n)
 # Convention: f0_bd = −h  (PETSc boundary residual sign)
@@ -192,7 +190,7 @@ function f0_bd_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     @inbounds for d in 1:dim_; xdotn += x[d]*n[d]; end
     f0[1] = -kap * 2.0 * xdotn
 end
-const f0_bd_field_ptr = PETSc.@bd_fn(f0_bd_field, Nf)
+const f0_bd_field_ptr = PETSc.@bd_fn(petsclib, f0_bd_field, Nf)
 
 # ── COEFF_NONLINEAR:  κ(u) = u + 1 ──────────────────────────────────────────
 #
@@ -214,28 +212,28 @@ function f0_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     u_val = u[Int(uOff[1])+1]
     f0[1] = (2*dim_+4)*u_val + 2*dim_
 end
-const f0_nl_ptr = PETSc.@residual_fn(f0_nl, Nf)
+const f0_nl_ptr = PETSc.@residual_fn(petsclib, f0_nl, Nf)
 
 function f1_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f1)
     kap = u[Int(uOff[1])+1] + 1.0
     @inbounds for d in 1:dim_; f1[d] = kap * u_x[d]; end
 end
-const f1_nl_ptr = PETSc.@residual_fn(f1_nl, dim_)
+const f1_nl_ptr = PETSc.@residual_fn(petsclib, f1_nl, dim_)
 
 # g0: ∂f0/∂u = 2·dim+4
 function g0_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g0)
     g0[1] = 2*dim_ + 4
 end
-const g0_nl_ptr = PETSc.@jacobian_fn(g0_nl, Nf*Nf)
+const g0_nl_ptr = PETSc.@jacobian_fn(petsclib, g0_nl, Nf*Nf)
 
 # g2: ∂f1[d]/∂u = u_x[d]   (κ′(u)=1, so ∂((u+1)·u_x[d])/∂u = u_x[d])
 function g2_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g2)
     @inbounds for d in 1:dim_; g2[d] = u_x[d]; end
 end
-const g2_nl_ptr = PETSc.@jacobian_fn(g2_nl, Nf*Nf*dim_)
+const g2_nl_ptr = PETSc.@jacobian_fn(petsclib, g2_nl, Nf*Nf*dim_)
 
 # g3: ∂f1[d]/∂(∂u/∂xd′) = (u+1)·δ_{dd′}
 function g3_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -244,7 +242,7 @@ function g3_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     fill!(g3, 0)
     @inbounds for d in 1:dim_; g3[(d-1)*dim_ + d] = kap; end
 end
-const g3_nl_ptr = PETSc.@jacobian_fn(g3_nl, dim_*dim_)
+const g3_nl_ptr = PETSc.@jacobian_fn(petsclib, g3_nl, dim_*dim_)
 
 # Neumann flux for COEFF_NONLINEAR:  h = (u+1)·2(x·n)
 # f0_bd = −h
@@ -255,7 +253,7 @@ function f0_bd_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     @inbounds for d in 1:dim_; xdotn += x[d]*n[d]; end
     f0[1] = -kap * 2.0 * xdotn
 end
-const f0_bd_nl_ptr = PETSc.@bd_fn(f0_bd_nl, Nf)
+const f0_bd_nl_ptr = PETSc.@bd_fn(petsclib, f0_bd_nl, Nf)
 
 # ── Build the mesh ────────────────────────────────────────────────────────────
 # 3-D simplex (TetGen) is broken in the current PETSc_jll; use hexahedra for 3D.
diff --git a/examples/dmplex/ex17.jl b/examples/dmplex/ex17.jl
index 946ce7508..8868242a9 100644
--- a/examples/dmplex/ex17.jl
+++ b/examples/dmplex/ex17.jl
@@ -123,8 +123,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -154,39 +152,39 @@ VALID_SOL = ("vlap_quad", "elas_quad", "vlap_trig", "elas_trig",
 function zero_u(t, x, u, ctx)
     fill!(u, 0)
 end
-const zero_u_ptr = PETSc.@simple_fn(zero_u)
+const zero_u_ptr = PETSc.@simple_fn(petsclib, zero_u)
 
 function ge_shift_u(t, x, u, ctx)
     u[1] = 0.1
     for d in 2:length(u); u[d] = 0.0; end
 end
-const ge_shift_u_ptr = PETSc.@simple_fn(ge_shift_u)
+const ge_shift_u_ptr = PETSc.@simple_fn(petsclib, ge_shift_u)
 
 function quadratic_2d_u(t, x, u, ctx)
     u[1] = x[1]^2
     u[2] = x[2]^2 - 2*x[1]*x[2]
 end
-const quadratic_2d_u_ptr = PETSc.@simple_fn(quadratic_2d_u)
+const quadratic_2d_u_ptr = PETSc.@simple_fn(petsclib, quadratic_2d_u)
 
 function quadratic_3d_u(t, x, u, ctx)
     u[1] = x[1]^2
     u[2] = x[2]^2 - 2*x[1]*x[2]
     u[3] = x[3]^2 - 2*x[2]*x[3]
 end
-const quadratic_3d_u_ptr = PETSc.@simple_fn(quadratic_3d_u)
+const quadratic_3d_u_ptr = PETSc.@simple_fn(petsclib, quadratic_3d_u)
 
 function trig_2d_u(t, x, u, ctx)
     u[1] = sin(2π * x[1])
     u[2] = sin(2π * x[2]) - 2*x[1]*x[2]
 end
-const trig_2d_u_ptr = PETSc.@simple_fn(trig_2d_u)
+const trig_2d_u_ptr = PETSc.@simple_fn(petsclib, trig_2d_u)
 
 function trig_3d_u(t, x, u, ctx)
     u[1] = sin(2π * x[1])
     u[2] = sin(2π * x[2]) - 2*x[1]*x[2]
     u[3] = sin(2π * x[3]) - 2*x[2]*x[3]
 end
-const trig_3d_u_ptr = PETSc.@simple_fn(trig_3d_u)
+const trig_3d_u_ptr = PETSc.@simple_fn(petsclib, trig_3d_u)
 
 function axial_disp_u(t, x, u, ctx)
     mu = mu_val; lambda = lambda_val; N = N_val
@@ -195,7 +193,7 @@ function axial_disp_u(t, x, u, ctx)
     u[2] = 0.25*lambda / (mu*(lambda+mu)) * N * x[2]
     for d in 3:length(u); u[d] = 0.0; end
 end
-const axial_disp_u_ptr = PETSc.@simple_fn(axial_disp_u)
+const axial_disp_u_ptr = PETSc.@simple_fn(petsclib, axial_disp_u)
 
 function uniform_strain_u(t, x, u, ctx)
     eps_xx = 0.1; eps_xy = 0.3; eps_yy = 0.25
@@ -203,7 +201,7 @@ function uniform_strain_u(t, x, u, ctx)
     u[2] = eps_xy * x[1] + eps_yy * x[2]
     for d in 3:length(u); u[d] = 0.0; end
 end
-const uniform_strain_u_ptr = PETSc.@simple_fn(uniform_strain_u)
+const uniform_strain_u_ptr = PETSc.@simple_fn(petsclib, uniform_strain_u)
 
 # ── Residual functions f0 (zeroth-order) ───────────────────────────────────────
 # f0 output size = Nc = dim (one entry per displacement component).
@@ -214,7 +212,7 @@ function f0_vlap_quadratic(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                            aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0)
     for d in 1:dim_; f0[d] += 2.0; end
 end
-const f0_vlap_quadratic_ptr = PETSc.@residual_fn(f0_vlap_quadratic, dim_)
+const f0_vlap_quadratic_ptr = PETSc.@residual_fn(petsclib, f0_vlap_quadratic, dim_)
 
 # Linear elasticity, quadratic solution: body force from −∇·σ.
 function f0_elas_quadratic(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -223,14 +221,14 @@ function f0_elas_quadratic(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     for d in 1:dim_-1; f0[d] += 2.0*mu; end
     f0[dim_] += 2.0*lambda + 4.0*mu
 end
-const f0_elas_quadratic_ptr = PETSc.@residual_fn(f0_elas_quadratic, dim_)
+const f0_elas_quadratic_ptr = PETSc.@residual_fn(petsclib, f0_elas_quadratic, dim_)
 
 # Vector Laplacian, trig solution: body force = 4π² sin(2πx_d) per component.
 function f0_vlap_trig(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                       aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0)
     for d in 1:dim_; f0[d] += -4.0*π^2 * sin(2π * x[d]); end
 end
-const f0_vlap_trig_ptr = PETSc.@residual_fn(f0_vlap_trig, dim_)
+const f0_vlap_trig_ptr = PETSc.@residual_fn(petsclib, f0_vlap_trig, dim_)
 
 # Linear elasticity, trig solution: body force from −∇·σ.
 function f0_elas_trig(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -242,7 +240,7 @@ function f0_elas_trig(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                  (d < dim_ ? 2.0*(mu + lambda) : 0.0)
     end
 end
-const f0_elas_trig_ptr = PETSc.@residual_fn(f0_elas_trig, dim_)
+const f0_elas_trig_ptr = PETSc.@residual_fn(petsclib, f0_elas_trig, dim_)
 
 # ── Residual functions f1 (first-order, diffusion flux) ────────────────────────
 # f1 output size = Nc * dim = dim * dim (gradient of Nc-component field).
@@ -254,7 +252,7 @@ function f1_vlap(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_ + d + 1] += u_x[c*dim_ + d + 1]
     end
 end
-const f1_vlap_ptr = PETSc.@residual_fn(f1_vlap, dim_*dim_)
+const f1_vlap_ptr = PETSc.@residual_fn(petsclib, f1_vlap, dim_*dim_)
 
 function f1_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                  aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f1)
@@ -264,7 +262,7 @@ function f1_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_ + c + 1] += lambda * u_x[d*dim_ + d + 1]
     end
 end
-const f1_elas_ptr = PETSc.@residual_fn(f1_elas, dim_*dim_)
+const f1_elas_ptr = PETSc.@residual_fn(petsclib, f1_elas, dim_*dim_)
 
 # ── Jacobian functions ─────────────────────────────────────────────────────────
 # g3 output size = Nc*Nc*dim*dim = dim^4 (full stiffness tensor).
@@ -277,7 +275,7 @@ function g3_vlap(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((c*dim_ + c)*dim_ + d)*dim_ + d + 1] = 1.0
     end
 end
-const g3_vlap_ptr = PETSc.@jacobian_fn(g3_vlap, dim_*dim_*dim_*dim_)
+const g3_vlap_ptr = PETSc.@jacobian_fn(petsclib, g3_vlap, dim_*dim_*dim_*dim_)
 
 function g3_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                  aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g3)
@@ -289,7 +287,7 @@ function g3_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((c*dim_ + d)*dim_ + c)*dim_ + d + 1] += lambda
     end
 end
-const g3_elas_ptr = PETSc.@jacobian_fn(g3_elas, dim_*dim_*dim_*dim_)
+const g3_elas_ptr = PETSc.@jacobian_fn(petsclib, g3_elas, dim_*dim_*dim_*dim_)
 
 # ── VTK post-processing: displacement (3-vector) and stress (3×3 tensor) ──────
 #
@@ -307,7 +305,7 @@ function copy_displacement_3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     for c in 1:dim_; out[c] = u[c]; end
     for c in dim_+1:3; out[c] = 0.0; end
 end
-const copy_displacement_3_ptr = PETSc.@residual_fn(copy_displacement_3, 3)
+const copy_displacement_3_ptr = PETSc.@residual_fn(petsclib, copy_displacement_3, 3)
 
 # Compute Cauchy stress σ_{cd} = μ(∂u_c/∂x_d + ∂u_d/∂x_c) + λ tr(ε) δ_{cd}
 # and place it in a 3×3 output (9 components, row-major).
@@ -323,7 +321,7 @@ function compute_stress_3x3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     end
     for c in 0:dim_-1; out[c*3+c+1] += lambda * div_u; end
 end
-const compute_stress_3x3_ptr = PETSc.@residual_fn(compute_stress_3x3, 9)
+const compute_stress_3x3_ptr = PETSc.@residual_fn(petsclib, compute_stress_3x3, 9)
 
 # ── Neumann BC for elas_axial_disp: traction N on right wall ──────────────────
 # f0_bd output size = Nc = dim.
@@ -332,7 +330,7 @@ function f0_elas_axial_disp_bd(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                                 numConstants, constants, f0)
     f0[1] = real(constants[3])   # N — tension force, constants[3] = N (1-indexed)
 end
-const f0_elas_axial_disp_bd_ptr = PETSc.@bd_fn(f0_elas_axial_disp_bd, dim_)
+const f0_elas_axial_disp_bd_ptr = PETSc.@bd_fn(petsclib, f0_elas_axial_disp_bd, dim_)
 
 # ── Near-null space constructor for AMG ────────────────────────────────────────
 # Called by PETSc for each level when setting up the preconditioner.
diff --git a/examples/dmplex/ex62.jl b/examples/dmplex/ex62.jl
index bb4c5ef64..ebd658f87 100644
--- a/examples/dmplex/ex62.jl
+++ b/examples/dmplex/ex62.jl
@@ -104,8 +104,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -135,7 +133,7 @@ function f1_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_+c+1] -= u[dim_+1]   # subtract pressure (field 1)
     end
 end
-const f1_u_ptr = PETSc.@residual_fn(f1_u, dim_*dim_)
+const f1_u_ptr = PETSc.@residual_fn(petsclib, f1_u, dim_*dim_)
 
 # Continuity: f0 = −∇·u
 function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -145,7 +143,7 @@ function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[1] -= u_x[d*dim_+d+1]
     end
 end
-const f0_p_ptr = PETSc.@residual_fn(f0_p, 1)
+const f0_p_ptr = PETSc.@residual_fn(petsclib, f0_p, 1)
 
 # ── Body-force functions (one per MMS solution type) ──────────────────────────
 #
@@ -163,7 +161,7 @@ function f0_quadratic_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[d] = 4mu - 1.0
     end
 end
-const f0_quadratic_u_ptr = PETSc.@residual_fn(f0_quadratic_u, dim_)
+const f0_quadratic_u_ptr = PETSc.@residual_fn(petsclib, f0_quadratic_u, dim_)
 
 # Trigonometric MMS body force:
 #   f_x = 2π cos(2πx) + μ(dim−1)π² sin(πx) + μπ² Σ_{d>0} sin(πx_d)
@@ -178,7 +176,7 @@ function f0_trig_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[d]  = -2π*cos(2π*x[d]) + mu*π^3*cos(π*x[1])*x[d]
     end
 end
-const f0_trig_u_ptr = PETSc.@residual_fn(f0_trig_u, dim_)
+const f0_trig_u_ptr = PETSc.@residual_fn(petsclib, f0_trig_u, dim_)
 
 # ── Jacobians ─────────────────────────────────────────────────────────────────
 
@@ -191,28 +189,28 @@ function g3_uu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((c*dim_+d)*dim_+d)*dim_+c+1] += mu
     end
 end
-const g3_uu_ptr = PETSc.@jacobian_fn(g3_uu, dim_*dim_*dim_*dim_)
+const g3_uu_ptr = PETSc.@jacobian_fn(petsclib, g3_uu, dim_*dim_*dim_*dim_)
 
 # J_up (g2): −⟨∇·v, p⟩
 function g2_up(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g2)
     for d in 0:dim_-1; g2[d*dim_+d+1] = -1.0; end
 end
-const g2_up_ptr = PETSc.@jacobian_fn(g2_up, dim_*dim_)
+const g2_up_ptr = PETSc.@jacobian_fn(petsclib, g2_up, dim_*dim_)
 
 # J_pu (g1): ⟨q, −∇·u⟩
 function g1_pu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g1)
     for d in 0:dim_-1; g1[d*dim_+d+1] = -1.0; end
 end
-const g1_pu_ptr = PETSc.@jacobian_fn(g1_pu, dim_*dim_)
+const g1_pu_ptr = PETSc.@jacobian_fn(petsclib, g1_pu, dim_*dim_)
 
 # J_pp preconditioner (g0): (1/μ)⟨q, p⟩
 function g0_pp(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g0)
     g0[1] = 1.0 / real(cst[1])
 end
-const g0_pp_ptr = PETSc.@jacobian_fn(g0_pp, 1)
+const g0_pp_ptr = PETSc.@jacobian_fn(petsclib, g0_pp, 1)
 
 # ── Exact solutions ───────────────────────────────────────────────────────────
 #
@@ -228,12 +226,12 @@ function quadratic_vel(t, x, u, ctx)
         u[c]  = 2*x[1]^2 - 2*x[1]*x[c]
     end
 end
-const quadratic_vel_ptr = PETSc.@simple_fn(quadratic_vel)
+const quadratic_vel_ptr = PETSc.@simple_fn(petsclib, quadratic_vel)
 
 function quadratic_pres(t, x, u, ctx)
     u[1] = sum(x) - 0.5*length(x)
 end
-const quadratic_pres_ptr = PETSc.@simple_fn(quadratic_pres)
+const quadratic_pres_ptr = PETSc.@simple_fn(petsclib, quadratic_pres)
 
 #
 # Trigonometric MMS:
@@ -249,23 +247,23 @@ function trig_vel(t, x, u, ctx)
         u[c]  = -π*cos(π*x[1])*x[c]
     end
 end
-const trig_vel_ptr = PETSc.@simple_fn(trig_vel)
+const trig_vel_ptr = PETSc.@simple_fn(petsclib, trig_vel)
 
 function trig_pres(t, x, u, ctx)
     u[1] = sum(sin(2π*xi) for xi in x)
 end
-const trig_pres_ptr = PETSc.@simple_fn(trig_pres)
+const trig_pres_ptr = PETSc.@simple_fn(petsclib, trig_pres)
 
 # Helper functions for the pressure null space constructor
 function zero_vel(t, x, u, ctx)
     for c in 1:length(u); u[c] = 0.0; end
 end
-const zero_vel_ptr = PETSc.@simple_fn(zero_vel)
+const zero_vel_ptr = PETSc.@simple_fn(petsclib, zero_vel)
 
 function one_pres(t, x, u, ctx)
     u[1] = 1.0
 end
-const one_pres_ptr = PETSc.@simple_fn(one_pres)
+const one_pres_ptr = PETSc.@simple_fn(petsclib, one_pres)
 
 # ── Select exact solution functions based on sol_type ─────────────────────────
 exact_vel_ptr  = sol_type == "quadratic" ? quadratic_vel_ptr  : trig_vel_ptr
diff --git a/examples/dmplex/ex62b.jl b/examples/dmplex/ex62b.jl
index 1a4a8f199..28007dd3f 100644
--- a/examples/dmplex/ex62b.jl
+++ b/examples/dmplex/ex62b.jl
@@ -350,7 +350,6 @@ PETSc.initialize(petsclib)
 
 @isdefined(PetscInt)    || const PetscInt    = petsclib.PetscInt
 @isdefined(PetscScalar) || const PetscScalar = Float64
-@isdefined(PetscReal)   || const PetscReal   = Float64
 
 # Module-level globals for background strain rates — read by @cfunction-based BCs.
 bg_exx = 0.0
@@ -518,7 +517,7 @@ function f1_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     compute_strainrate!(f1, dim_, u_x)
     deviatoric_stress_pressure!(f1, dim_, f1, phase, cst, p, tau_old)
 end
-const f1_u_ptr = PETSc.@residual_fn(f1_u, dim_*dim_)
+const f1_u_ptr = PETSc.@residual_fn(petsclib, f1_u, dim_*dim_)
 
 # Continuity: f0 = −∇·u
 function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -528,7 +527,7 @@ function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[1] -= u_x[d*dim_+d+1]
     end
 end
-const f0_p_ptr = PETSc.@residual_fn(f0_p, 1)
+const f0_p_ptr = PETSc.@residual_fn(petsclib, f0_p, 1)
 
 # ── Body-force functions (one per MMS solution type) ──────────────────────────
 #
@@ -546,7 +545,7 @@ function f0_quadratic_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[d] = 4mu - 1.0
     end
 end
-const f0_quadratic_u_ptr = PETSc.@residual_fn(f0_quadratic_u, dim_)
+const f0_quadratic_u_ptr = PETSc.@residual_fn(petsclib, f0_quadratic_u, dim_)
 
 # Trigonometric MMS body force:
 #   f_x = 2π cos(2πx) + μ(dim−1)π² sin(πx) + μπ² Σ_{d>0} sin(πx_d)
@@ -561,14 +560,14 @@ function f0_trig_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[d]  = -2π*cos(2π*x[d]) + mu*π^3*cos(π*x[1])*x[d]
     end
 end
-const f0_trig_u_ptr = PETSc.@residual_fn(f0_trig_u, dim_)
+const f0_trig_u_ptr = PETSc.@residual_fn(petsclib, f0_trig_u, dim_)
 
 # Zero body force for background deformation mode (no MMS forcing)
 function f0_zero_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                    aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, f0)
     for d in 1:dim_; f0[d] = 0.0; end
 end
-const f0_zero_u_ptr = PETSc.@residual_fn(f0_zero_u, dim_)
+const f0_zero_u_ptr = PETSc.@residual_fn(petsclib, f0_zero_u, dim_)
 
 # Gravity body force: f_body = ρ g (−ŷ direction).
 # In PETSc weak-form sign convention, f0 = −f_body, so:
@@ -581,7 +580,7 @@ function f0_grav_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     for d in 1:dim_; f0[d] = 0.0; end
     f0[dim_] = rho * g
 end
-const f0_grav_u_ptr = PETSc.@residual_fn(f0_grav_u, dim_)
+const f0_grav_u_ptr = PETSc.@residual_fn(petsclib, f0_grav_u, dim_)
 
 # ── Jacobians ─────────────────────────────────────────────────────────────────
 
@@ -593,7 +592,7 @@ function g3_uu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     p         = real(u[dim_+1])
     deviatoric_stress_tangent!(g3, dim_, u_x, phase, cst, p, tau_old)
 end
-const g3_uu_ptr = PETSc.@jacobian_fn(g3_uu, dim_*dim_*dim_*dim_)
+const g3_uu_ptr = PETSc.@jacobian_fn(petsclib, g3_uu, dim_*dim_*dim_*dim_)
 
 # J_up (g2): ∂σ_{fc,df}/∂p via AD of deviatoric_stress_pressure! w.r.t. p.
 # Layout: g2[fc*dim+df + 1] = ∂σ_{fc,df}/∂p.
@@ -615,14 +614,14 @@ function g2_up(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     end
     for i in 1:n; g2[i] = dσ_dp[i]; end
 end
-const g2_up_ptr = PETSc.@jacobian_fn(g2_up, dim_*dim_)
+const g2_up_ptr = PETSc.@jacobian_fn(petsclib, g2_up, dim_*dim_)
 
 # J_pu (g1): ⟨q, −∇·u⟩ — kinematic incompressibility, independent of rheology.
 function g1_pu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g1)
     for d in 0:dim_-1; g1[d*dim_+d+1] = -1.0; end
 end
-const g1_pu_ptr = PETSc.@jacobian_fn(g1_pu, dim_*dim_)
+const g1_pu_ptr = PETSc.@jacobian_fn(petsclib, g1_pu, dim_*dim_)
 
 # J_pp preconditioner (g0): (1/η_eff)⟨q, p⟩ using the Maxwell effective viscosity
 function g0_pp(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -640,7 +639,7 @@ function g0_pp(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     eta_eff = eta * G_mod * dt / (eta + G_mod * dt)
     g0[1]   = 1.0 / eta_eff
 end
-const g0_pp_ptr = PETSc.@jacobian_fn(g0_pp, 1)
+const g0_pp_ptr = PETSc.@jacobian_fn(petsclib, g0_pp, 1)
 
 # ── Deviatoric stress computation (for VTK and viscoelastic history) ─────────
 #
@@ -668,7 +667,7 @@ function compute_tau(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         out[i] = 2 * eta_eff * eps_buf[i] + relax_fac * real(tau_old[i])
     end
 end
-const compute_tau_ptr = PETSc.@residual_fn(compute_tau, dim_*dim_)
+const compute_tau_ptr = PETSc.@residual_fn(petsclib, compute_tau, dim_*dim_)
 
 # Same as compute_tau but zero-padded to 3×3 row-major for VTK tensor output.
 function compute_tau_3x3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -691,7 +690,7 @@ function compute_tau_3x3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         out[c*3+d+1] = 2 * eta_eff * eps_buf[c*dim_+d+1] + relax_fac * real(tau_old[c*dim_+d+1])
     end
 end
-const compute_tau_3x3_ptr = PETSc.@residual_fn(compute_tau_3x3, 9)
+const compute_tau_3x3_ptr = PETSc.@residual_fn(petsclib, compute_tau_3x3, 9)
 
 # Second invariants (DG-P0, 1 DOF/cell) — used for domain averages and VTK output.
 # ε_II = sqrt(½ ε:ε),  τ_II = sqrt(½ τ:τ)
@@ -703,7 +702,7 @@ function compute_eps_II(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     compute_strainrate!(eps_buf, dim_, eps_buf)
     out[1] = real(second_invariant(eps_buf, dim_))
 end
-const compute_eps_II_ptr = PETSc.@residual_fn(compute_eps_II, 1)
+const compute_eps_II_ptr = PETSc.@residual_fn(petsclib, compute_eps_II, 1)
 
 function compute_tau_II(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                         aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out)
@@ -727,7 +726,7 @@ function compute_tau_II(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     end
     out[1] = sqrt(0.5 * s)
 end
-const compute_tau_II_ptr = PETSc.@residual_fn(compute_tau_II, 1)
+const compute_tau_II_ptr = PETSc.@residual_fn(petsclib, compute_tau_II, 1)
 
 # ── VTK post-processing: strain rate tensor ─────────────────────────────────
 #
@@ -741,7 +740,7 @@ function compute_strainrate_3x3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         out[c*3+d+1] = 0.5*(u_x[c*dim_+d+1] + u_x[d*dim_+c+1])
     end
 end
-const compute_strainrate_3x3_ptr = PETSc.@residual_fn(compute_strainrate_3x3, 9)
+const compute_strainrate_3x3_ptr = PETSc.@residual_fn(petsclib, compute_strainrate_3x3, 9)
 
 # VTK copy kernels: extract velocity (field 0) and pressure (field 1) from the
 # combined solution vector.  Must be at top-level scope for @cfunction.
@@ -749,13 +748,13 @@ function copy_vel(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out_)
     for c in 1:dim_; out_[c] = u[c]; end
 end
-const copy_vel_ptr = PETSc.@residual_fn(copy_vel, dim_)
+const copy_vel_ptr = PETSc.@residual_fn(petsclib, copy_vel, dim_)
 
 function copy_pres(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                    aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out_)
     out_[1] = u[dim_+1]
 end
-const copy_pres_ptr = PETSc.@residual_fn(copy_pres, 1)
+const copy_pres_ptr = PETSc.@residual_fn(petsclib, copy_pres, 1)
 
 # copy_mu reports the Maxwell effective viscosity η_eff = η·G·dt/(η+G·dt).
 function copy_mu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -772,13 +771,13 @@ function copy_mu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     dt      = real(cst[17])
     out_[1] = eta * G_mod * dt / (eta + G_mod * dt)
 end
-const copy_mu_ptr = PETSc.@residual_fn(copy_mu, 1)
+const copy_mu_ptr = PETSc.@residual_fn(petsclib, copy_mu, 1)
 
 function copy_phase(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                     aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out_)
     out_[1] = real(a[1])   # phase from aux field — correct after Lagrangian advection
 end
-const copy_phase_ptr = PETSc.@residual_fn(copy_phase, 1)
+const copy_phase_ptr = PETSc.@residual_fn(petsclib, copy_phase, 1)
 
 # Pass-through kernel: copies dim×dim values from the source field to the output.
 # Used by advance_tau! to project tau_vec (on dm_tau) into field 1 of aux_vec.
@@ -786,7 +785,7 @@ function identity_tau(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                       aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, out)
     for i in 1:dim_*dim_; out[i] = u[i]; end
 end
-const identity_tau_ptr = PETSc.@residual_fn(identity_tau, dim_*dim_)
+const identity_tau_ptr = PETSc.@residual_fn(petsclib, identity_tau, dim_*dim_)
 
 # ── Time-stepping helpers ─────────────────────────────────────────────────────
 
@@ -986,12 +985,12 @@ function quadratic_vel(t, x, u, ctx)
         u[c]  = 2*x[1]^2 - 2*x[1]*x[c]
     end
 end
-const quadratic_vel_ptr = PETSc.@simple_fn(quadratic_vel)
+const quadratic_vel_ptr = PETSc.@simple_fn(petsclib, quadratic_vel)
 
 function quadratic_pres(t, x, u, ctx)
     u[1] = sum(x) - 0.5*length(x)
 end
-const quadratic_pres_ptr = PETSc.@simple_fn(quadratic_pres)
+const quadratic_pres_ptr = PETSc.@simple_fn(petsclib, quadratic_pres)
 
 #
 # Trigonometric MMS:
@@ -1007,23 +1006,23 @@ function trig_vel(t, x, u, ctx)
         u[c]  = -π*cos(π*x[1])*x[c]
     end
 end
-const trig_vel_ptr = PETSc.@simple_fn(trig_vel)
+const trig_vel_ptr = PETSc.@simple_fn(petsclib, trig_vel)
 
 function trig_pres(t, x, u, ctx)
     u[1] = sum(sin(2π*xi) for xi in x)
 end
-const trig_pres_ptr = PETSc.@simple_fn(trig_pres)
+const trig_pres_ptr = PETSc.@simple_fn(petsclib, trig_pres)
 
 # Helper functions for the pressure null space constructor
 function zero_vel(t, x, u, ctx)
     for c in 1:length(u); u[c] = 0.0; end
 end
-const zero_vel_ptr = PETSc.@simple_fn(zero_vel)
+const zero_vel_ptr = PETSc.@simple_fn(petsclib, zero_vel)
 
 function one_pres(t, x, u, ctx)
     u[1] = 1.0
 end
-const one_pres_ptr = PETSc.@simple_fn(one_pres)
+const one_pres_ptr = PETSc.@simple_fn(petsclib, one_pres)
 
 # Background deformation BC: Vx = x*exx_bg, Vy = y*eyy_bg, Vz = z*ezz_bg.
 # ezz_bg defaults to -(exx_bg+eyy_bg) to satisfy incompressibility.
@@ -1032,7 +1031,7 @@ function bg_vel(t, x, u, ctx)
     u[2] = x[2] * bg_eyy
     length(u) >= 3 && (u[3] = x[3] * bg_ezz)
 end
-const bg_vel_ptr = PETSc.@simple_fn(bg_vel)
+const bg_vel_ptr = PETSc.@simple_fn(petsclib, bg_vel)
 
 """
     update_tau!(petsclib, dm_tau, u, tau_vec, t=0.0)
@@ -1055,7 +1054,7 @@ end
 function phase_fn(t, x, u, ctx)
     u[1] = sum(xi^2 for xi in x) < inclusion_R^2 ? 2.0 : 1.0
 end
-const phase_fn_ptr = PETSc.@simple_fn(phase_fn)
+const phase_fn_ptr = PETSc.@simple_fn(petsclib, phase_fn)
 
 exact_vel_ptr  = sol_type == "quadratic" ? quadratic_vel_ptr  :
                  sol_type == "trig"      ? trig_vel_ptr        :
diff --git a/examples/dmplex/ex69.jl b/examples/dmplex/ex69.jl
index fa4168c73..725e1f29e 100644
--- a/examples/dmplex/ex69.jl
+++ b/examples/dmplex/ex69.jl
@@ -144,8 +144,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -176,7 +174,7 @@ function f0_stokes_u(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     f0[1] = 0.0
     f0[2] = -sin(real(cst[2]) * pi * x[2]) * cos(real(cst[1]) * pi * x[1])
 end
-const f0_stokes_u_ptr = PETSc.@residual_fn(f0_stokes_u, dim_)
+const f0_stokes_u_ptr = PETSc.@residual_fn(petsclib, f0_stokes_u, dim_)
 
 # Momentum flux σ = 2μ ε(u) − pI  for SolKx (μ = exp(2Bx), cst[3]=B)
 function stokes_momentum_kx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -189,7 +187,7 @@ function stokes_momentum_kx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_+c+1] -= u[dim_+1]      # subtract pressure
     end
 end
-const stokes_momentum_kx_ptr = PETSc.@residual_fn(stokes_momentum_kx, dim_*dim_)
+const stokes_momentum_kx_ptr = PETSc.@residual_fn(petsclib, stokes_momentum_kx, dim_*dim_)
 
 # Momentum flux for SolCx (μ = etaA or etaB, cst[3]=etaA, cst[4]=etaB, cst[5]=xc)
 function stokes_momentum_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -202,7 +200,7 @@ function stokes_momentum_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_+c+1] -= u[dim_+1]
     end
 end
-const stokes_momentum_cx_ptr = PETSc.@residual_fn(stokes_momentum_cx, dim_*dim_)
+const stokes_momentum_cx_ptr = PETSc.@residual_fn(petsclib, stokes_momentum_cx, dim_*dim_)
 
 # Continuity: f0 = −∇·u
 function stokes_mass(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -212,14 +210,14 @@ function stokes_mass(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f0[1] -= u_x[d*dim_+d+1]
     end
 end
-const stokes_mass_ptr = PETSc.@residual_fn(stokes_mass, 1)
+const stokes_mass_ptr = PETSc.@residual_fn(petsclib, stokes_mass, 1)
 
 # Zero f1 for pressure field
 function f1_zero_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                    aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, f1)
     for d in 1:dim_; f1[d] = 0.0; end
 end
-const f1_zero_p_ptr = PETSc.@residual_fn(f1_zero_p, dim_)
+const f1_zero_p_ptr = PETSc.@residual_fn(petsclib, f1_zero_p, dim_)
 
 # ── Pointwise functions: Jacobians ────────────────────────────────────────────
 
@@ -232,7 +230,7 @@ function stokes_vel_J_kx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((cI*dim_+d)*dim_+d)*dim_+cI+1] += mu
     end
 end
-const stokes_vel_J_kx_ptr = PETSc.@jacobian_fn(stokes_vel_J_kx, dim_*dim_*dim_*dim_)
+const stokes_vel_J_kx_ptr = PETSc.@jacobian_fn(petsclib, stokes_vel_J_kx, dim_*dim_*dim_*dim_)
 
 function stokes_vel_J_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                           aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g3)
@@ -242,35 +240,35 @@ function stokes_vel_J_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((cI*dim_+d)*dim_+d)*dim_+cI+1] += mu
     end
 end
-const stokes_vel_J_cx_ptr = PETSc.@jacobian_fn(stokes_vel_J_cx, dim_*dim_*dim_*dim_)
+const stokes_vel_J_cx_ptr = PETSc.@jacobian_fn(petsclib, stokes_vel_J_cx, dim_*dim_*dim_*dim_)
 
 # J_up (g2): −⟨∇·v, p⟩  →  g2[d*dim+d] = -1
 function stokes_pres_J(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                         aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g2)
     for d in 0:dim_-1; g2[d*dim_+d+1] = -1.0; end
 end
-const stokes_pres_J_ptr = PETSc.@jacobian_fn(stokes_pres_J, dim_*dim_)
+const stokes_pres_J_ptr = PETSc.@jacobian_fn(petsclib, stokes_pres_J, dim_*dim_)
 
 # J_pu (g1): ⟨q, ∇·u⟩  →  g1[d*dim+d] = -1
 function stokes_mass_J(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                         aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g1)
     for d in 0:dim_-1; g1[d*dim_+d+1] = -1.0; end
 end
-const stokes_mass_J_ptr = PETSc.@jacobian_fn(stokes_mass_J, dim_*dim_)
+const stokes_mass_J_ptr = PETSc.@jacobian_fn(petsclib, stokes_mass_J, dim_*dim_)
 
 # J_pp preconditioner (g0): (1/μ)⟨q, p⟩
 function stokes_id_J_kx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                           aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g0)
     g0[1] = 1.0 / exp(2.0 * real(cst[3]) * x[1])
 end
-const stokes_id_J_kx_ptr = PETSc.@jacobian_fn(stokes_id_J_kx, 1)
+const stokes_id_J_kx_ptr = PETSc.@jacobian_fn(petsclib, stokes_id_J_kx, 1)
 
 function stokes_id_J_cx(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                           aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g0)
     mu = x[1] < real(cst[5]) ? real(cst[3]) : real(cst[4])
     g0[1] = 1.0 / mu
 end
-const stokes_id_J_cx_ptr = PETSc.@jacobian_fn(stokes_id_J_cx, 1)
+const stokes_id_J_cx_ptr = PETSc.@jacobian_fn(petsclib, stokes_id_J_cx, 1)
 
 # ── Exact solution / BC callbacks ─────────────────────────────────────────────
 # PetscSimplePointFn: f(t, x, u, ctx) → nothing
@@ -278,38 +276,38 @@ const stokes_id_J_cx_ptr = PETSc.@jacobian_fn(stokes_id_J_cx, 1)
 function zero_vel(t, x, u, ctx)
     for c in 1:length(u); u[c] = 0.0; end
 end
-const zero_vel_ptr = PETSc.@simple_fn(zero_vel)
+const zero_vel_ptr = PETSc.@simple_fn(petsclib, zero_vel)
 
 function one_pres(t, x, u, ctx)
     u[1] = 1.0
 end
-const one_pres_ptr = PETSc.@simple_fn(one_pres)
+const one_pres_ptr = PETSc.@simple_fn(petsclib, one_pres)
 
 # SolKx exact solution (captures parameters via module-level globals)
 function exact_vel_kx(t, x, u, ctx)
     vx, vz, _ = SolKxSolution(x, m_val, n_val, B_val)
     u[1] = vx; u[2] = vz
 end
-const exact_vel_kx_ptr = PETSc.@simple_fn(exact_vel_kx)
+const exact_vel_kx_ptr = PETSc.@simple_fn(petsclib, exact_vel_kx)
 
 function exact_pres_kx(t, x, u, ctx)
     _, _, p = SolKxSolution(x, m_val, n_val, B_val)
     u[1] = p
 end
-const exact_pres_kx_ptr = PETSc.@simple_fn(exact_pres_kx)
+const exact_pres_kx_ptr = PETSc.@simple_fn(petsclib, exact_pres_kx)
 
 # SolCx exact solution
 function exact_vel_cx(t, x, u, ctx)
     vx, vz, _ = SolCxSolution(x, m_val, n_val, xc_val, etaA_val, etaB_val)
     u[1] = vx; u[2] = vz
 end
-const exact_vel_cx_ptr = PETSc.@simple_fn(exact_vel_cx)
+const exact_vel_cx_ptr = PETSc.@simple_fn(petsclib, exact_vel_cx)
 
 function exact_pres_cx(t, x, u, ctx)
     _, _, p = SolCxSolution(x, m_val, n_val, xc_val, etaA_val, etaB_val)
     u[1] = p
 end
-const exact_pres_cx_ptr = PETSc.@simple_fn(exact_pres_cx)
+const exact_pres_cx_ptr = PETSc.@simple_fn(petsclib, exact_pres_cx)
 
 exact_vel_ptr  = sol_type == "solkx" ? exact_vel_kx_ptr  : exact_vel_cx_ptr
 exact_pres_ptr = sol_type == "solkx" ? exact_pres_kx_ptr : exact_pres_cx_ptr
diff --git a/scripts/renames.jl b/scripts/renames.jl
index 2c239d566..a7cba9363 100644
--- a/scripts/renames.jl
+++ b/scripts/renames.jl
@@ -741,6 +741,33 @@ const EXTRA_SHIMS = NamedTuple{(:name, :code, :test), Tuple{String, String, Stri
         end
         """,
     ),
+    (
+        name = "pointwise macros without a library",
+        code = """
+        # Warns in src/dmplex.jl: the one-argument-shorter forms of @simple_fn,
+        # @residual_fn, @jacobian_fn and @bd_fn take PetscInt, PetscScalar and
+        # PetscReal from the calling module.
+        """,
+        test = """
+        @testset "pointwise macros without a library" begin
+            m = Module()
+            Core.eval(m, quote
+                using PETSc
+                const PetscInt = Int64
+                const PetscScalar = Float64
+                const PetscReal = Float64
+                one_u(t, x, u, ctx) = (u[1] = x[1] + 1; nothing)
+            end)
+            ptr = warns(() -> Core.eval(m, :(PETSc.@simple_fn(one_u))))
+            x = [2.0]
+            u = [0.0]
+            err = ccall(ptr, Cint, (Int64, Float64, Ptr{Float64}, Int64, Ptr{Float64}, Ptr{Cvoid}),
+                        1, 0.0, x, 1, u, C_NULL)
+            @test err == 0
+            @test u[1] == 3
+        end
+        """,
+    ),
 ]
 
 # ---------------------------------------------------------------------------
@@ -752,6 +779,13 @@ const INTERNAL = Set{Symbol}([
     :checkout_local_arrays,
     # the per-library body of the MatShell constructor
     :create_matshell,
+    # the expansions behind the pointwise macros and their type lookup
+    :simple_fn_expr,
+    :residual_fn_expr,
+    :jacobian_fn_expr,
+    :bd_fn_expr,
+    :pointwise_types,
+    :pointwise_types_from_scope,
     # leading underscore dropped (naming.md "Internals")
     :petsc_link,
     :petsc_subst,
diff --git a/src/deprecations.jl b/src/deprecations.jl
index 04127261d..ea44a7b81 100644
--- a/src/deprecations.jl
+++ b/src/deprecations.jl
@@ -432,6 +432,11 @@ end
 # Warns in src/dm.jl and returns a copy; `with_local_coordinate_array` hands
 # the array back when its block returns.
 
+# pointwise macros without a library
+# Warns in src/dmplex.jl: the one-argument-shorter forms of @simple_fn,
+# @residual_fn, @jacobian_fn and @bd_fn take PetscInt, PetscScalar and
+# PetscReal from the calling module.
+
 
 # Renamed types: a `const` alias, so `isa` and dispatch keep working
 
diff --git a/src/dmplex.jl b/src/dmplex.jl
index f7ee17a91..962d0c922 100644
--- a/src/dmplex.jl
+++ b/src/dmplex.jl
@@ -13,20 +13,42 @@ Base.unsafe_convert(::Type{LibPETSc.PetscDS}, v::PetscDS) = v.ptr
 # These macros generate a thin C-callable wrapper that unwraps all C pointer
 # arguments into Julia Vectors, then calls the user's pure-Julia function.
 #
-# PetscInt, PetscReal, PetscScalar must be concrete types in the calling scope
-# (e.g. `const PetscInt = petsclib.PetscInt`).
+# The first argument names the library (a PetscLibType value or type). It is
+# evaluated when the macro expands, so it must be a global defined before the
+# macro call; the C signature uses that library's PetscInt, PetscReal and
+# PetscScalar.
 #
 # Usage:
+#   const petsclib = PETSc.getlib(; PetscScalar = Float64)
 #   function my_fn(t, x, u, ctx);  u[1] = x[1]^2 + x[2]^2;  end
-#   const my_fn_ptr = PETSc.@simple_fn(my_fn)
+#   const my_fn_ptr = PETSc.@simple_fn(petsclib, my_fn)
 #
 #   function f0(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
 #               aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, out)
 #       out .= 4
 #   end
-#   const f0_ptr = PETSc.@residual_fn(f0, Nf)        # outsz = Nf
-#   const f1_ptr = PETSc.@residual_fn(f1, dim_*Nf)   # outsz = dim*Nf
-#   const g3_ptr = PETSc.@jacobian_fn(g3, dim_*dim_) # outsz = dim² (Nf=1)
+#   const f0_ptr = PETSc.@residual_fn(petsclib, f0, Nf)        # outsz = Nf
+#   const f1_ptr = PETSc.@residual_fn(petsclib, f1, dim_*Nf)   # outsz = dim*Nf
+#   const g3_ptr = PETSc.@jacobian_fn(petsclib, g3, dim_*dim_) # outsz = dim² (Nf=1)
+#
+# The wrapper's C pointer comes from an `@eval` at run time: within a single
+# top-level block the `@cfunction` is compiled before the wrapper exists.
+
+# Look up the PetscInt, PetscScalar and PetscReal of the library `lib` names
+# in module `mod`.
+function pointwise_types(mod::Module, lib)
+    L = Core.eval(mod, lib)
+    L isa LibPETSc.PetscLibType || (L isa Type && L <: LibPETSc.PetscLibType && isconcretetype(L)) ||
+        throw(ArgumentError("the first argument of a pointwise macro must name a PETSc library or its type, got `$lib` = $L"))
+    return LibPETSc.inttype(L), LibPETSc.scalartype(L), LibPETSc.realtype(L)
+end
+
+# Old form: the types are the globals PetscInt, PetscScalar and PetscReal of `mod`.
+# TODO: remove in v0.6 with the other 0.5 deprecations.
+function pointwise_types_from_scope(mod::Module, name::Symbol, f)
+    @warn "$name without a library is deprecated, use $name(petsclib, $f, …)" maxlog = 1 _id = name
+    return Core.eval(mod, :PetscInt), Core.eval(mod, :PetscScalar), Core.eval(mod, :PetscReal)
+end
 
 # Helper: substitute gensymmed param names into an output-size expression.
 petsc_subst(ex, m) =
@@ -35,7 +57,7 @@ petsc_subst(ex, m) =
     ex
 
 """
-    @simple_fn(f)
+    @simple_fn(petsclib, f)
 
 Generate a C-callable wrapper for `f(t, x, u, ctx)` and return its C pointer,
 matching the `PetscSimplePointFn` signature.  Arguments seen by `f`:
@@ -45,14 +67,20 @@ matching the `PetscSimplePointFn` signature.  Arguments seen by `f`:
   - `ctx` — user context (`Ptr{Cvoid}`)
 
 Used by `project_function!`, `add_boundary!`, `set_exact_solution!`.
-`PetscInt`, `PetscReal`, `PetscScalar` must be in scope.
+`petsclib` is a PETSc library or its type, and must be a global defined before
+the macro call: it sets the C types of the arguments.
 """
+macro simple_fn(lib, f)
+    return simple_fn_expr(f, pointwise_types(__module__, lib)...)
+end
+
 macro simple_fn(f)
+    return simple_fn_expr(f, pointwise_types_from_scope(__module__, Symbol("@simple_fn"), f)...)
+end
+
+function simple_fn_expr(f, PI, PS, PR)
     cfn = Symbol(f, :__cfn__)
     gdim = gensym("dim"); gNc = gensym("Nc")
-    PI = Core.eval(__module__, :PetscInt)
-    PS = Core.eval(__module__, :PetscScalar)
-    PR = Core.eval(__module__, :PetscReal)
     cfn_ptr_expr = :(Base.@cfunction($cfn, Cint,
         ($PI, $PR, Ptr{$PR}, $PI, Ptr{$PS}, Ptr{Cvoid})))
     quote
@@ -70,7 +98,7 @@ macro simple_fn(f)
 end
 
 """
-    @residual_fn(f, outsz)
+    @residual_fn(petsclib, f, outsz)
 
 Generate a C-callable wrapper for the pure-Julia residual function `f` and
 return its C pointer, matching the `PetscPointFn` signature.  All pointer
@@ -85,16 +113,21 @@ User function signature:
 f(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
   aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, out)
 ```
-`PetscInt`, `PetscReal`, `PetscScalar` must be in scope.
+`petsclib` is as in `@simple_fn`.
 """
+macro residual_fn(lib, f, outsz)
+    return residual_fn_expr(f, outsz, pointwise_types(__module__, lib)...)
+end
+
 macro residual_fn(f, outsz)
+    return residual_fn_expr(f, outsz, pointwise_types_from_scope(__module__, Symbol("@residual_fn"), f)...)
+end
+
+function residual_fn_expr(f, outsz, PI, PS, PR)
     cfn = Symbol(f, :__cfn__)
     gdim = gensym("dim"); gNf = gensym("Nf"); gNfAux = gensym("NfAux"); gnC = gensym("nC")
     m = Dict{Symbol,Any}(:dim_ => gdim, :Nf => gNf, :NfAux => gNfAux, :numConstants => gnC)
     outsz_s = petsc_subst(outsz, m)
-    PI = Core.eval(__module__, :PetscInt)
-    PS = Core.eval(__module__, :PetscScalar)
-    PR = Core.eval(__module__, :PetscReal)
     cfn_ptr_expr = :(Base.@cfunction($cfn, Cvoid,
         ($PI, $PI, $PI, Ptr{$PI}, Ptr{$PI},
          Ptr{$PS}, Ptr{$PS}, Ptr{$PS}, Ptr{$PI}, Ptr{$PI},
@@ -132,7 +165,7 @@ macro residual_fn(f, outsz)
 end
 
 """
-    @jacobian_fn(f, outsz)
+    @jacobian_fn(petsclib, f, outsz)
 
 Generate a C-callable wrapper for the pure-Julia Jacobian function `f` and
 return its C pointer, matching the `PetscPointJacFn` signature.  Identical to
@@ -143,16 +176,21 @@ User function signature:
 f(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
   aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, out)
 ```
-`PetscInt`, `PetscReal`, `PetscScalar` must be in scope.
+`petsclib` is as in `@simple_fn`.
 """
+macro jacobian_fn(lib, f, outsz)
+    return jacobian_fn_expr(f, outsz, pointwise_types(__module__, lib)...)
+end
+
 macro jacobian_fn(f, outsz)
+    return jacobian_fn_expr(f, outsz, pointwise_types_from_scope(__module__, Symbol("@jacobian_fn"), f)...)
+end
+
+function jacobian_fn_expr(f, outsz, PI, PS, PR)
     cfn = Symbol(f, :__cfn__)
     gdim = gensym("dim"); gNf = gensym("Nf"); gNfAux = gensym("NfAux"); gnC = gensym("nC")
     m = Dict{Symbol,Any}(:dim_ => gdim, :Nf => gNf, :NfAux => gNfAux, :numConstants => gnC)
     outsz_s = petsc_subst(outsz, m)
-    PI = Core.eval(__module__, :PetscInt)
-    PS = Core.eval(__module__, :PetscScalar)
-    PR = Core.eval(__module__, :PetscReal)
     cfn_ptr_expr = :(Base.@cfunction($cfn, Cvoid,
         ($PI, $PI, $PI, Ptr{$PI}, Ptr{$PI},
          Ptr{$PS}, Ptr{$PS}, Ptr{$PS}, Ptr{$PI}, Ptr{$PI},
@@ -191,7 +229,7 @@ macro jacobian_fn(f, outsz)
 end
 
 """
-    @bd_fn(f, outsz)
+    @bd_fn(petsclib, f, outsz)
 
 Generate a C-callable wrapper for the pure-Julia boundary pointwise function `f`
 and return its C pointer, matching the `PetscBdPointFn` signature.  Identical to
@@ -203,16 +241,21 @@ User function signature:
 f(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
   aOff, aOff_x, a, a_t, a_x, t, x, n, numConstants, constants, out)
 ```
-`PetscInt`, `PetscReal`, `PetscScalar` must be in scope.
+`petsclib` is as in `@simple_fn`.
 """
+macro bd_fn(lib, f, outsz)
+    return bd_fn_expr(f, outsz, pointwise_types(__module__, lib)...)
+end
+
 macro bd_fn(f, outsz)
+    return bd_fn_expr(f, outsz, pointwise_types_from_scope(__module__, Symbol("@bd_fn"), f)...)
+end
+
+function bd_fn_expr(f, outsz, PI, PS, PR)
     cfn = Symbol(f, :__cfn__)
     gdim = gensym("dim"); gNf = gensym("Nf"); gNfAux = gensym("NfAux"); gnC = gensym("nC")
     m = Dict{Symbol,Any}(:dim_ => gdim, :Nf => gNf, :NfAux => gNfAux, :numConstants => gnC)
     outsz_s = petsc_subst(outsz, m)
-    PI = Core.eval(__module__, :PetscInt)
-    PS = Core.eval(__module__, :PetscScalar)
-    PR = Core.eval(__module__, :PetscReal)
     cfn_ptr_expr = :(Base.@cfunction($cfn, Cvoid,
         ($PI, $PI, $PI, Ptr{$PI}, Ptr{$PI},
          Ptr{$PS}, Ptr{$PS}, Ptr{$PS}, Ptr{$PI}, Ptr{$PI},
diff --git a/test/core/deprecations.jl b/test/core/deprecations.jl
index 69e418bde..fac6caa63 100644
--- a/test/core/deprecations.jl
+++ b/test/core/deprecations.jl
@@ -380,4 +380,22 @@ const SHIMS = [
         @test coord[1, 3] ≈ 1
     end
 
+    @testset "pointwise macros without a library" begin
+        m = Module()
+        Core.eval(m, quote
+            using PETSc
+            const PetscInt = Int64
+            const PetscScalar = Float64
+            const PetscReal = Float64
+            one_u(t, x, u, ctx) = (u[1] = x[1] + 1; nothing)
+        end)
+        ptr = warns(() -> Core.eval(m, :(PETSc.@simple_fn(one_u))))
+        x = [2.0]
+        u = [0.0]
+        err = ccall(ptr, Cint, (Int64, Float64, Ptr{Float64}, Int64, Ptr{Float64}, Ptr{Cvoid}),
+                    1, 0.0, x, 1, u, C_NULL)
+        @test err == 0
+        @test u[1] == 3
+    end
+
 end
diff --git a/test/dm/dmplex.jl b/test/dm/dmplex.jl
index f47a29d9e..e9fc1814c 100644
--- a/test/dm/dmplex.jl
+++ b/test/dm/dmplex.jl
@@ -10,18 +10,15 @@ using MPI
 
 MPI.Initialized() || MPI.Init()
 
-# ── PETSc library and scalar types ───────────────────────────────────────────
-# @simple_fn / @residual_fn / @jacobian_fn resolve
-# PetscInt / PetscScalar / PetscReal from the calling module at macro-expansion
-# time, so they MUST be module-level names before the macros are invoked.
-# We define them once for the Float64/Int64 lib that owns the FEM callbacks.
+# ── PETSc library ────────────────────────────────────────────────────────────
+# @simple_fn / @residual_fn / @jacobian_fn evaluate their library argument at
+# macro-expansion time, so it must be a module-level name before the macros
+# are invoked. The FEM callbacks use the Float64/Int64 lib.
 
 const _dmplex_petsclib = PETSc.petsclibs[1]
 PETSc.initialize(_dmplex_petsclib)
 
 const PetscInt    = _dmplex_petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 # Serial-safe communicator (shared across all lib loops)
 const _TC = MPI.COMM_SELF
@@ -46,7 +43,7 @@ function _dm_f0(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                 aOff, aOff_x, a, a_t, a_x, t, x, nC, cst, f0)
     f0[1] = 0.0
 end
-const _dm_f0_ptr = PETSc.@residual_fn(_dm_f0, 1)
+const _dm_f0_ptr = PETSc.@residual_fn(_dmplex_petsclib, _dm_f0, 1)
 
 # f1: diffusion flux ∇u
 function _dm_f1(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -55,7 +52,7 @@ function _dm_f1(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[d] = u_x[d]
     end
 end
-const _dm_f1_ptr = PETSc.@residual_fn(_dm_f1, dim_)
+const _dm_f1_ptr = PETSc.@residual_fn(_dmplex_petsclib, _dm_f1, dim_)
 
 # g3: identity stiffness tensor (Laplacian Jacobian)
 function _dm_g3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -65,7 +62,7 @@ function _dm_g3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[d * dim_ + d + 1] = 1.0
     end
 end
-const _dm_g3_ptr = PETSc.@jacobian_fn(_dm_g3, dim_ * dim_)
+const _dm_g3_ptr = PETSc.@jacobian_fn(_dmplex_petsclib, _dm_g3, dim_ * dim_)
 
 # exact solution: u(t, x) = x₁ + x₂  (linear → exact in P1/Q1)
 function _dm_exact(t, x, u, ctx)
@@ -74,13 +71,17 @@ function _dm_exact(t, x, u, ctx)
         u[1] += x[d]
     end
 end
-const _dm_exact_ptr = PETSc.@simple_fn(_dm_exact)
+const _dm_exact_ptr = PETSc.@simple_fn(_dmplex_petsclib, _dm_exact)
 
 # zero BC function
 function _dm_zero(t, x, u, ctx)
     fill!(u, 0.0)
 end
-const _dm_zero_ptr = PETSc.@simple_fn(_dm_zero)
+const _dm_zero_ptr = PETSc.@simple_fn(_dmplex_petsclib, _dm_zero)
+
+# The library's type names the C types as well as the library does.
+_dm_one(t, x, u, ctx) = (u[1] = x[1] + 1; nothing)
+const _dm_one_ptr = PETSc.@simple_fn(typeof(_dmplex_petsclib), _dm_one)
 
 # ─────────────────────────────────────────────────────────────────────────────
 @testset "DMPlex" begin
@@ -89,6 +90,18 @@ if isintelmac
     @test_skip "DMPlex+PETSc tests skipped on macOS Intel (known SIGSEGV in PETSc_jll)"
 else
 
+@testset "pointwise macro library argument" begin
+    @test _dmplex_petsclib.PetscScalar === Float64
+    x = [2.0]
+    u = [0.0]
+    err = ccall(_dm_one_ptr, Cint, (PetscInt, Float64, Ptr{Float64}, PetscInt, Ptr{Float64}, Ptr{Cvoid}),
+                1, 0, x, 1, u, C_NULL)
+    @test err == 0
+    @test u[1] == 3
+    @test_throws "must name a PETSc library" macroexpand(@__MODULE__, :(PETSc.@simple_fn(1, _dm_one)))
+end
+
+
 for petsclib in PETSc.petsclibs
     PETSc.initialize(petsclib)
     PetscInt_t    = petsclib.PetscInt