Merge remote-tracking branch 'origin/master' into jacobian-reuse-rosenbrock-w-methods

Author: ChrisRackauckas

.github/workflows/SublibraryCI.yml

@@ -33,7 +33,7 @@ jobs:
             CHANGED=$(git diff --name-only ${{ github.event.before }} ${{ github.sha }} 2>/dev/null || git diff --name-only HEAD~1 HEAD)
           else
             # For pull requests, compare against the base branch
-            git fetch origin ${{ github.event.pull_request.base.ref }} --depth=1
+            git fetch origin ${{ github.event.pull_request.base.ref }}
             CHANGED=$(git diff --name-only origin/${{ github.event.pull_request.base.ref }}...HEAD)
           fi
 

lib/OrdinaryDiffEqDifferentiation/test/modelingtoolkit/Project.toml

@@ -0,0 +1,15 @@
+[deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[sources]
+OrdinaryDiffEq = {path = "../../../.."}
+
+[compat]
+ADTypes = "1.16"
+ForwardDiff = "0.10, 1"
+ModelingToolkit = "10.10, 11"
+OrdinaryDiffEq = "6"

lib/OrdinaryDiffEqDifferentiation/test/modelingtoolkit/jacobian_tests.jl

@@ -0,0 +1,104 @@
+using OrdinaryDiffEq, ForwardDiff, Test, ADTypes
+
+function d_alembert(du, u, p, t)
+    return du[1] = p[1] - p[2] * u[1] + p[3] * t
+end
+
+function d_alembert_jac(J, u, p, t)
+    return J[1] = -p[2]
+end
+
+function d_alembert_analytic(u0, p, t::Number)
+    a, b, c = p
+    ebt = exp(b * t)
+    return @. exp(-b * t) * (-a * b + c + ebt * (a * b + c * (b * t - 1)) + b^2 * u0) / (b^2)
+end
+
+p = (1.0, 2.0, 3.0)
+u0 = [1.0]
+tspan = (0.0, 10.0)
+prob = ODEProblem(
+    ODEFunction(
+        d_alembert,
+        jac = d_alembert_jac,
+        analytic = d_alembert_analytic
+    ),
+    u0, tspan, p
+)
+
+sol = solve(prob, Tsit5(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 1.0e-7
+sol = solve(prob, Rosenbrock23(), abstol = 1.0e-8, reltol = 1.0e-8)
+@test sol.errors[:l2] < 1.0e-7
+sol = solve(prob, Rodas4(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 1.0e-7
+sol = solve(prob, Veldd4(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 1.0e-7
+sol = solve(prob, Rodas5(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 1.0e-7
+sol = solve(prob, TRBDF2(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 2.0e-6
+sol = solve(prob, Trapezoid(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 2.0e-6
+sol = solve(prob, KenCarp3(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 8.0e-4
+sol = solve(prob, KenCarp4(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 1.0e-7
+sol = solve(prob, KenCarp47(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 1.0e-7
+sol = solve(prob, KenCarp58(), abstol = 1.0e-10, reltol = 1.0e-10)
+@test sol.errors[:l2] < 1.0e-7
+
+using ModelingToolkit
+function lotka(du, u, p, t)
+    x = u[1]
+    y = u[2]
+    du[1] = p[1] * x - p[2] * x * y
+    return du[2] = -p[3] * y + p[4] * x * y
+end
+
+prob = ODEProblem(lotka, [1.0, 1.0], (0.0, 1.0), [1.5, 1.0, 3.0, 1.0])
+de = ModelingToolkit.modelingtoolkitize(prob) |> complete
+prob2 = ODEProblem(de, [], prob.tspan; jac = true)
+
+sol = solve(prob, TRBDF2())
+
+for Alg in [Rodas5, Rosenbrock23, TRBDF2, KenCarp4]
+    @test Array(
+        solve(
+            prob2,
+            Alg(),
+            tstops = sol.t,
+            adaptive = false
+        )
+    ) ≈ Array(
+        solve(
+            prob,
+            Alg(),
+            tstops = sol.t,
+            adaptive = false
+        )
+    ) atol = 1.0e-4
+end
+
+## check chunk_size handling in ForwardDiff Jacobians
+const chunksize = 1
+function rober(du, u, p, t)
+    y₁, y₂, y₃ = u
+    k₁, k₂, k₃, check = p
+    if check && eltype(u) <: ForwardDiff.Dual && ForwardDiff.npartials(u[1]) != chunksize
+        @show ForwardDiff.npartials(u[1]), chunksize
+        error("chunk_size is not as specified")
+    end
+
+    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
+    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
+    du[3] = k₂ * y₂^2
+    return nothing
+end
+prob1 = ODEProblem(rober, [1.0, 0.0, 0.0], (0.0, 1.0e5), (0.04, 3.0e7, 1.0e4, true))
+sol1 = solve(prob1, TRBDF2(autodiff = AutoForwardDiff(chunksize = chunksize)))
+prob = ODEProblem(rober, [1.0, 0.0, 0.0], (0.0, 1.0e5), (0.04, 3.0e7, 1.0e4, false))
+sol = solve(prob, TRBDF2())
+@test sol.u[end] == sol1.u[end]
+@test length(sol.t) == length(sol1.t)

lib/OrdinaryDiffEqDifferentiation/test/runtests.jl

@@ -1,17 +1,42 @@
 using SafeTestsets
+using Pkg
 
 const TEST_GROUP = get(ENV, "ODEDIFFEQ_TEST_GROUP", "ALL")
 
+function activate_sparse_env()
+    Pkg.activate(joinpath(@__DIR__, "sparse"))
+    Pkg.develop(PackageSpec(path = dirname(dirname(dirname(@__DIR__)))))
+    return Pkg.instantiate()
+end
+
+function activate_modelingtoolkit_env()
+    Pkg.activate(joinpath(@__DIR__, "modelingtoolkit"))
+    Pkg.develop(PackageSpec(path = dirname(dirname(dirname(@__DIR__)))))
+    return Pkg.instantiate()
+end
+
 # Run QA tests (JET, Aqua)
-if TEST_GROUP != "Core" && isempty(VERSION.prerelease)
+if TEST_GROUP ∉ ("Core", "Sparse", "ModelingToolkit") && isempty(VERSION.prerelease)
     @time @safetestset "JET Tests" include("jet.jl")
     @time @safetestset "Aqua" include("qa.jl")
 end
 
 # Run functional tests
-if TEST_GROUP != "QA"
+if TEST_GROUP ∉ ("QA", "Sparse", "ModelingToolkit")
     @time @safetestset "OOP J_t Tracking" include("oop_jt_tracking_test.jl")
     @time @safetestset "Differentiation Trait Tests" include("differentiation_traits_tests.jl")
     @time @safetestset "Autodiff Error Tests" include("autodiff_error_tests.jl")
     @time @safetestset "No Jac Tests" include("nojac_tests.jl")
 end
+
+# Run sparse tests (separate environment due to ComponentArrays dep conflicts)
+if TEST_GROUP == "Sparse" || TEST_GROUP == "ALL"
+    activate_sparse_env()
+    @time @safetestset "Non-Full Diagonal Sparsity Tests" include("sparse/nonfulldiagonal_sparse_tests.jl")
+end
+
+# Run ModelingToolkit tests (separate environment due to heavy MTK dependency)
+if TEST_GROUP == "ModelingToolkit" && isempty(VERSION.prerelease)
+    activate_modelingtoolkit_env()
+    @time @safetestset "Jacobian Tests" include("modelingtoolkit/jacobian_tests.jl")
+end

lib/OrdinaryDiffEqDifferentiation/test/sparse/Project.toml

@@ -0,0 +1,15 @@
+[deps]
+ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[sources]
+OrdinaryDiffEq = {path = "../../../.."}
+
+[compat]
+ComponentArrays = "0.15, 1"
+LinearSolve = "3"
+OrdinaryDiffEq = "6"

lib/OrdinaryDiffEqDifferentiation/test/sparse/nonfulldiagonal_sparse_tests.jl

@@ -0,0 +1,154 @@
+using OrdinaryDiffEq, SparseArrays, LinearSolve, LinearAlgebra
+using ComponentArrays
+
+function enclosethetimedifferential(parameters::NamedTuple)::Function
+    @info "Enclosing the time differential"
+
+    (; Δr, r_space, countorderapprox) = parameters.compute
+    N = length(r_space)
+
+    function first_deriv(N)
+        dx = 1 / (N + 1)
+        du = -1 * ones(N - 1) # off diagonal
+        du2 = ones(N - 1) # off diagonal
+        diag = zeros(N)
+        lower = spzeros(Float64, N)
+        upper = spzeros(Float64, N)
+        lower[1] = -1.0
+        upper[end] = 1.0
+        M = hcat(lower, sparse(diagm(-1 => du, 0 => diag, 1 => du2)), upper)
+        return MatrixOperator(1 / dx * M)
+    end
+
+    function second_deriv(N)
+        dx = 1 / (N + 1)
+        du = ones(N - 1) # off diagonal
+        du2 = ones(N - 1) # off diagonal
+        diag = -2 * ones(N)
+        lower = spzeros(Float64, N)
+        upper = spzeros(Float64, N)
+        lower[1] = 1.0
+        upper[end] = 1.0
+        M = hcat(lower, sparse(diagm(-1 => du, 0 => diag, 1 => du2)), upper)
+        return MatrixOperator(1 / dx^2 * M)
+    end
+
+    function extender(N)
+        dx = 1 / (N + 1)
+        diag = ones(N)
+        lower = spzeros(Float64, N)
+        upper = spzeros(Float64, N)
+        lower[1] = 1.0
+        upper[end] = 1.0
+        M = vcat(
+            transpose(lower),
+            sparse(diagm(diag)),
+            transpose(upper)
+        )
+        return MatrixOperator(1 / dx^2 * M)
+    end
+
+    bc_handler = extender(N)
+
+    ∇ = first_deriv(N) * bc_handler
+    Δ = second_deriv(N) * bc_handler
+
+    bc_x = zeros(Real, N)
+    bc_xx = zeros(Real, N)
+
+    function timedifferentialclosure!(du, u, p, t)
+        (;
+            α, D, v, k_p, V_c, Q_l, Q_r, V_b,
+            S, Lm, Dm, V_v,
+        ) = p
+
+        c = u[1:(end - 3)]
+        c_v = u[end - 2]
+        c_c = u[end - 1]
+        c_b = u[end]
+
+        J_B0 = (Dm / Lm) * (α * c_v - c[1])
+        J_BL = (Dm / Lm) * (c[end] - α * c_c)
+        grad_0 = (v ./ D) .* c[1] .- J_B0 ./ D
+        grad_L = (v ./ D) .* c[end] .- J_BL ./ D
+
+        bc_x[1] = grad_0 / 2
+        bc_x[end] = grad_L / 2
+        grad_c = ∇ * c + bc_x
+
+        bc_xx[1] = -grad_0 / Δr
+        bc_xx[end] = grad_L / Δr
+        Lap_c = Δ * c + bc_xx
+
+        C = sum(Δr .* S * (k_p * (c .- c_b)))
+
+        dc_dt = D * Lap_c - v * grad_c .- k_p * (c .- c_b)
+        du[1:(end - 3)] = dc_dt[1:end]
+
+        dcv_dt = -S * J_B0 / V_v - (Q_l / V_v) * c_v
+        du[end - 2] = dcv_dt
+
+        dcc_dt = S * α * J_BL / V_c + (Q_l / V_c) * c_v - (Q_l / V_c) * c_c
+        du[end - 1] = dcc_dt
+
+        dcb_dt = (Q_l / V_b) * c_c + C / V_b
+        return du[end] = dcb_dt
+    end
+
+    return timedifferentialclosure!
+end
+
+prior = ComponentArray(;
+    α = 0.2,
+    D = 0.46,
+    v = 0.0,
+    k_p = 0.0,
+    V_c = 18,
+    Q_l = 20,
+    Q_r = 3.6,
+    V_b = 1490,
+    S = 52,
+    Lm = 0.05,
+    Dm = 0.046,
+    V_v = 18.0
+)
+
+r_space = collect(range(0.0, 2.0, length = 15))
+computeparams = (
+    Δr = r_space[2],
+    r_space = r_space,
+    countorderapprox = 2,
+)
+parameters = (
+    prior = prior,
+    compute = computeparams,
+)
+
+dudt = enclosethetimedifferential(parameters)
+IC = ones(length(r_space) + 3)
+odeprob = ODEProblem(
+    dudt,
+    IC,
+    (0, 2.1),
+    parameters.prior
+);
+du0 = copy(odeprob.u0);
+# Hardcoded sparsity pattern for 15 spatial points + 3 state variables (18x18 matrix)
+I = [1, 2, 16, 18, 1, 2, 3, 18, 2, 3, 4, 18, 3, 4, 5, 18, 4, 5, 6, 18, 5, 6, 7, 18, 6, 7, 8, 18, 7, 8, 9, 18, 8, 9, 10, 18, 9, 10, 11, 18, 10, 11, 12, 18, 11, 12, 13, 18, 12, 13, 14, 18, 13, 14, 15, 18, 14, 15, 17, 18, 1, 16, 17, 15, 17, 18, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 18]
+J = [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18]
+jac_sparsity = sparse(I, J, ones(Bool, length(I)), 18, 18);
+f = ODEFunction(
+    dudt;
+    jac_prototype = float.(jac_sparsity)
+);
+sparseodeprob = ODEProblem(
+    f,
+    odeprob.u0,
+    (0, 2.1),
+    parameters.prior
+);
+
+solve(odeprob, TRBDF2());
+solve(sparseodeprob, TRBDF2());
+solve(sparseodeprob, Rosenbrock23(linsolve = KLUFactorization()));
+solve(sparseodeprob, KenCarp47(linsolve = KrylovJL_GMRES()));

lib/OrdinaryDiffEqDifferentiation/test/test_groups.toml

@@ -0,0 +1,11 @@
+[Core]
+versions = ["lts", "1.11", "1", "pre"]
+
+[QA]
+versions = ["1"]
+
+[Sparse]
+versions = ["1"]
+
+[ModelingToolkit]
+versions = ["1"]

lib/OrdinaryDiffEqNonlinearSolve/test/linear_solver_split_ode_tests.jl

@@ -1,7 +1,6 @@
 using Test
 using OrdinaryDiffEqBDF, OrdinaryDiffEqSDIRK
 using LinearAlgebra, LinearSolve
-using OrdinaryDiffEqCore: dolinsolve
 
 n = 8
 dt = 1 / 1000

lib/OrdinaryDiffEqNonlinearSolve/test/modelingtoolkit/Project.toml

@@ -0,0 +1,29 @@
+[deps]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
+OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[sources]
+OrdinaryDiffEq = {path = "../../../.."}
+OrdinaryDiffEqBDF = {path = "../../../OrdinaryDiffEqBDF"}
+OrdinaryDiffEqNonlinearSolve = {path = "../.."}
+OrdinaryDiffEqSDIRK = {path = "../../../OrdinaryDiffEqSDIRK"}
+
+[compat]
+DiffEqDevTools = "2"
+IncompleteLU = "0.2"
+LinearSolve = "3"
+ModelingToolkit = "10.10, 11"
+NonlinearSolve = "4.10"
+OrdinaryDiffEq = "6"
+OrdinaryDiffEqBDF = "1"
+OrdinaryDiffEqNonlinearSolve = "1"
+OrdinaryDiffEqSDIRK = "1"