Draft: ImplicitDiscreteSolve

.github/workflows/CI_SimpleImplicitDiscreteSolve.yml

@@ -0,0 +1,120 @@
+name: CI (SimpleImplicitDiscreteSolve)
+
+on:
+  pull_request:
+    branches:
+      - master
+    paths:
+      - "lib/SimpleImplicitDiscreteSolve/**"
+      - ".github/workflows/CI_SimpleImplicitDiscreteSolve.yml"
+      - "lib/SimpleNonlinearSolve/**"
+  push:
+    branches:
+      - master
+
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - "1.10"
+          - "1"
+        os:
+          - ubuntu-latest
+          - macos-latest
+          - windows-latest
+        group:
+          - core
+          - adjoint
+          - alloc_check
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+      - uses: actions/cache@v4
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - name: "Install Dependencies and Run Tests"
+        run: |
+          import Pkg
+          Pkg.Registry.update()
+          # Install packages present in subdirectories
+          dev_pks = Pkg.PackageSpec[]
+          for path in ("lib/SimpleNonlinearSolve",)
+              push!(dev_pks, Pkg.PackageSpec(; path))
+          end
+          Pkg.develop(dev_pks)
+          Pkg.instantiate()
+          Pkg.test(; coverage="user")
+        shell: julia --color=yes --code-coverage=user --depwarn=yes --project=lib/SimpleImplicitDiscreteSolve {0}
+        env:
+          GROUP: ${{ matrix.group }}
+      - uses: julia-actions/julia-processcoverage@v1
+        with:
+          directories: lib/SimpleNonlinearSolve/src,lib/SimpleNonlinearSolve/ext,lib/SimpleImplicitDiscreteSolve/src
+      - uses: codecov/codecov-action@v5
+        with:
+          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          verbose: true
+          fail_ci_if_error: false
+
+  downgrade:
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - "1.10"
+        group:
+          - core
+          - adjoint
+          - alloc_check
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+      - uses: julia-actions/julia-downgrade-compat@v1
+        with:
+          skip: SimpleImplicitDiscreteSolve
+      - name: "Install Dependencies and Run Tests"
+        run: |
+          import Pkg
+          Pkg.Registry.update()
+          # Install packages present in subdirectories
+          dev_pks = Pkg.PackageSpec[]
+          for path in ("lib/SimpleNonlinearSolve",)
+              push!(dev_pks, Pkg.PackageSpec(; path))
+          end
+          Pkg.develop(dev_pks)
+          Pkg.instantiate()
+          Pkg.test(; coverage="user")
+        shell: julia --color=yes --code-coverage=user --depwarn=yes --project=lib/SimpleImplicitDiscreteSolve {0}
+        env:
+          GROUP: ${{ matrix.group }}
+      - uses: julia-actions/julia-processcoverage@v1
+        with:
+          directories: lib/SimpleImplicitDiscreteSolve/src,lib/SimpleNonlinearSolve/ext,lib/SimpleNonlinearSolve/src
+      - uses: codecov/codecov-action@v5
+        with:
+          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          verbose: true
+          fail_ci_if_error: false

Project.toml

@@ -25,6 +25,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SimpleImplicitDiscreteSolve = "8b67ef88-54bd-43ff-aca0-8be02588656a"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"

lib/SimpleImplicitDiscreteSolve/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2024 SciML
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

lib/SimpleImplicitDiscreteSolve/Project.toml

