Variations of the Trapping theorem – different Lyapunov functions
This example will be used to demonstrate that both the energy and enstrophy can be used as Lyapunov functions for promoting local and global stability with the trapping theorems used for quadratically nonlinear systems of ODEs.
We will demonstrate this using a canonical 2D fluid flow (2D so that the enstrophy is conserved). In many cases, the wake behind a bluff body is characterized by a periodic vortex shedding phenomenon known as a von Karman street. The two-dimensional incompressible flow past a cylinder is a stereotypical example of such behavior. The transient energy growth and saturation amplitude of this instability mode is of particular interest and has historically posed a significant modeling challenge. Noack et al. (2003) used an 8-mode POD basis that was augmented with a ninth “shift mode” parameterizing a mean flow deformation. The 9-mode quadratic Galerkin model does resolve the transient dynamics, nonlinear stability mechanism, and post-transient oscillation, accurately reproducing all of the key physical features of the vortex street. In general, the totally symmetric structure in the quadratic coefficients \(Q_{ijk}\) is weakly broken in models for the von Karman street, but it can be enforced to hold exactly.
This is precisely what is done in Schlegel and Noack (2015), and in this perfectly-skew-symmetric case, the global stability of the quadratic model was proven with \(m_9 = m_\text{shift} = \epsilon\), \(\epsilon > 1\), and \(m_i = 0\) for \(i = \{1,...,8\}\). Indeed, we do similarly for data-driven models obtained with the SINDy method in the trapping_sindy_paper_examples.ipynb Jupyter notebook corresponding to the Trapping SINDy paper (https://journals.aps.org/prfluids/abstract/10.1103/PhysRevFluids.6.094401).
### This notebook will show that both the energy and enstrophy can be used with the trapping theorem to promote global stability.
[1]:
import warnings
import matplotlib.gridspec as gridspec
import numpy as np
import scipy.io as sio
from matplotlib import pyplot as plt
from scipy.integrate import solve_ivp
import pysindy as ps
# ignore warnings
warnings.filterwarnings("ignore")
import pymech.neksuite as nek
from trapping_utils import (
galerkin_model,
integrator_keywords,
nel,
nGLL,
interp,
get_velocity,
get_vorticity,
obj_function,
sindy_library_no_bias,
check_local_stability,
make_trap_progress_plots,
make_bar,
nx,
ny,
plot_field,
make_progress_plots,
)
We have pre-loaded some useful functions, and now we can load in the von Karman DNS data. For simplicity (and speed of the code), we will limit ourselves to 5D models, using the first 4 POD modes + the shift mode.
[2]:
# define parameters and load in POD modes obtained from DNS
a = np.loadtxt("../data/vonKarman_pod/vonKarman_a.dat")
t = a[:, 0]
r = 5
a_temp = a[:, 1:r]
a_temp = np.hstack((a_temp, a[:, -1].reshape(3000, 1)))
a = a_temp
# optionally reduce the resolution slightly for more algorithm speed later
tbegin = 0
tend = 3000
skip = 1
t = t[tbegin:tend:skip]
a = a[tbegin:tend:skip, :]
dt = t[1] - t[0]
# define the POD-Galerkin models from Noack (2003)
galerkin9 = sio.loadmat("../data/vonKarman_pod/galerkin9.mat")
# Build two Galerkin models, one in which the nonlinearity is as
# calculated, and the other enforced to be exactly skew-symmetric.
gQ = 0.5 * (galerkin9["Q"] + np.transpose(galerkin9["Q"], [0, 2, 1]))
galerkin9["Q_ep"] = (
gQ
- (
gQ
+ np.transpose(gQ, [1, 0, 2])
+ np.transpose(gQ, [2, 1, 0])
+ np.transpose(gQ, [0, 2, 1])
+ np.transpose(gQ, [2, 0, 1])
+ np.transpose(gQ, [1, 2, 0])
)
/ 6.0
)
model9 = lambda t, a: galerkin_model(a, galerkin9["L"], galerkin9["Q"]) # noqa: E731
model9_ep = lambda t, a: galerkin_model( # noqa: E731
a, galerkin9["L"], galerkin9["Q_ep"]
)
# Generate initial condition from unstable eigenvectors
# lamb, Phi = np.linalg.eig(galerkin9['L'])
# idx = np.argsort(-np.real(lamb))
# lamb, Phi = lamb[idx], Phi[:, idx]
a0 = np.zeros(5)
a0[0] = 1e-3
# get the 5D POD-Galerkin coefficients
inds5 = np.ix_([0, 1, 2, 3, -1], [0, 1, 2, 3, -1])
galerkin5 = {}
galerkin5["L"] = galerkin9["L"][inds5]
inds5 = np.ix_([0, 1, 2, 3, -1], [0, 1, 2, 3, -1], [0, 1, 2, 3, -1])
galerkin5["Q"] = galerkin9["Q"][inds5]
galerkin5["Q_ep"] = galerkin9["Q_ep"][inds5]
model5 = lambda t, a: galerkin_model(a, galerkin5["L"], galerkin5["Q"]) # noqa: E731
# make the 5D POD-Galerkin model trajectories
t_span = (t[0], t[-1])
a_galerkin5 = solve_ivp(model5, t_span, a0, t_eval=t, **integrator_keywords).y.T
adot_galerkin5 = np.gradient(a_galerkin5, axis=0) / (t[1] - t[0])
# plot the first 4 POD modes + the shift mode
mode_numbers = range(10)
plt.figure(figsize=(16, 16))
for i in range(r):
plt.subplot(r, 1, i + 1)
if i == 0:
plt.title(
"DNS and POD-Galerkin models on first 4 POD modes + shift mode", fontsize=16
)
plt.plot(t, a[:, mode_numbers[i]], "r", label="POD from DNS")
plt.plot(t, a_galerkin5[:, mode_numbers[i]], "b", linewidth=2, label="POD-5 model")
ax = plt.gca()
plt.ylabel(r"$a_{" + str(mode_numbers[i]) + "}$", fontsize=20)
plt.grid(True)
if i == r - 1:
plt.xlabel("t", fontsize=18)
plt.legend(loc="upper left", fontsize=16)
else:
ax.set_xticklabels([])
plt.yticks(fontsize=16)
plt.xticks(fontsize=16)
Compute the velocity, vorticity, and enstrophy
The first thing we will do is promote globally stable models via the trapping SINDy algorithm, as in the trapping_sindy_paper_examples.ipynb notebook. We can actually do this by calculating a stability matrix \(\mathbf{A}^S\) with either the energy or the enstrophy (or any other positive definite, quadratic quantity that plays an important role in the dynamics). So we begin by constructing globally stable trapping models of the von Karman street in these two ways.
First, we load in the data and compute the energy and enstrophy.
[3]:
# path to POD mode files
field_path = "../data/vonKarman_pod/cyl0.snapshot"
mode_path = "../data/vonKarman_pod/pod_modes/"
# Read limit cycle flow field for grid points
field = nek.readnek(field_path)
n = nel * nGLL**2
# define cell values needed for the vorticity interpolation
Cx = np.array(
[
field.elem[i].pos[0, 0, j, k]
for i in range(nel)
for j in range(nGLL)
for k in range(nGLL)
]
)
Cy = np.array(
[
field.elem[i].pos[1, 0, j, k]
for i in range(nel)
for j in range(nGLL)
for k in range(nGLL)
]
)
filename = lambda t_idx: "cyl0.f{0:05d}".format(t_idx) # noqa: E731
# plot mean + leading POD modes
clim = [-1, 1]
file_order = [1, 2, 3, 4, 5, 10]
file_labels = ["Mean field", "Mode 1", "Mode 2", "Mode 3", "Mode 4", "Shift mode"]
# Plot the vorticity fields as we load them in
fig = plt.figure(figsize=(10, 8))
spec = gridspec.GridSpec(ncols=2, nrows=3, figure=fig, hspace=0.0, wspace=0.0)
u_list = []
v_list = []
vorticities = []
vorticities_flat = []
for i in range(len(file_order)):
plt.subplot(spec[i])
u, v = get_velocity(mode_path + filename(file_order[i]))
u_list.append(u)
v_list.append(v)
vort = interp(get_vorticity(mode_path + filename(file_order[i])), Cx, Cy)
vorticities.append(np.reshape(vort, [nx, ny], order="F").T)
vorticities_flat.append(get_vorticity(mode_path + filename(file_order[i])))
plot_field(np.reshape(vort, [nx, ny], order="F").T, clim=clim, label=file_labels[i])
plt.xlim([-1, 9])
plt.ylim([-2, 2])
ax = plt.gca()
ax.set_xticklabels([])
ax.set_xticks([])
ax.set_yticks([])
plt.title(file_labels[i], fontsize=16)
plt.show()
We have loaded in the vorticity fields of each POD mode, now need to calculate the energy and the enstrophy.
[4]:
mass_matrix = np.loadtxt("../data/vonKarman_pod/pod_modes/mass_matrix.dat")
ip1 = lambda a, b: np.dot(mass_matrix * a, b) # noqa: E731
ip2 = lambda a, b, c, d: np.dot(a * mass_matrix, c) + np.dot( # noqa: E731
b * mass_matrix, d
)
energy_integrals = np.zeros((6, 6))
enstrophy_integrals = np.zeros((6, 6))
for i, wi in enumerate(vorticities_flat):
for j, wj in enumerate(vorticities_flat):
if i == 0:
enstrophy_integrals[i, j] = ip2(u_list[i], v_list[i], wj, wj)
else:
enstrophy_integrals[i, j] = ip1(wi, wj)
energy_integrals[i, j] = ip2(u_list[i], v_list[i], u_list[j], v_list[j])
Do some checks to make sure energy eigenvalues and enstrophy eigenvalues make sense (energy eigenvalues should be identitity because we are using the eigenbasis of energy, and enstrophy eigenvalues should be positive since enstrophy is by construction a positive definite quantity).
[5]:
# Compute the energy eigenvalues
P_energy = energy_integrals[1:, 1:]
eigs_energy, eigvecs_energy = np.linalg.eigh(P_energy)
print(eigs_energy)
# Compute the enstrophy eigenvalues
P_enstrophy = enstrophy_integrals[1:, 1:]
eigs_enstrophy, eigvecs_enstrophy = np.linalg.eigh(P_enstrophy)
print(eigs_enstrophy)
# Define the linear part of the model,
# rotated into the eigenbasis of enstrophy
L_enstrophy = np.dot(P_enstrophy, galerkin5["L"])
[1. 1. 1. 1. 1.]
[ 2.67263208 2.77593793 2.99681091 10.56822815 10.72500109]
Check global stability of the POD-Galerkin models
Okay, so we have loaded in some DNS data from the von Karman Street and generated (analytic) 5D POD-Galerkin models for this system. The skew-symmetric models below are globally stable if and only if there exists a vector \(\mathbf{m}\) such that following matrix is negative definite:
\[A^S_{ij} = L^S_{ij} + (Q_{ijk} + Q_{jik})m_k.\]
Note that if the quadratic term \(Q_{ijk}\) has no totally-symmetric part this is equal to
\[A^S_{ij} = L^S_{ij} - Q_{kij}m_k.\]
A negative definite \(\mathbf{A}^S\) turns out to also be necessary when \(Q_{ijk}\) does have nontrivial totally symmetric component, but in this case is not sufficient for global boundedness and we can promote local stability as in the trapping_extended.ipynb notebook.
Next we check with a simple nonlinear algorithm (simulated annealing) that our analytic models can be shown to be globally stable (there is an \(\mathbf{m}\) such that \(\mathbf{A}^S\) is negative definite) using both the energy or the enstrophy to construct \(\mathbf{A}^S\).
[6]:
# Import simulated annealing algorithm from scipy
from scipy.optimize import dual_annealing as anneal_algo
# Search between -5000, 5000 for each component of m
boundvals = np.zeros((r, 2))
boundmax = 10
boundmin = -10
boundvals[:, 0] = boundmin
boundvals[:, 1] = boundmax
Ls_enstrophy = 0.5 * (L_enstrophy + L_enstrophy.T)
Ls = 0.5 * (galerkin5["L"] + galerkin5["L"].T)
# Run simulated annealing for the enstrophy eigenbasis
algo_sol = anneal_algo(
obj_function,
bounds=boundvals,
args=(galerkin5["L"], galerkin5["Q_ep"], P_enstrophy), # Factors of P_enstrophy taken care of in obj_func
maxiter=200,
)
opt_m = algo_sol.x
opt_energy = algo_sol.fun
opt_result = algo_sol.message
print("Enstrophy model result:")
print("Optimal m = ", opt_m)
print(
"Algorithm managed to reduce the largest eigenvalue of A^S to eig1 = ",
opt_energy,
"\n",
)
# Repeat using the energy
algo_sol = anneal_algo(
obj_function, bounds=boundvals, args=(galerkin5["L"], galerkin5["Q_ep"], P_energy), maxiter=1000
)
opt_m = algo_sol.x
opt_energy = algo_sol.fun
opt_result = algo_sol.message
print("Energy model result:")
print("Optimal m = ", opt_m)
print(
"Algorithm managed to reduce the largest eigenvalue of A^S to eig1 = ",
opt_energy,
"\n",
)
Enstrophy model result:
Optimal m = [3.52215579e-04 3.02271286e-04 1.16831772e-01 2.96243060e-02
7.79084155e+00]
Algorithm managed to reduce the largest eigenvalue of A^S to eig1 = -0.04389666571454352
Energy model result:
Optimal m = [ 6.29713134e-04 -1.29718141e-04 1.18528961e-01 1.97190841e-02
8.28766982e+00]
Algorithm managed to reduce the largest eigenvalue of A^S to eig1 = -0.04388638908485432
We have proven that both the models, with totally symmetric quadratic tensors, are globally stable.
We now fix the coefficients of \(\mathbf{m}\) except for the direction of the shift-mode, and scan this value to see how the largest eigenvalue of \(\mathbf{A}^S\) (and corresponding trapping region radius) varies.
[7]:
N = 500
alphas = np.linspace(0, 20, N)
r = 5
m = np.zeros(r)
# m[:-1] = [4.40596938e-04, 2.68248063e-04, 1.44130186e-01, 4.11438479e-02]
m_enstrophy = np.zeros(r)
# m_enstrophy[:-1] = [5.99535702e-04, -1.32512857e-04, 1.15691624e-01, 1.80317529e-02]
obj_energy = np.zeros(N)
obj_enstrophy = np.zeros(N)
Rm_energy = np.zeros(N)
Rm_enstrophy = np.zeros(N)
# Extract maximum and minimum eigenvalues,
# and compute radius of the trapping region
max_eig_energy = np.sort(eigs_energy)[-1]
max_eig_enstrophy = np.sort(eigs_enstrophy)[-1]
for i, alpha in enumerate(alphas):
m[-1] = alpha
m_enstrophy[-1] = alpha
obj_energy[i] = obj_function(m, Ls, galerkin5["Q_ep"], P_energy)
obj_enstrophy[i] = obj_function(
m_enstrophy, galerkin5["L"], galerkin5["Q_ep"], P_enstrophy
)
d_energy = np.dot(galerkin5["L"], m) + np.dot(
np.tensordot(galerkin5["Q_ep"], m, axes=([2], [0])), m
)
lsv, sing_vals, rsv = np.linalg.svd(P_enstrophy)
P_rt = lsv @ np.diag(np.sqrt(sing_vals)) @ rsv
d_enstrophy = np.dot(Ls, m_enstrophy) + np.dot(
np.tensordot(galerkin5["Q_ep"], m_enstrophy, axes=([2], [0])), m_enstrophy
)
d_enstrophy = P_rt @ d_enstrophy
Rm_energy[i] = np.linalg.norm(d_energy) / np.abs(obj_energy[i])
Rm_enstrophy[i] = np.linalg.norm(d_enstrophy) / np.abs(obj_enstrophy[i])
plt.figure(figsize=(6, 5))
plt.plot(alphas, obj_energy * 1e2, label=r"$\lambda_1^{energy}$ x $10^2$")
plt.plot(
alphas[obj_energy < 0],
Rm_energy[obj_energy < 0],
label="$R_m^{energy}$",
)
plt.plot(alphas, obj_enstrophy * 1e2, label=r"$\lambda_1^{enstrophy}$ x $10^2$")
plt.plot(
alphas[obj_enstrophy < 0],
Rm_enstrophy[obj_enstrophy < 0],
label="$R_m^{enstrophy}$",
)
plt.legend(fontsize=12, loc="upper right", framealpha=1.0)
plt.ylim(-5, 50)
# plt.xlim(0, 10000)
plt.grid(True)
plt.xlabel(r"$m_{shift}$", fontsize=18)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
Data-driven 5D models for the von Karman street using Trapping SINDy
We have investigated a number of 5D (analytic) POD-Galerkin models and shown that the analytic, constrained models are globally stable. Now we show that trapping SINDy can be used with the energy or enstrophy to build data-driven models that are also globally stable.
[8]:
galerkin_ep = (
gQ
+ np.transpose(gQ, [1, 0, 2])
+ np.transpose(gQ, [2, 1, 0])
+ np.transpose(gQ, [0, 2, 1])
+ np.transpose(gQ, [2, 0, 1])
+ np.transpose(gQ, [1, 2, 0])
) / 6.0
print(np.max(abs(galerkin_ep)))
gQ_enstrophy = enstrophy_integrals[1:, 1:] @ gQ[inds5]
galerkin_ep = (
gQ_enstrophy
+ np.transpose(gQ_enstrophy, [1, 0, 2])
+ np.transpose(gQ_enstrophy, [2, 1, 0])
+ np.transpose(gQ_enstrophy, [0, 2, 1])
+ np.transpose(gQ_enstrophy, [2, 0, 1])
+ np.transpose(gQ_enstrophy, [1, 2, 0])
) / 6.0
print(np.max(abs(galerkin_ep)))
0.0003931843531373586
0.04179936940675318
[ ]:
# define hyperparameters
max_iter = 10000
eta = 1.0
# don't need a reg_weight_lam if eta is sufficiently small
# which is good news because CVXPY is much slower
reg_weight_lam = 0
alpha_m = 1e-1 * eta
# run trapping SINDy
sindy_opt = ps.TrappingSR3(
_n_tgts=5,
_include_bias=False,
reg_weight_lam=reg_weight_lam,
eta=eta,
alpha_m=alpha_m,
max_iter=max_iter,
verbose=True,
)
model = ps.SINDy(
optimizer=sindy_opt,
feature_library=sindy_library_no_bias,
differentiation_method=ps.FiniteDifference(drop_endpoints=True),
)
model.fit(a, t=t)
Xi = model.coefficients().T
PL_tensor = sindy_opt.PL_
PQ_tensor = sindy_opt.PQ_
L = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))
Q = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))
Q_sum = np.max(np.abs((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1]))))
print("Max deviation from the constraints = ", Q_sum)
The previous model finds an \(\mathbf{m}\) such that \(\mathbf{A}^S\) is negative definite, while also fitting the data. Now we can repeat but in the eigenbasis of enstrophy. If the enstrophy is \(H = \mathbf{y}^T\mathcal{P}\mathbf{A}^S \mathbf{y}\), now we are searching for an \(\mathbf{m}\) such that \(\mathcal{P}\mathbf{A}^S\) is negative definite.
[ ]:
max_iter = 2000
eta = 1.0
reg_weight_lam = 0
alpha_m = 1e-2 * eta
mod_matrix = enstrophy_integrals[1:, 1:]
sindy_opt = ps.TrappingSR3(
method="global",
_n_tgts=5,
_include_bias=False,
reg_weight_lam=reg_weight_lam,
eta=eta,
alpha_m=alpha_m,
max_iter=max_iter,
mod_matrix=mod_matrix,
verbose=True,
)
model = ps.SINDy(
optimizer=sindy_opt,
feature_library=sindy_library_no_bias,
)
model.fit(a, t=t)
Xi = model.coefficients().T
Lenstrophy = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))
Qenstrophy = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))
mean_val = np.mean(a, axis=0)
mean_val = np.sqrt(np.sum(mean_val**2))
check_local_stability(Xi, sindy_opt, mean_val)
enstrophy_model = model
# Q = np.tensordot(sindy_opt.PQ_, Xi, axes=([4, 3], [0, 1]))
# Q_sum = np.max(np.abs((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1]))))
# print("Maximum deviation from the constraints = ", Q_sum)
make_progress_plots(r, sindy_opt)
Enstrophy model was successful! #### We have built two data-driven models, one using the energy as a Lyapunov function for trapping SINDy, and the other using the enstrophy. Now we compare the distribution of coefficients in each model.
[ ]:
make_bar(galerkin5, L, Q, Lenstrophy, Qenstrophy)
Compare the models
Below, we compare the 5D POD-Galerkin with trapping SINDy models obtained by considering the energy or by the enstrophy. Both trapping models outperform the POD-Galerkin model.
[ ]:
# Interpolate onto better time base
t_traj = np.linspace(t[0], t[-1], len(t) * 1)
# simulate trapping SINDy results
xtraj_energy = model.simulate(a0, t_traj)
xtraj_enstrophy = enstrophy_model.simulate(a0, t_traj)
# simulate and plot 5D von Karman trajectory results
t_span = (t_traj[0], t_traj[-1])
xtraj_pod9 = solve_ivp(model5, t_span, a0, t_eval=t_traj, **integrator_keywords).y.T
# Make awesome plot
fig, ax = plt.subplots(3, 3, subplot_kw={"projection": "3d"}, figsize=(16, 16))
data = [
a[:, [0, 1, -1]],
xtraj_energy[:, [0, 1, -1]],
a[:, [0, 1, -1]],
xtraj_pod9[:, [0, 1, -1]],
a[:, [0, 1, -1]],
xtraj_enstrophy[:, [0, 1, -1]],
a[:, [0, 1, 2]],
xtraj_energy[:, [0, 1, 2]],
a[:, [0, 1, 2]],
xtraj_pod9[:, [0, 1, 2]],
a[:, [0, 1, 2]],
xtraj_enstrophy[:, [0, 1, 2]],
a[:, [2, 3, -1]],
xtraj_energy[:, [2, 3, -1]],
a[:, [2, 3, -1]],
xtraj_pod9[:, [2, 3, -1]],
a[:, [2, 3, -1]],
xtraj_enstrophy[:, [2, 3, -1]],
]
data_labels = [
[r"$a_1$", r"$a_2$", r"$a_{shift}$"],
[r"$a_1$", r"$a_2$", r"$a_3$"],
[r"$a_2$", r"$a_3$", r"$a_{shift}$"],
]
for i in range(3):
ax[i, 0].plot(
data[6 * i][:, 0],
data[6 * i][:, 1],
data[6 * i][:, 2],
color="r",
label="POD trajectory from DNS",
)
ax[i, 0].plot(
data[6 * i + 1][:, 0],
data[6 * i + 1][:, 1],
data[6 * i + 1][:, 2],
color="k",
label="Trapping SINDy model (energy)",
)
ax[i, 1].plot(
data[6 * i + 2][:, 0],
data[6 * i + 2][:, 1],
data[6 * i + 2][:, 2],
color="r",
label="POD trajectory from DNS",
)
ax[i, 1].plot(
data[6 * i + 3][:, 0],
data[6 * i + 3][:, 1],
data[6 * i + 3][:, 2],
color="b",
label="(analytic) POD-Galerkin model",
)
ax[i, 2].plot(
data[6 * i + 4][:, 0],
data[6 * i + 4][:, 1],
data[6 * i + 4][:, 2],
color="r",
label="POD trajectory from DNS",
)
ax[i, 2].plot(
data[6 * i + 5][:, 0],
data[6 * i + 5][:, 1],
data[6 * i + 5][:, 2],
color="g",
label="Trapping SINDy model (enstrophy)",
)
ax[i, 0].legend(fontsize=10)
ax[i, 0].set_xticklabels([])
ax[i, 0].set_yticklabels([])
ax[i, 0].set_zticklabels([])
ax[i, 0].set_xlabel(data_labels[i][0], fontsize=14)
ax[i, 0].set_ylabel(data_labels[i][1], fontsize=14)
ax[i, 0].set_zlabel(data_labels[i][2], fontsize=14)
ax[i, 1].legend(fontsize=10)
ax[i, 1].set_xticklabels([])
ax[i, 1].set_yticklabels([])
ax[i, 1].set_zticklabels([])
ax[i, 1].set_xlabel(data_labels[i][0], fontsize=14)
ax[i, 1].set_ylabel(data_labels[i][1], fontsize=14)
ax[i, 1].set_zlabel(data_labels[i][2], fontsize=14)
ax[i, 2].legend(fontsize=10)
ax[i, 2].set_xticklabels([])
ax[i, 2].set_yticklabels([])
ax[i, 2].set_zticklabels([])
ax[i, 2].set_xlabel(data_labels[i][0], fontsize=14)
ax[i, 2].set_ylabel(data_labels[i][1], fontsize=14)
ax[i, 2].set_zlabel(data_labels[i][2], fontsize=14)
plt.show()
Rerun with enstrophy using local trapping method
[12]:
max_iter = 5000
eta = 1.0e4
reg_weight_lam = 0
alpha_m = 1e-1 * eta
# alpha = 1e10
beta = 1e-10
mod_matrix = enstrophy_integrals[1:, 1:]
sindy_opt = ps.TrappingSR3(
method="local",
_n_tgts=5,
_include_bias=False,
reg_weight_lam=reg_weight_lam,
eta=eta,
gamma=-1,
alpha_m=alpha_m,
# alpha=alpha,
beta=beta,
max_iter=max_iter,
verbose=True,
mod_matrix=mod_matrix,
)
model = ps.SINDy(
optimizer=sindy_opt,
feature_library=sindy_library_no_bias,
)
model.fit(a, t=t)
Xi = model.coefficients().T
PL_tensor = sindy_opt.PL_
PQ_tensor = sindy_opt.PQ_
Lenstrophy = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))
Qenstrophy = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))
mean_val = np.mean(a, axis=0)
mean_val = np.sqrt(np.sum(mean_val**2))
check_local_stability(Xi, sindy_opt, mean_val)
enstrophy_model = model
Q = np.tensordot(sindy_opt.PQ_, Xi, axes=([4, 3], [0, 1]))
Q = np.tensordot(mod_matrix, Q, axes=([1], [0]))
Q_sum = np.max(np.abs((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1]))))
print("Maximum deviation from the constraints = ", Q_sum)
rhos_minus, rhos_plus = make_trap_progress_plots(r, sindy_opt)
Iter ... |y-Xw|^2 ... |Pw-A|^2/eta ... |w|_1 ... |Qijk|/a ... |Qijk+...|/b ... Total:
0 ... 4.375e+00 ... 4.454e-04 ... 0.00e+00 ... 3.82e-19 ... 6.97e-09 ... 4.38e+00
500 ... 4.375e+00 ... 1.998e-04 ... 0.00e+00 ... 3.82e-19 ... 6.97e-09 ... 4.38e+00
1000 ... 4.375e+00 ... 1.387e-04 ... 0.00e+00 ... 3.82e-19 ... 6.98e-09 ... 4.38e+00
1500 ... 4.375e+00 ... 1.108e-04 ... 0.00e+00 ... 3.82e-19 ... 6.98e-09 ... 4.38e+00
2000 ... 4.375e+00 ... 9.291e-05 ... 0.00e+00 ... 3.82e-19 ... 6.98e-09 ... 4.38e+00
2500 ... 4.375e+00 ... 8.194e-05 ... 0.00e+00 ... 3.82e-19 ... 6.98e-09 ... 4.38e+00
3000 ... 4.375e+00 ... 7.462e-05 ... 0.00e+00 ... 3.82e-19 ... 6.98e-09 ... 4.38e+00
3500 ... 4.375e+00 ... 6.917e-05 ... 0.00e+00 ... 3.82e-19 ... 6.98e-09 ... 4.38e+00
4000 ... 4.375e+00 ... 6.510e-05 ... 0.00e+00 ... 3.82e-19 ... 6.98e-09 ... 4.37e+00
4500 ... 4.375e+00 ... 6.207e-05 ... 0.00e+00 ... 3.82e-19 ... 6.98e-09 ... 4.37e+00
optimal m: [ 2.12560595 -1.82937938 0.76799457 2.45305516 47.61642396]
As eigvals: [-3.15773393 -1.5835335 -0.92314318 -0.36068861 -0.1157516 ]
0.5 * |tilde{H}_0|_F = 5.906654037631929e-10
0.5 * |tilde{H}_0|_F^2 / beta = 6.977712384054714e-09
Estimate of trapping region size, Rm = 4005.17475735993
Normalized trapping region size, Reff = 2382.67546579448
Local stability size, R_ls= 293948204.118023
Maximum deviation from the constraints = 3.811306881207344e-09
[13]:
make_progress_plots(r, sindy_opt)
[ ]:
fs = 18
t_traj = np.linspace(t[0], t[-1], len(t)) # * 2)
fig, ax = plt.subplots(3, 3, subplot_kw={"projection": "3d"}, figsize=(16, 16))
a0s = [(np.random.rand(5) - 0.5) * 20 for i in range(2)]
for a0 in a0s:
print(a0)
# simulate trapping SINDy results
xtraj_energy = model.simulate(a0, t_traj)
xtraj_enstrophy = enstrophy_model.simulate(a0, t_traj)
# simulate and plot 5D von Karman trajectory results
t_span = (t_traj[0], t_traj[-1])
xtraj_pod9 = solve_ivp(model5, t_span, a0, t_eval=t_traj, **integrator_keywords).y.T
# Make awesome plot
data = [
a[:, [0, 1, -1]],
xtraj_energy[:, [0, 1, -1]],
a[:, [0, 1, -1]],
xtraj_pod9[:, [0, 1, -1]],
a[:, [0, 1, -1]],
xtraj_enstrophy[:, [0, 1, -1]],
a[:, [0, 1, 2]],
xtraj_energy[:, [0, 1, 2]],
a[:, [0, 1, 2]],
xtraj_pod9[:, [0, 1, 2]],
a[:, [0, 1, 2]],
xtraj_enstrophy[:, [0, 1, 2]],
a[:, [2, 3, -1]],
xtraj_energy[:, [2, 3, -1]],
a[:, [2, 3, -1]],
xtraj_pod9[:, [2, 3, -1]],
a[:, [2, 3, -1]],
xtraj_enstrophy[:, [2, 3, -1]],
]
data_labels = [
[r"$a_1$", r"$a_2$", r"$a_{shift}$"],
[r"$a_1$", r"$a_2$", r"$a_3$"],
[r"$a_2$", r"$a_3$", r"$a_{shift}$"],
]
ax[0, 0].set_title("Global Trapping SINDy, energy", fontsize=22)
ax[0, 1].set_title("POD-Galerkin model", fontsize=22)
ax[0, 2].set_title("Local Trapping SINDy, enstrophy", fontsize=22)
for i in range(3):
ax[i, 0].plot(
data[6 * i][:, 0],
data[6 * i][:, 1],
data[6 * i][:, 2],
"k--",
label="POD trajectory from DNS",
)
ax[i, 0].plot(
data[6 * i + 1][:, 0],
data[6 * i + 1][:, 1],
data[6 * i + 1][:, 2],
linestyle="--",
label="Trapping SINDy model (energy)",
)
ax[i, 1].plot(
data[6 * i + 2][:, 0],
data[6 * i + 2][:, 1],
data[6 * i + 2][:, 2],
"k--",
label="POD trajectory from DNS",
)
ax[i, 1].plot(
data[6 * i + 3][:, 0],
data[6 * i + 3][:, 1],
data[6 * i + 3][:, 2],
linestyle="--",
label="(analytic) POD-Galerkin model",
)
ax[i, 2].plot(
data[6 * i + 4][:, 0],
data[6 * i + 4][:, 1],
data[6 * i + 4][:, 2],
"k--",
label="POD trajectory from DNS",
)
ax[i, 2].plot(
data[6 * i + 5][:, 0],
data[6 * i + 5][:, 1],
data[6 * i + 5][:, 2],
linestyle="--",
label="Trapping SINDy model (enstrophy)",
)
if i == 2:
xlim = 1
ylim = 1
else:
xlim = 2
ylim = 2
if i == 1:
zlim = 0.5
zlim_minus = -0.5
else:
zlim = 2.5
zlim_minus = 0
# ax[i, 0].legend(fontsize=10)
ax[i, 0].set_xticklabels([])
ax[i, 0].set_yticklabels([])
ax[i, 0].set_zticklabels([])
ax[i, 0].set_xlabel(data_labels[i][0], fontsize=fs)
ax[i, 0].set_ylabel(data_labels[i][1], fontsize=fs)
ax[i, 0].set_zlabel(data_labels[i][2], fontsize=fs)
ax[i, 0].set_xlim(-xlim, xlim)
ax[i, 0].set_ylim(-ylim, ylim)
ax[i, 0].set_zlim(zlim_minus, zlim)
# ax[i, 1].legend(fontsize=10)
ax[i, 1].set_xticklabels([])
ax[i, 1].set_yticklabels([])
ax[i, 1].set_zticklabels([])
ax[i, 1].set_xlabel(data_labels[i][0], fontsize=fs)
ax[i, 1].set_ylabel(data_labels[i][1], fontsize=fs)
ax[i, 1].set_zlabel(data_labels[i][2], fontsize=fs)
ax[i, 1].set_xlim(-xlim, xlim)
ax[i, 1].set_ylim(-ylim, ylim)
ax[i, 1].set_zlim(zlim_minus, zlim)
# ax[i, 2].legend(fontsize=10)
ax[i, 2].set_xticklabels([])
ax[i, 2].set_yticklabels([])
ax[i, 2].set_zticklabels([])
ax[i, 2].set_xlabel(data_labels[i][0], fontsize=fs)
ax[i, 2].set_ylabel(data_labels[i][1], fontsize=fs)
ax[i, 2].set_zlabel(data_labels[i][2], fontsize=fs)
ax[i, 2].set_xlim(-xlim, xlim)
ax[i, 2].set_ylim(-ylim, ylim)
ax[i, 2].set_zlim(zlim_minus, zlim)
# plt.savefig("enstrophy_results.pdf")
[ 3.4093502 -1.65390395 1.17379657 -7.19226123 -6.03797022]
[ 6.01489137 9.36523151 -3.73151644 3.84645231 7.52778305]
capi_return is NULL
Call-back cb_f_in_lsoda__user__routines failed.
Traceback (most recent call last):
File "/Users/akaptanoglu/anaconda3/envs/pysindy_trapping/lib/python3.11/site-packages/scipy/integrate/_ivp/base.py", line 154, in fun
Fatal Python error: F2PySwapThreadLocalCallbackPtr: F2PySwapThreadLocalCallbackPtr: PyLong_AsVoidPtr failed
Python runtime state: initialized
return self.fun_single(t, y)
^^^^^^^^^^^^^^^^^^^^^
File "/Users/akaptanoglu/anaconda3/envs/pysindy_trapping/lib/python3.11/site-packages/scipy/integrate/_ivp/base.py", line 23, in fun_wrapped
return np.asarray(fun(t, y), dtype=dtype)
^^^^^^^^^
File "/Users/akaptanoglu/pysindy/pysindy/pysindy.py", line 684, in rhs
return self.predict(x[np.newaxis, :])[0]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/akaptanoglu/pysindy/pysindy/pysindy.py", line 316, in predict
result = [self.model.predict([xi]) for xi in x]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/akaptanoglu/pysindy/pysindy/pysindy.py", line 316, in <listcomp>
result = [self.model.predict([xi]) for xi in x]
^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/akaptanoglu/anaconda3/envs/pysindy_trapping/lib/python3.11/site-packages/sklearn/pipeline.py", line 602, in predict
Xt = transform.transform(Xt)
^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/akaptanoglu/anaconda3/envs/pysindy_trapping/lib/python3.11/site-packages/sklearn/utils/_set_output.py", line 295, in wrapped
data_to_wrap = f(self, X, *args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/akaptanoglu/pysindy/pysindy/feature_library/base.py", line 148, in func
result = wrapped_func(self, xs, *args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/akaptanoglu/pysindy/pysindy/feature_library/polynomial_library.py", line 253, in transform
xp[..., i] = x[..., combo].prod(-1)
~^^^^^^^^^^^^
File "/Users/akaptanoglu/pysindy/pysindy/utils/axes.py", line 316, in __getitem__
key, adv_inds = _standardize_indexer(self, key)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/akaptanoglu/pysindy/pysindy/utils/axes.py", line 700, in _standardize_indexer
ax_key = np.array(ax_key)
^^^^^^^^^^^^^^^^
KeyboardInterrupt
Extension modules: zmq.backend.cython._zmq
[ ]:
# define energies of the DNS, and both the 5D and 9D models
# for POD-Galerkin and the trapping SINDy models
E = np.sum(a**2, axis=1)
E_galerkin5 = np.sum(xtraj_pod9**2, axis=1)
E_sindy5 = np.sum(xtraj_energy**2, axis=1)
E_sindy5_local = np.sum(xtraj_enstrophy**2, axis=1)
# plot the energies
plt.figure(figsize=(16, 4))
plt.plot(t, E, "r", label="DNS")
plt.plot(t_traj, E_galerkin5, "m", label="POD-5")
plt.plot(t_traj, E_sindy5, "k", label=r"SINDy-energy")
plt.plot(t_traj, E_sindy5_local, "b", label=r"SINDy-enstrophy")
# do some formatting and save
plt.legend(fontsize=22, loc=2)
plt.grid()
plt.xlim([0, 300])
plt.ylim(0, 20)
ax = plt.gca()
ax.set_yticks([0, 10, 20])
ax.tick_params(axis="x", labelsize=20)
ax.tick_params(axis="y", labelsize=20)
plt.ylabel("Total energy", fontsize=20)
plt.xlabel("t", fontsize=20)
plt.show()