Modularize assembly internals and add SparseMatrixCSR extension

[termi-official]

Towards #848 and #628 . This should show how we could assemble into custom sparse matrix formats.

TODO

• Fix inhomogeneous boundary conditions
• Add devdocs
• Add test

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 98.42520% with 2 lines in your changes missing coverage. Please review.

Project coverage is 94.11%. Comparing base (c718c73) to head (0cce299).
Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/Dofs/ConstraintHandler.jl	97.50%	1 Missing ⚠️
src/assembler.jl	90.90%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #864      +/-   ##
==========================================
+ Coverage   94.03%   94.11%   +0.07%     
==========================================
  Files          39       40       +1     
  Lines        6555     6643      +88     
==========================================
+ Hits         6164     6252      +88     
  Misses        391      391              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[KnutAM]

Just took a look as I was curious, nice work on factoring out some parts - this will make it easier to extend!

[termi-official]

@fredrikekre Do you have some time to review?

[termi-official]

Thanks @KnutAM for the pointer on the canonical keyword. PR should be good to go now.

[KnutAM]

Really nice to make this more flexible!

Do I understand this PR correctly that the main purpose is to introduce the necessary abstractions to allow extensions to different (CSR) matrix formats, and that the SparseMatrixCSR extension is just an example?

While testing to understand, I realized the following behavior,

julia> K = sparsecsr([2], [2], [2.])
2×2 SparseMatrixCSR{1, Float64, Int64} with 1 stored entry:
  [2, 2]  =  2.0

julia> K+K
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  4.0

Which could lead to quite severe user issues if used in e.g. the transient heat tutorial here. So perhaps, instead of adding an extension that we have to maintain, could we add the implementation to the test suite instead?

And when JuliaSparse/SparseArrays.jl#546 is finalized, we could directly support most ops for abstract CSR inside Ferrite (but unless SparseMatricesCSR adopts the interface, we probably don't want to keep supporting that)?

[KnutAM]

Since this only modifies f.

Suggested change

	function add_inhomogeneities!(K, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping)
	function add_inhomogeneities!(f::AbstractVector, K::AbstractMatrix, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping)

I would prefer the interface
add_inhomogeneities!(f::AbstractVector, ch::ConstraintHandler, K::AbstractMatrix) for extra clarity (that only f is modified), and easier extensibility (following the change for zero_out_columns! and zero_out_rows!).

[termi-official]

Good catch. That makes sense from an interface point of view.

[KnutAM]

Out of curiosity, any reason not passing the constraint handler instead of inhomogeneities, prescribed_dofs, and dofmapping? (Esp. since dofmapping isn't always used).

[termi-official]

No specific reason. The internals are not consistent here. I changed it to take the ch for now to simplify forward compatibility.

[termi-official]

Thanks for the review Knut!

Do I understand this PR correctly that the main purpose is to introduce the necessary abstractions to allow extensions to different (CSR) matrix formats, and that the SparseMatrixCSR extension is just an example?

Almost. It should contain all changes to allow users to add support for any of their favorite serial matrix formats in their packages (BSR,CSR,CSC,...) and it adds helpers for CSR/CSC to do some of the stuff. The main blocker on master is right now that there is no easy way to add support for new matrix formats in the constraint handler and this PR introduces necessary hooks. I am not sure if I missed somthing more advanced than Dirichlet tho.

So perhaps, instead of adding an extension that we have to maintain, could we add the implementation to the test suite instead?

Can elaborate why adding an extension should make us responsible for fixing bugs in upstream packages? If I would not add the CSR stuff here then there would be again the question why the refactor which I propose is necessary in first place...

And when JuliaSparse/SparseArrays.jl#546 is finalized, we could directly support most ops for abstract CSR inside Ferrite (but unless SparseMatricesCSR adopts the interface, we probably don't want to keep supporting that)?

I have essentially given up on the effort of establishing an interface for different matrix types centrally right now, as it does not seem to be something people actually want, if you look at all the discussion in e.g. slack or discourse. But to answer the question: Yes, that was my initial plan, to have interfaces to different matrix formats which we can adopt. I just started adding a custom intermediate layer in Thunderbolt to circumvent the lack of a well-defined interface.

[KnutAM]

Can elaborate why adding an extension should make us responsible for fixing bugs in upstream packages?

No, I meant that it makes us responsible for maintaining the extension. And IIUC we don't actually have a use case to support SparseMatricesCSR.SparseMatrixCSR specifically (since AFAIU there are no solvers for that matrix type), it is an extra maintenance burden to support that one. Isn't our goal to support (or allow extensibility to support) formats with solvers, e.g. Ginkgo or Hypre?

Another option could be to add a how-to on how to support a new matrix format, if the interface is considered stable enough? (Or just mark as experimental/devdocs)

[termi-official]

No, I meant that it makes us responsible for maintaining the extension

Right, as with all other extensions.

AFAIU there are no solvers for that matrix type

IterativeSolvers.jl, Krylov.jl, KrylovKit.jl to name a few and we can implement how to cheaply convert the format to most direct solvers like Pardiso.jl , which require CSR. I also use it in Thunderbolt to wrap for my parallel CSR format and for parallel preconditioners therein.

Isn't our goal to support (or allow extensibility to support) formats with solvers, e.g. Ginkgo or Hypre?

Yes. Should I add these also to the PR?

Another option could be to add a how-to on how to support a new matrix format, if the interface is considered stable enough?

I wanted to do this later after PRs to support the big packages are done, so we actually can figure out if it is stable enough.

[KnutAM]

IterativeSolvers.jl, Krylov.jl, KrylovKit.jl to name a few

Aha, then I misunderstood the conversations here, #864 (comment), #864 (comment), #864 (comment), #864 (comment).

Should I add these also to the PR?

No, that wasn't what I meant. Related to above, I just didn't know that SparseMatricesCSR.jl had that wide solver support as you say.

[termi-official]

Probably to elaborate here, most of the iterative solvers packages literally just need to define matrix-vector products to be defined for the thing you shove into them.

[KnutAM]

Some details, but otherwise LGTM.

[termi-official]

Thanks for the detailed review Knut!

[KnutAM]

Just some formatting (and we should run CI as it includes the formatting check) and the changelog is included in the docs.