@@ -0,0 +1,25 @@
+name = "SimpleImplicitDiscreteSolve"
+uuid = "8b67ef88-54bd-43ff-aca0-8be02588656a"
+authors = ["vyudu <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
+
+[compat]
+DiffEqBase = "6.164.1"
+OrdinaryDiffEqSDIRK = "1.2.0"
+Reexport = "1.2.2"
+SciMLBase = "2.74.1"
+SimpleNonlinearSolve = "2.1.0"
+Test = "1.10"
+
+[extras]
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["OrdinaryDiffEqSDIRK", "Test"]

lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl

@@ -0,0 +1,95 @@
+module SimpleImplicitDiscreteSolve
+
+using SciMLBase
+using SimpleNonlinearSolve
+using Reexport
+@reexport using DiffEqBase
+
+"""
+    SimpleIDSolve()
+
+Simple solver for `ImplicitDiscreteSystems`. Uses `SimpleNewtonRaphson` to solve for the next state at every timestep.
+"""
+struct SimpleIDSolve <: SciMLBase.AbstractODEAlgorithm end
+
+function DiffEqBase.__init(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; dt = 1)
+    u0 = prob.u0
+    p = prob.p
+    f = prob.f
+    t = prob.tspan[1]
+
+    nlf = isinplace(f) ? (out, u, p) -> f(out, u, u0, p, t) : (u, p) -> f(u, u0, p, t)
+    prob = NonlinearProblem{isinplace(f)}(nlf, u0, p)
+    sol = solve(prob, SimpleNewtonRaphson())
+    sol, (sol.retcode != ReturnCode.Success)
+end
+
+function DiffEqBase.solve(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; 
+    dt = 1,
+    save_everystep = true,
+    save_start = true,
+    adaptive = false,
+    dense = false,
+    save_end = true,
+    kwargs...)
+
+    @assert !adaptive
+    @assert !dense
+    (initsol, initfail) = DiffEqBase.__init(prob, alg; dt)
+    if initfail
+        sol = DiffEqBase.build_solution(prob, alg, prob.tspan[1], u0, k = nothing, stats = nothing, calculate_error = false) 
+        return SciMLBase.solution_new_retcode(sol, ReturnCode.InitialFailure)
+    end
+
+    u0 = initsol.u
+    tspan = prob.tspan
+    f = prob.f
+    p = prob.p
+    t = tspan[1]
+    tf = prob.tspan[2]
+    ts = tspan[1]:dt:tspan[2]
+
+    if save_everystep && save_start
+        us = Vector{typeof(u0)}(undef, length(ts))
+        us[1] = u0
+    elseif save_everystep
+        us = Vector{typeof(u0)}(undef, length(ts) - 1)
+    elseif save_start
+        us = Vector{typeof(u0)}(undef, 2)
+        us[1] = u0
+    else
+        us = Vector{typeof(u0)}(undef, 1) # for interface compatibility
+    end
+
+    u = u0
+    convfail = false
+    for i in 2:length(ts)
+        uprev = u
+        t = ts[i]
+        nlf = isinplace(f) ? (out, u, p) -> f(out, u, uprev, p, t) : (u, p) -> f(u, uprev, p, t)
+        nlprob = NonlinearProblem{isinplace(f)}(nlf, uprev, p)
+        nlsol = solve(nlprob, SimpleNewtonRaphson())
+        u = nlsol.u
+        save_everystep && (us[i] = u)
+        convfail = (nlsol.retcode != ReturnCode.Success)
+
+        if convfail
+            sol = DiffEqBase.build_solution(prob, alg, ts[1:i], us[1:i], k = nothing, stats = nothing, calculate_error = false)
+            sol = SciMLBase.solution_new_retcode(sol, ReturnCode.ConvergenceFailure)
+            return sol
+        end
+    end
+
+    !save_everystep && save_end && (us[end] = u)
+    sol = DiffEqBase.build_solution(prob, alg, ts, us,
+        k = nothing, stats = nothing,
+        calculate_error = false)
+
+    DiffEqBase.has_analytic(prob.f) &&
+        DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true, dense_errors = false)
+    sol
+end
+
+export SimpleIDSolve
+
+end

lib/SimpleImplicitDiscreteSolve/test/runtests.jl

