From 4cc77e36bf20e697b97aafcce45656f5d9aacc40 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 9 Apr 2024 09:11:59 +0200 Subject: [PATCH] use linear interpolations (#62) * use linear interpolation for MPE * format * fix undefined reference * test also MVector and out-of-place arrays * in-place MPE * update MPRK22 - needs to be checked * format * fix convergence tests * bump version --- Project.toml | 2 +- src/PositiveIntegrators.jl | 24 ++-- src/mprk.jl | 225 ++++++++++++++++++++++++++++--------- test/Project.toml | 2 + test/runtests.jl | 19 +++- 5 files changed, 207 insertions(+), 65 deletions(-) diff --git a/Project.toml b/Project.toml index c5266064..535766c8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PositiveIntegrators" uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5" authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"] -version = "0.1.6-pre" +version = "0.1.6" [deps] FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 93abac27..7cbe494e 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -24,20 +24,24 @@ using SymbolicIndexingInterface using LinearSolve: LinearSolve, LinearProblem, LUFactorization -import SciMLBase: __has_mass_matrix, __has_analytic, __has_tgrad, +using SciMLBase: DEFAULT_OBSERVED +import SciMLBase: interp_summary, + __has_mass_matrix, __has_analytic, __has_tgrad, __has_jac, __has_jvp, __has_vjp, __has_jac_prototype, __has_sparsity, __has_Wfact, __has_Wfact_t, __has_paramjac, __has_syms, __has_indepsym, __has_paramsyms, - __has_observed, DEFAULT_OBSERVED, __has_colorvec, - __has_sys - -import OrdinaryDiffEq: @cache, - OrdinaryDiffEqAdaptiveAlgorithm, calculate_residuals, - calculate_residuals!, False, - OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + __has_observed, __has_colorvec, __has_sys + +using OrdinaryDiffEq: @cache, + DEFAULT_PRECS, + OrdinaryDiffEqAdaptiveAlgorithm, + OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, + False, + recursivefill!, _vec, wrapprecs, dolinsolve +import OrdinaryDiffEq: alg_order, isfsal, + calculate_residuals, calculate_residuals!, alg_cache, initialize!, perform_step!, - recursivefill!, _vec, DEFAULT_PRECS, wrapprecs, - dolinsolve, alg_order + _ode_interpolant, _ode_interpolant! # 2. Export functionality defining the public API export PDSFunction, PDSProblem diff --git a/src/mprk.jl b/src/mprk.jl index bc59161b..5528757a 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -176,6 +176,49 @@ function sum_destruction_terms!(D, P::Tridiagonal) return D end +##################################################################### +# 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 + ### MPE ##################################################################################### """ MPE([linsolve = ...]) @@ -221,7 +264,8 @@ function Base.getproperty(alg::MPE, f::Symbol) end end -alg_order(alg::MPE) = 1 +alg_order(::MPE) = 1 +isfsal(::MPE) = false struct MPECache{uType, rateType, PType, F, uNoUnitsType} <: OrdinaryDiffEqMutableCache u::uType @@ -241,9 +285,9 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} tmp = zero(u) + P = p_prototype(u, f) linsolve_tmp = zero(u) - weight = similar(u, uEltypeNoUnits) recursivefill!(weight, false) @@ -271,15 +315,23 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits}, end function initialize!(integrator, cache::MPEConstantCache) - integrator.kshortsize = 2 + integrator.kshortsize = 1 integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) - integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal - integrator.stats.nf += 1 # Avoid undefined entries if k is an array of arrays - integrator.fsallast = zero(integrator.fsalfirst) - integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast + integrator.fsalfirst = zero(integrator.u) + integrator.fsallast = integrator.fsalfirst + integrator.k[1] = integrator.fsallast + + # TODO: Do we need to set fsalfirst here? The other non-FSAL caches + # in OrdinaryDiffEq.jl use something like + # integrator.fsalfirst = integrator.f(integrator.uprev, integrator, + # integrator.t) # Pre-start fsal + # integrator.stats.nf += 1 + # integrator.fsallast = zero(integrator.fsalfirst) + # integrator.k[1] = integrator.fsalfirst + # Do we need something similar here to get a cache for k values + # with the correct units? end function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false) @@ -287,11 +339,13 @@ function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false) @unpack small_constant = cache # Attention: Implementation assumes that the pds is conservative, - # i.e. f.p[i,i] == 0 for all i + # i.e., P[i, i] == 0 for all i - P = f.p(uprev, p, t) # evaluate production matrix + # evaluate production matrix + P = f.p(uprev, p, t) + integrator.stats.nf += 1 - # avoid division by zero due to zero patankar weights + # avoid division by zero due to zero Patankar weights σ = add_small_constant(uprev, small_constant) # build linear system matrix @@ -303,46 +357,74 @@ function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false) alias_A = false, alias_b = false, assumptions = LinearSolve.OperatorAssumptions(true)) u = sol.u + integrator.stats.nsolve += 1 - k = f(u, p, t + dt) # For the interpolation, needs k at the updated point - integrator.stats.nf += 1 - integrator.fsallast = k - integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast integrator.u = u end function initialize!(integrator, cache::MPECache) - integrator.kshortsize = 2 @unpack k, fsalfirst = cache integrator.fsalfirst = fsalfirst integrator.fsallast = k + integrator.kshortsize = 1 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast - integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point - integrator.stats.nf += 1 end function perform_step!(integrator, cache::MPECache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator @unpack P, D, weight = cache - #@muladd @.. broadcast=false u=0*uprev + 123.0# + dt * integrator.fsalfirst - P .= 0.0 + # TODO: Shall we require the users to set unused entries to zero? + fill!(P, zero(eltype(P))) + f.p(P, uprev, p, t) # evaluate production terms sum_destruction_terms!(D, P) # store destruction terms in D + integrator.stats.nf += 1 + build_mprk_matrix!(P, 1, P, D, uprev, dt) - # linres = P\uprev # TODO: needs to be implemented without allocations + # Same as linres = P \ uprev linres = dolinsolve(integrator, cache.linsolve; A = P, b = _vec(uprev), du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight) - u .= linres + integrator.stats.nsolve += 1 +end - f(integrator.fsallast, u, p, t + dt) # For the interpolation, needs k at the updated point - integrator.stats.nf += 1 +# interpolation specializations +interp_summary(::MPE) = "Linear interpolation" + +function _ode_interpolant(Θ, dt, u0, u1, k, + cache::Union{MPEConstantCache, MPECache}, + 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::Union{MPEConstantCache, MPECache}, + 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::Union{MPEConstantCache, MPECache}, + 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::Union{MPEConstantCache, MPECache}, + idxs, # Optionally specialize for ::Nothing and others + T::Type{Val{1}}, + differential_vars::Nothing) + linear_interpolant!(out, Θ, dt, u0, u1, idxs, T) end ### MPRK ##################################################################################### @@ -407,7 +489,8 @@ function Base.getproperty(alg::MPRK22, f::Symbol) end end -alg_order(alg::MPRK22) = 2 +alg_order(::MPRK22) = 2 +isfsal(::MPRK22) = false struct MPRK22Cache{uType, rateType, PType, tabType, Thread, F, uNoUnitsType} <: OrdinaryDiffEqMutableCache @@ -482,15 +565,23 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, end function initialize!(integrator, cache::MPRK22ConstantCache) - integrator.kshortsize = 2 + integrator.kshortsize = 1 integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) - integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal - integrator.stats.nf += 1 # Avoid undefined entries if k is an array of arrays - integrator.fsallast = zero(integrator.fsalfirst) - integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast + integrator.fsalfirst = zero(integrator.u) + integrator.fsallast = integrator.fsalfirst + integrator.k[1] = integrator.fsallast + + # TODO: Do we need to set fsalfirst here? The other non-FSAL caches + # in OrdinaryDiffEq.jl use something like + # integrator.fsalfirst = integrator.f(integrator.uprev, integrator, + # integrator.t) # Pre-start fsal + # integrator.stats.nf += 1 + # integrator.fsallast = zero(integrator.fsalfirst) + # integrator.k[1] = integrator.fsalfirst + # Do we need something similar here to get a cache for k values + # with the correct units? end function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = false) @@ -498,10 +589,12 @@ function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = fal @unpack a21, b1, b2, small_constant = cache # Attention: Implementation assumes that the pds is conservative, - # i.e. f.p[i,i] == 0 for all i + # i.e. , P[i, i] == 0 for all i - P = f.p(uprev, p, t) # evaluate production matrix + # evaluate production matrix + P = f.p(uprev, p, t) Ptmp = a21 * P + integrator.stats.nf += 1 # avoid division by zero due to zero patankar weights σ = add_small_constant(uprev, small_constant) @@ -515,9 +608,10 @@ function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = fal alias_A = false, alias_b = false, assumptions = LinearSolve.OperatorAssumptions(true)) u = sol.u + integrator.stats.nsolve += 1 # compute Patankar weight denominator - if a21 == 1.0 + if isone(a21) σ = u else # σ = σ .* (u ./ σ) .^ (1 / a21) # generated Infs when solving brusselator @@ -528,6 +622,7 @@ function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = fal P2 = f.p(u, p, t + a21 * dt) Ptmp = b1 * P + b2 * P2 + integrator.stats.nf += 1 # build linear system matrix M = build_mprk_matrix(Ptmp, σ, dt) @@ -538,10 +633,7 @@ function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = fal alias_A = false, alias_b = false, assumptions = LinearSolve.OperatorAssumptions(true)) u = sol.u - - k = f(u, p, t + dt) # For the interpolation, needs k at the updated point - integrator.stats.nf += 1 - integrator.fsallast = k + integrator.stats.nsolve += 1 # copied from perform_step for HeunConstantCache # If a21 = 1.0, then σ is the MPE approximation and thus suited for stiff problems. @@ -551,21 +643,16 @@ function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = fal integrator.opts.reltol, integrator.opts.internalnorm, t) integrator.EEst = integrator.opts.internalnorm(atmp, t) - integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast integrator.u = u end function initialize!(integrator, cache::MPRK22Cache) - integrator.kshortsize = 2 @unpack k, fsalfirst = cache integrator.fsalfirst = fsalfirst integrator.fsallast = k + integrator.kshortsize = 1 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast - integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point - integrator.stats.nf += 1 end function perform_step!(integrator, cache::MPRK22Cache, repeat_step = false) @@ -577,34 +664,72 @@ function perform_step!(integrator, cache::MPRK22Cache, repeat_step = false) f.p(P, uprev, p, t) # evaluate production terms sum_destruction_terms!(D, P) # store destruction terms in D + integrator.stats.nf += 1 + build_mprk_matrix!(M, a21, P, D, uprev, dt) - # linres = M \ uprev #TODO: needs to be implemented without allocations. + # Same as linres = M \ uprev linres = dolinsolve(integrator, cache.linsolve; A = M, b = _vec(uprev), du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight) u .= linres + integrator.stats.nsolve += 1 u .= u .+ small_constant σ .= uprev .* (u ./ uprev) .^ (1 / a21) .+ small_constant f.p(P2, u, p, t + a21 * dt) # evaluate production terms + sum_destruction_terms!(D, P) # store destruction terms in D sum_destruction_terms!(D2, P2) # store destruction terms in D2 + build_mprk_matrix!(M, b1, P, D, b2, P2, D2, σ, dt) - # linres = M \ uprev #TODO: needs to be implemented without allocations. + # Same as linres = M \ uprev linres = dolinsolve(integrator, cache.linsolve; A = M, b = _vec(uprev), du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight) u .= linres + integrator.stats.nsolve += 1 tmp .= u .- σ calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t, thread) integrator.EEst = integrator.opts.internalnorm(atmp, t) +end - f(integrator.fsallast, u, p, t + dt) # For the interpolation, needs k at the updated point - integrator.stats.nf += 1 +# interpolation specializations +interp_summary(::MPRK22) = "Linear interpolation" + +function _ode_interpolant(Θ, dt, u0, u1, k, + cache::Union{MPRK22ConstantCache, MPRK22Cache}, + 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::Union{MPRK22ConstantCache, MPRK22Cache}, + 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::Union{MPRK22ConstantCache, MPRK22Cache}, + 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::Union{MPRK22ConstantCache, MPRK22Cache}, + idxs, # Optionally specialize for ::Nothing and others + T::Type{Val{1}}, + differential_vars::Nothing) + linear_interpolant!(out, Θ, dt, u0, u1, idxs, T) end diff --git a/test/Project.toml b/test/Project.toml index bc83d7a2..f8cca8d9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -11,3 +12,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Aqua = "0.7, 0.8" LinearSolve = "2" OrdinaryDiffEq = "6" +StaticArrays = "1.5" diff --git a/test/runtests.jl b/test/runtests.jl index 6afb7918..fe1b63a1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,8 @@ using LinearAlgebra using SparseArrays using Statistics: mean +using StaticArrays: MVector + using OrdinaryDiffEq using PositiveIntegrators @@ -70,6 +72,13 @@ function experimental_order_of_convergence(errors, dts) return mean(orders) end +const prob_pds_linmod_array = ConservativePDSProblem(prob_pds_linmod.f, + Array(prob_pds_linmod.u0), + prob_pds_linmod.tspan) +const prob_pds_linmod_mvector = ConservativePDSProblem(prob_pds_linmod_inplace.f, + MVector(prob_pds_linmod.u0), + prob_pds_linmod.tspan) + @testset "PositiveIntegrators.jl tests" begin @testset "Aqua.jl" begin # We do not test ambiguities since we get a lot of @@ -359,8 +368,9 @@ end @testset "Convergence tests" begin alg = MPE() - dts = 0.5 .^ (5:10) - problems = (prob_pds_linmod, prob_pds_linmod_inplace) + dts = 0.5 .^ (6:11) + problems = (prob_pds_linmod, prob_pds_linmod_array, + prob_pds_linmod_mvector, prob_pds_linmod_inplace) for prob in problems eoc = experimental_order_of_convergence(prob, alg, dts) @test isapprox(eoc, PositiveIntegrators.alg_order(alg); atol = 0.2) @@ -455,8 +465,9 @@ end end @testset "Convergence tests" begin - dts = 0.5 .^ (5:10) - problems = (prob_pds_linmod, prob_pds_linmod_inplace) + dts = 0.5 .^ (6:11) + problems = (prob_pds_linmod, prob_pds_linmod_array, + prob_pds_linmod_mvector, prob_pds_linmod_inplace) for alpha in (0.5, 1.0, 2.0), prob in problems alg = MPRK22(alpha) eoc = experimental_order_of_convergence(prob, alg, dts)