Factorize matrix for SvaerdKalischEquations1D

[JoshuaLampert]

Closes #107. Let's see if CI is happy with the cholesky decomposition this time. It is significantly faster and less allocating than lu. Since we know that hmD1betaD1 is sparse since D1betaD1 is constructed to be sparse, I only implemented the first branch in your proposal #107 (comment), @ranocha.

Result of julia> include("examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl"); for main:

──────────────────────────────────────────────────────────────────────────────
        DispersiveSWE                 Time                    Allocations      
                             ───────────────────────   ────────────────────────
      Tot / % measured:           2.34s /  99.4%           3.61GiB /  99.8%    

 Section             ncalls     time    %tot     avg     alloc    %tot      avg
 ──────────────────────────────────────────────────────────────────────────────
 rhs!                 5.16k    2.31s   99.2%   448μs   3.59GiB   99.7%   730KiB
   dv elliptic        5.16k    1.38s   59.4%   268μs   2.77GiB   76.9%   563KiB
   dv hyperbolic      5.16k    865ms   37.1%   168μs    590MiB   16.0%   117KiB
   deta hyperbolic    5.16k   54.8ms    2.3%  10.6μs    250MiB    6.8%  49.5KiB
   ~rhs!~             5.16k   8.03ms    0.3%  1.55μs    816KiB    0.0%     162B
   source terms       5.16k    381μs    0.0%  73.7ns     0.00B    0.0%    0.00B
 analyze solution        86   17.6ms    0.8%   205μs   10.0MiB    0.3%   119KiB
 ──────────────────────────────────────────────────────────────────────────────

and this branch (note that I added some @timeit statements to have a better overview):

────────────────────────────────────────────────────────────────────────────────────
           DispersiveSWE                    Time                    Allocations      
                                   ───────────────────────   ────────────────────────
         Tot / % measured:              1.61s /  99.2%           2.43GiB /  99.6%    

 Section                   ncalls     time    %tot     avg     alloc    %tot      avg
 ────────────────────────────────────────────────────────────────────────────────────
 rhs!                       5.16k    1.58s   99.0%   306μs   2.41GiB   99.6%   489KiB
   dv hyperbolic            5.16k    873ms   54.7%   169μs    590MiB   23.8%   117KiB
   dv elliptic              5.16k    641ms   40.1%   124μs   1.59GiB   65.6%   322KiB
     cholesky!              5.16k    478ms   29.9%  92.6μs   0.97GiB   40.1%   197KiB
     solve linear system    5.16k   82.3ms    5.2%  15.9μs    144MiB    5.8%  28.6KiB
     set hmD1betaD1         5.16k   75.5ms    4.7%  14.6μs    488MiB   19.7%  96.9KiB
     set tmp1               5.16k   2.49ms    0.2%   483ns     0.00B    0.0%    0.00B
     ~dv elliptic~          5.16k   2.35ms    0.1%   456ns   2.94KiB    0.0%    0.58B
   deta hyperbolic          5.16k   60.7ms    3.8%  11.7μs    250MiB   10.1%  49.5KiB
   ~rhs!~                   5.16k   5.47ms    0.3%  1.06μs    810KiB    0.0%     161B
   source terms             5.16k    403μs    0.0%  78.1ns     0.00B    0.0%    0.00B
 analyze solution              86   15.8ms    1.0%   184μs   10.0MiB    0.4%   119KiB
 ────────────────────────────────────────────────────────────────────────────────────

cholesky! still allocates more than I would have hoped, but I think there is not much we can do against, right? But overall this is a nice speedup.

[JoshuaLampert]

There are some quite significant changes in the reference values in the CG case...

[ranocha]

Thanks!

[ranocha]

It looks like some reference values changed

[JoshuaLampert]

It looks like some reference values changed

It looks like they changed a bit too much. This, e.g., doesn't look good 😬

[JoshuaLampert]

It looks like they changed a bit too much. This, e.g., doesn't look good 😬

Ah, the system matrix is not symmetric in the CG case.

[ranocha]

Right, something like

M * (Diagonal(h) - D1betaD1)

should be symmetric (in the usual sense) since the operator we wrap right now is symmetric with respect to the mass matrix M, isn't it?

[JoshuaLampert]

Yes, exactly. $M(\textrm{diag}(h) - D_1\textrm{diag}(\beta) D_1)$ is guaranteed to by symmetric and also positive definite (since $M$ positive definite, $h > 0$, and $\beta > 0$ due to $\tilde\beta > 0$ and $D > 0$):

$$ x^T(M\textrm{diag}(h) - MD_1\textrm{diag}(\beta) D_1)x = \underbrace{x^TM\textrm{diag}(h)x}_{> 0} + \underbrace{(D_1x)^TM\textrm{diag}(\beta) D_1x}_{\ge 0}>0 $$

Thus, it makes more sense to always apply the Cholesky decomposition to $M(\textrm{diag}(h) - D_1\textrm{diag}(\beta) D_1)$. Thanks for the hint! I've pushed an adapted version in 0643d81. Let's hope CI passes now.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 96.83%. Comparing base (bda70af) to head (0643d81).

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #114      +/-   ##
==========================================
+ Coverage   96.82%   96.83%   +0.01%     
==========================================
  Files          17       17              
  Lines        1227     1234       +7     
==========================================
+ Hits         1188     1195       +7     
  Misses         39       39              

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

[ranocha]

Thanks!