diff --git a/flow_over_cylinder_re200/compute_shedding_frequency.jl b/flow_over_cylinder_re200/compute_shedding_frequency.jl new file mode 100644 index 0000000..999e15d --- /dev/null +++ b/flow_over_cylinder_re200/compute_shedding_frequency.jl @@ -0,0 +1,62 @@ +# Load necessary libraries +using FFTW, Plots, NPZ, Statistics, LaTeXStrings +include("plot_functions.jl") + +# Function to compute the vortex shedding frequency using FFT +function vortex_shedding_frequency(data, dt) + # Perform FFT on the time-series data + N = length(data) # Number of data points + freq_data = fft(data) + + # Compute the corresponding frequencies + freqs = (0:N-1) ./ (N * dt) # Frequency axis + + # Compute the power spectrum (magnitude of FFT) + power_spectrum = abs.(freq_data[1:div(N,2)]) # Keep positive frequencies only + freqs = freqs[1:div(N,2)] # Positive frequencies + + # Find the index of the maximum peak in the power spectrum + peak_idx = argmax(power_spectrum) + shedding_freq = freqs[peak_idx] + + # Plot power spectrum for visualization + gr(size=(570,450), xtickfontsize=12, ytickfontsize=12, xguidefontsize=20, yguidefontsize=14, legendfontsize=12, + dpi=300, grid=(:y, :gray, :solid, 1, 0.4), palette=cgrad(:plasma, 3, categorical = true)); + plt = plot(freqs, power_spectrum/maximum(power_spectrum), xlabel="Frequency (Hz)", xlims=(0, 1), ylabel="Normalized power", title="Power Spectrum of Vortex Shedding") + display(plt) + return shedding_freq +end + +# Example usage: +# Assume `vorticity_data` is a time-series array of vorticity or velocity measured at a point downstream of the cylinder +# and `dt` is the time step of the data. + +# Generate sample data (replace with actual vorticity/velocity data) +# vort_path = "/Users/l352947/mori_zwanzig/mori_zwanzig/mzmd_code_release/data/vort_all_re600_t5100.npy" +re_num = 200; +vort_path = "../data/vort_all_re200_t145.npy" +X = npzread(vort_path); +X = X .- mean(X, dims=2); + +T = size(X, 2); +t_skip = 1; +X = X[:, 1:t_skip:end] .- mean(X, dims=2); +dt = t_skip*0.2; #physical time dt +t = 0:dt:((T-1)*dt); # Time vector + +Y = copy(X); +F_select = reshape(Y, (199, 449, T))[100, 200, :]; #select a point in wake to compute psd + +# simple sanity check +# f_shedding = 2.0 # Example shedding frequency (Hz) for demonstration +# vorticity_data = sin.(2 * π * f_shedding .* t) # Simulated vorticity signal with shedding frequency + +vorticity_data = F_select[1:end]; +plt = plot(vorticity_data[1:100]) +display(plt) + +# Call the function to compute the shedding frequency +shedding_frequency = vortex_shedding_frequency(vorticity_data, dt) +println("Estimated vortex shedding frequency: $shedding_frequency Hz") + + diff --git a/flow_over_cylinder_re200/convergence_of_modes.jl b/flow_over_cylinder_re200/convergence_of_modes.jl new file mode 100644 index 0000000..5bdab55 --- /dev/null +++ b/flow_over_cylinder_re200/convergence_of_modes.jl @@ -0,0 +1,263 @@ +using Plots, LaTeXStrings, LinearAlgebra, NPZ +using Statistics +using Colors, ColorSchemes + + +""" +Comparing MZMD errors + - Measuring convergence rates of modes +""" + +include("./src/main_mz_algorithm.jl") +include("./src/compute_obervability_amplitudes_optimized.jl") +include("./plot_functions.jl") +include("./src/mzmd_modes.jl") + +# include("./hodmd_rtilde_comp.jl") +re_num = 600 +method="companion" +a_method = "x0"; +sample = 1; +println("running e.method $(method) and amp.method $(a_method) eigendecomp computing errors over r, k, d") + + +#--------- Load data +t_skip = 2; +dt = t_skip*0.2; +if re_num==600 + vort_path = "/Users/l352947/mori_zwanzig/modal_analysis_mzmd/mori_zwanzig/mzmd_code_release/data/vort_all_re600_t5100.npy" + # vort_path = "/Users/l352947/mori_zwanzig/mori_zwanzig/mzmd_code_release/data/vort_all_re600_t5100.npy" + X = npzread(vort_path)[:, 300:5000]; +else + vort_path = "/Users/l352947/mori_zwanzig/cylinder_mzmd_analysis/data/vort_all_re200.npy"; + X = npzread(vort_path) #T=1500 +end +X_mean = mean(X, dims=2); +X = X .- X_mean; +# X = X .+ 0.2*maximum(X) * randn(size(X)) + + +m = size(X, 1); +t_end = 4000 #total n steps is 751 +t_pred = 100 +T_train = ((sample):(t_end+sample-1))[1:t_skip:end]; +T_test = ((t_end+sample):(t_end + t_pred + sample))[1:t_skip:end]; +time_train = dt * T_train +time_test = dt * t_skip * T_test; +t_all = 0:dt:(dt*(t_end + t_pred +1)); +X_train = X[:,T_train]; +X_test = X[:,T_test]; +X_train = X_train .+ 0.15*maximum(X_train) .* randn(size(X_train)); + +n_ks = 8; +#number of snapshots used in fitting +t_win = size(T_train, 1) - n_ks - 1; +r = 20; +# num_modes = r*n_ks; #all +num_modes = 6; + + + +#-------- Compute MZ on POD basis observables +X1 = X_train[:, 1:t_win]; #Used to obtain Ur on 1:twin +X0 = X_train[:, 1]; + +# set this on first run +U, Sigma, V = svd(X_train); +sigmas_ = diagm(Sigma); +X_proj = sigmas_[1:r,1:r] * V[:,1:r]'; +Ur = U[:, 1:r]; + +function select_dominant_modes(a, Ur, Vc, r, n_ks, lam, num_modes) + Φ_mz = Ur*Vc[1:r,:]; + nmodes = size(Φ_mz, 2); + amp = zeros(nmodes); + for i in 1:nmodes + amp[i] = norm(Φ_mz[:,i]*a[i]); + end + #sort according to largest amplitude: + nmodes2 = minimum([size(Φ_mz, 2), num_modes]); + ind = sortperm(abs.(amp), rev=true) + a_dom = a[ind][1:nmodes2]; + Vc_dom = Vc[1:r, ind][:, 1:nmodes2]; + lam_dom = lam[ind][1:nmodes2]; + return a_dom, Vc_dom, lam_dom +end + + +function convergence_mzmd_modes(t_range, n_ks) + n_t = length(t_range) + modes_t = zeros(n_t, r, num_modes); + modesdmd_t = zeros(n_t, r, r); + function obtain_mz_operators(X_proj, t_win, n_ks) + #compute two time covariance matrices + Cov = obtain_C(X_proj, t_win, n_ks); + #compute MZ operators with Mori projection + M, Ω = obtain_ker(Cov, n_ks); + return Ω + end + t_idx = 1; + for t in t_range + T_train = 1:t + t_win = size(T_train, 1) - n_ks - 1; + X_proj_ = X_proj[:, T_train]; + println("t = ", t) + Ω = obtain_mz_operators(X_proj_, t_win, n_ks); + C = form_companion(Ω, r, n_ks); + Λc, Vc = eigen(C); + Z0 = Ur'*X_test[:, 1:(n_ks+1)]; + z0 = zeros(n_ks*r) + for k in 1:n_ks + # z0[((k-1)*r+1):k*r] = Z0[:,k] + z0[((k-1)*r+1):k*r] = Z0[:,(n_ks - k + 1)] + end + a = pinv(Vc)*z0; + Λc, a, Vc[1:r, :], Ur + a_dom, Vc_dom, lam_dom = select_dominant_modes(a, Ur, Vc[1:r,:], r, n_ks, Λc, num_modes); + nmodes = minimum([size(Vc_dom, 2), num_modes]) + for m in 1:nmodes + modes_t[t_idx, :, m] = abs.(a_dom[m].*Vc_dom[:, m]); + end + lam_dmd, phi_dmd = eigen(Ω[1,:,:]); + x0 = Ur'*X_test[:, 1]; + admd = pinv(phi_dmd)*x0; + admd_dom, Vdmd_dom, lamdmd_dom = select_dominant_modes(admd, Ur, phi_dmd, r, n_ks, Λc, r); + for m in 1:r + modesdmd_t[t_idx, :, m] = abs.(admd_dom[m].*Vdmd_dom[:, m]); + end + t_idx+=1; + end + return modes_t, modesdmd_t +end + + +function hodmd_dominant_modes(d, U, Sigma, T, J, K, T_train, r1, r2; subtractmean::Bool = false) + # %% STEP 1: SVD of the original data + sigmas=diagm(Sigma); #create diagonal matrix from vector + n = size(Sigma, 1); #number of singular values + NormS=norm(sigmas,2); + # %% Spatial complexity: kk + kk = r1 + U=U[:,1:kk]; + # %% reduced snapshots matrix + hatT=sigmas[1:kk,1:kk]*T[T_train,1:kk]'; + N=size(hatT, 1); + # %% Create the modified snapshot matrix + tildeT=zeros(d*N, K-d+1); + for ppp=1:d + tildeT[(ppp-1)*N+1:ppp*N,:]=hatT[:,ppp:ppp+K-d]; + end + # %% Dimension reduction 2 + U1,Sigma1,T1 = svd(tildeT); + sigmas1=diagm(Sigma1); + n = size(sigmas1, 1); + NormS1=norm(sigmas1,2); + # Second Spatial dimension reduction + kk1 = r2 #d*r1 + if r2 > size(U1, 2) + kk1 = size(U1, 2); #this will break companion structure, but is a limit of this method! + end + U1=U1[:,1:kk1]; + hatT1=sigmas1[1:kk1,1:kk1]*T1[:,1:kk1]'; + # %% Reduced modified snapshot matrix (SVD Reduction 3) + K1 = size(hatT1, 2); + tildeU1,tildeSigma,tildeU2 =svd(hatT1[:,1:K1-1]); + tildesigmas = diagm(tildeSigma) + tildeR=hatT1[:,2:K1]*tildeU2*inv(tildesigmas)*tildeU1'; + lam, e_vects = eigen(tildeR); + Q=U1*e_vects; + a, Q = compute_amplitudes_given_ic_z0_hodmd(Ur, X_test[:, 1:n_ks], Q, n_ks, r) + a_dom, Vc_dom, lam_dom = select_dominant_modes(a, Ur, Q, r, 1, lam, num_modes); + return a_dom, Vc_dom, lam_dom +end + +function convergence_hodmd_dom_modes(t_range, d, r) + n_t = length(t_range) + modes = zeros(n_t, r, num_modes); + t_idx = 1; + for t in t_range + T_train = 1:t + X_train = X[:,T_train]; + J,K = size(X_train); + println("t = ", t) + a_dom, phi, _ = hodmd_dominant_modes(d, U, Sigma, V, J, K, T_train, r, d*r); + nmodes = minimum([size(phi, 2), num_modes]) + for m in 1:nmodes + modes[t_idx, :, m] = abs.(a_dom[m].*phi[:, m]); + end + t_idx+=1; + end + return modes +end + +function compute_convergence_mzmd_modes(modes_t) + nts = size(modes_t, 1); + conv = zeros(nts-1); + for i in 1:(nts-1) + conv[i] = norm(modes_t[i+1,:,:] .- modes_t[i,:,:])/norm(modes_t[2,:,:]); + end + return conv +end + +function compute_convergence_hodmd_modes(modes) + nts = size(modes, 1); + conv = zeros(nts-1); + for i in 1:(nts-1) + conv[i] = norm(modes[i+1,:,:] - modes[i,:,:])/norm(modes[1,:,:]); + end + return conv +end + +t_range = n_ks*r:10:1800 +# t_range = r:10:1500 +t_range_mz = r:10:1800 + +modes_t, modesdmd_t = convergence_mzmd_modes(t_range_mz, n_ks); +modes = convergence_hodmd_dom_modes(t_range, n_ks, r) + +conv_mz = compute_convergence_mzmd_modes(modes_t) +conv_dmd = compute_convergence_mzmd_modes(modesdmd_t) +conv_hodmd = compute_convergence_hodmd_modes(modes) + + +function make_directory(path::String) + if !isdir(path) # Check if the directory already exists + mkdir(path) + println("Directory created at: $path") + else + println("Directory already exists at: $path") + end +end + +# dir_path = "convergence_modes_noise_dom_dt$(t_skip)_te$(t_end)_re$(re_num)_nm$(num_modes)" +# make_directory(dir_path) + +# npzwrite("$(dir_path)/conv_dom_modes_mz_r$(r)_k$(n_ks).npy", conv_mz) +# npzwrite("$(dir_path)/conv_dom_modes_hodmd_r$(r)_k$(n_ks).npy", conv_hodmd) +# npzwrite("$(dir_path)/conv_dom_modes_dmd_r$(r)_k$(n_ks).npy", conv_dmd) + + +function plot_convergence(t_range, t_rangemz, mz_oms, hodmd_R, dmd_conv) + gr(size=(570,450), xtickfontsize=12, ytickfontsize=12, xguidefontsize=20, yguidefontsize=14, legendfontsize=12, + dpi=300, grid=(:y, :gray, :solid, 1, 0.4), palette=cgrad(:plasma, 3, categorical = true)); + plt = plot(t_range, hodmd_R, fillalpha=.4, ms = 7.0, linestyle=:dash, linewidth=1.5, color="black", + yaxis=:log, legendfontsize=12, left_margin=2mm, bottom_margin=2mm, label=L"\textbf{HODMD}") + plot!(t_rangemz, mz_oms, fillalpha=.4, ms = 7.0, linewidth=1.5, color="blue", + yaxis=:log, legendfontsize=12, left_margin=2mm, bottom_margin=2mm, label=L"\textbf{MZMD}") + # plot!(t_rangemz, dmd_conv, fillalpha=.4, ms = 7.0, linestyle=:dashdot, linewidth=1.5, color="green", + # yaxis=:log, legendfontsize=15, left_margin=2mm, bottom_margin=2mm, label=L"\textbf{DMD}") + if num_modes==r*n_ks + title!(L"\textrm{Convergence ~ of ~ modes}", titlefont=20) + else + title!(L"\textrm{Convergence ~ of ~ dominant ~ modes}", titlefont=20) + end + xlabel!(L"\textrm{Number ~ of ~ Snapshots}", xtickfontsize=14, xguidefontsize=16) + ylabel!(L"||\Phi_n - \Phi_{n-1}||/||\Phi_1||", ytickfontsize=14, yguidefontsize=16) + savefig(plt, "./figures/convergence_of_modes_noise_r$(r)_k$(n_ks).png") + return plt +end + + +plt = plot_convergence(t_range[1:(size(conv_hodmd, 1)-2)], t_range_mz[1:(size(conv_mz, 1)-2)], conv_mz[2:end-1], conv_hodmd[1:end-2], conv_dmd[1:end-2]) +display(plt) + diff --git a/flow_over_cylinder_re200/figures/pod_mode_energy_ts1_r20_k6.png b/flow_over_cylinder_re200/figures/pod_mode_energy_ts1_r20_k6.png new file mode 100644 index 0000000..329287e Binary files /dev/null and b/flow_over_cylinder_re200/figures/pod_mode_energy_ts1_r20_k6.png differ diff --git a/flow_over_cylinder_re200/figures/pod_mode_energy_ts2_r20_k6.png b/flow_over_cylinder_re200/figures/pod_mode_energy_ts2_r20_k6.png new file mode 100644 index 0000000..0995512 Binary files /dev/null and b/flow_over_cylinder_re200/figures/pod_mode_energy_ts2_r20_k6.png differ diff --git a/flow_over_cylinder_re200/figures/pod_mode_energy_ts2_r5_k6.png b/flow_over_cylinder_re200/figures/pod_mode_energy_ts2_r5_k6.png new file mode 100644 index 0000000..751f367 Binary files /dev/null and b/flow_over_cylinder_re200/figures/pod_mode_energy_ts2_r5_k6.png differ diff --git a/flow_over_cylinder_re200/gen_x0_errors_sweep_k_dominant_modes_pred.jl b/flow_over_cylinder_re200/gen_x0_errors_sweep_k_dominant_modes_pred.jl new file mode 100644 index 0000000..de3d383 --- /dev/null +++ b/flow_over_cylinder_re200/gen_x0_errors_sweep_k_dominant_modes_pred.jl @@ -0,0 +1,210 @@ +using Plots, LaTeXStrings, LinearAlgebra, NPZ +using Statistics +using Colors, ColorSchemes + + +""" +Comparing MZMD errors + - Measuring reconstrucion and validation errors +""" + +include("./src/main_mz_algorithm.jl") +include("./src/compute_obervability_amplitudes_optimized.jl") +include("./plot_functions.jl") +include("./src/mzmd_modes.jl") + + +re_num = 600; +# re_num = 200; +method="companion" +a_method = "x0" +sample = parse(Int, ARGS[1]); +println(" Sample = ", sample) +println("running e.method $(method) and amp.method $(a_method) eigendecomp computing errors over r, k, d") + + +#--------- Load data +t_skip = 2; +nmodes_ = 6; +dt = t_skip*0.2; +if re_num==600 + # vort_path = "/Users/l352947/mori_zwanzig/mori_zwanzig/mzmd_code_release/data/vort_all_re600_t5100.npy" + vort_path = "/Users/l352947/mori_zwanzig/modal_analysis_mzmd/mori_zwanzig/mzmd_code_release/data/vort_all_re600_t5100.npy" + X = npzread(vort_path)[:, 300:5100]; +else + vort_path = "/Users/l352947/mori_zwanzig/cylinder_mzmd_analysis/data/vort_all_re200.npy"; + X = npzread(vort_path)[:, 200:end] #T=1500 +end + +X_mean = mean(X, dims=2); +X = X .- X_mean; +#check robustness to noise +m = size(X,1); + +# t_end = 3000 #total n steps is 751 +# t_pred = 220 +t_end = 200 #total n steps is 751 +t_pred = 100 + +T_train = ((1):(t_end+1-1))[1:t_skip:end]; +T_test = ((t_end+sample):(t_end + t_pred + sample))[1:t_skip:end]; +time_train = dt * T_train +time_test = dt * t_skip * T_test; +t_all = 0:dt:(dt*(t_end + t_pred +1)); +X_train = X[:,T_train]; +X_train = X_train .+ 0.15*maximum(X_train) .* randn(size(X_train)) +X_test = X[:,T_test]; + +t_end = round(Int, t_end/t_skip); +t_pred = round(Int, t_pred/t_skip); +println("t end = ", t_end); +U, S, V = svd(X_train); + +function mzmd_modes_reduced_amps_for_errors(twin, X, r, n_ks) + #compute modes and amplitudes from C + function obtain_mz_operators_(G, t_win, n_ks) + Cov = obtain_C(G, t_win, n_ks); + _, Ω = obtain_ker(Cov, n_ks); + println("----------------------------"); + return Ω + end + Ur = U[:, 1:r]; + X_proj = Ur' * X[:, T_train]; + # if n_ks==1 #hodmd (i.e. no mz memory) + # # #TODO SVD 2 + # t_g = twin + n_ks + 1; + # t_gtilde = length(1:(1+t_g - 1)); + # G_tilde=zeros(r, t_gtilde); + # for i in 1:1 + # G_tilde[(i-1)*r+1:i*r,:] = X_proj[:,i:(i+t_g-1)]; + # end + # U1, Sigma1, T1 = svd(G_tilde); + # sigmas1=diagm(Sigma1); + # U1=U1[:,1:r]; + # #svd low rank approximation of G_tilde: + # hatG=sigmas1[1:r,1:r]*T1[:,1:r]'; + # # #TODO SVD 3 + # K1 = size(hatG, 2); + # tildeU1, tildeSigma, tildeV1 = svd(hatG[:,1:K1-1]); + # tildesigmas = diagm(tildeSigma) + # Ω_tilde=hatG[:,2:K1]*tildeV1*inv(tildesigmas)*tildeU1'; + # Λc, Vc = eigen(Ω_tilde); + # Vc=U1*Vc; + # a, Vc = compute_amplitudes_given_ic_z0_hodmd(Ur, X_test[:, 1], Vc, 1, r) + # else + Ω = obtain_mz_operators_(X_proj, twin, n_ks); + C = form_companion(Ω, r, n_ks); + Λc, Vc = eigen(C); + Z0 = Ur'*X_test[:, 1:(n_ks+1)]; + z0 = zeros(n_ks*r) + for k in 1:n_ks + # z0[((k-1)*r+1):k*r] = Z0[:,k] + z0[((k-1)*r+1):k*r] = Z0[:,(n_ks - k + 1)] + end + a = pinv(Vc)*z0; + return Λc, a, Vc[1:r, :], Ur +end + +function pointwise_mse_errors(t1, gt_data, Xmz_pred) + mz_pred = (Xmz_pred .+ X_mean)[:, t1]; + # mz_pred = Xmz_pred[:, t1]; + # gt_data = X_test[:, t2]; + mz_err = mean(mean((mz_pred .- gt_data).^2, dims=2))/mean(gt_data.^2); + max_mz = maximum(mz_pred .- gt_data)/maximum(gt_data) + return mz_err, max_mz +end + +function simulate_modes_time_delay(Ur, evects, a, r, λ, dt, t) + a_length = length(a); + t_length = length(t); + time_dynamics = complex(zeros(a_length, t_length)); + ts = (0:t_length-1) * dt; + omega = log.(λ)/dt; + for i in 1 : t_length + time_dynamics[:,i] = (a.*exp.(ts[i] .* omega)); + end + X_pred = evects * time_dynamics; #in higher dim delay space + Xmz_pred = Ur * X_pred[1:r,:]; + return real.(Xmz_pred) +end + +function select_dominant_modes(a, Ur, Vc, r, n_ks, lam, num_modes) + Φ_mz = Ur*Vc[1:r,:]; + amp = zeros(r*n_ks); + for i in 1:r*n_ks + amp[i] = norm(Φ_mz[:,i]*a[i]); + end + #sort according to largest amplitude: + ind = sortperm(abs.(amp), rev=true) + a_dom = a[ind][1:num_modes]; + Vc_dom = Vc[1:r, ind][:, 1:num_modes]; + lam_dom = lam[ind][1:num_modes]; + return a_dom, Vc_dom, lam_dom +end + + +function compute_error_r_nk(nk_range, nd_range, r_range, nk_end) + t_pred = length(T_test) + err_mzmd = zeros(t_pred, size(nk_range,1), size(nd_range,1), size(r_range,1)); + err_mzmd_max = zeros(t_pred, size(nk_range,1), size(nd_range,1), size(r_range,1)); + err_recon_mzmd = zeros(size(nk_range,1), size(nd_range,1), size(r_range,1)); + err_recon_max = zeros(size(nk_range,1), size(nd_range,1), size(r_range,1)); + r_idx = 1; + for r in r_range + nd_idx = 1; + for d in nd_range + nk_idx = 1; + for k in nk_range + t_win = size(T_train, 1) - k - 1; + # lam, Q, a, Ur = compute_modes_amplitudes_time_delay_observables(X, r, k, d, t_win, method); + lam, a, Q, Ur = mzmd_modes_reduced_amps_for_errors(t_win, X, r, k) + num_modes = minimum([k*r, nmodes_]); + a_dom, Vc_dom, lam_dom = select_dominant_modes(a, Ur, Q, r, k, lam, num_modes) #dominant three modes + solution_modes = simulate_modes_time_delay(Ur, Vc_dom, a_dom, r, lam_dom, dt, time_test) + for t in 1:(t_pred-nk_end) + err_mzmd[t, nk_idx, nd_idx, r_idx], err_mzmd_max[t, nk_idx, nd_idx, r_idx] = + pointwise_mse_errors(t, X_test[:, t+k-1] .+ X_mean, solution_modes) + end + println(" r = ", r, " nk = ", k, " nd = ", d, " mz err = ", mean(err_mzmd[:, nk_idx, nd_idx, r_idx]), + " mz err max = ", mean(err_mzmd_max[:, nk_idx, nd_idx, r_idx]), " recon err = ", err_recon_mzmd[nk_idx, nd_idx, r_idx]) + nk_idx += 1; + end + nd_idx += 1; + end + r_idx += 1; + end + return err_mzmd, err_mzmd_max, err_recon_mzmd, err_recon_max +end + +# r_range = cat(3:2:9); +r_range = 3:2:9; +nk_range = cat(1:2, 4:2:24, dims=1) +nd_range = [1]; +nk_end = nk_range[end]; + +err_mzmd, err_mzmd_max, err_recon_mzmd, err_recon_max = compute_error_r_nk(nk_range, nd_range, r_range, nk_end); + +nr = length(r_range); +nk = length(nk_range); +nd = length(nd_range); +nk_end = nk_range[end]; +nd_end = nd_range[end]; + +t_pred_length = t_pred; + + +function make_directory(path::String) + if !isdir(path) # Check if the directory already exists + mkdir(path) + println("Directory created at: $path") + else + println("Directory already exists at: $path") + end +end + +# dir_path = "noise_dom_modes_gen_error_x0_test_dt$(t_skip)_te$(t_end)_re$(re_num)_nm$(nmodes_)" +dir_path = "noise_dom_modes_gen_error_x0_test_dt$(t_skip)_te$(t_end)_re$(re_num)_nm$(nmodes_)" +make_directory(dir_path) + +npzwrite("$(dir_path)/cyl_ax0_err_mzmd_delay_k1dmd_dts$(t_skip)_st$(sample)_$(method)_$(a_method)_re$(re_num)_modes_nr$(nr)_nk$(nk)_tp$(t_pred_length)_te$(t_end)_nke$(nk_end)_nde$(nd_end).npy", err_mzmd) +npzwrite("$(dir_path)/cyl_ax0_max_err_mzmd_delay_k1dmd_dts$(t_skip)_st$(sample)_$(method)_$(a_method)_re$(re_num)_modes_nr$(nr)_nk$(nk)_tp$(t_pred_length)_te$(t_end)_nke$(nk_end)_nde$(nd_end).npy", err_mzmd_max) diff --git a/flow_over_cylinder_re200/plot_functions.jl b/flow_over_cylinder_re200/plot_functions.jl new file mode 100644 index 0000000..ea86628 --- /dev/null +++ b/flow_over_cylinder_re200/plot_functions.jl @@ -0,0 +1,471 @@ +using Measures + +function circle_shape(h,k,r) + θ=LinRange(0,2*pi,500); + return h .+ r*sin.(θ), k .+ r*cos.(θ) +end + +function plot_field(X, t, method) + gsy = 600; gsx = round(Int, 2.25*gsy) + gr(size=(gsx,gsy), dpi=300) + t_dim = size(X, 2); + Y = copy(X); + Y[abs.(Y).>5] .= 5 + F = reshape(Y, (199, 449, t_dim)); + x_range = range(-1, 8, length=449); y_range = range(-2, 2, length=199); + + plt = plot(x_range, y_range, F[:,:, t]', st=:heatmap, + # plt = plot(x_range, y_range, F[:,:, t], st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + scatter!([x_range[200]], [y_range[100]], color="black", ms=4) + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1, clims=(-8.1, 8.1)) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + #title!(L"\textbf{Vorticity}", titlefont=22) + title!(L"\textrm{%$(method)}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + display(plt) + return plt +end + +function plot_field_comp(Xgt, Xpr, t) + gsy = 1250; gsx = round(Int, 1.2*gsy) + gr(size=(gsx,gsy)) + t_dim_pr = size(Xpr, 2); + t_dim_gt = size(Xgt, 2); + Ygt = copy(Xgt); + Ygt[abs.(Ygt).>4] .= 4 + Ypr = copy(Xpr); + Ypr[abs.(Ypr).>4] .= 4 + Fgt = reshape(Ygt, (199, 449, t_dim_gt)); + Fpr = reshape(Ypr, (199, 449, t_dim_pr)); + + x_range = range(-1, 8, length=449); y_range = range(-2, 2, length=199); + + plt1 = plot(x_range, y_range, Fgt[:,:, t]', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1, clims=(-4.1, 4.1)) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + title!(L"\textbf{Vorticity}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + + plt2 = plot(x_range, y_range, Fpr[:,:, t]', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1, clims=(-4.1, 4.1)) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + title!(L"\textbf{Vorticity ~ Prediction}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + + plt = plot(plt1, plt2, layout=(2,1)) + display(plt) +end + +function plot_field_diff_mse_mz_hodmd_dmd(Xgt, Xpr, Xdmd, t, c_max) + gsy = 1250; gsx = round(Int, 1.2*gsy) + gr(size=(gsx,gsy)) + t_dim = size(Xpr, 2); + Ygt = copy(Xgt); + Ypr = copy(Xpr); + Ydmd = copy(Xdmd); + Fgt = reshape(Ygt, (199, 449)); + Fpr = reshape(Ypr, (199, 449)); + Fdmd = reshape(Ydmd, (199, 449)); + + x_range = range(-1, 8, length=449); y_range = range(-2, 2, length=199); + + plt1 = plot(x_range, y_range, Fgt[:,:]', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance)), clims=(0, c_max)) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + title!(L"\textbf{Vorticity ~ MZMD ~ MSE}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + + plt2 = plot(x_range, y_range, Fpr[:,:]', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance)), clims=(0, c_max)) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + title!(L"\textbf{Vorticity ~ HODMD ~ MSE}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + + #DMD + plt3 = plot(x_range, y_range, Fdmd[:,:]', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance)), clims=(0, c_max)) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + title!(L"\textbf{Vorticity ~ DMD ~ MSE}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + + plt = plot(plt1, plt2, plt3, layout=(3,1)) + display(plt) + # return plt +end + + +function plot_field_diff_mse(Xgt, Xpr, t, c_max) + gsy = 1250; gsx = round(Int, 1.2*gsy) + gr(size=(gsx,gsy)) + t_dim = size(Xpr, 2); + Ygt = copy(Xgt); + Ypr = copy(Xpr); + Fgt = reshape(Ygt, (199, 449)); + Fpr = reshape(Ypr, (199, 449)); + + x_range = range(-1, 8, length=449); y_range = range(-2, 2, length=199); + + plt1 = plot(x_range, y_range, Fgt[:,:]', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance)), clims=(0, c_max)) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + title!(L"\textbf{Vorticity ~ MZMD ~ MSE}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + + plt2 = plot(x_range, y_range, Fpr[:,:]', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance)), clims=(0, c_max)) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + title!(L"\textbf{Vorticity ~ DMD ~ MSE}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + + plt = plot(plt1, plt2, layout=(2,1)) + display(plt) + # return plt +end + +function plot_field_diff(Xgt, Xpr, t) + gsy = 1250; gsx = round(Int, 1.2*gsy) + gr(size=(gsx,gsy)) + t_dim_gt = size(Xgt, 2); + t_dim_pr = size(Xpr, 2); + Ygt = copy(Xgt); + Ygt[abs.(Ygt).>4] .= 4 + Ypr = copy(Xpr); + Ypr[abs.(Ypr).>4] .= 4 + Fgt = reshape(Ygt, (199, 449, t_dim_gt)); + Fpr = reshape(Ypr, (199, 449, t_dim_pr)); + + x_range = range(-1, 8, length=449); y_range = range(-2, 2, length=199); + + plt1 = plot(x_range, y_range, Fgt[:,:, t]', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + title!(L"\textbf{Vorticity ~ MZMD ~ Difference}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + + plt2 = plot(x_range, y_range, Fpr[:,:, t]', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + title!(L"\textbf{Vorticity ~ DMD ~ Difference}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + + plt = plot(plt1, plt2, layout=(2,1)) + display(plt) + # return plt +end + + +function plotting_sing_val(r_max, S) + gr(size=(700,600)) + rs = 1 : r_max + plt = scatter(rs, S[1:r_max], yaxis =:log, label=false, + ms=7.5, legendfontsize=10, color="black") + title!(L"\textrm{Singular Values}", titlefont=22) + xlabel!(L"\textrm{r}", xtickfontsize=12, xguidefontsize=16) + ylabel!(L"\sigma_k", ytickfontsize=12, yguidefontsize=16) + display(plt) + # savefig(plt, "./figures/singular_values_r.png") +end + +function plot_mzmd_Mmodes(Φ, m) + gsy = 600; gsx = round(Int, 2.25*gsy) + gr(size=(gsx,gsy)); + nx = 199; ny = 449; + Φ = reshape(Φ, (199, 449, r)); + x_range = range(-1, 8, length=449); y_range = range(-2, 2, length=199); + + plt = plot(x_range, y_range, real.(Φ[:,:, m])', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + mode_n = ceil(Int, m/2) + title!(L"\textrm{Markovian ~ (DMD) ~ mode ~ } \Phi^{m}_{%$mode_n}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + # savefig(plt, "./figures/mzmd_Markovian_mode$(m)_nk$(n_ks)_r$(r).png") + return plt +end + +function plot_hodmd_modes(Φ, m) + gsy = 600; gsx = round(Int, 2.25*gsy) + gr(size=(gsx,gsy)); + nx = 199; ny = 449; + Φ = reshape(Φ, (nx, ny, r)); + r_2 = floor(Int, n_ks*r/2) + + x_range = range(-1, 8, length=ny); y_range = range(-2, 2, length=nx); + + plt = plot(x_range, y_range, real.(Φ[:,:, m])', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + mode_n = ceil(Int, m/2) + title!(L"\textrm{HODMD ~ mode ~ } \Phi^{mz}_{%$mode_n}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + # savefig(plt, "./figures/mzmd_Markovian_mode$(m)_nk$(n_ks)_r$(r).png") + return plt +end + + +function plot_energy_content(s,r,t_win, eps=0.99) + gr(size=(700,600)); + Λ = 1/t_win * s.^2; + energy_r = zeros(r); + for i in 1 : r + energy_r[i] = sum(Λ[1:i])/sum(Λ); + end + plt = scatter(energy_r, ms=5.0, label=false, ylims=(0.0, 1.1), + legendfontsize=12, color="black") + y_99 = eps * ones(r); + per_eps = 100*eps; + plot!(y_99, linewidth=3.0, label=L"%$per_eps \textbf{\% ~ variations}", ls=:dash, color="grey") + title!(L"\textbf{Total ~ variation ~ contained ~ in ~ POD ~ modes}", titlefont=20) + xlabel!(L"\textbf{r}", xtickfontsize=12, xguidefontsize=16) + ylabel!(L"\textbf{Total variation}", ytickfontsize=12, yguidefontsize=16) + # savefig(plt, "./figures/pod_mode_energy_ts$(t_skip)_nt$(nt_samp)_r$(r).png") + display(plt) + # return plt +end + + + +function plot_mzmd_proj_memory_modes(Φ, m) + gsy = 600; gsx = round(Int, 2.25*gsy) + gr(size=(gsx,gsy)); + nx = 199; ny = 449; + Φ = reshape(Φ, (199, 449, r)); + r_2 = floor(Int, n_ks*r/2) + + x_range = range(-1, 8, length=449); y_range = range(-2, 2, length=199); + + plt = plot(x_range, y_range, real.(Φ[:,:, m])', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + mode_n = ceil(Int, m/2) + title!(L"\textrm{MZMD ~ Memory ~ modes ~ } \Phi^{mz}_{%$mode_n}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + # savefig(plt, "./figures/mzmd_Markovian_mode$(m)_nk$(n_ks)_r$(r).png") + return plt +end + +function plot_mzmd_memory_modes(Φ, m) + gsy = 600; gsx = round(Int, 2.25*gsy) + gr(size=(gsx,gsy)); + nx = 199; ny = 449; + Φ = reshape(Φ, (199, 449, n_ks*r)); + r_2 = floor(Int, n_ks*r/2) + + x_range = range(-1, 8, length=449); y_range = range(-2, 2, length=199); + + plt = plot(x_range, y_range, real.(Φ[:,:, m])', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + mode_n = ceil(Int, m/2) + title!(L"\textrm{MZMD ~ Memory ~ modes ~ } \Phi^{mz}_{%$mode_n}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + # savefig(plt, "./figures/mzmd_Markovian_mode$(m)_nk$(n_ks)_r$(r).png") + return plt +end + +function plot_mzmd_modes(Φ, m) + gsy = 600; gsx = round(Int, 2.25*gsy) + gr(size=(gsx,gsy)); + nx = 199; ny = 449; + Φ = reshape(Φ, (199, 449, n_ks*r)); + r_2 = floor(Int, n_ks*r/2) + + x_range = range(-1, 8, length=449); y_range = range(-2, 2, length=199); + + plt = plot(x_range, y_range, real.(Φ[:,:, m])', st=:heatmap, + bottom_margin=6mm, left_margin=6mm, right_margin=6mm, + top_margin=6mm, fill=(true, cgrad(:balance))) + + plot!(circle_shape(0,0,.5), seriestype=[:shape], lw = 0.5, + c = :black, legend=false, fillalpha=0.6, aspectratio=1) + xlims!(-1.06, 8.06); ylims!(-2.06, 2.06); + + mode_n = ceil(Int, m/2) + title!(L"\textrm{MZMD ~ mode ~ } \Phi^{mz}_{%$mode_n}", titlefont=22) + xlabel!(L"\textbf{x}", xtickfontsize=12, xguidefontsize=20) + ylabel!(L"\textbf{y}", ytickfontsize=12, yguidefontsize=20) + # savefig(plt, "./figures/mzmd_Markovian_mode$(m)_nk$(n_ks)_r$(r).png") + return plt +end + +function plot_amplitude_vs_frequency_markovian(λ, b, k) + gr(size=(700,600)); + # r_2 = floor(Int, r/2) + r_2 = r + # plt = scatter(imag(λ)[r_2:end], abs.(b)[r_2:end], linewidth=2.2, legend=false, ms=6.5) + plt = scatter(imag(λ)[1:r_2], abs.(b)[1:r_2], linewidth=2.2, legend=false, ms=6.5) + title!(L"\textrm{ Amplitude ~ vs ~ Frequency ~ } \mathbf{\Omega}^{(0)}", titlefont=18) + xlabel!(L"\textrm{Im}(\lambda)", xtickfontsize=14, xguidefontsize=14) + ylabel!(L"a_n", ytickfontsize=14, yguidefontsize=14) + # savefig(plt, "./figures/markovian_cylinder_spectrum_r$(r)_k$(k).png") + return plt +end + + +function plot_amplitude_vs_frequency_mzmd(λ, a) + gr(size=(700,600)); + kr = size(λ,1); + r_2 = floor(Int, kr/2) + # plt = scatter(imag(λ)[r_2:end], abs.(a)[r_2:end], legend=false, ms=6.5) + # plt = scatter(imag(λ)[1:r_2], abs.(a)[1:r_2], legend=false, ms=6.5) + plt = scatter(imag(λ), abs.(a), legend=false, ms=6.5) + + title!(L"\textrm{ Amplitude ~ vs ~ Frequency ~ } C", titlefont=18) + xlabel!(L"\textrm{Im}(\lambda)", xtickfontsize=14, xguidefontsize=14) + ylabel!(L"a_n", ytickfontsize=14, yguidefontsize=14) + # savefig(plt, "./figures/markovian_cylinder_spectrum_r$(r)_k$(k).png") + return plt +end + +function plotting_mem_norm(Ω) + gr(size=(700,600)) + Ω_norms = zeros(n_ks) + for i in 1:n_ks + Ω_norms[i] = norm(Ω[i,:,:]) + end + y_vals = Ω_norms[1:end] / norm(M); + plt_omega = Plots.scatter(y_vals, yaxis=:log, label=false, ms=7.5, color="black") + title!(L"\textrm{Memory ~ Effects: } k = %$n_ks", titlefont=22) + xlabel!(L"k", xtickfontsize=12, xguidefontsize=16) + ylabel!(L"||\Omega^{(k)}||_F / ||\Omega^{(0)}||_F", ytickfontsize=12, yguidefontsize=16) + display(plt_omega) + # savefig(plt_omega, "./figures/mem_mz_norm_nd$(n_ks)_r$(r).png") +end + +function plot_c_spectrum(Λc) + gr(size=(570,570), xtickfontsize=12, ytickfontsize=12, xguidefontsize=20, yguidefontsize=14, legendfontsize=12, + dpi=300, grid=(:y, :gray, :solid, 1, 0.4), palette=cgrad(:plasma, 3, categorical = true)); + + # Λc, Vc = load_py_eigen() + function circle_shape(h,k,r) + θ=LinRange(0,2*pi,500); + return h .+ r*sin.(θ), k .+ r*cos.(θ) + end + plt = scatter(Λc, linewidth=2.0, ms = 4.0, color = "black", legend=false, grid = true, framestyle = :box) + plot!(circle_shape(0,0,1.0), linewidth=1.5, color="black", linestyle=:dash) + title!(L"\textrm{MZMD ~ Eigenvalues }", titlefont=22) + xlabel!(L"\textrm{Re}(\lambda)", xtickfontsize=12, xguidefontsize=16) + ylabel!(L"\textrm{Im}(\lambda)", ytickfontsize=12, yguidefontsize=16) + return plt +end + +function plot_dmd_spectrum(Λc) + gr(size=(570,570), xtickfontsize=12, ytickfontsize=12, xguidefontsize=20, yguidefontsize=14, legendfontsize=12, + dpi=300, grid=(:y, :gray, :solid, 1, 0.4), palette=cgrad(:plasma, 3, categorical = true)); + + # Λc, Vc = load_py_eigen() + function circle_shape(h,k,r) + θ=LinRange(0,2*pi,500); + return h .+ r*sin.(θ), k .+ r*cos.(θ) + end + plt = scatter(Λc, linewidth=2.0, ms = 4.0, color = "black", legend=false, grid = true, framestyle = :box) + plot!(circle_shape(0,0,1.0), linewidth=1.5, color="black", linestyle=:dash) + title!(L"\textrm{DMD ~ Eigenvalues }", titlefont=22) + xlabel!(L"\textrm{Re}(\lambda)", xtickfontsize=12, xguidefontsize=16) + ylabel!(L"\textrm{Im}(\lambda)", ytickfontsize=12, yguidefontsize=16) + return plt +end + + +function plot_energy_content(s,r) + gr(size=(570,300), xtickfontsize=12, ytickfontsize=12, xguidefontsize=20, yguidefontsize=14, legendfontsize=12, + dpi=300, grid=(:y, :gray, :solid, 1, 0.4), palette=cgrad(:plasma, 3, categorical = true)); + + Λ = 1/t_win * s.^2; + energy_r = zeros(r); + for i in 1 : r + energy_r[i] = sum(Λ[1:i])/sum(Λ); + end + plt = plot(energy_r, lw=3.0, label=false, ylims=(0.0, 1.0), + legendfontsize=12, color="black") + y_99 = 0.99 * ones(r); + plot!(y_99, linewidth=3.0, label=L"\textbf{99\% ~ variations}", ls=:dash, color="grey") + title!(L"\textrm{Total ~ variation ~ in ~ POD ~ modes}", titlefont=22) + xlabel!(L"\textrm{r}", xtickfontsize=12, xguidefontsize=16) + ylabel!(L"\textrm{Voriticy ~ variation}", ytickfontsize=12, yguidefontsize=16) + savefig(plt, "./figures/pod_mode_energy_ts$(t_skip)_r$(r)_k$(n_ks).png") + return plt +end \ No newline at end of file diff --git a/flow_over_cylinder_re200/predictions_analysis_flared_cone.jl b/flow_over_cylinder_re200/predictions_analysis_flared_cone.jl new file mode 100644 index 0000000..dacffdb --- /dev/null +++ b/flow_over_cylinder_re200/predictions_analysis_flared_cone.jl @@ -0,0 +1,264 @@ +using Plots, LaTeXStrings, LinearAlgebra, NPZ +using Statistics + + +#--------- Include files +include("plot_functions.jl") +include("./src/compute_obervability_amplitudes_optimized.jl") +include("./src/main_mz_algorithm.jl") +include("./src/mzmd_modes.jl") + + +#--------- Load data +re_num = 9; +method="companion" +# a_method = "lsa"; +a_method = "x0"; +# sample = parse(Int, ARGS[1]); +sample = 1; +println("running e.method $(method) and amp.method $(a_method) eigendecomp computing errors over r, k, d") + +#--------- Load data +t_skip = 3; +dt = t_skip*0.2; +# vort_path = "/Users/l352947/mori_zwanzig/modal_analysis_mzmd/mori_zwanzig/mzmd_code_release/data/vort_all_re600_t5100.npy" +vort_path = "/Users/l352947/mori_zwanzig/mori_zwanzig/mzmd_code_release/data/vort_all_re600_t5100.npy" +X = npzread(vort_path)[:, 1100:5000]; + +# gx = load_grid_file(nx) +m = size(X, 1); +t_end = 1200 #total n steps is 751 +t_pred = 120 +T_train = ((sample):(t_end+sample-1))[1:t_skip:end]; +T_test = ((t_end+sample):(t_end + t_pred + sample))[1:t_skip:end]; +time_train = dt * T_train +time_test = dt * t_skip * T_test; +t_all = 0:dt:(dt*(t_end + t_pred +1)); +X_train = X[:,T_train]; +X_test = X[:,T_test]; + + +t_ = 20 +t_r = round((dt*(t_-2)), digits=2) + +plt = plot_field(X, t_, "DNS ~ ~ t = $(t_r)s") +savefig(plt, "figures/dns_snapshot_$(t_).png") + +#-------- Set params +#number of memory terms: +n_ks = 21 +#number of snapshots used in fitting +t_win = size(T_train, 1) - n_ks - 1; +r = 30; + + +#-------- Compute MZ on POD basis observables +X1 = X_train[:, 1:t_win]; #Used to obtain Ur on 1:twin +X0 = X_train[:, 1]; + +#compute svd based observables (X_proj: states projected onto pod modes) +#method of snapshot more efficient for tall skinny X +S, Ur, X_proj = svd_method_of_snapshots(X_train, r, subtractmean=true) + +#initial condtions with memory for obtaining predictions +X0r = X_proj[:, 1:n_ks]; + +function obtain_mz_operators() + #compute two time covariance matrices + Cov = obtain_C(X_proj, t_win, n_ks); + #compute MZ operators with Mori projection + M, Ω = obtain_ker(Cov, n_ks); + println("----------------------------"); + println("Two time covariance matrices: size(C)) = ", size(Cov)); + println("Markovian operator: size(M) = ", size(M)); + println("All mz operators: size(Ω) = ", size(Ω)); + return Cov, M, Ω +end +Cov, M, Ω = obtain_mz_operators() +plotting_mem_norm(Ω) + +function obtain_omega_norms(Ω) + Ω_norms = zeros(n_ks) + r_ = size(Ω, 2); + I_ = I(r_); + for i in 1:n_ks + if i==1 + Ω_norms[i] = norm(Ω[i,:,:] - I_) + else + Ω_norms[i] = norm(Ω[i,:,:]) + end + end + return Ω_norms +end +Ω_norms = obtain_omega_norms(Ω); +# npzwrite("./omega_norms_r$(r)_k$(n_ksk).npy", Ω_norms) + +plt = plot(Ω_norms, yaxis=:log) +display(plt) + +# C = form_companion(Ω, r, n_ks); + +#---------------------------------------------------------------- + #future state predictions +#---------------------------------------------------------------- + +function simulate_modes(Ur, evects, a, r, λ, dt, t) + a_length = length(a); + t_length = length(t); + time_dynamics = complex(zeros(a_length, t_length)); + ts = (0:t_length-1) * dt; + omega = log.(λ)/dt; + for i in 1 : t_length + time_dynamics[:,i] = (a.*exp.(ts[i] .* omega)); + end + X_pred = evects[1:r,:] * time_dynamics; + Xmz_pred = Ur * X_pred; + return real.(Xmz_pred) +end + +function predict_mzmd(amp_method, t_all, X_proj, C, Ur, r, n_ks) + if amp_method=="x0" + Λc, a_mz, Vc = mzmd_modes_reduced_amps_pred(C, X, Ur, r, n_ks); + solution_modes = simulate_modes(Ur, Vc, a_mz, r, Λc, dt, t_all); + end + if amp_method=="lsa" + Λc, Vc = eigen(C); + e_vects = Vc[1:r,:]; + a, amps, Q = compute_amplitudes_observability_matrix(e_vects, Λc, X_proj, Ur, m); + solution_modes = simulate_modes(Ur, Q, a, r, Λc, dt, t_all); + end + return Λc, Vc, a, amps, Q, solution_modes +end + +function predict_dmd(amp_method, t_all, X0r, M, Ur, r) + if amp_method=="x0" + lambda, e_vects, a, amps = compute_markovian_modes_reduced(M, X0r[:, 1], Ur, r) + solution_modes = simulate_modes(Ur, e_vects, a, r, lambda, dt, t_all); + end + if amp_method=="lsa" + lambda, e_vects = eigen(M); + a, amps, Q = compute_amplitudes_observability_matrix(e_vects, Λ, X_proj, Ur, m); + solution_modes = simulate_modes(Ur, Q, a, r, Λ, dt, t_all); + end + return lambda, e_vects, a, amps, solution_modes +end + +# Λc, Vc, a, amps, Q, solution_modes = predict_mzmd("lsa", 1:751, X_proj, C, Ur, r, n_ks); +# lambda, e_vects, a_m, amps_m, X_dmd = predict_dmd("x0", 1:751, X0r, M, Ur, r); + +# #plot spectrum +# # plt = plot_amplitude_vs_frequency_select_amps(lambda, amps_m, dt, "dmd") +# # display(plt) + +# # plot_c_spectrum(lambda, "dmd") + +# #compute modes and spectrum of MZMD +# Λc_red, Φ_mz, a_mz_red, amp_mz_red, Vc = mzmd_modes_reduced_amps(C, X, Ur, r, n_ks); +# amp_norm = zeros(r*n_ks); +# for i in 1:r*n_ks +# amp_norm[i] = norm(Φ_mz[:,i]*amps[i]); +# end +# plt = plot_amplitude_vs_frequency_select_amps(Λc_red, norm(Q)*sqrt(m)*amp_norm/10, dt, "mzmd") +# display(plt) +# plot_c_spectrum(Λc, "mzmd") + + +# function compute_avg_pressure_sig(X) +# F = reshape(X, (nx, nz, T)); #F[:,:,t] +# F_avgz = mean(F, dims=2); +# p_sig = mean(F_avgz[:,1,:], dims=2); +# return p_sig[:, 1] +# end +# p_amp_dns = compute_avg_pressure_sig(X .+ X_mean); +# p_amp_mz = compute_avg_pressure_sig(solution_modes .+X_mean); +# p_amp_dmd = compute_avg_pressure_sig(X_dmd .+ X_mean); + + +# function plot_p_amps() +# gr(size=(900, 900)) +# plt = plot(p_amp_dns, color="black", label="DNS") +# plot!(p_amp_mz, color="blue", label="MZMD") +# plot!(p_amp_dmd, color="red", label="DMD") +# display(plt) +# end +# plot_p_amps() + +# sorted_indices_mz = sortperm(abs.(amp_norm), rev=true); +# for i in 1:2:15 +# i_mode = sorted_indices_mz[i] +# plt_m = plot_mzmd_modes_axisymmetric(Φ_mz, nx, nz, gx[1:nx], i_mode) +# # savefig(plt_m, "./figures/mzmd_dominant_modes_r$(r)_k$(n_ks)_i$(i_mode).png") +# display(plt_m) +# end + + +function save_dominant_modes() + sorted_indices_mz = sortperm(abs.(amp_norm), rev=true); + for i in 1:2:15 + i_mode = sorted_indices_mz[i] + #plot mode check: + # plt_m = plot_mzmd_modes_axisymmetric(Φ_mz, nx, nz, gx[1:nx], i_mode) + # savefig(plt_m, "./figures/mzmd_dominant_modes_r$(r)_k$(n_ks)_i$(i_mode).png") + # display(plt_m) + Φ = reshape(norm(Q)*sqrt(m)*amp_norm[i_mode]/10 * Φ_mz[:, i_mode]/norm(Φ_mz[:, i_mode]), (nx, nz)); + println("i = ", i_mode, " amp = ", norm(Q)*sqrt(m)*amp_norm[i_mode]/10) + Φ_full = permutedims(hcat(Φ, reverse(Φ, dims=2)), (2, 1)); + npzwrite("./modes/phi_mz_mode_i$(i)_r$(r)_k$(n_ks).npy", Φ_full) + end +end +# save_dominant_modes() + + + + +# #plot diff fields +function obtain_save_diffs() + diff_mz = mean((solution_modes[:, 1:701] .- X[:, 1:701]).^2, dims=2); + diff_dmd = mean((X_dmd[:, 1:701] .- X[:, 1:701]).^2, dims=2); + c_max = maximum(diff_mz/mean(X.^2)) + + dmd_diff = reshape(diff_dmd, (nx, nz)); + dmd_diff = permutedims(hcat(dmd_diff, reverse(dmd_diff, dims=2)), (2, 1)); + + mz_diff = reshape(diff_mz, (nx, nz)); + mz_diff = permutedims(hcat(mz_diff, reverse(mz_diff, dims=2)), (2, 1)); + + npzwrite("./diff_dmd.npy", dmd_diff/mean(X.^2)) + npzwrite("./diff_mz.npy", mz_diff/mean(X.^2)) + + plt = plot_diff_field(diff_dmd/mean(X.^2), nx, nz, gx[1:nx], c_max, "DMD ~ Pointwise ~ MSE") + display(plt) +end + + +function obtain_save_preds() + F2 = reshape(X[:, 751] .+ X_mean, (nx, nz)); + F2 = permutedims(hcat(F2, reverse(F2, dims=2)), (2, 1)); + npzwrite("./dns_tf.npy", F2) + + Fmz = reshape(solution_modes[:, 751] .+ X_mean, (nx, nz)); + Fmz = permutedims(hcat(Fmz, reverse(Fmz, dims=2)), (2, 1)); + npzwrite("./mz_tf.npy", Fmz) + + Fdmd = reshape(X_dmd[:, 751] .+ X_mean, (nx, nz)); + Fdmd = permutedims(hcat(Fdmd, reverse(Fdmd, dims=2)), (2, 1)); + npzwrite("./dmd_tf.npy", Fdmd) +end +# obtain_save_preds() + + +# savefig(plt, "./figures/dmd_diff_field_r$(r).png") + +# plt = plot_diff_field(diff_mz/mean(X.^2), nx, nz, gx[1:nx], c_max, "MZMD ~ Pointwise ~ MSE") +# display(plt) +# savefig(plt, "./figures/mzmd_diff_field_r$(r)_k$(n_ks).png") + +# println("total error mzmd = ", mean(diff_mz)/mean(X.^2)) +# println("total error dmd = ", mean(diff_dmd)/mean(X.^2)) + + +#simulating predictions +# out_path = "./prediction_compare_all_full_r$(r)_k$(n_ks).mov" +# animate_field_comparison(X .+ X_mean, solution_modes .+ X_mean, X_dmd .+ X_mean, nx, nz, gx[1:nx], out_path) + + diff --git a/flow_over_cylinder_re200/spectral_analysis_mzmd.jl b/flow_over_cylinder_re200/spectral_analysis_mzmd.jl new file mode 100644 index 0000000..aaf9dc8 --- /dev/null +++ b/flow_over_cylinder_re200/spectral_analysis_mzmd.jl @@ -0,0 +1,198 @@ +using Plots, LaTeXStrings, LinearAlgebra, NPZ +using Statistics + + +#--------- Include files +include("plot_functions.jl") +include("./src/compute_obervability_amplitudes_optimized.jl") +include("./src/main_mz_algorithm.jl") +include("./src/mzmd_modes.jl") + + +#--------- Load data +# re_num = 9; +re_num = 200 +method="companion" +a_method = "x0"; +sample = 1; +println("running e.method $(method) and amp.method $(a_method) eigendecomp computing errors over r, k, d") + +#--------- Load data +t_skip = 1; +dt = t_skip*0.2; +vort_path = "../data/vort_all_re200_t145.npy" +X = npzread(vort_path) #T=1500 + +X_mean = mean(X, dims=2); +X = X .- X_mean; + +m = size(X, 1); +t_end = 140 #total n steps is 751 +t_pred = 2 +T_train = ((sample):(t_end+sample-1))[1:t_skip:end]; +T_test = ((t_end+sample):(t_end + t_pred + sample))[1:t_skip:end]; +time_train = dt * T_train +time_test = dt * t_skip * T_test; +t_all = 0:dt:(dt*(t_end + t_pred +1)); +X_train = X[:,T_train]; +X_train = X_train .+ 0.25*maximum(X_train) .* randn(size(X_train)) +X_test = X[:,T_test]; + + +t_ = 20 +t_r = round((dt*(t_-2)), digits=2) + +#number of memory terms: +n_ks = 6 +#number of snapshots used in fitting +t_win = size(T_train, 1) - n_ks - 1; +r = 20; + + +#-------- Compute MZ on POD basis observables +X1 = X_train[:, 1:t_win]; #Used to obtain Ur on 1:twin +X0 = X_train[:, 1]; + +#compute svd based observables (X_proj: states projected onto pod modes) +#method of snapshot more efficient for tall skinny X +S, Ur, X_proj = svd_method_of_snapshots(X_train, r, subtractmean=true) +plt = plot_energy_content(S,r) +display(plt) + +#initial condtions with memory for obtaining predictions +X0r = X_proj[:, 1:n_ks]; + +function obtain_mz_operators() + #compute two time covariance matrices + Cov = obtain_C(X_proj, t_win, n_ks); + #compute MZ operators with Mori projection + M, Ω = obtain_ker(Cov, n_ks); + println("----------------------------"); + println("Two time covariance matrices: size(C)) = ", size(Cov)); + println("Markovian operator: size(M) = ", size(M)); + println("All mz operators: size(Ω) = ", size(Ω)); + return Cov, M, Ω +end +Cov, M, Ω = obtain_mz_operators() +plt = plotting_mem_norm(Ω) +display(plt) + +C = form_companion(Ω, r, n_ks); + + +function plot_c_spectrum2(Λc) + gr(size=(570,570), xtickfontsize=12, ytickfontsize=12, xguidefontsize=20, yguidefontsize=14, legendfontsize=12, + dpi=300, grid=(:y, :gray, :solid, 1, 0.4), palette=cgrad(:plasma, 3, categorical = true)); + function circle_shape(h,k,r) + θ=LinRange(0,2*pi,500); + return h .+ r*sin.(θ), k .+ r*cos.(θ) + end + plt = scatter(Λc, linewidth=2.0, ms = 4.0, color = "black", legend=false, grid = true, framestyle = :box) + scatter!([Λc[5]], ms = 6.0, color = "blue", markershape=:diamond) + # scatter!([Λc[6]], ms = 6.0, color = "blue", markershape=:diamond) + scatter!([Λc[1]], ms = 6.0, color = "red", markershape=:star5) + scatter!([Λc[2]], ms = 6.0, color = "red", markershape=:star5) + scatter!([Λc[3]], ms = 6.0, color = "green", markershape=:hexagon) + scatter!([Λc[4]], ms = 6.0, color = "green", markershape=:hexagon) + plot!(circle_shape(0,0,1.0), linewidth=1.5, color="black", linestyle=:dash) + if n_ks>1 + title!(L"\textrm{MZMD ~ Eigenvalues }", titlefont=22) + else + title!(L"\textrm{DMD ~ Eigenvalues }", titlefont=22) + end + xlabel!(L"\textrm{Re}(\lambda)", xtickfontsize=12, xguidefontsize=16) + ylabel!(L"\textrm{Im}(\lambda)", ytickfontsize=12, yguidefontsize=16) + return plt +end + + +function mzmd_modes_reduced_amps_compute(C, Ur, r, n_ks) + #compute modes and amplitudes from C + Λc, Vc = eigen(C); + #mzmd modes: + Φ_mz = Ur*Vc[1:r,:]; + #use initial conditions with k memory terms: + X0_mem = X[:, 1:n_ks]; + #compute amplitudes of modes from observable space + Z0 = Ur' * X0_mem; + z0 = zeros(n_ks*r) + for k in 1:n_ks + # z0[((k-1)*r+1):k*r] = Z0[:,k] + z0[((k-1)*r+1):k*r] = Z0[:,(n_ks - k + 1)] + end + a = pinv(Vc)*z0; + amp = zeros(r*n_ks); + for i in 1:r*n_ks + amp[i] = norm(Φ_mz[:,i]*a[i]); + end + #sort according to largest amplitude: + ind = sortperm(abs.(amp), rev=true) + return Λc[ind], Φ_mz[:, ind], a[ind], amp[ind], Vc +end + +Λc, Φ_mz, a, amp = mzmd_modes_reduced_amps_compute(C, Ur, r, n_ks) #fast and accurate algorithm + +plt = plot_c_spectrum2(Λc) +display(plt) + + +function plot_amplitude_vs_frequency_select_amps_mzmd(lam, a, dt, method) + gr(size=(570,450), xtickfontsize=12, ytickfontsize=12, xguidefontsize=20, yguidefontsize=14, legendfontsize=12, + dpi=300, grid=(:y, :gray, :solid, 1, 0.4), palette=cgrad(:plasma, 3, categorical = true)); + r_2 = floor(Int, r/2); + freqs = imag.(log.(lam))/(2*pi*dt); + # freqs = imag.(log.(lam))/(dt); + # plt = scatter(abs.(real(freqs)), abs.(a)/maximum(abs.(a)), ms=3.0, color="black", legend=false, ylims=(0, 1+0.1)) + if n_ks>1 + plt = plot([abs.(real(freqs))[1], abs.(real(freqs))[1]], [0, (abs.(a)/maximum(abs.(a)))[1]], lw=2, color="blue", xlims=(0, 1.3), ylims=(0, 1+0.15), label=L"\textrm{MZMD}") + else + plt = plot([abs.(real(freqs))[1], abs.(real(freqs))[1]], [0, (abs.(a)/maximum(abs.(a)))[1]], lw=2, color="blue", xlims=(0, 1.3), ylims=(0, 1+0.15), label=L"\textrm{DMD}") + end + scatter!([abs.(real(freqs))[1]], [(abs.(a)/maximum(abs.(a)))[1]], ms = 6.0, color = "red", markershape=:star5, label=false, legend=false) + scatter!([abs.(real(freqs))[5]], [(abs.(a)/maximum(abs.(a)))[5]], ms = 6.0, color = "blue", markershape=:diamond, label=false, legend=false) + scatter!([abs.(real(freqs))[3]], [(abs.(a)/maximum(abs.(a)))[3]], ms = 6.0, color = "green", markershape=:hexagon, label=false, legend=false) + vline!([0.1967], label=L"f_s: \textrm{shedding ~ freq}", lw=2, linestyle=:dash, color=:green, legend=true) + vline!([0.3934], label=L"1^{st} \textrm{higher ~ harmonic}", lw=1.5, linestyle=:dashdot, color=:green, legend=true) + vline!([0.590], label=L"2^{nd} \textrm{higher ~ harmonic}", lw=1.5, linestyle=:dashdotdot, color=:green, legend=true) + # Add vertical lines from each point to the x-axis + for i in 1:length(real(freqs)) + plot!([abs.(real(freqs))[i], abs.(real(freqs))[i]], [0, (abs.(a)/maximum(abs.(a)))[i]], lw=2, color="blue", label="") + end + if n_ks>1 + title!(L"\textrm{Amplitude ~ vs. ~ frequency ~ MZMD ~ }", titlefont=22) + else + title!(L"\textrm{Amplitude ~ vs. ~ frequency ~ DMD ~ }", titlefont=22) + end + xlabel!(L"\textrm{Frequency ~ (Hz)}", xtickfontsize=12, xguidefontsize=16) + ylabel!(L"\textrm{Normalized ~ Amplitude}", ytickfontsize=12, yguidefontsize=16) + display(plt) + # savefig(plt, "./figures/$(method)_spectrum_vs_ampr$(r)_k$(n_ks)_tend$(t_end)_re$(re_num).png") + return plt +end + +plt = plot_amplitude_vs_frequency_select_amps_mzmd(Λc, amp, dt, "mzmd") + + +a_norm = abs.(amp)/maximum(abs.(amp)); +# println("a_norm dominant = ", a_norm[1]) +# println("a_norm first hh = ", a_norm[5]) +# println("a_norm second hh = ", a_norm[3]) + +# function select_save_dominant_modes(idx_mode) +# Φ = reshape(Φ_mz[:, idx_mode]/norm(Φ_mz[:, idx_mode]), (199, 449)); +# println("i = ", idx_mode, " amp = ", a_norm[idx_mode]) +# plt = plot_mzmd_modes(Φ_mz, idx_mode) +# display(plt) +# npzwrite("./modes/phi_mz_mode_i$(idx_mode)_r$(r)_k$(n_ks).npy", Φ) +# end + +# select_save_dominant_modes(1) +# select_save_dominant_modes(5) +# select_save_dominant_modes(3) + + +#check the error (residual) of projecting onto pod modes +# println("checking residual of projecting onto pod modes") +# v_ = randn(size(Ur, 1)); +# error_check = norm((I - Ur * Ur') * v_)/norm(v_) +# println("||(I - Ur * Ur') * v_||/||v_|| = ", error_check) \ No newline at end of file diff --git a/flow_over_cylinder_re200/src/compute_obervability_amplitudes_optimized.jl b/flow_over_cylinder_re200/src/compute_obervability_amplitudes_optimized.jl new file mode 100644 index 0000000..150550b --- /dev/null +++ b/flow_over_cylinder_re200/src/compute_obervability_amplitudes_optimized.jl @@ -0,0 +1,193 @@ +using LinearAlgebra + +function compute_amplitudes_observability_matrix_aq(evects, evals, x_proj) + #inputs: e.vects and e.vals of companion matrix. + #outputs: amplitudes obtained from optimized dmd + Q = evects; + num_modes = size(Q, 2); + #normalize modes + for m=1:num_modes + NormQ=Q[:,m]; + Q[:,m]= Q[:,m]/norm(NormQ,2); + end + M = diagm(evals); #matrix of evals given as vector + r, rk = size(Q); + nobs, T = size(x_proj); + b = zeros(nobs*T); + L = complex.(zeros(r*T, rk)) + aa = I(rk) + for i in 1:T + b[1+(i-1)*nobs:i*nobs] = x_proj[:, i] + L[1+(i-1)*nobs:i*nobs,:] = Q*aa; + aa=aa*M; + end + Ur,Sigmar,Vr = svd(L); + Sigmar = diagm(Sigmar) + a=Vr*(Sigmar\(Ur'*b)); + return a, Q +end + + +function compute_amplitudes_observability_matrix(evects, evals, x_proj, U, num_states) + #inputs: e.vects and e.vals of companion matrix. + #outputs: amplitudes obtained from optimized dmd + Q = evects; + num_modes = size(Q, 2); + #normalize modes + for m=1:num_modes + NormQ=Q[:,m]; + Q[:,m]= Q[:,m]/norm(NormQ,2); + end + M = diagm(evals); #matrix of evals given as vector + r, rk = size(Q); + nobs, T = size(x_proj); + b = zeros(nobs*T); + L = complex.(zeros(r*T, rk)) + aa = I(rk) + for i in 1:T + b[1+(i-1)*nobs:i*nobs] = x_proj[:, i] + L[1+(i-1)*nobs:i*nobs,:] = Q*aa; + aa=aa*M; + end + Ur,Sigmar,Vr = svd(L); + Sigmar = diagm(Sigmar) + a=Vr*(Sigmar\(Ur'*b)); + # a = L\b; + #compute amplitudes for plotting: + # Φ = U*Q; + # amplitude = zeros(rk); + # for i in 1:rk + # amplitude[i] = norm(Φ[:,i]*a[i]); + # end + u=complex.(zeros(r, rk)); + for m=1:num_modes + u[:,m]=a[m]*Q[:,m]; + end + amplitude=zeros(num_modes); + for m=1:num_modes + aca=U*u[:,m]; + amplitude[m]=norm(aca,2)/sqrt(num_states); + end + return a, amplitude, Q +end + +function compute_amplitudes_given_ic_z0_hodmd(Ur, x0_mem, evects, n_td, r) + Z0 = Ur' * x0_mem; + z0 = zeros(n_td*r) + # println(size(Z0)) + # println(size(z0)) + for k in 1:n_td + # z0[((k-1)*r+1):k*r] = Z0[:,(n_td-k+1)] + z0[((k-1)*r+1):k*r] = Z0[:, k] #hodmd order + end + # println(size(pinv(evects))) + a = pinv(evects)*z0; + # a = inv(evects)*z0; + # Q = evects[1:r, :]; + # #TODO why does this cause larger errors? + Q = evects[(n_td-1)*r+1:n_td*r, :]; + return a, Q +end + + +# function compute_amplitudes_given_ic_z0_hodmd(g0, evects, n_td, r) +# #inputs: e.vects and e.vals of companion matrix. +# #outputs: amplitudes obtained from optimized dmd +# a = pinv(evects)*g0; +# # a = inv(evects)*z0; +# # Q = evects[1:r*n_td, :]; +# Q=evects[(n_td-1)*r+1:n_td*r,:]; +# return a, Q +# end + +function compute_amplitudes_given_ic_z0(g0_mem, evects, n_ks, n_td, r) + #inputs: e.vects and e.vals of companion matrix. + #outputs: amplitudes obtained from optimized dmd + #mzmd modes: + z0 = zeros(n_td*n_ks*r); + for k in 1:n_ks + z0[((k-1)*r*n_td+1):(k*r*n_td)] = g0_mem[:,(n_ks-k+1)]; + # z0[((k-1)*r*n_td+1):(k*r*n_td)] = g0_mem[:, k]; + end + a = pinv(evects)*z0; + # a = inv(evects)*z0; + Q = evects[1:r*n_td, :]; + return a, Q +end + +function compute_amplitudes_given_ic_x0(Ur, X0, evects, d, r) + #inputs: e.vects and e.vals of companion matrix. + #outputs: amplitudes obtained from optimized dmd + #mzmd modes: + x0 = Ur' * X0; + a = pinv(evects)*x0; + Q = evects[1:d*r, :]; + return a, Q +end + + + +function compute_amplitudes_given_ic_z0_lsa(z0, evects, evals, U, num_states) + #inputs: e.vects and e.vals of companion matrix. + #outputs: amplitudes obtained from optimized dmd + Q = evects[1:r, :]; #equivalent to P0*Vc + # Q = evects; + num_modes = size(Q, 2); + # Phi_mz = Ur*Q; #must be done before normalizing + #normalize modes + for m=1:num_modes + NormQ=Q[:,m]; + Q[:,m]= Q[:,m]/norm(NormQ,2); + end + M = diagm(evals); #matrix of evals given as vector + r_, rk = size(Q); + Z0 = U'*z0; + nobs, t_ic= size(Z0); + b = zeros(r_*t_ic); + L = complex.(zeros(r_*t_ic, rk)) + aa = I(rk) + println("size L = ", size(L)) + println("size b = ", size(b)) + println("size Z = ", size(Z0)) + for i in 1:t_ic + b[1+(i-1)*nobs:i*nobs] = Z0[:, i] + L[1+(i-1)*nobs:i*nobs,:] = Q*aa; + aa=aa*M; + end + a = L\b; + # Phi_mz = Ur*Q; + # a = Phi_mz\X0; + #compute amplitudes for plotting: + u=complex.(zeros(r_, rk)); + for m=1:num_modes + u[:,m]=a[m]*Q[:,m]; + end + amplitude=zeros(num_modes); + for m=1:num_modes + aca=U*u[:,m]; + amplitude[m]=norm(aca,2)/sqrt(num_states); + end + return a, amplitude, Q +end + +function compute_amplitudes_given_ic_x0(X0, evects, U, num_states) + #inputs: e.vects and e.vals of companion matrix. + #outputs: amplitudes obtained from optimized dmd + Q = evects; + num_modes = size(Q, 2); + Phi_mz = U*Q; #must be done before normalizing + println("size Phi =", size(Phi_mz)) + r_, rk = size(Q); + a = Phi_mz\X0; + #compute amplitudes for plotting: + u=complex.(zeros(r_, rk)); + for m=1:num_modes + u[:,m]=a[m]*Q[:,m]; + end + amplitude=zeros(num_modes); + for m=1:num_modes + aca=U*u[:,m]; + amplitude[m]=norm(aca,2)/sqrt(num_states); + end + return a, amplitude, Q +end diff --git a/flow_over_cylinder_re200/src/dmd_basic_algorithm.jl b/flow_over_cylinder_re200/src/dmd_basic_algorithm.jl new file mode 100644 index 0000000..1c9aa1c --- /dev/null +++ b/flow_over_cylinder_re200/src/dmd_basic_algorithm.jl @@ -0,0 +1,91 @@ +using LinearAlgebra + +#see DMD book by Brunton et al. for more details. +#This is Tu "exact" dmd approach + +function dmd_alg(X1, X2, r, dt) + #inputs: X1 data matrix + #X2 shifted data matrix + #r target rank of SVD + #dt time step + + #outputs: Φ DMD modes + #ω = cont time DMD e.values + #λ, discrete time dmd e.values + #b vector of magnitutes of modes Φ + + U, S, V = svd(X1); + r = minimum([r, size(U)[2]]) + + U_r = U[:, 1:r] + S = Diagonal(S) + S_r = S[1:r, 1:r]; + V_r = V[:, 1:r]; + + Atilde = U_r'*X2*V_r / S_r + lambda, W_r = eigen(Atilde); + Phi = X2 * V_r / S_r * W_r; + + omega = log.(lambda)/dt; + x1 = X1[:, 1]; + b = Phi\x1; + + mm1 = size(X1)[2]; time_dynamics = complex(zeros(r, mm1)); + ts = (0:mm1-1) * dt; + for i in 1 : mm1 + time_dynamics[:,i] = (b.*exp.(ts[i] .* omega)); + end + Xdmd = Phi * time_dynamics; + return Xdmd, Phi, omega, lambda, b +end + +function extract_dmd_modes(X1, X2, r, dt) + U, S, V = svd(X1); + r = minimum([r, size(U)[2]]) + + U_r = U[:, 1:r] + S = Diagonal(S) + S_r = S[1:r, 1:r]; + V_r = V[:, 1:r]; + + Atilde = U_r'*X2*V_r / S_r + lambda, W_r = eigen(Atilde); + Phi = X2 * V_r / S_r * W_r; + omega = log.(lambda)/dt; + + return Phi, omega, lambda +end + +function dmd_prediction(Phi, omega, x1, t, dt, r) + time_dynamics = complex(zeros(r, t)); + ts = (0:t-1) * dt; b = Phi\x1; + for i in 1 : t + time_dynamics[:,i] = (b.*exp.(ts[i] .* omega)); + end + Xdmd = Phi * time_dynamics; + return Xdmd +end + + +function dmd_prediction_given_X(X1, X2, x1, t, dt, r) + U, S, V = svd(X1); + r = minimum([r, size(U)[2]]) + + U_r = U[:, 1:r] + S = Diagonal(S) + S_r = S[1:r, 1:r]; + V_r = V[:, 1:r]; + + Atilde = U_r'*X2*V_r / S_r + lambda, W_r = eigen(Atilde); + Phi = X2 * V_r / S_r * W_r; + omega = log.(lambda)/dt; + + time_dynamics = complex(zeros(r, t)); + ts = (0:t-1) * dt; b = Phi\x1; + for i in 1 : t + time_dynamics[:,i] = (b.*exp.(ts[i] .* omega)); + end + Xdmd = Phi * time_dynamics; + return Xdmd +end diff --git a/flow_over_cylinder_re200/src/homzmd_algorithm.jl b/flow_over_cylinder_re200/src/homzmd_algorithm.jl new file mode 100644 index 0000000..262ad0b --- /dev/null +++ b/flow_over_cylinder_re200/src/homzmd_algorithm.jl @@ -0,0 +1,88 @@ + + +# construct time delay obserables in projected (x_proj) space. Use this to find operators and e.vects +function compute_modes_amplitudes_time_delay_observables(X, T_train, r, n_ks, n_td, T, method, a_method) + #inputs: r number of observables (pod modes) + #n_ks number of mz memory terms + #n_td number of time delay embeddings + S, Ur, X_proj = svd_method_of_snapshots(X[:,T_train], r, subtractmean=true) + # %% Create time delay observable matrix + t_g = T + n_ks + 1; + t_gtilde = length(1:(1+t_g - n_td)); + G_tilde=zeros(n_td*r, t_gtilde); + for i in 1:n_td + G_tilde[(i-1)*r+1:i*r,:] = X_proj[:,i:(i+t_g-n_td)]; + end + #initial conditions for z0 in companion + # g0_test = Ur'*X_test[:, 1:(n_ks+n_td+1)]; + # G0_tilde_test = zeros(n_td*r, n_ks); + # for i in 1:n_td + # G0_tilde_test[(i-1)*r+1:i*r,:] = g0_test[:,i:(i+n_ks-1)]; + # end + + function obtain_mz_operators_(G, t_win, n_ks) + #compute two time covariance matrices + Cov = obtain_C(G, t_win, n_ks); + #compute MZ operators with Mori projection + _, Ω = obtain_ker(Cov, n_ks); + println("----------------------------"); + return Ω + end + twin_tilde = t_gtilde - n_ks - 1; + if n_ks==1 #hodmd (i.e. no mz memory) + # #TODO SVD 2 + U1, Sigma1, T1 = svd(G_tilde); + sigmas1=diagm(Sigma1); + # Second Spatial dimension reduction + r2 = n_td*r + if r2 >= T_train[end] + r2 = size(U1, 2); #this will break companion structure, but is a limit of this method! + end + U1=U1[:,1:r2]; + #svd low rank approximation of G_tilde: + hatG=sigmas1[1:r2,1:r2]*T1[:,1:r2]'; + # #TODO SVD 3 + K1 = size(hatG, 2); + tildeU1, tildeSigma, tildeV1 = svd(hatG[:,1:K1-1]); + tildesigmas = diagm(tildeSigma) + Ω_tilde=hatG[:,2:K1]*tildeV1*inv(tildesigmas)*tildeU1'; + lam, e_vects = eigen(Ω_tilde); + Q=U1*e_vects; + if a_method=="lsa" + a, Q = compute_amplitudes_observability_matrix_aq(Q, lam, G_tilde); + else + a, Q = compute_amplitudes_given_ic_z0_hodmd(Ur, X_test[:, 1:n_td], Q, n_td, r) + end + else + Ω_tilde = obtain_mz_operators_(G_tilde, twin_tilde, n_ks); + if method=="nep" + lam, e_vects = compute_nonlinear_eigen(Ω_tilde); + else + r2 = size(Ω_tilde, 2); + C = form_companion(Ω_tilde, r2, n_ks); + lam, Vc = eigen(C); + end + if a_method=="lsa" + e_vects = Vc[1:r2,:]; + a, Q = compute_amplitudes_observability_matrix_aq(e_vects, lam, G_tilde); + else + a, Q = compute_amplitudes_given_ic_z0(G0_tilde_test, Vc, n_ks, n_td, r) + end + end + return lam, Q, a, Ur, Ω_tilde +end + + +function simulate_modes_time_delay(Ur, evects, a, r, λ, dt, t) + a_length = length(a); + t_length = length(t); + time_dynamics = complex(zeros(a_length, t_length)); + ts = (0:t_length-1) * dt; + omega = log.(λ)/dt; + for i in 1 : t_length + time_dynamics[:,i] = (a.*exp.(ts[i] .* omega)); + end + X_pred = evects * time_dynamics; #in higher dim delay space + Xmz_pred = Ur * X_pred[1:r,:]; + return real.(Xmz_pred) +end \ No newline at end of file diff --git a/flow_over_cylinder_re200/src/main_mz_algorithm.jl b/flow_over_cylinder_re200/src/main_mz_algorithm.jl new file mode 100644 index 0000000..b5898de --- /dev/null +++ b/flow_over_cylinder_re200/src/main_mz_algorithm.jl @@ -0,0 +1,157 @@ +using LinearAlgebra + +function svd_low_rank_proj_observables(X1, X, r) + U, S, V = svd(X1); + r = minimum([r, size(U)[2]]); + + U_r = U[:, 1:r]; + S = Diagonal(S); + S_r = S[1:r, 1:r]; + V_r = V[:, 1:r]; + X_proj = U_r' * X; + return U_r, X_proj +end + +function svd_method_of_snapshots(X, r; subtractmean::Bool = false) + if subtractmean X .-= mean(X,dims=2); end + #number of snapshots + m = size(X, 2); + #Correlation matrix + C = X'*X; + + #Eigen decomp + E = eigen!(C); + e_vals = E.values; + e_vecs = E.vectors; + #sort eigenvectors + sort_ind = sortperm(abs.(e_vals)/m, rev=true) + e_vecs = e_vecs[:,sort_ind]; + e_vals = e_vals[sort_ind]; + #singular values + S = sqrt.(abs.(e_vals)); + #modes and coefficients + U_r = X*e_vecs*Diagonal(1 ./S)[:, 1:r] + # a = Diagonal(S)*e_vecs' + return S, U_r, U_r' * X +end + + + + +function obtain_C(X, t_win, n_ks) + ty = typeof(X[1,1]); m = size(X)[1]; + C = zeros(ty, n_ks+1, m, m); + for δ in 0 : n_ks + C[δ+1,:,:] = X[:, 1+δ : t_win+δ] * X[:, 1 : t_win]'; + end + return C +end + + +function sum_term_C(X, t_win, i, j, k) + ty = typeof(X[1,1]); + S = zeros(ty, 1)[1]; + for l in 1 : (t_win-k) + S += X[i, k+l] * X[j, l] + end + return 1/(t_win - k) * S +end + + +function sum_term_ker(Ω, C, k) + # S = complex(zeros(size(Ω[1,:,:]))); + ty = typeof(Ω[1,1,1]); m, m2 = size(Ω[1,:,:]); + S = zeros(ty, m, m2); + for l in 1 : k - 1 + S += Ω[l,:,:] * C[k - l + 1, :, :]; + end + return S +end + +function obtain_ker(C, n_ks) + ty = typeof(C[1,1,1]); + Ω = zeros(ty, size(C)[1]-1, size(C)[2], size(C)[3]); + C0_inv = pinv(C[1,:,:]); + Ω[1,:,:] = C[2,:,:] * C0_inv; + M = Ω[1,:,:]; + if n_ks > 1 + for k in 2 : n_ks + S = sum_term_ker(Ω, C, k); + Ω[k, :, :] = (C[k+1,:,:] - S) * C0_inv; + end + end + return M, Ω +end + + + +function obtain_Ω_modes(Ω) + ty = typeof(Ω[1,1,1]); + Ω_modes = zeros(ty, size(Ω)[1], size(Ω)[2], size(Ω)[3]); + λ = complex(zeros(size(Ω)[1], size(Ω)[3])); + for n_m in 2 : size(Ω)[1] + λ[n_m, :], Ω_ = eigen(Ω[n_m, :, :]); #Eigen modes of M + Ω_modes[n_m,:,:] = real.(Ω_) + end + return Ω_modes, λ +end + +function obtain_lifted_Ω_modes(Ω_m, Ur) + ty = typeof(Ω_m[1,1,1]); + Ωm_lift = zeros(ty, size(Ω_m)[1], size(Ur)[1], size(Ω_m)[3]); + for n_m in 2 : size(Ω_m)[1] + Ωm_lift[n_m,:,:] = Ur * Ω_m[n_m,:,:]; + end + return Ωm_lift +end + + + +""" + Compute predictions + g((k+1)Δt) = Σ Ω(l) ⋅ g((k-l)Δt) + Wₖ + + For now assume Wₖ is zero +""" + +function sum_term_pred(X, Ω, k, n_ks) + ty = typeof(X[1,1]); + S = zeros(ty, size(X[:,1])[1]); + for l in 1 : minimum([k, n_ks]) + S += Ω[l,:,:] * X[:, k - l + 1]; + end + return S +end + + +function obtain_prediction(X, Ω, n_ks) + ty = typeof(X[1,1]); + X_pred = zeros(ty, size(X)[1], size(X)[2]) + X_pred[:, 1] = X[:, 1]; + for k in 1 : size(X)[2] - 1 + X_pred[:, k+1] = sum_term_pred(X, Ω, k, n_ks) #+ W; + end + return X_pred +end + + +function obtain_future_state_prediction(Ur, X1_gen, Ω, n_ks, T_pred) + ty = typeof(X1_gen[1,1]); + X_pred = zeros(ty, size(X1_gen)[1], T_pred+n_ks) + #initial condition including history: + X_pred[:, 1:(n_ks)] = X1_gen[:, 1:(n_ks)]; + for k in (n_ks) : (T_pred + n_ks - 1) + X_pred[:, k+1] = sum_term_pred(X_pred, Ω, k, n_ks) + end + # return X_pred[:, (n_ks):end] + return X_pred +end + + +function obtain_norm_mem(Ω, n_ks) + out = zeros(n_ks) + for i in 1 : n_ks + out[i] = norm(Ω[i,:,:]) + end + return out +end diff --git a/flow_over_cylinder_re200/src/mzlib.jl b/flow_over_cylinder_re200/src/mzlib.jl new file mode 100644 index 0000000..67b8227 --- /dev/null +++ b/flow_over_cylinder_re200/src/mzlib.jl @@ -0,0 +1,43 @@ +using LinearAlgebra +using Printf +# using ProgressBars + +function compute_mz(g,h) + + M = size(g,1) + N = size(g,2) + C = zeros(M,M,h+2) + + # @printf("number of BLAS threads: %i\n",LinearAlgebra.BLAS.get_num_threads()) + + for k=1:h+2 + # @printf("computing C(%i)...\n",k-1) + temp = zeros(M,M) + # for l in ProgressBar(1:N-k+1) + for l in 1:(N-k+1) + #@printf("\tloop iter #%i of %i)\n",l,N-k+1) + LinearAlgebra.BLAS.ger!(1.0,g[:,l+k-1],g[:,l],temp) + #temp+=g[:,l+k-1]*g[:,l]' + end + C[:,:,k] = temp./(N-k+1) + end + + # print("computing C0 inverse...\n") + Cinv = inv(C[:,:,1]) + + Om = zeros(M,M,h+1) + + # print("computing Om(0)...\n") + Om[:,:,1] = C[:,:,2]*Cinv + + for k=2:h+1 + # @printf("computing Om(%i)...\n",k-1) + temp = zeros(M,M) + for l=1:k-1 + temp += Om[:,:,l]*C[:,:,k-l+1] + end + Om[:,:,k] = (C[:,:,k+1]-temp)*Cinv + end + + return Om,C +end diff --git a/flow_over_cylinder_re200/src/mzmd_modes.jl b/flow_over_cylinder_re200/src/mzmd_modes.jl new file mode 100644 index 0000000..b166656 --- /dev/null +++ b/flow_over_cylinder_re200/src/mzmd_modes.jl @@ -0,0 +1,146 @@ +using LinearAlgebra + +""" +Computing markovian (DMD) modes and modes of MZ with memory: MZMD + -Forms block companion matrix and extracts modes and spectrum. + -Computes predictions with preselected modes +""" + +#Markovian modes: (equivalent to DMD modes) +function compute_markovian_modes(M, Ur, X0, r) + lambda, W_r = eigen(M); + Φ_markovian = Ur * W_r; + a_m = Φ_markovian\X0; #ampliltude of modes + amp = zeros(r) + for i in 1:r + amp[i] = norm(Φ_markovian[:,i]*a_m[i]); + end + ind = sortperm(abs.(amp), rev=true) + return lambda[ind], Φ_markovian[:, ind], a_m[ind] +end + + + +function form_companion(Ω, r, n_ks) + #form block companion matrix + C_r1 = zeros(r, n_ks*r) #row 1 of C + for i in 1:n_ks + C_r1[:, ((i-1)*r+1):i*r] = Ω[i,:,:] + end + if n_ks > 1 + eye = I((n_ks)*r) .+ zeros((n_ks)*r,(n_ks)*r) + eye_sub = eye[1:((n_ks-1)*r), 1:((n_ks)*r)] + end + if n_ks > 1 + C = vcat(C_r1, eye_sub) + else + C = C_r1 + end + return C +end + +function form_companion_memory_only(Ω, r, n_ks) + #form block companion matrix + C_r1 = zeros(r, n_ks*r) #row 1 of C + for i in 2:n_ks + C_r1[:, ((i-1)*r+1):i*r] = Ω[i,:,:] + end + if n_ks > 1 + eye = I((n_ks)*r) .+ zeros((n_ks)*r,(n_ks)*r) + eye_sub = eye[1:((n_ks-1)*r), 1:((n_ks)*r)] + end + if n_ks > 1 + C = vcat(C_r1, eye_sub) + else + C = C_r1 + end + return C +end + +function mzmd_modes(C, X0) + #compute modes and amplitudes from C + Λc, Vc = eigen(C); + #mzmd modes: + Φ_mz = Ur*Vc[1:r,:]; + #compute amplitudes of modes + a = Φ_mz\X0; + amp = zeros(r*n_ks) + for i in 1:r*n_ks + amp[i] = norm(Φ_mz[:,i]*a[i]); + end + #sort according to largest amplitude: + ind = sortperm(abs.(amp), rev=true) + return Λc[ind], Φ_mz[:, ind], a[ind], amp[ind] +end + +function mzmd_modes_reduced_amps(C, Ur, r, n_ks) + #compute modes and amplitudes from C + Λc, Vc = eigen(C); + #mzmd modes: + Φ_mz = Ur*Vc[1:r,:]; + #use initial conditions with k memory terms: + X0_mem = X[:, 1:n_ks]; + #compute amplitudes of modes from observable space + Z0 = Ur' * X0_mem; + z0 = zeros(n_ks*r) + for k in 1:n_ks + # z0[((k-1)*r+1):k*r] = Z0[:,k] + z0[((k-1)*r+1):k*r] = Z0[:,(n_ks - k + 1)] + end + a = pinv(Vc)*z0; + amp = zeros(r*n_ks); + for i in 1:r*n_ks + amp[i] = norm(Φ_mz[:,i]*a[i]); + end + #sort according to largest amplitude: + ind = sortperm(abs.(amp), rev=true) + return Λc[ind], Φ_mz[:, ind], a[ind], amp[ind], Vc +end + + + +function mzmd_modes_full_amps(C) + #NOTE this is more computationally expensive than + #mzmd_modes_reduced_amps and they are equaivlent. + + #X0 needs to include first n_ks snapshots in X + X0_mem = X[:, 1:n_ks]; + #compute modes and amplitudes from C_g + Λc, Vc = eigen(C); + # Λc, Vc = eigen(C, sortby=x->abs.(x)) + #mzmd modes: + Φ_mz = Ur*Vc[1:r,:]; + #set i.c. and compute amplitudes of modes + z0 = zeros(n_ks*m) + for k in 1:n_ks + z0[((k-1)*m+1):k*m] = X0_mem[:,k] + end + #construct: full Φ + ty = typeof(Φ_mz[1,1]); + Φ = zeros(ty, m*n_ks, r*n_ks) + for i in 1:n_ks + # Φ[((i-1)*m+1):i*m, :] = Ur*(Vc[((i-1)*r+1):i*r, :] * diagm(1 ./(Λc.^(i-1)))) + Φ[((i-1)*m+1):i*m, :] = Φ_mz * diagm(1 ./(Λc.^(i-1))) + end + a = (pinv(Φ)*z0) + amp = zeros(r*n_ks) + for i in 1:r*n_ks + amp[i] = norm(Φ_mz[:,i]*a[i]); + end + #sort according to largest amplitude: + ind = sortperm(abs.(amp), rev=true) + return Λc[ind], Φ_mz[:, ind], a[ind], amp[ind] +end + + +function compute_prediction_modes(Φi, ai, λi, tpred) + #eigen-reconstruction with dominant modes Φi + #amplitude: ai, e.val: λi, prediction time steps: tpred + solution = zeros(size(Φi,1), tpred) + for t in 1:tpred + solution[:,t] = real.(Φi * (ai.* (λi.^t))); + end + return solution +end + + diff --git a/flow_over_cylinder_re200/vary_k_dominant_modes.sh b/flow_over_cylinder_re200/vary_k_dominant_modes.sh new file mode 100644 index 0000000..71092e7 --- /dev/null +++ b/flow_over_cylinder_re200/vary_k_dominant_modes.sh @@ -0,0 +1,5 @@ +#!/bin/bash + +for ((i=10; i <= 220; i+=10)); do + julia ./gen_x0_errors_sweep_k_dominant_modes_pred.jl $i +done diff --git a/mzmd_predictions_analysis_cylinder.jl b/mzmd_predictions_analysis_cylinder.jl index 498fa40..0510c93 100644 --- a/mzmd_predictions_analysis_cylinder.jl +++ b/mzmd_predictions_analysis_cylinder.jl @@ -56,6 +56,7 @@ X0 = X_train[:, 1]; #compute svd based observables (X_proj: states projected onto pod modes) #method of snapshot more efficient for tall skinny X S, Ur, X_proj = svd_method_of_snapshots(X_train, r, subtractmean=true) +plot_energy_content(S,r,t_win, 0.99) #initial condtions with memory for obtaining predictions X0r = X_proj[:, 1:n_ks]; diff --git a/prediction_compare_all_full_r12_k5.mov b/prediction_compare_all_full_r12_k5.mov new file mode 100644 index 0000000..383b388 Binary files /dev/null and b/prediction_compare_all_full_r12_k5.mov differ diff --git a/src/mzmd_modes.jl b/src/mzmd_modes.jl index 7357a9b..e52702e 100644 --- a/src/mzmd_modes.jl +++ b/src/mzmd_modes.jl @@ -67,7 +67,8 @@ function mzmd_modes_reduced_amps(X, C, Ur, r, n_ks) Z0 = Ur' * X0_mem; z0 = zeros(n_ks*r) for k in 1:n_ks - z0[((k-1)*r+1):k*r] = Z0[:,k] + # z0[((k-1)*r+1):k*r] = Z0[:,k] + z0[((k-1)*r+1):k*r] = Z0[:,(n_ks - k + 1)]; end a = pinv(Vc)*z0; amp = zeros(r*n_ks);