Careful handling of data types (Float32 & units) (#108)

* show u and sigma in MPRK43 during perform_step

* more output

* more output

* more output

* changed MPRK43II small_constant

* removed addtional output

* format

* updated doc strings

* added conversion to Float32 in SSPMPRK43

* Added small_constant_function_MPRK43II

* revised docstrings

* format

* typo

* Evaluation of PDSFunction and ConservativePDSFunction can handle units

* structure runtest

* bugfix runtest

* additoinal tests

* Use Unitful in tests

* More functionality from StaticArrays

* bugfixes

* Reset evaluation of PDSFunctions with sparse prototyp

* bugfix

* Evaluation of PDSFunctions with sparse Prototype and units is working

* format

* Enable solution of PDS with unity by OrdinaryDiffEq solvers

* test more OrdinaryDiffEq algorithms solving problems with units

* Float32 test

* removed tests for solving PDSProblems with units by implicit solvers

* bugfix test Compatibility of OrdinarDiffEq solvers with (Conservative)PDSProblems and units

* Don't use implicit OrdinaryDiffEq solvers with units in tests

* Update src/mprk.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/mprk.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/proddest.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/sspmprk.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update test/runtests.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* use unitful.jl's ustirp in tests

* Compare unitful data without ustrip

* Use ustrip before comparing solutions with units in tests

* Added missing test (SSPMPRK43 and Float32)

* similar and oneunit

* similar(P[:,1])

* set version to v0.2.6

---------

Co-authored-by: Hendrik Ranocha <[email protected]>

Author: SKopecz

Project.toml

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

src/mprk.jl

@@ -987,8 +987,9 @@ You can optionally choose the linear solver to be used by passing an
 algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
-to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
-`floatmin` of the floating point type used.
+to avoid divisions by zero. To display the default value for data type `type` evaluate 
+`MPRK43II(gamma).small_constant_function(type)`, where `type` can be, e.g.,
+`Float64`.
 
 ## References
 
@@ -1004,10 +1005,18 @@ struct MPRK43II{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm
     small_constant_function::T2
 end
 
-function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = nothing)
-    if isnothing(small_constant)
-        small_constant_function = floatmin
-    elseif small_constant isa Number
+function small_constant_function_MPRK43II(type)
+    if type == Float64
+        small_constant = 1e-50
+    else
+        small_constant = floatmin(type)
+    end
+    return small_constant
+end
+
+function MPRK43II(gamma; linsolve = LUFactorization(),
+                  small_constant = small_constant_function_MPRK43II)
+    if small_constant isa Number
         small_constant_function = Returns(small_constant)
     else # assume small_constant isa Function
         small_constant_function = small_constant

src/proddest.jl

@@ -95,11 +95,11 @@ function PDSProblem{iip}(P, D, u0, tspan, p = NullParameters();
 
     # p_prototype is used to store evaluations of P, if P is in-place.
     if isnothing(p_prototype) && iip
-        p_prototype = zeros(eltype(u0), (length(u0), length(u0)))
+        p_prototype = zeros(eltype(u0), (length(u0), length(u0))) / oneunit(first(tspan))
     end
     # If a PDSFunction is to be evaluated and D is in-place, then d_prototype is used to store
     # evaluations of D.
-    d_prototype = similar(u0)
+    d_prototype = similar(u0 ./ oneunit(first(tspan)))
 
     PD = PDSFunction{iip}(P, D; p_prototype, d_prototype,
                           analytic, std_rhs)
@@ -121,7 +121,7 @@ end
 # Most specific constructor for PDSFunction
 function PDSFunction{iip, FullSpecialize}(P, D;
                                           p_prototype = nothing,
-                                          d_prototype,
+                                          d_prototype = nothing,
                                           analytic = nothing,
                                           std_rhs = nothing) where {iip}
     if std_rhs === nothing
@@ -143,11 +143,22 @@ function (PD::PDSFunction)(du, u, p, t)
 end
 
 # Default implementation of the standard right-hand side evaluation function
-struct PDSStdRHS{P, D, PrototypeP, PrototypeD} <: Function
+struct PDSStdRHS{P, D, PrototypeP, PrototypeD, TMP} <: Function
     p::P
     d::D
     p_prototype::PrototypeP
     d_prototype::PrototypeD
+    tmp::TMP
+end
+
+function PDSStdRHS(P, D, p_prototype, d_prototype)
+    if p_prototype isa AbstractSparseMatrix
+        tmp = zeros(eltype(p_prototype), (size(p_prototype, 1),)) /
+              oneunit(first(p_prototype)) # drop units
+    else
+        tmp = nothing
+    end
+    PDSStdRHS(P, D, p_prototype, d_prototype, tmp)
 end
 
 # Evaluation of a PDSStdRHS (out-of-place)
@@ -163,9 +174,10 @@ function (PD::PDSStdRHS)(du, u, p, t)
     PD.p(PD.p_prototype, u, p, t)
 
     if PD.p_prototype isa AbstractSparseMatrix
-        # Same result but more efficient - at least currently for SparseMatrixCSC
-        fill!(PD.d_prototype, one(eltype(PD.d_prototype)))
-        mul!(vec(du), PD.p_prototype, PD.d_prototype)
+        # row sum coded as matrix-vector product 
+        fill!(PD.tmp, one(eltype(PD.tmp)))
+        mul!(vec(du), PD.p_prototype, PD.tmp)
+
         for i in 1:length(u)  #vec(du) .+= diag(PD.p_prototype)
             du[i] += PD.p_prototype[i, i]
         end
@@ -272,7 +284,7 @@ function ConservativePDSProblem{iip}(P, u0, tspan, p = NullParameters();
 
     # p_prototype is used to store evaluations of P, if P is in-place.
     if isnothing(p_prototype) && iip
-        p_prototype = zeros(eltype(u0), (length(u0), length(u0)))
+        p_prototype = zeros(eltype(u0), (length(u0), length(u0))) / oneunit(first(tspan))
     end
 
     PD = ConservativePDSFunction{iip}(P; p_prototype, analytic, std_rhs)
@@ -314,27 +326,30 @@ function (PD::ConservativePDSFunction)(du, u, p, t)
 end
 
 # Default implementation of the standard right-hand side evaluation function
-struct ConservativePDSStdRHS{P, PrototypeP, TMP} <: Function
+struct ConservativePDSStdRHS{P, PrototypeP, TMP, TMP2} <: Function
     p::P
     p_prototype::PrototypeP
     tmp::TMP
+    tmp2::TMP2
 end
 
 function ConservativePDSStdRHS(P, p_prototype)
     if p_prototype isa AbstractSparseMatrix
         tmp = zeros(eltype(p_prototype), (size(p_prototype, 1),))
+        tmp2 = tmp / oneunit(first(tmp)) # drop units
     else
         tmp = nothing
+        tmp2 = nothing
     end
-    ConservativePDSStdRHS(P, p_prototype, tmp)
+    ConservativePDSStdRHS(P, p_prototype, tmp, tmp2)
 end
 
 # Evaluation of a ConservativePDSStdRHS (out-of-place)
 function (PD::ConservativePDSStdRHS)(u, p, t)
     #vec(sum(PD.p(u, p, t), dims = 2)) - vec(sum(PD.p(u, p, t), dims = 1))
     P = PD.p(u, p, t)
 
-    f = zero(u)
+    f = zeros(eltype(P), size(u))
     @fastmath @inbounds @simd for I in CartesianIndices(P)
         if !iszero(P[I])
             f[I[1]] += P[I]
@@ -347,8 +362,8 @@ end
 function (PD::ConservativePDSStdRHS)(u::SVector, p, t)
     P = PD.p(u, p, t)
 
-    f = similar(u) #constructs MVector
-    zeroT = zero(eltype(u))
+    f = similar(P[:, 1]) #constructs MVector
+    zeroT = zero(eltype(P))
     for i in eachindex(f)
         f[i] = zeroT
     end
@@ -366,32 +381,37 @@ end
 # Evaluation of a ConservativePDSStdRHS (in-place)
 function (PD::ConservativePDSStdRHS)(du, u, p, t)
     PD.p(PD.p_prototype, u, p, t)
-    sum_terms!(du, PD.tmp, PD.p_prototype)
+    sum_terms!(du, PD.tmp, PD.tmp2, PD.p_prototype)
     return nothing
 end
 
 # Generic fallback (for dense arrays)
 # This implementation does not need any auxiliary vectors
-@inline function sum_terms!(du, tmp, P)
-    for i in 1:length(du)
+@inline function sum_terms!(du, tmp, tmp2, P)
+    for i in eachindex(du)
         du[i] = zero(eltype(du))
-        for j in 1:length(du)
+        for j in eachindex(du)
             du[i] += P[i, j] - P[j, i]
         end
     end
     return nothing
 end
 
 # Same result but more efficient - at least currently for SparseMatrixCSC
-@inline function sum_terms!(du, tmp, P::AbstractSparseMatrix)
-    fill!(tmp, one(eltype(tmp)))
-    mul!(vec(du), P, tmp)
+@inline function sum_terms!(du, tmp, tmp2, P::AbstractSparseMatrix)
+    # row sum coded as matrix vector product
+    fill!(tmp2, one(eltype(tmp2)))
+    mul!(vec(du), P, tmp2)
+    #sum!(vec(du), P)
+
+    #column sum
     sum!(tmp', P)
+
     vec(du) .-= tmp
     return nothing
 end
 
-@inline function sum_terms!(du, tmp, P::Tridiagonal)
+@inline function sum_terms!(du, tmp, tmp2, P::Tridiagonal)
     Base.require_one_based_indexing(du, P.dl, P.du)
     @assert length(du) == length(P.dl) + 1 == length(P.du) + 1
 

src/sspmprk.jl

@@ -466,8 +466,9 @@ You can optionally choose the linear solver to be used by passing an
 algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
 as keyword argument `linsolve`.
 You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators
-to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
-`floatmin` of the floating point type used.
+to avoid divisions by zero. To display the default value for data type `type` evaluate 
+`SSPMPRK43. small_constant_function(type)`, where `type` can be, e.g.,
+`Float64`.
 
 The current implementation only supports fixed time steps.
 
@@ -488,10 +489,26 @@ struct SSPMPRK43{F, T} <: OrdinaryDiffEqAlgorithm
     small_constant_function::T
 end
 
-function SSPMPRK43(; linsolve = LUFactorization(), small_constant = 1e-50)
-    if isnothing(small_constant)
-        small_constant_function = floatmin
-    elseif small_constant isa Number
+function small_constant_function_SSPMPRK43(type)
+    if type == Float64
+        small_constant = 1e-50
+    elseif type == Float32
+        # small_constant is chosen such that the problem below
+        # (zero initial condition) can be solved
+        # P_linmod(u, p, t) = [0 u[2]; 5*u[1] 0]
+        # u0 = [1.0f0, 0.0f0]
+        # prob = ConservativePDSProblem(P_linmod, u0, (0.0f0, 2.0f0))
+        # sol = solve(prob, SSPMPRK43(); dt=0.1f0)
+        small_constant = 1.0f-8
+    else
+        small_constant = floatmin(type)
+    end
+    return small_constant
+end
+
+function SSPMPRK43(; linsolve = LUFactorization(),
+                   small_constant = small_constant_function_SSPMPRK43)
+    if small_constant isa Number
         small_constant_function = Returns(small_constant)
     else # assume small_constant isa Function
         small_constant_function = small_constant
@@ -534,8 +551,7 @@ function get_constant_parameters(alg::SSPMPRK43)
     c3 = β20 + α21 * β10 + β21
 
     return n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20,
-           β21, β30,
-           β31, β32, c3
+           β21, β30, β31, β32, c3
 end
 
 struct SSPMPRK43ConstantCache{T} <: OrdinaryDiffEqConstantCache
@@ -573,11 +589,12 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("SSPMPRK43 can only be applied to production-destruction systems"))
     end
-    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = get_constant_parameters(alg)
+    const_param = get_constant_parameters(alg)
+    const_param = convert.(uEltypeNoUnits, const_param)
+    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = const_param
+    small_constant = alg.small_constant_function(uEltypeNoUnits)
     SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31,
-                           α32, β10,
-                           β20, β21, β30,
-                           β31, β32, c3, alg.small_constant_function(uEltypeNoUnits))
+                           α32, β10, β20, β21, β30, β31, β32, c3, small_constant)
 end
 
 function initialize!(integrator, cache::SSPMPRK43ConstantCache)

test/Project.toml

@@ -9,6 +9,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 Aqua = "0.7, 0.8"
@@ -18,3 +19,4 @@ OrdinaryDiffEq = "6.59"
 Plots = "1.11.1"
 StaticArrays = "1.5"
 Statistics = "1"
+Unitful = "1.21"