@@ -0,0 +1,71 @@
+#runtests
+using Test
+using SimpleImplicitDiscreteSolve
+using OrdinaryDiffEqSDIRK
+
+# Test implicit Euler using ImplicitDiscreteProblem
+@testset "Implicit Euler" begin
+    function lotkavolterra(u, p, t)
+        [1.5*u[1] - u[1]*u[2], -3.0*u[2] + u[1]*u[2]]
+    end
+
+    function f!(resid, u_next, u, p, t)
+        lv = lotkavolterra(u_next, p, t)
+        resid[1] = u_next[1] - u[1] - 0.01*lv[1]
+        resid[2] = u_next[2] - u[2] - 0.01*lv[2]
+        nothing
+    end
+    u0 = [1., 1.]
+    tspan = (0., 0.5)
+
+    idprob = ImplicitDiscreteProblem(f!, u0, tspan, [])
+    idsol = solve(idprob, SimpleIDSolve(), dt = 0.01)
+
+    oprob = ODEProblem(lotkavolterra, u0, tspan)
+    osol = solve(oprob, ImplicitEuler())
+
+    @test isapprox(idsol[end-1], osol[end], atol = 0.1)
+
+    ### free-fall
+    # y, dy
+    function ff(u, p, t) 
+        [u[2], -9.8]
+    end
+
+    function g!(resid, u_next, u, p, t) 
+        f = ff(u_next, p, t) 
+        resid[1] = u_next[1] - u[1] - 0.01*f[1]
+        resid[2] = u_next[2] - u[2] - 0.01*f[2]
+        nothing
+    end
+    u0 = [10., 0.]
+    tspan = (0, 0.5)
+
+    idprob = ImplicitDiscreteProblem(g!, u0, tspan, [])
+    idsol = solve(idprob, SimpleIDSolve(); dt = 0.01)
+
+    oprob = ODEProblem(ff, u0, tspan)
+    osol = solve(oprob, ImplicitEuler())
+
+    @test isapprox(idsol[end-1], osol[end], atol = 0.1)
+end
+
+@testset "Solver initializes" begin
+    function periodic!(resid, u_next, u, p, t) 
+        resid[1] = u_next[1] - u[1] - sin(t*π/4)
+        resid[2] = 16 - u_next[2]^2 - u_next[1]^2 
+    end
+
+    tsteps = 15
+    u0 = [1., 3.]
+    idprob = ImplicitDiscreteProblem(periodic!, u0, (0, tsteps), [])
+    initsol, initfail = DiffEqBase.__init(idprob, SimpleIDSolve())
+    @test initsol.u[1]^2 + initsol.u[2]^2 ≈ 16
+
+    idsol = solve(idprob, SimpleIDSolve())
+
+    for ts in 1:tsteps
+        step = idsol.u[ts]
+        @test step[1]^2 + step[2]^2 ≈ 16
+    end
+end

src/NonlinearSolve.jl

@@ -42,6 +42,7 @@ using NonlinearSolveFirstOrder: NonlinearSolveFirstOrder, GeneralizedFirstOrderA
 using NonlinearSolveQuasiNewton: NonlinearSolveQuasiNewton, QuasiNewtonAlgorithm
 using NonlinearSolveSpectralMethods: NonlinearSolveSpectralMethods, GeneralizedDFSane
 using SimpleNonlinearSolve: SimpleNonlinearSolve
+using SimpleImplicitDiscreteSolve: SimpleIDSolve
 
 const SII = SymbolicIndexingInterface
 
@@ -109,7 +110,8 @@ end
 # Rexexports
 @reexport using SciMLBase, NonlinearSolveBase, LineSearch, ADTypes
 @reexport using NonlinearSolveFirstOrder, NonlinearSolveSpectralMethods,
-                NonlinearSolveQuasiNewton, SimpleNonlinearSolve, BracketingNonlinearSolve
+                NonlinearSolveQuasiNewton, SimpleNonlinearSolve, BracketingNonlinearSolve,
+                SimpleImplicitDiscreteSolve
 @reexport using LinearSolve
 
 # Poly Algorithms