This model version is a simple 2-layer channel QG atmosphere truncated at wavenumber 2 on a beta-plane with a simple orography (a montain and a valley).
More detail can be found in the articles:
First, setting the path and loading of some modules
import sys, os
sys.path.extend([os.path.abspath('../../')])
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import rc
rc('font',**{'family':'serif','sans-serif':['Times'],'size':14})
Initializing the random number generator (for reproducibility). -- Disable if needed.
np.random.seed(210217)
Importing the model's modules
from qgs.params.params import QgParams
from qgs.integrators.integrator import RungeKuttaIntegrator
from qgs.functions.tendencies import create_tendencies
from qgs.plotting.util import std_plot
Importing the Lyapunovs Estimators
from qgs.toolbox.lyapunov import LyapunovsEstimator, CovariantLyapunovsEstimator
from qgs.toolbox.lyapunov import _compute_backward_lyap_traj_jit, _compute_forward_lyap_traj_jit
General parameters
# Time parameters
dt = 0.1
# Saving the model state n steps
write_steps = 5
number_of_trajectories = 1
number_of_perturbed_trajectories = 10
Setting some model parameters
# Model parameters instantiation with some non-default specs
model_parameters = QgParams({'phi0_npi': np.deg2rad(50.)/np.pi, 'hd':0.3})
# Mode truncation at the wavenumber 2 in both x and y spatial coordinate
model_parameters.set_atmospheric_channel_fourier_modes(2, 2)
# Changing (increasing) the orography depth and the meridional temperature gradient
model_parameters.ground_params.set_orography(0.4, 1)
model_parameters.atemperature_params.set_thetas(0.2, 0)
# Printing the model's parameters
model_parameters.print_params()
Creating the tendencies function
f, Df = create_tendencies(model_parameters)
Defining an integrator
integrator = RungeKuttaIntegrator()
integrator.set_func(f)
Start on a random initial condition and integrate over a transient time to obtain an initial condition on the attractors
%%time
ic = np.random.rand(model_parameters.ndim)*0.1
integrator.integrate(0., 200000., dt, ic=ic, write_steps=0)
time, ic = integrator.get_trajectories()
Now integrate to obtain a trajectory on the attractor
%%time
integrator.integrate(0., 100000., dt, ic=ic, write_steps=write_steps)
reference_time, reference_traj = integrator.get_trajectories()
varx = 0
vary = 1
varz = 2
fig = plt.figure(figsize=(10, 8))
axi = fig.add_subplot(111, projection='3d')
axi.scatter(reference_traj[varx], reference_traj[vary], reference_traj[varz], s=0.2);
axi.set_xlabel('$'+model_parameters.latex_var_string[varx]+'$')
axi.set_ylabel('$'+model_parameters.latex_var_string[vary]+'$')
axi.set_zlabel('$'+model_parameters.latex_var_string[varz]+'$');
varx = 2
vary = 1
plt.figure(figsize=(10, 8))
plt.plot(reference_traj[varx], reference_traj[vary], marker='o', ms=0.07, ls='')
plt.xlabel('$'+model_parameters.latex_var_string[varx]+'$')
plt.ylabel('$'+model_parameters.latex_var_string[vary]+'$');
var = 1
plt.figure(figsize=(10, 8))
plt.plot(model_parameters.dimensional_time*reference_time, reference_traj[var])
plt.xlabel('time (days)')
plt.ylabel('$'+model_parameters.latex_var_string[var]+'$');
Here we compare the two methods used to compute the CLVs: the method 0 (Ginelli et al. algorithm) and the method 1 (subspaces intersection method). These methods are described in:
Covariant Lyapunovs Estimator
clvint = CovariantLyapunovsEstimator()
%%time
clvint.set_func(f, Df)
clvint.compute_clvs(0., 10000., 40000., 50000., 0.1, 0.1, ic, write_steps=1)
ctl0, ctraj0, cexp0, cvec0 = clvint.get_clvs()
clvint.terminate()
Plotting the spectrum for reference
plt.figure(figsize=(15, 4))
mean_exp = np.mean(cexp0, axis=-1)
x_pos = np.arange(1.,model_parameters.ndim+1,1)
plt.bar(x_pos, mean_exp)
plt.vlines(x_pos, -0.55, np.minimum(0.,mean_exp)-0.035, linestyles='dashdot', colors='tab:gray')
plt.xticks(x_pos, map(str,range(1, model_parameters.ndim+1,1)))
yt=[-0.5,-0.4,-0.3,-0.2,-0.1,0.,0.1]
plt.yticks(yt, map(str,yt))
plt.xlim(x_pos[0]-1., x_pos[-1]+1.)
plt.ylim(np.min(mean_exp)-0.1, np.max(mean_exp)+0.1)
plt.ylabel("Lyapunov exponent");
plt.xlabel("Index of the Lyapunov exponent");
Computing the BLVs and FLVs
pretime = ctl0[:100001]
time = ctl0[100000:200001]
posttime = ctl0[200000:]
backtraj = ctraj0[..., :200001][np.newaxis, ...]
forwtraj = ctraj0[..., 100000:][np.newaxis, ...]
cvec0 = cvec0[..., 100000:200001]
ftraj, fexp, fvec = _compute_forward_lyap_traj_jit(f, Df, time, posttime, forwtraj, 0.1, model_parameters.ndim, 1, False, 1, clvint.b, clvint.c, clvint.a)
btraj, bexp, bvec = _compute_backward_lyap_traj_jit(f, Df, pretime, time, backtraj, 0.1, model_parameters.ndim, 1, False, 1, clvint.b, clvint.c, clvint.a)
Computing the subspaces intersections
ctraj1 = forwtraj[..., :100001]
n_records = ctraj1.shape[-1]
cvec1 = np.zeros((model_parameters.ndim, model_parameters.ndim, n_records))
i_traj = 0
for ti in range(n_records):
for j in range(model_parameters.ndim):
u, z, w = np.linalg.svd(bvec[i_traj, :, :j+1, ti].T @ fvec[i_traj, :, :model_parameters.ndim-j, ti])
basis = bvec[i_traj, :, :j+1, ti] @ u
cvec1[:, j, ti] = basis[:, 0]
Obtained by method 0
cvec0[:,0,0]
Obtained by method 1
cvec1[:,0,0]
Each component is plotted with a different color
vars = slice(0, model_parameters.ndim)
fig = plt.figure(figsize=(20, int(model_parameters.ndim*8/2)), constrained_layout=False)
grid = fig.add_gridspec(int(model_parameters.ndim/2), 2)
axs = grid.subplots()
for vec, ax in enumerate(axs.flatten()):
ax.plot(model_parameters.dimensional_time*time, (np.abs(cvec0[vars,vec,:])-np.abs(cvec1[vars,vec,:])).T);
ax.set_xlabel('time (days)')
ax.set_title('CLV '+str(vec+1))