From 01dac5b8e0bb489954661e7c05c3b1706d696b3e Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 26 Oct 2024 17:23:09 +0100 Subject: [PATCH 01/15] feat: `to_lab_frame` for multiple variables ansatz (#287) --- src/ExprUtils/Symbolics_utils.jl | 11 ++++ src/HC_wrapper/homotopy_interface.jl | 2 +- src/HarmonicBalance.jl | 12 ++++- src/HarmonicEquation.jl | 4 +- src/transform_solutions.jl | 77 +++++++++++++++------------- test/API.jl | 19 ++++++- test/symbolics.jl | 9 ++++ test/transform_solutions.jl | 28 ++++++---- 8 files changed, 111 insertions(+), 51 deletions(-) diff --git a/src/ExprUtils/Symbolics_utils.jl b/src/ExprUtils/Symbolics_utils.jl index 719897e1..a5560a1f 100644 --- a/src/ExprUtils/Symbolics_utils.jl +++ b/src/ExprUtils/Symbolics_utils.jl @@ -126,3 +126,14 @@ is_harmonic(x, t) = is_harmonic(Num(x), Num(t)) "Return true if `f` is a function of `var`." is_function(f, var) = any(isequal.(get_variables(f), var)) + +""" +Counts the number of derivatives of a symbolic variable. +""" +function count_derivatives(x::Symbolics.BasicSymbolic) + (Symbolics.isterm(x) || Symbolics.issym(x)) || + error("The input is not a single term or symbol") + bool = Symbolics.is_derivative(x) + return bool ? 1 + count_derivatives(first(arguments(x))) : 0 +end +count_derivatives(x::Num) = count_derivatives(Symbolics.unwrap(x)) diff --git a/src/HC_wrapper/homotopy_interface.jl b/src/HC_wrapper/homotopy_interface.jl index a5417b77..cde9bc61 100644 --- a/src/HC_wrapper/homotopy_interface.jl +++ b/src/HC_wrapper/homotopy_interface.jl @@ -61,7 +61,7 @@ function HarmonicBalance.Problem( system = HomotopyContinuation.System(eqs_HC; variables=conv_vars, parameters=conv_para) J = HarmonicBalance.get_Jacobian(equations, variables) #all derivatives are assumed to be in the left hand side; return Problem(variables, parameters, system, J) -end # is this funciton still needed/used? +end # TODO is this funciton still needed/used? function System(eom::HarmonicEquation) eqs = expand_derivatives.(_remove_brackets(eom)) diff --git a/src/HarmonicBalance.jl b/src/HarmonicBalance.jl index 20f74549..d1a0dee1 100644 --- a/src/HarmonicBalance.jl +++ b/src/HarmonicBalance.jl @@ -23,11 +23,19 @@ using Plots: Plots, plot, plot!, savefig, heatmap, Plot using Latexify: Latexify, latexify using Symbolics: - Symbolics, Num, Equation, @variables, expand_derivatives, get_variables, Differential + Symbolics, + Num, + Equation, + @variables, + expand_derivatives, + get_variables, + Differential, + unwrap, + wrap using SymbolicUtils: SymbolicUtils include("ExprUtils/ExprUtils.jl") -using .ExprUtils: is_harmonic, substitute_all, drop_powers +using .ExprUtils: is_harmonic, substitute_all, drop_powers, count_derivatives include("extention_functions.jl") include("utils.jl") diff --git a/src/HarmonicEquation.jl b/src/HarmonicEquation.jl index d5846309..c67940cf 100644 --- a/src/HarmonicEquation.jl +++ b/src/HarmonicEquation.jl @@ -266,7 +266,9 @@ function get_harmonic_equations( slow_time = isnothing(slow_time) ? (@variables T; T) : slow_time fast_time = isnothing(fast_time) ? get_independent_variables(diff_eom)[1] : fast_time - all(isempty.(values(diff_eom.harmonics))) && error("No harmonics specified!") + for pair in diff_eom.harmonics + isempty(pair[2]) && error("No harmonics specified for the variable $(pair[1])!") + end eom = harmonic_ansatz(diff_eom, fast_time) # substitute trig functions into the differential equation eom = slow_flow(eom; fast_time=fast_time, slow_time=slow_time, degree=degree) # drop 2nd order time derivatives fourier_transform!(eom, fast_time) # perform averaging over the frequencies originally specified in dEOM diff --git a/src/transform_solutions.jl b/src/transform_solutions.jl index b047e93f..2fa4660d 100644 --- a/src/transform_solutions.jl +++ b/src/transform_solutions.jl @@ -96,39 +96,39 @@ end # TRANSFORMATIONS TO THE LAB frame ### -function to_lab_frame(soln, res, times)::Vector{AbstractFloat} - timetrace = zeros(length(times)) - - for var in res.problem.eom.variables - val = real(substitute_all(Symbolics.unwrap(_remove_brackets(var)), soln)) - ω = real(Symbolics.unwrap(substitute_all(var.ω, soln))) - if var.type == "u" - timetrace .+= val * cos.(ω * times) - elseif var.type == "v" - timetrace .+= val * sin.(ω * times) - elseif var.type == "a" - timetrace .+= val - end - end - return timetrace -end - """ - to_lab_frame(res::Result, times; index::Int, branch::Int) - to_lab_frame(soln::OrderedDict, res::Result, times) + to_lab_frame(res::Result, x::Num, times; index::Int, branch::Int) + to_lab_frame(soln::OrderedDict, res::Result, nat_var::Num, times) -Transform a solution into the lab frame (i.e., invert the harmonic ansatz) for `times`. +Transform a solution into the lab frame (i.e., invert the harmonic ansatz) for the natural variable `x` for `times`. You can also compute the velocity by passing `d(x,t)` for the natural variable `x`. Either extract the solution from `res::Result` by `index` and `branch` or input `soln::OrderedDict` explicitly. """ -to_lab_frame(res::Result, times; index::Int, branch::Int) = - to_lab_frame(res[index][branch], res, times) +function to_lab_frame(soln::OrderedDict, res::Result, nat_var::Num, times) + count_derivatives(nat_var) > 2 && + throw(ArgumentError("Only the first derivative is supported")) -function to_lab_frame_velocity(soln, res, times) - timetrace = zeros(length(times)) + var = if Symbolics.is_derivative(unwrap(nat_var)) + wrap(first(Symbolics.arguments(unwrap(nat_var)))) + else + nat_var + end + vars = filter(x -> isequal(x.natural_variable, var), res.problem.eom.variables) + + return if Symbolics.is_derivative(unwrap(nat_var)) + _to_lab_frame_velocity(soln, vars, times) + else + _to_lab_frame(soln, vars, times) + end +end +function to_lab_frame(res::Result, nat_var::Num, times; index::Int, branch::Int) + return to_lab_frame(res[index][branch], res, nat_var, times) +end - for var in res.problem.eom.variables - val = real(substitute_all(Symbolics.unwrap(_remove_brackets(var)), soln)) - ω = real(real(Symbolics.unwrap(substitute_all(var.ω, soln)))) +function _to_lab_frame_velocity(soln, vars, times) + timetrace = zeros(length(times)) + for var in vars + val = real(substitute_all(unwrap(_remove_brackets(var)), soln)) + ω = real(real(unwrap(substitute_all(var.ω, soln)))) if var.type == "u" timetrace .+= -ω * val * sin.(ω * times) elseif var.type == "v" @@ -138,12 +138,19 @@ function to_lab_frame_velocity(soln, res, times) return timetrace end -""" - to_lab_frame_velocity(res::Result, times; index::Int, branch::Int) - to_lab_frame_velocity(soln::OrderedDict, res::Result, times) +function _to_lab_frame(soln, vars, times)::Vector{AbstractFloat} + timetrace = zeros(length(times)) -Transform a solution's velocity into the lab frame (i.e., invert the harmonic ansatz for dx/dt ) for `times`. -Either extract the solution from `res::Result` by `index` and `branch` or input `soln::OrderedDict` explicitly. -""" -to_lab_frame_velocity(res::Result, times; index, branch) = - to_lab_frame_velocity(res[index][branch], res, times) + for var in vars + val = real(substitute_all(unwrap(_remove_brackets(var)), soln)) + ω = real(unwrap(substitute_all(var.ω, soln))) + if var.type == "u" + timetrace .+= val * cos.(ω * times) + elseif var.type == "v" + timetrace .+= val * sin.(ω * times) + elseif var.type == "a" + timetrace .+= val + end + end + return timetrace +end diff --git a/test/API.jl b/test/API.jl index 60fd06fe..02696e1e 100644 --- a/test/API.jl +++ b/test/API.jl @@ -1,6 +1,6 @@ using HarmonicBalance -@testset begin +@testset "Equation{Vector}" begin # define equation of motion @variables ω1, ω2, t, ω, F, γ, α1, α2, k, x(t), y(t) rhs = [ @@ -27,3 +27,20 @@ end @test_throws MethodError get_steady_states(harmonic_eq, varied, threading=true) @test_throws ArgumentError get_steady_states(harmonic_eq, Dict(varied), threading=true) end + +@testset "forgot variable" begin + @variables Ω γ λ F x θ η α ω0 ω t T ψ + @variables x(t) y(t) + + natural_equation = [ + d(d(x, t), t) + γ * d(x, t) + Ω^2 * x + α * x^3 ~ F * cos(ω * t), + d(d(y, t), t) + γ * d(y, t) + Ω^2 * y + α * y^3 ~ 0, + ] + dEOM = DifferentialEquation(natural_equation, [x, y]) + + @test_throws ErrorException get_harmonic_equations(dEOM) + + add_harmonic!(dEOM, x, ω) + # add_harmonic!(dEOM, y, ω) + @test_throws ErrorException get_harmonic_equations(dEOM) +end diff --git a/test/symbolics.jl b/test/symbolics.jl index 99e5943a..6c173bab 100644 --- a/test/symbolics.jl +++ b/test/symbolics.jl @@ -177,3 +177,12 @@ end @eqtest get_independent(cos(t), t) == 0 @eqtest get_independent(cos(t)^2 + 5, t) == 5 end + +@testset "count_derivatives" begin + using HarmonicBalance.ExprUtils: count_derivatives + @variables t x(t) y(t) + @test count_derivatives(x) == 0 + @test count_derivatives(d(x, t)) == 1 + @test count_derivatives(d(d(x, t), t)) == 2 + @test_throws ErrorException count_derivatives(d(d(5 * x, t), t)) +end diff --git a/test/transform_solutions.jl b/test/transform_solutions.jl index 2b58a025..0411d61b 100644 --- a/test/transform_solutions.jl +++ b/test/transform_solutions.jl @@ -5,14 +5,17 @@ const SEED = 0xd8e5d8df Random.seed!(SEED) @variables Ω γ λ F x θ η α ω0 ω t T ψ -@variables x(t) +@variables x(t) y(t) -natural_equation = d(d(x, t), t) + γ * d(x, t) + Ω^2 * x + α * x^3 -forces = F * cos(ω * t) -dEOM = DifferentialEquation(natural_equation ~ forces, x) +natural_equation = [ + d(d(x, t), t) + γ * d(x, t) + Ω^2 * x + α * x^3 ~ F * cos(ω * t), + d(d(y, t), t) + γ * d(y, t) + Ω^2 * y + α * y^3 ~ 0, +] +dEOM = DifferentialEquation(natural_equation, [x, y]) add_harmonic!(dEOM, x, ω) -harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); +add_harmonic!(dEOM, y, ω) +harmonic_eq = get_harmonic_equations(dEOM); fixed = (Ω => 1.0, γ => 1e-2, F => 1e-3, α => 1.0) varied = ω => range(0.9, 1.1, 10) @@ -21,9 +24,12 @@ res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SE transform_solutions(res, "u1^2+v1^2") transform_solutions(res, "√(u1^2+v1^2)"; realify=true) -# transform_solutions(res.solutions[1][1], "√(u1^2+v1^2)", harmonic_eq) - -times = 0:1:10 -soln = res[5][1] -HarmonicBalance.to_lab_frame(soln, res, times) -HarmonicBalance.to_lab_frame_velocity(soln, res, times) +@testset "to_lab_frame" begin + using HarmonicBalance: to_lab_frame + @variables z(t) + times = 0:1:10 + @test to_lab_frame(res, x, times; index=1, branch=1) != zeros(length(times)) + @test all(isapprox.(to_lab_frame(res, y, times; index=1, branch=1), 0.0, atol=1e-10)) + @test all(to_lab_frame(res, z, times; index=1, branch=1) .≈ zeros(length(times))) + @test to_lab_frame(res, d(x, t), times; index=1, branch=1) != zeros(length(times)) +end From 26cf081ad8cf302f1fdb73564e67cf2ca3d5dbb8 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 26 Oct 2024 17:34:25 +0100 Subject: [PATCH 02/15] style: rename `ParameterSweep` to `AdiabaticSweep` (#289) --- docs/src/manual/time_dependent.md | 4 ++-- docs/src/tutorials/limit_cycles.md | 2 +- docs/src/tutorials/time_dependent.md | 4 ++-- ext/TimeEvolution/ODEProblem.jl | 4 ++-- ext/TimeEvolution/TimeEvolution.jl | 2 +- ext/TimeEvolution/sweeps.jl | 20 ++++++++++---------- src/HarmonicBalance.jl | 2 +- src/types.jl | 18 +++++++++--------- test/time_evolution.jl | 2 +- 9 files changed, 29 insertions(+), 29 deletions(-) diff --git a/docs/src/manual/time_dependent.md b/docs/src/manual/time_dependent.md index 28c08089..e86e8632 100644 --- a/docs/src/manual/time_dependent.md +++ b/docs/src/manual/time_dependent.md @@ -2,10 +2,10 @@ Generally, solving the ODE of oscillatory systems in time requires numerically tracking the oscillations. This is a computationally expensive process; however, using the harmonic ansatz removes the oscillatory time-dependence. Simulating instead the harmonic variables of a `HarmonicEquation` is vastly more efficient - a steady state of the system appears as a fixed point in multidimensional space rather than an oscillatory function. -The Extention `TimeEvolution` is used to interface `HarmonicEquation` with the solvers contained in `OrdinaryDiffEq.jl`. Time-dependent parameter sweeps are defined using the object `ParameterSweep`. To use the `TimeEvolution` extension, one must first load the `OrdinaryDiffEq.jl` package. +The Extention `TimeEvolution` is used to interface `HarmonicEquation` with the solvers contained in `OrdinaryDiffEq.jl`. Time-dependent parameter sweeps are defined using the object `AdiabaticSweep`. To use the `TimeEvolution` extension, one must first load the `OrdinaryDiffEq.jl` package. ```@docs ODEProblem(::HarmonicEquation, ::Any; timespan::Tuple) -ParameterSweep +AdiabaticSweep ``` ## Plotting diff --git a/docs/src/tutorials/limit_cycles.md b/docs/src/tutorials/limit_cycles.md index 8b265777..bdeaab84 100644 --- a/docs/src/tutorials/limit_cycles.md +++ b/docs/src/tutorials/limit_cycles.md @@ -105,7 +105,7 @@ using OrdinaryDiffEqTsit5 initial_state = result[1][1] T = 2e6 -sweep = ParameterSweep(F0 => (0.002, 0.011), (0,T)) +sweep = AdiabaticSweep(F0 => (0.002, 0.011), (0,T)) # start from initial_state, use sweep, total time is 2*T time_problem = ODEProblem(harmonic_eq, initial_state, sweep=sweep, timespan=(0,2*T)) diff --git a/docs/src/tutorials/time_dependent.md b/docs/src/tutorials/time_dependent.md index d6614b09..638272a3 100644 --- a/docs/src/tutorials/time_dependent.md +++ b/docs/src/tutorials/time_dependent.md @@ -76,9 +76,9 @@ Clearly when evolving from `x0 = [0.,0.]`, the system ends up in the low-amplitu Experimentally, the primary means of exploring the steady state landscape is an adiabatic sweep one or more of the system parameters. This takes the system along a solution branch. If this branch disappears or becomes unstable, a jump occurs. -The object [`ParameterSweep`](@ref ParameterSweep) specifies a sweep, which is then used as an optional `sweep` keyword in the `ODEProblem` constructor. +The object [`AdiabaticSweep`](@ref AdiabaticSweep) specifies a sweep, which is then used as an optional `sweep` keyword in the `ODEProblem` constructor. ```@example time_dependent -sweep = ParameterSweep(ω => (0.9,1.1), (0, 2e4)) +sweep = AdiabaticSweep(ω => (0.9,1.1), (0, 2e4)) ``` The sweep linearly interpolates between $\omega = 0.9$ at time 0 and $\omega = 1.1$ at time 2e4. For earlier/later times, $\omega$ is constant. diff --git a/ext/TimeEvolution/ODEProblem.jl b/ext/TimeEvolution/ODEProblem.jl index 05388ac2..859c4f5c 100644 --- a/ext/TimeEvolution/ODEProblem.jl +++ b/ext/TimeEvolution/ODEProblem.jl @@ -5,7 +5,7 @@ using OrdinaryDiffEqTsit5: ODEProblem, solve, ODESolution eom::HarmonicEquation; fixed_parameters, x0::Vector, - sweep::ParameterSweep, + sweep::AdiabaticSweep, timespan::Tuple ) @@ -16,7 +16,7 @@ If `x0` is specified, it is used as an initial condition; otherwise the values f function OrdinaryDiffEqTsit5.ODEProblem( eom::HarmonicEquation, fixed_parameters; - sweep::ParameterSweep=ParameterSweep(), + sweep::AdiabaticSweep=AdiabaticSweep(), x0::Vector=[], timespan::Tuple, perturb_initial=0.0, diff --git a/ext/TimeEvolution/TimeEvolution.jl b/ext/TimeEvolution/TimeEvolution.jl index c496c19e..bfb2d9a8 100644 --- a/ext/TimeEvolution/TimeEvolution.jl +++ b/ext/TimeEvolution/TimeEvolution.jl @@ -28,7 +28,7 @@ include("ODEProblem.jl") include("hysteresis_sweep.jl") include("plotting.jl") -export ParameterSweep +export AdiabaticSweep export transform_solutions, plot, plot! export ODEProblem, solve export follow_branch diff --git a/ext/TimeEvolution/sweeps.jl b/ext/TimeEvolution/sweeps.jl index 79ee3aeb..e99ba1d4 100644 --- a/ext/TimeEvolution/sweeps.jl +++ b/ext/TimeEvolution/sweeps.jl @@ -1,6 +1,6 @@ -using HarmonicBalance: ParameterSweep +using HarmonicBalance: AdiabaticSweep -function HarmonicBalance.ParameterSweep(functions::Dict, timespan::Tuple) +function HarmonicBalance.AdiabaticSweep(functions::Dict, timespan::Tuple) t0, t1 = timespan[1], timespan[2] sweep_func = Dict{Num,Any}([]) for swept_p in keys(functions) @@ -8,11 +8,11 @@ function HarmonicBalance.ParameterSweep(functions::Dict, timespan::Tuple) tfunc = swept_function(bounds, timespan) setindex!(sweep_func, tfunc, swept_p) end - return ParameterSweep(sweep_func) + return AdiabaticSweep(sweep_func) end -function HarmonicBalance.ParameterSweep(functions, timespan::Tuple) - return ParameterSweep(Dict(functions), timespan) +function HarmonicBalance.AdiabaticSweep(functions, timespan::Tuple) + return AdiabaticSweep(Dict(functions), timespan) end function swept_function(bounds, timespan) @@ -29,15 +29,15 @@ function swept_function(bounds, timespan) return f end -# overload so that ParameterSweep can be accessed like a Dict -Base.keys(s::ParameterSweep) = keys(s.functions) -Base.getindex(s::ParameterSweep, i) = getindex(s.functions, i) +# overload so that AdiabaticSweep can be accessed like a Dict +Base.keys(s::AdiabaticSweep) = keys(s.functions) +Base.getindex(s::AdiabaticSweep, i) = getindex(s.functions, i) # overload + -function Base.:+(s1::ParameterSweep, s2::ParameterSweep) +function Base.:+(s1::AdiabaticSweep, s2::AdiabaticSweep) common_params = intersect(keys(s1), keys(s2)) !isempty(common_params) && error("cannot combine sweeps of the same parameter") - return ParameterSweep(merge(s1.functions, s2.functions)) + return AdiabaticSweep(merge(s1.functions, s2.functions)) # combine sweeps of the same parameter #interval_overlap(s1.timespan, s2.timespan) && error("cannot combine sweeps with overlapping timespans") diff --git a/src/HarmonicBalance.jl b/src/HarmonicBalance.jl index d1a0dee1..9d0cb4fd 100644 --- a/src/HarmonicBalance.jl +++ b/src/HarmonicBalance.jl @@ -61,7 +61,7 @@ export rearrange_standard export plot, plot!, plot_phase_diagram, savefig, plot_spaghetti -export ParameterSweep, steady_state_sweep +export AdiabaticSweep, steady_state_sweep export plot_1D_solutions_branch, follow_branch include("HC_wrapper/HC_wrapper.jl") diff --git a/src/types.jl b/src/types.jl index ec2b7f11..04276447 100644 --- a/src/types.jl +++ b/src/types.jl @@ -284,36 +284,36 @@ $(TYPEDFIELDS) ```julia-repl # create a sweep of parameter a from 0 to 1 over time 0 -> 100 julia> @variables a,b; -julia> sweep = ParameterSweep(a => [0., 1.], (0, 100)); +julia> sweep = AdiabaticSweep(a => [0., 1.], (0, 100)); julia> sweep[a](50) 0.5 julia> sweep[a](200) 1.0 # do the same, varying two parameters simultaneously -julia> sweep = ParameterSweep([a => [0.,1.], b => [0., 1.]], (0,100)) +julia> sweep = AdiabaticSweep([a => [0.,1.], b => [0., 1.]], (0,100)) ``` Successive sweeps can be combined, ```julia-repl -sweep1 = ParameterSweep(ω => [0.95, 1.0], (0, 2e4)) -sweep2 = ParameterSweep(λ => [0.05, 0.01], (2e4, 4e4)) +sweep1 = AdiabaticSweep(ω => [0.95, 1.0], (0, 2e4)) +sweep2 = AdiabaticSweep(λ => [0.05, 0.01], (2e4, 4e4)) sweep = sweep1 + sweep2 ``` multiple parameters can be swept simultaneously, ```julia-repl -sweep = ParameterSweep([ω => [0.95;1.0], λ => [5e-2;1e-2]], (0, 2e4)) +sweep = AdiabaticSweep([ω => [0.95;1.0], λ => [5e-2;1e-2]], (0, 2e4)) ``` and custom sweep functions may be used. ```julia-repl ωfunc(t) = cos(t) -sweep = ParameterSweep(ω => ωfunc) +sweep = AdiabaticSweep(ω => ωfunc) ``` """ -struct ParameterSweep +struct AdiabaticSweep """Maps each swept parameter to a function.""" functions::Dict{Num,Function} - ParameterSweep(functions...) = new(Dict(functions...)) - ParameterSweep() = ParameterSweep([]) + AdiabaticSweep(functions...) = new(Dict(functions...)) + AdiabaticSweep() = AdiabaticSweep([]) end diff --git a/test/time_evolution.jl b/test/time_evolution.jl index c0e14202..0b3cac86 100644 --- a/test/time_evolution.jl +++ b/test/time_evolution.jl @@ -22,7 +22,7 @@ varied = ω => range(0.9, 1.1, 10) res = get_steady_states(p, varied, fixed; show_progress=false, seed=SEED) tspan = (0.0, 10) -sweep = ParameterSweep(ω => (0.9, 1.1), tspan) # linearly interpolate between two values at two times +sweep = AdiabaticSweep(ω => (0.9, 1.1), tspan) # linearly interpolate between two values at two times ode_problem = ODEProblem(harmonic_eq, fixed; sweep=sweep, x0=[0.01; 0.0], timespan=tspan) time_soln = solve(ode_problem, Tsit5(); saveat=1); From d91edcb30e4357ccf5762ca1bbff818042992ee3 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 26 Oct 2024 17:35:23 +0100 Subject: [PATCH 03/15] style: rename `x0` to be `u0` in line with `SciML` (#290) --- docs/src/tutorials/time_dependent.md | 14 +++++++------- ext/TimeEvolution/ODEProblem.jl | 16 ++++++++-------- test/SteadyStateDiffEqExt.jl | 4 ++-- test/time_evolution.jl | 4 ++-- 4 files changed, 19 insertions(+), 19 deletions(-) diff --git a/docs/src/tutorials/time_dependent.md b/docs/src/tutorials/time_dependent.md index 638272a3..6065960a 100644 --- a/docs/src/tutorials/time_dependent.md +++ b/docs/src/tutorials/time_dependent.md @@ -43,10 +43,10 @@ Given $\mathbf{u}(T_0)$, what is $\mathbf{u}(T)$ at future times? For constant parameters, a [`HarmonicEquation`](@ref HarmonicBalance.HarmonicEquation) object can be fed into the constructor of [`ODEProblem`](@ref ODEProblem). The syntax is similar to DifferentialEquations.jl : ```@example time_dependent using OrdinaryDiffEqTsit5 -x0 = [0.; 0.] # initial condition +u0 = [0.; 0.] # initial condition fixed = (ω0 => 1.0, γ => 1e-2, λ => 5e-2, F => 1e-3, α => 1.0, η => 0.3, θ => 0, ω => 1.0) # parameter values -ode_problem = ODEProblem(harmonic_eq, fixed, x0 = x0, timespan = (0,1000)) +ode_problem = ODEProblem(harmonic_eq, fixed, u0 = u0, timespan = (0,1000)) ``` OrdinaryDiffEq.jl takes it from here - we only need to use `solve`. @@ -55,10 +55,10 @@ time_evo = solve(ode_problem, Tsit5(), saveat=1.0); plot(time_evo, ["u1", "v1"], harmonic_eq) ``` -Running the above code with `x0 = [0.2, 0.2]` gives the plots +Running the above code with `u0 = [0.2, 0.2]` gives the plots ```@example time_dependent -x0 = [0.2; 0.2] # initial condition -ode_problem = remake(ode_problem, u0 = x0) +u0 = [0.2; 0.2] # initial condition +ode_problem = remake(ode_problem, u0 = u0) time_evo = solve(ode_problem, Tsit5(), saveat=1.0); plot(time_evo, ["u1", "v1"], harmonic_eq) ``` @@ -70,7 +70,7 @@ result = get_steady_states(harmonic_eq, varied, fixed) plot(result, "sqrt(u1^2 + v1^2)") ``` -Clearly when evolving from `x0 = [0.,0.]`, the system ends up in the low-amplitude branch 2. With `x0 = [0.2, 0.2]`, the system ends up in branch 3. +Clearly when evolving from `u0 = [0., 0.]`, the system ends up in the low-amplitude branch 2. With `u0 = [0.2, 0.2]`, the system ends up in branch 3. ## Adiabatic parameter sweeps @@ -84,7 +84,7 @@ The sweep linearly interpolates between $\omega = 0.9$ at time 0 and $\omega = Let us now define a new `ODEProblem` which incorporates `sweep` and again use `solve`: ```@example time_dependent -ode_problem = ODEProblem(harmonic_eq, fixed, sweep=sweep, x0=[0.1;0.0], timespan=(0, 2e4)) +ode_problem = ODEProblem(harmonic_eq, fixed, sweep=sweep, u0=[0.1;0.0], timespan=(0, 2e4)) time_evo = solve(ode_problem, Tsit5(), saveat=100) plot(time_evo, "sqrt(u1^2 + v1^2)", harmonic_eq) ``` diff --git a/ext/TimeEvolution/ODEProblem.jl b/ext/TimeEvolution/ODEProblem.jl index 859c4f5c..ec07993a 100644 --- a/ext/TimeEvolution/ODEProblem.jl +++ b/ext/TimeEvolution/ODEProblem.jl @@ -4,20 +4,20 @@ using OrdinaryDiffEqTsit5: ODEProblem, solve, ODESolution ODEProblem( eom::HarmonicEquation; fixed_parameters, - x0::Vector, - sweep::AdiabaticSweep, + u0::Vector, + sweep::ParameterSweep, timespan::Tuple ) 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. +If `u0` is specified, it is used as an initial condition; otherwise the values from `fixed_parameters` are used. """ function OrdinaryDiffEqTsit5.ODEProblem( eom::HarmonicEquation, fixed_parameters; - sweep::AdiabaticSweep=AdiabaticSweep(), - x0::Vector=[], + sweep::ParameterSweep=ParameterSweep(), + u0::Vector=[], timespan::Tuple, perturb_initial=0.0, kwargs..., @@ -54,11 +54,11 @@ function OrdinaryDiffEqTsit5.ODEProblem( end end - # the initial condition is x0 if specified, taken from fixed_parameters otherwise - initial = if isempty(x0) + # the initial condition is u0 if specified, taken from fixed_parameters otherwise + initial = if isempty(u0) real.(collect(values(fixed_parameters))[1:length(vars)]) * (1 - perturb_initial) else - x0 + u0 end return ODEProblem(f!, initial, timespan; kwargs...) diff --git a/test/SteadyStateDiffEqExt.jl b/test/SteadyStateDiffEqExt.jl index d1b7cf40..641f10ba 100644 --- a/test/SteadyStateDiffEqExt.jl +++ b/test/SteadyStateDiffEqExt.jl @@ -60,8 +60,8 @@ using HarmonicBalance, model = structural_simplify(model) 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) + u0 = [1.0, 0.0] + prob_ss = SteadyStateProblem{true}(model, u0, param; jac=true) prob_np = NonlinearProblem(prob_ss) ω_span = (0.9, 1.5) diff --git a/test/time_evolution.jl b/test/time_evolution.jl index 0b3cac86..e076bf95 100644 --- a/test/time_evolution.jl +++ b/test/time_evolution.jl @@ -22,8 +22,8 @@ varied = ω => range(0.9, 1.1, 10) res = get_steady_states(p, varied, fixed; show_progress=false, seed=SEED) tspan = (0.0, 10) -sweep = AdiabaticSweep(ω => (0.9, 1.1), tspan) # linearly interpolate between two values at two times -ode_problem = ODEProblem(harmonic_eq, fixed; sweep=sweep, x0=[0.01; 0.0], timespan=tspan) +sweep = ParameterSweep(ω => (0.9, 1.1), tspan) # linearly interpolate between two values at two times +ode_problem = ODEProblem(harmonic_eq, fixed; sweep=sweep, u0=[0.01; 0.0], timespan=tspan) time_soln = solve(ode_problem, Tsit5(); saveat=1); transform_solutions(time_soln, "sqrt(u1^2+v1^2)", harmonic_eq) From d9b94b338c4d4f8b80ac17ccd10aedafeccfdb4a Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 26 Oct 2024 17:55:51 +0100 Subject: [PATCH 04/15] Fix: wrong auto-merge (#291) --- ext/TimeEvolution/ODEProblem.jl | 4 ++-- test/time_evolution.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ext/TimeEvolution/ODEProblem.jl b/ext/TimeEvolution/ODEProblem.jl index ec07993a..ab2364ef 100644 --- a/ext/TimeEvolution/ODEProblem.jl +++ b/ext/TimeEvolution/ODEProblem.jl @@ -5,7 +5,7 @@ using OrdinaryDiffEqTsit5: ODEProblem, solve, ODESolution eom::HarmonicEquation; fixed_parameters, u0::Vector, - sweep::ParameterSweep, + sweep::AdiabaticSweep, timespan::Tuple ) @@ -16,7 +16,7 @@ If `u0` is specified, it is used as an initial condition; otherwise the values f function OrdinaryDiffEqTsit5.ODEProblem( eom::HarmonicEquation, fixed_parameters; - sweep::ParameterSweep=ParameterSweep(), + sweep::AdiabaticSweep=AdiabaticSweep(), u0::Vector=[], timespan::Tuple, perturb_initial=0.0, diff --git a/test/time_evolution.jl b/test/time_evolution.jl index e076bf95..1369e3bd 100644 --- a/test/time_evolution.jl +++ b/test/time_evolution.jl @@ -22,7 +22,7 @@ varied = ω => range(0.9, 1.1, 10) res = get_steady_states(p, varied, fixed; show_progress=false, seed=SEED) tspan = (0.0, 10) -sweep = ParameterSweep(ω => (0.9, 1.1), tspan) # linearly interpolate between two values at two times +sweep = AdiabaticSweep(ω => (0.9, 1.1), tspan) # linearly interpolate between two values at two times ode_problem = ODEProblem(harmonic_eq, fixed; sweep=sweep, u0=[0.01; 0.0], timespan=tspan) time_soln = solve(ode_problem, Tsit5(); saveat=1); From b2068dbffa9302b34bb9ee8379a5d747c64292c8 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 26 Oct 2024 18:27:00 +0100 Subject: [PATCH 05/15] test: enable JET (#284) --- src/HarmonicEquation.jl | 2 +- src/HarmonicVariable.jl | 4 +++- src/KrylovBogoliubov/KrylovEquation.jl | 2 +- src/LinearResponse/LinearResponse.jl | 2 +- src/LinearResponse/Lorentzian_spectrum.jl | 3 ++- test/runtests.jl | 8 ++++---- 6 files changed, 12 insertions(+), 9 deletions(-) diff --git a/src/HarmonicEquation.jl b/src/HarmonicEquation.jl index c67940cf..7822e915 100644 --- a/src/HarmonicEquation.jl +++ b/src/HarmonicEquation.jl @@ -167,7 +167,7 @@ end "Simplify the equations in HarmonicEquation." function simplify!(eom::HarmonicEquation) - return eom.equations = [simplify(eq) for eq in eom.equations] + return eom.equations = [Symbolics.simplify(eq) for eq in eom.equations] end """ diff --git a/src/HarmonicVariable.jl b/src/HarmonicVariable.jl index 7a1702c5..149c8c8e 100644 --- a/src/HarmonicVariable.jl +++ b/src/HarmonicVariable.jl @@ -25,7 +25,9 @@ end # when HV is used for substitute, substitute its symbol function ExprUtils.substitute_all(eq::Union{Num,Equation}, rules::Dict{HarmonicVariable}) - return substitute(eq, Dict(zip(getfield.(keys(rules), :symbol), values(rules)))) + return Symbolics.substitute( + eq, Dict(zip(getfield.(keys(rules), :symbol), values(rules))) + ) end function ExprUtils.substitute_all(var::HarmonicVariable, rules) diff --git a/src/KrylovBogoliubov/KrylovEquation.jl b/src/KrylovBogoliubov/KrylovEquation.jl index 439b9231..eb308be0 100644 --- a/src/KrylovBogoliubov/KrylovEquation.jl +++ b/src/KrylovBogoliubov/KrylovEquation.jl @@ -116,7 +116,7 @@ function van_der_Pol(eom::DifferentialEquation, t::Num) D = Differential(t) nvar_t = diff2term(D(unwrap(nvar))) vdP_rules = Dict(D(hvar_u.symbol) => 0, D(hvar_v.symbol) => 0) - rules[nvar_t] = substitute(expand_derivatives(D(rule)), vdP_rules) + rules[nvar_t] = Symbolics.substitute(expand_derivatives(D(rule)), vdP_rules) uv_idx += 1 push!(vars, hvar_u, hvar_v) diff --git a/src/LinearResponse/LinearResponse.jl b/src/LinearResponse/LinearResponse.jl index 23fb42be..c3390fbe 100644 --- a/src/LinearResponse/LinearResponse.jl +++ b/src/LinearResponse/LinearResponse.jl @@ -8,7 +8,7 @@ using Plots: heatmap, theme_palette, scatter, RGB, cgrad using Latexify: Latexify, latexify, @L_str using Latexify.LaTeXStrings: LaTeXStrings -using Symbolics: Num, Equation, substitute, unwrap +using Symbolics: Symbolics, Num, Equation, unwrap using LinearAlgebra: norm, eigvals, eigen, eigvecs using OrderedCollections: OrderedDict diff --git a/src/LinearResponse/Lorentzian_spectrum.jl b/src/LinearResponse/Lorentzian_spectrum.jl index d83e7ee1..44c1ebce 100644 --- a/src/LinearResponse/Lorentzian_spectrum.jl +++ b/src/LinearResponse/Lorentzian_spectrum.jl @@ -90,7 +90,8 @@ function JacobianSpectrum(res::Result; index::Int, branch::Int, force=false) for pair in _get_uv_pairs(hvars) u, v = hvars[pair] eigvec_2d = eigvec[pair] # fetch the relevant part of the Jacobian eigenvector - ωnum = real(unwrap(substitute(u.ω, solution_dict))) # the harmonic (numerical now) associated to this harmonic variable + ωnum = real(unwrap(Symbolics.substitute(u.ω, solution_dict))) + # ^ the harmonic (numerical now) associated to this harmonic variable # eigvec_2d is associated to a natural variable -> this variable gets Lorentzian peaks peaks = norm(eigvec_2d) * _pair_to_peaks(λ, eigvec_2d; ω=ωnum) diff --git a/test/runtests.jl b/test/runtests.jl index e4883fae..1081be93 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,10 +5,10 @@ using Random const SEED = 0xd8e5d8df Random.seed!(SEED) -# @testset "Code linting" begin -# using JET -# JET.test_package(HarmonicBalance; target_defined_modules=true) -# end +@testset "Code linting" begin + using JET + JET.test_package(HarmonicBalance; target_defined_modules=true) +end @testset "Code quality" begin using ExplicitImports, Aqua From dbc22df76aa7d9e294bcafeeaecec2080b8a3b21 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 27 Oct 2024 10:30:20 +0000 Subject: [PATCH 06/15] fix: fix "mathtt{}" symbolic bug (#292) --- Project.toml | 4 +-- src/HarmonicVariable.jl | 7 ++-- test/HarmonicVariable.jl | 71 ++++------------------------------------ 3 files changed, 12 insertions(+), 70 deletions(-) diff --git a/Project.toml b/Project.toml index da1a3ea5..59a8a3da 100644 --- a/Project.toml +++ b/Project.toml @@ -58,8 +58,8 @@ Plots = "1.39" PrecompileTools = "1.2" ProgressMeter = "1.7.2" SteadyStateDiffEq = "2.3.2" -SymbolicUtils = "3.7" -Symbolics = "~6.14" +SymbolicUtils = "3.5" +Symbolics = "6.4" julia = "1.10.0" [extras] diff --git a/src/HarmonicVariable.jl b/src/HarmonicVariable.jl index 149c8c8e..332f987a 100644 --- a/src/HarmonicVariable.jl +++ b/src/HarmonicVariable.jl @@ -77,9 +77,10 @@ function declare_variable(name::String, independent_variable::Num) end "Return the name of a variable (excluding independent variables)" -function var_name(x::Num) +function var_name(x::Num)::String var = Symbolics._toexpr(x) - return var isa Expr ? String(var.args[1]) : String(var) + var = var isa Expr ? String(var.args[1]) : String(var) + return String(replace(var, r"\\mathtt\{([^}]*)\}" => s"\1")) + # ^ remove "\\mathtt{}" from the variable name coming from Symbolics since Symbolics v6.14.1 (Symbolics#1305) end -# var_name(x::Term) = String(Symbolics._toexpr(x).args[1]) var_name(x::SymbolicUtils.Sym) = String(x.name) diff --git a/test/HarmonicVariable.jl b/test/HarmonicVariable.jl index 58937dec..88797607 100644 --- a/test/HarmonicVariable.jl +++ b/test/HarmonicVariable.jl @@ -5,69 +5,10 @@ using HarmonicBalance diff_eq = DifferentialEquation( d(x, t, 2) + ω0 * x + α * x^3 + γ * d(x, t) + η * x^2 * d(x, t) ~ F * cos(ω * t), x ) # define ODE +add_harmonic!(diff_eq, x, ω) +harmonic_eq = get_harmonic_equations(diff_eq) -add_harmonic!(diff_eq, x, ω) # specify the ansatz x = u(T) cos(ωt) + v(T) sin(ωt) - -harmonic_eq = get_harmonic_equations(diff_eq) # implement ansatz to get harmonic equations -get_variables(harmonic_eq) - -# @testset "HarmonicVariable Tests" begin -# # Test display function -# @testset "Display" begin -# var = HarmonicVariable(:x, "x", "u", 1.0, 2.0) -# @test display(var) == :x -# vars = [ -# HarmonicVariable(:x, "x", "u", 1.0, 2.0), -# HarmonicVariable(:y, "y", "v", 2.0, 3.0), -# ] -# @test display(vars) == [:x, :y] -# end - -# # Test _coordinate_transform function -# @testset "_coordinate_transform" begin -# @test _coordinate_transform(2.0, 1.0, 0.0, "u") == 2.0 -# @test _coordinate_transform(2.0, 1.0, 0.0, "v") == 0.0 -# @test _coordinate_transform(2.0, 1.0, 0.0, "a") == 2.0 -# end - -# # Test _create_harmonic_variable function -# @testset "_create_harmonic_variable" begin -# rule, hvar = _create_harmonic_variable(2.0, 1.0, 0.0, "u"; new_symbol="x") -# @test rule == 2.0 -# @test hvar == HarmonicVariable(:x, "u_{2.0}", "u", 1.0, 2.0) -# end - -# # Test substitute_all function -# @testset "substitute_all" begin -# eq = 2.0 * :x + 3.0 * :y -# rules = Dict(:x => 1.0, :y => 2.0) -# @test substitute_all(eq, rules) == 2.0 + 6.0 -# var = HarmonicVariable(:x, "x", "u", 1.0, 2.0) -# @test substitute_all(var, rules) == HarmonicVariable(1.0, "x", "u", 1.0, 2.0) -# vars = [ -# HarmonicVariable(:x, "x", "u", 1.0, 2.0), -# HarmonicVariable(:y, "y", "v", 2.0, 3.0), -# ] -# @test substitute_all(vars, rules) == [ -# HarmonicVariable(1.0, "x", "u", 1.0, 2.0), -# HarmonicVariable(2.0, "y", "v", 2.0, 3.0), -# ] -# end - -# # Test Symbolics.get_variables function -# @testset "Symbolics.get_variables" begin -# vars = [1.0, 2.0, 3.0] -# @test Symbolics.get_variables(vars) == [1.0, 2.0, 3.0] -# var = HarmonicVariable(:x, "x", "u", 1.0, 2.0) -# @test Symbolics.get_variables(var) == [:x] -# end - -# # Test isequal function -# @testset "isequal" begin -# var1 = HarmonicVariable(:x, "x", "u", 1.0, 2.0) -# var2 = HarmonicVariable(:x, "x", "u", 1.0, 2.0) -# var3 = HarmonicVariable(:y, "y", "v", 2.0, 3.0) -# @test isequal(var1, var2) == true -# @test isequal(var1, var3) == false -# end -# end +@testset "" begin + vars = get_variables(harmonic_eq) + @test HarmonicBalance.var_name.(vars) == ["u1", "v1"] +end From 19292cefceeddc367df55359d85a98cc2c1c8f22 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 27 Oct 2024 11:55:19 +0000 Subject: [PATCH 07/15] fix: make 2nd order response work and add some tests (#293) --- src/ExprUtils/Symbolics_utils.jl | 26 +++---------- src/LinearResponse/response.jl | 14 ++++--- test/linear_response.jl | 64 ++++++++++++++++++++------------ 3 files changed, 55 insertions(+), 49 deletions(-) diff --git a/src/ExprUtils/Symbolics_utils.jl b/src/ExprUtils/Symbolics_utils.jl index a5560a1f..ff467066 100644 --- a/src/ExprUtils/Symbolics_utils.jl +++ b/src/ExprUtils/Symbolics_utils.jl @@ -36,8 +36,8 @@ $(TYPEDSIGNATURES) Perform substitutions in `rules` on `x`. `include_derivatives=true` also includes all derivatives of the variables of the keys of `rules`. """ -subtype = Union{Num,Equation,BasicSymbolic} -function substitute_all(x::subtype, rules::Dict; include_derivatives=true) +Subtype = Union{Num,Equation,BasicSymbolic} +function substitute_all(x::Subtype, rules::Dict; include_derivatives=true) if include_derivatives rules = merge( rules, @@ -54,24 +54,10 @@ function substitute_all(dict::Dict, rules::Dict)::Dict end Collections = Union{Dict,Pair,Vector,OrderedDict} substitute_all(v::AbstractArray, rules) = [substitute_all(x, rules) for x in v] -substitute_all(x::subtype, rules::Collections) = substitute_all(x, Dict(rules)) -# Collections = Union{Dict,OrderedDict} -# function substitute_all(x, rules::Collections; include_derivatives=true) -# if include_derivatives -# rules = merge( -# rules, -# Dict([Differential(var) => Differential(rules[var]) for var in keys(rules)]), -# ) -# end -# return substitute(x, rules) -# end -# "Variable substitution - dictionary" -# function substitute_all(dict::Dict, rules::Dict)::Dict -# new_keys = substitute_all.(keys(dict), rules) -# new_values = substitute_all.(values(dict), rules) -# return Dict(zip(new_keys, new_values)) -# end -# substitute_all(v::AbstractArray, rules::Collections) = [substitute_all(x, rules) for x in v] +substitute_all(x::Subtype, rules::Collections) = substitute_all(x, Dict(rules)) +function substitute_all(x::Complex{Num}, rules::Collections) + return substitute_all(x.re, rules) + im * substitute_all(x.im, rules) +end get_independent(x::Num, t::Num) = get_independent(x.val, t) function get_independent(x::Complex{Num}, t::Num) diff --git a/src/LinearResponse/response.jl b/src/LinearResponse/response.jl index ac854834..312dcb68 100644 --- a/src/LinearResponse/response.jl +++ b/src/LinearResponse/response.jl @@ -20,12 +20,12 @@ function get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2): # get the response matrix by summing the orders M = get_Jacobian(eom.equations, get_variables(eom)) for n in 1:order - M += (i * freq)^n * get_Jacobian(eom.equations, d(get_variables(eom), T, n)) + M += (im * freq)^n * get_Jacobian(eom.equations, d(get_variables(eom), T, n)) end M = substitute_all( M, [var => declare_variable(var_name(var)) for var in get_variables(eom)] ) - return substitute_all(expand_derivatives.(M), i => im) + return M end "Get the response matrix corresponding to `res`. @@ -37,9 +37,13 @@ function ResponseMatrix(res::Result; rules=Dict()) M = get_response_matrix(res.problem.eom.natural_equation, Δ; order=2) M = substitute_all(M, merge(res.fixed_parameters, rules)) symbols = _free_symbols(res) - compiled_M = [Symbolics.build_function(el, cat(symbols, [Δ]; dims=1)) for el in M] - - return ResponseMatrix(eval.(compiled_M), symbols, res.problem.eom.variables) + compiled_M = map(M) do el + args = cat(symbols, [Δ]; dims=1) + f_re = Symbolics.build_function(el.re, args; expression=Val{false}) + f_im = Symbolics.build_function(el.im, args; expression=Val{false}) + (args...) -> f_re(args...) + im * f_im(args...) + end + return ResponseMatrix(compiled_M, symbols, res.problem.eom.variables) end """Evaluate the response matrix `resp` for the steady state `s` at (lab-frame) frequency `Ω`.""" diff --git a/test/linear_response.jl b/test/linear_response.jl index c0d86f55..a862cb91 100644 --- a/test/linear_response.jl +++ b/test/linear_response.jl @@ -1,4 +1,5 @@ -using HarmonicBalance +using HarmonicBalance, Symbolics +HB = HarmonicBalance @variables α, ω, ω0, F, γ, t, x(t); @@ -13,31 +14,46 @@ varied = ω => range(0.9, 1.1, 10) result = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) -plot_linear_response( - result, x; branch=1, Ω_range=range(0.9, 1.1, 10), order=1, logscale=true -) - -plot_rotframe_jacobian_response( - result; Ω_range=range(0.01, 1.1, 10), branch=1, logscale=true -) +@testset "first order" begin + plot_linear_response( + result, x; branch=1, Ω_range=range(0.9, 1.1, 10), order=1, logscale=true + ) + plot_rotframe_jacobian_response( + result; Ω_range=range(0.01, 1.1, 10), branch=1, logscale=true + ) +end -plot_eigenvalues(result; branch=1) -plot_eigenvalues(result; branch=1, type=:re, class="all") +@testset "second order" begin + response_matrix = HB.LinearResponse.ResponseMatrix(result) + M = response_matrix.matrix + @test M[1](ones(4)) isa ComplexF64 -@testset begin - @variables α λ ω0 ω ωₚ F t x(t) - diff_eq = DifferentialEquation( - d(x, t, 2) + (ω0^2 - λ * cos(2 * ω * t)) * x + α * x^3 + γ * d(x, t) ~ - F * cos(ωₚ * t), - x, + plot_linear_response( + result, x; branch=1, Ω_range=range(0.9, 1.1, 10), order=2, logscale=true ) +end - add_harmonic!(diff_eq, x, ω) - add_harmonic!(diff_eq, x, ωₚ) - harmonic_eq = get_harmonic_equations(diff_eq) - - fixed = (γ => 0.008, ω0 => 1.0, α => 1.0, F => 0.0, ωₚ => 1.0, λ => 0.016) - varied = (ω => range(0.995, 1.005, 20)) - result_ω = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) - @test_throws ErrorException plot_eigenvalues(result_ω; branch=1) +@testset "eigenvalues" begin + plot_eigenvalues(result; branch=1) + plot_eigenvalues(result; branch=1, type=:re, class="all") + + @testset "NaN Exception Error" begin + @variables α λ ω0 ω ωₚ F t x(t) + diff_eq = DifferentialEquation( + d(x, t, 2) + (ω0^2 - λ * cos(2 * ω * t)) * x + α * x^3 + γ * d(x, t) ~ + F * cos(ωₚ * t), + x, + ) + + add_harmonic!(diff_eq, x, ω) + add_harmonic!(diff_eq, x, ωₚ) + harmonic_eq = get_harmonic_equations(diff_eq) + + fixed = (γ => 0.008, ω0 => 1.0, α => 1.0, F => 0.0, ωₚ => 1.0, λ => 0.016) + varied = (ω => range(0.995, 1.005, 20)) + result_ω = get_steady_states( + harmonic_eq, varied, fixed; show_progress=false, seed=SEED + ) + @test_throws ErrorException plot_eigenvalues(result_ω; branch=1) + end end From 8715ee54c1fbf2a6117327f1908f966d0fd7d1fd Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 27 Oct 2024 12:21:23 +0000 Subject: [PATCH 08/15] style: move some of the functions (#294) --- src/classification.jl | 8 ++++ src/plotting_Plots.jl | 39 ---------------- src/solve_homotopy.jl | 92 -------------------------------------- src/sorting.jl | 27 +++++++++++ src/transform_solutions.jl | 74 +++++++++++++++++++++++++++++- src/types.jl | 10 +++++ src/utils.jl | 16 +++++++ 7 files changed, 133 insertions(+), 133 deletions(-) diff --git a/src/classification.jl b/src/classification.jl index 9085fd7c..dde2d413 100644 --- a/src/classification.jl +++ b/src/classification.jl @@ -148,3 +148,11 @@ function filter_result!(res::Result, class::String) res.classes[c] = [s[bools] for s in res.classes[c]] end end + +function _classify_default!(result) + classify_solutions!(result, _is_physical, "physical") + classify_solutions!(result, _is_stable(result), "stable") + classify_solutions!(result, _is_Hopf_unstable(result), "Hopf") + order_branches!(result, ["physical", "stable"]) # shuffle the branches to have relevant ones first + return classify_binaries!(result) # assign binaries to solutions depending on which branches are stable +end diff --git a/src/plotting_Plots.jl b/src/plotting_Plots.jl index 572b6f1c..3f5bc161 100644 --- a/src/plotting_Plots.jl +++ b/src/plotting_Plots.jl @@ -67,42 +67,6 @@ Similar to `plot` but adds a plot onto an existing plot. function Plots.plot!(res::Result, varargs...; kwargs...)::Plots.Plot return plot(res, varargs...; add=true, _set_Plots_default..., kwargs...) end -""" -$(TYPEDSIGNATURES) - -Return an array of bools to mark solutions in `res` which fall into `classes` but not `not_classes`. -Only `branches` are considered. -""" -function _get_mask(res, classes, not_classes=[]; branches=1:branch_count(res)) - classes == "all" && return fill(trues(length(branches)), size(res.solutions)) - bools = vcat( - [res.classes[c] for c in _str_to_vec(classes)], - [map(.!, res.classes[c]) for c in _str_to_vec(not_classes)], - ) - #m = map( x -> [getindex(x, b) for b in [branches...]], map(.*, bools...)) - - return m = map(x -> x[[branches...]], map(.*, bools...)) -end - -""" -$(TYPEDSIGNATURES) - -Go over a solution and an equally-sized array (a "mask") of booleans. -true -> solution unchanged -false -> changed to NaN (omitted from plotting) -""" -function _apply_mask(solns::Array{Vector{ComplexF64}}, booleans) - factors = replace.(booleans, 0 => NaN) - return map(.*, solns, factors) -end -function _apply_mask(solns::Vector{Vector{Vector{ComplexF64}}}, booleans) - Nan_vector = NaN .* similar(solns[1][1]) - new_solns = [ - [booleans[i][j] ? solns[i][j] : Nan_vector for j in eachindex(solns[i])] for - i in eachindex(solns) - ] - return new_solns -end """ Project the array `a` into the real axis, warning if its contents are complex. """ function _realify(a::Array{T} where {T<:Number}; warning="") @@ -119,9 +83,6 @@ function _realify(a::Array{T} where {T<:Number}; warning="") return a_real end -_str_to_vec(s::Vector) = s -_str_to_vec(s) = [s] - # return true if p already has a label for branch index idx function _is_labeled(p::Plot, idx::Int64) return in(string(idx), [sub[:label] for sub in p.series_list]) diff --git a/src/solve_homotopy.jl b/src/solve_homotopy.jl index 762f495e..2bed7e03 100644 --- a/src/solve_homotopy.jl +++ b/src/solve_homotopy.jl @@ -1,49 +1,3 @@ -# assume this order of variables in all compiled function (transform_solutions, Jacobians) -function _free_symbols(res::Result) - return cat(res.problem.variables, collect(keys(res.swept_parameters)); dims=1) -end -function _free_symbols(p::Problem, varied) - return cat(p.variables, collect(keys(OrderedDict(varied))); dims=1) -end -_symidx(sym::Num, args...) = findfirst(x -> isequal(x, sym), _free_symbols(args...)) - -""" -$(TYPEDSIGNATURES) -Return an ordered dictionary specifying all variables and parameters of the solution -in `result` on `branch` at the position `index`. -""" -function get_single_solution(res::Result; branch::Int64, index)::OrderedDict{Num,ComplexF64} - - # check if the dimensionality of index matches the solutions - if length(size(res.solutions)) !== length(index) - # if index is a number, use linear indexing - index = if length(index) == 1 - CartesianIndices(res.solutions)[index] - else - error("Index ", index, " undefined for a solution of size ", size(res.solutions)) - end - else - index = CartesianIndex(index) - end - - vars = OrderedDict(zip(res.problem.variables, res.solutions[index][branch])) - - # collect the swept parameters required for this call - swept_params = OrderedDict( - key => res.swept_parameters[key][index[i]] for - (i, key) in enumerate(keys(res.swept_parameters)) - ) - full_solution = merge(vars, swept_params, res.fixed_parameters) - - return OrderedDict(zip(keys(full_solution), ComplexF64.(values(full_solution)))) -end - -function get_single_solution(res::Result, index) - return [ - get_single_solution(res; index=index, branch=b) for b in 1:length(res.solutions[1]) - ] -end - """ get_steady_states(prob::Problem, swept_parameters::ParameterRange, @@ -162,14 +116,6 @@ function get_steady_states( return result end -function _classify_default!(result) - classify_solutions!(result, _is_physical, "physical") - classify_solutions!(result, _is_stable(result), "stable") - classify_solutions!(result, _is_Hopf_unstable(result), "Hopf") - order_branches!(result, ["physical", "stable"]) # shuffle the branches to have relevant ones first - return classify_binaries!(result) # assign binaries to solutions depending on which branches are stable -end - function get_steady_states(p::Problem, swept, fixed; kwargs...) return get_steady_states(p, ParameterRange(swept), ParameterList(fixed); kwargs...) end @@ -218,33 +164,6 @@ function compile_matrix(mat, variables; rules=Dict(), postproc=x -> x) return m end -"Find a branch order according `classification`. Place branches where true occurs earlier first." -function find_branch_order(classification::Vector{BitVector}) - branches = [getindex.(classification, k) for k in 1:length(classification[1])] # array of branches - indices = replace(findfirst.(branches), nothing => Inf) - negative = findall(x -> x == Inf, indices) # branches not true anywhere - leave out - return order = setdiff(sortperm(indices), negative) -end - -find_branch_order(classification::Array) = collect(1:length(classification[1])) # no ordering for >1D - -"Order the solution branches in `res` such that close classified positively by `classes` are first." -function order_branches!(res::Result, classes::Vector{String}) - for class in classes - order_branches!(res, find_branch_order(res.classes[class])) - end -end - -order_branches!(res::Result, class::String) = order_branches!(res, [class]) - -"Reorder the solutions in `res` to match the index permutation `order`." -function order_branches!(res::Result, order::Vector{Int64}) - res.solutions = _reorder_nested(res.solutions, order) - for key in keys(res.classes) - res.classes[key] = _reorder_nested(res.classes[key], order) - end -end - "Reorder EACH ELEMENT of `a` to match the index permutation `order`. If length(order) < length(array), the remanining positions are kept." function _reorder_nested(a::Array, order::Vector{Int64}) a[1] isa Union{Array,BitVector} || return a @@ -387,8 +306,6 @@ function pad_solutions(solutions::Array{Vector{Vector{ComplexF64}}}; padding_val return padded_solutions end -tuple_to_vector(t::Tuple) = [i for i in t] - function newton(prob::Problem, soln::OrderedDict) vars = _convert_or_zero.(substitute_all(prob.variables, soln)) pars = _convert_or_zero.(substitute_all(prob.parameters, soln)) @@ -405,12 +322,3 @@ Any variables/parameters not present in `soln` are set to zero. """ newton(res::Result, soln::OrderedDict) = newton(res.problem, soln) newton(res::Result; branch, index) = newton(res, res[index][branch]) - -function _convert_or_zero(x, t=ComplexF64) - try - convert(t, x) - catch ArgumentError - @warn string(x) * " not supplied: setting to zero" - return 0 - end -end diff --git a/src/sorting.jl b/src/sorting.jl index 10694a58..f3ddeb30 100644 --- a/src/sorting.jl +++ b/src/sorting.jl @@ -194,3 +194,30 @@ function sort_2D( end return sorted_solns end + +"Find a branch order according `classification`. Place branches where true occurs earlier first." +function find_branch_order(classification::Vector{BitVector}) + branches = [getindex.(classification, k) for k in 1:length(classification[1])] # array of branches + indices = replace(findfirst.(branches), nothing => Inf) + negative = findall(x -> x == Inf, indices) # branches not true anywhere - leave out + return order = setdiff(sortperm(indices), negative) +end + +find_branch_order(classification::Array) = collect(1:length(classification[1])) # no ordering for >1D + +"Order the solution branches in `res` such that close classified positively by `classes` are first." +function order_branches!(res::Result, classes::Vector{String}) + for class in classes + order_branches!(res, find_branch_order(res.classes[class])) + end +end + +order_branches!(res::Result, class::String) = order_branches!(res, [class]) + +"Reorder the solutions in `res` to match the index permutation `order`." +function order_branches!(res::Result, order::Vector{Int64}) + res.solutions = _reorder_nested(res.solutions, order) + for key in keys(res.classes) + res.classes[key] = _reorder_nested(res.classes[key], order) + end +end diff --git a/src/transform_solutions.jl b/src/transform_solutions.jl index 2fa4660d..11ffae6a 100644 --- a/src/transform_solutions.jl +++ b/src/transform_solutions.jl @@ -1,4 +1,39 @@ -_parse_expression(exp) = exp isa String ? Num(eval(Meta.parse(exp))) : exp +""" +$(TYPEDSIGNATURES) +Return an ordered dictionary specifying all variables and parameters of the solution +in `result` on `branch` at the position `index`. +""" +function get_single_solution(res::Result; branch::Int64, index)::OrderedDict{Num,ComplexF64} + + # check if the dimensionality of index matches the solutions + if length(size(res.solutions)) !== length(index) + # if index is a number, use linear indexing + index = if length(index) == 1 + CartesianIndices(res.solutions)[index] + else + error("Index ", index, " undefined for a solution of size ", size(res.solutions)) + end + else + index = CartesianIndex(index) + end + + vars = OrderedDict(zip(res.problem.variables, res.solutions[index][branch])) + + # collect the swept parameters required for this call + swept_params = OrderedDict( + key => res.swept_parameters[key][index[i]] for + (i, key) in enumerate(keys(res.swept_parameters)) + ) + full_solution = merge(vars, swept_params, res.fixed_parameters) + + return OrderedDict(zip(keys(full_solution), ComplexF64.(values(full_solution)))) +end + +function get_single_solution(res::Result, index) + return [ + get_single_solution(res; index=index, branch=b) for b in 1:length(res.solutions[1]) + ] +end """ $(TYPEDSIGNATURES) @@ -90,7 +125,42 @@ function _similar(type, res::Result; branches=1:branch_count(res)) return [type(undef, length(branches)) for k in res.solutions] end -## TODO move masks here +""" +$(TYPEDSIGNATURES) + +Return an array of bools to mark solutions in `res` which fall into `classes` but not `not_classes`. +Only `branches` are considered. +""" +function _get_mask(res, classes, not_classes=[]; branches=1:branch_count(res)) + classes == "all" && return fill(trues(length(branches)), size(res.solutions)) + bools = vcat( + [res.classes[c] for c in _str_to_vec(classes)], + [map(.!, res.classes[c]) for c in _str_to_vec(not_classes)], + ) + #m = map( x -> [getindex(x, b) for b in [branches...]], map(.*, bools...)) + + return m = map(x -> x[[branches...]], map(.*, bools...)) +end + +""" +$(TYPEDSIGNATURES) + +Go over a solution and an equally-sized array (a "mask") of booleans. +true -> solution unchanged +false -> changed to NaN (omitted from plotting) +""" +function _apply_mask(solns::Array{Vector{ComplexF64}}, booleans) + factors = replace.(booleans, 0 => NaN) + return map(.*, solns, factors) +end +function _apply_mask(solns::Vector{Vector{Vector{ComplexF64}}}, booleans) + Nan_vector = NaN .* similar(solns[1][1]) + new_solns = [ + [booleans[i][j] ? solns[i][j] : Nan_vector for j in eachindex(solns[i])] for + i in eachindex(solns) + ] + return new_solns +end ### # TRANSFORMATIONS TO THE LAB frame diff --git a/src/types.jl b/src/types.jl index 04276447..3defa00f 100644 --- a/src/types.jl +++ b/src/types.jl @@ -220,6 +220,11 @@ function Base.show(io::IO, p::Problem) return println(io, "Symbolic Jacobian: ", !(p.jacobian == false)) end +# assume this order of variables in all compiled function (transform_solutions, Jacobians) +function _free_symbols(p::Problem, varied) + return cat(p.variables, collect(keys(OrderedDict(varied))); dims=1) +end + """ $(TYPEDEF) @@ -263,6 +268,11 @@ function Base.show(io::IO, r::Result) return println(io, "\nClasses: ", join(keys(r.classes), ", ")) end +# assume this order of variables in all compiled function (transform_solutions, Jacobians) +function _free_symbols(res::Result) + return cat(res.problem.variables, collect(keys(res.swept_parameters)); dims=1) +end + # overload to use [] for indexing Base.getindex(r::Result, idx::Int...) = get_single_solution(r, idx) Base.size(r::Result) = size(r.solutions) diff --git a/src/utils.jl b/src/utils.jl index 59fa7dbf..9872504a 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -2,3 +2,19 @@ is_real(x) = abs(imag(x)) / abs(real(x)) < IM_TOL::Float64 || abs(x) < 1e-70 is_real(x::Array) = is_real.(x) flatten(a) = collect(Iterators.flatten(a)) + +_parse_expression(exp) = exp isa String ? Num(eval(Meta.parse(exp))) : exp +_symidx(sym::Num, args...) = findfirst(x -> isequal(x, sym), _free_symbols(args...)) + +tuple_to_vector(t::Tuple) = [i for i in t] +_str_to_vec(s::Vector) = s +_str_to_vec(s) = [s] + +function _convert_or_zero(x, t=ComplexF64) + try + convert(t, x) + catch ArgumentError + @warn string(x) * " not supplied: setting to zero" + return 0 + end +end From af6450dc87cdfaf8285473847ff3ecb505b00303 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 27 Oct 2024 14:52:03 +0000 Subject: [PATCH 09/15] build: tag v0.11 (#295) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 59a8a3da..a2846859 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.11" +version = "0.11.0" [deps] BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54" From 889753c3c211baec85bbe7a92e038546f627f21c Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Tue, 29 Oct 2024 10:16:21 +0000 Subject: [PATCH 10/15] feat: add `ODEProblem(::DifferentialEquation) (#300) --- ext/ModelingToolkitExt.jl | 51 ++++++++++++++++++++++++---- test/ModelingToolkitExt.jl | 69 +++++++++++++++++++++++++++++++------- 2 files changed, 101 insertions(+), 19 deletions(-) diff --git a/ext/ModelingToolkitExt.jl b/ext/ModelingToolkitExt.jl index b5110701..de0cd3c7 100644 --- a/ext/ModelingToolkitExt.jl +++ b/ext/ModelingToolkitExt.jl @@ -3,8 +3,17 @@ module ModelingToolkitExt export ODESystem, ODEProblem, SteadyStateProblem, NonlinearProblem using HarmonicBalance: - HarmonicEquation, is_rearranged, rearrange_standard, get_variables, ParameterList -using Symbolics: simplify, Equation, substitute, Num, @variables, expand, unwrap, arguments + HarmonicEquation, + is_rearranged, + rearrange_standard, + get_variables, + ParameterList, + DifferentialEquation, + get_independent_variables +using HarmonicBalance.KrylovBogoliubov: + rearrange_standard!, is_rearranged_standard, first_order_transform! +using Symbolics: + simplify, Equation, substitute, Num, @variables, expand, unwrap, arguments, wrap using ModelingToolkit: ModelingToolkit, ODESystem, @@ -38,6 +47,7 @@ function ModelingToolkit.ODESystem(eom::HarmonicEquation) @assert isone(length(slow_time)) "The argument of the variables are not the same." slow_time_ivp = @eval @independent_variables $(Symbol(first(slow_time))) + # we have the replace the param made with @variables with the ones made with @parameters par_names = declare_parameter.(eom.parameters) eqs = deepcopy(eom.equations) @@ -45,14 +55,43 @@ function ModelingToolkit.ODESystem(eom::HarmonicEquation) eqs = simplify.(expand.(eqs)) eqs = substitute(eqs, Dict(zip(eom.parameters, par_names))) - # compute jacobian for performance # ∨ mtk v9 need @mtkbuild @mtkbuild sys = ODESystem(eqs, first(slow_time_ivp), vars, par_names) return sys end +function ModelingToolkit.ODESystem(diff_eq::DifferentialEquation) + if !is_rearranged_standard(diff_eq) + rearrange_standard!(diff_eq) + end + + times = get_independent_variables(diff_eq) + @assert isone(length(times)) "Only one independent variable allowed." + iv = first(@eval @independent_variables $(Symbol(first(times)))) + + first_order_transform!(diff_eq, iv) + + eqs = collect(values(diff_eq.equations)) + vars = get_variables(diff_eq) + + diff_eq_sym = collect(Iterators.flatten(get_variables.(eqs))) + param_undeclared = setdiff(setdiff(wrap.(diff_eq_sym), vars), iv) + params = declare_parameter.(param_undeclared) + + eqs = substitute(eqs, Dict(zip(param_undeclared, params))) + + @mtkbuild sys = ODESystem(eqs, first(iv), vars, params) + + return sys +end + function ModelingToolkit.ODEProblem( - eom::HarmonicEquation, u0, tspan::Tuple, p::ParameterList; in_place=true, kwargs... + eom::Union{HarmonicEquation,DifferentialEquation}, + u0, + tspan::Tuple, + p::ParameterList; + in_place=true, + kwargs..., ) sys = ODESystem(eom) param = varmap_to_vars(p, parameters(sys)) @@ -60,7 +99,7 @@ function ModelingToolkit.ODEProblem( prob = ODEProblem{false}(sys, u0, tspan, param; jac=true, kwargs...) else # in-place prob = ODEProblem{true}(sys, u0, tspan, param; jac=true, kwargs...) - end + end # compute jacobian for performance return prob end @@ -80,7 +119,7 @@ function ModelingToolkit.SteadyStateProblem( prob = SteadyStateProblem{false}(sys, u0, param; jac=true, kwargs...) else # in-place prob = SteadyStateProblem{true}(sys, u0, param; jac=true, kwargs...) - end + end # compute jacobian for performance return prob end diff --git a/test/ModelingToolkitExt.jl b/test/ModelingToolkitExt.jl index 22e05866..778253f1 100644 --- a/test/ModelingToolkitExt.jl +++ b/test/ModelingToolkitExt.jl @@ -1,21 +1,64 @@ using HarmonicBalance using ModelingToolkit -using ModelingToolkit: varmap_to_vars + +ModelingToolkitExt = Base.get_extension(HarmonicBalance, :ModelingToolkitExt) + using Test -@variables α ω ω0 F γ t x(t) -diff_eq = DifferentialEquation( - d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x -) -add_harmonic!(diff_eq, x, ω) # -harmonic_eq = get_harmonic_equations(diff_eq) +@testset "Utilities" begin + @variables α + check = ModelingToolkitExt.declare_parameter(α) + @test ModelingToolkit.PARAMETER ∈ values(check.val.metadata) +end + +@testset "DifferentialEquation" begin + @testset "ODESystem" begin + @variables α ω ω0 F γ t x(t) + diff_eq = DifferentialEquation( + d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x + ) + + fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) + param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) + sys = ODESystem(diff_eq) + + for p in string.([α, ω, ω0, F, γ]) + @test p ∈ string.(parameters(sys)) + end + end + @testset "ODEProblem" begin + @variables α ω ω0 F γ t x(t) + diff_eq = DifferentialEquation( + d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x + ) -sys = ODESystem(harmonic_eq) -fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) -param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) + fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) + param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) -for p in string.([α, ω, ω0, F, γ]) - @test p ∈ string.(parameters(sys)) + ODEProblem(diff_eq, [1.0, 0.0], (0, 100), param) + end end -ODEProblem(harmonic_eq, [1.0, 0.0], (0, 100), param) +@testset "HarmonicEquation" begin + @variables α ω ω0 F γ t x(t) + diff_eq = DifferentialEquation( + d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x + ) + + add_harmonic!(diff_eq, x, ω) # + harmonic_eq = get_harmonic_equations(diff_eq) + + fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) + param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) + @testset "ODESystem" begin + sys = ODESystem(harmonic_eq) + + for p in string.([α, ω, ω0, F, γ]) + @test p ∈ string.(parameters(sys)) + end + end + + @testset "ODEProblem" begin + ODEProblem(harmonic_eq, [1.0, 0.0], (0, 100), param) + end +end From 9ef69957fb89ba82b4167f27bad53cdb94550097 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Tue, 29 Oct 2024 10:18:33 +0000 Subject: [PATCH 11/15] build: tag new patch version (#301) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a2846859..ae458653 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.11.0" +version = "0.11.1" [deps] BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54" From 641fd6437003ebd7a2ab5ddb211ea5de4f70f15e Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Tue, 29 Oct 2024 10:50:41 +0000 Subject: [PATCH 12/15] feat: add `ODEProblem(::DIfferentialEquation) (#302) --- ext/ModelingToolkitExt.jl | 1 + test/ModelingToolkitExt.jl | 7 ++++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/ext/ModelingToolkitExt.jl b/ext/ModelingToolkitExt.jl index de0cd3c7..ccfbabaa 100644 --- a/ext/ModelingToolkitExt.jl +++ b/ext/ModelingToolkitExt.jl @@ -61,6 +61,7 @@ function ModelingToolkit.ODESystem(eom::HarmonicEquation) end function ModelingToolkit.ODESystem(diff_eq::DifferentialEquation) + diff_eq = deepcopy(diff_eq) if !is_rearranged_standard(diff_eq) rearrange_standard!(diff_eq) end diff --git a/test/ModelingToolkitExt.jl b/test/ModelingToolkitExt.jl index 778253f1..690bc7b3 100644 --- a/test/ModelingToolkitExt.jl +++ b/test/ModelingToolkitExt.jl @@ -1,11 +1,9 @@ using HarmonicBalance using ModelingToolkit - -ModelingToolkitExt = Base.get_extension(HarmonicBalance, :ModelingToolkitExt) - using Test @testset "Utilities" begin + ModelingToolkitExt = Base.get_extension(HarmonicBalance, :ModelingToolkitExt) @variables α check = ModelingToolkitExt.declare_parameter(α) @test ModelingToolkit.PARAMETER ∈ values(check.val.metadata) @@ -25,6 +23,9 @@ end for p in string.([α, ω, ω0, F, γ]) @test p ∈ string.(parameters(sys)) end + + # can run a second time without error; diff_eq unmutated + ODESystem(diff_eq) end @testset "ODEProblem" begin @variables α ω ω0 F γ t x(t) From 8fe062b34fe113e0b0ee0297b3e506c3e5d36cf4 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 3 Nov 2024 11:58:36 +0100 Subject: [PATCH 13/15] docs: change font to JuliaMono (#304) --- docs/make.jl | 2 +- docs/src/.vitepress/theme/style.css | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 38bf0a24..5ea4d394 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -52,6 +52,6 @@ if CI devbranch="master", target="build", branch="gh-pages", - push_preview=false, + push_preview=true, ) end diff --git a/docs/src/.vitepress/theme/style.css b/docs/src/.vitepress/theme/style.css index 4cbbfd37..4dffba33 100644 --- a/docs/src/.vitepress/theme/style.css +++ b/docs/src/.vitepress/theme/style.css @@ -28,7 +28,7 @@ https://github.com/vuejs/vitepress/blob/main/src/client/theme-default/styles/var Cantarell, "Fira Sans", "Droid Sans", "Helvetica Neue", sans-serif; /* Code Snippet font */ - --vp-font-family-mono: JuliaMono-Regular, monospace; + --vp-font-family-mono: "Julia Mono", Menlo, Monaco, Consolas, "Courier New", monospace; } From 9d1efd808cc8ed6a0b630d1c626791521b33c760 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 3 Nov 2024 12:57:48 +0100 Subject: [PATCH 14/15] add Invalidations CI (#309) --- .github/workflows/Invalidations.yml | 40 +++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 .github/workflows/Invalidations.yml diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml new file mode 100644 index 00000000..4616f1cb --- /dev/null +++ b/.github/workflows/Invalidations.yml @@ -0,0 +1,40 @@ +name: Invalidations + +on: + pull_request: + +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: always. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + evaluate: + # Only run on PRs to the default branch. + # In the PR trigger above branches can be specified only explicitly whereas this check should work for master, main, or any other default branch + if: github.base_ref == github.event.repository.default_branch + runs-on: ubuntu-latest + steps: + - uses: julia-actions/setup-julia@latest + with: + version: '1' + - uses: actions/checkout@v4 + - uses: julia-actions/julia-buildpkg@latest + - uses: julia-actions/julia-invalidations@v1 + id: invs_pr + + - uses: actions/checkout@v4 + with: + ref: ${{ github.event.repository.default_branch }} + - uses: julia-actions/julia-buildpkg@latest + - uses: julia-actions/julia-invalidations + id: invs_default + + - name: Report invalidation counts + run: | + echo "Invalidations on default branch: ${{ steps.invs_default.outputs.total }} (${{ steps.invs_default.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY + echo "This branch: ${{ steps.invs_pr.outputs.total }} (${{ steps.invs_pr.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY + - name: Check if the PR does increase number of invalidations + if: steps.invs_pr.outputs.total > steps.invs_default.outputs.total + run: exit 1 \ No newline at end of file From 9a610018d46ecc26f9487b903cee4852fe3f7254 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 3 Nov 2024 21:19:20 +0100 Subject: [PATCH 15/15] feat: Method structs (#298) --- Project.toml | 6 +- benchmark/parametron.jl | 12 +- docs/src/.vitepress/config.mts | 2 + .../parametric_via_three_wave_mixing.md | 13 +- docs/src/examples/parametron.md | 5 +- docs/src/examples/wave_mixing.md | 8 +- docs/src/manual/methods.md | 18 ++ docs/src/manual/plotting.md | 2 +- docs/src/manual/solving_harmonics.md | 4 +- docs/src/tutorials/classification.md | 2 +- docs/src/tutorials/limit_cycles.md | 7 +- docs/src/tutorials/time_dependent.md | 1 + examples/parametric_via_three_wave_mixing.jl | 11 +- examples/parametron.jl | 3 +- examples/wave_mixing.jl | 6 +- src/DifferentialEquation.jl | 78 +++++ src/HarmonicBalance.jl | 31 +- src/HarmonicEquation.jl | 57 +++- src/HarmonicVariable.jl | 41 +++ src/LimitCycles/LimitCycles.jl | 6 + src/LimitCycles/gauge_fixing.jl | 60 ++-- src/LinearResponse/LinearResponse.jl | 11 + src/LinearResponse/Lorentzian_spectrum.jl | 5 +- src/LinearResponse/jacobians.jl | 2 +- src/LinearResponse/response.jl | 4 +- src/Problem.jl | 70 +++++ src/Result.jl | 60 ++++ src/methods.jl | 224 ++++++++++++++ src/precompilation.jl | 3 +- src/solve_homotopy.jl | 234 +++++++-------- src/types.jl | 276 ------------------ test/API.jl | 17 +- test/ModelingToolkitExt.jl | 7 +- test/hysteresis_sweep.jl | 3 +- test/krylov.jl | 4 +- test/limit_cycle.jl | 8 +- test/linear_response.jl | 6 +- test/methods.jl | 43 +++ test/parametron.jl | 10 +- test/plotting.jl | 4 +- test/runtests.jl | 1 + test/time_evolution.jl | 3 - test/transform_solutions.jl | 2 +- 43 files changed, 854 insertions(+), 516 deletions(-) create mode 100644 docs/src/manual/methods.md create mode 100644 src/Problem.jl create mode 100644 src/Result.jl create mode 100644 src/methods.jl create mode 100644 test/methods.jl diff --git a/Project.toml b/Project.toml index ae458653..57273395 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +EndpointRanges = "340492b5-2a47-5f55-813d-aca7ddf97656" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" @@ -42,6 +43,7 @@ DelimitedFiles = "1.9" Distances = "0.10.11" DocStringExtensions = "0.9.3" Documenter = "1.4" +EndpointRanges = "0.2.2" ExplicitImports = "1.6" FFTW = "1.8" HomotopyContinuation = "2.9" @@ -51,8 +53,8 @@ Latexify = "0.16" ModelingToolkit = "9.34" NonlinearSolve = "3.14" OrderedCollections = "1.6" -OrdinaryDiffEqTsit5 = "1.1" OrdinaryDiffEqRosenbrock = "1.1" +OrdinaryDiffEqTsit5 = "1.1" Peaks = "0.5" Plots = "1.39" PrecompileTools = "1.2" @@ -69,8 +71,8 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/benchmark/parametron.jl b/benchmark/parametron.jl index 4b257a8c..d1d23ea3 100644 --- a/benchmark/parametron.jl +++ b/benchmark/parametron.jl @@ -1,4 +1,5 @@ using HarmonicBalance +using BenchmarkTools using Random const SEED = 0xd8e5d8df @@ -18,8 +19,13 @@ dEOM = DifferentialEquation(natural_equation + forces, x) add_harmonic!(dEOM, x, ω) @time harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); -@time prob = HarmonicBalance.Problem(harmonic_eq) fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 0, α => 1.0, η => 0.3, θ => 0, ψ => 0) -varied = ω => range(0.9, 1.1, 20) -@time res = get_steady_states(prob, varied, fixed; show_progress=false) +varied = ω => range(0.9, 1.1, 100) + +@btime res = get_steady_states(harmonic_eq, WarmUp(), varied, fixed; show_progress=false) +@btime res = get_steady_states(harmonic_eq, WarmUp(;compile=true), varied, fixed; show_progress=false) +@btime res = get_steady_states(harmonic_eq, Polyhedral(;only_non_zero=true), varied, fixed; show_progress=false) +@btime res = get_steady_states(harmonic_eq, Polyhedral(;only_non_zero=false), varied, fixed; show_progress=false) +@btime res = get_steady_states(harmonic_eq, TotalDegree(;compile=true), varied, fixed; show_progress=false) +@btime res = get_steady_states(harmonic_eq, TotalDegree(), varied, fixed; show_progress=false) diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts index b823cdcc..6fc8ed59 100644 --- a/docs/src/.vitepress/config.mts +++ b/docs/src/.vitepress/config.mts @@ -64,6 +64,7 @@ export default defineConfig({ text: 'Manual', items: [ { text: 'Entering equations of motion', link: '/manual/entering_eom.md' }, { text: 'Computing effective system', link: '/manual/extracting_harmonics' }, + { text: 'Computing steady states', link: '/manual/methods' }, { text: 'Krylov-Bogoliubov', link: '/manual/Krylov-Bogoliubov_method' }, { text: 'Time evolution', link: '/manual/time_dependent' }, { text: 'Linear response', link: '/manual/linear_response' }, @@ -100,6 +101,7 @@ export default defineConfig({ "/manual/": [ { text: 'Entering equations of motion', link: '/manual/entering_eom.md' }, { text: 'Computing effective system', link: '/manual/extracting_harmonics' }, + { text: 'Computing steady states', link: '/manual/methods' }, { text: 'Krylov-Bogoliubov', link: '/manual/Krylov-Bogoliubov_method' }, { text: 'Time evolution', link: '/manual/time_dependent' }, { text: 'Linear response', link: '/manual/linear_response' }, diff --git a/docs/src/examples/parametric_via_three_wave_mixing.md b/docs/src/examples/parametric_via_three_wave_mixing.md index 7908825c..3bc2b036 100644 --- a/docs/src/examples/parametric_via_three_wave_mixing.md +++ b/docs/src/examples/parametric_via_three_wave_mixing.md @@ -1,5 +1,5 @@ ```@meta -EditURL = "parametric_via_three_wave_mixing.jl" +EditURL = "../../../examples/parametric_via_three_wave_mixing.jl" ``` # Parametric Pumping via Three-Wave Mixing @@ -33,7 +33,7 @@ If we both have quadratic and cubic nonlineariy, we observe the normal duffing o varied = (ω => range(0.99, 1.1, 200)) # range of parameter values fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters -result = get_steady_states(harmonic_eq, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq, varied, fixed) plot(result; y="u1^2+v1^2") ```` @@ -43,7 +43,7 @@ If we set the cubic nonlinearity to zero, we recover the driven damped harmonic varied = (ω => range(0.99, 1.1, 100)) fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) -result = get_steady_states(harmonic_eq, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq, varied, fixed) plot(result; y="u1^2+v1^2") ```` @@ -65,7 +65,7 @@ harmonic_eq2 = get_krylov_equations(diff_eq; order=2) varied = (ω => range(0.4, 1.1, 500)) fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.001, F => 0.005) -result = get_steady_states(harmonic_eq2, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq2, varied, fixed) plot(result; y="v1") ```` @@ -73,9 +73,8 @@ plot(result; y="v1") varied = (ω => range(0.4, 0.6, 100), F => range(1e-6, 0.01, 50)) fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.01) -result = get_steady_states( - harmonic_eq2, varied, fixed; threading=true, method=:total_degree -) +method = TotalDegree() +result = get_steady_states(harmonic_eq2, method, varied, fixed) plot_phase_diagram(result; class="stable") ```` diff --git a/docs/src/examples/parametron.md b/docs/src/examples/parametron.md index 26b5c929..4db5b685 100644 --- a/docs/src/examples/parametron.md +++ b/docs/src/examples/parametron.md @@ -1,5 +1,5 @@ ```@meta -EditURL = "parametron.jl" +EditURL = "../../../examples/parametron.jl" ``` # [Parametrically driven resonator](@id parametron) @@ -64,7 +64,7 @@ varied = ω => range(0.9, 1.1, 100) result = get_steady_states(harmonic_eq, varied, fixed) ```` -In `get_steady_states`, the default value for the keyword `method=:random_warmup` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `:total_degree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower. The `threading` keyword enables parallel tracking of homotopy paths, and it's set to `false` simply because we are using a single core computer for now. +In `get_steady_states`, the default method `WarmUp()` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `TotalDegree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower. After solving the system, we can save the full output of the simulation and the model (e.g. symbolic expressions for the harmonic equations) into a file @@ -100,6 +100,7 @@ The parametrically driven oscillator boasts a stability diagram called "Arnold's To perform a 2D sweep over driving frequency $\omega$ and parametric drive strength $\lambda$, we keep `fixed` from before but include 2 variables in `varied` ````@example parametron +fixed = (ω₀ => 1.0, γ => 1e-2, F => 1e-3, α => 1.0, η => 0.3) varied = (ω => range(0.8, 1.2, 50), λ => range(0.001, 0.6, 50)) result_2D = get_steady_states(harmonic_eq, varied, fixed); nothing #hide diff --git a/docs/src/examples/wave_mixing.md b/docs/src/examples/wave_mixing.md index 95df2bad..2be86998 100644 --- a/docs/src/examples/wave_mixing.md +++ b/docs/src/examples/wave_mixing.md @@ -1,5 +1,5 @@ ```@meta -EditURL = "wave_mixing.jl" +EditURL = "../../../examples/wave_mixing.jl" ``` # Three Wave Mixing vs four wave mixing @@ -39,7 +39,7 @@ response with no response at $2\omega$. ````@example wave_mixing varied = (ω => range(0.9, 1.2, 200)) # range of parameter values fixed = (α => 1.0, β => 0.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters -result = get_steady_states(harmonic_eq, varied, fixed; threading=true)# compute steady states +result = get_steady_states(harmonic_eq, varied, fixed)# compute steady states p1 = plot(result; y="√(u1^2+v1^2)", legend=:best) p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1)) @@ -57,7 +57,7 @@ We would like to investigate the three-wave mixing of the driven Duffing oscilla ````@example wave_mixing varied = (ω => range(0.9, 1.2, 200)) fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) -result = get_steady_states(harmonic_eq, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq, varied, fixed) p1 = plot(result; y="√(u1^2+v1^2)", legend=:best) p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1)) @@ -75,7 +75,7 @@ We would like to investigate the three-wave mixing of the driven Duffing oscilla ````@example wave_mixing varied = (ω => range(0.9, 1.2, 200)) fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) -result = get_steady_states(harmonic_eq, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq, varied, fixed) p1 = plot(result; y="√(u1^2+v1^2)", legend=:best) p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1)) diff --git a/docs/src/manual/methods.md b/docs/src/manual/methods.md new file mode 100644 index 00000000..80e8b52a --- /dev/null +++ b/docs/src/manual/methods.md @@ -0,0 +1,18 @@ +# Methods + +We offer several methods for solving the nonlinear algebraic equations that arise from the harmonic balance procedure. Each method has different tradeoffs between speed, robustness, and completeness. + +## Total Degree Method +```@docs +TotalDegree +``` + +## Polyhedral Method +```@docs +Polyhedral +``` + +## Warm Up Method +```@docs +WarmUp +``` \ No newline at end of file diff --git a/docs/src/manual/plotting.md b/docs/src/manual/plotting.md index bb096aab..6db443aa 100644 --- a/docs/src/manual/plotting.md +++ b/docs/src/manual/plotting.md @@ -12,7 +12,7 @@ The function `plot` is multiple-dispatched to plot 1D and 2D datasets. In 1D, the solutions are colour-coded according to the branches obtained by `sort_solutions`. ```@docs -HarmonicBalance.plot(::Result, varags...) +HarmonicBalance.plot(::HarmonicBalance.Result, varags...) ``` diff --git a/docs/src/manual/solving_harmonics.md b/docs/src/manual/solving_harmonics.md index 607d878a..65057e61 100644 --- a/docs/src/manual/solving_harmonics.md +++ b/docs/src/manual/solving_harmonics.md @@ -8,9 +8,9 @@ having called `get_harmonic_equations`, we need to set all time-derivatives to z Once defined, a `Problem` can be solved for a set of input parameters using `get_steady_states` to obtain `Result`. ```@docs -Problem +HarmonicBalance.Problem get_steady_states -Result +HarmonicBalance.Result ``` diff --git a/docs/src/tutorials/classification.md b/docs/src/tutorials/classification.md index 8996b79a..b79178ec 100644 --- a/docs/src/tutorials/classification.md +++ b/docs/src/tutorials/classification.md @@ -20,7 +20,7 @@ We performe a 2d sweep in the driving frequency $\omega$ and driving strength $\ fixed = (ω₀ => 1.0, γ => 0.002, α => 1.0) varied = (ω => range(0.99, 1.01, 100), λ => range(1e-6, 0.03, 100)) -result_2D = get_steady_states(harmonic_eq, varied, fixed, threading=true) +result_2D = get_steady_states(harmonic_eq, varied, fixed) ``` By default the steady states of the system are classified by four different catogaries: * `physical`: Solutions that are physical, i.e., all variables are purely real. diff --git a/docs/src/tutorials/limit_cycles.md b/docs/src/tutorials/limit_cycles.md index bdeaab84..acd790a9 100644 --- a/docs/src/tutorials/limit_cycles.md +++ b/docs/src/tutorials/limit_cycles.md @@ -62,8 +62,8 @@ Here we reconstruct the results of [Zambon et al., Phys Rev. A 102, 023526 (2020 ```@example lc using HarmonicBalance @variables γ F α ω0 F0 η ω J t x(t) y(t); -eqs = [d(x,t,2) + γ*d(x,t) + ω0^2*x + α*x^3+ 2*J*ω0*(x-y) - F0*cos(ω*t), - d(y,t,2) + γ * d(y,t) + ω0^2 * y + α*y^3 + 2*J*ω0*(y-x) - η*F0*cos(ω*t)] +eqs = [d(x,t,2) + γ*d(x,t) + ω0^2*x + α*x^3 + 2*J*ω0*(x-y) - F0*cos(ω*t), + d(y,t,2) + γ*d(y,t) + ω0^2*y + α*y^3 + 2*J*ω0*(y-x) - η*F0*cos(ω*t)] diff_eq = DifferentialEquation(eqs, [x,y]) ``` @@ -82,8 +82,7 @@ fixed = ( J => 154.1e-6, # coupling term α => 3.867e-7, # Kerr nonlinearity ω => 1.4507941, # pump frequency, resonant with antisymmetric mode (in paper, ħω0 + J) - η => -0.08, # pumping leaking to site 2 (F2 = ηF1) - F0 => 0.002 # pump amplitude (overriden in sweeps) + η => -0.08 # pumping leaking to site 2 (F2 = ηF1) ) varied = F0 => range(0.002, 0.03, 50) diff --git a/docs/src/tutorials/time_dependent.md b/docs/src/tutorials/time_dependent.md index 6065960a..0bd13c9f 100644 --- a/docs/src/tutorials/time_dependent.md +++ b/docs/src/tutorials/time_dependent.md @@ -65,6 +65,7 @@ plot(time_evo, ["u1", "v1"], harmonic_eq) Let us compare this to the steady state diagram. ```@example time_dependent +fixed = (ω0 => 1.0, γ => 1e-2, λ => 5e-2, F => 1e-3, α => 1.0, η => 0.3, θ => 0) varied = ω => range(0.9, 1.1, 100) result = get_steady_states(harmonic_eq, varied, fixed) plot(result, "sqrt(u1^2 + v1^2)") diff --git a/examples/parametric_via_three_wave_mixing.jl b/examples/parametric_via_three_wave_mixing.jl index 3b86887c..4a084fea 100644 --- a/examples/parametric_via_three_wave_mixing.jl +++ b/examples/parametric_via_three_wave_mixing.jl @@ -22,7 +22,7 @@ harmonic_eq.equations varied = (ω => range(0.99, 1.1, 200)) # range of parameter values fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters -result = get_steady_states(harmonic_eq, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq, varied, fixed) plot(result; y="u1^2+v1^2") # If we set the cubic nonlinearity to zero, we recover the driven damped harmonic oscillator. Indeed, thefirst order the quadratic nonlinearity has no affect on the system. @@ -30,7 +30,7 @@ plot(result; y="u1^2+v1^2") varied = (ω => range(0.99, 1.1, 100)) fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) -result = get_steady_states(harmonic_eq, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq, varied, fixed) plot(result; y="u1^2+v1^2") # ## 2nd order Krylov expansion @@ -50,7 +50,7 @@ harmonic_eq2 = get_krylov_equations(diff_eq; order=2) varied = (ω => range(0.4, 1.1, 500)) fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.001, F => 0.005) -result = get_steady_states(harmonic_eq2, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq2, varied, fixed) plot(result; y="v1") # @@ -58,7 +58,6 @@ plot(result; y="v1") varied = (ω => range(0.4, 0.6, 100), F => range(1e-6, 0.01, 50)) fixed = (α => 1.0, β => 2.0, ω0 => 1.0, γ => 0.01) -result = get_steady_states( - harmonic_eq2, varied, fixed; threading=true, method=:total_degree -) +method = TotalDegree() +result = get_steady_states(harmonic_eq2, method, varied, fixed) plot_phase_diagram(result; class="stable") diff --git a/examples/parametron.jl b/examples/parametron.jl index e9cff811..0db5a953 100644 --- a/examples/parametron.jl +++ b/examples/parametron.jl @@ -50,7 +50,7 @@ varied = ω => range(0.9, 1.1, 100) result = get_steady_states(harmonic_eq, varied, fixed) -# In `get_steady_states`, the default value for the keyword `method=:random_warmup` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `:total_degree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower. The `threading` keyword enables parallel tracking of homotopy paths, and it's set to `false` simply because we are using a single core computer for now. +# In `get_steady_states`, the default method `WarmUp()` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `TotalDegree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower. # After solving the system, we can save the full output of the simulation and the model (e.g. symbolic expressions for the harmonic equations) into a file @@ -74,6 +74,7 @@ plot(result, "sqrt(u1^2 + v1^2)"; class="all") # The parametrically driven oscillator boasts a stability diagram called "Arnold's tongues" delineating zones where the oscillator is stable from those where it is exponentially unstable (if the nonlinearity was absence). We can retrieve this diagram by calculating the steady states as a function of external detuning $\delta=\omega_L-\omega_0$ and the parametric drive strength $\lambda$. # To perform a 2D sweep over driving frequency $\omega$ and parametric drive strength $\lambda$, we keep `fixed` from before but include 2 variables in `varied` +fixed = (ω₀ => 1.0, γ => 1e-2, F => 1e-3, α => 1.0, η => 0.3) varied = (ω => range(0.8, 1.2, 50), λ => range(0.001, 0.6, 50)) result_2D = get_steady_states(harmonic_eq, varied, fixed); diff --git a/examples/wave_mixing.jl b/examples/wave_mixing.jl index 2af140f2..db548ac0 100644 --- a/examples/wave_mixing.jl +++ b/examples/wave_mixing.jl @@ -30,7 +30,7 @@ harmonic_eq = get_harmonic_equations(diff_eq) varied = (ω => range(0.9, 1.2, 200)) # range of parameter values fixed = (α => 1.0, β => 0.0, ω0 => 1.0, γ => 0.005, F => 0.0025) # fixed parameters -result = get_steady_states(harmonic_eq, varied, fixed; threading=true)# compute steady states +result = get_steady_states(harmonic_eq, varied, fixed)# compute steady states p1 = plot(result; y="√(u1^2+v1^2)", legend=:best) p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1)) @@ -46,7 +46,7 @@ plot(p1, p2, p3; layout=(1, 3), size=(900, 300), margin=5mm) varied = (ω => range(0.9, 1.2, 200)) fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) -result = get_steady_states(harmonic_eq, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq, varied, fixed) p1 = plot(result; y="√(u1^2+v1^2)", legend=:best) p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1)) @@ -62,7 +62,7 @@ plot(p1, p2, p3; layout=(1, 3), size=(900, 300), margin=5mm) varied = (ω => range(0.9, 1.2, 200)) fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025) -result = get_steady_states(harmonic_eq, varied, fixed; threading=true) +result = get_steady_states(harmonic_eq, varied, fixed) p1 = plot(result; y="√(u1^2+v1^2)", legend=:best) p2 = plot(result; y="√(u2^2+v2^2)", legend=:best, ylims=(-0.1, 0.1)) diff --git a/src/DifferentialEquation.jl b/src/DifferentialEquation.jl index 01235ddc..01721e9e 100644 --- a/src/DifferentialEquation.jl +++ b/src/DifferentialEquation.jl @@ -1,3 +1,81 @@ +""" +$(TYPEDEF) + +Holds differential equation(s) of motion and a set of harmonics to expand each variable. +This is the primary input for `HarmonicBalance.jl` ; after inputting the equations, the harmonics + ansatz needs to be specified using `add_harmonic!`. + +# Fields +$(TYPEDFIELDS) + +## Example +```julia-repl +julia> @variables t, x(t), y(t), ω0, ω, F, k; + +# equivalent ways to enter the simple harmonic oscillator +julia> DifferentialEquation(d(x,t,2) + ω0^2 * x - F * cos(ω*t), x); +julia> DifferentialEquation(d(x,t,2) + ω0^2 * x ~ F * cos(ω*t), x); + +# two coupled oscillators, one of them driven +julia> DifferentialEquation([d(x,t,2) + ω0^2 * x - k*y, d(y,t,2) + ω0^2 * y - k*x] .~ [F * cos(ω*t), 0], [x,y]); +``` +""" +mutable struct DifferentialEquation + """Assigns to each variable an equation of motion.""" + equations::OrderedDict{Num,Equation} + """Assigns to each variable a set of harmonics.""" + harmonics::OrderedDict{Num,OrderedSet{Num}} + + function DifferentialEquation(eqs) + return new(eqs, OrderedDict(var => OrderedSet() for var in keys(eqs))) + end + + # uses the above constructor if no harmonics defined + function DifferentialEquation(eqs::Vector{Equation}, vars::Vector{Num}) + return DifferentialEquation(OrderedDict(zip(vars, eqs))) + end + + # if expressions are entered instead of equations, automatically set them = 0 + function DifferentialEquation(exprs::Vector{Num}, vars::Vector{Num}) + return DifferentialEquation(exprs .~ Int(0), vars) + end + + function DifferentialEquation(eq::Equation, var::Num) + typerhs = typeof(eq.rhs) + typelhs = typeof(eq.lhs) + if eq.rhs isa AbstractVector || eq.lhs isa AbstractVector + throw( + ArgumentError( + "The equation is of the form $(typerhs)~$(typelhs) is not supported. Commenly one forgot to broadcast the equation symbol `~`.", + ), + ) + end + return DifferentialEquation([eq], [var]) + end + function DifferentialEquation(eq::Equation, vars::Vector{Num}) + typerhs = typeof(eq.rhs) + typelhs = typeof(eq.lhs) + throw( + ArgumentError( + "The variables are of type $(typeof(vars)). Whereas you gave one equation of type $(typerhs)~$(typelhs). Commenly one forgot to broadcast the equation symbol `~`.", + ), + ) + end + DifferentialEquation(lhs::Num, var::Num) = DifferentialEquation([lhs ~ Int(0)], [var]) +end + +function Base.show(io::IO, diff_eq::DifferentialEquation) + println(io, "System of ", length(keys(diff_eq.equations)), " differential equations") + println(io, "Variables: ", join(keys(diff_eq.equations), ", ")) + print(io, "Harmonic ansatz: ") + for var in keys(diff_eq.harmonics) + print(io, string(var), " => ", join(string.(diff_eq.harmonics[var]), ", ")) + print(io, "; ") + end + println(io, "\n") + return [println(io, eq) for eq in values(diff_eq.equations)] +end + """ $(TYPEDSIGNATURES) Add the harmonic `ω` to the harmonic ansatz used to expand the variable `var` in `diff_eom`. diff --git a/src/HarmonicBalance.jl b/src/HarmonicBalance.jl index 9d0cb4fd..f8d22f61 100644 --- a/src/HarmonicBalance.jl +++ b/src/HarmonicBalance.jl @@ -13,6 +13,7 @@ using OrderedCollections: OrderedDict, OrderedSet using ProgressMeter: ProgressMeter, Progress using LinearAlgebra: eigvals using Random: Random # for setting seed +using EndpointRanges: EndpointRanges using Distances: Distances using BijectiveHilbert: BijectiveHilbert, Simple2D, decode_hilbert!, encode_hilbert @@ -40,21 +41,25 @@ using .ExprUtils: is_harmonic, substitute_all, drop_powers, count_derivatives include("extention_functions.jl") include("utils.jl") include("types.jl") - include("DifferentialEquation.jl") include("HarmonicVariable.jl") include("HarmonicEquation.jl") +include("Problem.jl") +include("Result.jl") +include("methods.jl") + include("solve_homotopy.jl") include("sorting.jl") include("classification.jl") + include("saving.jl") include("transform_solutions.jl") include("plotting_Plots.jl") -export show, *, @variables, d, ComplexF64, Float64, IM_TOL +export show, *, @variables, d, IM_TOL +export WarmUp, TotalDegree, Polyhedral -export ParameterRange, ParameterList, StateDict, SteadyState, ParameterVector -export DifferentialEquation, HarmonicVariable, HarmonicEquation, Problem, Result +export DifferentialEquation, HarmonicVariable, HarmonicEquation export get_steady_states, get_single_solution, get_harmonic_equations, add_harmonic! export get_variables, get_independent_variables, classify_branch, classify_solutions! export rearrange_standard @@ -86,14 +91,14 @@ using .FFTWExt using PrecompileTools: @setup_workload, @compile_workload -@setup_workload begin - # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the - # precompile file and potentially make loading faster. - @compile_workload begin - # all calls in this block will be precompiled, regardless of whether - # they belong to your package or not (on Julia 1.8 and higher) - include("precompilation.jl") - end -end +# @setup_workload begin +# # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the +# # precompile file and potentially make loading faster. +# @compile_workload begin +# # all calls in this block will be precompiled, regardless of whether +# # they belong to your package or not (on Julia 1.8 and higher) +# include("precompilation.jl") +# end +# end end # module diff --git a/src/HarmonicEquation.jl b/src/HarmonicEquation.jl index 7822e915..1a126851 100644 --- a/src/HarmonicEquation.jl +++ b/src/HarmonicEquation.jl @@ -1,3 +1,57 @@ + +""" +$(TYPEDEF) + +Holds a set of algebraic equations governing the harmonics of a `DifferentialEquation`. + +# Fields +$(TYPEDFIELDS) +""" +mutable struct HarmonicEquation + """A set of equations governing the harmonics.""" + equations::Vector{Equation} + """A set of variables describing the harmonics.""" + variables::Vector{HarmonicVariable} + """The parameters of the equation set.""" + parameters::Vector{Num} + "The natural equation (before the harmonic ansatz was used)." + natural_equation::DifferentialEquation + + # use a self-referential constructor with _parameters + function HarmonicEquation(equations, variables, nat_eq) + return (x = new(equations, variables, Vector{Num}([]), nat_eq); + x.parameters = _parameters(x); + x) + end + function HarmonicEquation(equations, variables, parameters, natural_equation) + return new(equations, variables, parameters, natural_equation) + end +end + +function Base.show(io::IO, eom::HarmonicEquation) + println(io, "A set of ", length(eom.equations), " harmonic equations") + println(io, "Variables: ", join(string.(get_variables(eom)), ", ")) + println(io, "Parameters: ", join(string.(eom.parameters), ", ")) + println(io, "\nHarmonic ansatz: ", _show_ansatz(eom)) + println(io, "\nHarmonic equations:") + return [println(io, "\n", eq) for eq in eom.equations] +end + +"""Gives the full harmonic ansatz used to construct `eom`.""" +function _show_ansatz(eom::HarmonicEquation) + output = "" + vars = unique(getfield.(eom.variables, :natural_variable)) + for nat_var in vars + # the Hopf variable (limit cycle frequency) does not contribute a term + harm_vars = filter( + x -> isequal(nat_var, x.natural_variable) && x.type !== "Hopf", eom.variables + ) + ansatz = join([_show_ansatz(var) for var in harm_vars], " + ") + output *= "\n" * string(nat_var) * " = " * ansatz + end + return output +end + Base.show(eom::HarmonicEquation) = show_fields(eom) """ @@ -140,9 +194,6 @@ function Symbolics.get_variables(eom::HarmonicEquation)::Vector{Num} return get_variables.(eom.variables) end -Symbolics.get_variables(p::Problem)::Vector{Num} = get_variables(p.eom) -Symbolics.get_variables(res::Result)::Vector{Num} = get_variables(res.problem) - "Get the parameters (not time nor variables) of a HarmonicEquation" function _parameters(eom::HarmonicEquation) all_symbols = flatten([ diff --git a/src/HarmonicVariable.jl b/src/HarmonicVariable.jl index 332f987a..ebed7618 100644 --- a/src/HarmonicVariable.jl +++ b/src/HarmonicVariable.jl @@ -1,3 +1,44 @@ +""" +$(TYPEDEF) + +Holds a variable stored under `symbol` describing the harmonic `ω` of `natural_variable`. + +# Fields +$(TYPEDFIELDS) +""" +mutable struct HarmonicVariable + """Symbol of the variable in the HarmonicBalance namespace.""" + symbol::Num + """Human-readable labels of the variable, used for plotting.""" + name::String + """Type of the variable (u or v for quadratures, a for a constant, Hopf for Hopf etc.)""" + type::String + """The harmonic being described.""" + ω::Num + """The natural variable whose harmonic is being described.""" + natural_variable::Num +end + +function Base.show(io::IO, hv::HarmonicVariable) + return println( + io, + "Harmonic variable ", + string.(hv.symbol) * " for harmonic ", + string(hv.ω), + " of ", + string(hv.natural_variable), + ) +end + +"""Gives the relation between `var` and the underlying natural variable.""" +function _show_ansatz(var::HarmonicVariable) + t = var.natural_variable.val.arguments + t = length(t) == 1 ? string(t[1]) : error("more than 1 independent variable") + ω = string(var.ω) + terms = Dict("u" => "*cos(" * ω * t * ")", "v" => "*sin(" * ω * t * ")", "a" => "") + return string(string(var.symbol) * terms[var.type]) +end + # pretty-printing Base.display(var::HarmonicVariable) = display(var.name) Base.display(var::Vector{HarmonicVariable}) = display.(getfield.(var, Symbol("name"))) diff --git a/src/LimitCycles/LimitCycles.jl b/src/LimitCycles/LimitCycles.jl index 1babaf59..8e6ba97a 100644 --- a/src/LimitCycles/LimitCycles.jl +++ b/src/LimitCycles/LimitCycles.jl @@ -6,6 +6,11 @@ using Symbolics: Symbolics, Num, expand_derivatives using HarmonicBalance using HarmonicBalance: + HarmonicBalanceMethod, + WarmUp, + Problem, + Result, + get_steady_states, order_branches!, find_branch_order, _remove_brackets, @@ -15,6 +20,7 @@ using HarmonicBalance: _is_physical, substitute_all, var_name + using HarmonicBalance.LinearResponse: get_implicit_Jacobian using HarmonicBalance.ExprUtils: get_all_terms diff --git a/src/LimitCycles/gauge_fixing.jl b/src/LimitCycles/gauge_fixing.jl index f84f3936..b579c742 100644 --- a/src/LimitCycles/gauge_fixing.jl +++ b/src/LimitCycles/gauge_fixing.jl @@ -85,16 +85,7 @@ function _choose_fixed(eom, ω_lc) return first(vars) # This is arbitrary; better would be to substitute with values end -""" - get_limit_cycles(eom::HarmonicEquation, swept, fixed, ω_lc; kwargs...) - -Variant of `get_steady_states` for a limit cycle problem characterised by a Hopf frequency (usually called ω_lc) - -Solutions with ω_lc = 0 are labelled unphysical since this contradicts the assumption of distinct harmonic variables corresponding to distinct harmonics. -""" -function get_limit_cycles( - eom::HarmonicEquation, swept, fixed, ω_lc; explicit_Jacobian=false, kwargs... -) +function limit_cycle_problem(eom::HarmonicEquation, swept, fixed, ω_lc, explicit_Jacobian) prob = _cycle_Problem(eom, ω_lc) prob.jacobian = _gaugefixed_Jacobian( eom, @@ -103,20 +94,53 @@ function get_limit_cycles( sym_order=_free_symbols(prob, swept), rules=fixed, ) + return prob +end - result = get_steady_states( - prob, swept, fixed; method=:warmup, threading=true, classify_default=true, kwargs... - ) +""" + get_limit_cycles( + eom::HarmonicEquation, method::HarmonicBalanceMethod, swept, fixed, ω_lc; kwargs...) - _classify_limit_cycles!(result, ω_lc) +Variant of `get_steady_states` for a limit cycle problem characterised by a Hopf frequency (usually called ω_lc) - return result +Solutions with ω_lc = 0 are labelled unphysical since this contradicts the assumption of distinct harmonic variables corresponding to distinct harmonics. +""" +function get_limit_cycles( + eom::HarmonicEquation, swept, fixed, ω_lc; explicit_Jacobian=false, kwargs... +) + prob = limit_cycle_problem(eom, swept, fixed, ω_lc, explicit_Jacobian) + return get_limit_cycles(prob, WarmUp(), swept, fixed, ω_lc; kwargs...) end - function get_limit_cycles( - eom::HarmonicEquation, swept, fixed; limit_cycle_harmonic, kwargs... + eom::HarmonicEquation, + method::HarmonicBalanceMethod, + swept, + fixed, + ω_lc; + explicit_Jacobian=false, + kwargs..., ) - return get_limit_cycles(eom, swept, fixed, limit_cycle_harmonic; kwargs...) + prob = limit_cycle_problem(eom, swept, fixed, ω_lc, explicit_Jacobian) + return get_limit_cycles(prob, method, swept, fixed, ω_lc; kwargs...) +end +function get_limit_cycles(eom::HarmonicEquation, pairs::Dict, ω_lc; kwargs...) + swept = filter(x -> length(x[2]) > 1, pairs) + fixed = filter(x -> length(x[2]) == 1, pairs) + return get_limit_cycles(eom, swept, fixed, ω_lc; kwargs...) +end +function get_limit_cycles( + eom::HarmonicEquation, method::HarmonicBalanceMethod, pairs::Dict, ω_lc; kwargs... +) + swept = filter(x -> length(x[2]) > 1, pairs) + fixed = filter(x -> length(x[2]) == 1, pairs) + return get_limit_cycles(eom, method, swept, fixed, ω_lc; kwargs...) +end +function get_limit_cycles( + prob::Problem, method::HarmonicBalanceMethod, swept, fixed, ω_lc; kwargs... +) + result = get_steady_states(prob, method, swept, fixed; kwargs...) + _classify_limit_cycles!(result, ω_lc) + return result end # if abs(ω_lc) < tol, set all classifications to false diff --git a/src/LinearResponse/LinearResponse.jl b/src/LinearResponse/LinearResponse.jl index c3390fbe..4d8a6eb5 100644 --- a/src/LinearResponse/LinearResponse.jl +++ b/src/LinearResponse/LinearResponse.jl @@ -13,8 +13,18 @@ using LinearAlgebra: norm, eigvals, eigen, eigvecs using OrderedCollections: OrderedDict using HarmonicBalance +using HarmonicBalance: + Result, + HarmonicVariable, + HarmonicEquation, + DifferentialEquation, + StateDict, + get_variables, + get_independent_variables + using HarmonicBalance: var_name, + d, rearrange_standard, _remove_brackets, expand_derivatives, @@ -30,6 +40,7 @@ using HarmonicBalance: fourier_transform, declare_variable, is_rearranged + using ..HC_wrapper include("types.jl") diff --git a/src/LinearResponse/Lorentzian_spectrum.jl b/src/LinearResponse/Lorentzian_spectrum.jl index 44c1ebce..2fab36c7 100644 --- a/src/LinearResponse/Lorentzian_spectrum.jl +++ b/src/LinearResponse/Lorentzian_spectrum.jl @@ -80,9 +80,6 @@ function JacobianSpectrum(res::Result; index::Int, branch::Int, force=false) solution_dict = get_single_solution(res; branch=branch, index=index) λs, vs = eigen(res.jacobian(solution_dict)) - # convert OrderedDict to Dict - see Symbolics issue #601 - solution_dict = Dict(get_single_solution(res; index=index, branch=branch)) - for (j, λ) in enumerate(λs) eigvec = vs[:, j] # the eigenvector @@ -90,7 +87,7 @@ function JacobianSpectrum(res::Result; index::Int, branch::Int, force=false) for pair in _get_uv_pairs(hvars) u, v = hvars[pair] eigvec_2d = eigvec[pair] # fetch the relevant part of the Jacobian eigenvector - ωnum = real(unwrap(Symbolics.substitute(u.ω, solution_dict))) + ωnum = real(unwrap(Symbolics.substitute(u.ω, Dict(solution_dict)))) # ^ the harmonic (numerical now) associated to this harmonic variable # eigvec_2d is associated to a natural variable -> this variable gets Lorentzian peaks diff --git a/src/LinearResponse/jacobians.jl b/src/LinearResponse/jacobians.jl index 15e5ae15..706d8158 100644 --- a/src/LinearResponse/jacobians.jl +++ b/src/LinearResponse/jacobians.jl @@ -21,7 +21,7 @@ end " Obtain a Jacobian from a `DifferentialEquation` by first converting it into a `HarmonicEquation`. " function get_Jacobian(diff_eom::DifferentialEquation) - @variables T + Symbolics.@variables T harmonic_eq = get_harmonic_equations( diff_eom; slow_time=T, fast_time=first(get_independent_variables(diff_eom)) ) diff --git a/src/LinearResponse/response.jl b/src/LinearResponse/response.jl index 312dcb68..376d5743 100644 --- a/src/LinearResponse/response.jl +++ b/src/LinearResponse/response.jl @@ -7,7 +7,7 @@ This routine cannot accept a `HarmonicEquation` since there, some time-derivativ """ function get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)::Matrix - @variables T, i + Symbolics.@variables T, i time = get_independent_variables(diff_eq)[1] eom = harmonic_ansatz(diff_eq, time) @@ -33,7 +33,7 @@ Any substitution rules not specified in `res` can be supplied in `rules`." function ResponseMatrix(res::Result; rules=Dict()) # get the symbolic response matrix - @variables Δ + Symbolics.@variables Δ M = get_response_matrix(res.problem.eom.natural_equation, Δ; order=2) M = substitute_all(M, merge(res.fixed_parameters, rules)) symbols = _free_symbols(res) diff --git a/src/Problem.jl b/src/Problem.jl new file mode 100644 index 00000000..3d47115e --- /dev/null +++ b/src/Problem.jl @@ -0,0 +1,70 @@ + +""" +$(TYPEDEF) + +Holds a set of algebraic equations describing the steady state of a system. + +# Fields +$(TYPEDFIELDS) + +# Constructors +```julia +Problem(eom::HarmonicEquation; Jacobian=true) # find and store the symbolic Jacobian +Problem(eom::HarmonicEquation; Jacobian="implicit") # ignore the Jacobian for now, compute implicitly later +Problem(eom::HarmonicEquation; Jacobian=J) # use J as the Jacobian (a function that takes a Dict) +Problem(eom::HarmonicEquation; Jacobian=false) # ignore the Jacobian +``` +""" +mutable struct Problem + "The harmonic variables to be solved for." + variables::Vector{Num} + "All symbols which are not the harmonic variables." + parameters::Vector{Num} + "The input object for HomotopyContinuation.jl solver methods." + system::HC.System + "The Jacobian matrix (possibly symbolic). + If `false`, the Jacobian is ignored (may be calculated implicitly after solving)." + jacobian + "The HarmonicEquation object used to generate this `Problem`." + eom::HarmonicEquation + + function Problem(variables, parameters, system, jacobian) + return new(variables, parameters, system, jacobian) + end #incomplete initialization for user-defined symbolic systems + function Problem(variables, parameters, system, jacobian, eom) + return new(variables, parameters, system, jacobian, eom) + end +end + +Symbolics.get_variables(p::Problem)::Vector{Num} = get_variables(p.eom) + +function Base.show(io::IO, p::Problem) + println(io, length(p.system.expressions), " algebraic equations for steady states") + println(io, "Variables: ", join(string.(p.variables), ", ")) + println(io, "Parameters: ", join(string.(p.parameters), ", ")) + return println(io, "Symbolic Jacobian: ", !(p.jacobian == false)) +end + +# assume this order of variables in all compiled function (transform_solutions, Jacobians) +function _free_symbols(p::Problem, varied) + return cat(p.variables, collect(keys(OrderedDict(varied))); dims=1) +end + +""" Compile the Jacobian from `prob`, inserting `fixed_parameters`. + Returns a function that takes a dictionary of variables and `swept_parameters` to give the Jacobian.""" +function _compile_Jacobian( + prob::Problem, swept_parameters::ParameterRange, fixed_parameters::ParameterList +) + if prob.jacobian isa Matrix + compiled_J = compile_matrix( + prob.jacobian, _free_symbols(prob, swept_parameters); rules=fixed_parameters + ) + elseif prob.jacobian == "implicit" + compiled_J = LinearResponse.get_implicit_Jacobian( + prob, swept_parameters, fixed_parameters + ) # leave implicit Jacobian as is + else + return prob.jacobian + end + return compiled_J +end diff --git a/src/Result.jl b/src/Result.jl new file mode 100644 index 00000000..f9705f97 --- /dev/null +++ b/src/Result.jl @@ -0,0 +1,60 @@ +""" +$(TYPEDEF) + +Stores the steady states of a HarmonicEquation. + +# Fields +$(TYPEDFIELDS) + +""" +mutable struct Result + "The variable values of steady-state solutions." + solutions::Array{Vector{SteadyState}} + "Values of all parameters for all solutions." + swept_parameters::ParameterRange + "The parameters fixed throughout the solutions." + fixed_parameters::ParameterList + "The `Problem` used to generate this." + problem::Problem + "Maps strings such as \"stable\", \"physical\" etc to arrays of values, classifying the solutions (see method `classify_solutions!`)." + classes::Dict{String,Array} + "The Jacobian with `fixed_parameters` already substituted. Accepts a dictionary specifying the solution. + If problem.jacobian is a symbolic matrix, this holds a compiled function. + If problem.jacobian was `false`, this holds a function that rearranges the equations to find J + only after numerical values are inserted (preferable in cases where the symbolic J would be very large)." + jacobian::Function + "Seed used for the solver" + seed::UInt32 + + function Result(sol, swept, fixed, problem, classes, J, seed) + return new(sol, swept, fixed, problem, classes, J, seed) + end + Result(sol, swept, fixed, problem, classes) = new(sol, swept, fixed, problem, classes) + Result(sol, swept, fixed, problem) = new(sol, swept, fixed, problem, Dict([])) +end + +Symbolics.get_variables(res::Result)::Vector{Num} = get_variables(res.problem) + +function Base.show(io::IO, r::Result) + println(io, "A steady state result for ", length(r.solutions), " parameter points") + println(io, "\nSolution branches: ", length(r.solutions[1])) + println( + io, " of which real: ", sum(push!(any.(classify_branch(r, "physical")), false)) + ) + println( + io, " of which stable: ", sum(push!(any.(classify_branch(r, "stable")), false)) + ) + return println(io, "\nClasses: ", join(keys(r.classes), ", ")) +end + +# assume this order of variables in all compiled function (transform_solutions, Jacobians) +function _free_symbols(res::Result) + return cat(res.problem.variables, collect(keys(res.swept_parameters)); dims=1) +end + +# overload to use [] for indexing +Base.getindex(r::Result, idx::Int...) = get_single_solution(r, idx) +Base.size(r::Result) = size(r.solutions) + +branch_count(r::Result) = length(r.solutions[1]) +get_branch(r::Result, idx) = getindex.(r.solutions, idx) diff --git a/src/methods.jl b/src/methods.jl new file mode 100644 index 00000000..7d29b8d1 --- /dev/null +++ b/src/methods.jl @@ -0,0 +1,224 @@ +""" + HarmonicBalanceMethod + +Abstract type for harmonic balance methods. +""" +abstract type HarmonicBalanceMethod end + +""" + TotalDegree + +The Total Degree homotopy method. Performs a homotopy ``H(x, t) = γ t G(x) + (1-t) F(x)`` +from the trivial polynomial system ``xᵢ^{dᵢ} +aᵢ`` with the maximal degree ``dᵢ`` determined +by the [Bezout bound](https://en.wikipedia.org/wiki/B%C3%A9zout%27s_theorem). The method +guarantees to find all solutions, however, it comes with a high computational cost. See +[HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/guides/totaldegree/) +for more information. + +# Fields +$(TYPEDFIELDS) +""" +Base.@kwdef struct TotalDegree <: HarmonicBalanceMethod + """Complex multiplying factor of the start system G(x) for the homotopy""" + gamma::Complex = cis(2π * rand(Random.MersenneTwister(0xd8e5d8df))) + """Boolean indicating if threading is enabled.""" + thread::Bool = Threads.nthreads() > 1 + """Options for the tracker.""" + tracker_options::HomotopyContinuation.TrackerOptions = HomotopyContinuation.TrackerOptions() + """Options for the endgame.""" + endgame_options::HomotopyContinuation.EndgameOptions = HomotopyContinuation.EndgameOptions() + """Compilation options.""" + compile::Union{Bool,Symbol} = HomotopyContinuation.COMPILE_DEFAULT[] + """Seed for random number generation.""" + seed::UInt32 = 0xd8e5d8df +end + +""" + Polyhedral + +The Polyhedral homotopy method. This method constructs a homotopy based on the polyhedral +structure of the polynomial system. It is more efficient than the Total Degree method for +sparse systems, meaning most of the coefficients are zero. It can be especially useful if +you don't need to find the zero solutions (`only_non_zero = true`), resulting in speed up. +See [HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/guides/polyhedral/) +for more information. + +# Fields +$(TYPEDFIELDS) +""" +Base.@kwdef struct Polyhedral <: HarmonicBalanceMethod + """Boolean indicating if only non-zero solutions are considered.""" + only_non_zero::Bool = false + """Boolean indicating if threading is enabled.""" + thread::Bool = Threads.nthreads() > 1 + """Options for the tracker.""" + tracker_options::HomotopyContinuation.TrackerOptions = HomotopyContinuation.TrackerOptions() + """Options for the endgame.""" + endgame_options::HomotopyContinuation.EndgameOptions = HomotopyContinuation.EndgameOptions() + """Compilation options.""" + compile::Union{Bool,Symbol} = HomotopyContinuation.COMPILE_DEFAULT[] + """Seed for random number generation.""" + seed::UInt32 = 0xd8e5d8df +end + +""" + WarmUp + +The Warm Up method. This method prepares a warmup system using the parameter at `index` +perturbed by `perturbation_size` and performs a homotopy using the warmup system to all other +systems in the parameter sweep. It is very efficient for systems with less bifurcation in the +parameter sweep. The Warm Up method does not guarantee to find all solutions. See +[HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/guides/many-systems/) +for more information. + +# Fields +$(TYPEDFIELDS) +""" +Base.@kwdef struct WarmUp <: HarmonicBalanceMethod + """Size of the perturbation.""" + perturbation_size::ComplexF64 = 1e-6 + 1e-6 * im + """Index for the endpoint.""" + index::Union{Int,EndpointRanges.Endpoint} = EndpointRanges.iend ÷ 2 + """Boolean indicating if threading is enabled.""" + thread::Bool = Threads.nthreads() > 1 + """Options for the tracker.""" + tracker_options::HomotopyContinuation.TrackerOptions = HomotopyContinuation.TrackerOptions() + """Options for the endgame.""" + endgame_options::HomotopyContinuation.EndgameOptions = HomotopyContinuation.EndgameOptions() + """Compilation options.""" + compile::Union{Bool,Symbol} = HomotopyContinuation.COMPILE_DEFAULT[] + """Seed for random number generation.""" + seed::UInt32 = 0xd8e5d8df +end + +""" + thread(method::HarmonicBalanceMethod) -> Bool + +Returns whether threading is enabled for the given method. The number of available threads +is controlled by the environment variable `JULIA_NUM_THREADS`. +""" +thread(method::HarmonicBalanceMethod) = method.thread + +""" + tracker(method::HarmonicBalanceMethod) -> HomotopyContinuation.TrackerOptions + +Returns the tracker options for the given method. See `HomotopyContinuation.TrackerOptions` +for the available options. +""" +tracker(method::HarmonicBalanceMethod) = method.tracker_options + +""" + endgame(method::HarmonicBalanceMethod) -> HomotopyContinuation.EndgameOptions + +Returns the endgame options for the given method. See `HomotopyContinuation.EndgameOptions` +for the available options. +""" +endgame(method::HarmonicBalanceMethod) = method.endgame_options + +""" + compile(method::HarmonicBalanceMethod) -> Union{Bool,Symbol} + +Returns the compile options for the given method. If `true` then a system is compiled to a +straight line program for evaluation. This induces a compilation overhead. If `false` then +the generated program is only interpreted. This is slower than the compiled version, +but does not introduce compilation overhead. +""" +compile(method::HarmonicBalanceMethod) = method.compile + +""" + seed(method::HarmonicBalanceMethod) -> UInt32 + +Returns the seed for random number generation for the given method. +""" +seed(method::HarmonicBalanceMethod) = method.seed + +""" + alg_default_options(method::HarmonicBalanceMethod) -> NamedTuple + +Returns a named tuple of default algorithm options for the given method. +""" +function alg_default_options(method::HarmonicBalanceMethod) + return ( + threading=thread(method), + tracker_options=tracker(method), + endgame_options=endgame(method), + compile=compile(method), + seed=seed(method), + ) +end + +""" + alg_specific_options(method::TotalDegree) -> NamedTuple + +Returns a named tuple of specific algorithm options for the Total Degree method. +""" +alg_specific_options(method::TotalDegree) = (gamma=method.gamma,) + +""" + alg_specific_options(method::Polyhedral) -> NamedTuple + +Returns a named tuple of specific algorithm options for the Polyhedral method. +""" +alg_specific_options(method::Polyhedral) = (only_non_zero=method.only_non_zero,) + +""" + alg_specific_options(method::WarmUp) -> NamedTuple + +Returns a named tuple of specific algorithm options for the Warm Up method. +""" +function alg_specific_options(method::WarmUp) + return (perturbation_size=method.perturbation_size, index=method.index) +end + +""" + method_symbol(m::Polyhedral) -> Symbol + +Returns the symbol for the Polyhedral method identified by HomotopyContinuation/jl. +""" +method_symbol(m::Polyhedral) = :polyhedral + +""" + method_symbol(m::TotalDegree) -> Symbol + +Returns the symbol for the Total Degree method identified by HomotopyContinuation.jl. +""" +method_symbol(m::TotalDegree) = :total_degree + +""" + Base.show(io::IO, m::WarmUp) + +Displays information about the Warm Up method. +""" +function Base.show(io::IO, m::WarmUp) + println(io, "Warm up method:") + println(io, "perturbation_size: ", m.perturbation_size) + println(io, "Threading: ", thread(m)) + println(io, "Compile: ", compile(m)) + return println(io, "Seed: ", seed(m)) +end + +""" + Base.show(io::IO, m::TotalDegree) + +Displays information about the Total Degree method. +""" +function Base.show(io::IO, m::TotalDegree) + println(io, "Total degree method:") + println(io, "Gamma: ", m.gamma) + println(io, "Threading: ", thread(m)) + println(io, "Compile: ", compile(m)) + return println(io, "Seed: ", seed(m)) +end + +""" + Base.show(io::IO, m::Polyhedral) + +Displays information about the Polyhedral method. +""" +function Base.show(io::IO, m::Polyhedral) + println(io, "Polyhedral method:") + println(io, "Zero solutions: ", !m.only_non_zero) + println(io, "Threading: ", thread(m)) + println(io, "Compile: ", compile(m)) + return println(io, "Seed: ", seed(m)) +end diff --git a/src/precompilation.jl b/src/precompilation.jl index dd269803..36c98a1e 100644 --- a/src/precompilation.jl +++ b/src/precompilation.jl @@ -12,8 +12,7 @@ dEOM = DifferentialEquation(natural_equation + forces, x) add_harmonic!(dEOM, x, ω) harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); -prob = HarmonicBalance.Problem(harmonic_eq) fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 0, α => 1.0, η => 0.3, θ => 0, ψ => 0) varied = ω => range(0.9, 1.1, 20) -res = get_steady_states(prob, varied, fixed; show_progress=false) +res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) diff --git a/src/solve_homotopy.jl b/src/solve_homotopy.jl index 2bed7e03..f5e198ff 100644 --- a/src/solve_homotopy.jl +++ b/src/solve_homotopy.jl @@ -1,25 +1,25 @@ """ - get_steady_states(prob::Problem, + get_steady_states(problem::HarmonicEquation, + method::HarmonicBalanceMethod, swept_parameters::ParameterRange, fixed_parameters::ParameterList; - method=:warmup, - threading = Threads.nthreads() > 1, show_progress=true, - sorting="nearest") + sorting="nearest", + classify_default=true) -Solves `prob` over the ranges specified by `swept_parameters`, keeping `fixed_parameters` constant. -`swept_parameters` accepts pairs mapping symbolic variables to arrays or `LinRange`. +Solves `problem` with the `method` over the ranges specified by `swept_parameters`, +keeping `fixed_parameters` constant. +`swept_parameters` accepts pairs mapping symbolic variables to arrays or ranges. `fixed_parameters` accepts pairs mapping symbolic variables to numbers. -Keyword arguments -- `method`: If `:warmup` (default), a problem similar to `prob` but with random complex parameters is first solved to find all non-singular paths. The subsequent tracking to find results for all `swept_parameters` is then much faster than the initial solving. If `method=:total_degree`, each parameter point is solved separately by tracking the maximum number of paths (employs a total degree homotopy). -This takes far longer but can be more reliable. -- `threading`: If `true`, multithreaded support is activated. The number of available threads is set by the environment variable `JULIA_NUM_THREADS`. -- `sorting`: the method used by `sort_solutions` to get continuous solutions branches. The current options are `"hilbert"` (1D sorting along a Hilbert curve), `"nearest"` (nearest-neighbor sorting) and `"none"`. +### Keyword arguments - `show_progress`: Indicate whether a progress bar should be displayed. +- `sorting`: the method used by `sort_solutions` to get continuous solutions branches. The current options are `"hilbert"` (1D sorting along a Hilbert curve), `"nearest"` (nearest-neighbor sorting) and `"none"`. +- `classify_default`: If `true`, the solutions will be classified using the default classification method. -Example: solving a simple harmonic oscillator ``m \\ddot{x} + γ \\dot{x} + ω_0^2 x = F \\cos(ωt)`` -to obtain the response as a function of ``ω`` +### Example +solving a simple harmonic oscillator +``m \\ddot{x} + γ \\dot{x} + ω_0^2 x = F \\cos(ωt)`` to obtain the response as a function of ``ω`` ```julia-repl # having obtained a Problem object, let's find steady states julia> range = ParameterRange(ω => LinRange(0.8,1.2,100) ) # 100 parameter sets to solve @@ -36,10 +36,10 @@ A steady state result for 100 parameter points ``` -It is also possible to create multi-dimensional solutions plots. +It is also possible to perfrom 2-dimensional sweeps. ```julia-repl # The swept parameters take precedence over fixed -> use the same fixed -julia> range = ParameterRange(ω => LinRange(0.8,1.2,100), F => LinRange(0.1,1.0,10) ) # 100x10 parameter sets +julia> range = ParameterRange(ω => range(0.8,1.2,100), F => range(0.1,1.0,10) ) # The swept parameters take precedence over fixed -> the F in fixed is now ignored julia> get_steady_states(problem, range, fixed) @@ -56,98 +56,77 @@ A steady state result for 1000 parameter points """ function get_steady_states( prob::Problem, + method::HarmonicBalanceMethod, swept_parameters::ParameterRange, fixed_parameters::ParameterList; - method=:warmup, - threading=Threads.nthreads() > 1, show_progress=true, sorting="nearest", classify_default=true, - seed=nothing, - kwargs..., ) - - # set seed if provided - !isnothing(seed) && Random.seed!(seed) + Random.seed!(seed(method)) # make sure the variables are in our namespace to make them accessible later declare_variable.(string.(cat(prob.parameters, prob.variables; dims=1))) - # prepare an array of vectors, each representing one set of input parameters - # an n-dimensional sweep uses an n-dimensional array - unique_fixed = filter_duplicate_parameters(swept_parameters, fixed_parameters) - variable_names = var_name.([keys(fixed_parameters)..., keys(swept_parameters)...]) any([var_name(var) ∈ variable_names for var in get_variables(prob)]) && error("Cannot fix one of the variables!") - input_array = _prepare_input_params(prob, swept_parameters, unique_fixed) - # feed the array into HomotopyContinuation, get back an similar array of solutions - raw = _get_raw_solution( - prob, - input_array; - sweep=swept_parameters, - method=method, - threading=threading, - show_progress=show_progress, - seed=seed, - kwargs..., + unique_fixed, input_array = _prepare_input_params( + prob, swept_parameters, fixed_parameters ) - - # extract all the information we need from results - #rounded_solutions = unique_points.(HomotopyContinuation.solutions.(getindex.(raw, 1)); metric = EuclideanNorm(), atol=1E-14, rtol=1E-8) - rounded_solutions = HC.solutions.(getindex.(raw, 1)) - all(isempty.(rounded_solutions)) ? error("No solutions found!") : nothing - solutions = pad_solutions(rounded_solutions) + solutions = get_solutions(prob, method, input_array; show_progress=show_progress) compiled_J = _compile_Jacobian(prob, swept_parameters, unique_fixed) - # a "raw" solution struct result = Result( - solutions, swept_parameters, unique_fixed, prob, Dict(), compiled_J, seed + solutions, swept_parameters, unique_fixed, prob, Dict(), compiled_J, seed(method) ) - # sort into branches if sorting != "no_sorting" sort_solutions!(result; sorting=sorting, show_progress=show_progress) - else - nothing end classify_default ? _classify_default!(result) : nothing return result end -function get_steady_states(p::Problem, swept, fixed; kwargs...) - return get_steady_states(p, ParameterRange(swept), ParameterList(fixed); kwargs...) -end -function get_steady_states(eom::HarmonicEquation, swept, fixed; kwargs...) - return get_steady_states(Problem(eom), swept, fixed; kwargs...) -end -function get_steady_states(p, pairs::Dict; kwargs...) +function get_steady_states( + p::Problem, method::HarmonicBalanceMethod, swept, fixed; kwargs... +) return get_steady_states( - p, - filter(x -> length(x[2]) > 1, pairs), - filter(x -> length(x[2]) == 1, pairs); - kwargs..., + p, method, ParameterRange(swept), ParameterList(fixed); kwargs... ) end - -""" Compile the Jacobian from `prob`, inserting `fixed_parameters`. - Returns a function that takes a dictionary of variables and `swept_parameters` to give the Jacobian.""" -function _compile_Jacobian( - prob::Problem, swept_parameters::ParameterRange, fixed_parameters::ParameterList +function get_steady_states( + eom::HarmonicEquation, method::HarmonicBalanceMethod, swept, fixed; kwargs... ) - if prob.jacobian isa Matrix - compiled_J = compile_matrix( - prob.jacobian, _free_symbols(prob, swept_parameters); rules=fixed_parameters - ) - elseif prob.jacobian == "implicit" - compiled_J = LinearResponse.get_implicit_Jacobian( - prob, swept_parameters, fixed_parameters - ) # leave implicit Jacobian as is + return get_steady_states(Problem(eom), method, swept, fixed; kwargs...) +end +function get_steady_states(eom::HarmonicEquation, pairs::Dict; kwargs...) + swept = filter(x -> length(x[2]) > 1, pairs) + fixed = filter(x -> length(x[2]) == 1, pairs) + return get_steady_states(eom, swept, fixed; kwargs...) +end +function get_steady_states( + eom::HarmonicEquation, method::HarmonicBalanceMethod, pairs::Dict; kwargs... +) + swept = filter(x -> length(x[2]) > 1, pairs) + fixed = filter(x -> length(x[2]) == 1, pairs) + return get_steady_states(eom, method, swept, fixed; kwargs...) +end +function get_steady_states(eom::HarmonicEquation, swept, fixed; kwargs...) + return get_steady_states(Problem(eom), WarmUp(), swept, fixed; kwargs...) +end + +function get_solutions(prob, method, input_array; show_progress) + raw = _get_raw_solution(prob, method, input_array; show_progress=show_progress) + + solutions = HC.solutions.(getindex.(raw, 1)) + if all(isempty.(solutions)) + @warn "No solutions found!" + return solutions else - return prob.jacobian + pad_solutions(solutions) end - return compiled_J end """ @@ -176,7 +155,7 @@ function _prepare_input_params( prob::Problem, sweeps::ParameterRange, fixed_parameters::ParameterList ) # sweeping takes precedence over fixed_parameters - fixed_parameters = filter_duplicate_parameters(sweeps, fixed_parameters) + unique_fixed = filter_duplicate_parameters(sweeps, fixed_parameters) # fixed order of parameters all_keys = cat(collect(keys(sweeps)), collect(keys(fixed_parameters)); dims=1) @@ -200,7 +179,7 @@ function _prepare_input_params( # order each parameter vector to match the order in prob input_array = getindex.(input_array, [permutation]) # HC wants arrays, not tuples - return tuple_to_vector.(input_array) + return unique_fixed, tuple_to_vector.(input_array) end "Remove occurrences of `sweeps` elements from `fixed_parameters`." @@ -213,77 +192,70 @@ function filter_duplicate_parameters(sweeps, fixed_parameters) end "A random warmup solution is computed to use as `start_parameters` in the homotopy." -function _solve_warmup(problem::Problem, params_1D, sweep; threading, show_progress) +function _solve_warmup(problem::Problem, method::WarmUp, params; show_progress) # complex perturbation of the warmup parameters - complex_pert = [1e-5 * randn(ComplexF64) for p in problem.parameters] - real_pert = ones(length(params_1D[1])) - warmup_parameters = params_1D[end ÷ 2] .* (real_pert + complex_pert) + options = alg_specific_options(method) + l = length(params[1]) + + perturbation = ones(l) + options[:perturbation_size] * randn(ComplexF64, l) + warmup_parameters = params[options[:index]] .* perturbation warmup_solution = HC.solve( problem.system; start_system=:total_degree, target_parameters=warmup_parameters, - threading=threading, show_progress=show_progress, + alg_default_options(method)..., ) return warmup_parameters, warmup_solution end "Uses HomotopyContinuation to solve `problem` at specified `parameter_values`." function _get_raw_solution( - problem::Problem, - parameter_values; - sweep=ParameterRange(), - method=:warmup, - threading=false, - show_progress=true, - seed=nothing, - kwargs..., + problem::Problem, method::WarmUp, parameter_values; show_progress ) - # HomotopyContinuation accepts 1D arrays of parameter sets - params_1D = reshape(parameter_values, :, 1) + warmup_parameters, warmup_solution = _solve_warmup( + problem, method, parameter_values; show_progress=show_progress + ) + result_full = HC.solve( + problem.system, + HC.solutions(warmup_solution); + start_parameters=warmup_parameters, + target_parameters=parameter_values, + show_progress=show_progress, + alg_default_options(method)..., + ) - if method == :warmup && !isempty(sweep) - warmup_parameters, warmup_solution = _solve_warmup( - problem, params_1D, sweep; threading=threading, show_progress=show_progress - ) - result_full = HC.solve( - problem.system, - HC.solutions(warmup_solution); - start_parameters=warmup_parameters, - target_parameters=parameter_values, - threading=threading, - show_progress=show_progress, - seed=seed, - kwargs..., + return reshape(result_full, size(parameter_values)...) +end + +"Uses HomotopyContinuation to solve `problem` at specified `parameter_values`." +function _get_raw_solution( + problem::Problem, method::Union{TotalDegree,Polyhedral}, parameter_values; show_progress +) + result_full = Array{Vector{Any},1}(undef, length(parameter_values)) + if show_progress + bar = Progress( + length(parameter_values); + dt=1, + desc="Solving via $method homotopy ...", + barlen=50, ) - elseif method == :total_degree || method == :polyhedral - result_full = Array{Vector{Any},1}(undef, length(parameter_values)) - if show_progress - bar = Progress( - length(parameter_values); - dt=1, - desc="Solving via $method homotopy ...", - barlen=50, - ) - end - for i in eachindex(parameter_values) # do NOT thread this - p = parameter_values[i] - show_progress ? ProgressMeter.next!(bar) : nothing - result_full[i] = [ - HC.solve( - problem.system; - start_system=method, - target_parameters=p, - threading=threading, - show_progress=false, - seed=seed, - ), - p, - ] - end - else - error("Unknown method: ", string(method)) + end + for i in eachindex(parameter_values) # do NOT thread this + p = parameter_values[i] + show_progress ? ProgressMeter.next!(bar) : nothing + result_full[i] = [ + HC.solve( + problem.system; + start_system=method_symbol(method), + target_parameters=p, + show_progress=false, + alg_default_options(method)..., + alg_specific_options(method)..., + ), + p, + ] end # reshape back to the original shape diff --git a/src/types.jl b/src/types.jl index 3defa00f..645a22ca 100644 --- a/src/types.jl +++ b/src/types.jl @@ -4,282 +4,6 @@ const StateDict = OrderedDict{Num,ComplexF64}; const SteadyState = Vector{ComplexF64}; const ParameterVector = Vector{Float64}; -""" -$(TYPEDEF) - -Holds differential equation(s) of motion and a set of harmonics to expand each variable. -This is the primary input for `HarmonicBalance.jl` ; after inputting the equations, the harmonics - ansatz needs to be specified using `add_harmonic!`. - -# Fields -$(TYPEDFIELDS) - -## Example -```julia-repl -julia> @variables t, x(t), y(t), ω0, ω, F, k; - -# equivalent ways to enter the simple harmonic oscillator -julia> DifferentialEquation(d(x,t,2) + ω0^2 * x - F * cos(ω*t), x); -julia> DifferentialEquation(d(x,t,2) + ω0^2 * x ~ F * cos(ω*t), x); - -# two coupled oscillators, one of them driven -julia> DifferentialEquation([d(x,t,2) + ω0^2 * x - k*y, d(y,t,2) + ω0^2 * y - k*x] .~ [F * cos(ω*t), 0], [x,y]); -``` -""" -mutable struct DifferentialEquation - """Assigns to each variable an equation of motion.""" - equations::OrderedDict{Num,Equation} - """Assigns to each variable a set of harmonics.""" - harmonics::OrderedDict{Num,OrderedSet{Num}} - - function DifferentialEquation(eqs) - return new(eqs, OrderedDict(var => OrderedSet() for var in keys(eqs))) - end - - # uses the above constructor if no harmonics defined - function DifferentialEquation(eqs::Vector{Equation}, vars::Vector{Num}) - return DifferentialEquation(OrderedDict(zip(vars, eqs))) - end - - # if expressions are entered instead of equations, automatically set them = 0 - function DifferentialEquation(exprs::Vector{Num}, vars::Vector{Num}) - return DifferentialEquation(exprs .~ Int(0), vars) - end - - function DifferentialEquation(eq::Equation, var::Num) - typerhs = typeof(eq.rhs) - typelhs = typeof(eq.lhs) - if eq.rhs isa AbstractVector || eq.lhs isa AbstractVector - throw( - ArgumentError( - "The equation is of the form $(typerhs)~$(typelhs) is not supported. Commenly one forgot to broadcast the equation symbol `~`.", - ), - ) - end - return DifferentialEquation([eq], [var]) - end - function DifferentialEquation(eq::Equation, vars::Vector{Num}) - typerhs = typeof(eq.rhs) - typelhs = typeof(eq.lhs) - throw( - ArgumentError( - "The variables are of type $(typeof(vars)). Whereas you gave one equation of type $(typerhs)~$(typelhs). Commenly one forgot to broadcast the equation symbol `~`.", - ), - ) - end - DifferentialEquation(lhs::Num, var::Num) = DifferentialEquation([lhs ~ Int(0)], [var]) -end - -function Base.show(io::IO, diff_eq::DifferentialEquation) - println(io, "System of ", length(keys(diff_eq.equations)), " differential equations") - println(io, "Variables: ", join(keys(diff_eq.equations), ", ")) - print(io, "Harmonic ansatz: ") - for var in keys(diff_eq.harmonics) - print(io, string(var), " => ", join(string.(diff_eq.harmonics[var]), ", ")) - print(io, "; ") - end - println(io, "\n") - return [println(io, eq) for eq in values(diff_eq.equations)] -end - -""" -$(TYPEDEF) - -Holds a variable stored under `symbol` describing the harmonic `ω` of `natural_variable`. - -# Fields -$(TYPEDFIELDS) -""" -mutable struct HarmonicVariable - """Symbol of the variable in the HarmonicBalance namespace.""" - symbol::Num - """Human-readable labels of the variable, used for plotting.""" - name::String - """Type of the variable (u or v for quadratures, a for a constant, Hopf for Hopf etc.)""" - type::String - """The harmonic being described.""" - ω::Num - """The natural variable whose harmonic is being described.""" - natural_variable::Num -end - -function Base.show(io::IO, hv::HarmonicVariable) - return println( - io, - "Harmonic variable ", - string.(hv.symbol) * " for harmonic ", - string(hv.ω), - " of ", - string(hv.natural_variable), - ) -end - -""" -$(TYPEDEF) - -Holds a set of algebraic equations governing the harmonics of a `DifferentialEquation`. - -# Fields -$(TYPEDFIELDS) -""" -mutable struct HarmonicEquation - """A set of equations governing the harmonics.""" - equations::Vector{Equation} - """A set of variables describing the harmonics.""" - variables::Vector{HarmonicVariable} - """The parameters of the equation set.""" - parameters::Vector{Num} - "The natural equation (before the harmonic ansatz was used)." - natural_equation::DifferentialEquation - - # use a self-referential constructor with _parameters - function HarmonicEquation(equations, variables, nat_eq) - return (x = new(equations, variables, Vector{Num}([]), nat_eq); - x.parameters = _parameters(x); - x) - end - function HarmonicEquation(equations, variables, parameters, natural_equation) - return new(equations, variables, parameters, natural_equation) - end -end - -function Base.show(io::IO, eom::HarmonicEquation) - println(io, "A set of ", length(eom.equations), " harmonic equations") - println(io, "Variables: ", join(string.(get_variables(eom)), ", ")) - println(io, "Parameters: ", join(string.(eom.parameters), ", ")) - println(io, "\nHarmonic ansatz: ", _show_ansatz(eom)) - println(io, "\nHarmonic equations:") - return [println(io, "\n", eq) for eq in eom.equations] -end - -"""Gives the relation between `var` and the underlying natural variable.""" -function _show_ansatz(var::HarmonicVariable) - t = var.natural_variable.val.arguments - t = length(t) == 1 ? string(t[1]) : error("more than 1 independent variable") - ω = string(var.ω) - terms = Dict("u" => "*cos(" * ω * t * ")", "v" => "*sin(" * ω * t * ")", "a" => "") - return string(string(var.symbol) * terms[var.type]) -end - -"""Gives the full harmonic ansatz used to construct `eom`.""" -function _show_ansatz(eom::HarmonicEquation) - output = "" - vars = unique(getfield.(eom.variables, :natural_variable)) - for nat_var in vars - # the Hopf variable (limit cycle frequency) does not contribute a term - harm_vars = filter( - x -> isequal(nat_var, x.natural_variable) && x.type !== "Hopf", eom.variables - ) - ansatz = join([_show_ansatz(var) for var in harm_vars], " + ") - output *= "\n" * string(nat_var) * " = " * ansatz - end - return output -end - -""" -$(TYPEDEF) - -Holds a set of algebraic equations describing the steady state of a system. - -# Fields -$(TYPEDFIELDS) - -# Constructors -```julia -Problem(eom::HarmonicEquation; Jacobian=true) # find and store the symbolic Jacobian -Problem(eom::HarmonicEquation; Jacobian="implicit") # ignore the Jacobian for now, compute implicitly later -Problem(eom::HarmonicEquation; Jacobian=J) # use J as the Jacobian (a function that takes a Dict) -Problem(eom::HarmonicEquation; Jacobian=false) # ignore the Jacobian -``` -""" -mutable struct Problem - "The harmonic variables to be solved for." - variables::Vector{Num} - "All symbols which are not the harmonic variables." - parameters::Vector{Num} - "The input object for HomotopyContinuation.jl solver methods." - system::HC.System - "The Jacobian matrix (possibly symbolic). - If `false`, the Jacobian is ignored (may be calculated implicitly after solving)." - jacobian - "The HarmonicEquation object used to generate this `Problem`." - eom::HarmonicEquation - - function Problem(variables, parameters, system, jacobian) - return new(variables, parameters, system, jacobian) - end #incomplete initialization for user-defined symbolic systems - function Problem(variables, parameters, system, jacobian, eom) - return new(variables, parameters, system, jacobian, eom) - end -end - -function Base.show(io::IO, p::Problem) - println(io, length(p.system.expressions), " algebraic equations for steady states") - println(io, "Variables: ", join(string.(p.variables), ", ")) - println(io, "Parameters: ", join(string.(p.parameters), ", ")) - return println(io, "Symbolic Jacobian: ", !(p.jacobian == false)) -end - -# assume this order of variables in all compiled function (transform_solutions, Jacobians) -function _free_symbols(p::Problem, varied) - return cat(p.variables, collect(keys(OrderedDict(varied))); dims=1) -end - -""" -$(TYPEDEF) - -Stores the steady states of a HarmonicEquation. - -# Fields -$(TYPEDFIELDS) - -""" -mutable struct Result - "The variable values of steady-state solutions." - solutions::Array{Vector{SteadyState}} - "Values of all parameters for all solutions." - swept_parameters::ParameterRange - "The parameters fixed throughout the solutions." - fixed_parameters::ParameterList - "The `Problem` used to generate this." - problem::Problem - "Maps strings such as \"stable\", \"physical\" etc to arrays of values, classifying the solutions (see method `classify_solutions!`)." - classes::Dict{String,Array} - "The Jacobian with `fixed_parameters` already substituted. Accepts a dictionary specifying the solution. - If problem.jacobian is a symbolic matrix, this holds a compiled function. - If problem.jacobian was `false`, this holds a function that rearranges the equations to find J - only after numerical values are inserted (preferable in cases where the symbolic J would be very large)." - jacobian::Function - "Seed used for the solver" - seed::Union{Nothing,UInt32} - - function Result(sol, swept, fixed, problem, classes, J, seed) - return new(sol, swept, fixed, problem, classes, J, seed) - end - Result(sol, swept, fixed, problem, classes) = new(sol, swept, fixed, problem, classes) - Result(sol, swept, fixed, problem) = new(sol, swept, fixed, problem, Dict([])) -end - -function Base.show(io::IO, r::Result) - println(io, "A steady state result for ", length(r.solutions), " parameter points") - println(io, "\nSolution branches: ", length(r.solutions[1])) - println(io, " of which real: ", sum(any.(classify_branch(r, "physical")))) - println(io, " of which stable: ", sum(any.(classify_branch(r, "stable")))) - return println(io, "\nClasses: ", join(keys(r.classes), ", ")) -end - -# assume this order of variables in all compiled function (transform_solutions, Jacobians) -function _free_symbols(res::Result) - return cat(res.problem.variables, collect(keys(res.swept_parameters)); dims=1) -end - -# overload to use [] for indexing -Base.getindex(r::Result, idx::Int...) = get_single_solution(r, idx) -Base.size(r::Result) = size(r.solutions) - -branch_count(r::Result) = length(r.solutions[1]) -get_branch(r::Result, idx) = getindex.(r.solutions, idx) - """ Represents a sweep of one or more parameters of a `HarmonicEquation`. diff --git a/test/API.jl b/test/API.jl index 02696e1e..f380b36c 100644 --- a/test/API.jl +++ b/test/API.jl @@ -14,18 +14,25 @@ using HarmonicBalance @test_throws ArgumentError DifferentialEquation(rhs ~ [F * cos(ω * t), 0], [x, y]) end -@testset begin - @variables ω1, ω2, ωₘ, t, ω, F, γ, λ, x(t), y(t) - eqs = [d(x, t, 2) + (ω1^2 - λ * cos(ωₘ * t)) * x + γ * d(x, t)] +@testset "get_steady_states API" begin + @variables ω1, t, ω, F, γ, λ, x(t), y(t) + eqs = [d(x, t, 2) + (ω1^2 - λ * cos(2 * ω * t)) * x + γ * d(x, t)] diff_eq = DifferentialEquation(eqs, [x]) add_harmonic!(diff_eq, x, ω) # drive frequency, close to ω1 harmonic_eq = get_harmonic_equations(diff_eq) + prob = HarmonicBalance.Problem(harmonic_eq) + varied = ω => range(0.7, 1.3, 100) - @test_throws MethodError get_steady_states(harmonic_eq, varied, threading=true) - @test_throws ArgumentError get_steady_states(harmonic_eq, Dict(varied), threading=true) + @test_throws MethodError get_steady_states(harmonic_eq, varied) + @test_throws ArgumentError get_steady_states(harmonic_eq, Dict(varied)) + + fixed = Dict(ω1 => 1.0, γ => 0.005, λ => 0.1) + @test_throws MethodError get_steady_states(prob, Dict()) + @test_throws MethodError get_steady_states(prob, varied, fixed) + r = get_steady_states(prob, HarmonicBalance.WarmUp(), varied, fixed) end @testset "forgot variable" begin diff --git a/test/ModelingToolkitExt.jl b/test/ModelingToolkitExt.jl index 690bc7b3..8f85bc6e 100644 --- a/test/ModelingToolkitExt.jl +++ b/test/ModelingToolkitExt.jl @@ -17,7 +17,7 @@ end ) fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) - param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) + param = HarmonicBalance.ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) sys = ODESystem(diff_eq) for p in string.([α, ω, ω0, F, γ]) @@ -33,8 +33,9 @@ end d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x ) + sys = ODESystem(harmonic_eq) fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) - param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) + param = HarmonicBalance.ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) ODEProblem(diff_eq, [1.0, 0.0], (0, 100), param) end @@ -50,7 +51,7 @@ end harmonic_eq = get_harmonic_equations(diff_eq) fixed = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01) - param = ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) + param = HarmonicBalance.ParameterList(merge(Dict(fixed), Dict(ω => 1.1))) @testset "ODESystem" begin sys = ODESystem(harmonic_eq) diff --git a/test/hysteresis_sweep.jl b/test/hysteresis_sweep.jl index b6a8bb58..e49117f4 100644 --- a/test/hysteresis_sweep.jl +++ b/test/hysteresis_sweep.jl @@ -11,7 +11,8 @@ harmonic_eq = get_harmonic_equations(diff_eq) # implement ansatz to get harmonic fixed = (α => 1, ω0 => 1.0, γ => 0.005, F => 0.005, η => 0.2) # fixed parameters varied = ω => range(0.95, 1.1, 10) # range of parameter values -result = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) +method = HarmonicBalance.WarmUp(; seed=SEED) +result = get_steady_states(harmonic_eq, method, varied, fixed; show_progress=false) followed_branches, _ = follow_branch(1, result) followed_branches, _ = follow_branch(1, result; y="√(u1^2+v1^2)") diff --git a/test/krylov.jl b/test/krylov.jl index 4f1cae94..69813368 100644 --- a/test/krylov.jl +++ b/test/krylov.jl @@ -27,5 +27,5 @@ end fixed = (ω0 => 1.0, γ => 0.005, α => 1.0, η => 0, F => 0.0, ψ => 0.0, θ => 0.0) varied = (ω => range(0.99, 1.01, 5), λ => range(1e-6, 0.05, 5)) -res1 = get_steady_states(harmonic_eq1, varied, fixed; show_progress=false, seed=SEED); -res2 = get_steady_states(harmonic_eq2, varied, fixed; show_progress=false, seed=SEED); +res1 = get_steady_states(harmonic_eq1, varied, fixed; show_progress=false); +res2 = get_steady_states(harmonic_eq2, varied, fixed; show_progress=false); diff --git a/test/limit_cycle.jl b/test/limit_cycle.jl index 8d498db7..60950b68 100644 --- a/test/limit_cycle.jl +++ b/test/limit_cycle.jl @@ -15,10 +15,8 @@ import HarmonicBalance.LinearResponse.plot_linear_response fixed = () varied = μ => range(2, 3, 2) - - result = get_limit_cycles( - harmonic_eq, varied, fixed, ω_lc; show_progress=false, seed=SEED - ) + method = HarmonicBalance.WarmUp(; seed=SEED) + result = get_limit_cycles(harmonic_eq, method, varied, fixed, ω_lc; show_progress=false) @test sum(any.(classify_branch(result, "stable"))) == 4 @test sum(any.(classify_branch(result, "unique_cycle"))) == 1 @@ -51,5 +49,5 @@ end # varied = (ω => range(0.992, 0.995, 2)) # # results -# result = get_limit_cycles(harmonic_eq, varied, fixed, ω_lc; threading=true) +# result = get_limit_cycles(harmonic_eq, varied, fixed, ω_lc) # end diff --git a/test/linear_response.jl b/test/linear_response.jl index a862cb91..605cff62 100644 --- a/test/linear_response.jl +++ b/test/linear_response.jl @@ -12,7 +12,7 @@ harmonic_eq = get_harmonic_equations(diff_eq) fixed = (α => 1, ω0 => 1.0, γ => 1e-2, F => 1e-6) varied = ω => range(0.9, 1.1, 10) -result = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) +result = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) @testset "first order" begin plot_linear_response( @@ -51,9 +51,7 @@ end fixed = (γ => 0.008, ω0 => 1.0, α => 1.0, F => 0.0, ωₚ => 1.0, λ => 0.016) varied = (ω => range(0.995, 1.005, 20)) - result_ω = get_steady_states( - harmonic_eq, varied, fixed; show_progress=false, seed=SEED - ) + result_ω = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) @test_throws ErrorException plot_eigenvalues(result_ω; branch=1) end end diff --git a/test/methods.jl b/test/methods.jl new file mode 100644 index 00000000..f23053bf --- /dev/null +++ b/test/methods.jl @@ -0,0 +1,43 @@ +using HarmonicBalance +using Test + +@testset "WarmUp" begin + @variables α ω ω0 λ γ t x(t) + diff_eq = DifferentialEquation( + d(d(x, t), t) + γ * d(x, t) + ω0^2 * (1 - λ * cos(2 * ω * t)) * x + α * x^3, x + ) + + add_harmonic!(diff_eq, x, ω) # + harmonic_eq = get_harmonic_equations(diff_eq) + + fixed = HarmonicBalance.ParameterList((α => 1.0, ω0 => 1.1, λ => 0.01, γ => 0.01)) + varied = HarmonicBalance.ParameterRange((ω => range(0.9, 1.1, 20),)) + prob = HarmonicBalance.Problem(harmonic_eq) + + unique_fixed, input_array = HarmonicBalance._prepare_input_params(prob, varied, fixed) + @test length.(input_array) == fill(5, 20) + + method = WarmUp() + warmup_parameters, warmup_solution = HarmonicBalance._solve_warmup( + prob, method, input_array; show_progress=false + ) + diff = abs.((input_array[method.index] .- warmup_parameters)) + @test diff ≈ zeros(5) atol = real(method.perturbation_size) * 5 +end + +@testset "Polyhedral" begin + @variables α ω ω0 λ γ t x(t) + diff_eq = DifferentialEquation( + d(d(x, t), t) + γ * d(x, t) + ω0^2 * (1 - λ * cos(2 * ω * t)) * x + α * x^3, x + ) + + add_harmonic!(diff_eq, x, ω) # + harmonic_eq = get_harmonic_equations(diff_eq) + + fixed = HarmonicBalance.ParameterList((α => 1.0, ω0 => 1.1, λ => 0.001, γ => 0.01)) + varied = HarmonicBalance.ParameterRange((ω => range(0.9, 1.1, 10),)) + result = get_steady_states( + harmonic_eq, Polyhedral(; only_non_zero=true), varied, fixed; show_progress=false + ) + @test sum(reduce(vcat, classify_branch(result, "stable"))) == 0 +end diff --git a/test/parametron.jl b/test/parametron.jl index 7c5483a0..6167326c 100644 --- a/test/parametron.jl +++ b/test/parametron.jl @@ -16,6 +16,8 @@ dEOM = DifferentialEquation(natural_equation + forces, x) add_harmonic!(dEOM, x, ω) harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); +method = HarmonicBalance.WarmUp(; seed=SEED) + @testset "undriven parametron" begin fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 0, α => 1.0, η => 0.3, θ => 0, ψ => 0) varied = ω => range(0.9, 1.1, 20) @@ -30,11 +32,11 @@ harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); @test length(prob.variables) == 2 @test length(prob.parameters) == 9 @test length(prob.system.parameters) == 9 - res = get_steady_states(prob, varied, fixed; show_progress=false, seed=SEED) + res = get_steady_states(prob, method, varied, fixed; show_progress=false) end @testset "steady states" begin - res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) + res = get_steady_states(harmonic_eq, method, varied, fixed; show_progress=false) @test length(res.solutions) == 20 @test length(res.solutions[1]) == 5 @test sum(any.(classify_branch(res, "physical"))) == 5 @@ -53,13 +55,13 @@ harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); @testset "implicit jacobian" begin p = HarmonicBalance.Problem(harmonic_eq; Jacobian="implicit") - res = get_steady_states(p, varied, fixed; show_progress=false, seed=SEED) + res = get_steady_states(p, method, varied, fixed; show_progress=false) end @testset "2d sweep" begin # try to run a 2D calculation fixed = (Ω => 1.0, γ => 1e-2, F => 0, α => 1.0, η => 0.1, θ => 0, ψ => 0) varied = (ω => range(0.9, 1.1, 20), λ => range(0.01, 0.05, 20)) - res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) + res = get_steady_states(harmonic_eq, method, varied, fixed; show_progress=false) end end diff --git a/test/plotting.jl b/test/plotting.jl index 86135686..5388c098 100644 --- a/test/plotting.jl +++ b/test/plotting.jl @@ -17,7 +17,7 @@ harmonic_eq = get_harmonic_equations(dEOM); fixed = (ω0 => 1.0, γ => 1e-2, λ => 5e-2, α => 1.0, η => 0.3) varied = ω => range(0.9, 1.1, 10) -res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) +res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) # plot 1D result plot(res; x="ω", y="u1"); @@ -27,7 +27,7 @@ plot_phase_diagram(res); fixed = (ω0 => 1.0, γ => 1e-2, α => 1.0, η => 0.3) varied = (ω => range(0.9, 1.1, 5), λ => range(0.01, 0.05, 5)) -res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED) +res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false) # plot 2D result plot_phase_diagram(res); diff --git a/test/runtests.jl b/test/runtests.jl index 1081be93..2a835ac9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,6 +61,7 @@ end @testset "Computing steady states" begin include("parametron.jl") include("krylov.jl") + include("methods.jl") end @testset "Processing solutions" begin diff --git a/test/time_evolution.jl b/test/time_evolution.jl index 1369e3bd..125f9a70 100644 --- a/test/time_evolution.jl +++ b/test/time_evolution.jl @@ -15,11 +15,8 @@ forces = F * cos(ω * t + θ) dEOM = DifferentialEquation(natural_equation + forces, x) add_harmonic!(dEOM, x, ω) harmonic_eq = get_harmonic_equations(dEOM; slow_time=T, fast_time=t); -p = HarmonicBalance.Problem(harmonic_eq); fixed = (Ω => 1.0, γ => 1e-2, λ => 5e-2, F => 1e-3, α => 1.0, η => 0.3, θ => 0, ψ => 0) -varied = ω => range(0.9, 1.1, 10) -res = get_steady_states(p, varied, fixed; show_progress=false, seed=SEED) tspan = (0.0, 10) sweep = AdiabaticSweep(ω => (0.9, 1.1), tspan) # linearly interpolate between two values at two times diff --git a/test/transform_solutions.jl b/test/transform_solutions.jl index 0411d61b..bf68be3b 100644 --- a/test/transform_solutions.jl +++ b/test/transform_solutions.jl @@ -19,7 +19,7 @@ harmonic_eq = get_harmonic_equations(dEOM); fixed = (Ω => 1.0, γ => 1e-2, F => 1e-3, α => 1.0) varied = ω => range(0.9, 1.1, 10) -res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false, seed=SEED); +res = get_steady_states(harmonic_eq, varied, fixed; show_progress=false); transform_solutions(res, "u1^2+v1^2") transform_solutions(res, "√(u1^2+v1^2)"; realify=true)