Skip to content

Commit 4910f4a

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #15 from dynamic-queries/main
Initial Structure
2 parents 020f782 + a06a9f7 commit 4910f4a

9 files changed

Lines changed: 205 additions & 6 deletions

File tree

‎Project.toml‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,18 @@
11
name = "ModelOrderReduction"
22
uuid = "207801d6-6cee-43a9-ad0c-f0c64933efa0"
3-
authors = ["Bowen S. Zhu"]
3+
authors = ["Bowen S. Zhu", "Rahul Manavalan"]
44
version = "0.1.0"
55

6+
[deps]
7+
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
8+
RandomizedLinAlg = "0448d7d9-159c-5637-8537-fd72090fea46"
9+
TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
10+
611
[compat]
712
julia = "1.6"
813

914
[extras]
1015
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
1116

1217
[targets]
13-
test = ["Test"]
18+
test = ["Test"]

‎examples/POD/lorenz.jl‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
using ModelOrderReduction
2+
using Plots
3+
using OrdinaryDiffEq
4+
5+
function lorenz_prob()
6+
function lorenz!(du, u, p, t)
7+
du[1] = p[1] * (u[2] - u[1])
8+
du[2] = u[1] * (p[2] - u[3]) - u[2]
9+
du[3] = u[1] * u[2] - p[3] * u[3]
10+
end
11+
12+
u0 = [1, 0, 0]
13+
p = [10, 28, 8 / 3]
14+
tspan = (0, 100)
15+
prob = ODEProblem(lorenz!, u0, tspan, p)
16+
sol = solve(prob, Tsit5())
17+
sol
18+
end
19+
20+
sol = lorenz_prob()
21+
solution = Array(sol)
22+
plot(solution[1, :], solution[2, :], solution[3, :])
23+
# savefig("lorenz_attractor.png")
24+
25+
## Two way POD
26+
reducer = POD(solution, 2)
27+
reduce!(reducer, SVD())
28+
bases = reducer.rbasis
29+
plot(bases[:, 1], bases[:, 2], label = "POD2")
30+
# savefig("pod2.png")
31+
32+
## One way POD
33+
reducer = POD(solution, 1)
34+
reduce!(reducer, SVD())
35+
bases = reducer.rbasis
36+
plot(bases[:, 1], label = "POD1")
37+
# savefig("pod1.png")
38+
plot(solution[:, 3], label = 'z')
39+
# savefig("z.png")

‎src/DataReduction/POD.jl‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
function matricize(VoV::Vector{Vector{FT}}) where {FT}
2+
Matrix(reduce(hcat, VoV))
3+
end
4+
5+
mutable struct POD{FT} <: AbstractDRProblem
6+
snapshots::Union{Vector{Vector{FT}}, Matrix{FT}}
7+
nmodes::Int
8+
rbasis::Matrix{FT}
9+
renergy::FT
10+
11+
function POD(snaps::Vector{Vector{FT}}, nmodes::Int) where {FT}
12+
errorhandle(matricize(snaps), nmodes)
13+
new{eltype(snaps[1])}(snaps, nmodes, Array{FT, 2}(undef, size(snaps, 1), nmodes),
14+
FT(0))
15+
end
16+
17+
function POD(snaps::Matrix{FT}, nmodes::Int) where {FT}
18+
errorhandle(snaps, nmodes)
19+
new{eltype(snaps)}(snaps, nmodes, Array{FT, 2}(undef, size(snaps, 1), nmodes),
20+
FT(0))
21+
end
22+
23+
function POD(snaps::Adjoint{FT, Matrix{FT}}, nmodes::Int) where {FT}
24+
POD(Matrix(snaps), nmodes, Array{FT, 2}(undef, size(snaps, 1), nmodes), FT(0))
25+
end
26+
end
27+
28+
function reduce!(pod::POD{FT}, ::SVD) where {FT}
29+
op_matrix = pod.snapshots
30+
if typeof(pod.snapshots) == Vector{Vector{FT}}
31+
op_matrix = matricize(pod.snapshots)
32+
end
33+
u, s, v = svd(op_matrix)
34+
pod.rbasis .= u[:, 1:(pod.nmodes)]
35+
sr = s[1:(pod.nmodes)]
36+
pod.renergy = sum(s) / (sum(s) + (size(op_matrix, 1) - pod.nmodes) * s[end])
37+
nothing
38+
end
39+
40+
function reduce!(pod::POD{FT}, ::TSVD) where {FT}
41+
op_matrix = pod.snapshots
42+
if typeof(pod.snapshots) == Vector{Vector{FT}}
43+
op_matrix = matricize(pod.snapshots)
44+
end
45+
u, s, v = tsvd(op_matrix, pod.nmodes)
46+
pod.renergy = sum(s) / (sum(s) + (size(op_matrix, 1) - pod.nmodes) * s[end])
47+
pod.rbasis .= u
48+
nothing
49+
end
50+
51+
function reduce!(pod::POD{FT}, ::RSVD) where {FT}
52+
op_matrix = pod.snapshots
53+
if typeof(pod.snapshots) == Vector{Vector{FT}}
54+
op_matrix = matricize(pod.snapshots)
55+
end
56+
u, s, v = rsvd(op_matrix, pod.nmodes)
57+
pod.renergy = sum(s) / (sum(s) + (size(op_matrix, 1) - pod.nmodes) * s[end])
58+
pod.rbasis .= u
59+
nothing
60+
end
61+
62+
function Base.show(io::IO, pod::POD)
63+
print(io, "POD \n")
64+
print(io, "Reduction Order = ", pod.nmodes, "\n")
65+
print(io, "Snapshot size = (", size(pod.snapshots, 1), ",", size(pod.snapshots[1], 2),
66+
")\n")
67+
print(io, "Relative Energy = ", pod.renergy, "\n")
68+
end

