From afbb0909fa07aba3bb50b5786cf6a7bd160179d8 Mon Sep 17 00:00:00 2001 From: oameye Date: Tue, 20 Feb 2024 11:29:11 +0100 Subject: [PATCH] add SteadyStateDiffEqExt --- .gitignore | 1 + Project.toml | 6 +- ext/ModelingToolkitExt.jl | 36 +++++++--- ext/SteadyStateDiffEqExt.jl | 59 +++++++++++++++ src/HarmonicBalance.jl | 134 ++++++++++++++++++----------------- src/modules/TimeEvolution.jl | 1 + test/SteadyStateDiffEqExt.jl | 60 ++++++++++++++++ test/runtests.jl | 7 +- 8 files changed, 225 insertions(+), 79 deletions(-) create mode 100644 ext/SteadyStateDiffEqExt.jl create mode 100644 test/SteadyStateDiffEqExt.jl diff --git a/.gitignore b/.gitignore index 039c7c40..03da54a6 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,4 @@ local_tests .DS_Store docs/build/** examples/jdp.mplstyle +.JuliaFormatter.toml \ No newline at end of file diff --git a/Project.toml b/Project.toml index f984252f..0f137b79 100644 --- a/Project.toml +++ b/Project.toml @@ -25,6 +25,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [weakdeps] OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" [compat] BijectiveHilbert = "0.3.0" @@ -47,12 +48,15 @@ julia = "1.8.2" [extensions] TimeEvolution = "OrdinaryDiffEq" ModelingToolkitExt = "ModelingToolkit" +SteadyStateDiffEqExt = "SteadyStateDiffEq" [extras] Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" [targets] -test = ["Pkg", "Test", "OrdinaryDiffEq", "ModelingToolkit"] +test = ["Pkg", "Test", "OrdinaryDiffEq", "ModelingToolkit", "SteadyStateDiffEq", "NonlinearSolve"] diff --git a/ext/ModelingToolkitExt.jl b/ext/ModelingToolkitExt.jl index 00c31fb5..2b10a04a 100644 --- a/ext/ModelingToolkitExt.jl +++ b/ext/ModelingToolkitExt.jl @@ -3,23 +3,22 @@ module ModelingToolkitExt export ODESystem, ODEProblem using HarmonicBalance: HarmonicEquation, is_rearranged, - rearrange_standard, get_variables, simplify, ParameterList + rearrange_standard, get_variables, simplify, ParameterList using ModelingToolkit -import ModelingToolkit: ODESystem, ODEProblem +import ModelingToolkit: ODESystem, ODEProblem, NonlinearProblem, SteadyStateProblem swapsides(eq::Equation) = Equation(eq.rhs, eq.lhs) function declare_parameter(var::Num) var_sym = Symbol(var) new_var = @parameters $var_sym - @eval($(var_sym) = first($new_var)) # store the variable under "name" in this namespace + @eval($(var_sym)=first($new_var)) # store the variable under "name" in this namespace return eval(var_sym) - end +end function ODESystem(eom::HarmonicEquation) - if !is_rearranged(eom) # check if time-derivatives of the variable are on the right hand side - eom = rearrange_standard(eom) + eom = rearrange_standard(eom) end slow_time = (@variables T; T) @@ -35,13 +34,32 @@ function ODESystem(eom::HarmonicEquation) return sys end -function ODEProblem(eom::HarmonicEquation, u0, tspan::Tuple, p::ParameterList; in_place=true, kwargs...) +function ODEProblem(eom::HarmonicEquation, u0, tspan::Tuple, + p::ParameterList; in_place = true, kwargs...) + sys = ODESystem(eom) + param = ModelingToolkit.varmap_to_vars(p, parameters(sys)) + if !in_place # out-of-place + prob = ODEProblem{false}(sys, u0, tspan, param; jac = true, kwargs...) + else # in-place + prob = ODEProblem{true}(sys, u0, tspan, param; jac = true, kwargs...) + end + return prob +end + +function NonlinearProblem( + eom::HarmonicEquation, u0, p::ParameterList; in_place = true, kwargs...) + ss_prob = SteadyStateProblem(eom, u0, p::ParameterList; in_place = in_place, kwargs...) + return NonlinearProblem(ss_prob) # gives warning of some internal deprication +end + +function SteadyStateProblem( + eom::HarmonicEquation, u0, p::ParameterList; in_place = true, kwargs...) sys = ODESystem(eom) param = ModelingToolkit.varmap_to_vars(p, parameters(sys)) if !in_place # out-of-place - prob = ODEProblem{false}(sys, u0, tspan, param; jac=true, kwargs...) + prob = SteadyStateProblem{false}(sys, u0, param; jac = true, kwargs...) else # in-place - prob = ODEProblem{true}(sys, u0, tspan, param; jac=true, kwargs...) + prob = SteadyStateProblem{true}(sys, u0, param; jac = true, kwargs...) end return prob end diff --git a/ext/SteadyStateDiffEqExt.jl b/ext/SteadyStateDiffEqExt.jl new file mode 100644 index 00000000..ab830729 --- /dev/null +++ b/ext/SteadyStateDiffEqExt.jl @@ -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 diff --git a/src/HarmonicBalance.jl b/src/HarmonicBalance.jl index 3984c970..cfc4f89d 100644 --- a/src/HarmonicBalance.jl +++ b/src/HarmonicBalance.jl @@ -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 diff --git a/src/modules/TimeEvolution.jl b/src/modules/TimeEvolution.jl index cdca9a25..ecf20c09 100644 --- a/src/modules/TimeEvolution.jl +++ b/src/modules/TimeEvolution.jl @@ -1,3 +1,4 @@ function follow_branch end function plot_1D_solutions_branch end function FFT end +function steady_state_sweep end diff --git a/test/SteadyStateDiffEqExt.jl b/test/SteadyStateDiffEqExt.jl new file mode 100644 index 00000000..b460f827 --- /dev/null +++ b/test/SteadyStateDiffEqExt.jl @@ -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 diff --git a/test/runtests.jl b/test/runtests.jl index e721b9ad..60bcc18f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,12 +16,13 @@ files = [ "parametron.jl", "transform_solutions.jl", "plotting.jl", - "time_evolution.jl", "krylov.jl", - "hysteresis_sweep.jl", "linear_response.jl", "limit_cycle.jl", - "ModelingToolkitExt.jl" + "time_evolution.jl", + "hysteresis_sweep.jl", + "ModelingToolkitExt.jl", + "SteadyStateDiffEqExt.jl" ] for file in files