NonlinearSolve.jl with Sparse Jacobian: J isn't square?!

General Usage
1

I’m trying to use NonlinearSolve.jl’s LevenbergMarquardt solver on a powerflow. The solve call yields an error: LoadError: ArgumentError: KLU only accepts square matrices. That seems nonsensical: J0 is square. I expect things to have dimension equal to twice the number of buses in the system: in this case, that’ll be 10. The print statements I’ve added bear that out. What’s the issue?

Possibly relevant: sparsity detection

using NonlinearSolve
using SciMLBase
using SparseArrays
using PowerFlows, PowerSystemCaseBuilder
const PF = PowerFlows
const PSB = PowerSystemCaseBuilder
# small system for demo purposes: actual target system has 10k buses. 
sys = PSB.build_system(PSB.MatpowerTestSystems, "matpower_case5_sys")
pf = PF.ACPowerFlow()
data = PF.PowerFlowData(pf, sys)
residual = PF.ACPowerFlowResidual(data, 1)
J = PF.ACPowerFlowJacobian(data, 1) # this initializes the sparse structure.
J0 = deepcopy(J.Jv) # take a copy of the sparse structure.
x0 = PF.calculate_x0(data, 1)
@assert size(J0, 1) == size(J0, 2) # yes, dimensions match. 
@assert size(J0, 1) == size(x0, 1)

function f!(du::Vector{Float64}, u::Vector{Float64}, ::Any)
    println(size(du))
    println(size(u))
    residual(du, u, 1) # in-place compute f at u
end

function jac!(A::SparseMatrixCSC{Float64, Int32}, u::Vector{Float64}, ::Any)
    println(size(u))
    println(size(A))
    PF.update_data!(data, u, 1)
    J(A, 1) # in-place compute J at the point stored in data (set to u above)
end

func = NonlinearFunction{true, SciMLBase.AutoSpecialize}(
    f!;
    jac = jac!,
    jac_prototype = J0,
)

prob = NonlinearProblem(func, x0)
solver = NonlinearSolveFirstOrder.LevenbergMarquardt()
solve(prob, solver)