Merge pull request #86 from SKopecz/SK/SSPMPRK

SSPMPRK schemes and small_constant as algorithm field

Project.toml

@@ -1,7 +1,7 @@
 name = "PositiveIntegrators"
 uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5"
 authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"]
-version = "0.1.14"
+version = "0.1.15"
 
 [deps]
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"

src/PositiveIntegrators.jl

@@ -48,6 +48,7 @@ export PDSFunction, PDSProblem
 export ConservativePDSFunction, ConservativePDSProblem
 
 export MPE, MPRK22, MPRK43I, MPRK43II
+export SSPMPRK22, SSPMPRK43
 
 export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod,
        prob_pds_robertson, prob_pds_brusselator, prob_pds_sir,
@@ -61,6 +62,11 @@ include("proddest.jl")
 # modified Patankar-Runge-Kutta (MPRK) methods
 include("mprk.jl")
 
+# modified Patankar-Runge-Kutta based on the SSP formulation of RK methods (SSPMPRK)
+include("sspmprk.jl")
+
+include("interpolation.jl")
+
 # predefined PDS problems
 include("PDSProblemLibrary.jl")
 

src/interpolation.jl

@@ -0,0 +1,86 @@
+# Linear interpolations
+@muladd @inline function linear_interpolant(Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{0}})
+    Θm1 = (1 - Θ)
+    @.. broadcast=false Θm1 * u0+Θ * u1
+end
+
+@muladd @inline function linear_interpolant(Θ, dt, u0, u1, idxs, T::Type{Val{0}})
+    Θm1 = (1 - Θ)
+    @.. broadcast=false Θm1 * u0[idxs]+Θ * u1[idxs]
+end
+
+@muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing,
+                                             T::Type{Val{0}})
+    Θm1 = (1 - Θ)
+    @.. broadcast=false out=Θm1 * u0 + Θ * u1
+    out
+end
+
+@muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{0}})
+    Θm1 = (1 - Θ)
+    @views @.. broadcast=false out=Θm1 * u0[idxs] + Θ * u1[idxs]
+    out
+end
+
+@inline function linear_interpolant(Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{1}})
+    @.. broadcast=false (u1 - u0)/dt
+end
+
+@inline function linear_interpolant(Θ, dt, u0, u1, idxs, T::Type{Val{1}})
+    @.. broadcast=false (u1[idxs] - u0[idxs])/dt
+end
+
+@inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{1}})
+    @.. broadcast=false out=(u1 - u0) / dt
+    out
+end
+
+@inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{1}})
+    @views @.. broadcast=false out=(u1[idxs] - u0[idxs]) / dt
+    out
+end
+
+#######################################################################################
+# interpolation specializations
+const MPRKCaches = Union{MPEConstantCache, MPECache, MPEConservativeCache,
+                         MPRK22ConstantCache, MPRK22Cache, MPRK22ConservativeCache,
+                         MPRK43ConstantCache, MPRK43Cache, MPRK43ConservativeCache,
+                         SSPMPRK22ConstantCache, SSPMPRK22Cache, SSPMPRK22ConservativeCache,
+                         SSPMPRK43ConstantCache, SSPMPRK43Cache, SSPMPRK43ConservativeCache}
+
+function interp_summary(::Type{cacheType},
+                        dense::Bool) where {cacheType <: MPRKCaches}
+    "1st order linear"
+end
+
+function _ode_interpolant(Θ, dt, u0, u1, k,
+                          cache::MPRKCaches,
+                          idxs, # Optionally specialize for ::Nothing and others
+                          T::Type{Val{0}},
+                          differential_vars::Nothing)
+    linear_interpolant(Θ, dt, u0, u1, idxs, T)
+end
+
+function _ode_interpolant!(out, Θ, dt, u0, u1, k,
+                           cache::MPRKCaches,
+                           idxs, # Optionally specialize for ::Nothing and others
+                           T::Type{Val{0}},
+                           differential_vars::Nothing)
+    linear_interpolant!(out, Θ, dt, u0, u1, idxs, T)
+end
+
+function _ode_interpolant(Θ, dt, u0, u1, k,
+                          cache::MPRKCaches,
+                          idxs, # Optionally specialize for ::Nothing and others
+                          T::Type{Val{1}},
+                          differential_vars::Nothing)
+    linear_interpolant(Θ, dt, u0, u1, idxs, T)
+end
+
+function _ode_interpolant!(out, Θ, dt, u0, u1, k,
+                           cache::MPRKCaches,
+                           idxs, # Optionally specialize for ::Nothing and others
+                           T::Type{Val{1}},
+                           differential_vars::Nothing)
+    linear_interpolant!(out, Θ, dt, u0, u1, idxs, T)
+end