Error with ForwardColorJacCache

[bmxam]

I usually use ForwardColorJacCache to evaluate my jacobians. But I'm struggling with the following MWE. The last line shows that, when using ForwardColorJacCache, NaN are obtained (run 2 or 3 times because sometimes you're lucky). @ghislainb do you see something wrong?

I have tried to eliminate Bcube from the MWE but did not reproduce the problem (although I doubt this is a Bcube problem...). I'm running it with Julia 1.11

module tmp
using Bcube
using StaticArrays
using LinearAlgebra
using SparseDiffTools
using ForwardDiff
using UnPack
using SparseArrays

function eval_rhs(dq, q, p)
    @unpack U, V, dΩ = p

    h, hu = FEFunction(U, q)
    l((v_h, v_hu)) = ∫(zero(h) * (v_h + v_hu))dΩ

    assemble_linear!(dq, l, V)
end

function run_fem()
    n_fem = 18
    degree = 1

    # Mesh
    mesh = line_mesh(n_fem; xmin=-1.0, xmax=1.0)
    dΩ = Measure(CellDomain(mesh), 2 * degree + 1)

    # FEspaces
    fs_h = FunctionSpace(:Lagrange, degree)
    U_h = TrialFESpace(fs_h, mesh) # continuous !
    fs_hu = FunctionSpace(:Lagrange, degree)
    U_hu = TrialFESpace(fs_hu, mesh) # continuous !
    V_h = TestFESpace(U_h)
    V_hu = TestFESpace(U_hu)
    U = MultiFESpace(U_h, U_hu)
    V = MultiFESpace(V_h, V_hu)
    q = FEFunction(U)
    h, hu = q

    # Init
    projection_l2!(h, PhysicalFunction(x -> x[1] < 0 ? 3.0 : 2.0), mesh)

    # Cache
    q0 = get_dof_values(q)
    p = (; U, V, dΩ)
    f!(_y, _x) = eval_rhs(_y, _x, p)

    sparsity = Bcube.build_jacobian_sparsity_pattern(U, mesh)

    J = copy(sparsity)

    colorvec = matrix_colors(sparsity)
    fcjc = ForwardColorJacCache(f!, copy(q0); colorvec, sparsity)

    forwarddiff_color_jacobian!(J, f!, q0; colorvec, sparsity)
    @show all(isfinite.(J.nzval)) # OK

    forwarddiff_color_jacobian!(J, f!, q0, fcjc)
    @show all(isfinite.(J.nzval)) # KO
end


run_fem()

end