CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download

Bálint Kaszás and George Haller, 'Capturing the Edge of Chaos as a Spectral Submanifold in Pipe Flows ', Journal of Fluid Mechanics, submitted, Figure 5 key words: spectral submanifold, reduced-order model, edge of chaos, pipe flow

Views: 6
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2204
Tweet
Share
Share
Kernel: Python 3 (system-wide)

logo

Visualization of the mixed-mode SSM in 3D

In this notebook we use the training data located in the folder data/ generated by Openpipeflow (developed by Ashley Willis, https://openpipeflow.org/) and the previously obtained SSM-reduced model to visualize the geometry of the phase space of the Navier-Stokes equations.

  • We visualize the two-dimensional SSM-reduced model as a vector field and explicitly calculate the stable manifold of the lower branch traveling wave. Since the edge of chaos intersects the mixed-mode SSM transversely, we can conclude that this intersection is formed as the stable manifold of the lower branch traveling wave in the SSM-reduced model. Using the parametrization map we then lift the intersection of the mixed-mode SSM and the edge of chaos to the full phase space.

We recall that the SSM-reduced model has already been selected in Fig. 4. The polynomial order of the parametrization

u(t)=W(J(t),K(t))\mathbf{u}(t) = \mathbf{W}(J(t), K(t))

and the SSM-reduced dynamics

ddt(JK)=f(J,K)=∑n+m≤Mr(Rnm(J)KmJnRnm(K)KmJn),\frac{d}{dt}\begin{pmatrix}J \\ K \end{pmatrix} = \mathbf{f}(J,K)= \sum_{n+m\leq M_r} \begin{pmatrix} R^{(J)}_{nm} K^m J^n \\ R^{(K)}_{nm} K^m J^n\end{pmatrix},

was Mr=5M_r=5.

from joblib import load
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

We load the previously calculated model

ssm = load('saved_models/ssm_model')
INFO 2023-05-30 22:10:21 utils NumExpr defaulting to 4 threads. /usr/local/lib/python3.10/dist-packages/sklearn/base.py:318: UserWarning: Trying to unpickle estimator PolynomialFeatures from version 1.2.1 when using version 1.2.2. This might lead to breaking code or invalid results. Use at your own risk. For more info please refer to: https://scikit-learn.org/stable/model_persistence.html#security-maintainability-limitations warnings.warn( /usr/local/lib/python3.10/dist-packages/sklearn/base.py:318: UserWarning: Trying to unpickle estimator Ridge from version 1.2.1 when using version 1.2.2. This might lead to breaking code or invalid results. Use at your own risk. For more info please refer to: https://scikit-learn.org/stable/model_persistence.html#security-maintainability-limitations warnings.warn( /usr/local/lib/python3.10/dist-packages/sklearn/base.py:318: UserWarning: Trying to unpickle estimator Pipeline from version 1.2.1 when using version 1.2.2. This might lead to breaking code or invalid results. Use at your own risk. For more info please refer to: https://scikit-learn.org/stable/model_persistence.html#security-maintainability-limitations warnings.warn( /usr/local/lib/python3.10/dist-packages/sklearn/base.py:318: UserWarning: Trying to unpickle estimator StandardScaler from version 1.2.1 when using version 1.2.2. This might lead to breaking code or invalid results. Use at your own risk. For more info please refer to: https://scikit-learn.org/stable/model_persistence.html#security-maintainability-limitations warnings.warn( 
import shapely


from shapely.geometry import Polygon, Point

linecolor = [0.5111196874841155, 0.573802174897223, 0.1737010491405684]
surfacecolor = [0.8798537965385758, 0.8910197814921311, 0.5854008624716268]

We briefly repeat here the calculation to obtain the stable manifold of the lower branch traveling wave

def computeSpectrum(fixedPoint, delta, dyn):
#    fixedPoint = np.array(res.x).reshape(1,-1)
    xpdx = fixedPoint + np.array([delta, 0])
    xmdx = fixedPoint - np.array([delta, 0])
    xpdy = fixedPoint + np.array([0, delta])
    xmdy = fixedPoint - np.array([0, delta])

    fxpdx = dyn(xpdx)[0]
    fxmdx = dyn(xmdx)[0]

    fxpdy = dyn(xpdy)[0]
    fxmdy = dyn(xmdy)[0]


    dfxdx = (fxpdx[0] - fxmdx[0])/(2*delta)
    dfxdx = (fxpdx[0] - fxmdx[0])/(2*delta)

    dfxdy = (fxpdy[0] - fxmdy[0])/(2*delta)
    dfxdy = (fxpdy[0] - fxmdy[0])/(2*delta)

    dfydx = (fxpdx[1] - fxmdx[1])/(2*delta)
    dfydx = (fxpdx[1] - fxmdx[1])/(2*delta)


    dfydy = (fxpdy[1] - fxmdy[1])/(2*delta)
    dfydy = (fxpdy[1] - fxmdy[1])/(2*delta)
    A = np.array([[dfxdx, dfxdy], [dfydx, dfydy]])
    w, v = np.linalg.eig(A)
    return A, w, v

lb_jk_coords = [0.101187, 0.1011201]

predictor = lambda x: ssm.reduced_dynamics.predict(x)
def tocallfunDynamics(t, x, dyn):
    
    dy = dyn(np.array(x).reshape(1,-1))
    return np.squeeze(np.array(dy))
def fixedPt(x, dyn):
    xx = np.array(x)
    odefun = lambda t, x : tocallfunDynamics(t, x, dyn)
    return odefun(0, xx)

fixedPoint = np.array(lb_jk_coords).reshape(1,-1)

Df, w, v = computeSpectrum(fixedPoint, 1e-6, predictor)
ic = lb_jk_coords + 1e-6*v[:,0]

ic2 = lb_jk_coords - 1e-6*v[:,0]
teval = np.linspace(0, -1000, 5000)


solStab1 = solve_ivp(tocallfunDynamics, [0, -1000], ic, method = 'DOP853', t_eval = teval, rtol = 1e-12, atol = 1e-12, args=(predictor, ))
solStab2 = solve_ivp(tocallfunDynamics, [0, -1000], ic2, method = 'DOP853', t_eval = teval, rtol = 1e-12, atol = 1e-12, args=(predictor, ))

Selecting the domain

We select a bounded region in the J-K plane, that is not too far away from the training data. This will be then lifted to the full phase space using the parametrization map.

This can be changed manually. 

points = [[0.1294, 0.1198], [0.1109, 0.0860], [0.0856, 0.0492], [0.0559, 0.0163], [0.0099, -0.0042],
          [0,0], [0.0058, 0.0064], [0.0339, 0.0284], [0.0464, 0.0553], [0.0537, 0.0863], [0.0692, 0.1253], [0.0929, 0.1348]
    
]

polygon = Polygon(points)
polygon_smooth = polygon.buffer(0.02, join_style=1)#.buffer(-0.1, join_style=1)


x_linear = np.linspace(0.0,  0.13, 60)
y_linear = np.linspace(0.0, 0.13, 60)

X, Y = np.meshgrid(x_linear, y_linear)
points_to_eval = np.vstack((X.ravel(), Y.ravel()))

U = ssm.reduced_dynamics.predict(points_to_eval.T) # call the reduced vector field on the grid of points


UU = U[:,0].reshape(X.shape)
VV = U[:,1].reshape(Y.shape)

f = plt.figure(figsize = (8,6))
ax = f.add_subplot(111) 
streamplotobject = ax.streamplot(X, Y, UU, VV, density = 3, linewidth = 0.4, color='grey')

ax.plot((solStab1.y[0,:]), (solStab1.y[1,:]), '--', linewidth = 1.5, c='black')
ax.plot((solStab2.y[0,:]), (solStab2.y[1,:]), '--',  linewidth = 1.5, c='black')


ax.plot((solStab1.y[0,:4500]), (solStab1.y[1,:4500]), '-', linewidth = 1.5, c='black', label ='Stable manifold (predicted)')
ax.plot((solStab2.y[0,:4500]), (solStab2.y[1,:4500]), '-',  linewidth = 1.5, c='black')


ax.set_xlabel('$J$')
ax.set_ylabel('$K$')
flag = True
for x in ssm.emb_data['reduced_coordinates']:
    
    if flag:
        xx = np.hstack((x, np.array([0,0]).reshape(-1,1)))

        ax.plot(xx[0,:], xx[1,:], '-', c=linecolor, linewidth = 2, label= 'Training data')
    else:
        ax.plot(x[0,:], x[1,:], '-', c=linecolor, linewidth = 2)
    flag = False


ax.plot(0,0, 'o', c='black', label = 'Laminar flow')


ax.set_xlim(-0.0001,  0.13)
ax.set_ylim(-0.005, 0.13)
ax.plot(*lb_jk_coords, 'o', c='Crimson', label = 'Lower branch traveling wave (saddle)')
ax.legend(loc = 'upper left')
#shapely.plotting.plot_polygon(polygon_smooth, add_points=False, color=surfacecolor, alpha = 0.6, ax = ax) shapely command does not work
x_polygon, y_polygon = polygon_smooth.exterior.xy
ax.fill(x_polygon, y_polygon, color = surfacecolor, alpha = 0.6)
[<matplotlib.patches.Polygon at 0x7fd953386800>]
Image in a Jupyter notebook

Visualizing using the parametrization map

First select the points from the grid that are inside the given domain: 

points_inside = []
for p in points_to_eval.T:
    if polygon_smooth.contains(Point(p)):
        points_inside.append(p)

points_inside = np.array(points_inside)

To display an observable quantity, we define the relative energy as ΔE=12∥uHP−W(J,K)∥2, \Delta E = \frac{1}{2}\left \Vert \mathbf{u}_{HP} - \mathbf{W}(J,K) \right \Vert^2,

where uHP=(0,0,1−r2)\mathbf{u}_{HP} = (0,0,1-r^2) is the laminar velocity profile.

def calculate_energy(jk_coord):
    lifted =  ssm.decoder.predict(jk_coord)
    return np.linalg.norm(lifted, axis = 1) **2 /2

energies = calculate_energy(points_inside)

import matplotlib.tri as mtri
tri = mtri.Triangulation(points_inside[:,0], points_inside[:,1])

Then select the line-segments that were calculated by the streamplot. 

segments = []

# get the segments generated by the streamplot
for l in streamplotobject.lines.get_paths():
    v = np.array(l.vertices)
    if polygon_smooth.contains(Point(v[0,:])) and polygon_smooth.contains(Point(v[1,:])):
        segments.append(v)


energies_segments = []

for s in segments:
    energies_segments.append(calculate_energy(s))

%matplotlib inline

f = plt.figure()
ax = f.add_subplot(111, projection = '3d')
ax.plot_trisurf(points_inside[:,0], points_inside[:,1],
                energies,
                triangles = tri.triangles, 
                color=surfacecolor,
                alpha = 0.3, antialiased = True)


for i in range(len(energies_segments)):
    s = segments[i]
    e = energies_segments[i]
    ax.plot(s[:,0], s[:,1], e,'-', c='black', linewidth = 0.3)
    
ax.plot((solStab1.y[0,:]),
        (solStab1.y[1,:]),
        calculate_energy(solStab1.y.T),
        '--', linewidth = 1.5, c='black')
ax.plot((solStab2.y[0,:]),
        (solStab2.y[1,:]),
        calculate_energy(solStab2.y.T),
        '--',  linewidth = 1.5, c='black')


ax.plot((solStab1.y[0,:4500]),
        (solStab1.y[1,:4500]),
        calculate_energy(solStab1.y[:,:4500].T), '-', linewidth = 1.5, c='black', label ='Stable manifold (predicted)')
ax.plot((solStab2.y[0,:4500]),
        (solStab2.y[1,:4500]), 
        calculate_energy(solStab2.y[:,:4500].T),
        '-',  linewidth = 1.5, c='black')


    
    
ax.plot(*lb_jk_coords, calculate_energy(fixedPoint)[0], 'o', c='Crimson', label = 'Lower branch traveling wave (saddle)')
ax.plot(0,0,0, 'o', c='black', label ='Laminar state')
ax.set_xlabel('$J$')
ax.set_ylabel('$K$')
ax.set_zlabel('$\Delta E$')


ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))

