From 3b0dd5182429f29a9179e16d8c3fafb0ee159af0 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Mon, 4 Aug 2025 23:16:40 +0900 Subject: [PATCH 1/5] switched from SparseDiffTools to SparseConnectivityTracer and DifferentiationInterface in jacobian preparation in NonlinearOperator --- Project.toml | 10 +++- src/ExtendableFEM.jl | 9 ++- src/common_operators/nonlinear_operator.jl | 68 +++++++++------------- 3 files changed, 41 insertions(+), 46 deletions(-) diff --git a/Project.toml b/Project.toml index 7a7fabc5..b63cc407 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ authors = ["Christian Merdon ", "Patrick Jaap ([], "subset of regions where operator should be assembly only"), :sparse_jacobians => (true, "use sparse jacobians"), :sparse_jacobians_pattern => (nothing, "user provided sparsity pattern for the sparse jacobians (in case automatic detection fails)"), + :autodiff_backend => (AutoForwardDiff(), "backend for automatic differentiation"), :time_dependent => (false, "operator is time-dependent ?"), :verbosity => (0, "verbosity level"), ) @@ -236,31 +235,15 @@ function build_assembler!(A::AbstractMatrix, b::AbstractVector, O::NonlinearOper ## prepare parameters QPj = QPInfos(xgrid; time = time, x = ones(Tv, size(xgrid[Coordinates], 1)), params = O.parameters[:params]) kernel_params = (result, input) -> (O.kernel(result, input, QPj)) - if sparse_jacobians - input_args = zeros(Tv, op_offsets_args[end] + O.parameters[:extra_inputsize]) - result_kernel = zeros(Tv, op_offsets_test[end]) - if sparsity_pattern === nothing - sparsity_pattern = Symbolics.jacobian_sparsity(kernel_params, result_kernel, input_args) - end - jac = Float64.(sparse(sparsity_pattern)) - value = zeros(Tv, op_offsets_test[end]) - colors = matrix_colors(jac) - Dresult = nothing - cfg = ForwardColorJacCache( - kernel_params, input_args, nothing; - dx = nothing, - colorvec = colors, - sparsity = sparsity_pattern - ) - else - input_args = zeros(Tv, op_offsets_args[end] + O.parameters[:extra_inputsize]) - result_kernel = zeros(Tv, op_offsets_test[end]) - Dresult = DiffResults.JacobianResult(result_kernel, input_args) - jac = DiffResults.jacobian(Dresult) - value = DiffResults.value(Dresult) - cfg = ForwardDiff.JacobianConfig(kernel_params, result_kernel, input_args, ForwardDiff.Chunk{op_offsets_args[end]}()) - end - push!(Kj, KernelEvaluator(input_args, QPj, result_kernel, jac, value, Dresult, cfg, kernel_params)) + input_args = zeros(Tv, op_offsets_args[end] + O.parameters[:extra_inputsize]) + result_kernel = zeros(Tv, op_offsets_test[end]) + sparse_forward_backend = AutoSparse( + O.parameters[:autodiff_backend]; + sparsity_detector = TracerSparsityDetector(), + coloring_algorithm = GreedyColoringAlgorithm() + ) + jac_prep = prepare_jacobian(kernel_params, result_kernel, sparse_forward_backend, input_args) + push!(Kj, KernelEvaluator(input_args, QPj, result_kernel, jac_prep, sparse_forward_backend, kernel_params)) end ## prepare parallel assembly @@ -301,11 +284,12 @@ function build_assembler!(A::AbstractMatrix, b::AbstractVector, O::NonlinearOper ## extract kernel properties params = K.params input_args = K.input_args - result_kernel = K.result_kernel - cfg = K.cfg - Dresult = K.Dresult - jac = K.jac - value = K.value + value = K.result_kernel + jac_prep = K.jac_prep + jac_backend = K.jac_backend + # todo: get sparse jacobians to work (need to extract sparsity pattern) + sparse_jacobians = false + jac = zeros(Tv, length(input_args), length(value)) kernel_params = K.kernel params.time = time @@ -359,12 +343,14 @@ function build_assembler!(A::AbstractMatrix, b::AbstractVector, O::NonlinearOper ## get global x for quadrature point eval_trafo!(params.x, L2G, xref[qp]) if use_autodiff - if sparse_jacobians - forwarddiff_color_jacobian!(jac, kernel_params, input_args, cfg) - kernel_params(value, input_args) - else - ForwardDiff.chunk_mode_jacobian!(Dresult, kernel_params, result_kernel, input_args, cfg) - end + DifferentiationInterface.value_and_jacobian!( + kernel_params, + value, + jac, + jac_prep, + jac_backend, + input_args + ) else O.jacobian(jac, input_args, params) O.kernel(value, input_args, params) From e1c9e433e43d2a9471234fb2be35376c965dbc41 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Mon, 4 Aug 2025 23:19:37 +0900 Subject: [PATCH 2/5] fixed wrong size of jacobian --- src/common_operators/nonlinear_operator.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/common_operators/nonlinear_operator.jl b/src/common_operators/nonlinear_operator.jl index dcd56530..90d93ae9 100644 --- a/src/common_operators/nonlinear_operator.jl +++ b/src/common_operators/nonlinear_operator.jl @@ -289,7 +289,7 @@ function build_assembler!(A::AbstractMatrix, b::AbstractVector, O::NonlinearOper jac_backend = K.jac_backend # todo: get sparse jacobians to work (need to extract sparsity pattern) sparse_jacobians = false - jac = zeros(Tv, length(input_args), length(value)) + jac = zeros(Tv, length(value), length(input_args)) kernel_params = K.kernel params.time = time From 6c87755ec852c8673da943dc6a08b36fa3c36531 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Mon, 4 Aug 2025 23:43:12 +0900 Subject: [PATCH 3/5] fixed sparsity pattern issue --- Project.toml | 4 +++- src/ExtendableFEM.jl | 4 ++-- src/common_operators/nonlinear_operator.jl | 20 ++++++++++++++++---- 3 files changed, 21 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index b63cc407..a2c47b56 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "ExtendableFEM" uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed" -version = "1.3.0" authors = ["Christian Merdon ", "Patrick Jaap "] +version = "1.3.0" [deps] +ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" @@ -25,6 +26,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] +ADTypes = "1.16.0" Aqua = "0.8" CommonSolve = "0.2" DiffResults = "1" diff --git a/src/ExtendableFEM.jl b/src/ExtendableFEM.jl index 9e5cafb1..1579918c 100644 --- a/src/ExtendableFEM.jl +++ b/src/ExtendableFEM.jl @@ -68,13 +68,13 @@ using LinearSolve: LinearSolve, LinearProblem, UMFPACKFactorization, deleteat!, using Printf: Printf, @printf, @sprintf using SparseArrays: SparseArrays, AbstractSparseArray, SparseMatrixCSC, findnz, nnz, nzrange, rowvals, sparse - +using ADTypes: ADTypes, KnownJacobianSparsityDetector using SparseConnectivityTracer: SparseConnectivityTracer, TracerSparsityDetector using DifferentiationInterface: DifferentiationInterface, AutoSparse, AutoForwardDiff, prepare_jacobian -using SparseMatrixColorings: GreedyColoringAlgorithm #, sparsity_pattern +using SparseMatrixColorings: GreedyColoringAlgorithm, sparsity_pattern using Symbolics: Symbolics using SciMLBase: SciMLBase using TimerOutputs: TimerOutput, print_timer, @timeit diff --git a/src/common_operators/nonlinear_operator.jl b/src/common_operators/nonlinear_operator.jl index 90d93ae9..8032c85d 100644 --- a/src/common_operators/nonlinear_operator.jl +++ b/src/common_operators/nonlinear_operator.jl @@ -229,7 +229,6 @@ function build_assembler!(A::AbstractMatrix, b::AbstractVector, O::NonlinearOper Kj = Array{KernelEvaluator, 1}([]) sparse_jacobians = O.parameters[:sparse_jacobians] - sparsity_pattern = O.parameters[:sparse_jacobians_pattern] use_autodiff = O.jacobian === nothing for EG in EGs ## prepare parameters @@ -237,9 +236,14 @@ function build_assembler!(A::AbstractMatrix, b::AbstractVector, O::NonlinearOper kernel_params = (result, input) -> (O.kernel(result, input, QPj)) input_args = zeros(Tv, op_offsets_args[end] + O.parameters[:extra_inputsize]) result_kernel = zeros(Tv, op_offsets_test[end]) + if O.parameters[:sparse_jacobians_pattern] === nothing + sparsity_detector = TracerSparsityDetector() + else + sparsity_detector = KnownJacobianSparsityDetector(O.parameters[:sparse_jacobians_pattern]) + end sparse_forward_backend = AutoSparse( O.parameters[:autodiff_backend]; - sparsity_detector = TracerSparsityDetector(), + sparsity_detector = sparsity_detector, coloring_algorithm = GreedyColoringAlgorithm() ) jac_prep = prepare_jacobian(kernel_params, result_kernel, sparse_forward_backend, input_args) @@ -287,9 +291,17 @@ function build_assembler!(A::AbstractMatrix, b::AbstractVector, O::NonlinearOper value = K.result_kernel jac_prep = K.jac_prep jac_backend = K.jac_backend - # todo: get sparse jacobians to work (need to extract sparsity pattern) sparse_jacobians = false - jac = zeros(Tv, length(value), length(input_args)) + 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(sparsity_pattern)) + else + jac = zeros(Tv, length(value), length(input_args)) + end kernel_params = K.kernel params.time = time From c5b24b68c5a344f0f7c700400029c58680ac1322 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Tue, 5 Aug 2025 00:18:07 +0900 Subject: [PATCH 4/5] sparsity patterns for kernels of BilinearOperator now also done by SparseConnectivityTracer, could removed Symbolics as a main dependence now, version bump + changelog --- CHANGELOG.md | 8 ++++++++ Project.toml | 6 +++--- src/ExtendableFEM.jl | 5 ++--- src/common_operators/bilinear_operator.jl | 6 ++++-- src/common_operators/bilinear_operator_dg.jl | 6 ++++-- test/runtests.jl | 1 + 6 files changed, 22 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f45b8c23..715a52cc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,13 @@ # CHANGES +## v1.5.0 + +### Added + - new parameter `autodiff_backend` for NonlinearOperator to change differentiation backend + +### Changed + - sparsity patterns of local jacobians are now resolved by SparseConnectivityTracer + - local jacobians are now prepared by DifferentiationInterface ## v1.4.0 diff --git a/Project.toml b/Project.toml index a2c47b56..b0ce2309 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableFEM" uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed" authors = ["Christian Merdon ", "Patrick Jaap "] -version = "1.3.0" +version = "1.5.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" @@ -21,7 +21,6 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" -Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" @@ -70,9 +69,10 @@ OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea" Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" [targets] -test = ["Aqua", "ExampleJuggler", "ExplicitImports", "IncompleteLU", "Metis", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "SimplexGridFactory", "StaticArrays", "Test", "TetGen", "Triangulate"] +test = ["Aqua", "ExampleJuggler", "ExplicitImports", "IncompleteLU", "Metis", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "SimplexGridFactory", "StaticArrays", "Symbolics", "Test", "TetGen", "Triangulate"] diff --git a/src/ExtendableFEM.jl b/src/ExtendableFEM.jl index 1579918c..ebba4556 100644 --- a/src/ExtendableFEM.jl +++ b/src/ExtendableFEM.jl @@ -66,16 +66,15 @@ using LinearAlgebra: LinearAlgebra, copyto!, isposdef, mul!, norm using LinearSolve: LinearSolve, LinearProblem, UMFPACKFactorization, deleteat!, init, solve using Printf: Printf, @printf, @sprintf -using SparseArrays: SparseArrays, AbstractSparseArray, SparseMatrixCSC, findnz, nnz, +using SparseArrays: SparseArrays, AbstractSparseArray, findnz, nnz, nzrange, rowvals, sparse using ADTypes: ADTypes, KnownJacobianSparsityDetector -using SparseConnectivityTracer: SparseConnectivityTracer, TracerSparsityDetector +using SparseConnectivityTracer: SparseConnectivityTracer, TracerSparsityDetector, jacobian_sparsity using DifferentiationInterface: DifferentiationInterface, AutoSparse, AutoForwardDiff, prepare_jacobian using SparseMatrixColorings: GreedyColoringAlgorithm, sparsity_pattern -using Symbolics: Symbolics using SciMLBase: SciMLBase using TimerOutputs: TimerOutput, print_timer, @timeit using UnicodePlots: UnicodePlots diff --git a/src/common_operators/bilinear_operator.jl b/src/common_operators/bilinear_operator.jl index ba6e4770..64e1e9a2 100644 --- a/src/common_operators/bilinear_operator.jl +++ b/src/common_operators/bilinear_operator.jl @@ -410,7 +410,8 @@ function build_assembler!(A::AbstractMatrix, O::BilinearOperator{Tv}, FE_test, F coupling_matrix::Matrix{Bool} = ones(Bool, nansatz, ntest) if use_sparsity_pattern kernel_params = (result, input) -> (O.kernel(result, input, O.QP_infos[1])) - sparsity_pattern = SparseMatrixCSC{Float64, Int}(Symbolics.jacobian_sparsity(kernel_params, zeros(Tv, op_offsets_test[end]), zeros(Tv, op_offsets_ansatz[end]))) + detector = TracerSparsityDetector() + sparsity_pattern = jacobian_sparsity(kernel_params, zeros(Tv, op_offsets_test[end]), zeros(Tv, op_offsets_ansatz[end]), detector) ## find out which test and ansatz functions couple for id in 1:nansatz @@ -779,7 +780,8 @@ function build_assembler!(A, O::BilinearOperator{Tv}, FE_test, FE_ansatz; time = coupling_matrix::Matrix{Bool} = ones(Bool, nansatz, ntest) if use_sparsity_pattern kernel_params = (result, input) -> (O.kernel(result, input, O.QP_infos[1])) - sparsity_pattern = SparseMatrixCSC{Float64, Int}(Symbolics.jacobian_sparsity(kernel_params, zeros(Tv, op_offsets_test[end]), zeros(Tv, op_offsets_ansatz[end]))) + detector = TracerSparsityDetector() + sparsity_pattern = jacobian_sparsity(kernel_params, zeros(Tv, op_offsets_test[end]), zeros(Tv, op_offsets_ansatz[end]), detector) ## find out which test and ansatz functions couple for id in 1:nansatz diff --git a/src/common_operators/bilinear_operator_dg.jl b/src/common_operators/bilinear_operator_dg.jl index a4235cf8..a532fbcc 100644 --- a/src/common_operators/bilinear_operator_dg.jl +++ b/src/common_operators/bilinear_operator_dg.jl @@ -367,7 +367,8 @@ function build_assembler!(A, O::BilinearOperatorDG{Tv}, FE_test, FE_ansatz, FE_a coupling_matrix::Matrix{Bool} = ones(Bool, nansatz, ntest) if use_sparsity_pattern kernel_params = (result, input) -> (O.kernel(result, input, O.QP_infos[1])) - sparsity_pattern = SparseMatrixCSC{Float64, Int}(Symbolics.jacobian_sparsity(kernel_params, zeros(Tv, op_offsets_test[end]), zeros(Tv, op_offsets_ansatz[end]))) + detector = TracerSparsityDetector() + sparsity_pattern = jacobian_sparsity(kernel_params, zeros(Tv, op_offsets_test[end]), zeros(Tv, op_offsets_ansatz[end]), detector) ## find out which test and ansatz functions couple for id in 1:nansatz @@ -782,7 +783,8 @@ function build_assembler!(A, O::BilinearOperatorDG{Tv}, FE_test, FE_ansatz; time coupling_matrix::Matrix{Bool} = ones(Bool, nansatz, ntest) if use_sparsity_pattern kernel_params = (result, input) -> (O.kernel(result, input, O.QP_infos[1])) - sparsity_pattern = SparseMatrixCSC{Float64, Int}(Symbolics.jacobian_sparsity(kernel_params, zeros(Tv, op_offsets_test[end]), zeros(Tv, op_offsets_ansatz[end]))) + detector = TracerSparsityDetector() + sparsity_pattern = jacobian_sparsity(kernel_params, zeros(Tv, op_offsets_test[end]), zeros(Tv, op_offsets_ansatz[end]), detector) ## find out which test and ansatz functions couple for id in 1:nansatz diff --git a/test/runtests.jl b/test/runtests.jl index 0954d6d3..6093fdf3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ using Metis using Aqua using Triangulate using SimplexGridFactory +using Symbolics include("test_dgblf.jl") From cfe4d7e32581bc0eecaf735d2c2e4ea228268b77 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Thu, 7 Aug 2025 10:22:49 +0900 Subject: [PATCH 5/5] fixed sparsity_pattern in NonlinearOperator, added news packages to dependencies section in doc --- docs/src/index.md | 3 ++- src/ExtendableFEM.jl | 4 +--- src/common_operators/nonlinear_operator.jl | 9 ++++----- 3 files changed, 7 insertions(+), 9 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index ccea96fd..033aa4b8 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -23,7 +23,8 @@ - [GridVisualize.jl](https://github.com/WIAS-PDELib/GridVisualize.jl) - [DocStringExtensions.jl](https://github.com/JuliaDocs/DocStringExtensions.jl) - [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) -- [DiffResults.jl](https://github.com/JuliaDiff/DiffResults.jl) +- [DifferentiationInterface.jl](https://github.com/JuliaDiff/DifferentiationInterface.jl) +- [SparseConnectivityTracer.jl](https://github.com/adrhill/SparseConnectivityTracer.jl) ## Getting Started diff --git a/src/ExtendableFEM.jl b/src/ExtendableFEM.jl index 9e5459da..11147ba8 100644 --- a/src/ExtendableFEM.jl +++ b/src/ExtendableFEM.jl @@ -71,9 +71,7 @@ using SparseArrays: SparseArrays, AbstractSparseArray, findnz, nnz, using ADTypes: ADTypes, KnownJacobianSparsityDetector using SparseConnectivityTracer: SparseConnectivityTracer, TracerSparsityDetector, jacobian_sparsity using DifferentiationInterface: DifferentiationInterface, - AutoSparse, - AutoForwardDiff, - prepare_jacobian + AutoSparse, AutoForwardDiff, prepare_jacobian using SparseMatrixColorings: GreedyColoringAlgorithm, sparsity_pattern using SciMLBase: SciMLBase using TimerOutputs: TimerOutput, print_timer, @timeit diff --git a/src/common_operators/nonlinear_operator.jl b/src/common_operators/nonlinear_operator.jl index 8032c85d..8ec0a445 100644 --- a/src/common_operators/nonlinear_operator.jl +++ b/src/common_operators/nonlinear_operator.jl @@ -291,18 +291,17 @@ function build_assembler!(A::AbstractMatrix, b::AbstractVector, O::NonlinearOper value = K.result_kernel jac_prep = K.jac_prep jac_backend = K.jac_backend - sparse_jacobians = false - if sparse_jacobians + kernel_params = K.kernel + if O.parameters[:sparse_jacobians] if O.parameters[:sparse_jacobians_pattern] === nothing - jac_sparsity_pattern = DifferentiationInterface.sparsity_pattern(jac_prep) + jac_sparsity_pattern = sparsity_pattern(jac_prep) else jac_sparsity_pattern = O.parameters[:sparse_jacobians_pattern] end - jac = Tv.(sparse(sparsity_pattern)) + jac = Tv.(sparse(jac_sparsity_pattern)) else jac = zeros(Tv, length(value), length(input_args)) end - kernel_params = K.kernel params.time = time ndofs_test::Array{Int, 1} = [get_ndofs(ON_CELLS, FE, EG) for FE in FETypes_test]