Manual setting of the basis

In this example, we show how to setup a user-defined basis. It is done here for the ocean, but the approach is the same for the other components. We will project the ocean equations on four modes proposed in

  • S. Pierini. Low-frequency variability, coherence resonance, and phase selection in a low-order model of the wind-driven ocean circulation. Journal of Physical Oceanography, 41(9):1585–1604, 2011. doi:10.1175/JPO-D-10-05018.1.

  • S. Vannitsem and L. De Cruz. A 24-variable low-order coupled ocean–atmosphere model: OA-QG-WS v2. Geoscientific Model Development, 7(2):649–662, 2014. doi:10.5194/gmd-7-649-2014.

These four modes are given by

  • \(\tilde\phi_1(x,y)=2 \, e^{−\alpha x} \sin(n x/2) \sin(y)\)

  • \(\tilde\phi_2(x,y)=2 \, e^{−\alpha x} \sin(n x) \sin(y)\)

  • \(\tilde\phi_3(x,y)=2 \, e^{−\alpha x} \sin(n x/2) \sin(2y)\)

  • \(\tilde\phi_4(x,y)=2 \, e^{−\alpha x} \sin(n x) \sin(2y)\)

This ocean is then connected to a channel atmosphere using a symbolic basis of functions.

Modules import

Loading of some modules…

import sys, os, warnings
sys.path.extend([os.path.abspath('../')])
import numpy as np
import matplotlib.pyplot as plt
from numba import njit

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.basis.base import SymbolicBasis
from qgs.inner_products.definition import StandardSymbolicInnerProductDefinition
from qgs.inner_products.symbolic import AtmosphericSymbolicInnerProducts, OceanicSymbolicInnerProducts
from qgs.tensors.qgtensor import QgsTensor
from qgs.functions.sparse_mul import sparse_mul3
from qgs.integrators.integrator import RungeKuttaIntegrator

and diagnostics

from qgs.diagnostics.streamfunctions import MiddleAtmosphericStreamfunctionDiagnostic, OceanicLayerStreamfunctionDiagnostic
from qgs.diagnostics.temperatures import MiddleAtmosphericTemperatureDiagnostic, OceanicLayerTemperatureDiagnostic
from qgs.diagnostics.variables import VariablesDiagnostic, GeopotentialHeightDifferenceDiagnostic
from qgs.diagnostics.multi import MultiDiagnostic

and some SymPy functions

from sympy import symbols, sin, exp, pi

Systems definition

General parameters

# Time parameters
dt = 0.1
# Saving the model state n steps
write_steps = 100

number_of_trajectories = 1

Setting some model parameters and setting the atmosphere basis

# Model parameters instantiation with some non-default specs
model_parameters = QgParams({'n': 1.5})

# Mode truncation at the wavenumber 2 in both x and y spatial
# coordinates for the atmosphere
model_parameters.set_atmospheric_channel_fourier_modes(2, 2, mode="symbolic")

Creating the ocean basis

ocean_basis = SymbolicBasis()
x, y = symbols('x y')  # x and y coordinates on the model's spatial domain
n, al = symbols('n al', real=True, nonnegative=True)  # aspect ratio and alpha coefficients
for i in range(1, 3):
    for j in range(1, 3):
        ocean_basis.functions.append(2 * exp(- al * x) * sin(j * n * x / 2) * sin(i * y))

Setting the value of the parameter α to a certain value (here α=1). Please note that the α is then an extrinsic parameter of the model that you have to specify through a substitution:

ocean_basis.substitutions.append((al, 1.))
ocean_basis.substitutions.append((n, model_parameters.scale_params.n))

Setting now the ocean basis

model_parameters.set_oceanic_modes(ocean_basis)

Additionally, for these particular ocean basis functions, a special inner product needs to be defined instead of the standard one proposed. We thus set up a new definition for the ocean component as in the publication linked above:

class UserInnerProductDefinition(StandardSymbolicInnerProductDefinition):

    def symbolic_inner_product(self, S, G, symbolic_expr=False, integrand=False):
        """Function defining the inner product to be computed symbolically:
        :math:`(S, G) = \\frac{n}{2\\pi^2}\\int_0^\\pi\\int_0^{2\\pi/n}  e^{2 \\alpha x} \\,  S(x,y)\\, G(x,y)\\, \\mathrm{d} x \\, \\mathrm{d} y`.

        Parameters
        ----------
        S: Sympy expression
            Left-hand side function of the product.
        G: Sympy expression
            Right-hand side function of the product.
        symbolic_expr: bool, optional
            If `True`, return the integral as a symbolic expression object. Else, return the integral performed symbolically.
        integrand: bool, optional
            If `True`, return the integrand of the integral and its integration limits as a list of symbolic expression object. Else, return the integral performed symbolically.

        Returns
        -------
        Sympy expression
            The result of the symbolic integration
        """
        expr = (n / (2 * pi ** 2)) * exp(2 * al * x) * S * G
        if integrand:
            return expr, (x, 0, 2 * pi / n), (y, 0, pi)
        else:
            return self.integrate_over_domain(self.optimizer(expr), symbolic_expr=symbolic_expr)