ax.xaxis._axinfo["grid"].update({"linewidth":0.1})
ax.yaxis._axinfo["grid"].update({"linewidth":0.1})
ax.zaxis._axinfo["grid"].update({"linewidth":0.1})

ax.legend()
ax.view_init(azim=-90)

plt.show()
/usr/local/lib/python3.10/dist-packages/IPython/core/pylabtools.py:151: UserWarning: Creating legend with loc="best" can be slow with large amounts of data. fig.canvas.print_figure(bytes_io, **kw)
Image in a Jupyter notebook

Making a figure with the 2D and the 3D plot side by side

f = plt.figure(figsize = (16,6))
ax = f.add_subplot(121) 
ax2 = f.add_subplot(122, projection = '3d') 


## left figure
streamplotobject = ax.streamplot(X, Y, UU, VV, density = 3, linewidth = 0.4, color='grey')

ax.plot((solStab1.y[0,:]), (solStab1.y[1,:]), '--', linewidth = 1.5, c='black')
ax.plot((solStab2.y[0,:]), (solStab2.y[1,:]), '--',  linewidth = 1.5, c='black')


ax.plot((solStab1.y[0,:4500]), (solStab1.y[1,:4500]), '-', linewidth = 1.5, c='black', label ='Stable manifold (predicted)')
ax.plot((solStab2.y[0,:4500]), (solStab2.y[1,:4500]), '-',  linewidth = 1.5, c='black')


