Split MIRK methods into subpackage

Project.toml

@@ -52,6 +52,7 @@ Logging = "1.10"
 NonlinearSolve = "3.8.1"
 ODEInterface = "0.5"
 OrdinaryDiffEq = "6.89.0"
+Pkg = "1.10.0"
 PreallocationTools = "0.4.24"
 PrecompileTools = "1.2"
 Preferences = "1.4"
@@ -74,11 +75,12 @@ JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "DiffEqDevTools", "JET", "LinearSolve", "ODEInterface", "OrdinaryDiffEq", "Random", "ReTestItems", "RecursiveArrayTools", "StaticArrays", "Test"]
+test = ["Aqua", "DiffEqDevTools", "JET", "LinearSolve", "ODEInterface", "OrdinaryDiffEq", "Pkg", "Random", "ReTestItems", "RecursiveArrayTools", "StaticArrays", "Test"]

lib/BoundaryValueDiffEqMIRK/LICENSE.md

@@ -0,0 +1,21 @@
+The BoundaryValueDiffEq.jl package is licensed under the MIT "Expat" License:
+
+> Copyright (c) 2017: ChrisRackauckas.
+> 
+> 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/BoundaryValueDiffEqMIRK/Project.toml

@@ -0,0 +1,77 @@
+name = "BoundaryValueDiffEqMIRK"
+uuid = "1a22d4ce-7765-49ea-b6f2-13c8438986a6"
+version = "0.1.0"
+
+[deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
+ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
+FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
+
+[compat]
+ADTypes = "1.2"
+Adapt = "4"
+Aqua = "0.8.7"
+ArrayInterface = "7.7"
+BandedMatrices = "1.4"
+BoundaryValueDiffEq = "5.10"
+ConcreteStructs = "0.2.3"
+DiffEqBase = "6.146"
+DiffEqDevTools = "2.44"
+FastAlmostBandedMatrices = "0.1.1"
+FastClosures = "0.3"
+ForwardDiff = "0.10.36"
+JET = "0.8"
+LinearAlgebra = "1.10"
+LinearSolve = "2.21"
+Logging = "1.10"
+NonlinearSolve = "3.8.1"
+OrdinaryDiffEq = "6.89.0"
+PreallocationTools = "0.4.24"
+PrecompileTools = "1.2"
+Preferences = "1.4"
+Random = "1.10"
+ReTestItems = "1.23.1"
+RecursiveArrayTools = "3.27.0"
+Reexport = "1.2"
+SciMLBase = "2.40"
+Setfield = "1"
+SparseArrays = "1.10"
+SparseDiffTools = "2.14"
+StaticArrays = "1.8.1"
+Test = "1.10"
+julia = "1.10"
+
+[extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Aqua", "DiffEqDevTools", "JET", "LinearSolve", "OrdinaryDiffEq", "Random", "ReTestItems", "RecursiveArrayTools", "StaticArrays", "Test"]

lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl

@@ -0,0 +1,160 @@
+module BoundaryValueDiffEqMIRK
+
+import PrecompileTools: @compile_workload, @setup_workload
+
+using ADTypes, Adapt, ArrayInterface, DiffEqBase, ForwardDiff, LinearAlgebra,
+      NonlinearSolve, Preferences, RecursiveArrayTools, Reexport, SciMLBase, Setfield,
+      SparseDiffTools
+
+using PreallocationTools: PreallocationTools, DiffCache
+
+# Special Matrix Types
+using BandedMatrices, FastAlmostBandedMatrices, SparseArrays
+
+import BoundaryValueDiffEq: BoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
+                            recursive_flatten, recursive_flatten!, recursive_unflatten!,
+                            __concrete_nonlinearsolve_algorithm, diff!,
+                            __FastShortcutBVPCompatibleNonlinearPolyalg,
+                            __FastShortcutBVPCompatibleNLLSPolyalg,
+                            concrete_jacobian_algorithm, eval_bc_residual,
+                            eval_bc_residual!, get_tmp, __maybe_matmul!, __append_similar!,
+                            __extract_problem_details, __initial_guess,
+                            __maybe_allocate_diffcache, __get_bcresid_prototype, __similar,
+                            __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
+                            recursive_flatten_twopoint!, __unsafe_nonlinearfunction,
+                            __internal_nlsolve_problem,
+                            __generate_sparse_jacobian_prototype, __extract_mesh,
+                            __extract_u0, __has_initial_guess, __initial_guess_length,
+                            __initial_guess_on_mesh, __flatten_initial_guess,
+                            __build_solution, __Fix3, __sparse_jacobian_cache,
+                            __sparsity_detection_alg, _sparse_like, ColoredMatrix
+
+import ADTypes: AbstractADType
+import ArrayInterface: matrix_colors, parameterless_type, undefmatrix, fast_scalar_indexing
+import ConcreteStructs: @concrete
+import DiffEqBase: solve
+import FastClosures: @closure
+import ForwardDiff: ForwardDiff, pickchunksize
+import Logging
+import RecursiveArrayTools: ArrayPartition, DiffEqArray
+import SciMLBase: AbstractDiffEqInterpolation, StandardBVProblem, __solve, _unwrap_val
+
+@reexport using ADTypes, DiffEqBase, NonlinearSolve, SparseDiffTools, SciMLBase
+
+include("types.jl")
+include("algorithms.jl")
+include("mirk.jl")
+include("adaptivity.jl")
+include("alg_utils.jl")
+include("collocation.jl")
+include("interpolation.jl")
+include("mirk_tableaus.jl")
+include("sparse_jacobians.jl")
+
+@setup_workload begin
+    function f1!(du, u, p, t)
+        du[1] = u[2]
+        du[2] = 0
+    end
+    f1 = (u, p, t) -> [u[2], 0]
+
+    function bc1!(residual, u, p, t)
+        residual[1] = u[:, 1][1] - 5
+        residual[2] = u[:, end][1]
+    end
+
+    bc1 = (u, p, t) -> [u[:, 1][1] - 5, u[:, end][1]]
+
+    bc1_a! = (residual, ua, p) -> (residual[1] = ua[1] - 5)
+    bc1_b! = (residual, ub, p) -> (residual[1] = ub[1])
+
+    bc1_a = (ua, p) -> [ua[1] - 5]
+    bc1_b = (ub, p) -> [ub[1]]
+
+    tspan = (0.0, 5.0)
+    u0 = [5.0, -3.5]
+    bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1))
+
+    probs = [BVProblem(f1!, bc1!, u0, tspan; nlls = Val(false)),
+        BVProblem(f1, bc1, u0, tspan; nlls = Val(false)),
+        TwoPointBVProblem(
+            f1!, (bc1_a!, bc1_b!), u0, tspan; bcresid_prototype, nlls = Val(false)),
+        TwoPointBVProblem(
+            f1, (bc1_a, bc1_b), u0, tspan; bcresid_prototype, nlls = Val(false))]
+
+    algs = []
+
+    jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))
+
+    if Preferences.@load_preference("PrecompileMIRK", true)
+        append!(algs, [MIRK2(; jac_alg), MIRK4(; jac_alg), MIRK6(; jac_alg)])
+    end
+
+    @compile_workload begin
+        @sync for prob in probs, alg in algs
+            Threads.@spawn solve(prob, alg; dt = 0.2)
+        end
+    end
+
+    f1_nlls! = (du, u, p, t) -> begin
+        du[1] = u[2]
+        du[2] = -u[1]
+    end
+
+    f1_nlls = (u, p, t) -> [u[2], -u[1]]
+
+    bc1_nlls! = (resid, sol, p, t) -> begin
+        solₜ₁ = sol[:, 1]
+        solₜ₂ = sol[:, end]
+        resid[1] = solₜ₁[1]
+        resid[2] = solₜ₂[1] - 1
+        resid[3] = solₜ₂[2] + 1.729109
+        return nothing
+    end
+    bc1_nlls = (sol, p, t) -> [sol[:, 1][1], sol[:, end][1] - 1, sol[:, end][2] + 1.729109]
+
+    bc1_nlls_a! = (resid, ua, p) -> (resid[1] = ua[1])
+    bc1_nlls_b! = (resid, ub, p) -> (resid[1] = ub[1] - 1;
+    resid[2] = ub[2] + 1.729109)
+
+    bc1_nlls_a = (ua, p) -> [ua[1]]
+    bc1_nlls_b = (ub, p) -> [ub[1] - 1, ub[2] + 1.729109]
+
+    tspan = (0.0, 100.0)
+    u0 = [0.0, 1.0]
+    bcresid_prototype1 = Array{Float64}(undef, 3)
+    bcresid_prototype2 = (Array{Float64}(undef, 1), Array{Float64}(undef, 2))
+
+    probs = [
+        BVProblem(BVPFunction(f1_nlls!, bc1_nlls!; bcresid_prototype = bcresid_prototype1),
+            u0, tspan, nlls = Val(true)),
+        BVProblem(BVPFunction(f1_nlls, bc1_nlls; bcresid_prototype = bcresid_prototype1),
+            u0, tspan, nlls = Val(true)),
+        TwoPointBVProblem(f1_nlls!, (bc1_nlls_a!, bc1_nlls_b!), u0, tspan;
+            bcresid_prototype = bcresid_prototype2, nlls = Val(true)),
+        TwoPointBVProblem(f1_nlls, (bc1_nlls_a, bc1_nlls_b), u0, tspan;
+            bcresid_prototype = bcresid_prototype2, nlls = Val(true))]
+
+    jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))
+
+    nlsolvers = [LevenbergMarquardt(; disable_geodesic = Val(true)), GaussNewton()]
+
+    algs = []
+
+    if Preferences.@load_preference("PrecompileMIRKNLLS", false)
+        for nlsolve in nlsolvers
+            append!(algs, [MIRK2(; jac_alg, nlsolve), MIRK6(; jac_alg, nlsolve)])
+        end
+    end
+
+    @compile_workload begin
+        @sync for prob in probs, alg in algs
+            Threads.@spawn solve(prob, alg; dt = 0.2, abstol = 1e-2)
+        end
+    end
+end
+
+export MIRK2, MIRK3, MIRK4, MIRK5, MIRK6
+export BVPJacobianAlgorithm
+
+end