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);

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.14194    0.0152349   -0.0800551  …  -0.0258886    -0.0152696
  0.140543   0.0153246   -0.0790431     -0.0180305    -0.00952855
  0.138571   0.0156234   -0.0783249     -0.00939853   -0.0037204
  0.136239   0.0164329   -0.077899      -0.000959092   0.0014893
  0.133733   0.0179383   -0.0775806      0.00635362    0.00564145
  0.131235   0.0202009   -0.077098   …   0.0119154     0.0085169
  0.128913   0.0231614   -0.076146       0.0154992     0.0100634
  0.126906   0.0266551   -0.0744303      0.0172009     0.0102969
  0.125316   0.0304381   -0.0717035      0.0172914     0.00925677
  0.124198   0.0342168   -0.067791       0.0160731     0.00703318
  ⋮                                  ⋱                
 -0.112304  -0.0344556   -0.0532147     -0.0120673     0.00853449
 -0.114349  -0.0217861   -0.0588577     -0.0090689     0.0097446
 -0.115442  -0.00828272  -0.0625237     -0.00497924    0.00960382
 -0.115547   0.00579294  -0.0640886  …  -0.00014938    0.00812758
 -0.11465    0.0201689   -0.0634731      0.00501764    0.00548013
 -0.112757   0.0345678   -0.0606473      0.0101173     0.00194107
 -0.109901   0.048711    -0.0556342      0.0147904    -0.0021433
 -0.106135   0.0623227   -0.0485125      0.0187581    -0.00641432
 -0.101538   0.0751335   -0.0394171  …   0.0218405    -0.0105507

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, 154, 104, Float64, Vector{Float64}} data
u (in grid orientation)
103×154 Matrix{Float64}:
 1.032     1.03397   1.03586   1.03773   …  0.994236  0.992556  0.991314
 1.03115   1.03309   1.03501   1.03691      0.993158  0.991552  0.990161
 1.03018   1.0321    1.03403   1.03596      0.991975  0.990395  0.988936
 1.02911   1.03104   1.03297   1.03491      0.990679  0.989092  0.987592
 1.02797   1.02989   1.03182   1.03377      0.989258  0.987656  0.986126
 1.02673   1.02866   1.0306    1.03255   …  0.987707  0.986089  0.984537
 1.02541   1.02734   1.02928   1.03124      0.986024  0.984393  0.982825
 1.024     1.02592   1.02787   1.02983      0.984205  0.982566  0.98099
 1.02249   1.02441   1.02635   1.02832      0.98225   0.980608  0.979029
 1.02088   1.0228    1.02474   1.0267       0.980157  0.978516  0.976941
 ⋮                                       ⋱                      
 0.950972  0.951146  0.951332  0.951533     1.03813   1.03596   1.03368
 0.952807  0.953047  0.953298  0.953563  …  1.03729   1.03525   1.0331
 0.954589  0.954895  0.955211  0.955538     1.03645   1.03453   1.0325
 0.956318  0.956691  0.957072  0.95746      1.0356    1.0338    1.03187
 0.957994  0.958436  0.958882  0.95933      1.03475   1.03305   1.03123
 0.959616  0.960131  0.960643  0.961148     1.0339    1.0323    1.03058
 0.961176  0.961772  0.96235   0.962904  …  1.03303   1.03154   1.02991
 0.962655  0.963345  0.963978  0.964566     1.03215   1.0308    1.02924
 0.963999  0.964798  0.965449  0.966044     1.03125   1.03007   1.02854
v (in grid orientation)
104×153 Matrix{Float64}:
  0.079378    0.0796791   0.0799135  …  -0.086456   -0.0852579  -0.0836666
  0.0813472   0.0815723   0.0817851     -0.0879538  -0.0869376  -0.0849086
  0.0832878   0.0834915   0.083691      -0.0894933  -0.0885432  -0.0862995
  0.0852161   0.0854194   0.0856175     -0.091066   -0.0901241  -0.0877583
  0.0871396   0.0873514   0.0875567     -0.0926647  -0.0917105  -0.0892583
  0.089062    0.0892871   0.0895049  …  -0.0942849  -0.0933124  -0.0907882
  0.0909849   0.0912261   0.0914604     -0.0959226  -0.09493    -0.0923406
  0.0929083   0.0931678   0.0934215     -0.0975739  -0.0965605  -0.0939084
  0.0948313   0.0951108   0.0953864     -0.0992345  -0.0981994  -0.0954849
  0.0967519   0.097053    0.0973528     -0.1009     -0.0998417  -0.0970637
  ⋮                                  ⋱   ⋮                      
 -0.0351887  -0.0356327  -0.0361486  …   0.0347907   0.0362171   0.0387756
 -0.0349488  -0.035381   -0.0358838      0.0328253   0.0341775   0.0366264
 -0.0346428  -0.0350646  -0.0355569      0.0309782   0.0322575   0.0345954
 -0.03427    -0.0346837  -0.035169       0.0292412   0.0304494   0.0326741
 -0.0338283  -0.0342378  -0.0347216      0.0276078   0.0287478   0.0308556
 -0.0333131  -0.0337257  -0.0342172  …   0.0260736   0.0271513   0.0291339
 -0.0327164  -0.0331483  -0.0336627      0.0246344   0.0256662   0.027502
 -0.0320263  -0.0325156  -0.033075       0.0232839   0.0243153   0.0259442
 -0.0312274  -0.0318645  -0.0324798      0.0220126   0.0231401   0.0244124

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)\$"))

The energy associated with this mode is

modes.lambda[1]

0.7533616864230657

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)\$"))

The energy associated with this mode is

modes.lambda[end]

0.006931914393839419

This page was generated using Literate.jl.

Proper orthogonal decomposition (POD)

Get the flow field data

Assemble snapshots of the velocity field from the solution data

Perform the POD