ax.set_xlabel('$J$')
ax.set_ylabel('$K$')
flag = True
for x in ssm.emb_data['reduced_coordinates']:
    
    if flag:
        xx = np.hstack((x, np.array([0,0]).reshape(-1,1)))

        ax.plot(xx[0,:], xx[1,:], '-', c=linecolor, linewidth = 2, label= 'Training data')
    else:
        ax.plot(x[0,:], x[1,:], '-', c=linecolor, linewidth = 2)
    flag = False


ax.plot(0,0, 'o', c='black', label = 'Laminar flow')


ax.set_xlim(-0.0001,  0.13)
ax.set_ylim(-0.005, 0.13)
ax.plot(*lb_jk_coords, 'o', c='Crimson', label = 'Lower branch traveling wave (saddle)')
ax.legend(loc = 'upper left')
#shapely.plotting.plot_polygon(polygon_smooth, add_points=False, color=surfacecolor, alpha = 0.6, ax = ax) shapely command does not work
x_polygon, y_polygon = polygon_smooth.exterior.xy
ax.fill(x_polygon, y_polygon, color = surfacecolor, alpha = 0.6)


## right figure
ax2.plot_trisurf(points_inside[:,0], points_inside[:,1],
                energies,
                triangles = tri.triangles, 
                color=surfacecolor,
                alpha = 0.3, antialiased = True)


