(Linear) Matrix operator for efficient IMEX

[DanielDoehring]

Motivation: Tell the IMEX schemes the diffusion operator is linear.
The recommended way of doing so is to construct a Splitfunction with a MatrixOperator (this is somewhat undocumented, I took this from NumericalMathematics/DispersiveShallowWater.jl#207 (comment)) and then a SplitProblem

In this way, we can actually use the fact that the diffusion operator is linear (for advection-diffusion only).

Runtimes:

IMEX with above procedure:

─────────────────────────────────────────────────────────────────────────────────
           Trixi.jl                     Time                    Allocations      
                               ───────────────────────   ────────────────────────
       Tot / % measured:            105ms /  11.3%           3.54MiB /   1.7%    

Section                ncalls     time    %tot     avg     alloc    %tot      avg
─────────────────────────────────────────────────────────────────────────────────
rhs!                      301   10.9ms   91.3%  36.2μs   9.33KiB   15.6%    31.7B
  volume integral         301   6.99ms   58.6%  23.2μs     0.00B    0.0%    0.00B
  interface flux          301   1.66ms   13.9%  5.52μs     0.00B    0.0%    0.00B
  prolong2interfaces      301    842μs    7.1%  2.80μs     0.00B    0.0%    0.00B
  surface integral        301    758μs    6.3%  2.52μs     0.00B    0.0%    0.00B
  ~rhs!~                  301    233μs    1.9%   773ns   9.33KiB   15.6%    31.7B
  reset ∂u/∂t             301    216μs    1.8%   719ns     0.00B    0.0%    0.00B
  Jacobian                301    152μs    1.3%   504ns     0.00B    0.0%    0.00B
  prolong2mortars         301   16.7μs    0.1%  55.4ns     0.00B    0.0%    0.00B
  prolong2boundaries      301   12.6μs    0.1%  41.9ns     0.00B    0.0%    0.00B
  mortar flux             301   7.85μs    0.1%  26.1ns     0.00B    0.0%    0.00B
  boundary flux           301   5.05μs    0.0%  16.8ns     0.00B    0.0%    0.00B
  source terms            301   4.82μs    0.0%  16.0ns     0.00B    0.0%    0.00B
analyze solution            2   1.04ms    8.7%   518μs   50.6KiB   84.4%  25.3KiB
─────────────────────────────────────────────────────────────────────────────────

Note that rhs_parabolic is never called!

Explicit method

───────────────────────────────────────────────────────────────────────────────────────
              Trixi.jl                        Time                    Allocations      
                                     ───────────────────────   ────────────────────────
          Tot / % measured:               685ms /  98.3%           99.2KiB /  80.7%    

Section                      ncalls     time    %tot     avg     alloc    %tot      avg
───────────────────────────────────────────────────────────────────────────────────────
parabolic rhs!                4.50k    506ms   75.1%   112μs   18.1KiB   22.7%    4.13B
  calculate gradient          4.50k    233ms   34.5%  51.7μs   7.34KiB    9.2%    1.67B
    volume integral           4.50k    146ms   21.7%  32.5μs     0.00B    0.0%    0.00B
    Jacobian                  4.50k   28.8ms    4.3%  6.39μs     0.00B    0.0%    0.00B
    surface integral          4.50k   17.6ms    2.6%  3.92μs     0.00B    0.0%    0.00B
    interface flux            4.50k   16.4ms    2.4%  3.65μs     0.00B    0.0%    0.00B
    prolong2interfaces        4.50k   14.2ms    2.1%  3.15μs     0.00B    0.0%    0.00B
    reset gradients           4.50k   6.02ms    0.9%  1.34μs     0.00B    0.0%    0.00B
    ~calculate gradient~      4.50k   2.56ms    0.4%   568ns   7.34KiB    9.2%    1.67B
    prolong2boundaries        4.50k    175μs    0.0%  39.0ns     0.00B    0.0%    0.00B
    prolong2mortars           4.50k    146μs    0.0%  32.4ns     0.00B    0.0%    0.00B
    mortar flux               4.50k    102μs    0.0%  22.7ns     0.00B    0.0%    0.00B
    boundary flux             4.50k    102μs    0.0%  22.6ns     0.00B    0.0%    0.00B
  volume integral             4.50k    129ms   19.1%  28.6μs     0.00B    0.0%    0.00B
  calculate viscous fluxes    4.50k   80.1ms   11.9%  17.8μs     0.00B    0.0%    0.00B
  interface flux              4.50k   16.6ms    2.5%  3.68μs     0.00B    0.0%    0.00B
  prolong2interfaces          4.50k   15.0ms    2.2%  3.33μs     0.00B    0.0%    0.00B
  transform variables         4.50k   14.7ms    2.2%  3.26μs     0.00B    0.0%    0.00B
  surface integral            4.50k   10.9ms    1.6%  2.42μs     0.00B    0.0%    0.00B
  ~parabolic rhs!~            4.50k   3.32ms    0.5%   737ns   10.8KiB   13.5%    2.46B
  reset ∂u/∂t                 4.50k   2.05ms    0.3%   456ns     0.00B    0.0%    0.00B
  Jacobian                    4.50k   1.74ms    0.3%   387ns     0.00B    0.0%    0.00B
  prolong2mortars             4.50k    146μs    0.0%  32.5ns     0.00B    0.0%    0.00B
  prolong2boundaries          4.50k    141μs    0.0%  31.3ns     0.00B    0.0%    0.00B
  boundary flux               4.50k    110μs    0.0%  24.5ns     0.00B    0.0%    0.00B
  mortar flux                 4.50k   98.7μs    0.0%  21.9ns     0.00B    0.0%    0.00B
rhs!                          4.50k    167ms   24.7%  37.0μs   9.33KiB   11.7%    2.12B
  volume integral             4.50k    112ms   16.7%  25.0μs     0.00B    0.0%    0.00B
  interface flux              4.50k   24.9ms    3.7%  5.54μs     0.00B    0.0%    0.00B
  prolong2interfaces          4.50k   11.5ms    1.7%  2.55μs     0.00B    0.0%    0.00B
  surface integral            4.50k   10.9ms    1.6%  2.42μs     0.00B    0.0%    0.00B
  ~rhs!~                      4.50k   2.68ms    0.4%   595ns   9.33KiB   11.7%    2.12B
  reset ∂u/∂t                 4.50k   1.95ms    0.3%   433ns     0.00B    0.0%    0.00B
  Jacobian                    4.50k   1.71ms    0.3%   379ns     0.00B    0.0%    0.00B
  prolong2boundaries          4.50k    133μs    0.0%  29.5ns     0.00B    0.0%    0.00B
  prolong2mortars             4.50k    121μs    0.0%  27.0ns     0.00B    0.0%    0.00B
  source terms                4.50k   97.8μs    0.0%  21.7ns     0.00B    0.0%    0.00B
  mortar flux                 4.50k   88.1μs    0.0%  19.6ns     0.00B    0.0%    0.00B
  boundary flux               4.50k   75.0μs    0.0%  16.7ns     0.00B    0.0%    0.00B
analyze solution                  2    897μs    0.1%   449μs   52.6KiB   65.7%  26.3KiB
───────────────────────────────────────────────────────────────────────────────────────

IMEX with sparsity detection:

(This is very slow since I needed to decrease the timestep to make the Newton converge)

───────────────────────────────────────────────────────────────────────────────────────
              Trixi.jl                        Time                    Allocations      
                                     ───────────────────────   ────────────────────────
          Tot / % measured:               68.4s /  25.3%           18.3GiB /   0.0%    

Section                      ncalls     time    %tot     avg     alloc    %tot      avg
───────────────────────────────────────────────────────────────────────────────────────
parabolic rhs!                93.5k    17.1s   98.8%   183μs   18.1KiB   22.6%    0.20B
  calculate gradient          93.5k    7.70s   44.5%  82.3μs   7.34KiB    9.1%    0.08B
    volume integral           93.5k    4.85s   28.1%  51.9μs     0.00B    0.0%    0.00B
    Jacobian                  93.5k    920ms    5.3%  9.84μs     0.00B    0.0%    0.00B
    surface integral          93.5k    595ms    3.4%  6.36μs     0.00B    0.0%    0.00B
    interface flux            93.5k    545ms    3.1%  5.82μs     0.00B    0.0%    0.00B
    prolong2interfaces        93.5k    476ms    2.7%  5.09μs     0.00B    0.0%    0.00B
    reset gradients           93.5k    202ms    1.2%  2.16μs     0.00B    0.0%    0.00B
    ~calculate gradient~      93.5k   90.4ms    0.5%   967ns   7.34KiB    9.1%    0.08B
    prolong2boundaries        93.5k   5.76ms    0.0%  61.6ns     0.00B    0.0%    0.00B
    prolong2mortars           93.5k   5.74ms    0.0%  61.4ns     0.00B    0.0%    0.00B
    mortar flux               93.5k   4.12ms    0.0%  44.1ns     0.00B    0.0%    0.00B
    boundary flux             93.5k   2.72ms    0.0%  29.1ns     0.00B    0.0%    0.00B
  volume integral             93.5k    4.23s   24.4%  45.2μs     0.00B    0.0%    0.00B
  calculate viscous fluxes    93.5k    2.94s   17.0%  31.4μs     0.00B    0.0%    0.00B
  interface flux              93.5k    541ms    3.1%  5.79μs     0.00B    0.0%    0.00B
  prolong2interfaces          93.5k    497ms    2.9%  5.31μs     0.00B    0.0%    0.00B
  transform variables         93.5k    490ms    2.8%  5.24μs     0.00B    0.0%    0.00B
  surface integral            93.5k    369ms    2.1%  3.95μs     0.00B    0.0%    0.00B
  reset ∂u/∂t                 93.5k    128ms    0.7%  1.37μs     0.00B    0.0%    0.00B
  ~parabolic rhs!~            93.5k    119ms    0.7%  1.28μs   10.8KiB   13.4%    0.12B
  Jacobian                    93.5k   69.6ms    0.4%   745ns     0.00B    0.0%    0.00B
  prolong2mortars             93.5k   6.31ms    0.0%  67.5ns     0.00B    0.0%    0.00B
  prolong2boundaries          93.5k   5.51ms    0.0%  58.9ns     0.00B    0.0%    0.00B
  mortar flux                 93.5k   3.93ms    0.0%  42.0ns     0.00B    0.0%    0.00B
  boundary flux               93.5k   2.25ms    0.0%  24.0ns     0.00B    0.0%    0.00B
rhs!                          3.00k    202ms    1.2%  67.4μs   9.33KiB   11.6%    3.18B
  volume integral             3.00k    124ms    0.7%  41.2μs     0.00B    0.0%    0.00B
  interface flux              3.00k   30.5ms    0.2%  10.2μs     0.00B    0.0%    0.00B
  prolong2interfaces          3.00k   14.6ms    0.1%  4.88μs     0.00B    0.0%    0.00B
  surface integral            3.00k   13.4ms    0.1%  4.47μs     0.00B    0.0%    0.00B
  reset ∂u/∂t                 3.00k   6.96ms    0.0%  2.32μs     0.00B    0.0%    0.00B
  ~rhs!~                      3.00k   6.75ms    0.0%  2.25μs   9.33KiB   11.6%    3.18B
  Jacobian                    3.00k   3.12ms    0.0%  1.04μs     0.00B    0.0%    0.00B
  prolong2boundaries          3.00k   1.20ms    0.0%   399ns     0.00B    0.0%    0.00B
  prolong2mortars             3.00k   1.16ms    0.0%   387ns     0.00B    0.0%    0.00B
  mortar flux                 3.00k    571μs    0.0%   190ns     0.00B    0.0%    0.00B
  boundary flux               3.00k    136μs    0.0%  45.4ns     0.00B    0.0%    0.00B
  source terms                3.00k    105μs    0.0%  34.8ns     0.00B    0.0%    0.00B
analyze solution                  2   1.93ms    0.0%   963μs   52.9KiB   65.8%  26.4KiB
───────────────────────────────────────────────────────────────────────────────────────

@cwittens @MarcoArtiano @JoshuaLampert

Question: Should we add a certain semidiscretize (similar to the one for the sparsity-aware Jacobians)? Or is this too specific for right now only advection-diffusion?

Related NumericalMathematics/DispersiveShallowWater.jl#207

[github-actions]

Review checklist

This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging.

Purpose and scope

• The PR has a single goal that is clear from the PR title and/or description.
• All code changes represent a single set of modifications that logically belong together.
• No more than 500 lines of code are changed or there is no obvious way to split the PR into multiple PRs.

Code quality

• The code can be understood easily.
• Newly introduced names for variables etc. are self-descriptive and consistent with existing naming conventions.
• There are no redundancies that can be removed by simple modularization/refactoring.
• There are no leftover debug statements or commented code sections.
• The code adheres to our conventions and style guide, and to the Julia guidelines.

Documentation

• New functions and types are documented with a docstring or top-level comment.
• Relevant publications are referenced in docstrings (see example for formatting).
• Inline comments are used to document longer or unusual code sections.
• Comments describe intent ("why?") and not just functionality ("what?").
• If the PR introduces a significant change or new feature, it is documented in NEWS.md with its PR number.

Testing

• The PR passes all tests.
• New or modified lines of code are covered by tests.
• New or modified tests run in less then 10 seconds.

Performance

• There are no type instabilities or memory allocations in performance-critical parts.
• If the PR intent is to improve performance, before/after time measurements are posted in the PR.

Verification

• The correctness of the code was verified using appropriate tests.
• If new equations/methods are added, a convergence test has been run and the results
  are posted in the PR.

Created with ❤️ by the Trixi.jl community.

[codecov]

Codecov Report

❌ Patch coverage is 96.15385% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 96.85%. Comparing base (fee6f96) to head (4c10ef5).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
...ization/semidiscretization_hyperbolic_parabolic.jl	92.86%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2674      +/-   ##
==========================================
- Coverage   96.85%   96.85%   -0.00%     
==========================================
  Files         546      547       +1     
  Lines       43276    43302      +26     
==========================================
+ Hits        41914    41939      +25     
- Misses       1362     1363       +1     

Flag	Coverage Δ
unittests	96.85% <96.15%> (-<0.01%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ 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.

[MarcoArtiano]

Thanks Daniel, that seems really interesting. Do you think it would be possible to provide a user defined linear operator for a generic set of equations (see IMEX user defined #2518), so that the user can define its own linearized operator (maybe we can choose to update it every linear_operator_interval)? Basically we let the user define how the linear operator should look like for a generic set of equations. Also, if we do that, than in-fact we are doing a Rosenbrock-like method.

[DanielDoehring]

Thanks Daniel, that seems really interesting. Do you think it would be possible to provide a user defined linear operator for a generic set of equations (see IMEX user defined #2518), so that the user can define its own linearized operator (maybe we can choose to update it every linear_operator_interval)? Basically we let the user define how the linear operator should look like for a generic set of equations. Also, if we do that, than in-fact we are doing a Rosenbrock-like method.

I guess you can supply whatever you want, as now only the diffusion part is handed over. So as long you can construct a linear operator (which could, at least according to the documentation, be updated from time to time) this should work. Question is of course again if you have to factorize multiple times how efficient the overall procedure then is.

[DanielDoehring]

Compat is a real joy again...

[ranocha]

Question: Should we add a certain semidiscretize (similar to the one for the sparsity-aware Jacobians)? Or is this too specific for right now only advection-diffusion?

It would be quite specific and would also not work with AMR or MPI...

[DanielDoehring]

Looks like I found settings that work - at the cost of massive updates