Skip to content

switch to SparseConnectivityTracer and DifferentiationInterface #78

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 6 commits into
base: master
Choose a base branch
from

Conversation

chmerdon
Copy link
Member

@chmerdon chmerdon commented Aug 4, 2025

sparsity detection for BilinearOperator and jacobian preparation for NonlinearOperator now is done with SparseConnectivityTracer and DifferentiationInterface (as suggested by issue #50); SparseDiffTools and Symbolics removed as a main dependency

chmerdon added 4 commits August 4, 2025 23:16
…ntiationInterface in jacobian preparation in NonlinearOperator
…arseConnectivityTracer, could removed Symbolics as a main dependence now, version bump + changelog
@chmerdon chmerdon force-pushed the chmerdon/diffinterface branch from f4f4a38 to c5b24b6 Compare August 4, 2025 15:29
@chmerdon chmerdon marked this pull request as ready for review August 4, 2025 15:47
@chmerdon chmerdon requested review from j-fu and pjaap August 6, 2025 06:07
@gdalle
Copy link

gdalle commented Aug 6, 2025

Hey there! Just popping by to ask how the transition went?

@@ -8,6 +8,7 @@ using Metis
using Aqua
using Triangulate
using SimplexGridFactory
using Symbolics

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why add symbolic here?

@j-fu j-fu requested a review from Copilot August 6, 2025 18:39
Copy link

@Copilot Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull Request Overview

This pull request modernizes the automatic differentiation infrastructure by switching from SparseDiffTools and Symbolics to SparseConnectivityTracer and DifferentiationInterface for sparsity detection and jacobian preparation.

  • Replaces Symbolics.jacobian_sparsity with SparseConnectivityTracer's jacobian_sparsity function
  • Refactors NonlinearOperator to use DifferentiationInterface for jacobian computation with sparse backends
  • Adds new autodiff_backend parameter for configurable automatic differentiation backends

Reviewed Changes

Copilot reviewed 6 out of 7 changed files in this pull request and generated 3 comments.

Show a summary per file
File Description
src/ExtendableFEM.jl Updates imports to use new AD packages and removes old dependencies
src/common_operators/nonlinear_operator.jl Major refactor of KernelEvaluator and jacobian computation using DifferentiationInterface
src/common_operators/bilinear_operator.jl Replaces Symbolics sparsity detection with SparseConnectivityTracer
src/common_operators/bilinear_operator_dg.jl Replaces Symbolics sparsity detection with SparseConnectivityTracer
test/runtests.jl Adds Symbolics import for test compatibility
CHANGELOG.md Documents the changes and new parameter

else
jac_sparsity_pattern = O.parameters[:sparse_jacobians_pattern]
end
jac = Tv.(sparse(sparsity_pattern))
Copy link
Preview

Copilot AI Aug 6, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The sparse_jacobians variable is hardcoded to false, which contradicts the parameter O.parameters[:sparse_jacobians] used earlier in the function. This will prevent sparse jacobian functionality from working correctly.

Suggested change
jac = Tv.(sparse(sparsity_pattern))
sparse_jacobians = O.parameters[:sparse_jacobians]
if sparse_jacobians
if O.parameters[:sparse_jacobians_pattern] === nothing
jac_sparsity_pattern = DifferentiationInterface.sparsity_pattern(jac_prep)
else
jac_sparsity_pattern = O.parameters[:sparse_jacobians_pattern]
end
jac = Tv.(sparse(jac_sparsity_pattern))

Copilot uses AI. Check for mistakes.

else
jac_sparsity_pattern = O.parameters[:sparse_jacobians_pattern]
end
jac = Tv.(sparse(sparsity_pattern))
Copy link
Preview

Copilot AI Aug 6, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The variable 'sparsity_pattern' is undefined in this scope. It should be 'jac_sparsity_pattern' which is defined above.

Suggested change
jac = Tv.(sparse(sparsity_pattern))
jac = Tv.(sparse(jac_sparsity_pattern))

Copilot uses AI. Check for mistakes.

@@ -8,6 +8,7 @@ using Metis
using Aqua
using Triangulate
using SimplexGridFactory
using Symbolics
Copy link
Preview

Copilot AI Aug 6, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Adding Symbolics as a test dependency contradicts the PR's purpose of removing Symbolics as a main dependency. Consider refactoring tests to use the new SparseConnectivityTracer instead.

Suggested change
using Symbolics

Copilot uses AI. Check for mistakes.

@j-fu
Copy link
Member

j-fu commented Aug 6, 2025

Bin im Urlaub. Für mich ok wenn die Tests laufen

@pjaap
Copy link
Member

pjaap commented Aug 7, 2025

It currently breaks one of my project examples, similar to Example330, which fails in the docs. I will look into it.

@chmerdon
Copy link
Member Author

chmerdon commented Aug 7, 2025

@pjaap thanks!

@gdalle thanks for popping by and the great packages! So far the transition went well, we just have a problem with some examples were we need to sparse-detect and differentiate functionals that are ForwardDiff jacobians by themself, see e.g. here. I also did not yet compare runtimes and allocations and other backends than ForwardDiff, but will do this later maybe.

@j-fu Symbolics is still needed in some of the examples (some have tests) where given functions are differentiated symbolically to construct analytic benchmarks. Dann erstmal schönen Urlaub! Ich bin auch noch im Urlaub.

@gdalle
Copy link

gdalle commented Aug 7, 2025

@chmerdon do you have an MWE demonstrating the failure you're talking about? What code should I run and in which environment? What's the error message?

@pjaap
Copy link
Member

pjaap commented Aug 11, 2025

Sorry for the delay, I fell sick last week....

Here is a MWE for the problem. It is a strip-down of Example330.

Basically, it should simply minimize

$$\int u^2 - 2u \, dx$$

with the optimal solution

$$u(x) = 1$$

In our FE Library, we solve weak formulations, hence we solve

$$\int DF(u) \cdot v \, dx = 0 \quad \forall v$$

with

$$DF(u) = D[u^2 - 2u] = 2u - 2$$

for the linear operator. This operator is computed automatically by ForwardDiff here.

module MWE

using ExtendableFEM
using ExtendableGrids
using DiffResults
using ForwardDiff


function nonlinear_F!(result, input, qpinfo)
    result[1] = input[1]^2 - 2*input[1] # => input = 1 is minimizer
end

## derivative of nonlinear_F via AD
function nonlinear_DF!()
    Dresult = nothing
    cfg = nothing
    result_dummy = nothing
    F(qpinfo) = (a, b) -> nonlinear_F!(a, b, qpinfo)

    return function closure(result, input, qpinfo)
        if Dresult === nothing
            ## first initialization of DResult when type of input = F is known
            result_dummy = zeros(eltype(input), 1)
            Dresult = DiffResults.JacobianResult(result_dummy, input)
            cfg = ForwardDiff.JacobianConfig(F(qpinfo), result_dummy, input, ForwardDiff.Chunk{length(input)}())
        end
        Dresult = ForwardDiff.vector_mode_jacobian!(Dresult, F(qpinfo), result_dummy, input, cfg)
        copyto!(result, DiffResults.jacobian(Dresult))
        return nothing
    end
end


function main()

    ## define unknowns
    u = Unknown("u"; name = "potential")

    ## define problem
    PD = ProblemDescription()
    assign_unknown!(PD, u)
    assign_operator!(PD, NonlinearOperator(nonlinear_DF!(), [id(u)]; sparse_jacobians = false))

    ## simple 2D grid
    xgrid = grid_unitsquare(Triangle2D)

    ## solve
    FES = FESpace{H1P1{1}}(xgrid)
    sol = solve(PD, FES; maxiterations = 20)

    # should be all ≈ 1
    @show sol[u].entries

    return sol
end

end # module

@pjaap
Copy link
Member

pjaap commented Aug 11, 2025

The MWE works fine on current master, but with SCT we get the following error:

julia> MWE.main()
┌ Info: SOLVING My problem @ time = 0.0
│                               unknowns = ["potential"]
│                               fetypes = ["H1P1{1}"]
└                               ndofs = [5]
 #IT    ------- RESIDUALS -------
        NONLINEAR       LINEAR
ERROR: MethodError: no method matching ForwardDiff.Dual{…}(::ForwardDiff.Dual{…}, ::ForwardDiff.Partials{…})
The type `ForwardDiff.Dual{ForwardDiff.Tag{Main.MWE.var"#nonlinear_DF!##1#nonlinear_DF!##2"{ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}, 1}` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  ForwardDiff.Dual{T, V, N}(::ForwardDiff.Dual{T, V, N}) where {T, V, N}
   @ ForwardDiff ~/.julia/packages/ForwardDiff/Wq9Wb/src/dual.jl:76
  ForwardDiff.Dual{T, V, N}(::Number) where {T, V, N}
   @ ForwardDiff ~/.julia/packages/ForwardDiff/Wq9Wb/src/dual.jl:78
  ForwardDiff.Dual{T, V, N}(::Any) where {T, V, N}
   @ ForwardDiff ~/.julia/packages/ForwardDiff/Wq9Wb/src/dual.jl:77
  ...

Stacktrace:
  [1] seed!
    @ ~/.julia/packages/ForwardDiff/Wq9Wb/src/apiutils.jl:84 [inlined]
  [2] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/Wq9Wb/src/apiutils.jl:29 [inlined]
  [3] vector_mode_jacobian!(result::DiffResults.MutableDiffResult{…}, f!::Main.MWE.var"#nonlinear_DF!##1#nonlinear_DF!##2"{…}, y::Vector{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/Wq9Wb/src/jacobian.jl:157
  [4] (::Main.MWE.var"#closure#nonlinear_DF!##3"{…})(result::Vector{…}, input::Vector{…}, qpinfo::ExtendableFEMBase.QPInfos{…})
    @ Main.MWE ~/git-stuff/julia/ExtendableFEM/examples/MWE.jl:27
  [5] #1501
    @ ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:236 [inlined]
  [6] compute_ydual_twoarg
    @ ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceForwardDiffExt/twoarg.jl:59 [inlined]
  [7] pushforward!
    @ ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceForwardDiffExt/twoarg.jl:122 [inlined]
  [8] _sparse_jacobian_aux!(::Tuple{…}, ::Matrix{…}, ::DifferentiationInterfaceSparseMatrixColoringsExt.SMCPushforwardSparseJacobianPrep{…}, ::ADTypes.AutoSparse{…}, ::Vector{…})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:308
  [9] jacobian!
    @ ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:237 [inlined]
 [10] value_and_jacobian!(::ExtendableFEM.var"#1501#1502"{…}, ::Vector{…}, ::Matrix{…}, ::DifferentiationInterfaceSparseMatrixColoringsExt.SMCPushforwardSparseJacobianPrep{…}, ::ADTypes.AutoSparse{…}, ::Vector{…})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:277
 [11] (::ExtendableFEM.var"#assembly_loop#1511"{…})(A::ExtendableSparse.ExtendableSparseMatrixCSC{…}, b::Vector{…}, sol::Vector{…}, items::SubArray{…}, EG::Type{…}, QF::ExtendableFEMBase.SQuadratureRule{…}, BE_test::Vector{…}, BE_args::Vector{…}, BE_test_vals::Vector{…}, BE_args_vals::Vector{…}, L2G::ExtendableGrids.L2GTransformer{…}, K::ExtendableFEM.KernelEvaluator{…}, part::Int64; time::Float64)
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:357
 [12] assembly_loop
    @ ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:271 [inlined]
 [13] macro expansion
    @ ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:464 [inlined]
 [14] macro expansion
    @ ./timing.jl:461 [inlined]
 [15] (::ExtendableFEM.var"#assembler#1517"{…})(A::ExtendableSparse.ExtendableSparseMatrixCSC{…}, b::Vector{…}, sol::Vector{…}; kwargs::@Kwargs{…})
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:440
 [16] assemble!(A::ExtendableFEMBase.FEMatrix{…}, b::ExtendableFEMBase.FEVector{…}, sol::ExtendableFEMBase.FEVector{…}, O::ExtendableFEM.NonlinearOperator{…}, SC::ExtendableFEM.SolverConfiguration{…}; assemble_matrix::Bool, assemble_rhs::Bool, time::Float64, kwargs::@Kwargs{…})
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:490
 [17] macro expansion
    @ ./timing.jl:645 [inlined]
 [18] assemble_system!(A::ExtendableFEMBase.FEMatrix{…}, b::ExtendableFEMBase.FEVector{…}, sol::ExtendableFEMBase.FEVector{…}, PD::ExtendableFEM.ProblemDescription, SC::ExtendableFEM.SolverConfiguration{…}, timer::TimerOutputs.TimerOutput; kwargs::@Kwargs{…})
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/solvers.jl:134
 [19] solve(PD::ExtendableFEM.ProblemDescription, FES::ExtendableFEMBase.FESpace{…}, SC::Nothing; unknowns::Vector{…}, kwargs::@Kwargs{…})
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/solvers.jl:533
 [20] solve (repeats 2 times)
    @ ~/git-stuff/julia/ExtendableFEM/src/solvers.jl:469 [inlined]
 [21] main()
    @ Main.MWE ~/git-stuff/julia/ExtendableFEM/examples/MWE.jl:49
 [22] top-level scope
    @ REPL[5]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> show(err)
1-element ExceptionStack:
MethodError: no method matching ForwardDiff.Dual{ForwardDiff.Tag{Main.MWE.var"#nonlinear_DF!##1#nonlinear_DF!##2"{ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}, 1}(::ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}, ::ForwardDiff.Partials{1, SparseConnectivityTracer.GradientTracer{Int64, BitSet}})
The type `ForwardDiff.Dual{ForwardDiff.Tag{Main.MWE.var"#nonlinear_DF!##1#nonlinear_DF!##2"{ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}, 1}` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  ForwardDiff.Dual{T, V, N}(::ForwardDiff.Dual{T, V, N}) where {T, V, N}
   @ ForwardDiff ~/.julia/packages/ForwardDiff/Wq9Wb/src/dual.jl:76
  ForwardDiff.Dual{T, V, N}(::Number) where {T, V, N}
   @ ForwardDiff ~/.julia/packages/ForwardDiff/Wq9Wb/src/dual.jl:78
  ForwardDiff.Dual{T, V, N}(::Any) where {T, V, N}
   @ ForwardDiff ~/.julia/packages/ForwardDiff/Wq9Wb/src/dual.jl:77
  ...

Stacktrace:
  [1] seed!
    @ ~/.julia/packages/ForwardDiff/Wq9Wb/src/apiutils.jl:84 [inlined]
  [2] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/Wq9Wb/src/apiutils.jl:29 [inlined]
  [3] vector_mode_jacobian!(result::DiffResults.MutableDiffResult{1, Vector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}}, Tuple{Matrix{SparseConnectivityTracer.GradientTracer{Int64, BitSet}}}}, f!::Main.MWE.var"#nonlinear_DF!##1#nonlinear_DF!##2"{ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, y::Vector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}}, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{Main.MWE.var"#nonlinear_DF!##1#nonlinear_DF!##2"{ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}, 1, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{Main.MWE.var"#nonlinear_DF!##1#nonlinear_DF!##2"{ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Main.MWE.var"#nonlinear_DF!##1#nonlinear_DF!##2"{ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}}, SparseConnectivityTracer.GradientTracer{Int64, BitSet}, 1}}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/Wq9Wb/src/jacobian.jl:157
  [4] (::Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"})(result::Vector{ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}}, input::Vector{ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}}, qpinfo::ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing})
    @ Main.MWE ~/git-stuff/julia/ExtendableFEM/examples/MWE.jl:27
  [5] #1501
    @ ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:236 [inlined]
  [6] compute_ydual_twoarg
    @ ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceForwardDiffExt/twoarg.jl:59 [inlined]
  [7] pushforward!
    @ ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceForwardDiffExt/twoarg.jl:122 [inlined]
  [8] _sparse_jacobian_aux!(::Tuple{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Vector{Float64}}, ::Matrix{Float64}, ::DifferentiationInterfaceSparseMatrixColoringsExt.SMCPushforwardSparseJacobianPrep{Tuple{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Vector{Float64}, ADTypes.AutoSparse{ADTypes.AutoForwardDiff{nothing, Nothing}, SparseConnectivityTracer.TracerSparsityDetector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}, SparseConnectivityTracer.HessianTracer{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}, SparseMatrixColorings.GreedyColoringAlgorithm{:direct, SparseMatrixColorings.NaturalOrder}}, Vector{Float64}, Tuple{}}, DifferentiationInterface.BatchSizeSettings{1, true, true}, SparseArrays.SparseMatrixCSC{Bool, Int64}, SparseMatrixColorings.ColumnColoringResult{SparseArrays.SparseMatrixCSC{Bool, Int64}, Int64, SparseMatrixColorings.BipartiteGraph{Int64}, Vector{Int64}, Vector{SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}}, Vector{Int64}, Nothing}, Matrix{Float64}, Vector{Tuple{Vector{Float64}}}, Vector{Tuple{Vector{Float64}}}, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgPushforwardPrep{Tuple{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Vector{Float64}, ADTypes.AutoForwardDiff{nothing, Nothing}, Vector{Float64}, Tuple{Vector{Float64}}, Tuple{}}, ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}}, Tuple{}}}, ::ADTypes.AutoSparse{ADTypes.AutoForwardDiff{nothing, Nothing}, SparseConnectivityTracer.TracerSparsityDetector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}, SparseConnectivityTracer.HessianTracer{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}, SparseMatrixColorings.GreedyColoringAlgorithm{:direct, SparseMatrixColorings.NaturalOrder}}, ::Vector{Float64})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:308
  [9] jacobian!
    @ ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:237 [inlined]
 [10] value_and_jacobian!(::ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, ::Vector{Float64}, ::Matrix{Float64}, ::DifferentiationInterfaceSparseMatrixColoringsExt.SMCPushforwardSparseJacobianPrep{Tuple{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Vector{Float64}, ADTypes.AutoSparse{ADTypes.AutoForwardDiff{nothing, Nothing}, SparseConnectivityTracer.TracerSparsityDetector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}, SparseConnectivityTracer.HessianTracer{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}, SparseMatrixColorings.GreedyColoringAlgorithm{:direct, SparseMatrixColorings.NaturalOrder}}, Vector{Float64}, Tuple{}}, DifferentiationInterface.BatchSizeSettings{1, true, true}, SparseArrays.SparseMatrixCSC{Bool, Int64}, SparseMatrixColorings.ColumnColoringResult{SparseArrays.SparseMatrixCSC{Bool, Int64}, Int64, SparseMatrixColorings.BipartiteGraph{Int64}, Vector{Int64}, Vector{SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}}, Vector{Int64}, Nothing}, Matrix{Float64}, Vector{Tuple{Vector{Float64}}}, Vector{Tuple{Vector{Float64}}}, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgPushforwardPrep{Tuple{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Vector{Float64}, ADTypes.AutoForwardDiff{nothing, Nothing}, Vector{Float64}, Tuple{Vector{Float64}}, Tuple{}}, ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}}, Tuple{}}}, ::ADTypes.AutoSparse{ADTypes.AutoForwardDiff{nothing, Nothing}, SparseConnectivityTracer.TracerSparsityDetector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}, SparseConnectivityTracer.HessianTracer{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}, SparseMatrixColorings.GreedyColoringAlgorithm{:direct, SparseMatrixColorings.NaturalOrder}}, ::Vector{Float64})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/D3LUI/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:277
 [11] (::ExtendableFEM.var"#assembly_loop#1511"{ExtendableFEM.var"#assembly_loop#1461#1512"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, Int64, Int64, Vector{Union{ExtendableGrids.SerialVariableTargetAdjacency{Int32}, ExtendableGrids.Adjacency{Int32}}}, Vector{Union{ExtendableGrids.SerialVariableTargetAdjacency{Int32}, ExtendableGrids.Adjacency{Int32}}}, Bool, Bool, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Bool}, Int64, Int64, Vector{DataType}, Vector{DataType}, Vector{Int32}, Vector{Float64}}})(A::ExtendableSparse.ExtendableSparseMatrixCSC{Float64, Int64}, b::Vector{Float64}, sol::Vector{ExtendableFEMBase.FEVectorBlock{Float64, Float64, Int32, Vector{Float64}, H1P1{1}, ON_CELLS}}, items::SubArray{Int32, 1, Vector{Int32}, Tuple{UnitRange{Int64}}, true}, EG::Type{ExtendableGrids.Triangle2D}, QF::ExtendableFEMBase.SQuadratureRule{Float64, ExtendableGrids.Triangle2D, 2, 3}, BE_test::Vector{ExtendableFEMBase.SingleFEEvaluator{Float64, Float64, Int32, id, H1P1{1}, ExtendableGrids.Triangle2D, ExtendableFEMBase.var"#closure#get_basis##32"{1, Int64}, ExtendableFEMBase.var"#20#21", ExtendableFEMBase.var"#20#21", ExtendableFEMBase.var"#20#21"}}, BE_args::Vector{ExtendableFEMBase.SingleFEEvaluator{Float64, Float64, Int32, id, H1P1{1}, ExtendableGrids.Triangle2D, ExtendableFEMBase.var"#closure#get_basis##32"{1, Int64}, ExtendableFEMBase.var"#20#21", ExtendableFEMBase.var"#20#21", ExtendableFEMBase.var"#20#21"}}, BE_test_vals::Vector{Array{Float64, 3}}, BE_args_vals::Vector{Array{Float64, 3}}, L2G::ExtendableGrids.L2GTransformer{Float64, Int32, ExtendableGrids.Triangle2D, ExtendableGrids.Cartesian2D}, K::ExtendableFEM.KernelEvaluator{Float64, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}, DifferentiationInterfaceSparseMatrixColoringsExt.SMCPushforwardSparseJacobianPrep{Tuple{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Vector{Float64}, ADTypes.AutoSparse{ADTypes.AutoForwardDiff{nothing, Nothing}, SparseConnectivityTracer.TracerSparsityDetector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}, SparseConnectivityTracer.HessianTracer{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}, SparseMatrixColorings.GreedyColoringAlgorithm{:direct, SparseMatrixColorings.NaturalOrder}}, Vector{Float64}, Tuple{}}, DifferentiationInterface.BatchSizeSettings{1, true, true}, SparseArrays.SparseMatrixCSC{Bool, Int64}, SparseMatrixColorings.ColumnColoringResult{SparseArrays.SparseMatrixCSC{Bool, Int64}, Int64, SparseMatrixColorings.BipartiteGraph{Int64}, Vector{Int64}, Vector{SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}}, Vector{Int64}, Nothing}, Matrix{Float64}, Vector{Tuple{Vector{Float64}}}, Vector{Tuple{Vector{Float64}}}, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgPushforwardPrep{Tuple{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Vector{Float64}, ADTypes.AutoForwardDiff{nothing, Nothing}, Vector{Float64}, Tuple{Vector{Float64}}, Tuple{}}, ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}, Float64}, Float64, 1}}, Tuple{}}}, ADTypes.AutoSparse{ADTypes.AutoForwardDiff{nothing, Nothing}, SparseConnectivityTracer.TracerSparsityDetector{SparseConnectivityTracer.GradientTracer{Int64, BitSet}, SparseConnectivityTracer.HessianTracer{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}, SparseMatrixColorings.GreedyColoringAlgorithm{:direct, SparseMatrixColorings.NaturalOrder}}, ExtendableFEM.var"#1501#1502"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, ExtendableFEMBase.QPInfos{Int32, Float64, Float64, Float64, Float64, Float64, Int32, Nothing}}}, part::Int64; time::Float64)
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:357
 [12] assembly_loop
    @ ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:271 [inlined]
 [13] macro expansion
    @ ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:464 [inlined]
 [14] macro expansion
    @ ./timing.jl:461 [inlined]
 [15] (::ExtendableFEM.var"#assembler#1517"{ExtendableFEM.var"#assembler#1464#1518"{ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, Vector{ExtendableFEM.KernelEvaluator}, Vector{DataType}, ExtendableGrids.ExtendableGrid{Float64, Int32}}})(A::ExtendableSparse.ExtendableSparseMatrixCSC{Float64, Int64}, b::Vector{Float64}, sol::Vector{ExtendableFEMBase.FEVectorBlock{Float64, Float64, Int32, Vector{Float64}, H1P1{1}, ON_CELLS}}; kwargs::@Kwargs{time::Float64})
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:440
 [16] assemble!(A::ExtendableFEMBase.FEMatrix{Float64, Int64, Float64, Int32, 1, 1, 1}, b::ExtendableFEMBase.FEVector{Float64, Float64, Int32, Vector{Float64}}, sol::ExtendableFEMBase.FEVector{Float64, Float64, Int32, Vector{Float64}}, O::ExtendableFEM.NonlinearOperator{Float64, ExtendableFEM.Unknown{Symbol}, Main.MWE.var"#closure#nonlinear_DF!##3"{Main.MWE.var"#F#nonlinear_DF!##0"}, Nothing}, SC::ExtendableFEM.SolverConfiguration{ExtendableFEMBase.FEMatrix{Float64, Int64, Float64, Int32, 1, 1, 1}, ExtendableFEMBase.FEVector{Float64, Float64, Int32, Vector{Float64}}, ExtendableFEMBase.FEVector{Float64, Float64, Int32, Vector{Float64}}}; assemble_matrix::Bool, assemble_rhs::Bool, time::Float64, kwargs::@Kwargs{maxiterations::Int64})
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/common_operators/nonlinear_operator.jl:490
 [17] macro expansion
    @ ./timing.jl:645 [inlined]
 [18] assemble_system!(A::ExtendableFEMBase.FEMatrix{Float64, Int64, Float64, Int32, 1, 1, 1}, b::ExtendableFEMBase.FEVector{Float64, Float64, Int32, Vector{Float64}}, sol::ExtendableFEMBase.FEVector{Float64, Float64, Int32, Vector{Float64}}, PD::ExtendableFEM.ProblemDescription, SC::ExtendableFEM.SolverConfiguration{ExtendableFEMBase.FEMatrix{Float64, Int64, Float64, Int32, 1, 1, 1}, ExtendableFEMBase.FEVector{Float64, Float64, Int32, Vector{Float64}}, ExtendableFEMBase.FEVector{Float64, Float64, Int32, Vector{Float64}}}, timer::TimerOutputs.TimerOutput; kwargs::@Kwargs{maxiterations::Int64})
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/solvers.jl:134
 [19] solve(PD::ExtendableFEM.ProblemDescription, FES::ExtendableFEMBase.FESpace{Float64, Int32, H1P1{1}, ON_CELLS}, SC::Nothing; unknowns::Vector{ExtendableFEM.Unknown}, kwargs::@Kwargs{maxiterations::Int64})
    @ ExtendableFEM ~/git-stuff/julia/ExtendableFEM/src/solvers.jl:533
 [20] solve (repeats 2 times)
    @ ~/git-stuff/julia/ExtendableFEM/src/solvers.jl:469 [inlined]
 [21] main()
    @ Main.MWE ~/git-stuff/julia/ExtendableFEM/examples/MWE.jl:49
 [22] top-level scope
    @ REPL[5]:1
 [23] eval(m::Module, e::Any)
    @ Core ./boot.jl:489
 [24] eval
    @ ./Base_compiler.jl:146 [inlined]
 [25] repleval(m::Module, code::Expr, ::String)
    @ VSCodeServer ~/.vscode-oss/extensions/julialang.language-julia-1.149.2-universal/scripts/packages/VSCodeServer/src/repl.jl:229
 [26] #evalrepl##2
    @ ~/.vscode-oss/extensions/julialang.language-julia-1.149.2-universal/scripts/packages/VSCodeServer/src/repl.jl:192 [inlined]
 [27] with_logstate(f::VSCodeServer.var"#evalrepl##2#evalrepl##3"{Module, Expr, REPL.LineEditREPL, REPL.LineEdit.Prompt}, logstate::Base.CoreLogging.LogState)
    @ Base.CoreLogging ./logging/logging.jl:540
 [28] with_logger
    @ ./logging/logging.jl:651 [inlined]
 [29] (::VSCodeServer.var"#evalrepl##0#evalrepl##1"{Module, Expr, REPL.LineEditREPL, REPL.LineEdit.Prompt})()
    @ VSCodeServer ~/.vscode-oss/extensions/julialang.language-julia-1.149.2-universal/scripts/packages/VSCodeServer/src/repl.jl:193
 [30] (::VSCodeServer.var"#start_eval_backend##0#start_eval_backend##1")()
    @ VSCodeServer ~/.vscode-oss/extensions/julialang.language-julia-1.149.2-universal/scripts/packages/VSCodeServer/src/eval.jl:34