for i in range(len(energies_segments)):
    s = segments[i]
    e = energies_segments[i]
    ax2.plot(s[:,0], s[:,1], e,'-', c='black', linewidth = 0.3)
    
ax2.plot((solStab1.y[0,:]),
        (solStab1.y[1,:]),
        calculate_energy(solStab1.y.T),
        '--', linewidth = 1.5, c='black')
ax2.plot((solStab2.y[0,:]),
        (solStab2.y[1,:]),
        calculate_energy(solStab2.y.T),
        '--',  linewidth = 1.5, c='black')


ax2.plot((solStab1.y[0,:4500]),
        (solStab1.y[1,:4500]),
        calculate_energy(solStab1.y[:,:4500].T), '-', linewidth = 1.5, c='black', label ='Stable manifold (predicted)')
ax2.plot((solStab2.y[0,:4500]),
        (solStab2.y[1,:4500]), 
        calculate_energy(solStab2.y[:,:4500].T),
        '-',  linewidth = 1.5, c='black')


    
    
ax2.plot(*lb_jk_coords, calculate_energy(fixedPoint)[0], 'o', c='Crimson', label = 'Lower branch traveling wave (saddle)')
ax2.plot(0,0,0, 'o', c='black', label ='Laminar state')
ax2.set_xlabel('$J$')
ax2.set_ylabel('$K$')
ax2.set_zlabel('$\Delta E$')


ax2.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax2.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax2.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))

ax2.xaxis._axinfo["grid"].update({"linewidth":0.1})
ax2.yaxis._axinfo["grid"].update({"linewidth":0.1})
ax2.zaxis._axinfo["grid"].update({"linewidth":0.1})
ax2.view_init(azim=-90)


## Also add the trajectories: 
for x in ssm.emb_data['reduced_coordinates']:
    
    if flag:
        xx = np.hstack((x, np.array([0,0]).reshape(-1,1)))
        
        ax2.plot(xx[0,:], xx[1,:], calculate_energy(xx.T), '-', c=linecolor, linewidth = 2, label= 'Training data')
    else:
        if x[0,-1] < 0.1: # do not plot trajectories that go to turbulence
            ax2.plot(x[0,:], x[1,:], calculate_energy(x.T),  '-', c=linecolor, linewidth = 2)
    flag = False

ax2.set_xlim(-0.05, 0.15)
ax2.set_ylim(-0.05, 0.15)
ax2.set_zlim(-0.05, 0.7)


ax2.legend()
<matplotlib.legend.Legend at 0x7fd9458652a0>
Image in a Jupyter notebook