MPDeC methods (#146)

* Added MPDeC for conservative PDS

* fix

* first tests

* fix error test

* added TODOs in tests

* added mpdec api reference

* Added MPDeC for non-autonomous PDS

* added adaptive MPDeC

* mpdec interpolation

* MPDeC for PDSProblem (out-of-place)

* added mpdec for non-conservative PDS (in-place)

* bugfix MPDeC for PDSProblem & additional tests

* fix testset 'zero intial values'

* code restructuring

* performance enhancement

* check theta < 0 only once

* Support Tridiagonal

* Support for sparse matrices

* fixed sparse matrix support

* minor changes and tests

* format

* fix - change in sparsity pattern

* Use POSITIVE PDS in tests

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* additiional changes after review

* Removed breaks in loops for sparse matrices

* fixed breaks in loops for sparse matrices

* use algebraic expressions for integration nodes and weights

* fix strato benchmark

* Convergence tests for high order schemes

* Updated NEWS.md

* Added DoubleFloats to test environment

* made SVector available in tests

* made SMatrix available in tests

* import only Double64 from DoubleFloats

* added suggested changes from reviews

* format

* use better variable name M_plus_1

---------

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

Author: SKopecz

NEWS.md

@@ -12,6 +12,9 @@ for human readability.
 
 - The minimum required Julia version was updated to v1.10
 
+### Added
+
+- MPDeC methods of Öffner and Torlo
 
 ## Changes when updating to v0.2 from v0.1.x
 

docs/src/api_reference.md

@@ -35,6 +35,7 @@ SSPMPRK22
 MPRK43I
 MPRK43II
 SSPMPRK43
+MPDeC
 ```
 
 ## Auxiliary functions

docs/src/stratospheric_reaction_benchmark.md

@@ -23,7 +23,7 @@ function stratreac_plot(sols, labels = fill("", length(sols)), sol_ref = nothing
         labels = [labels]
     end
 
-    tspan = prob.tspan
+    tspan = prob_pds_stratreac.tspan
     layout = (3, 2)
     linewidth = 2
     xticks = (range(first(tspan), last(tspan), 4), range(12.0, 84.0, 4))

src/PositiveIntegrators.jl

@@ -6,7 +6,7 @@ using Statistics: median
 
 using SparseArrays: SparseArrays, AbstractSparseMatrix,
                     issparse, nonzeros, nzrange, rowvals, spdiagm
-using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix
+using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix, MMatrix
 
 using FastBroadcast: @..
 using MuladdMacro: @muladd
@@ -48,6 +48,7 @@ export ConservativePDSFunction, ConservativePDSProblem
 
 export MPE, MPRK22, MPRK43I, MPRK43II
 export SSPMPRK22, SSPMPRK43
+export MPDeC
 
 export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod,
        prob_pds_robertson, prob_pds_brusselator, prob_pds_sir,
@@ -69,6 +70,9 @@ include("mprk.jl")
 # modified Patankar-Runge-Kutta based on the SSP formulation of RK methods (SSPMPRK)
 include("sspmprk.jl")
 
+# MPDeC methods
+include("mpdec.jl")
+
 # interpolation for dense output
 include("interpolation.jl")
 

src/interpolation.jl

@@ -46,7 +46,8 @@ const MPRKCaches = Union{MPEConstantCache, MPECache, MPEConservativeCache,
                          MPRK22ConstantCache, MPRK22Cache, MPRK22ConservativeCache,
                          MPRK43ConstantCache, MPRK43Cache, MPRK43ConservativeCache,
                          SSPMPRK22ConstantCache, SSPMPRK22Cache, SSPMPRK22ConservativeCache,
-                         SSPMPRK43ConstantCache, SSPMPRK43Cache, SSPMPRK43ConservativeCache}
+                         SSPMPRK43ConstantCache, SSPMPRK43Cache, SSPMPRK43ConservativeCache,
+                         MPDeCConstantCache, MPDeCCache, MPDeCConservativeCache}
 
 function interp_summary(::Type{cacheType},
                         dense::Bool) where {cacheType <: MPRKCaches}

src/mpdec.jl

@@ -1,6 +1,7 @@
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
@@ -21,6 +22,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 ADTypes = "1.12"
 Aqua = "0.8"
+DoubleFloats = "1.4.3"
 ExplicitImports = "1.0.1"
 LinearAlgebra = "1"
 LinearSolve = "3.7.1"