From 0d5d1109be43f43e2d57bc697431cb46bcdb682e Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 13 Oct 2024 10:37:29 +0200 Subject: [PATCH] Fix #259 --- .github/workflows/downgrade.yml | 4 +--- Project.toml | 16 +++++++------ docs/Project.toml | 2 +- docs/make.jl | 2 +- docs/src/manual/time_dependent.md | 2 +- docs/src/tutorials/limit_cycles.md | 2 +- examples/Project.toml | 5 +++- ext/ModelingToolkitExt.jl | 14 +++++++---- ext/TimeEvolution/ODEProblem.jl | 6 ++--- ext/TimeEvolution/TimeEvolution.jl | 2 +- ext/TimeEvolution/hysteresis_sweep.jl | 2 +- ext/TimeEvolution/plotting.jl | 8 +++++-- src/modules/FFTWExt.jl | 2 +- src/modules/LinearResponse/plotting.jl | 32 +++++++++++++------------- src/solve_homotopy.jl | 5 +++- test/ModelingToolkitExt.jl | 22 +++++++----------- test/SteadyStateDiffEqExt.jl | 13 ++++++----- test/hysteresis_sweep.jl | 2 +- test/runtests.jl | 3 --- test/time_evolution.jl | 2 +- 20 files changed, 77 insertions(+), 69 deletions(-) diff --git a/.github/workflows/downgrade.yml b/.github/workflows/downgrade.yml index 288d3a39..a5591d5c 100644 --- a/.github/workflows/downgrade.yml +++ b/.github/workflows/downgrade.yml @@ -19,9 +19,7 @@ jobs: strategy: matrix: version: - - '1.10' - # - '1.9' - # - 'nightly' + - 'min' os: - ubuntu-latest arch: diff --git a/Project.toml b/Project.toml index 0460dbed..dd607300 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "HarmonicBalance" uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e" authors = ["Jan Kosata ", "Javier del Pino ", "Orjan Ameye "] -version = "0.10.9" +version = "0.10.10" [deps] BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54" @@ -26,13 +26,13 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [weakdeps] ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" [extensions] ModelingToolkitExt = "ModelingToolkit" SteadyStateDiffEqExt = "SteadyStateDiffEq" -TimeEvolution = "OrdinaryDiffEq" +TimeEvolution = "OrdinaryDiffEqTsit5" [compat] Aqua = "0.8" @@ -48,10 +48,11 @@ HomotopyContinuation = "2.9" JET = "0.9.9" JLD2 = "0.4.48, 0.5" Latexify = "0.16" -ModelingToolkit = "9.34 - 9.40" +ModelingToolkit = "9.34" NonlinearSolve = "3.14" OrderedCollections = "1.6" -OrdinaryDiffEq = "v6.89" +OrdinaryDiffEqTsit5 = "1.1" +OrdinaryDiffEqRosenbrock = "1.1" Peaks = "0.5" Plots = "1.39" PrecompileTools = "1.2" @@ -68,10 +69,11 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Pkg", "Test", "OrdinaryDiffEq", "ModelingToolkit", "SteadyStateDiffEq", "NonlinearSolve", "ExplicitImports", "Aqua", "JET", "Documenter"] +test = ["Pkg", "Test", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqRosenbrock", "ModelingToolkit", "SteadyStateDiffEq", "NonlinearSolve", "ExplicitImports", "Aqua", "JET", "Documenter"] \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml index 3912bfc0..b6621bba 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,7 +5,7 @@ DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365" HarmonicBalance = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" diff --git a/docs/make.jl b/docs/make.jl index 8aebbbdb..38bf0a24 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,7 +7,7 @@ using DocumenterCitations # extentions using ModelingToolkit -using OrdinaryDiffEq +using OrdinaryDiffEqTsit5 using SteadyStateDiffEq bib = CitationBibliography( diff --git a/docs/src/manual/time_dependent.md b/docs/src/manual/time_dependent.md index e46ab381..28c08089 100644 --- a/docs/src/manual/time_dependent.md +++ b/docs/src/manual/time_dependent.md @@ -11,7 +11,7 @@ ParameterSweep ## Plotting ```@docs -HarmonicBalance.plot(::OrdinaryDiffEq.ODESolution, ::Any, ::HarmonicEquation) +HarmonicBalance.plot(::OrdinaryDiffEqTsit5.ODESolution, ::Any, ::HarmonicEquation) ``` ## Miscellaneous diff --git a/docs/src/tutorials/limit_cycles.md b/docs/src/tutorials/limit_cycles.md index e694a0dc..9c342545 100644 --- a/docs/src/tutorials/limit_cycles.md +++ b/docs/src/tutorials/limit_cycles.md @@ -101,7 +101,7 @@ According to Zambon et al., a limit cycle solution exists around $F_0 \cong 0.01 Let us try and simulate the limit cycle. We could in principle run a time-dependent simulation with a fixed value of $F_0$, but this would require a suitable initial condition. Instead, we will sweep $F_0$ upwards from a low starting value. To observe the dynamics just after the jump has occurred, we follow the sweep by a time interval where the system evolves under fixed parameters. ```@example lc -using OrdinaryDiffEq +using OrdinaryDiffEqTsit5 initial_state = result[1][1] T = 2e6 diff --git a/examples/Project.toml b/examples/Project.toml index 3912bfc0..451d5d4c 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -3,9 +3,12 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365" HarmonicBalance = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" diff --git a/ext/ModelingToolkitExt.jl b/ext/ModelingToolkitExt.jl index 1866f7ed..b5110701 100644 --- a/ext/ModelingToolkitExt.jl +++ b/ext/ModelingToolkitExt.jl @@ -4,7 +4,7 @@ export ODESystem, ODEProblem, SteadyStateProblem, NonlinearProblem using HarmonicBalance: HarmonicEquation, is_rearranged, rearrange_standard, get_variables, ParameterList -using Symbolics: simplify, Equation, substitute, Num, @variables, expand +using Symbolics: simplify, Equation, substitute, Num, @variables, expand, unwrap, arguments using ModelingToolkit: ModelingToolkit, ODESystem, @@ -31,9 +31,14 @@ function ModelingToolkit.ODESystem(eom::HarmonicEquation) eom = rearrange_standard(eom) end - slow_time = (@independent_variables T; T) - par_names = declare_parameter.(eom.parameters) vars = get_variables(eom) + slow_times = arguments.(unwrap.(vars)) + @assert all(isone.(length.(slow_times))) "Only one argument for the variables are allowed." + slow_time = unique(first.(slow_times)) + @assert isone(length(slow_time)) "The argument of the variables are not the same." + slow_time_ivp = @eval @independent_variables $(Symbol(first(slow_time))) + + par_names = declare_parameter.(eom.parameters) eqs = deepcopy(eom.equations) eqs = swapsides.(eqs) @@ -41,7 +46,8 @@ function ModelingToolkit.ODESystem(eom::HarmonicEquation) eqs = substitute(eqs, Dict(zip(eom.parameters, par_names))) # compute jacobian for performance - @mtkbuild sys = ODESystem(eqs, slow_time, vars, par_names) #mtk v9 need @mtkbuild + # ∨ mtk v9 need @mtkbuild + @mtkbuild sys = ODESystem(eqs, first(slow_time_ivp), vars, par_names) return sys end diff --git a/ext/TimeEvolution/ODEProblem.jl b/ext/TimeEvolution/ODEProblem.jl index 9306622e..05388ac2 100644 --- a/ext/TimeEvolution/ODEProblem.jl +++ b/ext/TimeEvolution/ODEProblem.jl @@ -1,4 +1,4 @@ -using OrdinaryDiffEq: ODEProblem, solve, ODESolution +using OrdinaryDiffEqTsit5: ODEProblem, solve, ODESolution """ ODEProblem( @@ -9,11 +9,11 @@ using OrdinaryDiffEq: ODEProblem, solve, ODESolution timespan::Tuple ) -Creates an ODEProblem object used by OrdinaryDiffEq.jl from the equations in `eom` to simulate time-evolution within `timespan`. +Creates an ODEProblem object used by OrdinaryDiffEqTsit5.jl from the equations in `eom` to simulate time-evolution within `timespan`. `fixed_parameters` must be a dictionary mapping parameters+variables to numbers (possible to use a solution index, e.g. solutions[x][y] for branch y of solution x). If `x0` is specified, it is used as an initial condition; otherwise the values from `fixed_parameters` are used. """ -function OrdinaryDiffEq.ODEProblem( +function OrdinaryDiffEqTsit5.ODEProblem( eom::HarmonicEquation, fixed_parameters; sweep::ParameterSweep=ParameterSweep(), diff --git a/ext/TimeEvolution/TimeEvolution.jl b/ext/TimeEvolution/TimeEvolution.jl index bf61d750..2e9e185e 100644 --- a/ext/TimeEvolution/TimeEvolution.jl +++ b/ext/TimeEvolution/TimeEvolution.jl @@ -2,7 +2,7 @@ module TimeEvolution using DocStringExtensions using Symbolics: Num, substitute, unwrap -using OrdinaryDiffEq: OrdinaryDiffEq +using OrdinaryDiffEqTsit5: OrdinaryDiffEqTsit5 using HarmonicBalance: HarmonicBalance, diff --git a/ext/TimeEvolution/hysteresis_sweep.jl b/ext/TimeEvolution/hysteresis_sweep.jl index b691d628..9034d523 100644 --- a/ext/TimeEvolution/hysteresis_sweep.jl +++ b/ext/TimeEvolution/hysteresis_sweep.jl @@ -59,7 +59,7 @@ function HarmonicBalance.follow_branch( end problem_t = ODEProblem(res.problem.eom, sol_dict; timespan=(0, tf)) - res_t = solve(problem_t, OrdinaryDiffEq.Tsit5(); saveat=tf) + res_t = solve(problem_t, OrdinaryDiffEqTsit5.Tsit5(); saveat=tf) # closest branch to final state followed_branch[i] = _closest_branch_index(res, res_t.u[end], next_index) diff --git a/ext/TimeEvolution/plotting.jl b/ext/TimeEvolution/plotting.jl index 46b2c932..5f40ee28 100644 --- a/ext/TimeEvolution/plotting.jl +++ b/ext/TimeEvolution/plotting.jl @@ -21,7 +21,11 @@ Parametric plot of f[1] against f[2] Also callable as plot! """ function Plots.plot( - soln::OrdinaryDiffEq.ODESolution, funcs, harm_eq::HarmonicEquation; add=false, kwargs... + soln::OrdinaryDiffEqTsit5.ODESolution, + funcs, + harm_eq::HarmonicEquation; + add=false, + kwargs..., ) # start a new plot if needed @@ -51,7 +55,7 @@ function Plots.plot( end end -function Plots.plot!(soln::OrdinaryDiffEq.ODESolution, varargs...; kwargs...) +function Plots.plot!(soln::OrdinaryDiffEqTsit5.ODESolution, varargs...; kwargs...) return plot(soln, varargs...; add=true, kwargs...) end diff --git a/src/modules/FFTWExt.jl b/src/modules/FFTWExt.jl index 95e5ba21..370cc4c2 100644 --- a/src/modules/FFTWExt.jl +++ b/src/modules/FFTWExt.jl @@ -24,7 +24,7 @@ function FFT(soln_u, soln_t; window=DSP.Windows.hanning) return (fft_u / length(fft_f), 2 * pi * fft_f) end -# function HarmonicBalance.FFT(soln::OrdinaryDiffEq.ODESolution; window=DSP.Windows.hanning) +# function HarmonicBalance.FFT(soln::OrdinaryDiffEqTsit5.ODESolution; window=DSP.Windows.hanning) # return HarmonicBalance.FFT(soln.u, soln.t; window=window) # end diff --git a/src/modules/LinearResponse/plotting.jl b/src/modules/LinearResponse/plotting.jl index 0feb4c77..444703cd 100644 --- a/src/modules/LinearResponse/plotting.jl +++ b/src/modules/LinearResponse/plotting.jl @@ -9,10 +9,10 @@ function get_jacobian_response( if show_progress bar = Progress( - length(CartesianIndices(C)), - 1, - "Diagonalizing the Jacobian for each solution ... ", - 50, + length(CartesianIndices(C)); + dt=1, + desc="Diagonalizing the Jacobian for each solution ... ", + barlen=50, ) end # evaluate the Jacobians for the different values of noise frequency Ω @@ -38,10 +38,10 @@ function get_jacobian_response( if show_progress bar = Progress( - length(CartesianIndices(C)), - 1, - "Diagonalizing the Jacobian for each solution ... ", - 50, + length(CartesianIndices(C)); + dt=1, + desc="Diagonalizing the Jacobian for each solution ... ", + barlen=50, ) end # evaluate the Jacobians for the different values of noise frequency Ω @@ -64,10 +64,10 @@ function get_linear_response( # note: this could be optimized by not grabbing the entire huge dictionary every time if show_progress bar = Progress( - length(C), - 1, - "Solving the linear response ODE for each solution and input frequency ... ", - 50, + length(C); + dt=1, + desc="Solving the linear response ODE for each solution and input frequency ... ", + barlen=50, ) end for j in findall(stable) @@ -92,10 +92,10 @@ function get_rotframe_jacobian_response( if show_progress bar = Progress( - length(C), - 1, - "Solving the linear response ODE for each solution and input frequency ... ", - 50, + length(C); + dt=1, + desc="Solving the linear response ODE for each solution and input frequency ... ", + barlen=50, ) end diff --git a/src/solve_homotopy.jl b/src/solve_homotopy.jl index 0b4fe854..e820c37e 100644 --- a/src/solve_homotopy.jl +++ b/src/solve_homotopy.jl @@ -342,7 +342,10 @@ function _get_raw_solution( result_full = Array{Vector{Any},1}(undef, length(parameter_values)) if show_progress bar = Progress( - length(parameter_values), 1, "Solving via $method homotopy ...", 50 + length(parameter_values); + dt=1, + desc="Solving via $method homotopy ...", + barlen=50, ) end for i in eachindex(parameter_values) # do NOT thread this diff --git a/test/ModelingToolkitExt.jl b/test/ModelingToolkitExt.jl index 33b37020..22e05866 100644 --- a/test/ModelingToolkitExt.jl +++ b/test/ModelingToolkitExt.jl @@ -1,5 +1,6 @@ using HarmonicBalance using ModelingToolkit +using ModelingToolkit: varmap_to_vars using Test @variables α ω ω0 F γ t x(t) @@ -9,19 +10,12 @@ diff_eq = DifferentialEquation( add_harmonic!(diff_eq, x, ω) # harmonic_eq = get_harmonic_equations(diff_eq) -ODESystem(harmonic_eq) +sys = ODESystem(harmonic_eq) +fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) +param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) -force = 0.01 -omega0 = 1.1 -alpha = 1.0 -gamma = 0.01 -fixed_nonlin = (α => alpha, ω0 => omega0, F => force, γ => gamma) -ω_span = (0.9, 1.5) -ω_range = range(ω_span..., 100) -varied = ω => ω_range +for p in string.([α, ω, ω0, F, γ]) + @test p ∈ string.(parameters(sys)) +end -param = ParameterList(merge(Dict(fixed_nonlin), Dict(ω => 1.1))) -varied = 1 => ω_range -x0 = [1.0, 0.0] - -ODEProblem(harmonic_eq, x0, (0, 100), param) +ODEProblem(harmonic_eq, [1.0, 0.0], (0, 100), param) diff --git a/test/SteadyStateDiffEqExt.jl b/test/SteadyStateDiffEqExt.jl index 42eddfa6..d1b7cf40 100644 --- a/test/SteadyStateDiffEqExt.jl +++ b/test/SteadyStateDiffEqExt.jl @@ -1,5 +1,10 @@ -using ModelingToolkit, SteadyStateDiffEq, OrdinaryDiffEq, LinearAlgebra, NonlinearSolve +using HarmonicBalance, + ModelingToolkit, + SteadyStateDiffEq, + OrdinaryDiffEqRosenbrock, + LinearAlgebra, + NonlinearSolve @testset "Steady state sweeps" begin @testset "one variable ODE" begin @@ -54,11 +59,7 @@ using ModelingToolkit, SteadyStateDiffEq, OrdinaryDiffEq, LinearAlgebra, Nonline @named model = ODESystem(eqs, t, [u1, v1], [α, ω, ω0, F, γ]) model = structural_simplify(model) - force = 0.01 - omega0 = 1.0 - alpha = 1.0 - gamma = 0.01 - param = [alpha, 1.2, omega0, force, gamma] + param = [α, ω, ω0, F, γ] .=> [1.0, 1.2, 1.0, 0.01, 0.01] x0 = [1.0, 0.0] prob_ss = SteadyStateProblem{true}(model, x0, param; jac=true) prob_np = NonlinearProblem(prob_ss) diff --git a/test/hysteresis_sweep.jl b/test/hysteresis_sweep.jl index 7987af6a..b6a8bb58 100644 --- a/test/hysteresis_sweep.jl +++ b/test/hysteresis_sweep.jl @@ -1,4 +1,4 @@ -using HarmonicBalance, OrdinaryDiffEq +using HarmonicBalance, OrdinaryDiffEqTsit5 @variables α ω ω0 F γ η t x(t); # declare constant variables and a function x(t) diff --git a/test/runtests.jl b/test/runtests.jl index 1716d5b8..c92402d8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,6 @@ using HarmonicBalance using Test -# WARNING: Method definition (::Type{Symbolics.Num})(Base.Complex{Symbolics.Num}) in module HarmonicBalance at E:\HarmonicBalance.jl\src\Symbolics_customised.jl:165 overwritten in module HarmonicBalance on the same line (check for duplicate calls to `include`). -# WARNING: Method definition (::Type{HomotopyContinuation.ModelKit.Variable})(Symbolics.Num) in module HC_wrapper at E:\HarmonicBalance.jl\src\modules\HC_wrapper\homotopy_interface.jl:2 overwritten in module HC_wrapper on the same line (check for duplicate calls to `include`). - using Random const SEED = 0xd8e5d8df Random.seed!(SEED) diff --git a/test/time_evolution.jl b/test/time_evolution.jl index e9163ef3..c0e14202 100644 --- a/test/time_evolution.jl +++ b/test/time_evolution.jl @@ -1,5 +1,5 @@ using HarmonicBalance -using Symbolics, OrdinaryDiffEq +using OrdinaryDiffEqTsit5 using Test @variables Ω, γ, λ, F, x, θ, η, α, ω0, ω, t, T, ψ