WIP: Linear multistep Patankar schemes

[SKopecz]

This implements #107.

• Out-of-place implementations
• In-place implementations
• Check stats count
• Tests
• Docs
• Additional test problems
• Benchmark diffusion test problem
• Ensure that dt remains constant for all times (no changes by callbacks allowed), see AB3 (and other multistep methods) fails to reinitialize after callback / discontinuity because integrator.derivative_discontinuity is reset too early SciML/DifferentialEquations.jl#1154

[SKopecz]

All the onestep MPRK/SSPMPRK and also the multistep Patankar schemes can be interpreted as sequences of MPE steps with slightly modified input data. Before continuing this PR we should implement a function like _permform_step_MPE! and rewrite the implementations of MPRK and SSPMPRK schemes. This will make the implementation of the multistep Patankar schemes easier and cleaner. The corresponding PR is #183 .

[codecov]

Codecov Report

❌ Patch coverage is 16.35389% with 312 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/mplm.jl	0.00%	296 Missing ⚠️
src/mprk.jl	27.27%	16 Missing ⚠️

📢 Thoughts on this report? Let us know!

[coveralls]

Coverage Report for CI Build 28241095053

Coverage increased (+0.9%) to 98.233%

Details

• Coverage increased (+0.9%) from the base build.
• Patch coverage: 17 uncovered changes across 1 file (1844 of 1861 lines covered, 99.09%).
• No coverage regressions found.

Uncovered Changes

File	Changed	Covered	%
src/mplm.jl	1765	1748	99.04%
Total (3 files)	1861	1844	99.09%

Coverage Regressions

No coverage regressions found.

---

Coverage Stats

Relevant Lines:	3566
Covered Lines:	3503
Line Coverage:	98.23%
Coverage Strength:	61897055.86 hits per line

---

💛 - Coveralls

[SKopecz]

The new problems need to be added here, to be found in the tests.

[SKopecz]

The authors of the original paper kindly shared their code; consequently, the current out-of-place implementations of MPLM22, MPLM33, and MPLM43 yield the same results as those in the paper.

However, the implementation of the startup phase for these m-step methods is quite cumbersome and difficult to maintain. To achieve the required accuracy for the starting values, the authors utilize a sub-stepping strategy (dt/4) with nested lower-order MPLM schemes. This currently necessitates a redundant implementation: one for the standard perform_step! and another for a perform_substeps_MPLMXX function, which is used within the perform_step! of higher-order methods to compute initial values. This "double implementation" is error-prone and leads to significant code duplication.

Should we keep this sub-stepping strategy? If so, is there a way to eliminate the "double implementation"? Alternatively, we could use MPRK or MPDeC one-step schemes to compute the initial values. Since these are available for all orders, the reduced time step would not be necessary. However, this would also require implementations of these methods that can be integrated into the MPLM perform_step!.

@ranocha @JoshuaLampert Do you have any suggestions?

[SKopecz]

Currently, the multistep coefficients of all MPLM schemes are based on the original paper and are hardcoded. We should decide between two directions:

Simplify code by keeping original coefficients

If we restrict the implementation strictly to the coefficients of the original paper, the code can be significantly simplified since many of these coefficients are zero. For instance, for MPLM106 we have α1 = ... = α9 = 0. This means the current implementation

    @.. broadcast=false linsolve.b=α1 * uprev + α2 * uprevprev + α3 * uprev3 + α4 * uprev4 +
                                   α5 * uprev5 + α6 * uprev6 + α7 * uprev7 + α8 * uprev8 +
                                   α9 * uprev9 + α10 * uprev10

could be simplified to

@.. broadcast=false linsolve.b = uprev10

Allow for user-defined coefficients

Alternatively, we could allow users to specify custom multistep coefficients. In this case, we would need to introduce a constructor likeMPLM106(α, β).

Technical detail: In each step, the denominators of the Patankar weights are computed via recursive use of lower-order MPLM schemes. The coefficients of these underlying schemes would remain fixed.