-
Notifications
You must be signed in to change notification settings - Fork 8
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
8 changed files
with
225 additions
and
79 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -9,3 +9,4 @@ local_tests | |
| .DS_Store | ||
| docs/build/** | ||
| examples/jdp.mplstyle | ||
| .JuliaFormatter.toml | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,59 @@ | ||
| module SteadyStateDiffEqExt | ||
|
|
||
| export steady_state_sweep | ||
| import HarmonicBalance: steady_state_sweep | ||
| using SteadyStateDiffEq: solve, NonlinearProblem, SteadyStateProblem, DynamicSS, remake | ||
| using LinearAlgebra: norm, eigvals | ||
|
|
||
| function steady_state_sweep( | ||
| prob::SteadyStateProblem, alg::DynamicSS; | ||
| varied::Pair, kwargs...) | ||
| varied_idx, sweep_range = varied | ||
|
|
||
| # if p is dual number (AD) result is dual number | ||
| result = [similar(prob.p, 2) for _ in sweep_range] | ||
|
|
||
| foreach(pairs(sweep_range)) do (i, value) | ||
| u0 = i == 1 ? [0.0, 0.0] : result[i - 1] | ||
| # make type-stable: FD.Dual or Float64 | ||
| parameters = eltype(prob_np.p)[i == varied_idx ? value : x | ||
| for (i, x) in enumerate(prob_np.p)] | ||
| sol = solve(remake(prob, p = parameters, u0 = u0), alg; kwargs...) | ||
| result[i] = sol.u | ||
| end | ||
| return result | ||
| end | ||
|
|
||
| function steady_state_sweep( | ||
| prob_np::NonlinearProblem, prob_ss::SteadyStateProblem, | ||
| alg_np, alg_ss::DynamicSS; | ||
| varied::Pair, kwargs...) | ||
| varied_idx, sweep_range = varied | ||
| # if p is dual number (AD) result is dual number | ||
| result = [similar(prob_np.p) for _ in sweep_range] | ||
|
|
||
| foreach(pairs(sweep_range)) do (i, value) | ||
| u0 = i == 1 ? Base.zeros(length(prob_np.u0)) : result[i - 1] | ||
| # make type-stable: FD.Dual or Float64 | ||
| parameters = eltype(prob_np.p)[i == varied_idx ? value : x | ||
| for (i, x) in enumerate(prob_np.p)] | ||
| sol_nn = solve(remake(prob_np, p = parameters, u0 = u0), alg_np; kwargs...) | ||
|
|
||
| # last argument is time but does not matter | ||
| zeros = norm(prob_np.f.f.f.f.f_oop(sol_nn.u, parameters, 0)) | ||
| jac = prob_np.f.jac.f.f.f_oop(sol_nn.u, parameters, 0) | ||
| eigval = jac isa Vector ? jac : eigvals(jac) # eigvals favourable supports FD.Dual | ||
|
|
||
| if !isapprox(zeros, 0, atol = 1e-5) || any(λ -> λ > 0, real.(eigval)) | ||
| sol_ss = solve(remake(prob_ss, p = parameters, u0 = u0), | ||
| alg_ss, abstol = 1e-5, reltol = 1e-5) | ||
| result[i] = sol_ss.u | ||
| else | ||
| result[i] = sol_nn.u | ||
| end | ||
| end | ||
|
|
||
| return result | ||
| end | ||
|
|
||
| end # module |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,83 +1,85 @@ | ||
| module HarmonicBalance | ||
|
|
||
| using Printf | ||
| using OrderedCollections | ||
| using Symbolics | ||
| using ProgressMeter | ||
| using DocStringExtensions | ||
| using BijectiveHilbert | ||
| using LinearAlgebra | ||
| using Plots, Latexify | ||
| using Random | ||
| import HomotopyContinuation; | ||
| const HC = HomotopyContinuation | ||
| import Distances | ||
| # using SnoopPrecompile | ||
| using Printf | ||
| using OrderedCollections | ||
| using Symbolics | ||
| using ProgressMeter | ||
| using DocStringExtensions | ||
| using BijectiveHilbert | ||
| using LinearAlgebra | ||
| using Plots, Latexify | ||
| using Random | ||
| import HomotopyContinuation | ||
| const HC = HomotopyContinuation | ||
| import Distances | ||
| # using SnoopPrecompile | ||
|
|
||
| import Base: show, display; export show | ||
| export * | ||
| export @variables | ||
| export d | ||
| import Base: show, display | ||
| export show | ||
| export * | ||
| export @variables | ||
| export d | ||
|
|
||
| import Base: ComplexF64, Float64; export ComplexF64, Float64 | ||
| ComplexF64(x::Complex{Num}) = ComplexF64(Float64(x.re) + im*Float64(x.im)) | ||
| Float64(x::Complex{Num}) = Float64(ComplexF64(x)) | ||
| Float64(x::Num) = Float64(x.val) | ||
| import Base: ComplexF64, Float64 | ||
| export ComplexF64, Float64 | ||
| ComplexF64(x::Complex{Num}) = ComplexF64(Float64(x.re) + im * Float64(x.im)) | ||
| Float64(x::Complex{Num}) = Float64(ComplexF64(x)) | ||
| Float64(x::Num) = Float64(x.val) | ||
|
|
||
| # default global settings | ||
| export IM_TOL | ||
| IM_TOL::Float64 = 1E-6 | ||
| function set_imaginary_tolerance(x::Float64) | ||
| @eval(IM_TOL::Float64 = $x) | ||
| end | ||
| # default global settings | ||
| export IM_TOL | ||
| IM_TOL::Float64 = 1E-6 | ||
| function set_imaginary_tolerance(x::Float64) | ||
| @eval(IM_TOL::Float64=$x) | ||
| end | ||
|
|
||
| export is_real | ||
| is_real(x) = abs(imag(x)) / abs(real(x)) < IM_TOL::Float64 || abs(x) < 1e-70 | ||
| is_real(x::Array) = is_real.(x) | ||
| export is_real | ||
| is_real(x) = abs(imag(x)) / abs(real(x)) < IM_TOL::Float64 || abs(x) < 1e-70 | ||
| is_real(x::Array) = is_real.(x) | ||
|
|
||
| # Symbolics does not natively support complex exponentials of variables | ||
| import Base: exp | ||
| exp(x::Complex{Num}) = x.re.val == 0 ? exp(im*x.im.val) : exp(x.re.val + im*x.im.val) | ||
| # Symbolics does not natively support complex exponentials of variables | ||
| import Base: exp | ||
| exp(x::Complex{Num}) = x.re.val == 0 ? exp(im * x.im.val) : exp(x.re.val + im * x.im.val) | ||
|
|
||
| include("types.jl") | ||
| include("types.jl") | ||
|
|
||
| include("Symbolics_customised.jl") | ||
| include("Symbolics_utils.jl") | ||
| include("DifferentialEquation.jl") | ||
| include("HarmonicVariable.jl") | ||
| include("HarmonicEquation.jl") | ||
| include("solve_homotopy.jl") | ||
| include("sorting.jl") | ||
| include("classification.jl") | ||
| include("saving.jl") | ||
| include("transform_solutions.jl") | ||
| include("plotting_Plots.jl") | ||
| include("Symbolics_customised.jl") | ||
| include("Symbolics_utils.jl") | ||
| include("DifferentialEquation.jl") | ||
| include("HarmonicVariable.jl") | ||
| include("HarmonicEquation.jl") | ||
| include("solve_homotopy.jl") | ||
| include("sorting.jl") | ||
| include("classification.jl") | ||
| include("saving.jl") | ||
| include("transform_solutions.jl") | ||
| include("plotting_Plots.jl") | ||
|
|
||
| include("modules/HC_wrapper.jl") | ||
| using .HC_wrapper | ||
| include("modules/HC_wrapper.jl") | ||
| using .HC_wrapper | ||
|
|
||
| include("modules/LinearResponse.jl") | ||
| using .LinearResponse | ||
| export plot_linear_response, plot_rotframe_jacobian_response | ||
| include("modules/LinearResponse.jl") | ||
| using .LinearResponse | ||
| export plot_linear_response, plot_rotframe_jacobian_response | ||
|
|
||
| include("modules/TimeEvolution.jl") | ||
| if !isdefined(Base, :get_extension) | ||
| include("ext/TimeEvolution.jl") | ||
| using .TimeEvolution | ||
| end | ||
| export ParameterSweep, ODEProblem, solve, ODESystem | ||
| export plot_1D_solutions_branch, follow_branch | ||
| include("modules/TimeEvolution.jl") | ||
| if !isdefined(Base, :get_extension) | ||
| include("ext/TimeEvolution.jl") | ||
| using .TimeEvolution | ||
| end | ||
| export ParameterSweep, ODEProblem, solve, ODESystem, steady_state_sweep | ||
| export plot_1D_solutions_branch, follow_branch | ||
|
|
||
| include("modules/LimitCycles.jl") | ||
| using .LimitCycles | ||
| include("modules/LimitCycles.jl") | ||
| using .LimitCycles | ||
|
|
||
| include("modules/KrylovBogoliubov.jl") | ||
| using .KrylovBogoliubov | ||
| export first_order_transform!, is_rearranged_standard, rearrange_standard!, get_equations | ||
| export get_krylov_equations | ||
| include("modules/KrylovBogoliubov.jl") | ||
| using .KrylovBogoliubov | ||
| export first_order_transform!, is_rearranged_standard, rearrange_standard!, get_equations | ||
| export get_krylov_equations | ||
|
|
||
| # precomp_path = (@__DIR__) * "/../test/" | ||
| # @precompile_all_calls include(precomp_path * "parametron.jl") | ||
| # @precompile_all_calls include(precomp_path * "plotting.jl") | ||
| # precomp_path = (@__DIR__) * "/../test/" | ||
| # @precompile_all_calls include(precomp_path * "parametron.jl") | ||
| # @precompile_all_calls include(precomp_path * "plotting.jl") | ||
|
|
||
| end # module |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,3 +1,4 @@ | ||
| function follow_branch end | ||
| function plot_1D_solutions_branch end | ||
| function FFT end | ||
| function steady_state_sweep end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,60 @@ | ||
| using ModelingToolkit, SteadyStateDiffEq, OrdinaryDiffEq, LinearAlgebra, NonlinearSolve | ||
|
|
||
| @testset "Steady state sweeps" begin | ||
|
|
||
| @testset "one variable ODE" begin | ||
| @variables t v(t)=0 | ||
| @parameters g=9.8 k=0.2 | ||
| D = Differential(t) | ||
| eqs = [D(v) ~ g - k * v] | ||
| @named model = ODESystem(eqs) | ||
|
|
||
| model = structural_simplify(model) | ||
|
|
||
| prob_ss = SteadyStateProblem{false}(model, [], [], jac = true) | ||
| prob_np = NonlinearProblem(prob_ss) | ||
|
|
||
| varied = 2 => range(0, 1, 100) | ||
|
|
||
| swept = steady_state_sweep(prob_np, prob_ss, NewtonRaphson(), DynamicSS(Rodas5()); varied=varied) | ||
| end | ||
|
|
||
| @testset "two variable ODE (duffing)" begin | ||
| @variables t u1(t) v1(t); | ||
| @parameters α ω ω0 F γ | ||
|
|
||
| eqs = [ | ||
| Differential(t)(u1) ~ (F*γ - u1*γ*(ω^2) - u1*γ*(ω0^2) - v1*(γ^2)*ω - 2*v1*(ω^3) + 2*v1*ω*(ω0^2) - (3//4)*(u1^3)*α*γ + (3//2)*(u1^2)*v1*α*ω - (3//4)*u1*(v1^2)*α*γ + (3//2)*(v1^3)*α*ω) / (γ^2 + (4)*(ω^2)), | ||
| Differential(t)(v1) ~ (-2*F*ω - u1*(γ^2)*ω - (2)*u1*(ω^3) + 2*u1*ω*(ω0^2) + v1*γ*(ω^2) + v1*γ*(ω0^2) + (3//2)*(u1^3)*α*ω + (3//4)*(u1^2)*v1*α*γ + (3//2)*u1*(v1^2)*α*ω + (3//4)*(v1^3)*α*γ) / (-(γ^2) - (4)*(ω^2)) | ||
| ] | ||
|
|
||
| @named model = ODESystem(eqs, t, [u1, v1], [α, ω, ω0, F, γ]) | ||
| model = structural_simplify(model) | ||
|
|
||
| force = 0.01; omega0 = 1.1; alpha = 1.0; gamma = 0.01; | ||
| param = [alpha, omega0, omega0, force, gamma] | ||
| x0 = [1.0, 0.0]; | ||
| prob_ss = SteadyStateProblem{true}(model, x0, param, jac = true) | ||
| prob_np = NonlinearProblem(prob_ss) | ||
|
|
||
| ω_span = (0.9, 1.5); ω_range = range(ω_span..., 100) | ||
| varied = 2 => ω_range | ||
| swept = steady_state_sweep(prob_np, prob_ss, NewtonRaphson(), DynamicSS(Rodas5()); varied=varied) | ||
|
|
||
| @test length(swept) == 100 | ||
| @test length(swept[1]) == 2 | ||
| @test norm.(swept) isa Array{Float64, 1} | ||
| # using Plots, LinearAlgebra; plot(norm.(swept)) | ||
|
|
||
| function has_discontinuity(v::Vector{Float64}) | ||
| threshold = 1.0e-1 # Define a threshold for the discontinuity | ||
| for i in 2:length(v) | ||
| abs(v[i] - v[i-1]) > threshold && return true | ||
| end | ||
| return false | ||
| end | ||
|
|
||
| @test has_discontinuity(norm.(swept)) | ||
| end | ||
|
|
||
| end # Steady state sweeps |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters