WIP: Simplified and cleaner code using the function basic_patankar_step.

[SKopecz]

All MPRK and SSPMPRK schemes are build upon a basic Patankar step with varying Patankar matrices, destruction vectors and denominators. This PR introduces the function basic_patankar_step (better name suggestions are welcome), which makes the implementations cleaner and can also be used to implement the multistep Patankar schemes in #107.

1. Out-of-place implementations

• MPE
• MPRK22
• MPRK43I & MPRK43II
• SSPMPRK22
• SSPMPRK43

2. In-place implementations

• MPE
• MPRK22
• MPRK43I & MPRK43II
• SSPMPRK22
• SSPMPRK43

[codecov]

Codecov Report

❌ Patch coverage is 97.05882% with 2 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/mprk.jl	95.55%	2 Missing ⚠️

📢 Thoughts on this report? Let us know!

[coveralls]

Pull Request Test Coverage Report for Build 20692027876

Warning: This coverage report may be inaccurate.

This pull request's base commit is no longer the HEAD commit of its target branch. This means it includes changes from outside the original pull request, including, potentially, unrelated coverage changes.

• For more information on this, see Tracking coverage changes with pull request builds.
• To avoid this issue with future PRs, see these Recommended CI Configurations.
• For a quick fix, rebase this PR at GitHub. Your next report should be accurate.

Details

• 66 of 68 (97.06%) changed or added relevant lines in 2 files are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage decreased (-0.2%) to 97.612%

---

Changes Missing Coverage	Covered Lines	Changed/Added Lines	%
src/mprk.jl	43	45	95.56%

Totals
Change from base Build 20199034964:	-0.2%
Covered Lines:	1962
Relevant Lines:	2010

---

💛 - Coveralls

[SKopecz]

Introduced a function to compute linear combinations for out-of-place computations. Although implemented recursively, it's internally unrolled.

using BenchmarkTools, StaticArrays 

# --- Base cases (End of recursion) ---

# For PDSFunctions (with destruction vectors)
@inline lincomb(c1::Number, P1, d1::AbstractArray) = (c1 .* P1, c1 .* d1)

# For ConservativePDSFunctions (without destruction vectors)
@inline lincomb(c1::Number, P1, d1::Nothing) = (c1 .* P1, nothing)


# --- Recursive steps ---

# For PDSFunctions: Processes triplets (coeff, P, d)
@inline function lincomb(c1::Number, P1, d1::AbstractArray, c2, P2, d2, args...)
    P_tail, d_tail = lincomb(c2, P2, d2, args...)
    return (c1 .* P1 .+ P_tail, c1 .* d1 .+ d_tail)
end

# For ConservativePDSFunctions: Processes triplets (coeff, P, nothing)
@inline function lincomb(c1::Number, P1, d1::Nothing, c2, P2, d2, args...)
    P_tail, _ = lincomb(c2, P2, d2, args...)
    return (c1 .* P1 .+ P_tail, nothing)
end

function run_comprehensive_benchmark(N)
    println("\n" * "█"^70)
    println("  BENCHMARKING SYSTEM SIZE N = $N")
    println("█"^70)

    # Setup test data (3 stages)
    P1, P2, P3 = rand(SMatrix{N,N}), rand(SMatrix{N,N}), rand(SMatrix{N,N})
    d1, d2, d3 = rand(SVector{N}), rand(SVector{N}), rand(SVector{N})
    c1, c2, c3 = 0.5, 0.3, 0.2

    # --- CASE A: Standard PDS (with destruction vectors) ---
    println("\n>>> CASE A: Standard PDS (d is AbstractArray)")
    
    println("\n[A.1] Direct Input (Manual):")
    b_direct_a = @benchmark ($c1 .* $P1 .+ $c2 .* $P2 .+ $c3 .* $P3, 
                             $c1 .* $d1 .+ $c2 .* $d2 .+ $c3 .* $d3)
    display(b_direct_a)

    println("\n[A.2] lincomb Function:")
    b_func_a = @benchmark lincomb($c1, $P1, $d1, $c2, $P2, $d2, $c3, $P3, $d3)
    display(b_func_a)


    # --- CASE B: Conservative PDS (d is nothing) ---
    println("\n" * "-"^70)
    println(">>> CASE B: Conservative PDS (d is nothing)")
    
    println("\n[B.1] Direct Input (Manual):")
    # In manual case, we only calculate P
    b_direct_b = @benchmark ($c1 .* $P1 .+ $c2 .* $P2 .+ $c3 .* $P3, nothing)
    display(b_direct_b)

    println("\n[B.2] lincomb Function:")
    b_func_b = @benchmark lincomb($c1, $P1, nothing, $c2, $P2, nothing, $c3, $P3, nothing)
    display(b_func_b)
    
    return nothing
