From 1eeecc6fd34f052df1a9bfd627969de924fb8d2e Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Mon, 25 Dec 2023 18:00:12 -0500 Subject: [PATCH] Add JET Tests and disable Multiple Shooting tests for now --- Manifest.toml | 16 +- Project.toml | 5 +- src/algorithms.jl | 81 ++--- src/solve/single_shooting.jl | 9 +- test/misc/odeinterface.jl | 100 ------ test/misc/odeinterface_wrapper.jl | 69 +++-- test/runtests.jl | 2 +- test/shooting/basic_problems.jl | 378 +++++++++++++++++++++++ test/shooting/basic_tests.jl | 324 ------------------- test/shooting/nonlinear_least_squares.jl | 2 +- test/shooting/orbital.jl | 2 +- 11 files changed, 492 insertions(+), 496 deletions(-) delete mode 100644 test/misc/odeinterface.jl create mode 100644 test/shooting/basic_problems.jl delete mode 100644 test/shooting/basic_tests.jl diff --git a/Manifest.toml b/Manifest.toml index 6d6134b8..2464fce9 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,12 +2,12 @@ julia_version = "1.10.0-rc2" manifest_format = "2.0" -project_hash = "caa3a84a29ae467d2583d2084e7f50dd07a4148f" +project_hash = "0ac16cb78a3540d2a0e82de32da74fda24f340a4" [[deps.ADTypes]] -git-tree-sha1 = "332e5d7baeff8497b923b730b994fa480601efc7" +git-tree-sha1 = "41c37aa88889c171f1300ceac1313c06e891d245" uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b" -version = "0.2.5" +version = "0.2.6" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] @@ -587,10 +587,10 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" [[deps.NonlinearSolve]] -deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "EnumX", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LazyArrays", "LineSearches", "LinearAlgebra", "LinearSolve", "MaybeInplace", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SimpleNonlinearSolve", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"] -git-tree-sha1 = "e61a283ef89110668b83db24cd7de8868fb8133e" +deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "EnumX", "FastBroadcast", "FastClosures", "FiniteDiff", "ForwardDiff", "LazyArrays", "LineSearches", "LinearAlgebra", "LinearSolve", "MaybeInplace", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SimpleNonlinearSolve", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"] +git-tree-sha1 = "72b036b728461272ae1b1c3f7096cb4c319d8793" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -version = "3.3.0" +version = "3.4.0" [deps.NonlinearSolve.extensions] NonlinearSolveBandedMatricesExt = "BandedMatrices" @@ -599,6 +599,7 @@ version = "3.3.0" NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim" NonlinearSolveMINPACKExt = "MINPACK" NonlinearSolveNLsolveExt = "NLsolve" + NonlinearSolveSIAMFANLEquationsExt = "SIAMFANLEquations" NonlinearSolveSpeedMappingExt = "SpeedMapping" NonlinearSolveSymbolicsExt = "Symbolics" NonlinearSolveZygoteExt = "Zygote" @@ -610,6 +611,7 @@ version = "3.3.0" LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891" MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" + SIAMFANLEquations = "084e46ad-d928-497d-ad5e-07fa361a48c4" SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" @@ -775,7 +777,7 @@ version = "0.6.42" [[deps.SciMLBase]] deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FillArrays", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables", "TruncatedStacktraces"] -git-tree-sha1 = "50f2e6905c201e1212871824da13a802804835d8" +git-tree-sha1 = "b8f7a0807314cce87bf846ba5fd12c1b0ef512b7" repo-rev = "ap/nlls_bvp" repo-url = "https://github.com/SciML/SciMLBase.jl.git" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" diff --git a/Project.toml b/Project.toml index 4dbb133c..a3132353 100644 --- a/Project.toml +++ b/Project.toml @@ -44,9 +44,11 @@ ArrayInterface = "7" BandedMatrices = "1" ConcreteStructs = "0.2" DiffEqBase = "6.145" +DiffEqDevTools = "2.44" FastAlmostBandedMatrices = "0.1" FastClosures = "0.3" ForwardDiff = "0.10" +JET = "0.8" LinearAlgebra = "1.9" LinearSolve = "2.20" NonlinearSolve = "2.6.1, 3" @@ -74,6 +76,7 @@ julia = "1.9" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -83,4 +86,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["StaticArrays", "Random", "DiffEqDevTools", "OrdinaryDiffEq", "Test", "SafeTestsets", "ODEInterface", "Aqua", "LinearSolve", "RecursiveArrayTools"] +test = ["StaticArrays", "Random", "DiffEqDevTools", "OrdinaryDiffEq", "Test", "SafeTestsets", "ODEInterface", "Aqua", "LinearSolve", "RecursiveArrayTools", "JET"] diff --git a/src/algorithms.jl b/src/algorithms.jl index 81794496..12aaa6fc 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -1,10 +1,24 @@ # Algorithms abstract type BoundaryValueDiffEqAlgorithm <: SciMLBase.AbstractBVPAlgorithm end +abstract type AbstractShooting <: BoundaryValueDiffEqAlgorithm end abstract type AbstractMIRK <: BoundaryValueDiffEqAlgorithm end ## Disable the ugly verbose printing by default +@inline __modifier_text!(list, fieldname, field) = push!(list, "$fieldname = $(field)") +@inline __modifier_text!(list, fieldname, ::Nothing) = list +@inline __modifier_text!(list, fieldname, ::Missing) = list +@inline function __modifier_text!(list, fieldname, field::SciMLBase.AbstractODEAlgorithm) + push!(list, "$fieldname = $(__nameof(field))()") +end + function Base.show(io::IO, alg::BoundaryValueDiffEqAlgorithm) - print(io, "$(nameof(typeof(alg)))()") + print(io, "$(__nameof(alg))(") + modifiers = String[] + for field in fieldnames(typeof(alg)) + __modifier_text!(modifiers, field, getfield(alg, field)) + end + print(io, join(modifiers, ", ")) + print(io, ")") end """ @@ -28,7 +42,7 @@ Single shooting method, reduces BVP to an initial value problem and solves the I and problem type. If `BVPJacobianAlgorithm` is provided, only `diffmode` is used (defaults to `AutoForwardDiff` if possible else `AutoFiniteDiff`). """ -@concrete struct Shooting{J <: BVPJacobianAlgorithm} <: BoundaryValueDiffEqAlgorithm +@concrete struct Shooting{J <: BVPJacobianAlgorithm} <: AbstractShooting ode_alg nlsolve jac_alg::J @@ -40,25 +54,17 @@ end @inline Shooting(ode_alg; kwargs...) = Shooting(; ode_alg, kwargs...) @inline Shooting(ode_alg, nlsolve; kwargs...) = Shooting(; ode_alg, nlsolve, kwargs...) -function Base.show(io::IO, alg::Shooting) - print(io, "Shooting(") - modifiers = String[] - alg.nlsolve !== nothing && push!(modifiers, "nlsolve = $(alg.nlsolve)") - alg.jac_alg !== nothing && push!(modifiers, "jac_alg = $(alg.jac_alg)") - alg.ode_alg !== nothing && push!(modifiers, "ode_alg = $(__nameof(alg.ode_alg))()") - print(io, join(modifiers, ", ")) - print(io, ")") -end - @inline function concretize_jacobian_algorithm(alg::Shooting, prob) alg.jac_alg.diffmode === nothing && (return @set alg.jac_alg.diffmode = __default_nonsparse_ad(prob.u0)) return alg end - """ - MultipleShooting(nshoots::Int, ode_alg = nothing; nlsolve = nothing, - grid_coarsening = true, jac_alg = BVPJacobianAlgorithm()) + MultipleShooting(; nshoots::Int, ode_alg = nothing, nlsolve = nothing, + grid_coarsening = true, jac_alg = nothing) + MultipleShooting(nshoots::Int; kwargs...) + MultipleShooting(nshoots::Int, ode_alg; kwargs...) + MultipleShooting(nshoots::Int, ode_alg, nlsolve; kwargs...) Multiple Shooting method, reduces BVP to an initial value problem and solves the IVP. Significantly more stable than Single Shooting. @@ -66,16 +72,12 @@ Significantly more stable than Single Shooting. ## Arguments - `nshoots`: Number of shooting points. + - `ode_alg`: ODE algorithm to use for solving the IVP. Any solver which conforms to the SciML `ODEProblem` interface can be used! (Defaults to `nothing` which will use poly-algorithm if `DifferentialEquations.jl` is loaded else this must be supplied) - -## Keyword Arguments - - `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. - + `NonlinearProblem` 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. @@ -83,8 +85,8 @@ Significantly more stable than Single Shooting. + For `TwoPointBVProblem`, only `diffmode` is used (defaults to `AutoSparseForwardDiff` if possible else `AutoSparseFiniteDiff`). + For `BVProblem`, `bc_diffmode` and `nonbc_diffmode` are used. For `nonbc_diffmode` - defaults to `AutoSparseForwardDiff` if possible else `AutoSparseFiniteDiff`. For - `bc_diffmode`, defaults to `AutoForwardDiff` if possible else `AutoFiniteDiff`. + we default to `AutoSparseForwardDiff` if possible else `AutoSparseFiniteDiff`. For + `bc_diffmode`, we default to `AutoForwardDiff` if possible else `AutoFiniteDiff`. - `grid_coarsening`: Coarsening the multiple-shooting grid to generate a stable IVP solution. Possible Choices: @@ -97,13 +99,8 @@ Significantly more stable than Single Shooting. + `Function`: Takes the current number of shooting points and returns the next number of shooting points. For example, if `nshoots = 10` and `grid_coarsening = n -> n ÷ 2`, then the grid will be coarsened to `[5, 2]`. - -!!! note - - For type-stability, the chunksizes for ForwardDiff ADTypes in `BVPJacobianAlgorithm` - must be provided. """ -@concrete struct MultipleShooting{J <: BVPJacobianAlgorithm} +@concrete struct MultipleShooting{J <: BVPJacobianAlgorithm} <: AbstractShooting ode_alg nlsolve jac_alg::J @@ -111,6 +108,16 @@ Significantly more stable than Single Shooting. grid_coarsening end +# function Base.show(io::IO, alg::MultipleShooting) +# print(io, "MultipleShooting(") +# modifiers = String[] +# alg.nlsolve !== nothing && push!(modifiers, "nlsolve = $(alg.nlsolve)") +# alg.jac_alg !== nothing && push!(modifiers, "jac_alg = $(alg.jac_alg)") +# alg.ode_alg !== nothing && push!(modifiers, "ode_alg = $(__nameof(alg.ode_alg))()") +# print(io, join(modifiers, ", ")) +# print(io, ")") +# 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, @@ -122,18 +129,22 @@ function update_nshoots(alg::MultipleShooting, nshoots::Int) alg.grid_coarsening) end -function MultipleShooting(nshoots::Int, ode_alg = nothing; nlsolve = nothing, - grid_coarsening = true, jac_alg = BVPJacobianAlgorithm()) - @assert grid_coarsening isa Bool || grid_coarsening isa Function || - grid_coarsening isa AbstractVector{<:Integer} || - grid_coarsening isa NTuple{N, <:Integer} where {N} +function MultipleShooting(; nshoots::Int, ode_alg = nothing, nlsolve = nothing, + grid_coarsening::Union{Bool, Function, <:AbstractVector{<:Integer}, + Tuple{Vararg{Integer}}} = true, jac_alg = nothing) grid_coarsening isa Tuple && (grid_coarsening = Vector(grid_coarsening...)) if grid_coarsening isa AbstractVector sort!(grid_coarsening; rev = true) @assert all(grid_coarsening .> 0) && 1 ∉ grid_coarsening end - return MultipleShooting(ode_alg, nlsolve, jac_alg, nshoots, grid_coarsening) + return MultipleShooting(ode_alg, nlsolve, + __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...) +@inline MultipleShooting(nshoots::Int, ode_alg, nlsolve; kwargs...) = MultipleShooting(; + nshoots, ode_alg, nlsolve, kwargs...) for order in (2, 3, 4, 5, 6) alg = Symbol("MIRK$(order)") diff --git a/src/solve/single_shooting.jl b/src/solve/single_shooting.jl index 74f3d6f9..46fffb30 100644 --- a/src/solve/single_shooting.jl +++ b/src/solve/single_shooting.jl @@ -72,16 +72,15 @@ function __solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;), nlf = __unsafe_nonlinearfunction{iip}(loss_fn; jac_prototype, resid_prototype, jac = jac_fn) nlprob = __internal_nlsolve_problem(prob, resid_prototype, u0, nlf, vec(u0), prob.p) - opt = __solve(nlprob, alg.nlsolve; nlsolve_kwargs..., verbose, kwargs...) + nlsol = __solve(nlprob, alg.nlsolve; nlsolve_kwargs..., verbose, 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 - internal_prob_final = ODEProblem{iip}(prob.f, reshape(opt.u, u0_size), prob.tspan, + internal_prob_final = ODEProblem{iip}(prob.f, reshape(nlsol.u, u0_size), prob.tspan, prob.p) - sol = __solve(internal_prob_final, alg.ode_alg; actual_ode_kwargs...) + odesol = __solve(internal_prob_final, alg.ode_alg; actual_ode_kwargs...) - retcode = SciMLBase.successful_retcode(opt) ? sol.retcode : opt.retcode - return BVPSolution(sol; original = opt, retcode) + return SciMLBase.build_solution(prob, odesol, nlsol) end function __single_shooting_loss!(resid_, u0_, p, cache, bc::BC, u0_size, diff --git a/test/misc/odeinterface.jl b/test/misc/odeinterface.jl deleted file mode 100644 index 42e61f86..00000000 --- a/test/misc/odeinterface.jl +++ /dev/null @@ -1,100 +0,0 @@ -using Test, BoundaryValueDiffEq, LinearAlgebra, ODEInterface, Random, RecursiveArrayTools - -# Adaptation of https://github.com/luchr/ODEInterface.jl/blob/958b6023d1dabf775033d0b89c5401b33100bca3/examples/BasicExamples/ex7.jl -function ex7_f!(du, u, p, t) - ϵ = p[1] - u₁, λ = u - du[1] = (sin(t)^2 + λ * sin(t)^4 / u₁) / ϵ - du[2] = 0 - return nothing -end - -function ex7_2pbc1!(resa, ua, p) - resa[1] = ua[1] - 1 - return nothing -end - -function ex7_2pbc2!(resb, ub, p) - resb[1] = ub[1] - 1 - return nothing -end - -u0 = [0.5, 1.0] -p = [0.1] -tspan = (-π / 2, π / 2) - -tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan, p; - bcresid_prototype = (zeros(1), zeros(1))) -sol_bvpm2 = solve(tpprob, BVPM2(); dt = π / 20) - -@testset "BVPM2" begin - @info "Testing BVPM2" - - sol_bvpm2 = solve(tpprob, BVPM2(); dt = π / 20) - @test SciMLBase.successful_retcode(sol_bvpm2) - resid_f = (Array{Float64, 1}(undef, 1), Array{Float64, 1}(undef, 1)) - ex7_2pbc1!(resid_f[1], sol_bvpm2(tspan[1]), nothing) - ex7_2pbc2!(resid_f[2], sol_bvpm2(tspan[2]), nothing) - @test norm(resid_f, Inf) < 1e-6 -end - -# Just test that it runs. BVPSOL only works with linearly separable BCs. -@testset "BVPSOL" begin - @info "Testing BVPSOL" - - @info "BVPSOL with Vector{<:AbstractArray}" - - initial_u0 = [sol_bvpm2(t) .+ rand() for t in tspan[1]:(π / 20):tspan[2]] - tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan; - bcresid_prototype = (zeros(1), zeros(1))) - - # Just test that it runs. BVPSOL only works with linearly separable BCs. - sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) - - @info "BVPSOL with VectorOfArray" - - initial_u0 = VectorOfArray([sol_bvpm2(t) .+ rand() for t in tspan[1]:(π / 20):tspan[2]]) - tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan; - bcresid_prototype = (zeros(1), zeros(1))) - - # Just test that it runs. BVPSOL only works with linearly separable BCs. - sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) - - @info "BVPSOL with DiffEqArray" - - ts = collect(tspan[1]:(π / 20):tspan[2]) - initial_u0 = DiffEqArray([sol_bvpm2(t) .+ rand() for t in ts], ts) - tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan; - bcresid_prototype = (zeros(1), zeros(1))) - - sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) - - @info "BVPSOL with initial guess function" - - initial_u0 = (p, t) -> sol_bvpm2(t) .+ rand() - tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan, p; - bcresid_prototype = (zeros(1), zeros(1))) - sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) -end - -#= -@info "COLNEW" - -function f!(du, u, p, t) - du[1] = u[2] - du[2] = u[1] -end -function bca!(resid_a, u_a, p) - resid_a[1] = u_a[1] - 1 -end -function bcb!(resid_b, u_b, p) - resid_b[1] = u_b[1] -end - -fun = BVPFunction(f!, (bca!, bcb!), bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true)) -tspan = (0.0, 1.0) - -prob = TwoPointBVProblem(fun, [1.0, 0.0], tspan) -sol_colnew = solve(prob, COLNEW(), dt = 0.01) -@test SciMLBase.successful_retcode(sol_colnew) -=# diff --git a/test/misc/odeinterface_wrapper.jl b/test/misc/odeinterface_wrapper.jl index 2dfa85ff..42e61f86 100644 --- a/test/misc/odeinterface_wrapper.jl +++ b/test/misc/odeinterface_wrapper.jl @@ -1,4 +1,4 @@ -using Test, BoundaryValueDiffEq, LinearAlgebra, ODEInterface, Random +using Test, BoundaryValueDiffEq, LinearAlgebra, ODEInterface, Random, RecursiveArrayTools # Adaptation of https://github.com/luchr/ODEInterface.jl/blob/958b6023d1dabf775033d0b89c5401b33100bca3/examples/BasicExamples/ex7.jl function ex7_f!(du, u, p, t) @@ -25,33 +25,59 @@ tspan = (-π / 2, π / 2) tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan, p; bcresid_prototype = (zeros(1), zeros(1))) - -@info "BVPM2" - sol_bvpm2 = solve(tpprob, BVPM2(); dt = π / 20) -@test SciMLBase.successful_retcode(sol_bvpm2) -resid_f = (Array{Float64, 1}(undef, 1), Array{Float64, 1}(undef, 1)) -ex7_2pbc1!(resid_f[1], sol_bvpm2(tspan[1]), nothing) -ex7_2pbc2!(resid_f[2], sol_bvpm2(tspan[2]), nothing) -@test norm(resid_f, Inf) < 1e-6 -function ex7_f2!(du, u, p, t) - u₁, λ = u - du[1] = (sin(t)^2 + λ * sin(t)^4 / u₁) / 0.1 - du[2] = 0 - return nothing +@testset "BVPM2" begin + @info "Testing BVPM2" + + sol_bvpm2 = solve(tpprob, BVPM2(); dt = π / 20) + @test SciMLBase.successful_retcode(sol_bvpm2) + resid_f = (Array{Float64, 1}(undef, 1), Array{Float64, 1}(undef, 1)) + ex7_2pbc1!(resid_f[1], sol_bvpm2(tspan[1]), nothing) + ex7_2pbc2!(resid_f[2], sol_bvpm2(tspan[2]), nothing) + @test norm(resid_f, Inf) < 1e-6 end -@info "BVPSOL" +# Just test that it runs. BVPSOL only works with linearly separable BCs. +@testset "BVPSOL" begin + @info "Testing BVPSOL" -initial_u0 = [sol_bvpm2(t) .+ rand() for t in tspan[1]:(π / 20):tspan[2]] -tpprob = TwoPointBVProblem(ex7_f2!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan; - bcresid_prototype = (zeros(1), zeros(1))) + @info "BVPSOL with Vector{<:AbstractArray}" -# Just test that it runs. BVPSOL only works with linearly separable BCs. -# TODO: Implement appolo reentry example from ODEInterface.jl -sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) + initial_u0 = [sol_bvpm2(t) .+ rand() for t in tspan[1]:(π / 20):tspan[2]] + tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan; + bcresid_prototype = (zeros(1), zeros(1))) + + # Just test that it runs. BVPSOL only works with linearly separable BCs. + sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) + + @info "BVPSOL with VectorOfArray" + + initial_u0 = VectorOfArray([sol_bvpm2(t) .+ rand() for t in tspan[1]:(π / 20):tspan[2]]) + tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan; + bcresid_prototype = (zeros(1), zeros(1))) + + # Just test that it runs. BVPSOL only works with linearly separable BCs. + sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) + + @info "BVPSOL with DiffEqArray" + + ts = collect(tspan[1]:(π / 20):tspan[2]) + initial_u0 = DiffEqArray([sol_bvpm2(t) .+ rand() for t in ts], ts) + tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan; + bcresid_prototype = (zeros(1), zeros(1))) + + sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) + + @info "BVPSOL with initial guess function" + + initial_u0 = (p, t) -> sol_bvpm2(t) .+ rand() + tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan, p; + bcresid_prototype = (zeros(1), zeros(1))) + sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) +end +#= @info "COLNEW" function f!(du, u, p, t) @@ -71,3 +97,4 @@ tspan = (0.0, 1.0) prob = TwoPointBVProblem(fun, [1.0, 0.0], tspan) sol_colnew = solve(prob, COLNEW(), dt = 0.01) @test SciMLBase.successful_retcode(sol_colnew) +=# diff --git a/test/runtests.jl b/test/runtests.jl index 4b9d674a..3366d4eb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,4 +61,4 @@ const GROUP = uppercase(get(ENV, "GROUP", "ALL")) end end end -end \ No newline at end of file +end diff --git a/test/shooting/basic_problems.jl b/test/shooting/basic_problems.jl new file mode 100644 index 00000000..fb4950ca --- /dev/null +++ b/test/shooting/basic_problems.jl @@ -0,0 +1,378 @@ +using BoundaryValueDiffEq, LinearAlgebra, LinearSolve, OrdinaryDiffEq, Test, JET + +# FIXME: Reenable the Multiple Shooting Tests + +@testset "Basic Shooting Tests" begin + SOLVERS = [ + Shooting(Tsit5(), NewtonRaphson(; autodiff = AutoFiniteDiff())), + Shooting(Tsit5(), NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 2))), + Shooting(Tsit5()), + # MultipleShooting(10, Tsit5(), + # NewtonRaphson(; autodiff = AutoFiniteDiff())), + # MultipleShooting(10, Tsit5(), + # NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 2))), + # MultipleShooting(10, Tsit5()), + ] + JET_SKIP = [false, false, true, false, false, true] + JET_BROKEN = [false, false, false, false, false, false] + + tspan = (0.0, 100.0) + u0 = [0.0, 1.0] + + # Inplace + function f1!(du, u, p, t) + du[1] = u[2] + du[2] = -u[1] + return nothing + end + + function bc1!(resid, sol, p, t) + t₀, t₁ = first(t), last(t) + resid[1] = sol(t₀)[1] + resid[2] = sol(t₁)[1] - 1 + return nothing + end + + @info "Basic MP Shooting IIP" + + bvp1 = BVProblem(f1!, bc1!, u0, tspan) + @test SciMLBase.isinplace(bvp1) + + for (i, solver) in enumerate(SOLVERS) + @info "Solver: $solver" + sol = @time solve(bvp1, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3), maxiters = 10000) + + @test SciMLBase.successful_retcode(sol) + @test norm(sol.resid, Inf) < 1e-8 + + JET_SKIP[i] && continue + @test_opt target_modules=(SciMLBase, DiffEqBase, NonlinearSolve, + BoundaryValueDiffEq) solve(bvp1, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) broken=JET_BROKEN[i] + @test_call target_modules=(SciMLBase, DiffEqBase, NonlinearSolve, + BoundaryValueDiffEq) solve(bvp1, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) broken=JET_BROKEN[i] + end + + # Out of Place + f1(u, p, t) = [u[2], -u[1]] + + function bc1(sol, p, t) + t₀, t₁ = first(t), last(t) + return [sol(t₀)[1], sol(t₁)[1] - 1] + end + + @test_throws SciMLBase.NonconformingFunctionsError BVProblem(f1!, bc1, u0, tspan) + @test_throws SciMLBase.NonconformingFunctionsError BVProblem(f1, bc1!, u0, tspan) + + @info "Basic MP Shooting OOP" + + bvp2 = BVProblem(f1, bc1, u0, tspan) + @test !SciMLBase.isinplace(bvp2) + + for (i, solver) in enumerate(SOLVERS) + @info "Solver: $solver" + sol = @time solve(bvp2, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3), maxiters = 10000) + + @test SciMLBase.successful_retcode(sol) + @test norm(sol.resid, Inf) < 1e-8 + + JET_SKIP[i] && continue + @test_opt target_modules=(SciMLBase, DiffEqBase, NonlinearSolve, + BoundaryValueDiffEq) solve(bvp2, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) broken=JET_BROKEN[i] + @test_call target_modules=(SciMLBase, DiffEqBase, NonlinearSolve, + BoundaryValueDiffEq) solve(bvp2, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) broken=JET_BROKEN[i] + end + + # Inplace + bc2a!(resid, ua, p) = (resid[1] = ua[1]) + bc2b!(resid, ub, p) = (resid[1] = ub[1] - 1) + + @info "Basic TP Shooting IIP" + + bvp3 = TwoPointBVProblem(f1!, (bc2a!, bc2b!), u0, tspan; + bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1))) + @test SciMLBase.isinplace(bvp3) + + for (i, solver) in enumerate(SOLVERS) + @info "Solver: $solver" + sol = @time solve(bvp3, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3), maxiters = 10000) + + @test SciMLBase.successful_retcode(sol) + @test norm(sol.resid, Inf) < 1e-8 + + JET_SKIP[i] && continue + @test_opt target_modules=(SciMLBase, DiffEqBase, NonlinearSolve, + BoundaryValueDiffEq) solve(bvp3, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) broken=JET_BROKEN[i] + @test_call target_modules=(SciMLBase, DiffEqBase, NonlinearSolve, + BoundaryValueDiffEq) solve(bvp3, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) broken=JET_BROKEN[i] + end + + # Out of Place + bc2a(ua, p) = [ua[1]] + bc2b(ub, p) = [ub[1] - 1] + + @info "Basic TP Shooting OOP" + + bvp4 = TwoPointBVProblem(f1, (bc2a, bc2b), u0, tspan) + @test !SciMLBase.isinplace(bvp4) + + for (i, solver) in enumerate(SOLVERS) + @info "Solver: $solver" + sol = @time solve(bvp4, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3), maxiters = 10000) + + @test SciMLBase.successful_retcode(sol) + @test norm(sol.resid, Inf) < 1e-8 + + JET_SKIP[i] && continue + @test_opt target_modules=(SciMLBase, DiffEqBase, NonlinearSolve, + BoundaryValueDiffEq) solve(bvp4, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) broken=JET_BROKEN[i] + @test_call target_modules=(SciMLBase, DiffEqBase, NonlinearSolve, + BoundaryValueDiffEq) solve(bvp4, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) broken=JET_BROKEN[i] + end +end + +@testset "Shooting with Complex Values" begin + SOLVERS = [ + Shooting(Vern7(), NewtonRaphson(; autodiff = AutoFiniteDiff())), + Shooting(Vern7()), + # MultipleShooting(10, Vern7(), + # NewtonRaphson(; autodiff = AutoFiniteDiff())), + # MultipleShooting(10, Vern7()), + ] + JET_SKIP = [false, true, false, true] + JET_BROKEN = [false, false, false, false] + + function f1!(du, u, p, t) + du[1] = u[2] + du[2] = -u[1] + return nothing + end + + function bc1!(resid, sol, p, t) + t₀, t₁ = first(t), last(t) + resid[1] = sol(t₀)[1] + resid[2] = sol(t₁)[1] - 1 + return nothing + end + + tspan = (0.0, 100.0) + u0 = [0.0, 1.0] .+ 1im + bvp = BVProblem(f1!, bc1!, u0, tspan) + + @info "Shooting with Complex Values" + + for (i, solver) in enumerate(SOLVERS) + @info "Solver: $solver" + sol = @time solve(bvp, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3), maxiters = 10000) + + @test SciMLBase.successful_retcode(sol) + @test norm(sol.resid, Inf) < 1e-8 + end +end + +@testset "Flow In a Channel" begin + function flow_in_a_channel!(du, u, p, t) + R, P = p + A, f′′, f′, f, h′, h, θ′, θ = u + du[1] = 0 + du[2] = R * (f′^2 - f * f′′) - R * A + du[3] = f′′ + du[4] = f′ + du[5] = -R * f * h′ - 1 + du[6] = h′ + du[7] = -P * f * θ′ + du[8] = θ′ + end + + function bc_flow!(resid, sol, p, tspan) + t₁, t₂ = extrema(tspan) + solₜ₁ = sol(t₁) + solₜ₂ = sol(t₂) + resid[1] = solₜ₁[4] + resid[2] = solₜ₁[3] + resid[3] = solₜ₂[4] - 1 + resid[4] = solₜ₂[3] + resid[5] = solₜ₁[6] + resid[6] = solₜ₂[6] + resid[7] = solₜ₁[8] + resid[8] = solₜ₂[8] - 1 + end + + tspan = (0.0, 1.0) + p = [10.0, 7.0] + u0 = zeros(8) + + flow_bvp = BVProblem{true}(flow_in_a_channel!, bc_flow!, u0, tspan, p) + + @info "Flow in a Channel" + + for solver in [ + Shooting(AutoTsit5(Rosenbrock23())), + # MultipleShooting(10, AutoTsit5(Rosenbrock23())), + ] + @info "Solver: $solver" + sol = @time solve(flow_bvp, solver; abstol = 1e-8, reltol = 1e-8, + odesolve_kwargs = (; abstol = 1e-8, reltol = 1e-8)) + + @test SciMLBase.successful_retcode(sol) + @test norm(sol.resid, Inf) < 1e-8 + end +end + +# @testset "Ray Tracing" begin +# @inline v(x, y, z, p) = 1 / (4 + cos(p[1] * x) + sin(p[2] * y) - cos(p[3] * z)) +# @inline ux(x, y, z, p) = -p[1] * sin(p[1] * x) +# @inline uy(x, y, z, p) = p[2] * cos(p[2] * y) +# @inline uz(x, y, z, p) = p[3] * sin(p[3] * z) + +# function ray_tracing(u, p, t) +# du = similar(u) +# ray_tracing!(du, u, p, t) +# return du +# end + +# function ray_tracing!(du, u, p, t) +# x, y, z, ξ, η, ζ, T, S = u + +# nu = v(x, y, z, p) # Velocity of a sound wave, function of space; +# μx = ux(x, y, z, p) # ∂(slowness)/∂x, function of space +# μy = uy(x, y, z, p) # ∂(slowness)/∂y, function of space +# μz = uz(x, y, z, p) # ∂(slowness)/∂z, function of space + +# du[1] = S * nu * ξ +# du[2] = S * nu * η +# du[3] = S * nu * ζ + +# du[4] = S * μx +# du[5] = S * μy +# du[6] = S * μz + +# du[7] = S / nu +# du[8] = 0 + +# return nothing +# end + +# function ray_tracing_bc(sol, p, t) +# res = similar(first(sol)) +# ray_tracing_bc!(res, sol, p, t) +# return res +# end + +# function ray_tracing_bc!(res, sol, p, t) +# ua = sol(0.0) +# ub = sol(1.0) +# nu = v(ua[1], ua[2], ua[3], p) # Velocity of a sound wave, function of space; + +# res[1] = ua[1] - p[4] +# res[2] = ua[2] - p[5] +# res[3] = ua[3] - p[6] +# res[4] = ua[7] # T(0) = 0 +# res[5] = ua[4]^2 + ua[5]^2 + ua[6]^2 - 1 / nu^2 +# res[6] = ub[1] - p[7] +# res[7] = ub[2] - p[8] +# res[8] = ub[3] - p[9] +# return nothing +# end + +# function ray_tracing_bc_a(ua, p) +# resa = similar(ua, 5) +# ray_tracing_bc_a!(resa, ua, p) +# return resa +# end + +# function ray_tracing_bc_a!(resa, ua, p) +# nu = v(ua[1], ua[2], ua[3], p) # Velocity of a sound wave, function of space; + +# resa[1] = ua[1] - p[4] +# resa[2] = ua[2] - p[5] +# resa[3] = ua[3] - p[5] +# resa[4] = ua[7] +# resa[5] = ua[4]^2 + ua[5]^2 + ua[6]^2 - 1 / nu^2 + +# return nothing +# end + +# function ray_tracing_bc_b(ub, p) +# resb = similar(ub, 3) +# ray_tracing_bc_b!(resb, ub, p) +# return resb +# end + +# function ray_tracing_bc_b!(resb, ub, p) +# resb[1] = ub[1] - p[7] +# resb[2] = ub[2] - p[8] +# resb[3] = ub[3] - p[9] +# return nothing +# end + +# p = [0, 1, 2, 0, 0, 0, 4, 3, 2.0] + +# dx = p[7] - p[4] +# dy = p[8] - p[5] +# dz = p[9] - p[6] + +# u0 = zeros(8) +# u0[1:3] .= 0 # position +# u0[4] = dx / v(p[4], p[5], p[6], p) +# u0[5] = dy / v(p[4], p[5], p[6], p) +# u0[6] = dz / v(p[4], p[5], p[6], p) +# u0[8] = 1 + +# tspan = (0.0, 1.0) + +# prob_oop = BVProblem(BVPFunction{false}(ray_tracing, ray_tracing_bc), u0, tspan, p; +# nlls = Val(false)) +# prob_iip = BVProblem(BVPFunction{true}(ray_tracing!, ray_tracing_bc!), u0, tspan, p; +# nlls = Val(true)) +# prob_tp_oop = BVProblem(BVPFunction{false}(ray_tracing, +# (ray_tracing_bc_a, ray_tracing_bc_b); twopoint = Val(true)), u0, tspan, p; +# nlls = Val(true)) +# prob_tp_iip = BVProblem(BVPFunction{true}(ray_tracing!, +# (ray_tracing_bc_a!, ray_tracing_bc_b!); +# bcresid_prototype = (zeros(5), zeros(3)), twopoint = Val(true)), u0, tspan, p; +# nlls = Val(true)) + +# @info "Ray Tracing: Multiple Shooting" + +# alg_sp = MultipleShooting(10, AutoVern7(Rodas4P()); grid_coarsening = true, +# nlsolve = TrustRegion(), +# jac_alg = BVPJacobianAlgorithm(; bc_diffmode = AutoForwardDiff(; chunksize = 8), +# nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 8))) +# alg_dense = MultipleShooting(10, AutoVern7(Rodas4P()); grid_coarsening = true, +# nlsolve = TrustRegion(), +# jac_alg = BVPJacobianAlgorithm(; bc_diffmode = AutoForwardDiff(; chunksize = 8), +# nonbc_diffmode = AutoForwardDiff(; chunksize = 8))) +# alg_default = MultipleShooting(10, AutoVern7(Rodas4P()); grid_coarsening = true) + +# for (prob, alg) in Iterators.product((prob_oop, prob_iip, prob_tp_oop, prob_tp_iip), +# (alg_sp, alg_dense, alg_default)) +# @info "Solver: $alg" +# @time sol = solve(prob, alg; abstol = 1e-6, reltol = 1e-6, maxiters = 1000, +# odesolve_kwargs = (; abstol = 1e-8, reltol = 1e-5)) +# @test SciMLBase.successful_retcode(sol.retcode) + +# if prob.problem_type isa TwoPointBVProblem +# resida, residb = zeros(5), zeros(3) +# ray_tracing_bc_a!(resida, sol.u[1], p) +# ray_tracing_bc_b!(residb, sol.u[end], p) +# @test norm(vcat(resida, residb), Inf) < 1e-6 +# else +# resid = zeros(8) +# ray_tracing_bc!(resid, sol, p, sol.t) +# @test norm(resid, Inf) < 1e-6 +# end +# end +# end diff --git a/test/shooting/basic_tests.jl b/test/shooting/basic_tests.jl deleted file mode 100644 index 8f99d115..00000000 --- a/test/shooting/basic_tests.jl +++ /dev/null @@ -1,324 +0,0 @@ -using BoundaryValueDiffEq, LinearAlgebra, LinearSolve, OrdinaryDiffEq, Test - -@testset "Basic Shooting Tests" begin - SOLVERS = [Shooting(Tsit5()), MultipleShooting(10, Tsit5())] - - tspan = (0.0, 100.0) - u0 = [0.0, 1.0] - - # Inplace - function f1!(du, u, p, t) - du[1] = u[2] - du[2] = -u[1] - return nothing - end - - function bc1!(resid, sol, p, t) - t₀, t₁ = first(t), last(t) - resid[1] = sol(t₀)[1] - resid[2] = sol(t₁)[1] - 1 - return nothing - end - - @info "Basic MP Shooting IIP" - - bvp1 = BVProblem(f1!, bc1!, u0, tspan) - @test SciMLBase.isinplace(bvp1) - for solver in SOLVERS - @info "Solver: $solver" - resid_f = Array{Float64}(undef, 2) - sol = @time solve(bvp1, solver; abstol = 1e-13, reltol = 1e-13, - odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) - @test SciMLBase.successful_retcode(sol) - bc1!(resid_f, sol, nothing, sol.t) - @test norm(resid_f, Inf) < 1e-12 - end - - # Out of Place - f1(u, p, t) = [u[2], -u[1]] - - function bc1(sol, p, t) - t₀, t₁ = first(t), last(t) - return [sol(t₀)[1], sol(t₁)[1] - 1] - end - - @test_throws SciMLBase.NonconformingFunctionsError BVProblem(f1!, bc1, u0, tspan) - @test_throws SciMLBase.NonconformingFunctionsError BVProblem(f1, bc1!, u0, tspan) - - @info "Basic MP Shooting OOP" - - bvp2 = BVProblem(f1, bc1, u0, tspan) - @test !SciMLBase.isinplace(bvp2) - for solver in SOLVERS - @info "Solver: $solver" - sol = @time solve(bvp2, solver; abstol = 1e-13, reltol = 1e-13, - odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) - @test SciMLBase.successful_retcode(sol) - resid_f = bc1(sol, nothing, sol.t) - @test norm(resid_f, Inf) < 1e-12 - end - - # Inplace - bc2a!(resid, ua, p) = (resid[1] = ua[1]) - bc2b!(resid, ub, p) = (resid[1] = ub[1] - 1) - - @info "Basic TP Shooting IIP" - - bvp3 = TwoPointBVProblem(f1!, (bc2a!, bc2b!), u0, tspan; - bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1))) - @test SciMLBase.isinplace(bvp3) - for solver in SOLVERS - @info "Solver: $solver" - sol = @time solve(bvp3, solver; abstol = 1e-13, reltol = 1e-13, - odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) - @test SciMLBase.successful_retcode(sol) - resid_f = (Array{Float64, 1}(undef, 1), Array{Float64, 1}(undef, 1)) - bc2a!(resid_f[1], sol(tspan[1]), nothing) - bc2b!(resid_f[2], sol(tspan[2]), nothing) - @test norm(reduce(vcat, resid_f), Inf) < 1e-12 - end - - # Out of Place - bc2a(ua, p) = [ua[1]] - bc2b(ub, p) = [ub[1] - 1] - - @info "Basic TP Shooting OOP" - - bvp4 = TwoPointBVProblem(f1, (bc2a, bc2b), u0, tspan) - @test !SciMLBase.isinplace(bvp4) - for solver in SOLVERS - @info "Solver: $solver" - sol = @time solve(bvp4, solver; abstol = 1e-13, reltol = 1e-13, - odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-3)) - @test SciMLBase.successful_retcode(sol) - resid_f = reduce(vcat, (bc2a(sol(tspan[1]), nothing), bc2b(sol(tspan[2]), nothing))) - @test norm(resid_f, Inf) < 1e-12 - end -end - -@testset "Shooting with Complex Values" begin - function f1!(du, u, p, t) - du[1] = u[2] - du[2] = -u[1] - return nothing - end - - function bc1!(resid, sol, p, t) - t₀, t₁ = first(t), last(t) - resid[1] = sol(t₀)[1] - resid[2] = sol(t₁)[1] - 1 - return nothing - end - - tspan = (0.0, 100.0) - u0 = [0.0, 1.0] .+ 1im - bvp = BVProblem(f1!, bc1!, u0, tspan) - resid_f = Array{ComplexF64}(undef, 2) - - @info "Shooting with Complex Values" - - for solver in [Shooting(Vern7()), MultipleShooting(10, Vern7())] - @info "Solver: $solver" - sol = @time solve(bvp, solver; abstol = 1e-13, reltol = 1e-13, - odesolve_kwargs = (; abstol = 1e-8, reltol = 1e-8)) - @test SciMLBase.successful_retcode(sol) - bc1!(resid_f, sol, nothing, sol.t) - @test norm(resid_f, Inf) < 1e-12 - end -end - -@testset "Flow In a Channel" begin - function flow_in_a_channel!(du, u, p, t) - R, P = p - A, f′′, f′, f, h′, h, θ′, θ = u - du[1] = 0 - du[2] = R * (f′^2 - f * f′′) - R * A - du[3] = f′′ - du[4] = f′ - du[5] = -R * f * h′ - 1 - du[6] = h′ - du[7] = -P * f * θ′ - du[8] = θ′ - end - - function bc_flow!(resid, sol, p, tspan) - t₁, t₂ = extrema(tspan) - solₜ₁ = sol(t₁) - solₜ₂ = sol(t₂) - resid[1] = solₜ₁[4] - resid[2] = solₜ₁[3] - resid[3] = solₜ₂[4] - 1 - resid[4] = solₜ₂[3] - resid[5] = solₜ₁[6] - resid[6] = solₜ₂[6] - resid[7] = solₜ₁[8] - resid[8] = solₜ₂[8] - 1 - end - - tspan = (0.0, 1.0) - p = [10.0, 7.0] - u0 = zeros(8) - - flow_bvp = BVProblem{true}(flow_in_a_channel!, bc_flow!, u0, tspan, p) - - @info "Flow in a Channel" - - for solver in [ - Shooting(AutoTsit5(Rosenbrock23())), - MultipleShooting(10, AutoTsit5(Rosenbrock23())), - ] - @info "Solver: $solver" - sol = @time solve(flow_bvp, solver; abstol = 1e-13, reltol = 1e-13, - odesolve_kwargs = (; abstol = 1e-8, reltol = 1e-8)) - @test SciMLBase.successful_retcode(sol) - resid = zeros(8) - bc_flow!(resid, sol, p, sol.t) - @test norm(resid, Inf) < 1e-6 - end -end - -@testset "Ray Tracing" begin - @inline v(x, y, z, p) = 1 / (4 + cos(p[1] * x) + sin(p[2] * y) - cos(p[3] * z)) - @inline ux(x, y, z, p) = -p[1] * sin(p[1] * x) - @inline uy(x, y, z, p) = p[2] * cos(p[2] * y) - @inline uz(x, y, z, p) = p[3] * sin(p[3] * z) - - function ray_tracing(u, p, t) - du = similar(u) - ray_tracing!(du, u, p, t) - return du - end - - function ray_tracing!(du, u, p, t) - x, y, z, ξ, η, ζ, T, S = u - - nu = v(x, y, z, p) # Velocity of a sound wave, function of space; - μx = ux(x, y, z, p) # ∂(slowness)/∂x, function of space - μy = uy(x, y, z, p) # ∂(slowness)/∂y, function of space - μz = uz(x, y, z, p) # ∂(slowness)/∂z, function of space - - du[1] = S * nu * ξ - du[2] = S * nu * η - du[3] = S * nu * ζ - - du[4] = S * μx - du[5] = S * μy - du[6] = S * μz - - du[7] = S / nu - du[8] = 0 - - return nothing - end - - function ray_tracing_bc(sol, p, t) - res = similar(first(sol)) - ray_tracing_bc!(res, sol, p, t) - return res - end - - function ray_tracing_bc!(res, sol, p, t) - ua = sol(0.0) - ub = sol(1.0) - nu = v(ua[1], ua[2], ua[3], p) # Velocity of a sound wave, function of space; - - res[1] = ua[1] - p[4] - res[2] = ua[2] - p[5] - res[3] = ua[3] - p[6] - res[4] = ua[7] # T(0) = 0 - res[5] = ua[4]^2 + ua[5]^2 + ua[6]^2 - 1 / nu^2 - res[6] = ub[1] - p[7] - res[7] = ub[2] - p[8] - res[8] = ub[3] - p[9] - return nothing - end - - function ray_tracing_bc_a(ua, p) - resa = similar(ua, 5) - ray_tracing_bc_a!(resa, ua, p) - return resa - end - - function ray_tracing_bc_a!(resa, ua, p) - nu = v(ua[1], ua[2], ua[3], p) # Velocity of a sound wave, function of space; - - resa[1] = ua[1] - p[4] - resa[2] = ua[2] - p[5] - resa[3] = ua[3] - p[5] - resa[4] = ua[7] - resa[5] = ua[4]^2 + ua[5]^2 + ua[6]^2 - 1 / nu^2 - - return nothing - end - - function ray_tracing_bc_b(ub, p) - resb = similar(ub, 3) - ray_tracing_bc_b!(resb, ub, p) - return resb - end - - function ray_tracing_bc_b!(resb, ub, p) - resb[1] = ub[1] - p[7] - resb[2] = ub[2] - p[8] - resb[3] = ub[3] - p[9] - return nothing - end - - p = [0, 1, 2, 0, 0, 0, 4, 3, 2.0] - - dx = p[7] - p[4] - dy = p[8] - p[5] - dz = p[9] - p[6] - - u0 = zeros(8) - u0[1:3] .= 0 # position - u0[4] = dx / v(p[4], p[5], p[6], p) - u0[5] = dy / v(p[4], p[5], p[6], p) - u0[6] = dz / v(p[4], p[5], p[6], p) - u0[8] = 1 - - tspan = (0.0, 1.0) - - prob_oop = BVProblem(BVPFunction{false}(ray_tracing, ray_tracing_bc), u0, tspan, p; - nlls = Val(false)) - prob_iip = BVProblem(BVPFunction{true}(ray_tracing!, ray_tracing_bc!), u0, tspan, p; - nlls = Val(true)) - prob_tp_oop = BVProblem(BVPFunction{false}(ray_tracing, - (ray_tracing_bc_a, ray_tracing_bc_b); twopoint = Val(true)), u0, tspan, p; - nlls = Val(true)) - prob_tp_iip = BVProblem(BVPFunction{true}(ray_tracing!, - (ray_tracing_bc_a!, ray_tracing_bc_b!); - bcresid_prototype = (zeros(5), zeros(3)), twopoint = Val(true)), u0, tspan, p; - nlls = Val(true)) - - @info "Ray Tracing: Multiple Shooting" - - alg_sp = MultipleShooting(10, AutoVern7(Rodas4P()); grid_coarsening = true, - nlsolve = TrustRegion(), - jac_alg = BVPJacobianAlgorithm(; bc_diffmode = AutoForwardDiff(; chunksize = 8), - nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 8))) - alg_dense = MultipleShooting(10, AutoVern7(Rodas4P()); grid_coarsening = true, - nlsolve = TrustRegion(), - jac_alg = BVPJacobianAlgorithm(; bc_diffmode = AutoForwardDiff(; chunksize = 8), - nonbc_diffmode = AutoForwardDiff(; chunksize = 8))) - alg_default = MultipleShooting(10, AutoVern7(Rodas4P()); grid_coarsening = true) - - for (prob, alg) in Iterators.product((prob_oop, prob_iip, prob_tp_oop, prob_tp_iip), - (alg_sp, alg_dense, alg_default)) - @info "Solver: $alg" - @time sol = solve(prob, alg; abstol = 1e-6, reltol = 1e-6, maxiters = 1000, - odesolve_kwargs = (; abstol = 1e-8, reltol = 1e-5)) - @test SciMLBase.successful_retcode(sol.retcode) - - if prob.problem_type isa TwoPointBVProblem - resida, residb = zeros(5), zeros(3) - ray_tracing_bc_a!(resida, sol.u[1], p) - ray_tracing_bc_b!(residb, sol.u[end], p) - @test norm(vcat(resida, residb), Inf) < 1e-6 - else - resid = zeros(8) - ray_tracing_bc!(resid, sol, p, sol.t) - @test norm(resid, Inf) < 1e-6 - end - end -end \ No newline at end of file diff --git a/test/shooting/nonlinear_least_squares.jl b/test/shooting/nonlinear_least_squares.jl index 7fac9eb3..304cc805 100644 --- a/test/shooting/nonlinear_least_squares.jl +++ b/test/shooting/nonlinear_least_squares.jl @@ -91,4 +91,4 @@ using BoundaryValueDiffEq, LinearAlgebra, OrdinaryDiffEq, Test sol = @time solve(bvp4, solver; verbose = false) @test norm(sol.resid, Inf) < 1e-4 end -end \ No newline at end of file +end diff --git a/test/shooting/orbital.jl b/test/shooting/orbital.jl index 552e809c..8e2d5067 100644 --- a/test/shooting/orbital.jl +++ b/test/shooting/orbital.jl @@ -124,4 +124,4 @@ using BoundaryValueDiffEq, OrdinaryDiffEq, LinearAlgebra, Test @test SciMLBase.successful_retcode(sol) @test norm(reduce(vcat, resid_f_2p), Inf) < 1e-6 end -end \ No newline at end of file +end