diff --git a/README.md b/README.md index 6c015b4c..f61d47fe 100644 --- a/README.md +++ b/README.md @@ -37,6 +37,21 @@ prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan) sol = solve(prob, MIRK4(), dt = 0.05) ``` +```julia +tspan = (0.0, 1.0) +function lotka_volterra!(du, u, p, t) + du[1] = 1.5 * u[1] - 1.0 * u[1] * u[2] + du[2] = -3.0 * u[2] + 1.0 * u[1] * u[2] +end +function bc!(residual, u, p, t) + residual[1] = u[1] - 1.0 + residual[2] = u[2] - 1.0 +end +prob = BVProblem( + lotka_volterra!, bc!, [4.0, 2.0], tspan, lcons = [0.0, 0.0], ucons = [Inf, Inf]) +sol = solve(prob, MIRK4(; optimizer = Ipopt.Optimizer()), dt = 0.05) +``` + ## Available Solvers For the list of available solvers, please refer to the [DifferentialEquations.jl BVP Solvers page](https://docs.sciml.ai/DiffEqDocs/stable/solvers/bvp_solve/). For options for the `solve` command, see the [common solver options page](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). diff --git a/docs/make.jl b/docs/make.jl index f0ca989e..07383973 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -27,7 +27,7 @@ makedocs(; sitename = "BoundaryValueDiffEq.jl", clean = true, doctest = false, checkdocs = :exports, - warnonly = [:missing_docs], + warnonly = [:missing_docs, :cross_references], plugins = [bib, interlinks], format = Documenter.HTML(assets = ["assets/favicon.ico"], canonical = "https://docs.sciml.ai/BoundaryValueDiffEq/stable/"), diff --git a/docs/src/basics/solve.md b/docs/src/basics/solve.md index d99d5944..4d48df59 100644 --- a/docs/src/basics/solve.md +++ b/docs/src/basics/solve.md @@ -5,6 +5,7 @@ - `controller`: Error controller for collocation methods, default as `DefectControl()`, more controller options in [Error Control Adaptivity](@ref error_control). - `defect_threshold`: Monitor of the size of defect norm. Defaults to `0.1`. - `odesolve_kwargs`: OrdinaryDiffEq.jl solvers kwargs for passing to ODE solving in shooting methods. For more information, see the documentation for OrdinaryDiffEq: [Common Solver Options](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/). - - `nlsolve_kwargs`: NonlinearSolve.jl solvers kwargs for passing to nonlinear solving in collocation methods and shooting methods. For more information, see the documentation for NonlinearSolve: [Common Solver Options](https://docs.sciml.ai/NonlinearSolve/stable/basics/solve/). The default absolute tolerance of nonlinear solving in collocaio + - `nlsolve_kwargs`: NonlinearSolve.jl solvers kwargs for passing to nonlinear solving in collocation methods and shooting methods. For more information, see the documentation for NonlinearSolve: [Common Solver Options](https://docs.sciml.ai/NonlinearSolve/stable/basics/solve/). The default internal nonlinear solver is [customized polyalgorithm](https://github.com/SciML/BoundaryValueDiffEq.jl/blob/master/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl) and the default absolute tolerance of nonlinear solving in collocation and shooting methods is `1e-6`. + - `optimize_kwargs`: Optimization.jl solvers kwargs for passing to optimization problem solving in collocation methods and shooting methods. For more information, see the documentation for Optimization: [Common Solver Options](https://docs.sciml.ai/Optimization/stable/API/solve/). The internal optimization solver should be specified and the default absolute tolerance of optimization problem solving in collocation and shooting methods is `1e-6`. - `verbose`: Toggles whether warnings are thrown when the solver exits early. Defaults to `true`. - `ensemblealg`: Whether `MultipleShooting` uses multithreading, default as `EnsembleThreads()`. For more information, see the documentation for OrdinaryDiffEq: [EnsembleAlgorithms](https://docs.sciml.ai/DiffEqDocs/latest/features/ensemble/#EnsembleAlgorithms). diff --git a/docs/src/tutorials/unknown_parameters.md b/docs/src/tutorials/unknown_parameters.md index b064d5ba..5bb31530 100644 --- a/docs/src/tutorials/unknown_parameters.md +++ b/docs/src/tutorials/unknown_parameters.md @@ -54,7 +54,7 @@ sol = solve(bvp, MIRK4(), dt = 0.05) plot(sol) ``` -after solving the boundary value problem, the estimated unknown parameters can be accessed with +after solving the boundary value problem, the estimated unknown parameters can be accessed in the solution ```@example unknown sol.prob.p diff --git a/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl b/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl index 1c2725ea..38752727 100644 --- a/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl +++ b/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl @@ -6,11 +6,11 @@ using AlmostBlockDiagonals: AlmostBlockDiagonals, IntermediateAlmostBlockDiagona using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, AbstractBoundaryValueDiffEqCache, BVPJacobianAlgorithm, __extract_problem_details, concrete_jacobian_algorithm, - __Fix3, __concrete_nonlinearsolve_algorithm, + __Fix3, __concrete_solve_algorithm, __internal_nlsolve_problem, __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!, __extract_mesh, get_dense_ad, - __get_bcresid_prototype, __split_kwargs, - __default_nonsparse_ad + __get_bcresid_prototype, __split_kwargs, __concrete_kwargs, + __default_nonsparse_ad, __construct_internal_problem using ConcreteStructs: @concrete using DiffEqBase: DiffEqBase diff --git a/lib/BoundaryValueDiffEqAscher/src/algorithms.jl b/lib/BoundaryValueDiffEqAscher/src/algorithms.jl index 7e9d522d..784579ff 100644 --- a/lib/BoundaryValueDiffEqAscher/src/algorithms.jl +++ b/lib/BoundaryValueDiffEqAscher/src/algorithms.jl @@ -14,6 +14,9 @@ for stage in (1, 2, 3, 4, 5, 6, 7) - `nlsolve`: Internal Nonlinear solver. Any solver which conforms to the SciML `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. Note that any autodiff argument for + the solver will be ignored and a custom jacobian algorithm will be used. - `max_num_subintervals`: Number of maximal subintervals, default as 3000. - `zeta`: side condition points, should always be provided. @@ -46,8 +49,9 @@ for stage in (1, 2, 3, 4, 5, 6, 7) } ``` """ - @kwdef struct $(alg){N, J <: BVPJacobianAlgorithm} <: AbstractAscher + @kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm} <: AbstractAscher nlsolve::N = nothing + optimize::O = nothing zeta::Vector{Float64} = nothing jac_alg::J = BVPJacobianAlgorithm() max_num_subintervals::Int = 3000 diff --git a/lib/BoundaryValueDiffEqAscher/src/ascher.jl b/lib/BoundaryValueDiffEqAscher/src/ascher.jl index d442106a..41804f46 100644 --- a/lib/BoundaryValueDiffEqAscher/src/ascher.jl +++ b/lib/BoundaryValueDiffEqAscher/src/ascher.jl @@ -39,6 +39,7 @@ TU valstr nlsolve_kwargs + optimize_kwargs kwargs end @@ -60,7 +61,8 @@ end function SciMLBase.__init( prob::BVProblem, alg::AbstractAscher; dt = 0.0, controller = GlobalErrorControl(), - adaptive = true, abstol = 1e-4, nlsolve_kwargs = (; abstol = abstol), kwargs...) + adaptive = true, abstol = 1e-4, nlsolve_kwargs = (; abstol = abstol), + optimize_kwargs = (; abstol = abstol), kwargs...) (; tspan, p) = prob _, T, ncy, n, u0 = __extract_problem_details(prob; dt, check_positive_dt = true) t₀, t₁ = tspan @@ -147,9 +149,9 @@ function SciMLBase.__init( g = build_almost_block_diagonals(zeta, ncomp, mesh, T) cache = AscherCache{iip, T}( prob, f, jac, bc, bcjac, k, copy(mesh), mesh, mesh_dt, ncomp, ny, p, zeta, - fixpnt, alg, prob.problem_type, bcresid_prototype, residual, zval, yval, - gval, err, g, w, v, lz, ly, dmz, delz, deldmz, dqdmz, dmv, pvtg, pvtw, TU, - valst, nlsolve_kwargs, (; abstol, dt, adaptive, controller, kwargs...)) + fixpnt, alg, prob.problem_type, bcresid_prototype, residual, zval, yval, gval, + err, g, w, v, lz, ly, dmz, delz, deldmz, dqdmz, dmv, pvtg, pvtw, TU, valst, + nlsolve_kwargs, optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...)) return cache end @@ -176,8 +178,10 @@ function __perform_ascher_iteration(cache::AscherCache{iip, T}, abstol, adaptive iip, T} info::ReturnCode.T = ReturnCode.Success nlprob = __construct_nlproblem(cache) - nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve) - nlsol = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs...) + solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize) + kwargs = __concrete_kwargs( + cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs) + nlsol = __solve(nlprob, solve_alg; kwargs...) error_norm = 2 * abstol info = nlsol.retcode @@ -203,7 +207,7 @@ function __perform_ascher_iteration(cache::AscherCache{iip, T}, abstol, adaptive __expand_cache_for_error!(cache) _nlprob = __construct_nlproblem(cache) - nlsol = __solve(_nlprob, nlsolve_alg; cache.nlsolve_kwargs...) + nlsol = solve(_nlprob, solve_alg; kwargs...) error_norm = error_estimate!(cache) if norm(error_norm) > abstol @@ -346,9 +350,9 @@ function __construct_nlproblem(cache::AscherCache{iip, T}) where {iip, T} jac_prototype, u, diffmode, jac_cache, loss, cache.p) end - nlf = NonlinearFunction{iip}( - loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype) - return __internal_nlsolve_problem(cache.prob, similar(lz), lz, nlf, lz, cache.p) + return __construct_internal_problem( + cache.prob, alg, loss, jac, jac_prototype, resid_prototype, + lz, cache.p, cache.ncomp, length(cache.mesh)) end function __ascher_mpoint_jacobian!(J, x, diffmode, diffcache, loss, resid, p) diff --git a/lib/BoundaryValueDiffEqAscher/src/utils.jl b/lib/BoundaryValueDiffEqAscher/src/utils.jl index 0c5efb5e..3772efc1 100644 --- a/lib/BoundaryValueDiffEqAscher/src/utils.jl +++ b/lib/BoundaryValueDiffEqAscher/src/utils.jl @@ -112,7 +112,7 @@ end return nothing end -@views function recursive_flatten!(y::Vector, x::Vector{Vector{T}}) where {T} +@views function recursive_flatten!(y::AbstractArray, x::Vector{Vector{T}}) where {T} i = 0 for xᵢ in x copyto!(y[(i + 1):(i + length(xᵢ))], xᵢ) diff --git a/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl b/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl index b4335aa3..f3ff338a 100644 --- a/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl +++ b/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl @@ -29,7 +29,7 @@ include("utils.jl") include("algorithms.jl") include("abstract_types.jl") include("alg_utils.jl") -include("default_nlsolve.jl") +include("default_internal_solve.jl") include("calc_errors.jl") function SciMLBase.__solve(prob::AbstractBVProblem, diff --git a/lib/BoundaryValueDiffEqCore/src/algorithms.jl b/lib/BoundaryValueDiffEqCore/src/algorithms.jl index 0001712f..7d6cdf38 100644 --- a/lib/BoundaryValueDiffEqCore/src/algorithms.jl +++ b/lib/BoundaryValueDiffEqCore/src/algorithms.jl @@ -26,3 +26,12 @@ function Base.show(io::IO, alg::AbstractBoundaryValueDiffEqAlgorithm) print(io, join(modifiers, ", ")) print(io, ")") end + +# Check what's the internal solver, nonlinear or optimization? +function __internal_solver(alg::AbstractBoundaryValueDiffEqAlgorithm) + # We don't allow both `nlsolve` and `optimize` to be specified at the same time + (isnothing(alg.nlsolve) && isnothing(alg.optimize)) && + error("Either `nlsolve` or `optimize` must be specified in the algorithm, but not both.") + isnothing(alg.nlsolve) && return alg.optimize + isnothing(alg.optimize) && return alg.nlsolve +end diff --git a/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl b/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl similarity index 62% rename from lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl rename to lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl index f57a1324..800b6e95 100644 --- a/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl +++ b/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl @@ -46,11 +46,37 @@ function __FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = noth return NonlinearSolvePolyAlgorithm(algs) end -@inline __concrete_nonlinearsolve_algorithm(prob, alg) = alg -@inline function __concrete_nonlinearsolve_algorithm(prob, ::Nothing) +""" + __concrete_solve_algorithm(prob, nlsolve_alg, optimize_alg) + +Automatic solver choosing according to the input solver. +If none of the solvers are specified, we use nonlinear solvers from NonlinearSolve.jl. +If both of the nonlinear solver and optimization solver are specified, we throw an error. +If only one of the nonlinear solver and optimization solver is specified, we use that solver. +""" +@inline __concrete_solve_algorithm(prob, alg) = alg +@inline __concrete_solve_algorithm(prob, alg, ::Nothing) = alg +@inline __concrete_solve_algorithm(prob, ::Nothing, alg) = alg +@inline __concrete_solve_algorithm(prob, + alg1, + alg2) = error("Both `nlsolve` and `optimize` are specified in the algorithm, but only one of them is allowed. Please specify only one of them.") +@inline function __concrete_solve_algorithm(prob, ::Nothing) + if prob isa NonlinearLeastSquaresProblem + return __FastShortcutBVPCompatibleNLLSPolyalg(eltype(prob.u0)) + else + return __FastShortcutBVPCompatibleNonlinearPolyalg(eltype(prob.u0)) + end +end +@inline function __concrete_solve_algorithm(prob, ::Nothing, ::Nothing) if prob isa NonlinearLeastSquaresProblem return __FastShortcutBVPCompatibleNLLSPolyalg(eltype(prob.u0)) else return __FastShortcutBVPCompatibleNonlinearPolyalg(eltype(prob.u0)) end end + +# Some optimization algorithms (solvers from interfacing packages) don't support the __solve(prob) interface +@inline __internal_solve( + prob::Union{SciMLBase.NonlinearProblem, SciMLBase.NonlinearLeastSquaresProblem}, + alg; kwargs...) = __solve(prob, alg; kwargs...) +@inline __internal_solve(prob::SciMLBase.OptimizationProblem, alg; kwargs...) = solve(prob, alg; kwargs...) diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl index af518e7e..b89cfff3 100644 --- a/lib/BoundaryValueDiffEqCore/src/utils.jl +++ b/lib/BoundaryValueDiffEqCore/src/utils.jl @@ -467,6 +467,19 @@ end return NonlinearProblem(args...; kwargs...) end +# Construct the internal OptimizationProblem +@inline function __internal_optimization_problem( + ::BVProblem{uType, tType, iip}, args...; kwargs...) where {uType, tType, iip} + prob = OptimizationProblem(args...; kwargs...) + return prob +end + +@inline function __internal_optimization_problem(::SecondOrderBVProblem{uType, tType, iip}, + args...; kwargs...) where {uType, tType, iip} + prob = OptimizationProblem(args...; kwargs...) + return prob +end + # Handling Initial Guesses """ __extract_u0(u₀, t₀) @@ -560,10 +573,14 @@ end end # Construct BVP Solution -function __build_solution(prob::AbstractBVProblem, odesol, nlsol) +function __build_solution(prob::AbstractBVProblem, odesol, nlsol::SciMLBase.NonlinearSolution) retcode = ifelse(SciMLBase.successful_retcode(nlsol), odesol.retcode, nlsol.retcode) return SciMLBase.solution_new_original_retcode(odesol, nlsol, retcode, nlsol.resid) end +function __build_solution(prob::AbstractBVProblem, odesol, optsol::SciMLBase.OptimizationSolution) + retcode = ifelse(SciMLBase.successful_retcode(optsol), odesol.retcode, optsol.retcode) + return SciMLBase.solution_new_original_retcode(odesol, optsol, retcode, zeros(length(first(odesol)))) # Need a patch in SciMLBase +end # Fix3 @concrete struct __Fix3 @@ -597,3 +614,124 @@ end function __split_kwargs(; abstol, adaptive, controller, kwargs...) return ((abstol, adaptive, controller), (; abstol, adaptive, kwargs...)) end + +@inline __concrete_kwargs(nlsolve, ::Nothing, nlsolve_kwargs, optimize_kwargs) = (; + nlsolve_kwargs...) +@inline __concrete_kwargs(::Nothing, optimize, nlsolve_kwargs, optimize_kwargs) = (;) # Doesn't support for now +@inline __concrete_kwargs(::Nothing, ::Nothing, nlsolve_kwargs, optimize_kwargs) = (; + nlsolve_kwargs...) + +## Optimization solver related utils ## + +@inline __default_cost(::Nothing) = (x, p) -> 0.0 +@inline __default_cost(f) = f +@inline __default_cost(fun::BVPFunction) = __default_cost(fun.cost) + +@inline function __extract_lcons_ucons(prob::AbstractBVProblem, ::Type{T}, M, N) where {T} + lcons = if isnothing(prob.lcons) + zeros(T, N*M) + else + lcons_length = length(prob.lcons) + vcat(prob.lcons, zeros(T, N*M - lcons_length)) + end + ucons = if isnothing(prob.ucons) + zeros(T, N*M) + else + ucons_length = length(prob.ucons) + vcat(prob.ucons, zeros(T, N*M - ucons_length)) + end + return lcons, ucons +end + +""" + __construct_internal_problem + +Constructs the internal problem based on the type of the boundary value problem and the +algorithm used. It returns either a `NonlinearProblem` or an `OptimizationProblem`. +""" +function __construct_internal_problem(prob::AbstractBVProblem, alg, loss, jac, + jac_prototype, resid_prototype, y, p, M::Int, N::Int) + T = eltype(y) + # multiple shooting always use iip + iip = SciMLBase.isinplace(prob) + if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize)) + !isnothing(prob.lcons) || + !isnothing(prob.ucons) && + throw(ArgumentError("Inequality constraints can only be handled by optimization solvers, please specify e.g. `MIRK4(; optimize = Ipopt.Optimizer())`.")) + nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype, + jac_prototype = jac_prototype) + return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p) + else + optf = OptimizationFunction{true}(__default_cost(prob.f), AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem + cons = loss, + cons_j = jac, cons_jac_prototype = jac_prototype) + lcons, ucons = __extract_lcons_ucons(prob, T, M, N) + return __internal_optimization_problem( + prob, optf, y, p; lcons = lcons, ucons = ucons) + end +end + +function __construct_internal_problem(prob::TwoPointBVProblem, alg, loss, jac, + jac_prototype, resid_prototype, y, p, M::Int, N::Int) + T = eltype(y) + iip = SciMLBase.isinplace(prob) + if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize)) + !isnothing(prob.lcons) || + !isnothing(prob.ucons) && + throw(ArgumentError("Inequality constraints can only be handled by optimization solvers, please specify e.g. `MIRK4(; optimize = Ipopt.Optimizer())`.")) + nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype, + jac_prototype = jac_prototype) + return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p) + else + optf = OptimizationFunction{true}( + __default_cost(prob.f), get_dense_ad(alg.jac_alg.diffmode), + cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype) + lcons, ucons = __extract_lcons_ucons(prob, T, M, N) + + return __internal_optimization_problem( + prob, optf, y, p; lcons = lcons, ucons = ucons) + end +end +# Multiple shooting only use inplace version internal problem constructor +function __construct_internal_problem(prob, alg, loss, jac, jac_prototype, + resid_prototype, y, p, M::Int, N::Int, ::Nothing) + T = eltype(y) + if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize)) + !isnothing(prob.lcons) || + !isnothing(prob.ucons) && + throw(ArgumentError("Inequality constraints can only be handled by optimization solvers, please specify e.g. `MIRK4(; optimize = Ipopt.Optimizer())`.")) + nlf = NonlinearFunction{true}(loss; jac = jac, resid_prototype = resid_prototype, + jac_prototype = jac_prototype) + return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p) + else + optf = OptimizationFunction{true}( + __default_cost(prob.f), get_dense_ad(alg.jac_alg.diffmode), + cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype) + lcons, ucons = __extract_lcons_ucons(prob, T, M, N) + + return __internal_optimization_problem( + prob, optf, y, p; lcons = lcons, ucons = ucons) + end +end +function __construct_internal_problem( + prob::TwoPointBVProblem, alg, loss, jac, jac_prototype, + resid_prototype, y, p, M::Int, N::Int, ::Nothing) + T = eltype(y) + iip = SciMLBase.isinplace(prob) + if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize)) + !isnothing(prob.lcons) || + !isnothing(prob.ucons) && + throw(ArgumentError("Inequality constraints can only be handled by optimization solvers, please specify e.g. `MIRK4(; optimize = Ipopt.Optimizer())`.")) + nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype, + jac_prototype = jac_prototype) + return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p) + else + optf = OptimizationFunction{true}( + __default_cost(prob.f), get_dense_ad(alg.jac_alg.nonbc_diffmode), + cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype) + lcons, ucons = __extract_lcons_ucons(prob, T, M, N) + + return __internal_optimization_problem( + prob, optf, y, p; lcons = lcons, ucons = ucons) + end +end diff --git a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl index 824df0bb..7d6b6f91 100644 --- a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl +++ b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl @@ -6,13 +6,13 @@ using BandedMatrices: BandedMatrix, Ones using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, AbstractBoundaryValueDiffEqCache, BVPJacobianAlgorithm, recursive_flatten, recursive_flatten!, recursive_unflatten!, - __concrete_nonlinearsolve_algorithm, diff!, EvalSol, + __concrete_solve_algorithm, diff!, EvalSol, concrete_jacobian_algorithm, eval_bc_residual, interval, eval_bc_residual!, get_tmp, __maybe_matmul!, __resize!, __extract_problem_details, __initial_guess, nodual_value, __maybe_allocate_diffcache, __restructure_sol, __get_bcresid_prototype, __vec, __vec_f, __vec_f!, __vec_bc, - __vec_bc!, recursive_flatten_twopoint!, + __vec_bc!, recursive_flatten_twopoint!, __concrete_kwargs, __internal_nlsolve_problem, __extract_mesh, __extract_u0, __default_coloring_algorithm, __maybe_allocate_diffcache, __restructure_sol, __get_bcresid_prototype, safe_similar, @@ -20,9 +20,10 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, recursive_flatten_twopoint!, __internal_nlsolve_problem, __extract_mesh, __extract_u0, DiffCacheNeeded, NoDiffCacheNeeded, __has_initial_guess, - __initial_guess_length, __initial_guess_on_mesh, - __flatten_initial_guess, __build_solution, __Fix3, - __split_kwargs, _sparse_like, get_dense_ad + __construct_internal_problem, __initial_guess_length, + __initial_guess_on_mesh, __flatten_initial_guess, + __build_solution, __Fix3, __split_kwargs, _sparse_like, + get_dense_ad, __internal_optimization_problem using ConcreteStructs: @concrete using DiffEqBase: DiffEqBase diff --git a/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl b/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl index 11870451..d8d80e09 100644 --- a/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl +++ b/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl @@ -106,7 +106,7 @@ end end τ = (t - mesh[j]) - nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve) + nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve) nestprob_p = zeros(T, cache.M + 2) yᵢ = copy(cache.y[j].du) @@ -152,7 +152,7 @@ end end τ = (t - mesh[j]) - nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve) + nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve) nestprob_p = zeros(T, cache.M + 2) yᵢ = copy(cache.y[j]) @@ -457,7 +457,7 @@ end (; f, mesh, mesh_dt, defect, ITU, nest_prob, alg) = cache (; q_coeff, τ_star) = ITU - nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, cache.alg.nlsolve) + nlsolve_alg = __concrete_solve_algorithm(nest_prob, cache.alg.nlsolve) nestprob_p = zeros(T, cache.M + 2) for i in 1:(length(mesh) - 1) @@ -509,7 +509,7 @@ end (; f, mesh, mesh_dt, defect, ITU, nest_prob, alg) = cache (; q_coeff, τ_star) = ITU - nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, cache.alg.nlsolve) + nlsolve_alg = __concrete_solve_algorithm(nest_prob, cache.alg.nlsolve) nestprob_p = zeros(T, cache.M + 2) for i in 1:(length(mesh) - 1) diff --git a/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl b/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl index c8156a86..1220506c 100644 --- a/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl +++ b/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl @@ -17,6 +17,9 @@ for stage in (1, 2, 3, 5, 7) `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. Note that any autodiff argument for + the solver will be ignored and a custom jacobian algorithm will be used. - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to use based on the input types and problem type. @@ -87,23 +90,15 @@ for stage in (1, 2, 3, 5, 7) } ``` """ - Base.@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK + Base.@kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK nlsolve::N = nothing + optimize::O = nothing jac_alg::J = BVPJacobianAlgorithm() nested_nlsolve::Bool = false nested_nlsolve_kwargs::NamedTuple = (;) defect_threshold::T = 0.1 max_num_subintervals::Int = 3000 end - $(alg)(nlsolve::N, - jac_alg::J; - nested = false, - nested_nlsolve_kwargs::NamedTuple = (;), - defect_threshold::T = 0.1, - max_num_subintervals::Int = 3000) where {N, - J, - T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs, - defect_threshold, max_num_subintervals) end end @@ -123,6 +118,10 @@ for stage in (2, 3, 4, 5) `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. Note that any autodiff argument for + the solver will be ignored and a custom jacobian algorithm will be used. + - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to use based on the input types and problem type. @@ -192,23 +191,15 @@ for stage in (2, 3, 4, 5) } ``` """ - Base.@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK + Base.@kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK nlsolve::N = nothing + optimize::O = nothing jac_alg::J = BVPJacobianAlgorithm() nested_nlsolve::Bool = false nested_nlsolve_kwargs::NamedTuple = (;) defect_threshold::T = 0.1 max_num_subintervals::Int = 3000 end - $(alg)(nlsolve::N, - jac_alg::J; - nested = false, - nested_nlsolve_kwargs::NamedTuple = (;), - defect_threshold::T = 0.1, - max_num_subintervals::Int = 3000) where {N, - J, - T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs, - defect_threshold, max_num_subintervals) end end @@ -228,6 +219,9 @@ for stage in (2, 3, 4, 5) `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. Note that any autodiff argument for + the solver will be ignored and a custom jacobian algorithm will be used. - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to use based on the input types and problem type. @@ -299,23 +293,15 @@ for stage in (2, 3, 4, 5) } ``` """ - Base.@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK + Base.@kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK nlsolve::N = nothing + optimize::O = nothing jac_alg::J = BVPJacobianAlgorithm() nested_nlsolve::Bool = false nested_nlsolve_kwargs::NamedTuple = (;) defect_threshold::T = 0.1 max_num_subintervals::Int = 3000 end - $(alg)(nlsolve::N, - jac_alg::J; - nested = false, - nested_nlsolve_kwargs::NamedTuple = (;), - defect_threshold::T = 0.1, - max_num_subintervals::Int = 3000) where {N, - J, - T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs, - defect_threshold, max_num_subintervals) end end @@ -335,6 +321,9 @@ for stage in (2, 3, 4, 5) `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. Note that any autodiff argument for + the solver will be ignored and a custom jacobian algorithm will be used. - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to use based on the input types and problem type. @@ -406,23 +395,15 @@ for stage in (2, 3, 4, 5) } ``` """ - Base.@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK + Base.@kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK nlsolve::N = nothing + optimize::O = nothing jac_alg::J = BVPJacobianAlgorithm() nested_nlsolve::Bool = false nested_nlsolve_kwargs::NamedTuple = (;) defect_threshold::T = 0.1 max_num_subintervals::Int = 3000 end - $(alg)(nlsolve::N, - jac_alg::J; - nested = false, - nested_nlsolve_kwargs::NamedTuple = (;), - defect_threshold::T = 0.1, - max_num_subintervals::Int = 3000) where {N, - J, - T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs, - defect_threshold, max_num_subintervals) end end diff --git a/lib/BoundaryValueDiffEqFIRK/src/collocation.jl b/lib/BoundaryValueDiffEqFIRK/src/collocation.jl index 0c60b7df..bbdcf135 100644 --- a/lib/BoundaryValueDiffEqFIRK/src/collocation.jl +++ b/lib/BoundaryValueDiffEqFIRK/src/collocation.jl @@ -121,7 +121,7 @@ end T = eltype(u) nestprob_p = vcat(T(mesh[1]), T(mesh_dt[1]), get_tmp(y[1], u)) - nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve) + nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve) for i in eachindex(k_discrete) residᵢ = residual[i] @@ -151,7 +151,7 @@ end T = eltype(u) nestprob_p = vcat(T(mesh[1]), T(mesh_dt[1]), y[1]) - nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, cache.alg.nlsolve) + nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, cache.alg.nlsolve) for i in eachindex(k_discrete) residᵢ = residual[i] @@ -265,7 +265,7 @@ end T = eltype(u) nestprob_p = vcat(T(mesh[1]), T(mesh_dt[1]), get_tmp(y[1], u)) - nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve) + nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve) for i in eachindex(k_discrete) residᵢ = residuals[i] @@ -296,7 +296,7 @@ end T = eltype(u) nestprob_p = vcat(T(mesh[1]), T(mesh_dt[1]), y[1]) - nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve) + nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve) for i in eachindex(k_discrete) residᵢ = residuals[i] diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl index 3a58727b..a6d7abb5 100644 --- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl +++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl @@ -27,6 +27,7 @@ nest_prob resid_size nlsolve_kwargs + optimize_kwargs kwargs end @@ -60,6 +61,7 @@ Base.eltype(::FIRKCacheNested{iip, T}) where {iip, T} = T defect resid_size nlsolve_kwargs + optimize_kwargs kwargs end @@ -85,22 +87,27 @@ end function SciMLBase.__init( prob::BVProblem, alg::AbstractFIRK; dt = 0.0, abstol = 1e-6, adaptive = true, - controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), kwargs...) + controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), + optimize_kwargs = (; abstol = abstol), kwargs...) if alg.nested_nlsolve return init_nested(prob, alg; dt = dt, abstol = abstol, adaptive = adaptive, - controller = controller, nlsolve_kwargs = nlsolve_kwargs, kwargs...) + controller = controller, nlsolve_kwargs = nlsolve_kwargs, + optimize_kwargs = optimize_kwargs, kwargs...) else return init_expanded(prob, alg; dt = dt, abstol = abstol, adaptive = adaptive, - controller = controller, nlsolve_kwargs = nlsolve_kwargs, kwargs...) + controller = controller, nlsolve_kwargs = nlsolve_kwargs, + optimize_kwargs = optimize_kwargs, kwargs...) end end function init_nested( prob::BVProblem, alg::AbstractFIRK; dt = 0.0, abstol = 1e-6, adaptive = true, - controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), kwargs...) + controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), + optimize_kwargs = (; abstol = abstol), kwargs...) @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg) iip = isinplace(prob) + @assert (iip || isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers " if adaptive && isa(alg, FIRKNoAdaptivity) error("Algorithm doesn't support adaptivity. Please choose a higher order algorithm.") end @@ -197,25 +204,25 @@ function init_nested( end return FIRKCacheNested{iip, T, typeof(diffcache), fit_parameters}( - alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type, - prob.p, alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, - y, y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, nestprob, resid₁_size, - nlsolve_kwargs, (; abstol, dt, adaptive, controller, kwargs...)) + alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type, prob.p, + alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, y, y₀, residual, + fᵢ_cache, fᵢ₂_cache, defect, nestprob, resid₁_size, nlsolve_kwargs, + optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...)) end function init_expanded( prob::BVProblem, alg::AbstractFIRK; dt = 0.0, abstol = 1e-6, adaptive = true, - controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), kwargs...) + controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), + optimize_kwargs = (; abstol = abstol), kwargs...) @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg) - + iip = isinplace(prob) + @assert (iip || isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers " if adaptive && isa(alg, FIRKNoAdaptivity) error("Algorithm $(alg) doesn't support adaptivity. Please choose a higher order algorithm.") end diffcache = __cache_trait(alg.jac_alg) fit_parameters = haskey(prob.kwargs, :fit_parameters) - iip = isinplace(prob) - t₀, t₁ = prob.tspan ig, T, M, @@ -297,9 +304,9 @@ function init_expanded( return FIRKCacheExpand{iip, T, typeof(diffcache), fit_parameters}( alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type, - prob.p, alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, - y, y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, resid₁_size, - nlsolve_kwargs, (; abstol, dt, adaptive, controller, kwargs...)) + prob.p, alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, y, + y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, resid₁_size, nlsolve_kwargs, + optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...)) end """ @@ -398,9 +405,11 @@ function SciMLBase.solve!(cache::FIRKCacheNested{ end function __perform_firk_iteration(cache::Union{FIRKCacheExpand, FIRKCacheNested}, abstol, adaptive::Bool) - nlprob = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀)) - nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve) - sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true) + nlprob = __construct_problem(cache, vec(cache.y₀), copy(cache.y₀)) + solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize) + kwargs = __concrete_kwargs( + cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs) + sol_nlprob = __solve(nlprob, solve_alg; kwargs...) recursive_unflatten!(cache.y₀, sol_nlprob.u) defect_norm = 2 * abstol @@ -445,7 +454,7 @@ function __perform_firk_iteration(cache::Union{FIRKCacheExpand, FIRKCacheNested} end # Constructing the Nonlinear Problem -function __construct_nlproblem(cache::Union{FIRKCacheNested{iip}, FIRKCacheExpand{iip}}, +function __construct_problem(cache::Union{FIRKCacheNested{iip}, FIRKCacheExpand{iip}}, y::AbstractVector, y₀::AbstractVectorOfArray) where {iip} pt = cache.problem_type (; jac_alg) = cache.alg @@ -485,10 +494,10 @@ function __construct_nlproblem(cache::Union{FIRKCacheNested{iip}, FIRKCacheExpan u, p, cache.y, pt, cache.bc, cache.mesh, cache, eval_sol, trait) end - return __construct_nlproblem(cache, y, loss_bc, loss_collocation, loss, pt) + return __construct_problem(cache, y, loss_bc, loss_collocation, loss, pt) end -function __construct_nlproblem( +function __construct_problem( cache::FIRKCacheExpand{iip}, y, loss_bc::BC, loss_collocation::C, loss::LF, ::StandardBVProblem) where {iip, BC, C, LF} (; alg, stage) = cache @@ -569,13 +578,11 @@ function __construct_nlproblem( end resid_prototype = vcat(resid_bc, resid_collocation) - nlf = NonlinearFunction{iip}( - loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype) - - return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p) + return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype, + resid_prototype, y, cache.p, cache.M, (N - 1) * (stage + 1) + 1) end -function __construct_nlproblem( +function __construct_problem( cache::FIRKCacheExpand{iip}, y, loss_bc::BC, loss_collocation::C, loss::LF, ::TwoPointBVProblem) where {iip, BC, C, LF} (; jac_alg) = cache.alg @@ -630,12 +637,11 @@ function __construct_nlproblem( end resid_prototype = copy(resid) - nlf = NonlinearFunction{iip}( - loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype) - return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p) + return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype, + resid_prototype, y, cache.p, cache.M, (N - 1) * (stage + 1) + 1) end -function __construct_nlproblem( +function __construct_problem( cache::FIRKCacheNested{iip}, y, loss_bc::BC, loss_collocation::C, loss::LF, ::StandardBVProblem) where {iip, BC, C, LF} (; jac_alg) = cache.alg @@ -710,13 +716,11 @@ function __construct_nlproblem( end resid_prototype = vcat(resid_bc, resid_collocation) - nlf = NonlinearFunction{iip}( - loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype) - - return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p) + return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype, + resid_prototype, y, cache.p, cache.M, N) end -function __construct_nlproblem( +function __construct_problem( cache::FIRKCacheNested{iip}, y, loss_bc::BC, loss_collocation::C, loss::LF, ::TwoPointBVProblem) where {iip, BC, C, LF} (; jac_alg) = cache.alg @@ -761,9 +765,8 @@ function __construct_nlproblem( end resid_prototype = copy(resid) - nlf = NonlinearFunction{iip}( - loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype) - return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p) + return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype, + resid_prototype, y, cache.p, cache.M, N) end @views function __firk_loss!(resid, u, p, y, pt::StandardBVProblem, bc!::BC, residual, @@ -825,7 +828,6 @@ end @views function __firk_loss(u, p, y::AbstractVector, pt::TwoPointBVProblem, bc::Tuple{BC1, BC2}, mesh, cache, _, trait) where {BC1, BC2} y_ = recursive_unflatten!(y, u) - soly_ = VectorOfArray(y_) resid_bca, resid_bcb = eval_bc_residual(pt, bc, y_, p, mesh) resid_co = Φ(cache, y_, u, trait) return vcat(resid_bca, mapreduce(vec, vcat, resid_co), resid_bcb) diff --git a/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl b/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl index 2adc0ab4..56b5e9b0 100644 --- a/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl +++ b/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl @@ -95,7 +95,7 @@ end τ = (t - mesh[j]) length_z = length(z) - nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve) + nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve) nestprob_p = zeros(T, cache.M + 2) yᵢ = copy(cache.y[j].du) @@ -141,7 +141,7 @@ end τ = (t - mesh[j]) length_dz = length(dz) - nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve) + nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve) nestprob_p = zeros(T, cache.M + 2) yᵢ = copy(cache.y[j].du) @@ -413,7 +413,7 @@ function (s::EvalSol{C})(tval::Number) where {C <: FIRKCacheNested} h = mesh_dt[j] τ = tval - t[j] - nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve) + nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve) nestprob_p = zeros(cache.M + 2) yᵢ = u[j] diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl index 88e496ed..68b72b69 100644 --- a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl +++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl @@ -27,7 +27,8 @@ AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))] for jac_alg in jac_algs - sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1) + sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested); + trajectories = 10, dt = 0.1) @test sol.converged end end @@ -39,7 +40,8 @@ AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))] for jac_alg in jac_algs - sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1) + sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested); + trajectories = 10, dt = 0.1) @test sol.converged end end @@ -50,7 +52,8 @@ AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))] for jac_alg in jac_algs - sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1) + sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested); + trajectories = 10, dt = 0.1) @test sol.converged end end @@ -61,7 +64,8 @@ AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))] for jac_alg in jac_algs - sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1) + sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested); + trajectories = 10, dt = 0.1) @test sol.converged end end diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl index 04cf30dd..e6a421cd 100644 --- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl +++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl @@ -232,33 +232,36 @@ end jac_alg = BVPJacobianAlgorithm(; bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff())) - nl_solve = NewtonRaphson() + nlsolve = NewtonRaphson() nested = false # Using ForwardDiff might lead to Cache expansion warnings - @test_nowarn solve(bvp1, LobattoIIIa2(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIa3(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIa4(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIa5(nl_solve, jac_alg; nested); dt = 0.005) - - @test_nowarn solve(bvp1, LobattoIIIb2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false) - @test_nowarn solve(bvp1, LobattoIIIb3(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIb4(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIb5(nl_solve, jac_alg; nested); dt = 0.005) - - @test_nowarn solve(bvp1, LobattoIIIc2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false) - @test_nowarn solve(bvp1, LobattoIIIc3(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIc4(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIc5(nl_solve, jac_alg; nested); dt = 0.005) - - @test_nowarn solve(bvp1, RadauIIa1(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false) - @test_nowarn solve(bvp1, RadauIIa2(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, RadauIIa3(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, RadauIIa5(nl_solve, jac_alg; nested); dt = 0.05) - @test_nowarn solve(bvp1, RadauIIa7(nl_solve, jac_alg; nested); dt = 0.05) + @test_nowarn solve(bvp1, LobattoIIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIa4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + + @test_nowarn solve( + bvp1, LobattoIIIb2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false) + @test_nowarn solve(bvp1, LobattoIIIb3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIb4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIb5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + + @test_nowarn solve( + bvp1, LobattoIIIc2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false) + @test_nowarn solve(bvp1, LobattoIIIc3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIc4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIc5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + + @test_nowarn solve( + bvp1, RadauIIa1(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false) + @test_nowarn solve(bvp1, RadauIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, RadauIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, RadauIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05) + @test_nowarn solve(bvp1, RadauIIa7(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05) end -@testitem "Interpolation" begin +@testitem "Interpolation" setup=[FIRKExpandedConvergenceTests] begin using LinearAlgebra λ = 1 @@ -291,26 +294,6 @@ end testTol = 1e-6 nested = false - for stage in (2, 3, 5, 7) - s = Symbol("RadauIIa$(stage)") - @eval radau_solver(::Val{$stage}) = $(s)(NewtonRaphson(), BVPJacobianAlgorithm(); nested) - end - - for stage in (3, 4, 5) - s = Symbol("LobattoIIIa$(stage)") - @eval lobattoIIIa_solver(::Val{$stage}) = $(s)(NewtonRaphson(), BVPJacobianAlgorithm(); nested) - end - - for stage in (3, 4, 5) - s = Symbol("LobattoIIIb$(stage)") - @eval lobattoIIIb_solver(::Val{$stage}) = $(s)(NewtonRaphson(), BVPJacobianAlgorithm(); nested) - end - - for stage in (3, 4, 5) - s = Symbol("LobattoIIIc$(stage)") - @eval lobattoIIIc_solver(::Val{$stage}) = $(s)(NewtonRaphson(), BVPJacobianAlgorithm(); nested) - end - @testset "Radau interpolations" begin @testset "Interpolation tests for RadauIIa$stage" for stage in (2, 3, 5, 7) @time sol = solve(prob_bvp_linear, radau_solver(Val(stage)); dt = 0.001) @@ -336,9 +319,11 @@ end for (id, lobatto_solver) in zip(("a", "b", "c"), ( lobattoIIIa_solver, lobattoIIIb_solver, lobattoIIIc_solver)) begin - @testset "Interpolation tests for LobattoIII$(id)$stage" for stage in - (3, 4, 5) - @time sol = solve(prob_bvp_linear, lobatto_solver(Val(stage)); dt = 0.001) + @testset "Interpolation tests for LobattoIII$(id)$stage" for stage in ( + 2, 3, 4, 5) + adaptive = ifelse(stage == 2, false, true) # LobattoIIIa2 is not adaptive + @time sol = solve( + prob_bvp_linear, lobatto_solver(Val(stage)); dt = 0.001, adaptive = adaptive) @test sol(0.001)≈[0.998687464, -1.312035941] atol=testTol @test sol(0.001; idxs = [1, 2])≈[0.998687464, -1.312035941] atol=testTol @test sol(0.001; idxs = 1)≈0.998687464 atol=testTol @@ -347,8 +332,10 @@ end @testset "Derivative Interpolation tests for lobatto$(id)$stage" for stage in ( - 3, 4, 5) - @time sol = solve(prob_bvp_linear, lobatto_solver(Val(stage)); dt = 0.001) + 2, 3, 4, 5) + adaptive = ifelse(stage == 2, false, true) # LobattoIIIa2 is not adaptive + @time sol = solve( + prob_bvp_linear, lobatto_solver(Val(stage)); dt = 0.001, adaptive = adaptive) sol_analytic = prob_bvp_linear_analytic(nothing, λ, 0.04) dsol_analytic = prob_bvp_linear_analytic_derivative(nothing, λ, 0.04) diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl index 5a4ba948..6f565b2d 100644 --- a/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl +++ b/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl @@ -27,7 +27,8 @@ AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))] for jac_alg in jac_algs - sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1) + sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested); + trajectories = 10, dt = 0.1) @test sol.converged end end @@ -39,7 +40,8 @@ AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))] for jac_alg in jac_algs - sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1) + sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested); + trajectories = 10, dt = 0.1) @test sol.converged end end @@ -50,7 +52,8 @@ AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))] for jac_alg in jac_algs - sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1) + sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested); + trajectories = 10, dt = 0.1) @test sol.converged end end @@ -61,7 +64,8 @@ AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))] for jac_alg in jac_algs - sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1) + sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested); + trajectories = 10, dt = 0.1) @test sol.converged end end diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl index 31ba3f1b..118fcf47 100644 --- a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl +++ b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl @@ -252,33 +252,36 @@ end jac_alg = BVPJacobianAlgorithm(; bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff())) - nl_solve = NewtonRaphson() + nlsolve = NewtonRaphson() nested = true # Using ForwardDiff might lead to Cache expansion warnings - @test_nowarn solve(bvp1, LobattoIIIa2(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIa3(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIa4(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIa5(nl_solve, jac_alg; nested); dt = 0.005) - - @test_nowarn solve(bvp1, LobattoIIIb2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false) - @test_nowarn solve(bvp1, LobattoIIIb3(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIb4(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIb5(nl_solve, jac_alg; nested); dt = 0.005) - - @test_nowarn solve(bvp1, LobattoIIIc2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false) - @test_nowarn solve(bvp1, LobattoIIIc3(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIc4(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, LobattoIIIc5(nl_solve, jac_alg; nested); dt = 0.005) - - @test_nowarn solve(bvp1, RadauIIa1(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false) - @test_nowarn solve(bvp1, RadauIIa2(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, RadauIIa3(nl_solve, jac_alg; nested); dt = 0.005) - @test_nowarn solve(bvp1, RadauIIa5(nl_solve, jac_alg; nested); dt = 0.05) - @test_nowarn solve(bvp1, RadauIIa7(nl_solve, jac_alg; nested); dt = 0.05) + @test_nowarn solve(bvp1, LobattoIIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIa4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + + @test_nowarn solve( + bvp1, LobattoIIIb2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false) + @test_nowarn solve(bvp1, LobattoIIIb3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIb4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIb5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + + @test_nowarn solve( + bvp1, LobattoIIIc2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false) + @test_nowarn solve(bvp1, LobattoIIIc3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIc4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, LobattoIIIc5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + + @test_nowarn solve( + bvp1, RadauIIa1(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false) + @test_nowarn solve(bvp1, RadauIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, RadauIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005) + @test_nowarn solve(bvp1, RadauIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05) + @test_nowarn solve(bvp1, RadauIIa7(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05) end -@testitem "Interpolation" begin +@testitem "Interpolation" setup=[FIRKNestedConvergenceTests] begin using LinearAlgebra λ = 1 @@ -311,26 +314,6 @@ end testTol = 1e-6 nested = true - for stage in (2, 3, 5, 7) - s = Symbol("RadauIIa$(stage)") - @eval radau_solver(::Val{$stage}) = $(s)(NewtonRaphson(), BVPJacobianAlgorithm(); nested) - end - - for stage in (3, 4, 5) - s = Symbol("LobattoIIIa$(stage)") - @eval lobattoIIIa_solver(::Val{$stage}) = $(s)(NewtonRaphson(), BVPJacobianAlgorithm(); nested) - end - - for stage in (3, 4, 5) - s = Symbol("LobattoIIIb$(stage)") - @eval lobattoIIIb_solver(::Val{$stage}) = $(s)(NewtonRaphson(), BVPJacobianAlgorithm(); nested) - end - - for stage in (3, 4, 5) - s = Symbol("LobattoIIIc$(stage)") - @eval lobattoIIIc_solver(::Val{$stage}) = $(s)(NewtonRaphson(), BVPJacobianAlgorithm(); nested) - end - @testset "Radau interpolations" begin @testset "Interpolation tests for RadauIIa$stage" for stage in (2, 3, 5, 7) @time sol = solve(prob_bvp_linear, radau_solver(Val(stage)); dt = 0.001) diff --git a/lib/BoundaryValueDiffEqMIRK/Project.toml b/lib/BoundaryValueDiffEqMIRK/Project.toml index a51b63e4..9477595b 100644 --- a/lib/BoundaryValueDiffEqMIRK/Project.toml +++ b/lib/BoundaryValueDiffEqMIRK/Project.toml @@ -23,8 +23,8 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -[sources.BoundaryValueDiffEqCore] -path = "../BoundaryValueDiffEqCore" +[sources] +BoundaryValueDiffEqCore = {path = "../BoundaryValueDiffEqCore"} [compat] ADTypes = "1.14" @@ -42,10 +42,12 @@ FastClosures = "0.3.2" ForwardDiff = "0.10.38, 1" Hwloc = "3" InteractiveUtils = "<0.0.1, 1" +Ipopt = "1.10.6" JET = "0.9.18" LinearAlgebra = "1.10" LinearSolve = "2.36.2, 3" Mooncake = "0.4.146" +OptimizationMOI = "0.5.5" OrdinaryDiffEqRosenbrock = "1" PreallocationTools = "0.4.24" PrecompileTools = "1.2" @@ -67,9 +69,11 @@ DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" Hwloc = "0e44f5e4-bd66-52a0-8798-143a42290a1d" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" +OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823" @@ -78,4 +82,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "DiffEqDevTools", "Enzyme", "Hwloc", "InteractiveUtils", "JET", "LinearSolve", "Mooncake", "OrdinaryDiffEqRosenbrock", "Random", "ReTestItems", "RecursiveArrayTools", "StaticArrays", "Test"] +test = ["Aqua", "DiffEqDevTools", "Enzyme", "Hwloc", "InteractiveUtils", "Ipopt", "JET", "LinearSolve", "Mooncake", "OptimizationMOI", "OrdinaryDiffEqRosenbrock", "Random", "ReTestItems", "RecursiveArrayTools", "StaticArrays", "Test"] diff --git a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl index 9d6bb6d7..9801f9b3 100644 --- a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl +++ b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl @@ -6,7 +6,7 @@ using BandedMatrices: BandedMatrix, Ones using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, AbstractBoundaryValueDiffEqCache, BVPJacobianAlgorithm, recursive_flatten, recursive_flatten!, recursive_unflatten!, - __concrete_nonlinearsolve_algorithm, diff!, EvalSol, + __concrete_solve_algorithm, diff!, EvalSol, concrete_jacobian_algorithm, eval_bc_residual, eval_bc_residual!, get_tmp, __maybe_matmul!, __resize!, __extract_problem_details, __initial_guess, interval, @@ -14,15 +14,16 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, __restructure_sol, __cache_trait, __get_bcresid_prototype, safe_similar, __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!, recursive_flatten_twopoint!, __internal_nlsolve_problem, - __extract_mesh, __extract_u0, __has_initial_guess, - __initial_guess_length, __initial_guess_on_mesh, - __flatten_initial_guess, __build_solution, __Fix3, - get_dense_ad, _sparse_like, AbstractErrorControl, - DefectControl, GlobalErrorControl, SequentialErrorControl, - HybridErrorControl, HOErrorControl, __use_both_error_control, - __default_coloring_algorithm, DiffCacheNeeded, - NoDiffCacheNeeded, __split_kwargs, - __FastShortcutNonlinearPolyalg + __internal_optimization_problem, __extract_mesh, + __extract_u0, __has_initial_guess, __initial_guess_length, + __initial_guess_on_mesh, __flatten_initial_guess, + __build_solution, __Fix3, get_dense_ad, _sparse_like, + AbstractErrorControl, DefectControl, GlobalErrorControl, + SequentialErrorControl, HybridErrorControl, HOErrorControl, + __use_both_error_control, __default_coloring_algorithm, + DiffCacheNeeded, NoDiffCacheNeeded, __split_kwargs, + __concrete_kwargs, __FastShortcutNonlinearPolyalg, + __construct_internal_problem, __internal_solve using ConcreteStructs: @concrete using DiffEqBase: DiffEqBase diff --git a/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl b/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl index 5bcfdcf0..750bdd26 100644 --- a/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl +++ b/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl @@ -496,7 +496,7 @@ end new_prob = remake(prob, u0 = high_sol) high_cache = SciMLBase.__init(new_prob, alg, adaptive = false) - high_nlprob = __construct_nlproblem(high_cache, vec(high_sol), VectorOfArray(high_sol.u)) + high_nlprob = __construct_problem(high_cache, vec(high_sol), VectorOfArray(high_sol.u)) high_sol_original = __solve(high_nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true) recursive_unflatten!(high_sol, high_sol_original.u) error_norm = global_error(VectorOfArray(copy(high_sol.u[1:2:end])), copy(cache.y₀), errors) @@ -513,7 +513,7 @@ end new_prob = remake(prob, u0 = high_sol) high_cache = SciMLBase.__init(new_prob, __high_order_method(alg), adaptive = false) - high_nlprob = __construct_nlproblem(high_cache, sol.u, high_sol) + high_nlprob = __construct_problem(high_cache, sol.u, high_sol) high_sol_nlprob = __solve(high_nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true) recursive_unflatten!(high_sol, high_sol_nlprob) error_norm = global_error(VectorOfArray(high_sol.u), cache.y₀, errors) diff --git a/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl b/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl index 996d161c..eae0b547 100644 --- a/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl +++ b/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl @@ -17,6 +17,10 @@ for order in (2, 3, 4, 5, 6) `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. Note that any autodiff argument for + the solver will be ignored and a custom jacobian algorithm will be used. Optimization + solvers should first be loaded to allow this functionality. - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to use based on the input types and problem type. @@ -48,8 +52,9 @@ for order in (2, 3, 4, 5, 6) } ``` """ - @kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractMIRK + @kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractMIRK nlsolve::N = nothing + optimize::O = nothing jac_alg::J = BVPJacobianAlgorithm() defect_threshold::T = 0.1 max_num_subintervals::Int = 3000 @@ -73,6 +78,9 @@ for order in (6) `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. Note that any autodiff argument for + the solver will be ignored and a custom jacobian algorithm will be used. - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to use based on the input types and problem type. @@ -104,8 +112,9 @@ for order in (6) } ``` """ - @kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractMIRK + @kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractMIRK nlsolve::N = nothing + optimize::O = nothing jac_alg::J = BVPJacobianAlgorithm() defect_threshold::T = 0.1 max_num_subintervals::Int = 3000 diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl index 266b9b8a..98d89646 100644 --- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl +++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl @@ -28,6 +28,7 @@ new_stages resid_size nlsolve_kwargs + optimize_kwargs kwargs end @@ -35,10 +36,12 @@ Base.eltype(::MIRKCache{iip, T, use_both}) where {iip, T, use_both} = T function SciMLBase.__init( prob::BVProblem, alg::AbstractMIRK; dt = 0.0, abstol = 1e-6, adaptive = true, - controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), kwargs...) + controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), + optimize_kwargs = (; abstol = abstol), kwargs...) @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg) iip = isinplace(prob) diffcache = __cache_trait(alg.jac_alg) + @assert (iip || isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers " fit_parameters = haskey(prob.kwargs, :fit_parameters) t₀, t₁ = prob.tspan @@ -126,9 +129,9 @@ function SciMLBase.__init( return MIRKCache{iip, T, use_both, typeof(diffcache), fit_parameters}( alg_order(alg), stage, N, size(X), f, bc, prob_, prob.problem_type, prob.p, - alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, k_interp, - y, y₀, residual, fᵢ_cache, fᵢ₂_cache, errors, new_stages, resid₁_size, - nlsolve_kwargs, (; abstol, dt, adaptive, controller, fit_parameters, kwargs...)) + alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, k_interp, y, y₀, + residual, fᵢ_cache, fᵢ₂_cache, errors, new_stages, resid₁_size, nlsolve_kwargs, + optimize_kwargs, (; abstol, dt, adaptive, controller, fit_parameters, kwargs...)) end """ @@ -185,9 +188,11 @@ function SciMLBase.solve!(cache::MIRKCache{iip, T, use_both, diffcache, end function __perform_mirk_iteration(cache::MIRKCache, abstol, adaptive::Bool, controller::AbstractErrorControl) - nlprob = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀)) - nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve) - sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true) + nlprob = __construct_problem(cache, vec(cache.y₀), copy(cache.y₀)) + solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize) + kwargs = __concrete_kwargs( + cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs) + sol_nlprob = __internal_solve(nlprob, solve_alg; kwargs...) recursive_unflatten!(cache.y₀, sol_nlprob.u) error_norm = 2 * abstol @@ -200,7 +205,7 @@ function __perform_mirk_iteration(cache::MIRKCache, abstol, adaptive::Bool, cont if info == ReturnCode.Success # Nonlinear Solve was successful error_norm, info = error_estimate!( - cache, controller, cache.errors, sol_nlprob, nlsolve_alg, abstol) + cache, controller, cache.errors, sol_nlprob, solve_alg, abstol) end if info == ReturnCode.Success # Nonlinear Solve Successful and defect norm is acceptable @@ -233,7 +238,7 @@ function __perform_mirk_iteration(cache::MIRKCache, abstol, adaptive::Bool, cont end # Constructing the Nonlinear Problem -function __construct_nlproblem(cache::MIRKCache{iip}, y::AbstractVector, y₀::AbstractVectorOfArray) where {iip} +function __construct_problem(cache::MIRKCache{iip}, y::AbstractVector, y₀::AbstractVectorOfArray) where {iip} pt = cache.problem_type (; jac_alg) = cache.alg @@ -272,7 +277,7 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y::AbstractVector, y₀::A u, p, cache.y, pt, cache.bc, cache.mesh, cache, eval_sol, trait) end - return __construct_nlproblem(cache, y, loss_bc, loss_collocation, loss, pt) + return __construct_problem(cache, y, loss_bc, loss_collocation, loss, pt) end @views function __mirk_loss!(resid, u, p, y, pt::StandardBVProblem, bc!::BC, residual, @@ -378,7 +383,7 @@ end return mapreduce(vec, vcat, resids) end -function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collocation::C, +function __construct_problem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collocation::C, loss::LF, ::StandardBVProblem) where {iip, BC, C, LF} (; jac_alg) = cache.alg (; bc_diffmode) = jac_alg @@ -454,10 +459,8 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo end resid_prototype = vcat(resid_bc, resid_collocation) - nlf = NonlinearFunction{iip}( - loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype) - - return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p) + return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype, + resid_prototype, y, cache.p, cache.M, N) end function __mirk_mpoint_jacobian!( @@ -502,9 +505,9 @@ function __mirk_mpoint_jacobian( return J end -function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collocation::C, +function __construct_problem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collocation::C, loss::LF, ::TwoPointBVProblem) where {iip, BC, C, LF} - (; nlsolve, jac_alg) = cache.alg + (; jac_alg) = cache.alg N = length(cache.mesh) resid = vcat(@view(cache.bcresid_prototype[1:prod(cache.resid_size[1])]), @@ -545,9 +548,8 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo end resid_prototype = copy(resid) - nlf = NonlinearFunction{iip}( - loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype) - return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p) + return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype, + resid_prototype, y, cache.p, cache.M, N) end function __mirk_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid, p) where {L} diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl index 72043ca7..5c0a7247 100644 --- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl +++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl @@ -457,3 +457,43 @@ end @test sol.prob.p≈[17.09658] atol=1e-5 end + +@testitem "Convergence with optimization based solver" setup=[MIRKConvergenceTests] begin + using LinearAlgebra, DiffEqDevTools, OptimizationMOI, Ipopt + + # Only test on inplace problems + @testset "Problem: $i" for i in (3, 5, 9) + prob = probArr[i] + @testset "MIRK$order" for (_, order) in enumerate((2, 3, 4, 5, 6)) + sim = test_convergence( + dts, prob, mirk_solver(Val(order), optimize = Ipopt.Optimizer()); + abstol = 1e-8, reltol = 1e-8) + @test sim.𝒪est[:final]≈order atol=testTol + end + + @testset "MIRK$(order)I" for (_, order) in enumerate((6,)) + sim = test_convergence(dts, prob, MIRK6I(; optimize = Ipopt.Optimizer()); + abstol = 1e-8, reltol = 1e-8) + @test sim.𝒪est[:final]≈order atol=testTol + end + end +end + +@testitem "BVP with inequality constraints" begin + using BoundaryValueDiffEqMIRK, OptimizationMOI, Ipopt + + tspan = (0.0, pi / 2) + function simplependulum!(du, u, p, t) + θ = u[1] + dθ = u[2] + du[1] = dθ + du[2] = -9.81 * sin(θ) + end + function bc!(residual, u, p, t) + residual[1] = u(pi / 4)[1] + pi / 2 + residual[2] = u(pi / 2)[1] - pi / 2 + end + prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan, + lcons = [0.0, 0.0], ucons = [Inf, Inf]) + @test_nowarn sol = solve(prob, MIRK4(; optimize = Ipopt.Optimizer()), dt = 0.05) +end diff --git a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl index 234b3f05..fd90fa1d 100644 --- a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl +++ b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl @@ -6,9 +6,9 @@ using BandedMatrices: BandedMatrix, Ones using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, AbstractBoundaryValueDiffEqCache, BVPJacobianAlgorithm, recursive_flatten, recursive_flatten!, recursive_unflatten!, - __concrete_nonlinearsolve_algorithm, diff!, EvalSol, - eval_bc_residual, eval_bc_residual!, get_tmp, - __maybe_matmul!, __extract_problem_details, __initial_guess, + __concrete_solve_algorithm, diff!, EvalSol, eval_bc_residual, + eval_bc_residual!, get_tmp, __maybe_matmul!, + __extract_problem_details, __initial_guess, __maybe_allocate_diffcache, __restructure_sol, __get_bcresid_prototype, safe_similar, __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!, __vec_so_bc!, __vec_so_bc, @@ -19,7 +19,8 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, __default_sparse_ad, __default_nonsparse_ad, get_dense_ad, concrete_jacobian_algorithm, __default_coloring_algorithm, __default_sparsity_detector, interval, __split_kwargs, - NoErrorControl + NoErrorControl, __construct_internal_problem, + __concrete_kwargs using ConcreteStructs: @concrete using DiffEqBase: DiffEqBase diff --git a/lib/BoundaryValueDiffEqMIRKN/src/algorithms.jl b/lib/BoundaryValueDiffEqMIRKN/src/algorithms.jl index ede6be63..de36e6cf 100644 --- a/lib/BoundaryValueDiffEqMIRKN/src/algorithms.jl +++ b/lib/BoundaryValueDiffEqMIRKN/src/algorithms.jl @@ -6,7 +6,7 @@ for order in (4, 6) @eval begin """ - $($alg)(; nlsolve = NewtonRaphson(), jac_alg = BVPJacobianAlgorithm(), + $($alg)(; nlsolve = NewtonRaphson(), optimize = nothing, jac_alg = BVPJacobianAlgorithm(), defect_threshold = 0.1, max_num_subintervals = 3000) $($order)th order Monotonic Implicit Runge Kutta Nyström method. @@ -16,6 +16,10 @@ for order in (4, 6) - `nlsolve`: Internal Nonlinear solver. Any solver which conforms to the SciML `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. Note that any autodiff argument for + the solver will be ignored and a custom jacobian algorithm will be used. Optimization + solvers should first be loaded to allow this functionality. - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to use based on the input types and problem type. @@ -46,8 +50,9 @@ for order in (4, 6) } ``` """ - @kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractMIRKN + @kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractMIRKN nlsolve::N = nothing + optimize::O = nothing jac_alg::J = BVPJacobianAlgorithm() defect_threshold::T = 0.1 max_num_subintervals::Int = 3000 diff --git a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl index 18d825dd..05bdb2cd 100644 --- a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl +++ b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl @@ -21,16 +21,19 @@ fᵢ₂_cache resid_size nlsolve_kwargs + optimize_kwargs kwargs end Base.eltype(::MIRKNCache{iip, T}) where {iip, T} = T -function SciMLBase.__init(prob::SecondOrderBVProblem, alg::AbstractMIRKN; dt = 0.0, - adaptive = false, abstol = 1e-6, controller = NoErrorControl(), - nlsolve_kwargs = (; abstol = abstol), kwargs...) +function SciMLBase.__init(prob::SecondOrderBVProblem, alg::AbstractMIRKN; + dt = 0.0, adaptive = false, abstol = 1e-6, + controller = NoErrorControl(), nlsolve_kwargs = (; abstol = abstol), + optimize_kwargs = (; abstol = abstol), kwargs...) @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg) iip = isinplace(prob) + @assert (iip || isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers " t₀, t₁ = prob.tspan ig, T, M, Nig, X = __extract_problem_details(prob; dt, check_positive_dt = true) mesh = __extract_mesh(prob.u0, t₀, t₁, Nig) @@ -93,9 +96,10 @@ function SciMLBase.__init(prob::SecondOrderBVProblem, alg::AbstractMIRKN; dt = 0 prob_ = !(prob.u0 isa AbstractArray) ? remake(prob; u0 = X) : prob return MIRKNCache{iip, T}( - alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type, prob.p, alg, TU, - bcresid_prototype, mesh, mesh_dt, k_discrete, y, y₀, residual, fᵢ_cache, fᵢ₂_cache, - resid_size, nlsolve_kwargs, (; abstol, dt, adaptive, controller, kwargs...)) + alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type, + prob.p, alg, TU, bcresid_prototype, mesh, mesh_dt, k_discrete, + y, y₀, residual, fᵢ_cache, fᵢ₂_cache, resid_size, nlsolve_kwargs, + optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...)) end function SciMLBase.solve!(cache::MIRKNCache{iip, T}) where {iip, T} @@ -111,9 +115,11 @@ function SciMLBase.solve!(cache::MIRKNCache{iip, T}) where {iip, T} end function __perform_mirkn_iteration(cache::MIRKNCache) - nlprob::NonlinearProblem = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀)) - nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve) - sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true) + nlprob = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀)) + solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize) + kwargs = __concrete_kwargs( + cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs) + sol_nlprob = __solve(nlprob, solve_alg; kwargs...) recursive_unflatten!(cache.y₀, sol_nlprob.u) return sol_nlprob, sol_nlprob.retcode @@ -224,9 +230,8 @@ function __construct_nlproblem(cache::MIRKNCache{iip}, y, loss_bc::BC, loss_coll cache_collocation, loss_bc, loss_collocation, L, cache.p) end resid_prototype = vcat(resid_bc, resid_collocation) - nlf = NonlinearFunction{iip}( - loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype) - __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p) + return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype, + resid_prototype, y, cache.p, cache.M, 2 * N) end function __construct_nlproblem(cache::MIRKNCache{iip}, y, loss_bc::BC, loss_collocation::C, @@ -267,9 +272,8 @@ function __construct_nlproblem(cache::MIRKNCache{iip}, y, loss_bc::BC, loss_coll end resid_prototype = copy(resid) - nlf = NonlinearFunction{iip}( - loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype) - return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p) + return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype, + resid_prototype, y, cache.p, cache.M, 2 * N) end function __mirkn_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid, p) where {L} diff --git a/lib/BoundaryValueDiffEqShooting/Project.toml b/lib/BoundaryValueDiffEqShooting/Project.toml index ab0a608e..c7aa7966 100644 --- a/lib/BoundaryValueDiffEqShooting/Project.toml +++ b/lib/BoundaryValueDiffEqShooting/Project.toml @@ -10,13 +10,10 @@ BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a" ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" -FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e" FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" 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" @@ -36,7 +33,6 @@ ConcreteStructs = "0.2.3" DiffEqBase = "6.167" DiffEqDevTools = "2.44" DifferentiationInterface = "0.6.42, 0.7" -FastAlmostBandedMatrices = "0.1.4" FastClosures = "0.3.2" ForwardDiff = "0.10.38, 1" Hwloc = "3" @@ -49,8 +45,6 @@ OrdinaryDiffEqTsit5 = "1" OrdinaryDiffEqVerner = "1" Pkg = "1.10.0" PreallocationTools = "0.4.24" -PrecompileTools = "1.2" -Preferences = "1.4" Random = "1.10" ReTestItems = "1.23.1" RecursiveArrayTools = "3.27.0" diff --git a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl index 9292cc9f..d8d7faed 100644 --- a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl +++ b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl @@ -5,10 +5,11 @@ using ArrayInterface: fast_scalar_indexing using BandedMatrices: BandedMatrix, Ones using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm, recursive_flatten, recursive_flatten!, recursive_unflatten!, - __concrete_nonlinearsolve_algorithm, diff!, __any_sparse_ad, + __concrete_solve_algorithm, diff!, __any_sparse_ad, __cache_trait, concrete_jacobian_algorithm, eval_bc_residual, eval_bc_residual!, get_tmp, __maybe_matmul!, - __extract_problem_details, __initial_guess, + __concrete_kwargs, __extract_problem_details, + __initial_guess, __construct_internal_problem, __default_coloring_algorithm, __default_sparsity_detector, __maybe_allocate_diffcache, __get_bcresid_prototype, safe_similar, __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!, @@ -23,13 +24,9 @@ using ConcreteStructs: @concrete using DiffEqBase: DiffEqBase, solve using DifferentiationInterface: DifferentiationInterface, Constant, prepare_jacobian, overloaded_input_type -using FastAlmostBandedMatrices: AlmostBandedMatrix, fillpart, exclusive_bandpart, - finish_part_setindex! using FastClosures: @closure using ForwardDiff: ForwardDiff, pickchunksize using LinearAlgebra -using PrecompileTools: @compile_workload, @setup_workload -using Preferences: Preferences using Reexport: @reexport using RecursiveArrayTools: ArrayPartition, DiffEqArray, VectorOfArray using SciMLBase: SciMLBase, AbstractDiffEqInterpolation, StandardBVProblem, __solve, diff --git a/lib/BoundaryValueDiffEqShooting/src/algorithms.jl b/lib/BoundaryValueDiffEqShooting/src/algorithms.jl index baa20464..bef614c6 100644 --- a/lib/BoundaryValueDiffEqShooting/src/algorithms.jl +++ b/lib/BoundaryValueDiffEqShooting/src/algorithms.jl @@ -4,7 +4,7 @@ abstract type AbstractShooting <: AbstractBoundaryValueDiffEqAlgorithm end """ Shooting(ode_alg; kwargs...) Shooting(ode_alg, nlsolve; kwargs...) - Shooting(; ode_alg = nothing, nlsolve = nothing, jac_alg = nothing) + Shooting(; ode_alg = nothing, nlsolve = nothing, optimize = nothing, jac_alg = nothing) Single shooting method, reduces BVP to an initial value problem and solves the IVP. @@ -16,6 +16,9 @@ Single shooting method, reduces BVP to an initial value problem and solves the I - `nlsolve`: Internal Nonlinear solver. Any solver which conforms to the SciML `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver will be ignored and a custom jacobian algorithm will be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. Note that any autodiff argument for the solver + will be ignored and a custom jacobian algorithm will be used. - `jac_alg`: Jacobian Algorithm used for the Nonlinear Solver. If this is not set, we check if `nlsolve.ad` exists and is not nothing. If it is, we use that to construct the jacobian. If not, we try to use the best algorithm based on the input types @@ -25,11 +28,12 @@ Single shooting method, reduces BVP to an initial value problem and solves the I @concrete struct Shooting{J <: BVPJacobianAlgorithm} <: AbstractShooting ode_alg nlsolve + optimize jac_alg::J end -function Shooting(; ode_alg = nothing, nlsolve = nothing, jac_alg = nothing) - return Shooting(ode_alg, nlsolve, __materialize_jacobian_algorithm(nlsolve, jac_alg)) +function Shooting(; ode_alg = nothing, nlsolve = nothing, optimize = nothing, jac_alg = nothing) + return Shooting(ode_alg, nlsolve, optimize, __materialize_jacobian_algorithm(nlsolve, jac_alg)) end @inline Shooting(ode_alg; kwargs...) = Shooting(; ode_alg, kwargs...) @inline Shooting(ode_alg, nlsolve; kwargs...) = Shooting(; ode_alg, nlsolve, kwargs...) @@ -42,7 +46,7 @@ end """ MultipleShooting(; nshoots::Int, ode_alg = nothing, nlsolve = nothing, - grid_coarsening = true, jac_alg = nothing) + optimize = nothing, grid_coarsening = true, jac_alg = nothing) MultipleShooting(nshoots::Int; kwargs...) MultipleShooting(nshoots::Int, ode_alg; kwargs...) MultipleShooting(nshoots::Int, ode_alg, nlsolve; kwargs...) @@ -59,6 +63,8 @@ Significantly more stable than Single Shooting. poly-algorithm if `DifferentialEquations.jl` is loaded else this must be supplied) - `nlsolve`: Internal Nonlinear solver. Any solver which conforms to the SciML `NonlinearProblem` interface can be used. + - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML + `OptimizationProblem` interface can be used. - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to use based on the input types and problem type. @@ -85,6 +91,7 @@ Significantly more stable than Single Shooting. @concrete struct MultipleShooting{J <: BVPJacobianAlgorithm} <: AbstractShooting ode_alg nlsolve + optimize jac_alg::J nshoots::Int grid_coarsening @@ -93,17 +100,18 @@ end function concretize_jacobian_algorithm(alg::MultipleShooting, prob) jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg) return MultipleShooting( - alg.ode_alg, alg.nlsolve, jac_alg, alg.nshoots, alg.grid_coarsening) + alg.ode_alg, alg.nlsolve, alg.optimize, jac_alg, alg.nshoots, alg.grid_coarsening) end function update_nshoots(alg::MultipleShooting, nshoots::Int) return MultipleShooting( - alg.ode_alg, alg.nlsolve, alg.jac_alg, nshoots, alg.grid_coarsening) + alg.ode_alg, alg.nlsolve, alg.optimize, alg.jac_alg, nshoots, alg.grid_coarsening) end function MultipleShooting(; nshoots::Int, ode_alg = nothing, nlsolve = nothing, + optimize = nothing, grid_coarsening::Union{ Bool, Function, <:AbstractVector{<:Integer}, Tuple{Vararg{Integer}}} = true, jac_alg = nothing) @@ -112,9 +120,8 @@ function MultipleShooting(; nshoots::Int, sort!(grid_coarsening; rev = true) @assert all(grid_coarsening .> 0) && 1 ∉ grid_coarsening end - return MultipleShooting( - ode_alg, nlsolve, __materialize_jacobian_algorithm(nlsolve, jac_alg), - nshoots, grid_coarsening) + return MultipleShooting(ode_alg, nlsolve, optimize, + __materialize_jacobian_algorithm(nlsolve, jac_alg), nshoots, grid_coarsening) end @inline MultipleShooting(nshoots::Int; kwargs...) = MultipleShooting(; nshoots, kwargs...) @inline MultipleShooting(nshoots::Int, ode_alg; kwargs...) = MultipleShooting(; nshoots, ode_alg, kwargs...) diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl index 940aa299..a384f4ed 100644 --- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl +++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl @@ -1,5 +1,7 @@ -function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwargs = (;), - nlsolve_kwargs = (;), ensemblealg = EnsembleThreads(), verbose = true, kwargs...) +function SciMLBase.__solve( + prob::BVProblem, _alg::MultipleShooting; abstol = 1e-6, odesolve_kwargs = (;), + nlsolve_kwargs = (; abstol = abstol), optimize_kwargs = (; abstol = abstol), + ensemblealg = EnsembleThreads(), verbose = true, kwargs...) (; f, tspan) = prob if !(ensemblealg isa EnsembleSerial) && !(ensemblealg isa EnsembleThreads) @@ -14,6 +16,7 @@ function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwa bcresid_prototype, resid_size = __get_bcresid_prototype(prob, u0) iip, bc, u0, u0_size = isinplace(prob), prob.f.bc, deepcopy(u0), size(u0) + @assert (iip || isnothing(_alg.optimize)) "Out-of-place constraints don't allow optimization solvers " __alg = concretize_jacobian_algorithm(_alg, prob) alg = if has_initial_guess && Nig != __alg.nshoots @@ -64,14 +67,14 @@ function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwa if prob.problem_type isa TwoPointBVProblem __solve_nlproblem!(prob.problem_type, alg, bcresid_prototype, u_at_nodes, nodes, - cur_nshoot, M, N, resida_len, residb_len, solve_internal_odes!, - bc[1], bc[2], prob, u0, ode_cache_loss_fn, ensemblealg, - internal_ode_kwargs; verbose, kwargs..., nlsolve_kwargs...) + cur_nshoot, M, N, resida_len, residb_len, solve_internal_odes!, bc[1], + bc[2], prob, u0, ode_cache_loss_fn, ensemblealg, internal_ode_kwargs; + verbose, kwargs..., nlsolve_kwargs, optimize_kwargs...) else __solve_nlproblem!(prob.problem_type, alg, bcresid_prototype, u_at_nodes, nodes, - cur_nshoot, M, N, prod(resid_size), solve_internal_odes!, - bc, prob, f, u0_size, u0, ode_cache_loss_fn, ensemblealg, - internal_ode_kwargs; verbose, kwargs..., nlsolve_kwargs...) + cur_nshoot, M, N, prod(resid_size), solve_internal_odes!, bc, prob, + f, u0_size, u0, ode_cache_loss_fn, ensemblealg, internal_ode_kwargs; + verbose, kwargs..., nlsolve_kwargs, optimize_kwargs...) end end @@ -80,11 +83,12 @@ function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwa else diffmode_shooting = __get_non_sparse_ad(alg.jac_alg.bc_diffmode) end - shooting_alg = Shooting(alg.ode_alg, alg.nlsolve, BVPJacobianAlgorithm(diffmode_shooting)) + + shooting_alg = Shooting(alg.ode_alg, alg.nlsolve, alg.optimize, BVPJacobianAlgorithm(diffmode_shooting)) single_shooting_prob = remake(prob; u0 = reshape(u_at_nodes[1:N], u0_size)) return __solve(single_shooting_prob, shooting_alg; odesolve_kwargs, - nlsolve_kwargs, verbose, kwargs...) + nlsolve_kwargs, optimize_kwargs, verbose, kwargs...) end # TODO: We can save even more memory by hoisting the preallocated caches for the ODEs @@ -139,9 +143,12 @@ function __solve_nlproblem!( jac_prototype = jac_prototype) # NOTE: u_at_nodes is updated inplace - nlprob = __internal_nlsolve_problem(prob, M, N, loss_function!, u_at_nodes, prob.p) - nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, alg.nlsolve) - __solve(nlprob, nlsolve_alg; kwargs..., alias_u0 = true) + nlprob = __construct_internal_problem( + prob, alg, loss_fn, jac_fn, jac_prototype, resid_prototype, + u_at_nodes, prob.p, M, length(nodes), nothing) + + nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize) + __solve(nlprob, nlsolve_alg; kwargs...) return nothing end @@ -211,13 +218,12 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_ J, u, p, similar(bcresid_prototype), resid_nodes, ode_jac_cache, bc_jac_cache, ode_fn, bc_fn, nonbc_diffmode, bc_diffmode, N, M, __cache_trait(alg.jac_alg)) - loss_function! = NonlinearFunction{true}(loss_fn; resid_prototype = resid_prototype, - jac_prototype = jac_prototype, jac = jac_fn) - # NOTE: u_at_nodes is updated inplace - nlprob = __internal_nlsolve_problem(prob, M, N, loss_function!, u_at_nodes, prob.p) - nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, alg.nlsolve) - __solve(nlprob, nlsolve_alg; kwargs..., alias_u0 = true) + nlprob = __construct_internal_problem( + prob, alg, loss_fn, jac_fn, jac_prototype, resid_prototype, + u_at_nodes, prob.p, M, length(nodes), nothing) + nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize) + __solve(nlprob, nlsolve_alg; kwargs...) return nothing end diff --git a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl index 8290531e..c8e9a567 100644 --- a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl +++ b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl @@ -1,5 +1,6 @@ -function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;), - nlsolve_kwargs = (;), verbose = true, kwargs...) +function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; abstol = 1e-6, + odesolve_kwargs = (;), nlsolve_kwargs = (; abstol = abstol), + optimize_kwargs = (; abstol = abstol), verbose = true, kwargs...) # Setup the problem if prob.u0 isa AbstractArray{<:Number} u0 = prob.u0 @@ -13,6 +14,7 @@ function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (; bcresid_prototype, resid_size = __get_bcresid_prototype(prob, u0) iip, bc, u0, u0_size = isinplace(prob), prob.f.bc, deepcopy(u0), size(u0) + @assert (iip || isnothing(alg_.optimize)) "Out-of-place constraints don't allow optimization solvers " resid_prototype = __vec(bcresid_prototype) # Construct the residual function @@ -73,12 +75,12 @@ function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (; jac_prototype, u, jac_cache, diffmode, loss_fnₚ) end - nlf = NonlinearFunction{iip}(loss_fn; jac_prototype = jac_prototype, - resid_prototype = resid_prototype, jac = jac_fn) - nlprob = __internal_nlsolve_problem(prob, resid_prototype, u0, nlf, vec(u0), prob.p) - nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, alg.nlsolve) - nlsol::SciMLBase.NonlinearSolution = __solve( - nlprob, nlsolve_alg; nlsolve_kwargs..., verbose, kwargs...) + nlprob = __construct_internal_problem(prob, alg, loss_fn, jac_fn, jac_prototype, + resid_prototype, u0, prob.p, length(u0), 1) + solve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize) + kwargs = __concrete_kwargs(alg.nlsolve, alg.optimize, nlsolve_kwargs, optimize_kwargs) + #TODO: add verbose kwarg + nlsol = __solve(nlprob, solve_alg; kwargs...) # There is no way to reinit with the same cache with different cache. But not saving # the internal values gives a significant speedup. So we just create a new cache diff --git a/test/misc/default_solvers.jl b/test/misc/default_solvers.jl index 91cf9543..7c936b88 100644 --- a/test/misc/default_solvers.jl +++ b/test/misc/default_solvers.jl @@ -1,24 +1,24 @@ @testitem "Default Solvers" begin using BoundaryValueDiffEq, Test - + function f(du, u, p, t) (x, v) = u du[1] = v du[2] = -x end - + function bc!(resid, sol, p, t) resid[1] = sol[1][1] resid[2] = sol[end][1] - 1 end - + tspan = (0.0, 100.0) u0 = [0.0, 1.0] bvp = BVProblem(f, bc!, u0, tspan) resid_f = Array{Float64}(undef, 2) sol = solve(bvp, Shooting(Tsit5())) sol2 = solve(bvp) - + @test sol2.alg == Tsit5() @test all(sol.u .== sol2.u) -end \ No newline at end of file +end diff --git a/test/misc/manifolds_tests.jl b/test/misc/manifolds_tests.jl index 9ae55eb3..8bb061fa 100644 --- a/test/misc/manifolds_tests.jl +++ b/test/misc/manifolds_tests.jl @@ -65,7 +65,7 @@ for alg in algs if alg isa Shooting || alg isa MultipleShooting - sol = solve(bvp, alg) + sol = solve(bvp, alg; abstol = 1e-8) else sol = solve(bvp, alg; dt, abstol = 1e-8) end diff --git a/test/misc/qa_tests.jl b/test/misc/qa_tests.jl index e406bc0c..78dc3275 100644 --- a/test/misc/qa_tests.jl +++ b/test/misc/qa_tests.jl @@ -1,10 +1,8 @@ @testitem "Quality Assurance" begin using Aqua - Aqua.test_all(BoundaryValueDiffEq; - ambiguities = false, - piracies = (broken = false, - treat_as_own = [SciMLBase.BVProblem])) + Aqua.test_all(BoundaryValueDiffEq; ambiguities = false, + piracies = (broken = false, treat_as_own = [SciMLBase.BVProblem])) end @testitem "JET Package Test" begin