Finally setting some other model’s parameters:

# Setting MAOOAM parameters according to the publication linked above
model_parameters.set_params({'kd': 0.029, 'kdp': 0.029, 'r': 1.0e-7,
                             'h': 136.5, 'd': 1.1e-7})
model_parameters.atemperature_params.set_params({'eps': 0.76, 'T0': 289.3,
                                                 'hlambda': 15.06})
model_parameters.gotemperature_params.set_params({'gamma': 5.6e8, 'T0': 301.46})

Setting the short-wave radiation component as in the publication above: \(C_{\text{a},1}\) and \(C_{\text{o},1}\)

model_parameters.atemperature_params.set_insolation(103.3333, 0)
model_parameters.gotemperature_params.set_insolation(310, 0)

Printing the model’s parameters

model_parameters.print_params()
Qgs parameters summary
======================

General Parameters:
'time_unit': days,
'rr': 287.058  [J][kg^-1][K^-1]  (gas constant of dry air),
'sb': 5.67e-08  [J][m^-2][s^-1][K^-4]  (Stefan-Boltzmann constant),

Scale Parameters:
'scale': 5000000.0  [m]  (characteristic space scale (L*pi)),
'f0': 0.0001032  [s^-1]  (Coriolis parameter at the middle of the domain),
'n': 1.5    (aspect ratio (n = 2 L_y / L_x)),
'rra': 6370000.0  [m]  (earth radius),
'phi0_npi': 0.25    (latitude expressed in fraction of pi),
'deltap': 50000.0  [Pa]  (pressure difference between the two atmospheric layers),

Atmospheric Parameters:
'kd': 0.029  [nondim]  (atmosphere bottom friction coefficient),
'kdp': 0.029  [nondim]  (atmosphere internal friction coefficient),
'sigma': 0.2  [nondim]  (static stability of the atmosphere),

Atmospheric Temperature Parameters:
'gamma': 10000000.0  [J][m^-2][K^-1]  (specific heat capacity of the atmosphere),
'C[1]': 103.3333  [W][m^-2]  (spectral component 1 of the short-wave radiation of the atmosphere),
'C[2]': 0.0  [W][m^-2]  (spectral component 2 of the short-wave radiation of the atmosphere),
'C[3]': 0.0  [W][m^-2]  (spectral component 3 of the short-wave radiation of the atmosphere),
'C[4]': 0.0  [W][m^-2]  (spectral component 4 of the short-wave radiation of the atmosphere),
'C[5]': 0.0  [W][m^-2]  (spectral component 5 of the short-wave radiation of the atmosphere),
'C[6]': 0.0  [W][m^-2]  (spectral component 6 of the short-wave radiation of the atmosphere),
'C[7]': 0.0  [W][m^-2]  (spectral component 7 of the short-wave radiation of the atmosphere),
'C[8]': 0.0  [W][m^-2]  (spectral component 8 of the short-wave radiation of the atmosphere),
'C[9]': 0.0  [W][m^-2]  (spectral component 9 of the short-wave radiation of the atmosphere),
'C[10]': 0.0  [W][m^-2]  (spectral component 10 of the short-wave radiation of the atmosphere),
'eps': 0.76    (emissivity coefficient for the grey-body atmosphere),
'T0': 289.3  [K]  (stationary solution for the 0-th order atmospheric temperature),
'sc': 1.0    (ratio of surface to atmosphere temperature),
'hlambda': 15.06  [W][m^-2][K^-1]  (sensible+turbulent heat exchange between ocean and atmosphere),

Oceanic Parameters:
'gp': 0.031  [m][s^-2]  (reduced gravity),
'r': 1e-07  [s^-1]  (frictional coefficient at the bottom of the ocean),
'h': 136.5  [m]  (depth of the water layer of the ocean),
'd': 1.1e-07  [s^-1]  (strength of the ocean-atmosphere mechanical coupling),

Oceanic Temperature Parameters:
'gamma': 560000000.0  [J][m^-2][K^-1]  (specific heat capacity of the ocean),
'C[1]': 310.0  [W][m^-2]  (spectral component 1 of the short-wave radiation of the ocean),
'C[2]': 0.0  [W][m^-2]  (spectral component 2 of the short-wave radiation of the ocean),
'C[3]': 0.0  [W][m^-2]  (spectral component 3 of the short-wave radiation of the ocean),
'C[4]': 0.0  [W][m^-2]  (spectral component 4 of the short-wave radiation of the ocean),
'C[5]': 0.0  [W][m^-2]  (spectral component 5 of the short-wave radiation of the ocean),
'C[6]': 0.0  [W][m^-2]  (spectral component 6 of the short-wave radiation of the ocean),
'C[7]': 0.0  [W][m^-2]  (spectral component 7 of the short-wave radiation of the ocean),
'C[8]': 0.0  [W][m^-2]  (spectral component 8 of the short-wave radiation of the ocean),
'C[9]': 0.0  [W][m^-2]  (spectral component 9 of the short-wave radiation of the ocean),
'C[10]': 0.0  [W][m^-2]  (spectral component 10 of the short-wave radiation of the ocean),
'T0': 301.46  [K]  (stationary solution for the 0-th order oceanic temperature),

