Proper Orthogonal Decomposition (POD) method

.github/workflows/documentation.yml

@@ -0,0 +1,23 @@
+name: Documentation
+
+on:
+ push:
+ branches:
+ - main
+ pull_request:
+
+jobs:
+ build:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v2
+ - uses: julia-actions/setup-julia@latest
+ with:
+ version: '1'
+ - name: Install dependencies
+ run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+ - name: Build and deploy
+ env:
+ GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
+ DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
+ run: julia --project=docs/ docs/make.jl

Project.toml

@@ -3,11 +3,15 @@ uuid = "207801d6-6cee-43a9-ad0c-f0c64933efa0"
 authors = ["Bowen S. Zhu"]
 version = "0.1.0"
 
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
 [compat]
 julia = "1.6"
 
 [extras]
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Test", "SafeTestsets"]

docs/Project.toml

@@ -0,0 +1,9 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+
+[compat]
+Documenter = "0.27"

docs/make.jl

@@ -0,0 +1,8 @@
+using Documenter, ModelOrderReduction
+
+include("pages.jl")
+
+makedocs(sitename = "ModelOrderReduction.jl",
+ pages = pages)
+
+deploydocs(repo = "github.com/SciML/ModelOrderReduction.jl.git")

docs/pages.jl

@@ -0,0 +1,5 @@
+pages = [
+ "Home" => "index.md",
+ "Functions" => "functions.md",
+ "Tutorials" => Any["tutorials/proper_orthogonal_decomposition.md"],
+]

docs/src/functions.md

@@ -0,0 +1,6 @@
+```@index
+```
+
+```@docs
+proper_orthogonal_decomposition(snapshots::AbstractMatrix, dim::Integer)
+```

docs/src/index.md

@@ -0,0 +1,10 @@
+# ModelOrderReduction.jl
+
+## Installation
+
+To install ModelOrderReduction.jl, use the Julia package manager:
+
+```julia
+using Pkg
+Pkg.add(url = "https:/SciML/ModelOrderReduction.jl")
+```

docs/src/tutorials/proper_orthogonal_decomposition.md