@gdalle
Copy link

gdalle commented Aug 14, 2025

@adrhill can you maybe take a look?

@adrhill
Copy link

adrhill commented Aug 14, 2025

The stack trace is odd, ForwardDiff's Dual and SparseConnectivityTracer's GradientTracer should never interact. They are part of two separate passes: first sparsity detection with SCT, then AD with ForwardDiff.
Yet we somehow end up with the following type in L9 of the stacktrace (shortened for legibility):

Dual{
    Tag{
        Main.MWE.var"#nonlinear_DF!##1#nonlinear_DF!##2"{QPInfos{...}},
        GradientTracer{...}
    }, 
    GradientTracer{...}, 
    1
}

Where and how do you call SCT?

@adrhill
Copy link

adrhill commented Aug 14, 2025

In case you are attempting to obtain a second-order sparsity pattern by using SCT to trace through a ForwardDiff computation, you might just want to use SCT's hessian_sparsity function instead.

@gdalle
Copy link

gdalle commented Aug 14, 2025

Unless the two differentiations are with respect to different variables?

@pjaap
Copy link
Member

pjaap commented Aug 14, 2025

Thanks for your input. I think we'll continue the discussion once @chmerdon and me are back from vacation in September 😄

@adrhill
Copy link

adrhill commented Aug 14, 2025

Same here, so that works well for me 😂

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants