GitHub

Proper orthogonal decomposition (POD)

In this example, we will demonstrate the use of proper orthogonal decomposition (POD) for decomposing a flow field into basis modes.

using ILMPostProcessing
using Plots

Get the flow field data

First, we need some flow field data to analyze. For this purpose, we will use ViscousFlow.jl to get snapshots of the flow for a flat plate at 30 degrees angle of attack at Reynolds number 100.

using ViscousFlow

my_params = Dict()
my_params["Re"] = 100
my_params["freestream speed"] = 1.0 # in x-dir
my_params["freestream angle"] = 0.0 # relative to horizontal

xlim = (-1.0, 5.0)
ylim = (-2.0, 2.0)
my_params["grid Re"] = 4.0

g = setup_grid(xlim, ylim, my_params)

Δs = surface_point_spacing(g, my_params)
body = Plate(1.0, Δs)

cent = [0.5,0.0]
α = -30π/180
X = MotionTransform(cent, α)
joint = Joint(X)

m = RigidBodyMotion(joint, body)
x = init_motion_state(body, m)
update_body!(body, x, m)
sys = viscousflow_system(g, body, phys_params = my_params, motions = m);

u0 = init_sol(sys)
tspan = (0.0, 20.0)
integrator = init(u0, tspan, sys)

# Solve to 10 convective time units
step!(integrator, 10)

Assemble snapshots of the velocity field from the solution data

Here, we make use of the capability of the velocity function to generate an array of velocity fields at a range of times. We will save every 5th time step in this array.

sol = integrator.sol
tpod = sol.t[2:5:end]
X = velocity(sol, sys, tpod);

Perform the POD

The POD is simply performed with the PODModes function. This provides a structure containing the modes (phi), the expansion coefficients (a), and the modal energies (lambda). By default, PODModes retains 99% of the energy. This can be changed with the optional argument tolerance.

modes = pod(X);

The a array is of size $N_t \times r$, where $N_t$ is the number of time values, and $r$ is the number of modes. The modes are ordered from highest energy to lowest energy.

modes.a
101×10 Matrix{Float64}:
  0.139488    0.00736536  -0.0768184  …  -0.0256712    -0.0123668
  0.138101    0.00750601  -0.0758466     -0.0178839    -0.00779976
  0.136146    0.00779498  -0.075173      -0.00937233   -0.00320754
  0.133823    0.00853107  -0.0748546     -0.00108265    0.000901443
  0.131309    0.00992053  -0.0747563      0.00607686    0.00419214
  0.128777    0.0120599   -0.0746398  …   0.0115034     0.0065311
  0.12639     0.0149328   -0.0742144      0.014987      0.00791158
  0.124287    0.0184193   -0.0731801      0.0166409     0.00834968
  0.122573    0.0223155   -0.0712662      0.0167553     0.00783233
  0.121308    0.0263596   -0.0682605      0.0156501     0.00634241
  ⋮                                   ⋱                
 -0.109734   -0.0348703   -0.0499863     -0.0118352     0.00881401
 -0.111545   -0.0231225   -0.0565372     -0.00877868    0.0102615
 -0.112448   -0.0104652   -0.0612092     -0.00464563    0.0102618
 -0.112411    0.0028515   -0.0638602  …   0.000203529   0.00879315
 -0.111426    0.0165666   -0.064389       0.00535478    0.00599874
 -0.1095      0.0304125   -0.0627401      0.010392      0.00215723
 -0.106666    0.0441202   -0.0589079      0.0149453    -0.00236466
 -0.102978    0.0574227   -0.0529397      0.0187283    -0.00716857
 -0.0985084   0.070059    -0.0449373  …   0.0215588    -0.0118793

In this case, 7 modes were retained, at 51 times.

If we wanted to re-assemble the modes and coefficients to recover the flow at some time instant, we could use the mapreduce function, e.g.,

vel_assemble = mapreduce((aj, phi_j) -> aj .* phi_j, +, modes.a[end,:], modes.phi) + modes.Xmean
CartesianGrids.Edges{CartesianGrids.Primal, 160, 104, Float64, Vector{Float64}} data
u (in grid orientation)
103×160 Matrix{Float64}:
 1.02786   1.02968   1.03141   1.03312   …  0.987704  0.986192  0.985118
 1.02696   1.02874   1.03049   1.03222      0.986571  0.985142  0.983924
 1.02595   1.0277    1.02945   1.03119      0.985336  0.983939  0.98266
 1.02485   1.02658   1.02832   1.03007      0.98399   0.982594  0.981281
 1.02366   1.02539   1.02712   1.02886      0.982521  0.981118  0.979782
 1.0224    1.02411   1.02584   1.02758   …  0.980927  0.979515  0.978164
 1.02106   1.02276   1.02448   1.02621      0.979204  0.977787  0.976428
 1.01964   1.02133   1.02303   1.02476      0.97735   0.975934  0.974573
 1.01813   1.0198    1.0215    1.02321      0.975366  0.973954  0.972599
 1.01653   1.01819   1.01987   1.02157      0.97325   0.971849  0.970506
 ⋮                                       ⋱                      
 0.951463  0.95149   0.951528  0.951578     1.0403    1.03811   1.03584
 0.95317   0.953264  0.953368  0.953482  …  1.03951   1.03745   1.03531
 0.954834  0.954996  0.955165  0.955342     1.0387    1.03676   1.03474
 0.956457  0.956686  0.95692   0.95716      1.03788   1.03606   1.03414
 0.958035  0.958334  0.958635  0.958936     1.03705   1.03533   1.03352
 0.959567  0.959941  0.96031   0.960669     1.03621   1.0346    1.03289
 0.961047  0.961503  0.961939  0.96235   …  1.03537   1.03387   1.03224
 0.962455  0.963005  0.963497  0.963945     1.03452   1.03315   1.03159
 0.963738  0.964395  0.96491   0.965368     1.03364   1.03244   1.03091
