Fix #259

NonlinearOscillations · Oct 13, 2024 · 0d5d110 · 0d5d110
1 parent f12e4c9
commit 0d5d110
Show file tree

Hide file tree

Showing 20 changed files with 77 additions and 69 deletions.
diff --git 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
@@ -1,7 +1,7 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Jan Kosata <[email protected]>", "Javier del Pino <[email protected]>", "Orjan Ameye <[email protected]>"]
-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"]
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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,17 +31,23 @@ 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)
     eqs = simplify.(expand.(eqs))
     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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -1,5 +1,5 @@
 using HarmonicBalance
-using Symbolics, OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 using Test
 
 @variables Ω, γ, λ, F, x, θ, η, α, ω0, ω, t, T, ψ