improve performance with sparse matrices

[ranocha]

Closes #98

[codecov-commenter]

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

Attention: Patch coverage is 96.85864% with 6 lines in your changes missing coverage. Please review.

Project coverage is 98.09%. Comparing base (4fdc654) to head (1d157ca).
Report is 1 commits behind head on main.

Files	Patch %	Lines
src/mprk.jl	95.20%	6 Missing ⚠️

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #101      +/-   ##
==========================================
- Coverage   98.32%   98.09%   -0.23%     
==========================================
  Files           6        6              
  Lines        1250     1415     +165     
==========================================
+ Hits         1229     1388     +159     
- Misses         21       27       +6     

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

[ranocha]

There seem to be some issues with floating point rounding errors on macOS CI. The other tests seem to be fine

[ranocha]

@SKopecz Do you have an idea why MPRK43II may be so sensitive here - or how we should fix it?

[SKopecz]

@SKopecz Do you have an idea why MPRK43II may be so sensitive here - or how we should fix it?

If possible, I'd exclude MPRK43II from this specific test on macOS for now and create an issue to deal with this later.

I think we need much more information to see what the actual problem is. It looks to me as if the result of linear solve has a negative element.

[ranocha]

@SKopecz Do you have an idea why MPRK43II may be so sensitive here - or how we should fix it?

If possible, I'd exclude MPRK43II from this specific test on macOS for now and create an issue to deal with this later.

👍

I think we need much more information to see what the actual problem is. It looks to me as if the result of linear solve has a negative element.

Yes, that's what I think, too - computing the powers like σ^(1 - q1) * u^q1 seems to be the problem when σ becomes slightly negative

[ranocha]

This PR is ready for a review. Let's keep fingers crossed that CI passes this time 🤞

[SKopecz]

Are you saying that this is the reason for the large number of allocations we see in the linear advection tutorial?

There's also something else we should look at. When we initialize the linear solvers we use something like

linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
                        assumptions = LinearSolve.OperatorAssumptions(true))

The statement alias_A=true allows the linear solver to reuse the memory of the system matrix and modify it. The idea behind this was that in an MPRK scheme we no longer need the system matrix once the system is solved, so the linear solver can change it at will.

For a sparse matrix this modification of the system matrix could become problematic, since the sparsity pattern of L and U does not match that of A, especially in the linear advection tutorial, where A is a banded matrix with an entry in the upper right corner, for which in genaral L or U could become a full triangular matrix.

All I want to say is that we should try alias_A=false and check if we see any difference.

[ranocha]

Shall we use alias_A = !issparse(the_matrix_we_put_in) for the general case?

[SKopecz]

Looks good to me, thanks! I just had a minor comment on the FAQ and the concern regarding alias_A=true could also be moved to an issue.