From 41416e6bae4b61440eaecb52fd52e74712993b2d Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Mon, 10 Jun 2024 11:16:02 +0200 Subject: [PATCH] Fix Limit cycle bug (#186) * Fix symbolic version * increase timout-minutes of CI * format limit_cycle module * get rid of unsued kwarg `explicit_Jacobian` in `_cycle_Problem` * change `cycle_harmonic` kwarg to `limit_cycle_harmonic` * instead of sorting always take the first in the list * clean test file --- src/modules/LimitCycles/analysis.jl | 4 +- src/modules/LimitCycles/gauge_fixing.jl | 68 +++++++++++++++---------- test/limit_cycle.jl | 29 ++++++----- 3 files changed, 58 insertions(+), 43 deletions(-) diff --git a/src/modules/LimitCycles/analysis.jl b/src/modules/LimitCycles/analysis.jl index a2bc049e..b6a7423b 100644 --- a/src/modules/LimitCycles/analysis.jl +++ b/src/modules/LimitCycles/analysis.jl @@ -1,6 +1,6 @@ import HarmonicBalance: classify_solutions, _free_symbols, _symidx, _is_physical -function classify_unique!(res::Result, Δω; class_name="unique_cycle") +function classify_unique!(res::Result, Δω; class_name = "unique_cycle") # 1st degeneracy: arbitrary sign of Δω i1 = _symidx(Δω, res) @@ -12,4 +12,4 @@ function classify_unique!(res::Result, Δω; class_name="unique_cycle") c2 = classify_solutions(res, soln -> _is_physical(soln) && real(soln[i2]) >= 0) res.classes[class_name] = map(.*, c1, c2) -end \ No newline at end of file +end diff --git a/src/modules/LimitCycles/gauge_fixing.jl b/src/modules/LimitCycles/gauge_fixing.jl index 5212dba7..eea32644 100644 --- a/src/modules/LimitCycles/gauge_fixing.jl +++ b/src/modules/LimitCycles/gauge_fixing.jl @@ -3,8 +3,8 @@ export add_pairs! using HarmonicBalance: is_rearranged, rearrange_standard, _remove_brackets using HarmonicBalance.LinearResponse: get_implicit_Jacobian, get_Jacobian -import HarmonicBalance: is_stable, is_physical, is_Hopf_unstable, order_branches!, classify_binaries!, find_branch_order - +import HarmonicBalance: is_stable, is_physical, is_Hopf_unstable, order_branches!, + classify_binaries!, find_branch_order function add_pairs!(eom::DifferentialEquation, ω_lc::Num) for var in get_variables(eom), ω in eom.harmonics[var] @@ -13,27 +13,24 @@ function add_pairs!(eom::DifferentialEquation, ω_lc::Num) end end - """ add_pairs!(eom::DifferentialEquation; ω_lc::Num, n=1) Add a limit cycle harmonic `ω_lc` to the system Equivalent to adding `n` pairs of harmonics ω +- ω_lc for each existing ω. """ -add_pairs!(eom::DifferentialEquation; ω_lc::Num, n::Int) = [add_pairs!(eom, ω_lc) for k in 1:n] - +add_pairs!(eom::DifferentialEquation; ω_lc::Num, n::Int) = [add_pairs!(eom, ω_lc) + for k in 1:n] """ $(TYPEDSIGNATURES) Return the harmonic variables which participate in the limit cycle labelled by `ω_lc`. """ -function get_cycle_variables(eom::HarmonicEquation, ω_lc::Num; sorted=false) +function get_cycle_variables(eom::HarmonicEquation, ω_lc::Num) vars = filter(x -> any(isequal.(ω_lc, get_all_terms(x.ω))), eom.variables) - sorted ? sort(vars, by = x -> x.ω / ω_lc) : vars end - """ Obtain the Jacobian of `eom` with a gauge-fixed variable `fixed_var`. `fixed_var` marks the variable which is fixed to zero due to U(1) symmetry. @@ -42,17 +39,17 @@ end For limit cycles, we always use an 'implicit' Jacobian - a function which only returns the numerical Jacobian when a numerical solution is inserted. Finding the analytical Jacobian is usually extremely time-consuming. """ -function _gaugefixed_Jacobian(eom::HarmonicEquation, fixed_var::HarmonicVariable; explicit=false, sym_order, rules) +function _gaugefixed_Jacobian(eom::HarmonicEquation, fixed_var::HarmonicVariable; + explicit = false, sym_order, rules) if explicit _fix_gauge!(get_Jacobian(eom), fixed_var) # get a symbolic explicit J, compile later else rules = Dict(rules) setindex!(rules, 0, _remove_brackets(fixed_var)) - get_implicit_Jacobian(eom, rules=rules, sym_order=sym_order) + get_implicit_Jacobian(eom, rules = rules, sym_order = sym_order) end end - """ $(TYPEDSIGNATURES) @@ -60,11 +57,11 @@ Construct a `Problem` from `eom` in the case where U(1) symmetry is present due to having added a limit cycle frequency `ω_lc`. `explicit_Jacobian=true` attempts to derive a symbolic Jacobian (usually not possible). """ -function _cycle_Problem(eom::HarmonicEquation, ω_lc::Num; explicit_Jacobian=false) - +function _cycle_Problem(eom::HarmonicEquation, ω_lc::Num) eom = deepcopy(eom) # do not mutate eom isempty(get_cycle_variables(eom, ω_lc)) ? error("No Hopf variables found!") : nothing - !any(isequal.(eom.parameters, ω_lc)) ? error(ω_lc, " is not a parameter of the harmonic equation!") : nothing + !any(isequal.(eom.parameters, ω_lc)) ? + error(ω_lc, " is not a parameter of the harmonic equation!") : nothing # eliminate one of the u,v variables (gauge fixing) fixed_var = _choose_fixed(eom, ω_lc) # remove the HV of the lowest subharmonic @@ -74,11 +71,14 @@ function _cycle_Problem(eom::HarmonicEquation, ω_lc::Num; explicit_Jacobian=fal _fix_gauge!(eom, ω_lc, fixed_var) # define Problem as usual but with the Hopf Jacobian (always computed implicitly) - p = Problem(eom; Jacobian=nothing) + p = Problem(eom; Jacobian = nothing) return p end -_choose_fixed(eom, ω_lc) = get_cycle_variables(eom, ω_lc, sorted=true)[1] +function _choose_fixed(eom, ω_lc) + vars = get_cycle_variables(eom, ω_lc) + first(vars) # This is arbitrary; better would be to substitute with values +end """ get_limit_cycles(eom::HarmonicEquation, swept, fixed, ω_lc; kwargs...) @@ -87,44 +87,56 @@ Variant of `get_steady_states` for a limit cycle problem characterised by a Hopf 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 = _cycle_Problem(eom, ω_lc, explicit_Jacobian=explicit_Jacobian); - prob.jacobian = _gaugefixed_Jacobian(eom, _choose_fixed(eom, ω_lc), explicit=explicit_Jacobian, sym_order=_free_symbols(prob, swept), rules=fixed) - result = get_steady_states(prob, swept, fixed; method=:warmup, threading=true, classify_default=true, kwargs...) +function get_limit_cycles( + eom::HarmonicEquation, swept, fixed, ω_lc; explicit_Jacobian = false, kwargs...) + + prob = _cycle_Problem(eom, ω_lc) + prob.jacobian = _gaugefixed_Jacobian( + eom, _choose_fixed(eom, ω_lc), explicit = explicit_Jacobian, + sym_order = _free_symbols(prob, swept), rules = fixed) + + result = get_steady_states(prob, swept, fixed; method = :warmup, + threading = true, classify_default = true, kwargs...) _classify_limit_cycles!(result, ω_lc) result end -get_limit_cycles(eom::HarmonicEquation,swept,fixed; cycle_harmonic, kwargs...) = get_limit_cycles(eom, swept, fixed, cycle_harmonic; kwargs...) +function get_limit_cycles(eom::HarmonicEquation, swept, fixed; limit_cycle_harmonic, kwargs...) + get_limit_cycles(eom, swept, fixed, limit_cycle_harmonic; kwargs...) +end # if abs(ω_lc) < tol, set all classifications to false # TOLERANCE HARDCODED FOR NOW function _classify_limit_cycles!(res::Result, ω_lc::Num) ω_lc_idx = findfirst(x -> isequal(x, ω_lc), res.problem.variables) - for idx in CartesianIndices(res.solutions), c in filter(x -> x != "binary_labels", keys(res.classes)) - res.classes[c][idx] .*= abs.(getindex.(res.solutions[idx], ω_lc_idx)) .> 1E-10 + for idx in CartesianIndices(res.solutions), + c in filter(x -> x != "binary_labels", keys(res.classes)) + + res.classes[c][idx] .*= abs.(getindex.(res.solutions[idx], ω_lc_idx)) .> 1e-10 end classify_unique!(res, ω_lc) - unique_stable = find_branch_order(map(.*, res.classes["stable"], res.classes["unique_cycle"])) + unique_stable = find_branch_order(map( + .*, res.classes["stable"], res.classes["unique_cycle"])) # branches which are unique but never stable - unique_unstable = setdiff(find_branch_order(map(.*, res.classes["unique_cycle"], res.classes["physical"])), unique_stable) + unique_unstable = setdiff( + find_branch_order(map(.*, res.classes["unique_cycle"], res.classes["physical"])), + unique_stable) order_branches!(res, vcat(unique_stable, unique_unstable)) # shuffle to have relevant branches first end - """ $(TYPEDSIGNATURES) Fix the gauge in `eom` where `ω_lc` is the limit cycle frequency by constraining `fixed_var` to zero and promoting `ω_lc` to a variable. """ function _fix_gauge!(eom::HarmonicEquation, ω_lc::Num, fixed_var::HarmonicVariable) - - new_symbol = HarmonicBalance.declare_variable(string(ω_lc), first(get_independent_variables(eom))) + new_symbol = HarmonicBalance.declare_variable( + string(ω_lc), first(get_independent_variables(eom))) rules = Dict(ω_lc => new_symbol, fixed_var.symbol => Num(0)) eom.equations = expand_derivatives.(substitute_all(eom.equations, rules)) eom.parameters = setdiff(eom.parameters, [ω_lc]) # ω_lc is now NOT a parameter anymore diff --git a/test/limit_cycle.jl b/test/limit_cycle.jl index 0e54df0d..e0c0d548 100644 --- a/test/limit_cycle.jl +++ b/test/limit_cycle.jl @@ -2,24 +2,27 @@ using HarmonicBalance import HarmonicBalance.LinearResponse.plot_linear_response import HarmonicBalance.LimitCycles.get_limit_cycles -@variables ω_lc, t, ω0, x(t), μ -natural_equation = d(d(x, t), t) - μ * (1 - x^2) * d(x, t) + x -dEOM = DifferentialEquation(natural_equation, x) +@testset "van der Pol oscillator " begin + @variables ω_lc, t, ω0, x(t), μ -# order should NOT matter if correct gauge is fixed (that corresponding to 1*ω_lc) -add_harmonic!(dEOM, x, [ω_lc, 3 * ω_lc]) + natural_equation = d(d(x, t), t) - μ * (1 - x^2) * d(x, t) + x + dEOM = DifferentialEquation(natural_equation, x) -harmonic_eq = get_harmonic_equations(dEOM) + # order matters for 1*ω_lc gauge to be fixed + add_harmonic!(dEOM, x, [ω_lc, 3 * ω_lc]) -fixed = (); -varied = μ => range(1, 5, 5) + harmonic_eq = get_harmonic_equations(dEOM) + HarmonicBalance.LimitCycles._choose_fixed(harmonic_eq, ω_lc) -result = get_limit_cycles(harmonic_eq, varied, fixed, cycle_harmonic=ω_lc, show_progress=false, seed=SEED) + fixed = (); + varied = μ => range(1, 5, 5) -@test sum(any.(classify_branch(result, "stable"))) == 4 -@test sum(any.(classify_branch(result, "unique_cycle"))) == 1 + result = get_limit_cycles(harmonic_eq, varied, fixed, ω_lc; show_progress=false, seed=SEED) -plot(result, y="ω_lc") + @test sum(any.(classify_branch(result, "stable"))) == 4 + @test sum(any.(classify_branch(result, "unique_cycle"))) == 1 -plot_linear_response(result, x, branch=1, Ω_range=range(0.9, 1.1, 2), order=1) + plot(result, y="ω_lc") + plot_linear_response(result, x, branch=1, Ω_range=range(0.9, 1.1, 2), order=1) +end