Add cache-based solve test for multiple right-hand sides

[Da-Be-Ru]

For my application, I need to extract the state matrix for a linear convection-diffusion system, factorize it, and solve the system for multiple right-hand sides.

It'd be nice to have a test case for this that checks for regressions in the interface. It seems to me that such a regression has been introduced some time after the move to VoronoiFVM v2 as it used to be that I could simply feed the ExtendableSparseMatrix to the LinearSolve package to repeatedly solve! with updated right-hand sides in the cache.
This seems to not be possible anymore as of the latest version of LinearSolve v3.7.0 - is that intended?

Furthermore, the test environment of VoronoiFVM is again somehow broken and outdated, partly due to a related issue for ExtendableSparse and also because LinearSolve is fixed to v2.39.x by some compatibility constraints I have been unable to figure out so far.

Thus, to reproduce the regression in this test, one has to run it within an environment that includes LinearSolve v3.7.0 which will produce a similar error as in the related issue from ExtendableSparse.

I could also live with having to extract the sparse matrix directly. It'd just be nice to settle on some interface and have a test included to guard against downstream breakage.

[j-fu]

I think we need to invent a good public API call signature for this. eval_and_assemble is not the right one.
EDIT:

Could this be sometinh like this ?

cache=linearize(state, u)

or

cache=linearize(system, u)

And what tools do you need for assembling your rhs ?

[Da-Be-Ru]

cache=linearize(system, u)

That'd already be quite helpful.

Currently, what I settled on is to define a system that discretizes a convection-diffusion operator for some solution $Y_0$. My basis elements $Y_i$ satisfy the same operator except that they have an additional Neumann condition on one single boundary segment associated with a respective boundary node.
So I extract the state matrix $A$ ("the linearization") associated with $Y_0$ and assemble the appropriate right-hand side $F_i$ such that $A \cdot Y_i = F_i$ in some separate loop like this (see also your merge request in my project which has the latest working state in compute_basis_elements!):

    F=unknowns(rbssys.basis_system;inival=0.)

    state=VoronoiFVM.SystemState(rbssys.basis_system)
    VoronoiFVM.evaluate_residual_and_jacobian!(state, unknowns(rbssys.basis_system, inival=0))
    
    mtx = SparseMatrixCSC(state.matrix)
    p = LinearProblem(mtx, VoronoiFVM.dofs(F))
    linear_cache = init(p,LinearSolve.UMFPACKFactorization())
    
    for (i, ireg) in enumerate(rbssys.breaction_face_regions)
        ibasis=basis[i]
        for (k,inode) in enumerate(Nodes[i])
            # extract bfacenodefactors for inode 
            # to write into RHS

            # 1. get bfaces for inode
            faceindices=findall(bfacenodes.==inode)
            # 2. calculate contributions of bfaces to inode
            γ=0.
            for faceind in faceindices
                # if the bface is in ireg
                # (violated e.g. in inodes on very left or right of regions)
                if grid[BFaceRegions][faceind[2]] == ireg
                    γ+=bfacenodefactors[faceind[1],faceind[2]]
                end
            end
            # 3. iterate through reactions and build a basis element per reaction
            for jreac in 1:nreac
                for ispec in 1:nspecies
                    idof=dof(F,ispec,inode)
                    if chemtest
                        F[idof]+=data.S[jreac,ispec]*data.molar_weights[ispec]*γ
                    else
                        F[idof]+=γ
                    end
                end
                # reuse LU factorization
                linear_cache.b=VoronoiFVM.dofs(F)
                sol=LinearSolve.solve!(linear_cache)
                VoronoiFVM.dofs(ibasis[k,jreac]).=sol.u
                F.=0
            end         
        end
    end

Since the assembly machinery of VoronoiFVM always assembles $A$ as well as the residual simultaneously, this seemed to run through much faster (and also clearer) than re-defining the system for each $Y_i$ and letting VoronoiFVM take care of assembling F_i. But this setup is not a hill I'd die on - if we could e.g. re-arrange the assembly such that I could trigger the assembly of only the boundary nodes into only the residual to construct $F_i$, that'd also be adequate.

[j-fu]

Ah we already have eval_residual_and_jacobian in the API, see
https://wias-pdelib.github.io/VoronoiFVM.jl/stable/solver/#Matrix-extraction

This indeed lacks a test. May be you can update the test to this API ?
Also it seems the dofs method is not exported. I would add this to this PR, probably just as public.
Are there any other functions from internal API (not exported/public) which you use ?

[j-fu]

Also see #187: I want to downgrade dof from export to public in 3.0.

[Da-Be-Ru]

eval_residual_and_jacobian contained a bug I took the liberty to fix.

I also removed the _initialize_dirichlet! call in there since that modifies the argument u causing incorrect entries in the residual (we want both the matrix and residual evaluated exactly in the vector u passed by the user).

[Da-Be-Ru]

I also incorporated your suggestion to make dofs public.

And I took the liberty to bump the gitleaks version in the pre-commit hooks due to an upstream-induced panic on my OS. Wasn't sure if I should open up a separate PR for that or not. Feel free to re-arrange ofc.

[Da-Be-Ru]

Concerning the constantly failing tests, it looks to me that NonlinearSolve is pinned to some older and faulty version which is limited by ModelingToolkit which should be upgradable, but simply won't update.

Similarly, LinearSolve is pinned to an older version by AlgebraicMultigrid which should also be upgradable, but won't update.

I think we should open a separate PR to fix the testing environment.

[j-fu]

AMG seems to be held back by ExtendableSparse which gets a new version soon...

As for ModelingToolkit I'll have a look...

Could you design your project code in a branch of your project such that it works an we can see that the changes in this PR are sufficient for your project ?

[Da-Be-Ru]

Could you design your project code in a branch of your project

That's already available. I rebased the branch voronoifvm-v2 from your old MR which works perfectly fine in its current state with the latest VoronoiFVM release.

The changes from this PR probably won't have much of an impact, but I'll check the downstream ramifications nonetheless. The only thing that'd be useful would be some assembly method for only the residuals without going through the effort of computing the linearization $A$. But things are still fine the way they currently are.

AMG seems to be held back by ExtendableSparse

Oh, yes, that definitely seems correct. But then I'm confused as to why this is not actually displayed by checking the test environment with the status interface:

(test) pkg> st --outdated -m
Status `.../VoronoiFVM.jl/test/Manifest.toml`
⌃ [2169fc97] AlgebraicMultigrid v0.6.0 (<v1.0.0)
[...]
⌃ [961ee093] ModelingToolkit v9.32.0 (<v9.69.0)

I would have expected here an output like ⌅ [2169fc97] AlgebraicMultigrid v0.6.0 (<v1.0.0): ExtendableSparse and a similar listing of the culprit for ModelingToolkit. Instead, they're displayed as upgradable without restriction, but clearly they are restricted. This seems very odd.

[j-fu]

I think after Chris' intervention we also should move test dependencies to [extras] @pjaap .
It seems that this way, we also will be forced to add compat entries for test deps. I think this should be another PR.

[Da-Be-Ru]

I incorporated some changes suggested by your review.
And also

It seems that this way, we also will be forced to add compat entries for test deps. I think this should be another PR.

I opened #188 to address this, but wasn't sure which [compats] to put in there so we should discuss this there.

I'd suggest to rebase this branch onto the version of #188 once we're done with that one to get a CI check for the changes here.

[Da-Be-Ru]

Tests seem fine now. Since we patched a buggy function, I'd propose a patch version bump and release once someone has had a look at the current state.

[j-fu]

looks good to me