Sandu projection

Project.toml

@@ -18,8 +18,15 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
+[weakdeps]
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+
+[extensions]
+JuMPExt = "JuMP"
+
 [compat]
 FastBroadcast = "0.3.5"
+JuMP = "1.28"
 LinearAlgebra = "1"
 LinearSolve = "3.7.1"
 MuladdMacro = "0.2.4"

docs/Project.toml

@@ -1,11 +1,13 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEqFIRK = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
@@ -22,12 +24,14 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 BenchmarkTools = "1"
+Clarabel = "0.11"
 DiffEqBase = "6.160"
 DiffEqCallbacks = "4"
 DiffEqDevTools = "2.45.1"
 Documenter = "1"
 DoubleFloats = "1.4"
 InteractiveUtils = "1"
+JuMP = "1.28"
 LinearAlgebra = "1"
 LinearSolve = "3.7.1"
 OrdinaryDiffEqFIRK = "1.7"

docs/make.jl

@@ -9,6 +9,7 @@ if (get(ENV, "CI", nothing) != "true") &&
 end
 
 using PositiveIntegrators
+using JuMP # load JuMPExt
 
 # Define module-wide setups such that the respective modules are available in doctests
 DocMeta.setdocmeta!(PositiveIntegrators,
@@ -68,7 +69,7 @@ EditURL = "https://github.com/NumericalMathematics/PositiveIntegrators.jl/blob/m
 end
 
 # Make documentation
-makedocs(modules = [PositiveIntegrators],
+makedocs(modules = [PositiveIntegrators, Base.get_extension(PositiveIntegrators, :JuMPExt)],
          sitename = "PositiveIntegrators.jl",
          format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                                   canonical = "https://NumericalMathematics.github.io/PositiveIntegrators.jl/stable"),
@@ -82,7 +83,8 @@ makedocs(modules = [PositiveIntegrators],
                  "Linear Advection" => "linear_advection.md",
                  "Heat Equation, Neumann BCs" => "heat_equation_neumann.md",
                  "Heat Equation, Dirichlet BCs" => "heat_equation_dirichlet.md",
-                 "Scalar equation" => "scalar_pds.md"
+                 "Scalar equation" => "scalar_pds.md",
+                 "Sandu Projection" => "sandu_projection.md"
              ],
              "Benchmarks" => [
                  "Experimental order of convergence" => "convergence.md",

docs/src/api_reference.md

@@ -24,6 +24,7 @@ prob_pds_npzd
 prob_pds_robertson
 prob_pds_sir
 prob_pds_stratreac
+prob_ode_stratreac_scaled
 ```
 
 ## Algorithms
@@ -38,6 +39,12 @@ SSPMPRK43
 MPDeC
 ```
 
+## Callbacks
+
+```@docs
+SanduProjection
+```
+
 ## Auxiliary functions
 
 ```@docs
@@ -51,4 +58,5 @@ work_precision_adaptive
 work_precision_adaptive!
 work_precision_fixed
 work_precision_fixed!
+get_numsteps_SanduProjection
 ```

docs/src/sandu_projection.md

@@ -0,0 +1,78 @@
+# [Tutorial: Positive-projection method](@id tutorial-sandu)
+
+This tutorial is about solving an ODE using the projection method introduced by Adrian Sandu in [Positive Numerical Integration Methods for Chemical Kinetic Systems](https://doi.org/10.1006%2Fjcph.2001.6750). It guarantees positivity by solving an optimization problem while preserving all linear invariants.
+
+The Sandu projection is a post-processing technique that can be used in combination with any ODE solver.
+If the ODE solver computes a negative approximation at any time step, the projection method calculates a positive approximation, also taking into account the linear invariants.
+
+## Solution of the ODE system
+
+As an example we want to the solve the NPZD problem [`prob_pds_npzd`](@ref), which is an ODE system in which negative approximations quickly lead to unacceptable solutions. First, we solve the problem without Sandu projection and select `ROS2` form [`OrdinaryDiffEq.jl`](https://docs.sciml.ai/OrdinaryDiffEq/stable/) as ODE solver.
+
+```@example Sandu_NPZD
+using PositiveIntegrators
+using OrdinaryDiffEqRosenbrock
+using Plots
+
+prob = prob_pds_npzd
+
+ref_sol = solve(prob, ROS2(); abstol = 1e-8, reltol = 1e-6); # reference solution for plotting
+
+sol = solve(prob, ROS2(); abstol = 5e-2, reltol = 1e-1)
+
+plot(ref_sol, linestyle = :dash, label = "", color = palette(:default)[1:4]')
+plot!(sol, ylims = (-2.5, 12.5), denseplot = false,  markers = :circle, linewidth = 2, color = palette(:default)[1:4]', label = ["N" "P" "Z" "D"], legend = :right)
+```
+
+The plot shows the numerical solution obtained with `ROS2` compared to a reference solution (dashed lines).
+We see that the `ROS2` method produces negative approximations, which can occur because Rosenbrock methods are not positivity-preserving. For the NPZD problem, however, this is fatal and leads to a completely unacceptable numerical solution. It is therefore particularly important to use techniques that guarantee positivity of the numerical approximations for this problem. We achieve this below with the [`SanduProjection`](@ref).
+
+To apply the [`SanduProjection`](@ref) we need to choose an [optimization solver](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers) which is supported by [JuMP.jl](https://jump.dev/JuMP.jl/stable/) and can handle quadratic optimization problems (QP). In this tutorial we select [Clarabel.jl](https://clarabel.org/stable/) as optimization solver.
+
+In addition, we need to specify the linear invariants of the problem. 
+The only linear invariant of the NPZD problem is ``N(t)+P(t)+Z(t)+D(t)=N(0)+P(0)+Z(0)+D(0)=15`` for all times ``t≥0``.
+This can be written in the form 
+```math
+\mathbf{A}^T \begin{pmatrix} N(t)\\ P(t)\\ Z(t)\\ D(t) \end{pmatrix} = \mathbf{b}
+```
+with ``\mathbf{A}^T = [1.0,\  1.0,\  1.0,\  1.0]`` and ``\mathbf{b} = [15]``.
+
+The projection method [`SanduProjection`](@ref) is implemented as a callback and hence, must be passed as an argument to the keyword `callback`. In addition, we must also use `save_everystep = false`.
+
+```@example Sandu_NPZD
+using JuMP, Clarabel
+
+AT = [1.0 1.0 1.0 1.0]
+b = [15.0]
+proj = SanduProjection(Model(Clarabel.Optimizer), AT, b)
+
+sol_proj = solve(prob, ROS2(); abstol = 5e-2, reltol = 1e-1,
+                 save_everystep = false, callback = proj);
+
+plot(ref_sol, linestyle = :dash, label = "", color = palette(:default)[1:4]')
+plot!(sol_proj, ylims = (-2.5, 12.5), denseplot = false,  markers = :circle, linewidth = 2, color = palette(:default)[1:4]', label = ["N" "P" "Z" "D"], legend = :right)            
+```
+
+As intended, negative approximations no longer occur and we obtain an acceptable approximation.
+
+The [`SanduProjection`](@ref) is implemented as a [`DiscreteCallback`](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#SciMLBase.DiscreteCallback) and we can display the number of projection steps in the following way.
+
+```@example Sandu_NPZD
+@show get_numsteps_SanduProjection(proj) 
+```
+
+We can see that in this example, a single projection step was already sufficient.
+
+## Package versions
+
+These results were obtained using the following versions.
+```@example NPZD
+using InteractiveUtils
+versioninfo()
+println()
+
+using Pkg
+Pkg.status(["PositiveIntegrators", "JuMP", "Clarabel", "OrdinaryDiffEqRosenbrock", "Plots"],
+           mode=PKGMODE_MANIFEST)
+nothing # hide
+```

ext/JuMPExt.jl

@@ -0,0 +1,142 @@
+module JuMPExt
+
+using StaticArrays: StaticArray, SVector
+using JuMP: @variable, @objective, @constraint, print, set_silent,
+            optimize!, is_solved_and_feasible, value, set_string_names_on_creation,
+            set_objective_coefficient, @expression
+using SciMLBase: DiscreteCallback
+using PositiveIntegrators
+
+mutable struct SanduProjection{M} <: PositiveIntegrators.SanduProjection
+    model::M
+    cnt::Int
+end
+
+"""
+    SanduProjection(model, AT, b, eps = nothing; [save = true, verbose = false])
+
+A projection method which ensures conservation of prescribed linear invariants and positivity.
+If the current approximation ``\\mathbf{u}`` has negative components then a projection ``\\mathbf{z}`` is computed such that
+```math
+\\min \\lVert \\mathbf{z} - \\mathbf{u} \\rVert_{\\mathbf{G}},\\quad \\mathbf{A}^T\\mathbf{z}=\\mathbf{b},\\quad \\mathbf{z}≥ \\mathbf{0},
+```
+is satisfied, where the matrix ``\\mathbf{A^T}`` and the vector ``\\mathbf{b}`` define the linear invariants.
+The ``G``-norm is defined by ``\\lVert \\mathbf{u}\\rVert_{\\mathbf{G}} = \\sqrt{\\mathbf{u}^T\\mathbf{G}\\mathbf{u}}``,
+and assuming ``\\mathbf{u} = (u_1,\\dots, u_s)^T``,
+the positive definite diagonal matrix ``\\mathbf{G}`` is given by
+```math
+\\mathbf{G(\\mathbf{u})}=\\operatorname{diag}\\biggl(\\frac{1}{s(\\mathtt{abstol}+\\mathtt{reltol}\\lvert u_i \\rvert^2)}\\biggr),
+``` 
+where `abstol` and `reltol` denote the absolute and relative tolerances of the adaptive step size control, respectively.
+See Sandu (2001) for details.
+
+To use this callback one must also specify `save_everystep = false`.
+
+## Arguments
+
+  - `model`: A [`JuMP Model`](https://jump.dev/JuMP.jl/stable/api/JuMP/#JuMP.Model) to solve the minimization problem.
+  - `AT`: The matrix ``\\mathbf{A}^T`` defining the linear invariants.
+  - `b`: The vector ``\\mathbf{b}`` defining the linear invariants.
+  - `eps`: It may be helpful for the optimization solver that feasible solutions are bounded away from 0. To achieve this one can specify the optional parameter `eps`. The positivity constraint is then replaced by ``\\mathbf{z}≥```eps`, where `eps` can either be a scalar or a vector. 
+
+## Keyword Arguments
+
+  - `save`: If the keyword argument `save` is set to `false` only the initial value and the last approximation will be saved.
+            The default value is `true`.
+  - `verbose`: Enables additional output of the optimization solver. The default value is `false`.
+
+## References
+
+- Adrian Sandu.
+  "Positive numerical integration methods for chemical kinetic systems."
+  Journal of Computational Physics 170 (2001): 589-602.
+  [DOI: 10.1006/jcph.2001.6750](https://doi.org/10.1006/jcph.2001.6750)
+"""
+function PositiveIntegrators.SanduProjection(args...; kwargs...)
+    SanduProjection(args...; kwargs...)
+end
+
+function SanduProjection(model, AT, b, eps = nothing; save = true, verbose = false)
+    if isnothing(eps) || eps isa Number
+        epsv = zeros(eltype(AT), size(AT, 2))
+        if eps isa Number
+            fill!(epsv, eps)
+        end
+    else
+        epsv = eps
+    end
+
+    # Set up optimization problem
+    s = size(AT, 2)
+    if !verbose
+        set_silent(model)
+    end
+    set_string_names_on_creation(model, false)
+
+    @variable(model, z[i = 1:s]>=epsv[i])
+    @constraint(model, AT * z.==b)
+    # This just initializes the objective. The correct coefficients will be set later
+    @expression(model, obj_exp, sum(z .^ 2)+sum(z))
+    @objective(model, Min, obj_exp)
+
+    if verbose
+        print(model)
+    end
+
+    affect! = SanduProjection(model, 0)
+
+    return DiscreteCallback(Returns(true), affect!; save_positions = (false, save),
+                            initialize = initialize_sandu_projection!)
+end
+
+function initialize_sandu_projection!(c, u, t, integrator)
+    return initialize_sandu_projection!(c.affect!)
+end
+
+function initialize_sandu_projection!(proj::SanduProjection)
+    proj.cnt = 0
+end
+
+function (proj::SanduProjection)(integrator)
+    u = integrator.u
+
+    if isnegative(u)
+        proj.cnt += 1
+
+        rtol = integrator.opts.reltol
+        atol = integrator.opts.abstol
+        model = proj.model
+
+        s = length(u)
+        g = @. 1 / (s * (atol + rtol * abs(u))^2)
+
+        # set coefficients of quadratic terms
+        set_objective_coefficient(model, model[:z], model[:z], Vector(1 / 2 .* g))
+        # set coefficients of linear terms
+        set_objective_coefficient(model, model[:z], Vector(-g .* u))
+
+        # solve optimization problem
+        optimize!(model)
+        if !is_solved_and_feasible(model)
+            error("Solver did not find an optimal solution")
+        end
+
+        if integrator.u isa StaticArray
+            integrator.u = SVector{length(integrator.u)}(value.(model[:z]))
+        else
+            integrator.u = value.(model[:z])
+        end
+    end
+    return nothing
+end
+
+"""
+    get_numsteps_SanduProjection(proj)
+
+For a `SanduProjection` `proj`, this function returns the number of the performed projection steps.
+"""
+function PositiveIntegrators.get_numsteps_SanduProjection(proj)
+    return proj.affect!.cnt
+end
+
+end

src/PDSProblemLibrary.jl

@@ -450,6 +450,71 @@ prob_pds_stratreac = PDSProblem(P_stratreac, d_stratreac, u0_stratreac, (4.32e4,
                                 linear_invariants = @SMatrix[1.0 1.0 3.0 2.0 1.0 2.0;
                                                              0.0 0.0 0.0 0.0 1.0 1.0])
 
+function f_stratreac_scaled(u, p, t)
+    uc = [9.906e1, 6.624e8, 5.326e11, 1.697e16, 4e6, 1.093e9]
+
+    Tr = 4.5
+    Ts = 19.5
+    T = mod(t / 3600, 24)
+    if (Tr <= T) && (T <= Ts)
+        Tfrac = (2 * T - Tr - Ts) / (Ts - Tr)
+        sigma = 0.5 + 0.5 * cos(pi * abs(Tfrac) * Tfrac)
+    else
+        sigma = zero(t)
+    end
+
+    M = 8.120e16
+
+    k1 = 2.643e-10 * sigma^3
+    k2 = 8.018e-17
+    k3 = 6.120e-4 * sigma
+    k4 = 1.567e-15
+    k5 = 1.070e-3 * sigma^2
+    k6 = 7.110e-11
+    k7 = 1.200e-10
+    k8 = 6.062e-15
+    k9 = 1.069e-11
+    k10 = 1.289e-2 * sigma
+    k11 = 1.0e-8
+
+    r1 = k1 * u[4] * uc[4]
+    r2 = k2 * u[2] * uc[2] * u[4] * uc[4]
+    r3 = k3 * u[3] * uc[3]
+    r4 = k4 * u[3] * uc[3] * u[2] * uc[2]
+    r5 = k5 * u[3] * uc[3]
+    r6 = k6 * M * u[1] * uc[1]
+    r7 = k7 * u[1] * uc[1] * u[3] * uc[3]
+    r8 = k8 * u[3] * uc[3] * u[5] * uc[5]
+    r9 = k9 * u[6] * uc[6] * u[2] * uc[2]
+    r10 = k10 * u[6] * uc[6]
+    r11 = k11 * u[5] * uc[5] * u[2] * uc[2]
+
+    return @SVector [(r5 - r6 - r7) / uc[1];
+                     (2 * r1 - r2 + r3 - r4 + r6 - r9 + r10 - r11) / uc[2];
+                     (r2 - r3 - r4 - r5 - r7 - r8) / uc[3];
+                     (-r1 - r2 + r3 + 2 * r4 + r5 + 2 * r7 + r8 + r9) / uc[4];
+                     (-r8 + r9 + r10 - r11) / uc[5];
+                     (r8 - r9 - r10 + r11) / uc[6]]
+end
+u0_stratreac_scaled = @SVector ones(6)
+"""
+    prob_ode_stratreac_scaled
+
+Scaled version of the stratosperic reaction problem [`prob_pds_stratreac`](@ref).
+Each component is scaled by its corresponding original initial value.  
+
+The initial value is ``\\mathbf{u}_0 = (1,1,1,1,1,1)^T`` and the time domain ``(4.32⋅10^{4}, 3.024⋅10^5)``.
+
+There are two independent linear invariants. The function `linear_invariants_stratreac_scaled` returns the invariance matrix.
+"""
+prob_ode_stratreac_scaled = ODEProblem(f_stratreac_scaled, u0_stratreac_scaled,
+                                       (4.32e4, 3.024e5))
+
+function linear_invariants_stratreac_scaled()
+    return @SMatrix [99.06 6.624e8 1.5978e12 3.394e16 4.0e6 2.186e9;
+                     0.0 0.0 0.0 0.0 4.0e6 1.093e9]
+end
+
 # mapk problem
 function f_minmapk(u, p, t)
     k1 = 100 / 3