@@ -0,0 +1,145 @@
+# Proper Orthogonal Decomposition
+
+The following fluid dynamics example is adapted from this 
+[YouTube video](https://youtu.be/F7rWoxeGrko) which demonstrates the the Stable Fluids 
+algorithm for a 2d fluid flow. 
+
+```@example pod
+using FFTW, Plots, Interpolations, LinearAlgebra
+N_POINTS = 250
+KINEMATIC_VISCOSITY = 0.0001
+TIME_STEP_LENGTH = 0.01
+N_TIME_STEPS = 100
+ELEMENT_LENGTH = 1.0 / (N_POINTS - 1)
+X_INTERVAL = 0.0:ELEMENT_LENGTH:1.0
+Y_INTERVAL = 0.0:ELEMENT_LENGTH:1.0
+function backtrace!(backtraced_positions, original_positions, direction)
+ backtraced_positions[:] = mod1.(original_positions - TIME_STEP_LENGTH * direction, 1.0)
+end
+function interpolate_positions!(field_interpolated, field, interval_x, interval_y,
+ query_points_x, query_points_y)
+ interpolator = LinearInterpolation((interval_x, interval_y), field)
+ field_interpolated[:] = interpolator.(query_points_x, query_points_y)
+end
+function stable_fluids_fft()
+ coordinates_x = [x for x in X_INTERVAL, y in Y_INTERVAL]
+ coordinates_y = [y for x in X_INTERVAL, y in Y_INTERVAL]
+ wavenumbers_1d = fftfreq(N_POINTS) .* N_POINTS
+ wavenumbers_x = [k_x for k_x in wavenumbers_1d, k_y in wavenumbers_1d]
+ wavenumbers_y = [k_y for k_x in wavenumbers_1d, k_y in wavenumbers_1d]
+ wavenumbers_norm = [norm([k_x, k_y]) for k_x in wavenumbers_1d, k_y in wavenumbers_1d]
+ decay = exp.(-TIME_STEP_LENGTH .* KINEMATIC_VISCOSITY .* wavenumbers_norm .^ 2)
+ wavenumbers_norm[iszero.(wavenumbers_norm)] .= 1.0
+ normalized_wavenumbers_x = wavenumbers_x ./ wavenumbers_norm
+ normalized_wavenumbers_y = wavenumbers_y ./ wavenumbers_norm
+ force_x = 100.0 .* (exp.(-1.0 / (2 * 0.005) *
+ ((coordinates_x .- 0.2) .^ 2 + (coordinates_y .- 0.45) .^ 2))
+ -
+ exp.(-1.0 / (2 * 0.005) *
+ ((coordinates_x .- 0.8) .^ 2 + (coordinates_y .- 0.55) .^ 2)))
+ backtraced_coordinates_x = similar(coordinates_x)
+ backtraced_coordinates_y = similar(coordinates_y)
+ velocity_x = similar(coordinates_x)
+ velocity_y = similar(coordinates_y)
+ velocity_x_prev = similar(velocity_x)
+ velocity_y_prev = similar(velocity_y)
+ velocity_x_fft = similar(velocity_x)
+ velocity_y_fft = similar(velocity_y)
+ pressure_fft = similar(coordinates_x)
+ curls = similar(velocity_x, length(X_INTERVAL) - 1, length(Y_INTERVAL) - 1,
+ N_TIME_STEPS)
+ for iter in 1:N_TIME_STEPS
+ time_current = (iter - 1) * TIME_STEP_LENGTH
+ pre_factor = max(1 - time_current, 0.0)
+ velocity_x_prev += TIME_STEP_LENGTH * pre_factor * force_x
+ backtrace!(backtraced_coordinates_x, coordinates_x, velocity_x_prev)
+ backtrace!(backtraced_coordinates_y, coordinates_y, velocity_y_prev)
+ interpolate_positions!(velocity_x, velocity_x_prev, X_INTERVAL, Y_INTERVAL,
+ backtraced_coordinates_x, backtraced_coordinates_y)
+ interpolate_positions!(velocity_y, velocity_y_prev, X_INTERVAL, Y_INTERVAL,
+ backtraced_coordinates_x, backtraced_coordinates_y)
+ velocity_x_fft = fft(velocity_x)
+ velocity_y_fft = fft(velocity_y)
+ velocity_x_fft .*= decay
+ velocity_y_fft .*= decay
+ pressure_fft = (velocity_x_fft .* normalized_wavenumbers_x
+ +
+ velocity_y_fft .* normalized_wavenumbers_y)
+ velocity_x_fft -= pressure_fft .* normalized_wavenumbers_x
+ velocity_y_fft -= pressure_fft .* normalized_wavenumbers_y
+ velocity_x = real(ifft(velocity_x_fft))
+ velocity_y = real(ifft(velocity_y_fft))
+ velocity_x_prev = velocity_x
+ velocity_y_prev = velocity_y
+ d_u__d_y = diff(velocity_x, dims = 2)[2:end, :]
+ d_v__d_x = diff(velocity_y, dims = 1)[:, 2:end]
+ curls[:, :, iter] = d_u__d_y - d_v__d_x
+ end
+ return curls
+end
+curls = stable_fluids_fft(); nothing
+```
+
+We got the data of curl at a few time steps. The time-averaged field is visualized as 
+follows.
+
+```@example pod
+using Statistics
+time_averaged = dropdims(mean(curls; dims = 3); dims = 3)
+heatmap(X_INTERVAL, Y_INTERVAL, time_averaged', aspect_ratio = :equal, size = (1700, 1600),
+ xlabel = "x", ylabel = "y")
+```
+
+We subtract the time-averaged field from each individual snapshot and change the shape to 
+our snapshot matrix:
+
+```@example pod
+snapshots = curls .- time_averaged
+snapshots = reshape(snapshots, :, N_TIME_STEPS)
+```
+
+```math
+D =
+\begin{pmatrix}
+ d_{11} & d_{12} & \cdots & d_{1m} \\
+ \vdots & \vdots & & \vdots \\
+ d_{nm} & d_{n2} & \cdots & d_{nm}
+\end{pmatrix}
+=
+\begin{pmatrix}
+ d(x_1,y_1,t_1) & d(x_1,y_1,t_2) & \cdots & d(x_1,y_1,t_m) \\
+ \vdots & \vdots & & \vdots \\
+ d(x_{N_x},y_{N_y},t_1) & d(x_{N_x},y_{N_y},t_2) & \cdots & d(x_{N_x},y_{N_y},t_m)
+\end{pmatrix}
+```
+
+where each column represents one snapshot. ``n`` is the number of discrete points and ``m``
+is the number of snapshots. Each element ``d_{ij}`` of ``D`` is the scalar value at point 
+``i`` measured at time ``j``. 
+
+Then we compute the POD with dimension, say, 3.
+
+```@example pod
+using ModelOrderReduction
+pod_dim = 3
+pod_basis, singular_vals = proper_orthogonal_decomposition(snapshots, pod_dim); nothing
+```
+
+As an illustration, let's see the first 3 POD modes.
+
+```@example pod
+mode_plots = Any[]
+for dim in 1:pod_dim
+ mode = reshape((@view pod_basis[:, dim]),
+ (length(X_INTERVAL) - 1, length(Y_INTERVAL) - 1))
+ push!(mode_plots,
+ heatmap(X_INTERVAL, Y_INTERVAL, mode', aspect_ratio = :equal,
+ title = "mode $dim", xlabel = "x", ylabel = "y"))
+ push!(mode_plots,
+ plot(vec((@view pod_basis[:, dim])' * snapshots), ylim = (-15, 15),
+ title = "mode $dim", xlabel = "t", ylabel = "a", lw = 3))
+end
+plot(mode_plots..., layout = (pod_dim, 2), legend = false, size = (3000, 2000))
+```
+
+The figures on the right show the time coefficient of the first 3 POD modes.

src/ModelOrderReduction.jl

@@ -1,5 +1,5 @@
 module ModelOrderReduction
 
-# Write your package code here.
+include("proper_orthogonal_decomposition.jl")
 
 end

src/proper_orthogonal_decomposition.jl

@@ -0,0 +1,46 @@
+using LinearAlgebra
+
+"""
+ proper_orthogonal_decomposition(snapshots::AbstractMatrix, dim::Integer)
+
+Compute the Proper Orthogonal Decomposition (POD) for a model with a snapshot matrix and a
+POD dimension. Return the POD basis and singular values.
+
+The results are sorted in descending order of singular values.
+
+The POD dimension should be smaller than the number of snapshots.
+
+# Arguments
+- `snapshots::AbstractMatrix`: a matrix with snapshots in the columns.
+- `dim::Integer`: the POD dimension.
+
+# Examples
+```jldoctest
+julia> n_eq = 10; # number of equations
+
+julia> n_snapshot = 6; # number of snapshots
+
+julia> snapshots = rand(n_eq, n_snapshot);
+
+julia> dim = 3; # POD dimension
+
+julia> pod_basis, singular_vals = proper_orthogonal_decomposition(snapshots, dim);
+
+julia> size(pod_basis)
+(10, 3)
+
+julia> size(singular_vals)
+(3,)
+```
+"""
+function proper_orthogonal_decomposition(snapshots::AbstractMatrix, dim::Integer)
+ eigen_vecs, singuler_vals = svd(snapshots)
+ if size(snapshots, 2) < dim
+ @warn "The POD dimension is larger than the number of snapshots"
+ return eigen_vecs, singuler_vals
+ else
+ return (@view eigen_vecs[:, 1:dim]), (@view singuler_vals[1:dim])
+ end
+end
+
+export proper_orthogonal_decomposition

test/proper_orthogonal_decomposition_test.jl

@@ -0,0 +1,22 @@
+using ModelOrderReduction, Test
+
+@testset "POD $T" for T in (Float32, Float64)
+ n_eq = 10 # number of equations
+ n_snapshot = 6 # number of snapshots
+ snapshots = rand(T, n_eq, n_snapshot)
+
+ dim₁ = 3 # POD dimension
+ pod_basis₁, singular_vals₁ = proper_orthogonal_decomposition(snapshots, dim₁)
+ @test size(pod_basis₁) == (n_eq, dim₁)
+ @test size(singular_vals₁, 1) == dim₁
+ @test eltype(pod_basis₁) == T
+ @test eltype(singular_vals₁) == T
+
+ dim₂ = 8 # larger than the number of snapshots
+ pod_basis₂, singular_vals₂ = @test_logs (:warn,) proper_orthogonal_decomposition(snapshots,
+ dim₂)
+ @test size(pod_basis₂) == (n_eq, n_snapshot)
+ @test size(singular_vals₂, 1) == n_snapshot
+ @test eltype(pod_basis₂) == T
+ @test eltype(singular_vals₂) == T
+end

test/runtests.jl

@@ -1,6 +1,3 @@
-using ModelOrderReduction
-using Test
+using SafeTestsets
 
-@testset "ModelOrderReduction.jl" begin
- # Write your tests here.
-end
+@safetestset "Proper Orthogonal Decomposition Test" begin include("proper_orthogonal_decomposition_test.jl") end