v (in grid orientation)
104×159 Matrix{Float64}:
  0.0752184   0.0755549   0.075833   …  -0.0818985  -0.0806515  -0.0790153
  0.0770362   0.0772896   0.0775429     -0.0832243  -0.0821639  -0.0800897
  0.0788121   0.0790406   0.0792777     -0.0845858  -0.0835934  -0.0813074
  0.0805632   0.0807891   0.0810229     -0.0859737  -0.0849907  -0.0825866
  0.0822974   0.0825304   0.0827697     -0.0873806  -0.0863863  -0.0839
  0.0840187   0.0842636   0.0845142  …  -0.0888015  -0.0877898  -0.0852362
  0.0857289   0.0859886   0.0862544     -0.0902324  -0.0892017  -0.0865875
  0.0874281   0.0877046   0.0879885     -0.0916688  -0.0906184  -0.0879467
  0.0891155   0.0894103   0.0897146     -0.0931061  -0.0920353  -0.0893069
  0.090789    0.0911036   0.0914303     -0.0945392  -0.093447   -0.0906617
  ⋮                                  ⋱                          
 -0.0316417  -0.0320594  -0.0325492  …   0.0359338   0.0372454   0.0395536
 -0.031547   -0.0319554  -0.0324353      0.0339406   0.03519     0.0374097
 -0.0313856  -0.031786   -0.032258       0.032066    0.0332539   0.0353837
 -0.0311566  -0.0315511  -0.0320186      0.0303021   0.03143     0.0334677
 -0.0308577  -0.0312498  -0.0317179      0.0286428   0.0297128   0.0316549
 -0.0304842  -0.030881   -0.0313584  …   0.0270831   0.0281006   0.0299387
 -0.0300284  -0.0304453  -0.0309466      0.0256185   0.0265982   0.0283116
 -0.0294786  -0.0299525  -0.030499       0.0242424   0.0252262   0.0267573
 -0.0288208  -0.0294384  -0.0300405      0.0229449   0.0240227   0.0252288

In this last line, modes.a[end,:] obtains the expansion coefficients at the last time available.

Let's print the first mode, and the corresponding history of the modal coefficient in the decomposition

plot(layout=[2;1],plot(modes.phi[1].u,sys,title="u"),
   plot(modes.phi[1].v,sys,title="v"),
   plot(tpod,modes.a[:,1],xlim=(0,Inf),xlabel="\$t\$",ylabel="\$a_1(t)\$"))
Example block output

The energy associated with this mode is

modes.lambda[1]
0.719789956020543

Now let's print the $r$th mode, and the history of the coefficient in the decomposition

plot(layout=[2;1],plot(modes.phi[end].u,sys,title="u"),
   plot(modes.phi[end].v,sys,title="v"),
   plot(tpod,modes.a[:,end],xlim=(0,Inf),xlabel="\$t\$",ylabel="\$a_r(t)\$"))
Example block output

The energy associated with this mode is

modes.lambda[end]
0.006852880686168134

POD functions

ILMPostProcessing.pod — Function
pod(X::Vector{T}[; tolerance=0.99])

Calculate the POD modes and associated time-varying coefficients from an array of snapshot data X. This X plays the role of a snapshot matrix, whose columns are snapshots of the data. However, it is actually to be stored as a type Vector{T} where T<:GridData. It can be generated with a function like velocity(sol,sys,trange), where sol is a ODESolution, sys is an ILMSystem, and trange is an array of times, e.g., trange=range(1,10,length=100). The number of POD modes retained in the decomposition is set by tolerance: this specifies the fraction of the total energy to keep, and defaults to 99 percent.

The output of PODModes is a structure with the following fields

  • Xmean: temporal mean of the data. type T
  • Xnorm: original X vector with mean removed. Each element is of type T
  • phi: vector of POD modes. Each element is of type T
  • a: matrix of POD coefficients. Number of columns is same as number of entries in phi. Column k constitutes the time-varying coefficient for mode k in phi.
  • lambda: vector of modal energies, arranged in decreasing order, corresponding to the modes in phi
  • psi: matrix of
source
pod(X::Vector{T},r::Int)

Perform POD on snapshot data X and truncate to r modes

source

This page was generated using Literate.jl.