Ways to calculate sparse Jacobians
by Jürgen Fuhrmann, 2022-10-21
package | version | |
---|---|---|
Any | Any | |
1 | VoronoiFVM | "0.18.3" |
2 | Symbolics | "4.13.0" |
3 | ForwardDiff | "0.10.32" |
4 | SparseDiffTools | "1.27.0" |
We investigate different ways to assemble sparse Jacobians for a finite difference operator $A_h(u)$ discretizing the differential operator
$$ A(u)=-\Delta u^2 + u^2 -1 $$
in $\Omega=(0,1)$ with homogeneous Neumann boundary conditions on an equidistant grid of $n$ points. We only discuss approaches which can be generalized to higher space dimensions and unstructured grids, so e.g. banded and tridiagonal matrix structures are not discussed.
This is the corresponding nonlinear finite difference operator:
Number of discretization points for initial testing:
The methods
Straightforward Jacobian calculation
Just take the operator and calculate the jacobian using ForwardDiff.jl and use a sparse matrix to collect the derivatives. Here we use the ExtendableSparse.jl package in order to avoid cumbersome handling of intermediate coordinate format data.
As we will see in the timing test, the complexity of this operation is $O(n^2)$. Without further structural information, automatic differencing needs to investigate the partial derivatives with respect to all unknowns.
5×5 SparseMatrixCSC{Float64, Int64} with 13 stored entries: 8.25 -8.0 ⋅ ⋅ ⋅ -8.0 16.5 -8.0 ⋅ ⋅ ⋅ -8.0 16.5 -8.0 ⋅ ⋅ ⋅ -8.0 16.5 -8.0 ⋅ ⋅ ⋅ -8.0 8.25
1.6042e-5
Symbolics.jl for sparsity detection
While in the particular case the sparsity pattern of the operator is known – it is tridiagonal – we assume it is not known and we use Symbolics.jl to detect the sparsity pattern and convey the sparsity information to ForwardDiff.jl.
This approach uses the multiple dispatch feature of Julia and calls A_h!
with symbolic data, resulting in a symbolic representation of the operator which can be investigated for sparsity in $O(n)$ time. For a description of this approach see this paper.
In addition, for calculating the Jacobian efficiently it is useful to provide a matrix coloring obtained using SparseDiffTools.jl.
5×5 SparseMatrixCSC{Float64, Int64} with 13 stored entries: 8.25 -8.0 ⋅ ⋅ ⋅ -8.0 16.5 -8.0 ⋅ ⋅ ⋅ -8.0 16.5 -8.0 ⋅ ⋅ ⋅ -8.0 16.5 -8.0 ⋅ ⋅ ⋅ -8.0 8.25
0.000336296
3.0492e-5
SparsityDetection.jl for sparsity detection
Before the advent of Symbolics.jl, SparsityDetection.jl (now deprecated) provided a similar method for sparsity detection of nonlinear operators in $O(n)$ time which we test here in order to be comprehensive.
5×5 SparseMatrixCSC{Float64, Int64} with 13 stored entries: 8.25 -8.0 ⋅ ⋅ ⋅ -8.0 16.5 -8.0 ⋅ ⋅ ⋅ -8.0 16.5 -8.0 ⋅ ⋅ ⋅ -8.0 16.5 -8.0 ⋅ ⋅ ⋅ -8.0 8.25
0.000830708
3.5003e-5
Explored path: SparsityDetection.Path(Bool[], 1)
Assembly from local jacobians
Usually, for finite difference, finite volume and finite element methods, the sparsity structure is defined by the discretization grid topology, and the operator is essentially a superposition of local contributions translated in space according to the grid structure. So we can assemble the global Jacobi matrix from local contributions, very much like in classical finite element assembly. The local Jacobi matrix can be obtained without the need to detect its sparsity. However, for larger coupled systems of PDEs this possibility appears to be worth to be investigated.
For this purpose, we define the local contributions to the nonlinear difference operator as mutating Julia functions, ready to be generalized for systems of PDEs.