end

Results:

julia> run_comprehensive_benchmark(2)

██████████████████████████████████████████████████████████████████████
  BENCHMARKING SYSTEM SIZE N = 2
██████████████████████████████████████████████████████████████████████

>>> CASE A: Standard PDS (d is AbstractArray)

[A.1] Direct Input (Manual):
BenchmarkTools.Trial: 10000 samples with 999 evaluations per sample.
 Range (min … max):  6.135 ns … 183.600 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     7.063 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   7.266 ns ±   3.170 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

    ▅ █ ▇ ▄  ▁                                                 
  ▇██▃█▃█▃█▆▃█▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂ ▃
  6.14 ns         Histogram: frequency by time        14.8 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[A.2] lincomb Function:
BenchmarkTools.Trial: 10000 samples with 999 evaluations per sample.
 Range (min … max):  5.700 ns … 152.339 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     6.569 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   6.755 ns ±   2.887 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

     █ █     ▁ ▂ ▅ ▆                                           
  ▇▅▅█▄█▄█▆▇▅█▅█▃█▄█▄▃▇▃▇▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▂▁▂▂▁▁▂▂▂▂▂▂▁▁▂▂ ▃
  5.7 ns          Histogram: frequency by time        10.4 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

----------------------------------------------------------------------
>>> CASE B: Conservative PDS (d is nothing)

[B.1] Direct Input (Manual):
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  4.823 ns … 174.098 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     5.560 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.677 ns ±   2.636 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅  ▆  ▅  ▂   █                                               
  █▃▄█▃▃█▄▃█▆▃▃█▇▃▃▅▇▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▂▁▂▂▂▂▂▂▂▁▁▂▂▁▂▂▁▂ ▃
  4.82 ns         Histogram: frequency by time        9.31 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[B.2] lincomb Function:
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  4.823 ns … 70.316 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     5.562 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.758 ns ±  1.852 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▂█▄▇▄▃▇▃▆█▄▃▆▃▂▁▁▁▁▁                                      ▃
  █████████████████████▇▇█▇▇██▇██████▇▇▇▇▇▇▆▇▆▅▆▆▁▅▅▅▅▄▄▄▃▃▄ █
  4.82 ns      Histogram: log(frequency) by time     10.8 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> run_comprehensive_benchmark(3)

██████████████████████████████████████████████████████████████████████
  BENCHMARKING SYSTEM SIZE N = 3
██████████████████████████████████████████████████████████████████████

>>> CASE A: Standard PDS (d is AbstractArray)

[A.1] Direct Input (Manual):
BenchmarkTools.Trial: 10000 samples with 999 evaluations per sample.
 Range (min … max):   9.181 ns … 153.478 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     10.074 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   10.712 ns ±   3.368 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▇▃▆▃▆▂▆▄▁▅▁                                                 ▂
  █████████████▇▇▆▅▆▅▄▆▆▆▇▇▇▅▆▇▇▇▇▆▇▆▆▅▅▆▂▆▅▆▅▅▆▄▅▅▅▅▄▃▃▄▂▄▃▂▃ █
  9.18 ns       Histogram: log(frequency) by time      24.1 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[A.2] lincomb Function:
BenchmarkTools.Trial: 10000 samples with 999 evaluations per sample.
 Range (min … max):  8.308 ns … 175.174 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     8.814 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   9.438 ns ±   3.640 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▂▇▃▂▆▂▁▅▁  ▇▂▁▄▅                                           ▂
  ███████████▆██████▇▇▆▆▆▄▄▅▄▅▅▃▅▅▃▄▄▄▂▅▅▄▄▅▄▅▃▃▄▄▄▄▃▄▃▄▂▄▂▄▄ █
  8.31 ns      Histogram: log(frequency) by time      16.9 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

----------------------------------------------------------------------
>>> CASE B: Conservative PDS (d is nothing)

[B.1] Direct Input (Manual):
BenchmarkTools.Trial: 10000 samples with 999 evaluations per sample.
 Range (min … max):  6.126 ns … 166.082 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     6.728 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   7.217 ns ±   3.308 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▂▇▁▆▁▅▅▁▇▁▁▅▁▁                                             ▂
  ███████████████▇▇▇▆▅▆▅▅▆▇▇▅▆▅▆▆▆▆▇▇▆▇▇▆▇▆▆▇▅▅▅▅▆▆▅▅▆▇▅▃▄▅▅▄ █
  6.13 ns      Histogram: log(frequency) by time      14.1 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[B.2] lincomb Function:
BenchmarkTools.Trial: 10000 samples with 999 evaluations per sample.
 Range (min … max):  6.127 ns … 149.530 ns  ┊ GC (min … max): 0.00% … 0.00%