We now construct the tendencies of the model by first constructing the ocean and atmosphere inner products objects. In addition, a inner product definition instance defined above must be passed to the ocean inner products object:

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    ip = UserInnerProductDefinition()

    aip = AtmosphericSymbolicInnerProducts(model_parameters)
    oip = OceanicSymbolicInnerProducts(model_parameters, inner_product_definition=ip)

and finally we create manually the tendencies function, first by creating the tensor object:

aotensor = QgsTensor(model_parameters, aip, oip)

and then the Python-Numba callable for the model’s tendencies \(\boldsymbol{f}\) :

coo = aotensor.tensor.coords.T
val = aotensor.tensor.data

@njit
def f(t, x):
    xx = np.concatenate((np.full((1,), 1.), x))
    xr = sparse_mul3(coo, val, xx, xx)

    return xr[1:]

Time integration

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

## Might take several minutes, depending on your cpu computational power.
ic = np.random.rand(model_parameters.ndim)*0.0001
integrator.integrate(0., 10000000., dt, ic=ic, write_steps=0)
time, ic = integrator.get_trajectories()

Now integrate to obtain a trajectory on the attractor

integrator.integrate(0., 1000000., dt, ic=ic, write_steps=write_steps)
reference_time, reference_traj = integrator.get_trajectories()

Plotting the result in 3D and 2D

varx = 21
vary = 25
varz = 0

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]+'$');
../../_images/manual_basis_51_0.png
varx = 21
vary = 25
plt.figure(figsize=(10, 8))

plt.plot(reference_traj[varx], reference_traj[vary], marker='o', ms=0.1, ls='')

plt.xlabel('$'+model_parameters.latex_var_string[varx]+'$')
plt.ylabel('$'+model_parameters.latex_var_string[vary]+'$');
../../_images/manual_basis_52_0.png
var = 21
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]+'$');
../../_images/manual_basis_53_0.png

Showing the resulting fields (animation)

Here, we want to show that the diagnostics adapt to the manually set basis.

One can use objects called Diagnostic to plot movies of the model fields. Here we show how to use them and plot the simultaneous time evolution of the variables \(\psi_{{\rm a}, 1}\), \(\psi_{{\rm o}, 2}\) and \(\theta_{{\rm o}, 2}\), with the corresponding atmospheric and oceanic streamfunctions and temperature at 500 hPa.

Creating the diagnostics (for field plots, we must specify the grid step):

  • For the 500hPa geopotential height:

psi_a = MiddleAtmosphericStreamfunctionDiagnostic(model_parameters, delta_x=0.1, delta_y=0.1, geopotential=True)
  • For the 500hPa atmospheric temperature:

theta_a = MiddleAtmosphericTemperatureDiagnostic(model_parameters, delta_x=0.1, delta_y=0.1)
  • For the ocean streamfunction:

psi_o = OceanicLayerStreamfunctionDiagnostic(model_parameters, delta_x=0.1, delta_y=0.1, conserved=False)
  • For the ocean temperature:

theta_o = OceanicLayerTemperatureDiagnostic(model_parameters, delta_x=0.1, delta_y=0.1)
  • For the nondimensional variables \(\psi_{{\rm a}, 1}\), \(\psi_{{\rm o}, 2}\) and \(\theta_{{\rm o}, 2}\):

variable_nondim = VariablesDiagnostic([21, 25, 0], model_parameters, False)
  • For the geopotential height difference between North and South:

geopot_dim = GeopotentialHeightDifferenceDiagnostic([[[np.pi/model_parameters.scale_params.n, np.pi/4], [np.pi/model_parameters.scale_params.n, 3*np.pi/4]]],
                                                    model_parameters, True)
# setting also the background
background = VariablesDiagnostic([21, 25, 0], model_parameters, False)
background.set_data(reference_time, reference_traj)

Selecting a subset of the data to plot:

stride = 10
time = reference_time[10000:10000+5200*stride:stride]
traj = reference_traj[:, 10000:10000+5200*stride:stride]

Creating a multi diagnostic with both:

m = MultiDiagnostic(3,2)
m.add_diagnostic(geopot_dim, diagnostic_kwargs={'style':'moving-timeserie'})
m.add_diagnostic(variable_nondim, diagnostic_kwargs={'show_time':False, 'background': background, 'style':'3Dscatter'}, plot_kwargs={'ms': 0.2})
m.add_diagnostic(theta_a, diagnostic_kwargs={'show_time':False})
m.add_diagnostic(psi_a, diagnostic_kwargs={'show_time':False})
m.add_diagnostic(theta_o, diagnostic_kwargs={'show_time':False})
m.add_diagnostic(psi_o, diagnostic_kwargs={'show_time':False})

m.set_data(time, traj)

and showing finally a movie:

m.movie(figsize=(15,18), anim_kwargs={'interval': 100, 'frames':1000})