diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml index 68fc360..901321a 100644 --- a/.github/workflows/Tests.yml +++ b/.github/workflows/Tests.yml @@ -24,7 +24,6 @@ jobs: matrix: version: - "1" - - "1.6" uses: "SciML/.github/.github/workflows/tests.yml@v1" with: julia-version: "${{ matrix.version }}" diff --git a/Project.toml b/Project.toml index 10c093f..9a967c1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelOrderReduction" uuid = "207801d6-6cee-43a9-ad0c-f0c64933efa0" -authors = ["Bowen S. Zhu and contributors"] -version = "0.1.1" +authors = ["Bowen S. Zhu and contributors"] +version = "0.1.2" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -10,18 +10,14 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" RandomizedLinAlg = "0448d7d9-159c-5637-8537-fd72090fea46" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" -Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e" [compat] DocStringExtensions = "0.8, 0.9" LinearAlgebra = "1" -ModelingToolkit = "8.21" +ModelingToolkit = "9" RandomizedLinAlg = "0.1" Setfield = "0.8, 1" SparseArrays = "1" -SymbolicUtils = "0.19" -Symbolics = "4.10" TSVD = "0.4" -julia = "1.6" +julia = "1.10" diff --git a/docs/Project.toml b/docs/Project.toml index 5e2e801..702b756 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -12,7 +12,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" DifferentialEquations = "7.6" Documenter = "1" LaTeXStrings = "1.3" -MethodOfLines = "0.6, 0.7" +MethodOfLines = "0.11" ModelOrderReduction = "0.1" -ModelingToolkit = "8.33" -Plots = "1.36" +ModelingToolkit = "9" +Plots = "1.40" diff --git a/docs/make.jl b/docs/make.jl index fc3171f..ecb4609 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -6,12 +6,12 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true) include("pages.jl") makedocs(sitename = "ModelOrderReduction.jl", - authors = "Bowen S. Zhu", - modules = [ModelOrderReduction], - clean = true, doctest = false, linkcheck = true, - warnonly = [:missing_docs, :example_block], - format = Documenter.HTML(assets = ["assets/favicon.ico"], - canonical = "https://docs.sciml.ai/ModelOrderReduction/stable/"), - pages = pages) + authors = "Bowen S. Zhu", + modules = [ModelOrderReduction], + clean = true, doctest = false, linkcheck = true, + warnonly = [:missing_docs, :example_block], + format = Documenter.HTML(assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/ModelOrderReduction/stable/"), + pages = pages) deploydocs(repo = "github.com/SciML/ModelOrderReduction.jl.git") diff --git a/docs/pages.jl b/docs/pages.jl index 2a7c950..8ddcb50 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,5 +1,5 @@ pages = [ "Home" => "index.md", "Functions" => "functions.md", - "Tutorials" => Any["tutorials/deim_FitzHugh-Nagumo.md"], + "Tutorials" => Any["tutorials/deim_FitzHugh-Nagumo.md"] ] diff --git a/docs/src/tutorials/deim_FitzHugh-Nagumo.md b/docs/src/tutorials/deim_FitzHugh-Nagumo.md index 41e17cc..e53fef7 100644 --- a/docs/src/tutorials/deim_FitzHugh-Nagumo.md +++ b/docs/src/tutorials/deim_FitzHugh-Nagumo.md @@ -107,7 +107,7 @@ svdval_w = svdvals(snapshot_w) svdval_fv = svdvals(snapshot_fv) using Plots, LaTeXStrings svd_plt = plot(yscale = :log10, xticks = eachindex(svdval_v), titlefont = 11, - legendfont = 10, title = "Singular values of the snapshots") + legendfont = 10, title = "Singular values of the snapshots") plot!(svd_plt, svdval_v, markershape = :circle, label = L"Singular Val of $v$") plot!(svd_plt, svdval_w, markershape = :circle, label = L"Singular Val of $w$") plot!(svd_plt, svdval_fv, markershape = :circle, label = L"Singular Val of $f(v)$") @@ -118,9 +118,9 @@ points ``x`` from the full-order system. ```@example deim_FitzHugh_Nagumo full_plt = plot(xlabel = L"v(x,t)", ylabel = L"x", zlabel = L"w(x,t)", xlims = (-0.5, 2.0), - ylims = (0.0, L), zlims = (0.0, 0.25), legend = false, xflip = true, - camera = (50, 30), titlefont = 10, - title = "Phase−Space diagram of full $(nameof(pde_sys)) system") + ylims = (0.0, L), zlims = (0.0, 0.25), legend = false, xflip = true, + camera = (50, 30), titlefont = 10, + title = "Phase−Space diagram of full $(nameof(pde_sys)) system") @views for i in 1:nₓ plot!(full_plt, snapshot_v[i, :], _ -> sol_x[i], snapshot_w[i, :]) end @@ -149,9 +149,9 @@ And plot the result from the POD-DEIM reduced system. ```@example deim_FitzHugh_Nagumo deim_plt = plot(xlabel = L"v(x,t)", ylabel = L"x", zlabel = L"w(x,t)", xlims = (-0.5, 2.0), - ylims = (0.0, L), zlims = (0.0, 0.25), legend = false, xflip = true, - camera = (50, 30), titlefont = 10, - title = "Phase−Space diagram of reduced $(nameof(pde_sys)) system") + ylims = (0.0, L), zlims = (0.0, 0.25), legend = false, xflip = true, + camera = (50, 30), titlefont = 10, + title = "Phase−Space diagram of reduced $(nameof(pde_sys)) system") @views for i in 1:nₓ plot!(deim_plt, sol_deim_v[i, :], _ -> sol_deim_x[i], sol_deim_w[i, :]) end @@ -178,12 +178,12 @@ function unconnected(v::AbstractVector, nₜ::Integer) vec(data) end plt_2 = plot(xlabel = L"v(x,t)", ylabel = L"x", zlabel = L"w(x,t)", xlims = (-0.5, 2.0), - ylims = (0.0, L), zlims = (0.0, 0.25), xflip = true, camera = (50, 30), - titlefont = 10, title = "Comparison of full and reduced systems") + ylims = (0.0, L), zlims = (0.0, 0.25), xflip = true, camera = (50, 30), + titlefont = 10, title = "Comparison of full and reduced systems") plot!(plt_2, unconnected(snapshot_v), unconnected(sol_x, nₜ), unconnected(snapshot_w), - label = "Full$(length(ode_sys.eqs))") + label = "Full$(length(ode_sys.eqs))") plot!(plt_2, unconnected(sol_deim_v), unconnected(sol_deim_x, nₜ_deim), - unconnected(sol_deim_w), label = "POD$(pod_dim)/DEIM$(deim_dim)") + unconnected(sol_deim_w), label = "POD$(pod_dim)/DEIM$(deim_dim)") ``` As we can see, the reduced-order system captures the limit cycle of the original full-order diff --git a/src/DataReduction/POD.jl b/src/DataReduction/POD.jl index 854b706..2abd01a 100644 --- a/src/DataReduction/POD.jl +++ b/src/DataReduction/POD.jl @@ -39,9 +39,9 @@ mutable struct POD <: AbstractDRProblem spectrum::Any # constructors function POD(snaps; - min_renergy = 1.0, - min_nmodes::Int = 1, - max_nmodes::Int = length(snaps[1])) + min_renergy = 1.0, + min_nmodes::Int = 1, + max_nmodes::Int = length(snaps[1])) nmodes = min_nmodes errorhandle(snaps, nmodes, min_renergy, min_nmodes, max_nmodes) new(snaps, min_renergy, min_nmodes, max_nmodes, nmodes, missing, 1.0, missing) @@ -66,7 +66,7 @@ end function reduce!(pod::POD, alg::SVD) u, s, v = _svd(pod.snapshots; alg.kwargs...) pod.nmodes, pod.renergy = determine_truncation(s, pod.min_nmodes, pod.max_nmodes, - pod.min_renergy) + pod.min_renergy) pod.rbasis = u[:, 1:(pod.nmodes)] pod.spectrum = s nothing @@ -94,6 +94,6 @@ function Base.show(io::IO, pod::POD) print(io, "POD \n") print(io, "Reduction Order = ", pod.nmodes, "\n") print(io, "Snapshot size = (", size(pod.snapshots, 1), ",", size(pod.snapshots[1], 2), - ")\n") + ")\n") print(io, "Relative Energy = ", pod.renergy, "\n") end diff --git a/src/ErrorHandle.jl b/src/ErrorHandle.jl index 2994f42..a867ecb 100644 --- a/src/ErrorHandle.jl +++ b/src/ErrorHandle.jl @@ -7,7 +7,7 @@ function errorhandle(data::AbstractMatrix, nmodes::Int, max_energy, min_nmodes, end function errorhandle(data::AbstractVector{T}, nmodes::Int, max_energy, min_nmodes, - max_nmodes) where {T <: AbstractVector} + max_nmodes) where {T <: AbstractVector} @assert size(data[1], 1)>1 "State vector is expected to be vector valued." s = min(size(data, 1), size(data[1], 1)) @assert 0 eq.lhs, old_observed); dvs; ModelingToolkit.get_states(sys)] + fullstates = [map(eq -> eq.lhs, old_observed); dvs; ModelingToolkit.get_unknowns(sys)] new_observed = [old_observed; linear_projection_eqs] new_sorted_observed = ModelingToolkit.topsort_equations(new_observed, fullstates; - kwargs...) + kwargs...) @set! sys.observed = new_sorted_observed inv_dict = Dict(Symbolics.scalarize(ŷ .=> V' * dvs)) # reduced vars to original vars diff --git a/src/utils.jl b/src/utils.jl index 69b2b54..bab1754 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -27,16 +27,41 @@ end """ $(SIGNATURES) -Given a vector of expressions `exprs` and variables `vars`, -returns a sparse coefficient matrix `A`, constant terms `c` and nonlinear terms `n`, -such that `exprs = A * vars + c + n`, -where the constant terms do not contain any variables in `vars`. +Returns `true` is `expr` contains variables in `dvs` only and does not contain `iv`. + +""" +function only_dvs(expr, dvs, iv) + if isequal(expr, iv) + return false + elseif expr in dvs + return true + elseif SymbolicUtils.istree(expr) + args = arguments(expr) + for arg in args + if only_dvs(arg, dvs, iv) + return true + end + end + end + return false +end + +""" +$(SIGNATURES) + +Given a vector of expressions `exprs`, variables `vars` and a single variable `iv`, +where `vars(iv)` are dependent variables of `iv`, +returns a sparse coefficient matrix `A`, other terms `g` and nonlinear terms `F`, +such that `exprs = A * vars(iv) + g(iv) + F(vars(iv))`, +where the nonlinear terms are functions of `vars` only and do not contain `iv`. Variables in `vars` must be unique. """ -function linear_terms(exprs::AbstractVector, vars) +function separate_terms(exprs::AbstractVector, vars, iv) vars = Symbolics.unwrap.(vars) exprs = Symbolics.unwrap.(exprs) + # expand is helpful for separating terms but is harmful for generating efficient runtime functions + exprs = expand.(exprs) linear_I = Int[] # row idx for sparse matrix linear_J = Int[] # col idx for sparse matrix linear_V = Float64[] # values @@ -55,26 +80,26 @@ function linear_terms(exprs::AbstractVector, vars) nothing end - const_terms = similar(exprs, Num) # create a vector of the same size - const_terms .= 0 # manually set to Int 0 because `Num` doesn't have a corresponding zero + other_terms = similar(exprs, Num) # create a vector of the same size + other_terms .= 0 # manually set to Int 0 because `Num` doesn't have a corresponding zero nonlinear_terms = similar(exprs, Num) nonlinear_terms .= 0 - # check if the expr is a constant or nolinear term about vars + # check if the expr is a nolinear term about vars only # and add it to the corresponding collection - @inline function const_nonlinear(i, expr) - # expr is nonlinear if it contains any variable in vars - if Symbolics.has_vars(expr, vars) + @inline function other_nonlinear(i, expr) + # expr is nonlinear if it contains vars only + if only_dvs(expr, vars, iv) nonlinear_terms[i] += expr - else # expr is constant if it doesn't have vars - const_terms[i] += expr + else + other_terms[i] += expr end nothing end for (i, expr) in enumerate(exprs) if expr isa Number # just a number, e.g. Int, Float64 - const_terms[i] = expr + other_terms[i] = expr elseif expr in vars # expr is a variables in vars push_sparse_coeff(i, expr, 1) elseif SymbolicUtils.ismul(expr) && length(expr.dict) == 1 @@ -82,21 +107,21 @@ function linear_terms(exprs::AbstractVector, vars) if base in vars && exp == 1 # a var with a coeff push_sparse_coeff(i, base, expr.coeff) else - const_nonlinear(i, expr) + other_nonlinear(i, expr) end elseif SymbolicUtils.isadd(expr) - const_terms[i] += expr.coeff + other_terms[i] += expr.coeff for (term, coeff) in expr.dict if term in vars push_sparse_coeff(i, term, coeff) else - const_nonlinear(i, term * coeff) + other_nonlinear(i, term * coeff) end end else - const_nonlinear(i, expr) + other_nonlinear(i, expr) end end linear = sparse(linear_I, linear_J, linear_V, length(exprs), length(vars)) - return linear, const_terms, nonlinear_terms + return linear, other_terms, nonlinear_terms end diff --git a/test/Project.toml b/test/Project.toml index 8a8f06b..3a568ac 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,13 +4,11 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" -DifferentialEquations = "7.6" -MethodOfLines = "0.6, 0.7" -ModelingToolkit = "8.33" -SafeTestsets = "0.0.1" -Symbolics = "4.13" +DifferentialEquations = "7.13" +MethodOfLines = "0.11" +ModelingToolkit = "9" +SafeTestsets = "0.1" diff --git a/test/deim.jl b/test/deim.jl index 05d61d7..d8d2d3d 100644 --- a/test/deim.jl +++ b/test/deim.jl @@ -40,7 +40,7 @@ pod_dim = 3 deim_sys = @test_nowarn deim(simp_sys, snapshot_simpsys, pod_dim) # check the number of dependent variables in the new system -@test length(ModelingToolkit.get_states(deim_sys)) == pod_dim +@test length(ModelingToolkit.get_unknowns(deim_sys)) == pod_dim deim_prob = ODEProblem(deim_sys, nothing, tspan) diff --git a/test/runtests.jl b/test/runtests.jl index 821f7d6..4953946 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,14 @@ using SafeTestsets -@safetestset "Quality Assurance" begin include("qa.jl") end -@safetestset "POD" begin include("DataReduction.jl") end -@safetestset "utils" begin include("utils.jl") end -@safetestset "DEIM" begin include("deim.jl") end +@safetestset "Quality Assurance" begin + include("qa.jl") +end +@safetestset "POD" begin + include("DataReduction.jl") +end +@safetestset "utils" begin + include("utils.jl") +end +@safetestset "DEIM" begin + include("deim.jl") +end diff --git a/test/utils.jl b/test/utils.jl index 2515cf8..03a3ef7 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,32 +1,32 @@ using Test, ModelOrderReduction -using Symbolics +using ModelingToolkit @variables t w(t) x(t) y(t) z(t) -@testset "linear_terms" begin - @testset "linear_terms full" begin - vars = [x, y, z] - exprs = [3.0x + 4.5y + 6.0 - 2.0z + 3.4w + 7.0 + sin(x) - 9.8 + x * (1.0 - y) - 5.6y + 1.3z^2] - A, c, n = ModelOrderReduction.linear_terms(exprs, vars) - @test size(A) == (length(exprs), length(vars)) - @test A == [3.0 4.5 0.0 - 0.0 0.0 2.0 - 0.0 0.0 0.0 - 0.0 5.6 0.0] - @test length(c) == length(exprs) - @test isequal(c, [6.0, 3.4w + 7.0, 9.8, 0.0]) - @test length(n) == length(exprs) - @test isequal(n, [0.0, sin(x), x * (1.0 - y), 1.3z^2]) - end +@testset "separate_terms" begin + # @testset "separate_terms full" begin + # vars = [x, y, z] + # exprs = [3.0x + 4.5y + 6.0 + # 2.0z + 3.4w + 7.0 + sin(x) + # 9.8 + x * (1.0 - y) + # 5.6y + 1.3z^2] + # A, c, n = ModelOrderReduction.separate_terms(exprs, vars, t) + # @test size(A) == (length(exprs), length(vars)) + # @test A == [3.0 4.5 0.0 + # 0.0 0.0 2.0 + # 0.0 0.0 0.0 + # 0.0 5.6 0.0] + # @test length(c) == length(exprs) + # @test isequal(c, [6.0, 3.4w + 7.0, 9.8, 0.0]) + # @test length(n) == length(exprs) + # @test isequal(n, [0.0, sin(x), x * (1.0 - y), 1.3z^2]) + # end - @testset "linear_terms empty exprs" begin + @testset "separate_terms empty exprs" begin vars = [x, y, z] exprs = Vector{Num}(undef, 4) fill!(exprs, false) - A, c, n = ModelOrderReduction.linear_terms(exprs, vars) + A, c, n = ModelOrderReduction.separate_terms(exprs, vars, t) @test size(A) == (length(exprs), length(vars)) @test iszero(A) @test length(c) == length(exprs) @@ -35,10 +35,10 @@ using Symbolics @test iszero(n) end - @testset "linear_terms diagonal" begin + @testset "separate_terms diagonal" begin vars = [x, y, z] exprs = [x, 2y, 3z, 4w] - A, c, n = ModelOrderReduction.linear_terms(exprs, vars) + A, c, n = ModelOrderReduction.separate_terms(exprs, vars, t) @test size(A) == (length(exprs), length(vars)) @test A == [1.0 0.0 0.0 0.0 2.0 0.0 @@ -50,12 +50,12 @@ using Symbolics @test iszero(n) end - @testset "linear_terms nonunique vars" begin + @testset "separate_terms nonunique vars" begin vars = [x, y, y] exprs = [3.0x + 4.5y + 6.0 2.0z + 3.4w + 7.0 + sin(x) 9.8 + x * (1.0 - y) 5.6y + 1.3z^2] - @test_throws ArgumentError ModelOrderReduction.linear_terms(exprs, vars) + @test_throws ArgumentError ModelOrderReduction.separate_terms(exprs, vars, t) end end