Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 0 additions & 2 deletions .JuliaFormatter.toml

This file was deleted.

14 changes: 10 additions & 4 deletions .github/workflows/FormatCheck.yml
Original file line number Diff line number Diff line change
@@ -1,13 +1,19 @@
name: "Format Check"
name: format-check

on:
push:
branches:
- 'master'
- 'main'
- 'release-'
tags: '*'
pull_request:

jobs:
format-check:
name: "Format Check"
uses: "SciML/.github/.github/workflows/format-check.yml@v1"
runic:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- uses: fredrikekre/runic-action@v1
with:
version: '1'
12 changes: 8 additions & 4 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,17 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)

include("pages.jl")

makedocs(sitename = "ModelOrderReduction.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)
format = Documenter.HTML(
assets = ["assets/favicon.ico"],
canonical = "https://docs.sciml.ai/ModelOrderReduction/stable/"
),
pages = pages
)

deploydocs(repo = "github.com/SciML/ModelOrderReduction.jl.git")
32 changes: 19 additions & 13 deletions src/DataReduction/POD.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ using TSVD: tsvd
using RandomizedLinAlg: rsvd

function matricize(VoV::Vector{Vector{T}}) where {T}
reduce(hcat, VoV)
return reduce(hcat, VoV)
end

function _svd(data::Vector{Vector{T}}; kwargs...) where {T}
Expand Down Expand Up @@ -38,17 +38,19 @@ mutable struct POD <: AbstractDRProblem
renergy::Any
spectrum::Any
# constructors
function POD(snaps;
function POD(
snaps;
min_renergy = 1.0,
min_nmodes::Int = 1,
max_nmodes::Int = length(snaps[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)
return new(snaps, min_renergy, min_nmodes, max_nmodes, nmodes, missing, 1.0, missing)
end
function POD(snaps, nmodes::Int)
errorhandle(snaps, nmodes, 0.0, nmodes, nmodes)
new(snaps, 0.0, nmodes, nmodes, nmodes, missing, 1.0, missing)
return new(snaps, 0.0, nmodes, nmodes, nmodes, missing, 1.0, missing)
end
end

Expand All @@ -66,11 +68,13 @@ 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.renergy = determine_truncation(
s, pod.min_nmodes, pod.max_nmodes,
pod.min_renergy
)
pod.rbasis = u[:, 1:(pod.nmodes)]
pod.spectrum = s
nothing
return nothing
end

function reduce!(pod::POD, alg::TSVD)
Expand All @@ -79,7 +83,7 @@ function reduce!(pod::POD, alg::TSVD)
pod.renergy = sum(s) / (sum(s) + (n_max - pod.nmodes) * s[end])
pod.rbasis = u
pod.spectrum = s
nothing
return nothing
end

function reduce!(pod::POD, alg::RSVD)
Expand All @@ -88,13 +92,15 @@ function reduce!(pod::POD, alg::RSVD)
pod.renergy = sum(s) / (sum(s) + (n_max - pod.nmodes) * s[end])
pod.rbasis = u
pod.spectrum = s
nothing
return nothing
end

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")
print(io, "Relative Energy = ", pod.renergy, "\n")
print(
io, "Snapshot size = (", size(pod.snapshots, 1), ",", size(pod.snapshots[1], 2),
")\n"
)
return print(io, "Relative Energy = ", pod.renergy, "\n")
end
22 changes: 12 additions & 10 deletions src/ErrorHandle.jl
Original file line number Diff line number Diff line change
@@ -1,16 +1,18 @@
function errorhandle(data::AbstractMatrix, nmodes::Int, max_energy, min_nmodes, max_nmodes)
@assert size(data, 1)>1 "State vector is expected to be vector valued."
@assert size(data, 1) > 1 "State vector is expected to be vector valued."
s = minimum(size(data))
@assert 0<nmodes<=s "Number of modes should be in {1,2,...,$s}."
@assert min_nmodes<=max_nmodes "Minimum number of modes must lie below maximum number of modes"
@assert 0.0<=max_energy<=1.0 "Maximum relative energy must be in [0,1]"
@assert 0 < nmodes <= s "Number of modes should be in {1,2,...,$s}."
@assert min_nmodes <= max_nmodes "Minimum number of modes must lie below maximum number of modes"
return @assert 0.0 <= max_energy <= 1.0 "Maximum relative energy must be in [0,1]"
end

function errorhandle(data::AbstractVector{T}, nmodes::Int, max_energy, min_nmodes,
max_nmodes) where {T <: AbstractVector}
@assert size(data[1], 1)>1 "State vector is expected to be vector valued."
function errorhandle(
data::AbstractVector{T}, nmodes::Int, max_energy, min_nmodes,
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<nmodes<=s "Number of modes should be in {1,2,...,$s}."
@assert min_nmodes<=max_nmodes "Minimum number of modes must lie below maximum number of modes"
@assert 0.0<=max_energy<=1.0 "Maximum relative energy must be in [0,1]"
@assert 0 < nmodes <= s "Number of modes should be in {1,2,...,$s}."
@assert min_nmodes <= max_nmodes "Minimum number of modes must lie below maximum number of modes"
return @assert 0.0 <= max_energy <= 1.0 "Maximum relative energy must be in [0,1]"
end
4 changes: 2 additions & 2 deletions src/ModelOrderReduction.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@ module ModelOrderReduction
using DocStringExtensions: DocStringExtensions, FUNCTIONNAME, SIGNATURES, TYPEDSIGNATURES

using ModelingToolkit: ModelingToolkit, @variables, Differential, Equation, Num, ODESystem,
SymbolicUtils, Symbolics, arguments, build_function, complete,
expand, substitute, tearing_substitution
SymbolicUtils, Symbolics, arguments, build_function, complete,
expand, substitute, tearing_substitution
using LinearAlgebra: LinearAlgebra, /, \, mul!, qr, svd

using Setfield: Setfield, @set!
Expand Down
6 changes: 3 additions & 3 deletions src/Types.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,20 +7,20 @@ abstract type AbstractSVD end
struct SVD <: AbstractSVD
kwargs::Any
function SVD(; kwargs...)
new(kwargs)
return new(kwargs)
end
end

struct TSVD <: AbstractSVD
kwargs::Any
function TSVD(; kwargs...)
new(kwargs)
return new(kwargs)
end
end

struct RSVD <: AbstractSVD
p::Int
function RSVD(p::Int = 0)
new(p)
return new(p)
end
end
22 changes: 14 additions & 8 deletions src/deim.jl
Original file line number Diff line number Diff line change
Expand Up @@ -64,10 +64,12 @@ the ``\\rho_i``-th column of the identity matrix ``I_n\\in\\mathbb R^{n\\times n
- `reduced_rhss`: the right-hand side of ROM.
- `linear_projection_eqs`: the linear projection mapping ``\\mathbf y=V\\hat{\\mathbf y}``.
"""
function deim(full_vars::AbstractVector, linear_coeffs::AbstractMatrix,
function deim(
full_vars::AbstractVector, linear_coeffs::AbstractMatrix,
constant_part::AbstractVector, nonlinear_part::AbstractVector,
reduced_vars::AbstractVector, linear_projection_matrix::AbstractMatrix,
nonlinear_projection_matrix::AbstractMatrix; kwargs...)
nonlinear_projection_matrix::AbstractMatrix; kwargs...
)
# rename variables for convenience
y = full_vars
A = linear_coeffs
Expand All @@ -91,7 +93,7 @@ function deim(full_vars::AbstractVector, linear_coeffs::AbstractMatrix,
 = V' * A * V
ĝ = V' * g
reduced_rhss = Â * ŷ + ĝ + F̂
reduced_rhss, linear_projection_eqs
return reduced_rhss, linear_projection_eqs
end
"""
$(FUNCTIONNAME)(
Expand All @@ -118,9 +120,11 @@ The POD basis used for DEIM interpolation is obtained from the snapshot matrix o
nonlinear terms, which is computed by executing the runtime-generated function for
nonlinear expressions.
"""
function deim(sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer;
function deim(
sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer;
deim_dim::Integer = pod_dim, name::Symbol = Symbol(nameof(sys), :_deim),
kwargs...)::ODESystem
kwargs...
)::ODESystem
sys = deepcopy(sys)
@set! sys.name = name

Expand Down Expand Up @@ -166,15 +170,17 @@ function deim(sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer;
old_observed = ModelingToolkit.get_observed(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...)
new_sorted_observed = ModelingToolkit.topsort_equations(
new_observed, fullstates;
kwargs...
)
@set! sys.observed = new_sorted_observed

inv_dict = Dict(Symbolics.scalarize(ŷ .=> V' * dvs)) # reduced vars to original vars
@set! sys.defaults = merge(ModelingToolkit.defaults(sys), inv_dict)

# CRITICAL: Call complete() on the system before returning to ensure all subsystems,
# variables, and parameters are properly registered and namespaced. Without this,
# variables, and parameters are properly registered and namespaced. Without this,
# attempting to create an ODEProblem from the DEIM system will fail with errors about
# missing initial conditions for variables that should exist in the system.
# This is required due to changes in ModelingToolkit.jl's internal structure handling.
Expand Down
8 changes: 4 additions & 4 deletions src/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ function get_deqs(sys::ODESystem)::Tuple{Vector{Equation}, Vector{Equation}}
push!(others, eq)
end
end
deqs, others
return deqs, others
end

"""
Expand Down Expand Up @@ -78,7 +78,7 @@ function separate_terms(exprs::AbstractVector, vars, iv)
push!(linear_I, row_index)
push!(linear_J, idxmap[term])
push!(linear_V, value)
nothing
return nothing
end

other_terms = similar(exprs, Num) # create a vector of the same size
Expand All @@ -95,7 +95,7 @@ function separate_terms(exprs::AbstractVector, vars, iv)
else
other_terms[i] += expr
end
nothing
return nothing
end

for (i, expr) in enumerate(exprs)
Expand Down Expand Up @@ -125,4 +125,4 @@ function separate_terms(exprs::AbstractVector, vars, iv)
end
linear = sparse(linear_I, linear_J, linear_V, length(exprs), length(vars))
return linear, other_terms, nonlinear_terms
end
end
4 changes: 2 additions & 2 deletions test/DataReduction.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,15 @@ function lorenz_prob()
function lorenz!(du, u, p, t)
du[1] = p[1] * (u[2] - u[1])
du[2] = u[1] * (p[2] - u[3]) - u[2]
du[3] = u[1] * u[2] - p[3] * u[3]
return du[3] = u[1] * u[2] - p[3] * u[3]
end

u0 = [1, 0, 0]
p = [10, 28, 8 / 3]
tspan = (0, 100)
prob = ODEProblem(lorenz!, u0, tspan, p)
sol = solve(prob, Tsit5())
sol
return sol
end

@testset "POD-Utils" begin
Expand Down
18 changes: 12 additions & 6 deletions test/deim.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,20 @@ const γ = 2.0
const c = 0.05
f(v) = v * (v - 0.1) * (1.0 - v)
i₀(t) = 50000.0t^3 * exp(-15.0t)
eqs = [ε * Dt(v(x, t)) ~ ε^2 * Dxx(v(x, t)) + f(v(x, t)) - w(x, t) + c,
Dt(w(x, t)) ~ b * v(x, t) - γ * w(x, t) + c]
bcs = [v(x, 0.0) ~ 0.0,
eqs = [
ε * Dt(v(x, t)) ~ ε^2 * Dxx(v(x, t)) + f(v(x, t)) - w(x, t) + c,
Dt(w(x, t)) ~ b * v(x, t) - γ * w(x, t) + c,
]
bcs = [
v(x, 0.0) ~ 0.0,
w(x, 0) ~ 0.0,
Dx(v(0, t)) ~ -i₀(t),
Dx(v(L, t)) ~ 0.0]
domains = [x ∈ (0.0, L),
t ∈ (0.0, 14.0)]
Dx(v(L, t)) ~ 0.0,
]
domains = [
x ∈ (0.0, L),
t ∈ (0.0, 14.0),
]
ivs = [x, t]
dvs = [v(x, t), w(x, t)]
pde_sys = PDESystem(eqs, bcs, domains, ivs, dvs; name = :FN)
Expand Down
6 changes: 4 additions & 2 deletions test/qa.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,10 @@ using ModelOrderReduction, Aqua
#Aqua.find_persistent_tasks_deps(ModelOrderReduction)
Aqua.test_ambiguities(ModelOrderReduction, recursive = false)
Aqua.test_deps_compat(ModelOrderReduction)
Aqua.test_piracies(ModelOrderReduction,
treat_as_own = [])
Aqua.test_piracies(
ModelOrderReduction,
treat_as_own = []
)
Aqua.test_project_extras(ModelOrderReduction)
Aqua.test_stale_deps(ModelOrderReduction)
Aqua.test_unbound_args(ModelOrderReduction)
Expand Down
16 changes: 9 additions & 7 deletions test/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,12 @@ using ModelingToolkit
vars = [x, y, z]
exprs = [x, 2y, 3z, 4w]
A, c, n = ModelOrderReduction.separate_terms(exprs, vars, t)

# Test dimensions are correct
@test size(A) == (length(exprs), length(vars))
@test length(c) == length(exprs)
@test length(n) == length(exprs)

# Test that function doesn't error and returns expected types
end

Expand All @@ -32,10 +32,12 @@ using ModelingToolkit

@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]
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.separate_terms(exprs, vars, t)
end
end
end
Loading