‎src/ErrorHandle.jl‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
function errorhandle(data::Matrix{FT}, modes::IT) where {FT, IT}
2+
@assert size(data, 1)>1 "State vector is expected to be vector valued."
3+
s = size(data, 2)
4+
@assert (modes > 0)&(modes < s) "Number of modes should be in {1,2,...,$(s-1)}."
5+
end

‎src/ModelOrderReduction.jl‎

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,17 @@
11
module ModelOrderReduction
2+
#========================Data Reduction=========================================#
3+
include("Types.jl")
4+
include("ErrorHandle.jl")
25

3-
# Write your package code here.
6+
using LinearAlgebra
7+
using TSVD
8+
using RandomizedLinAlg
49

10+
include("DataReduction/POD.jl")
11+
12+
export SVD, TSVD, RSVD
13+
export POD, reduce!, matricize
14+
#========================Model Reduction========================================#
15+
16+
#===============================================================================#
517
end

‎src/Types.jl‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
abstract type AbstractReductionProblem end
2+
abstract type AbstractMORProblem <: AbstractReductionProblem end
3+
abstract type AbstractDRProblem <: AbstractReductionProblem end
4+
5+
abstract type AbstractSVD end
6+
struct SVD <: AbstractSVD end
7+
struct TSVD <: AbstractSVD end
8+
struct RSVD <: AbstractSVD end

‎test/DataReduction.jl‎

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
#--------- Data Reduction -----------------#
2+
3+
function lorenz_prob()
4+
function lorenz!(du, u, p, t)
5+
du[1] = p[1] * (u[2] - u[1])
6+
du[2] = u[1] * (p[2] - u[3]) - u[2]
7+
du[3] = u[1] * u[2] - p[3] * u[3]
8+
end
9+
10+
u0 = [1, 0, 0]
11+
p = [10, 28, 8 / 3]
12+
tspan = (0, 100)
13+
prob = ODEProblem(lorenz!, u0, tspan, p)
14+
sol = solve(prob, Tsit5())
15+
sol
16+
end
17+
18+
@testset "POD-Utils" begin
19+
solution = lorenz_prob()
20+
VoV = solution.u
21+
M = matricize(VoV)
22+
@test size(M, 1) == length(VoV[1]) # Parameters
23+
@test size(M, 2) == length(VoV) # Time
24+
end
25+
26+
@testset "POD - Attractor Test" begin
27+
sol = lorenz_prob()
28+
solution = Array(sol)
29+
30+
order = 2
31+
solver = SVD()
32+
reducer = POD(solution, order)
33+
reduce!(reducer, solver)
34+
35+
reducer.renergy
36+
37+
# Ad-hoc tests. To be checked with Chris.
38+
@test size(reducer.rbasis, 2) == reducer.nmodes
39+
@test size(reducer.rbasis, 1) == size(solution, 1)
40+
@test reducer.renergy > 0.9
41+
42+
order = 2
43+
solver = TSVD()
44+
reducer = POD(solution, order)
45+
reduce!(reducer, solver)
46+
47+
# Ad-hoc tests. To be checked with Chris.
48+
@test size(reducer.rbasis, 2) == reducer.nmodes
49+
@test size(reducer.rbasis, 1) == size(solution, 1)
50+
51+
order = 2
52+
solver = RSVD()
53+
reducer = POD(solution, order)
54+
reduce!(reducer, solver)
55+
56+
# Ad-hoc tests. To be checked with Chris.
57+
@test size(reducer.rbasis, 2) == reducer.nmodes
58+
@test size(reducer.rbasis, 1) == size(solution, 1)
59+
end

‎test/Project.toml‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
[deps]
2+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
3+
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

‎test/runtests.jl‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
using ModelOrderReduction
22
using Test
3+
using OrdinaryDiffEq
34

4-
@testset "ModelOrderReduction.jl" begin
5-
# Write your tests here.
6-
end
5+
include("DataReduction.jl")
6+
#---------- Model Reduction ----------------#

0 commit comments

Comments
 (0)